// Swarm library. Copyright (C) 1996 Santa Fe Institute.
//  Swarm library. Copyright (C) 1996 Santa Fe Institute.
// This library is distributed without any warranty; 
// without even the implied warranty of merchantability 
// or fitness for a particular purpose.
// See file LICENSE for details and terms of copying.

/*
Name:            NormalDistribution.m
Description:     Normal (Gaussian) distribution returning double values
Library:         random
Original Author: Manor Askenazi
Modified by:     Sven Thommesen
Date:            1997-01-15
*/

/*
123456789|123456789|123456789|123456789|123456789|123456789|123456789|123456789|
*/

// 
// Methods to generate the next pseudo-random number in a stream:
// 
// // Using the bit generator to create Uniform numbers.
//
// // Bit generators return values in the interval [0,maxValue].
// // For the implemented generators, maxValue is:
// //    PMMLCG:  2,147,483,645 = 2^31 - 3
// //    SWB:     4,294,967,295 = 2^32 - 1
// //    LCG:     4,294,967,295 = 2^32 - 1
// //    ACG:     4,294,967,295 = 2^32 - 1
// //    SCG:       999,999,999 = 10^9 - 1 (29.9 bits)
// 


#import <collections.h>
#import <math.h>
#import <random/NormalDistribution.h>


@implementation NormalDistribution

// data struct used by getState and setState:
//
typedef struct {
   // Object identification:
   unsigned magic;
   // unsigned generator_magic;	// later
   // ProbabilityDistribution variables:
   unsigned stateSize;
   BOOL optionsInitialized;
   unsigned currentCount;
   unsigned generatorMax;
   double maxDivisor;
   // Distribution specific data:
   double theMean;
   double theVariance;
   double theStdDev;
   BOOL stored;
   double stored_double;
} state_struct_t;

// ----- internal methods -----

-resetState {

// Called by setGenerator in the superclass 
// and setMean:Variance below

   stored = NO;
   stored_double = 0.0;
   currentCount = 0;

   return self;
}

-initState {
   UnsignedOverlayDouble duuTemp;
   // volatile double x;
   // double y;

// Called from createBegin in the superclass

   stateSize = sizeof(state_struct_t);

   theMean       = 0.0;
   theVariance   = 0.0;
   theStdDev     = 0.0;
  
   stored        = NO;
   stored_double = 0.0;
   currentCount  = 0;

// State variables are allocated as instance variables
// and initialized in resetState above.

// ----------------------------
   // Now for the libg++ magic:

   if ( sizeof(double) != 2*sizeof(unsigned) )
   [InvalidCombination raiseEvent:
   " double <> 2 unsigneds: math will fail!!! \n"];

// #if_IEEE == 1
// #ifdef _IEEE

   //   printf(" initState: IEEE is defined \n");

   duuTemp.d = 1.5;
   if (duuTemp.u[1] == 0) {	// sun word order?
     //printf(" initState: sun word order \n");
     duuTemp.u[0] = 0x3fffffff;
     duuTemp.u[1] = 0xffffffff;
   } else {			// encore word order?
     //printf(" initState: encore word order \n");
     duuTemp.u[0] = 0xffffffff;
     duuTemp.u[1] = 0x3fffffff;
   }

/*
#else IEEE

   printf(" initState: IEEE is *not* defined \n");

   x = 1.0;
   y = 0.5;

   do {
     duuTemp.d = x;
     x += y;
     y *= 0.5;
   } while ( x != duuTemp.d && x < 2.0);

#endif IEEE
*/

   duuMask.d = 1.0;
   duuMask.u[0] ^= duuTemp.u[0];
   duuMask.u[1] ^= duuTemp.u[1];

// ----------------------------

   return self;
}

// This method returns values in the interval [0.0,1.0):

-(double) rDouble {
  double rdValue;
  unsigned bitValue;

  bitValue = [randomGenerator getUnsignedSample];
  rdValue = (double) (bitValue / maxDivisor);

  return rdValue;
}

// This method returns values in the interval (0.0,1.0):

