/*
Background:
: I'm looking for a formula that mutates a bit string (vector of 1's and 0's)
: given that the probability of mutating an individual bit is pmut.
: ...........
: I've seen references to a faster formula, that does not require generating a
: random number at each bit.  But I've never seen the formula itself or the
: methodology.

I use what I think is a relatively efficient method based on Poisson
distribution, which uses one call to a random number per chromosome
(this decides how many mutations there are) followed by one call to
a random number for each (if any) such mutations (this decides where
they are -- actually where there is a potential clash, there is an extra
call to random, but this is rare).

There is an additional once-only overhead of calculating various numbers
relating to a Poisson distribution --- if these calculations are
stored somewhere they do not have to be repeated.

Some C code, here using drand48() to produce random numbers in range [0,1].

Inman Harvey
inmanh@cogs.susx.ac.uk
*/

/************************************************************************/

#include <math.h>

/* Poisson distribution
 * poisson(x,m)
 * Inman Harvey Feb 1995
 * This gives probability that a longish genotype, with expected m
 * mutns per genotype, actually has x mutns exactly.
 * Taken from 'Mathematical Statistics' J.E. Freund, Prentice-Hall 1992
 * page 204. "Excellent approx for gen. length >=100, m<10"
 */
double poisson(int x, double m)
{
    int i;
    double lambda, result;
    lambda=m;
    result=pow(lambda,(double)x) * exp(-lambda);
    for (i=2;i<=x;i++)          /* divide by factorial(x) */
        result/=(double)i;
    return(result);
}

/* Find-poisson
 * find_poisson(n)
 * This assumes a global array of doubles
 * double pois[LEN+1]
 * This array will be initialised [ELSEWHERE**] to all -1.0, but then
 * when the probability that there are exactly n mutns
 * in a genotype of length LEN is needed (0<= n <= LEN)
 * it is calculated and stored here. Calculations below are based
 * on an expected MUT mutns, alter to suit.
 */
double find_poisson(int n)
{
    if (pois[n]<0.0) pois[n]=poisson(n,MUT); /*calculate first time only*/
    return(pois[n]);
}

/* Poisson mutation
 * Inman Harvey Feb 1995
 * First calc actual number of mutns, in a genotype length LEN
 * (Here assumed an array of ints, alter to suit)
 * Then do each at a different spot
 */
void poisson-mutate(int *genotype)
{
    int i,nm,pos;
    double r,tot=0.0;
    int mut_pos[LEN];       /* keep track of mutn posns*/

/* number of Poisson mutations nm
 * Spin the roulette wheel for r, 0.0<r<1.0
 * Then sum in turn the probability of 0 mutns, of 1 mutn,
 * of 2 mutns...until the sum exceeds r
 */
    r=drand48();
    for (nm=0;nm<LEN;nm++) {
        tot+=find_poisson(nm);
        if (tot>r) break;   /* exit with nm as number of mutns*/
    }
/* cheap exit if no actual mutations, only needed one call to drand48*/
    if (nm==0) return;

/* Now assume there were one or more actual mutations, nm>0
 * initialise list of positions where mutations are placed*/
    for (i=0;i<LEN;i++) mut_pos[i]=0;

/* set the nm positions for mutations at different positions*/
    for (i=0;i<nm;i++) {
        do (pos=(int) (drand48()*(double)LEN);}
            while (mut_pos[pos]==1);  /*avoid 2 at same place*/
        mut_pos[pos]=1;           /*keep track of mutn posns*/
        genotype[pos]=1-genotype[pos];    /* Substitute here your
                                        favourite mutation */
    }
}
/*******************************************************************/