// 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:            ExponentialDistribution.m
Description:     Exponential distribution returning double values
Library:         random
Original Author: 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/ExponentialDistribution.h>


@implementation ExponentialDistribution

// 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;
   double duuMask;
   // Distribution specific data:
   double theMean;
} state_struct_t;

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

-resetState {

// Called by setGenerator in the superclass
// and by setMean below.

   currentCount = 0;

   return self;
}

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

// Called by createBegin in the superclass

   stateSize = sizeof(state_struct_t);

   theMean       = 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 Exponential -----

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

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

   aDistribution->theMean = mean;

   // This object is now fixed:

   aDistribution->optionsInitialized = YES;

   [ aDistribution resetState ];

   return aDistribution;
}

-(void) setMean: (double) mean {

   if (optionsInitialized)
   [InvalidCombination raiseEvent:
   "ExponentialDistribution: setting parameters more than once not allowed\n"];

   theMean = mean;

   // This object is now fixed:

   optionsInitialized = YES;

   [self resetState];

   //   return self;
}

-(double) getMean {
   return theMean;
}

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

-(double) getSampleWithMean: (double) mean {
   double xpon;
   double rdValue;

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

   currentCount++ ;

   // Generate a value in [0.0,1.0)
   // rdValue = [self rDouble];

   rdValue = [self duuDouble];

   // Transform the uniform value:
   // (Code from Russell 1992)

   xpon = -mean*log(rdValue);	// need rdValue > 0 !!!

   return xpon;

}

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

-(double) getDoubleSample {
   double xpon;
   double rdValue;

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

   currentCount++ ;

   // Generate a value in [0.0,1.0)
   // rdValue = [self rDouble];
   rdValue = [self duuDouble];

   // Transform the uniform value:
   // (Code from Russell 1992)

   xpon = -theMean*log(rdValue);	// need rdValue > 0 !!!

   return xpon;

}

// ----- 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 = EXPONENTIALDISTMAGIC;
  // internalState->generator_magic = generator_magic;
  // ProbabilityDistribution data:
  internalState->stateSize = stateSize;
  internalState->optionsInitialized = optionsInitialized;
  internalState->currentCount = currentCount;
  internalState->generatorMax = generatorMax;
  internalState->maxDivisor = maxDivisor;
  internalState->duuMask = duuMask.d;
  // distribution specific data:
  internalState->theMean = theMean;

  // 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     != EXPONENTIALDISTMAGIC)
       || (internalState->stateSize != stateSize)
//     || (internalState->generator_magic != generator_magic)
     )
  [InvalidCombination raiseEvent:
  "ExponentialDistribution: 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;
    // duuMask.d = internalState->duuMask;
  // distribution specific data:
  theMean       = internalState->theMean;

  // nothing is returned from a (void) function

}

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

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

  (void)sprintf(buffer," ExponentialDistribution: \n");
  (void)sprintf(buffer,"       theMean = %f\n", theMean);

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

  //  return self;
}

-(int) verifySelf {

   return 1;
}

@end