-(double) nzrDouble {
  double rdValue;
  unsigned bitValue;

  do { bitValue = [randomGenerator getUnsignedSample]; }
  while (bitValue == 0);		// So we don't return 0.0

  rdValue = (double) (bitValue / maxDivisor);

  return rdValue;
}

// This method fills each double with two successive generator 
// values (64 bits). It will never return 0.0.

-(double) duuDouble {
   unsigned bitValue0, bitValue1;
   UnsignedOverlayDouble duuValue;

   bitValue0 = [randomGenerator getUnsignedSample];
   bitValue1 = [randomGenerator getUnsignedSample];

//   printf(" duuValue: generator inputs were %u and %u \n", 
//	bitValue0, bitValue1);

   duuValue.d = 1.0;
   duuValue.u[0] |= (bitValue0 & duuMask.u[0]);
   duuValue.u[1] |= (bitValue1 & duuMask.u[1]);
   duuValue.d -= 1.0;

   if ((duuValue.d >= 1.0) || (duuValue.d < 0))
   [InvalidCombination raiseEvent:
   "duuDouble: result was not in [0.0,1.0) -- sigh... \n"];

//   printf(" duuValue: d0=%60.50f\n", duuValue.d);

   return duuValue.d;
}

// ----- protocol Normal -----

+create: (id) aZone setGenerator: (id) generator
	setMean: (double) mean setVariance: (double) variance {
   NormalDistribution * aDistribution;

   if (variance == 0.0)
   [InvalidCombination raiseEvent:
   "NormalDistribution: setting Variance = 0 not allowed !\n"];

   aDistribution = [ super create: aZone setGenerator: generator ];

   aDistribution->theMean = mean;
   aDistribution->theVariance = variance;
   aDistribution->theStdDev = sqrt(variance);

   // This object is now fixed:

   aDistribution->optionsInitialized = YES;

   [ aDistribution resetState ];

   return aDistribution;
}

-(void) setMean: (double) mean setVariance: (double) variance {

   if (optionsInitialized)
   [InvalidCombination raiseEvent:
   "NormalDistribution: resetting parameters not allowed\n"];

   if (variance == 0.0)
   [InvalidCombination raiseEvent:
   "NormalDistribution: setting Variance = 0 not allowed !\n"];

   theMean = mean;
   theVariance = variance;
   theStdDev = sqrt(variance);

   // This object is now fixed:

   optionsInitialized = YES;

   [self resetState];

   //   return self;
}

-(double) getMean {
   return theMean;
}

-(double) getVariance {
   return theVariance;
}

-(double) getStdDev {
  return theStdDev;
}

// ==================================================

// Now do the transformation from Uniform to Normal.
// Both versions do all calculations using
// doubles, returning doubles. Input values
// for Mean and Standard Deviation are doubles.

// Source for the algorithm used:
// 'Numerical recipes in C' by W.H.Press et al.

-(double) getSampleWithMean: (double) mean withVariance: (double) variance {
   double fac,radius,v1,v2 ;
   double rd1Value, rd2Value;
   double stdDev;

   if (optionsInitialized)
   [InvalidCombination raiseEvent:
   "NormalDistribution: getSampleWithMean:withVariance: options already initialized\n"];

   if (variance == 0.0)
   [InvalidCombination raiseEvent:
   "NormalDistribution: setting Variance = 0 not allowed !\n"];

   stdDev = sqrt(variance);

   currentCount++ ;

  if (stored) {
    stored = NO ;
    return stored_double*stdDev + mean ;
  } else {
    do {
      // rd1Value = [self rDouble];
      // rd2Value = [self rDouble];
      rd1Value = [self duuDouble];
      rd2Value = [self duuDouble];
      v1 = (2.0 * rd1Value) - 1.0 ;
      v2 = (2.0 * rd2Value) - 1.0 ;
      radius = v1*v1 + v2*v2 ; 
        } while (radius >= 1.0) ;
    fac = sqrt(-2.0*log(radius)/radius) ;
    stored_double = v1*fac ;
    stored = YES;
    return v2*fac*stdDev + mean ;
  }
}

