/* Mutation routines.

   09/20/93 AW  Created. (from Melanie Mitchell and Stephanie Forrest's lisp code)
   
   I adapted this to my Lewontin model, 07/25/94 CCM.
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "random.h"
#include "poisson.h"




#ifdef NOTDEFINED
static double **poissons;

/*This stuff was crashing after the 4th calculation of the first lambda.*/

/*This initializes the poisson values so they don't have to be recalculated each
  time.*/
  
void setup_poissons()
{
	int i, j;
	double lambda_term, p;
	
	poissons = (double **) calloc(MAXLAMBDA*MAXRANGE, sizeof(double));
	
	printf("Okay, I've allocated the space.\n");
	
	if (poissons == NULL) printf("Ran out of memory.\n");
	else
	{
	for(i=1; i < MAXLAMBDA; i++)
	{
		printf("\nIterating i = %d\nJ: ", i);
		lambda_term = exp(-i);
		p = lambda_term;
		printf("here we go...");
		poissons[i][0] = p;
		printf(" there we went with (%f).\n", poissons[i][0]);
		
		for(j=1; j < MAXRANGE; j++)
		{
			printf("*** %d ", j);
		      if (p <= 0.00005) p = 0.0;
		      else p = p * ((double) i / j);
		      printf("-");
			fflush(stdout);
		      poissons[i][j] = p + poissons[i][j-1];
		}
	}
	}

}

  /*This is a faster version of "poisson" below.  Instead of doing any calculating,
  it does a binary search on an array of summed poisson values.*/

int get_poisson(int lambda)
{
	double rand;
	int front, back;
	
	rand = knuth_random();
	
	for(front=-1, back=MAXRANGE-1; front+1 < back;)
	{
		if (rand < poissons[lambda][(front + back) / 2])
			back = (front + back) / 2;
		else front = (front + back) / 2;
	}
	
	return back;
}
	
#endif  /*NOTDEFINED*/

/********** poisson **********/
/* parameters:	lambda
   called by:	
   actions:	given a mutation probability (lambda), 
   			returns a random number somewhere around lambda
   			using the Poisson distribution.
		This algorithm was copied directly from the Lisp
		royal road code of Mitchell and Forrest.
*/
int poisson(double lambda)
{
   double lambda_term;
   double p;
   double sum;
   double unif_rand;
   int i;



   if (lambda < 0.0)
      {
      printf(" Error(poisson): bad (neg.) lambda value.\n");
      return -1;
      }  /* if */

   unif_rand = knuth_random();
   lambda_term = exp(-lambda);
   sum = p = lambda_term;


   i = 0;
   while ((sum <= unif_rand) && i < MAXRANGE /*(p > 0.0000005)*/)
      {
      i++;
      p = p * lambda / i;
      sum += p;
      }  /* while */



   return i;
}  /* poisson */