/* Mutation routines.

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

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


static int mean = 0;

double int_poissons[INTRANGE];
double frac_poissons[FRACRANGE];



/*This initializes the poisson values so they don't have to be recalculated each
  time.*/
  
void setup_poissons(int lambda, double mut_rate)
{
	int i;
	double lambda_term, p;
	
	if (lambda <= 0) printf("Lambda cannot be <= 0 : lambda = %d\n", lambda);
	else 
	{
		mean = lambda;
		lambda_term = exp(-lambda);
		
		for(i=1, p=lambda_term, int_poissons[0]=p; i < INTRANGE; i++)
		{		
			  if (i >= 2 * lambda) int_poissons[i] = 1.0; /*prevents skew*/
			  else
			  {
		      	if (p <= MIN_P) p = 0.0;
		      	else p = p * ((double) lambda / i);
		      	int_poissons[i] = p + int_poissons[i-1];
		      }
		}
	}

	if (mut_rate <= 0) printf("Mutation rate cannot be <= 0 : mut_rate = %f\n", mut_rate);
	else 
	{
		lambda_term = exp(-1 * mut_rate);
		
		for(i=1, p=lambda_term, frac_poissons[0]=p; i < FRACRANGE; i++)
		{		
		      if (p <= MIN_P) p = 0.0;
		      else p = p * (mut_rate / i);
		      frac_poissons[i] = p + frac_poissons[i-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.  And returns a value
  with a mean of 0.  Note that is only is designed to work with integer lambdas.*/

int get_int_poisson()
{
	double rand;
	int front, back;
	
	rand = knuth_random();
	
	if (rand >= int_poissons[INTRANGE-1]) return INTRANGE-1;
	
	for(front=-1, back=2*mean; front+1 < back;)
	{
		if (rand < int_poissons[(front + back) / 2])
			back = (front + back) / 2;
		else front = (front + back) / 2;
	}
	
	return back - mean;
}

/*This works with fractional lambdas*/
int get_frac_poisson()
{
	double rand;
	int front, back;
	
	rand = knuth_random();
	
	if (rand >= frac_poissons[FRACRANGE-1]) return FRACRANGE - 1;
	
	for(front=-1, back=FRACRANGE-1; front+1 < back;)
	{
		if (rand < frac_poissons[(front + back) / 2])
			back = (front + back) / 2;
		else front = (front + back) / 2;
	}
	
	return back;
}	
	

/********** 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 > MIN_P)*/)
      {
      i++;
      p = p * lambda / i;
      sum += p;
      }  /* while */



   return i;
}  /* poisson */