// ----- protocol DoubleDistribution -----

-(double) getDoubleSample {
   double fac,radius,v1,v2 ;
   double rd1Value, rd2Value;

   if (!optionsInitialized)
   [InvalidCombination raiseEvent:
   "NormalDistribution: getDoubleSample: options have not been initialized\n"];

   currentCount++ ;

  if (stored) {
    stored = NO ;
    return stored_double*theStdDev + theMean ;
  } else {
    do {
      // rd1Value = [self rDouble];
      // rd2Value = [self rDouble];
      rd1Value = [self duuDouble];
      rd2Value = [self duuDouble];
      v1 = (2.0 * rd1Value) - 1.0 ;
      v2 = (2.0 * rd2Value) - 1.0 ;
      radius = v1*v1 + v2*v2 ; 
        } while (radius >= 1.0) ;
    fac = sqrt(-2.0*log(radius)/radius) ;
    stored_double = v1*fac ;
    stored = YES;
    return v2*fac*theStdDev + theMean ;
  }
}

// ----- protocol InternalState -----

-(void) getState: (void *) stateBuf {		// override superclass method
   state_struct_t * internalState;
   // unsigned generator_magic;

  // obtain my generator's magic number:
  // (no way to do that yet)

  // recast the caller's pointer:
  internalState = (state_struct_t *) stateBuf;

  // fill the caller's buffer with state data:
  // object identification:
  internalState->magic = NORMALDISTMAGIC;
  // internalState->generator_magic = generator_magic;
  // ProbabilityDistribution data:
  internalState->stateSize = stateSize;
  internalState->optionsInitialized = optionsInitialized;
  internalState->currentCount = currentCount;
  internalState->generatorMax = generatorMax;
  internalState->maxDivisor = maxDivisor;
  // distribution specific data:
  internalState->theMean = theMean;
  internalState->theVariance = theVariance;
  internalState->theStdDev = theStdDev;
  internalState->stored = stored;
  internalState->stored_double = stored_double;

  // nothing is returned from a (void) function

}

-(void) setState: (void *) stateBuf {		// override superclass method
   state_struct_t * internalState;

  // obtain my generator's magic number:
  // (no way to do that yet)

  // recast the caller's pointer:
  internalState = (state_struct_t *) stateBuf;

  // TEST the integrity of the external data:
  if (    (internalState->magic     != NORMALDISTMAGIC)
       || (internalState->stateSize != stateSize)
//     || (internalState->generator_magic != generator_magic)
     )
  [InvalidCombination raiseEvent:
  "NormalDistribution: you are passing bad data to setState!\n %u %u\n",
   internalState->magic, internalState->stateSize];

  // set internal state from data in caller's buffer:
  // ProbabilityDistribution data:
    // stateSize = internalState->stateSize;
  optionsInitialized = internalState->optionsInitialized;
  currentCount = internalState->currentCount;
    // generatorMax = internalState->generatorMax;
    // maxDivisor = internalState->maxDivisor;
  // distribution specific data:
  theMean       = internalState->theMean;
  theVariance   = internalState->theVariance;
  theStdDev     = internalState->theStdDev;
  stored        = internalState->stored;
  stored_double = internalState->stored_double;

  // nothing is returned from a (void) function

}

// ----- temporary methods -----

- (void) describe: outStream { 
  char buffer[200];

  (void)sprintf(buffer," NormalDistribution describe: outStream: \n");
  (void)sprintf(buffer,"       theMean = %f\n", theMean);
  (void)sprintf(buffer,"   theVariance = %f\n", theVariance);
  (void)sprintf(buffer,"     theStdDev = %f\n", theStdDev);
  (void)sprintf(buffer,"       stored? = %d\n", stored);
  (void)sprintf(buffer," stored_double = %f\n", stored_double);

  [ super describe: outStream ];
  [outStream catC: buffer];
  [outStream catC: "\n"];

  //return self;
}

-(int) verifySelf {

   return 1;
}

@end
