// 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:            GammaDistribution.m
Description:     Gamma 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/GammaDistribution.h>


@implementation GammaDistribution

// 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 theAlpha;
   double theBeta;
} state_struct_t;

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

-resetState {

// Called by setGenerator in the superclass
// and by setAlpha:Beta: below.

   currentCount = 0;

   return self;
}

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

// Called from createBegin in the superclass.

   stateSize = sizeof(state_struct_t);

   theAlpha      = 0.0;
   theBeta       = 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 Gamma -----

+create: (id) aZone setGenerator: (id) generator
	setAlpha: (double) alpha setBeta: (double) beta {
   GammaDistribution * aDistribution;

   if (alpha <= 0.0)
   [InvalidCombination raiseEvent:
   "GammaDistribution: setting alpha <= 0.0 not supported\n"];

   if (beta <= 0.0)
   [InvalidCombination raiseEvent:
   "GammaDistribution: setting beta <= 0.0 not supported\n"];

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

   aDistribution->theAlpha = alpha;
   aDistribution->theBeta  = beta;

   // This object is now fixed:

   aDistribution->optionsInitialized = YES;

   [ aDistribution resetState ];

   return aDistribution;
}

-(void) setAlpha: (double) alpha setBeta: (double) beta {

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

   if (alpha <= 0.0)
   [InvalidCombination raiseEvent:
   "GammaDistribution: setting alpha <= 0.0 not supported\n"];

   if (beta <= 0.0)
   [InvalidCombination raiseEvent:
   "GammaDistribution: setting beta <= 0.0 not supported\n"];

   theAlpha = alpha;
   theBeta  = beta;

   // This object is now fixed:

   optionsInitialized = YES;

   [self resetState];

   //   return self;
}

-(double) getAlpha {
   return theAlpha;
}

-(double) getBeta {
   return theBeta;
}

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

// Method that returns an Exponential:
-(double) getExponentialWithMean: (double) theMean 
{
  // return theMean * (-log([self nzrDouble]));
  return theMean * (-log([self duuDouble]));
}

-(double) getSampleWithAlpha: (double) alpha withBeta: (double) beta {
   // unsigned bitValue;
   // double rdValue;
   double  x, avg, am, e, s, v1, v2, y;
   // int i;
   // double u1, u2;

   if (optionsInitialized)
   [InvalidCombination raiseEvent:
"GammaDistribution: getSampleWithAlpha:withBeta: options already initialized\n"
   ];

   if (alpha <= 0.0)
   [InvalidCombination raiseEvent:
   "GammaDistribution: setting alpha <= 0.0 not supported\n"];

   if (beta <= 0.0)
   [InvalidCombination raiseEvent:
   "GammaDistribution: setting beta <= 0.0 not supported\n"];

   currentCount++ ;

  // Now transform Uniform to Gamma:
  // code from  Watkins (1994)

  // For alpha == 1, just use exponential:

  if (alpha == 1.0) 
  {
    return ([self getExponentialWithMean: 1.0] / beta);
  }

  // For alpha < 1, rejection method does not work,
  // so use alternative:
 
  if (alpha < 1.0) 
  {
    do 
    {
      x = pow([self duuDouble], 1.0/alpha);
      y = pow([self duuDouble], 1.0/(1.0-alpha));
    } while (x+y > 1.0);

    x = x/(x+y);
    y = [self getExponentialWithMean: 1.0];

    return (x*y/beta);
  }

  // For alpha > 1.0, use rejection method:
   
  do 
  {
    do 
    {
      do 
      {
        v1 = 2.0 * [self duuDouble] - 1.0;
        v2 = 2.0 * [self duuDouble] - 1.0;
      } while (v1 * v1 + v2 * v2 > 1.0);
      y = v2 / v1;
      am = alpha - 1.0;
      s = sqrt(2.0 * am + 1.0);
      avg = s * y + am;
    } while (avg <= 0.0);
    e = (1.0 + y * y) * exp(am * log(avg / am) - s * y);
  } while ([self duuDouble] > e);

  return (avg / beta) ;


}

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

-(double) getDoubleSample {
   // unsigned bitValue;
   // double rdValue;
   double  x, avg, am, e, s, v1, v2, y;
   // int i;
   // double u1, u2;

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

   currentCount++ ;

  // Now transform Uniform to Gamma:
  // code from Watkins (1994)

  // For alpha == 1, just use exponential:

  if (theAlpha == 1.0) 
  {
    return ([self getExponentialWithMean: 1.0] / theBeta);
  }

  // For alpha < 1, rejection method does not work,
  // so use alternative:
 
  if (theAlpha < 1.0) 
  {
    do 
    {
      x = pow([self duuDouble], 1.0 / theAlpha);
      y = pow([self duuDouble], 1.0 / (1.0-theAlpha));
    } while (x+y > 1.0);

    x = x / (x+y);
    y = [self getExponentialWithMean: 1.0];

    return (x*y / theBeta);
  }

  // For alpha > 1.0, use rejection method:
   
  do 
  {
    do 
    {
      do 
      {
        v1 = 2.0 * [self duuDouble] - 1.0;
        v2 = 2.0 * [self duuDouble] - 1.0;
      } while (v1 * v1 + v2 * v2 > 1.0);
      y = v2 / v1;
      am = theAlpha - 1.0;
      s = sqrt(2.0 * am + 1.0);
      avg = s * y + am;
    } while (avg <= 0.0);
    e = (1.0 + y * y) * exp(am * log(avg / am) - s * y);
  } while ([self duuDouble] > e);

  return (avg / theBeta) ;

}

// ----- 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 = GAMMADISTMAGIC;
  // 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->theAlpha = theAlpha;
  internalState->theBeta  = theBeta;

  // 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     != GAMMADISTMAGIC)
       || (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;
    // duuMask.d = internalState->duuMask;
  // distribution specific data:
  theAlpha      = internalState->theAlpha;
  theBeta       = internalState->theBeta;

  // nothing is returned from a (void) function

}

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

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

  (void)sprintf(buffer," GammaDistribution printSelf: \n");
  (void)sprintf(buffer,"      theAlpha = %f\n", theAlpha);
  (void)sprintf(buffer,"       theBeta = %f\n", theBeta);

  [ super describe: outStream ];
  [outStream catC: buffer];
  [outStream catC: "\n"];
  //  return self;
}

-(int) verifySelf {

   return 1;
}

@end
