/*-------------------------------------------------------------

Copyright (C) 2000 Peter Clote
All Rights Reserved.

Permission to use, copy, modify, and distribute this
software and its documentation for NON-COMMERCIAL purposes
and without fee is hereby granted provided that this
copyright notice appears in all copies.


THE AUTHOR AND PUBLISHER MAKE NO REPRESENTATIONS OR
WARRANTIES ABOUT THE SUITABILITY OF THE SOFTWARE, EITHER
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE, OR NON-INFRINGEMENT. THE AUTHORS
AND PUBLISHER SHALL NOT BE LIABLE FOR ANY DAMAGES SUFFERED
BY LICENSEE AS A RESULT OF USING, MODIFYING OR DISTRIBUTING
THIS SOFTWARE OR ITS DERIVATIVES.

-------------------------------------------------------------*/


/* Program: genetic.c  */


/* Sample program using the Population Genetics Algorithm */
/* Set includes {0,1,2,...,2^n-1} where n = NumBits */
/* Written by Rich Fiorito, modified by P. Clote, Feb 10, 1998 */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define MAX_STEPS 10000  /* Arbitrary upper bound for convergence */
#define NumBits 16
#define MAX ((unsigned) (1<<NumBits)-1) /* equals 2^NumBits  - 1 */
#define P 32		/* population size roughly log(MAX) */
#define cprob 0.25	/* crossover prob */
#define mprob 0.001	/* mutation prob  */
#define change 0.00001	/* termination condition */
#define F(x) (sin(pow(sqrt(x),1.5)) + 2)



/* Goal: find maximum of F in integer domain {0,...,2^NumBits -1} 
Must assume that F(x) always nonnegative. Here
       F(x) = -(x-1000)^2 + MAX^2  
Earlier, I tried F(x) = (double)(x) for debugging purposes. 
A subtle defect with this program is that for even simple functions,
the values pfitness and nfitness (population fitness values) overflow
and form nonsense values (i.e. negative values, though all individual
fitness values are positive. One should rather keep track of average
population fitness, and even in this case, one cannot just compute
pfitness/P. In any case, the population fitness (as a sum) plays
no role in the algorithm (except for termination). Here, we simply
wanted to display how the population on the average becomes fitter.

However on functions whose range is contained in small interval, as
in F(x) above, or as in the Unger-Moult application maximizing unit contacts,
there is no problem with overflow, and convergence (to a stable population)
is very rapid.  */

 
#define POW2(x) ((unsigned int) 1<<x)


void initialize (unsigned int []);
void crossover(unsigned int [], unsigned int [], double);
void popfitness (unsigned int [], double*, double*, unsigned int*);
void relativefitness (unsigned int [], double [], double);
void mutate (unsigned int []);
void refresh (unsigned int [], unsigned int []);

main()
{
  /* variable declarations */
  unsigned int pop[P], npop[P], best, bestOffspring;
    /* npop is new population, best is individual with largest F(x) */
  double pfitness, nfitness, cfitness;
  int i, step;

  /* initialize the population at random */
  srand((unsigned)time(NULL));
  initialize(pop);
  
  /* initialize population fitness nfitness, fitness of best individuum
     cfitness and index of this best individumm. */
  popfitness(pop,&nfitness,&cfitness,&best);

  /* the program loop */
  step = 0;
  do {
    step++;
    pfitness = nfitness;
    crossover(pop,npop,pfitness); 	/* perform crossover */
    mutate(npop); 			/* perform pointwise mutations */
    popfitness(npop,&nfitness,&cfitness,&bestOffspring);
    	/* determine fitness of new population */
    if(cfitness > F(best)) best = bestOffspring;
    refresh(pop,npop);
  }
  while ((step < MAX_STEPS) && (abs(pfitness-nfitness)>change)); 


  /* output */
  printf("Last population fitness values\n\n");
  printf("  INDIVIDUAL\t FITNESS\n\n");
  printf("-----------------------------------------------\n");
  for(i=0;i<P;i++) printf("%4d  %u \t%f\n", i,pop[i],F(pop[i]));
  printf("\n\n");
  printf("Average fitness in last population: %f\n", nfitness/P);
  printf("Number of iteration steps: %d\n", step);
  printf("Difference in pop fitness: %f\n", abs(pfitness-nfitness));
  printf("F achieves max of %f at %u\n", F(best), best);
}

/* initialize population pop of size P with random integers
*/
void initialize(unsigned int pop[P])
{
  int i;
  for(i=0; i<P; i++)
    pop[i] = rand()%MAX;	
	/* pop consists of random integers in domain {0,...,2^NumBits -1} */
}

/* perform crossovers
   generates new population D by crossovers from A, gets the population 
   fitness of pop A <pfitness>.

   first, the individuals are selected with probs according to their
   fitness. Therefore, a roulette wheele method is applied. the
   selected individuals are copied to population B
   
   second, choose individuals to crossover with crossover prob <cprob>
   and copy them to the lower part of population C. Individuals, which
   do not crossover are copied to the upper part of C.

   third, individuals of the lower part of C are crossed over. Result
   is copied to D.
*/
void crossover(unsigned int A[P], unsigned int D[P], double pfitness)
{
  static double p[P]; /* p[i] is relative fitness of i-th individual,
			 i.e. the fitness of i divided by the population
			 fitness */
  static unsigned int B[P], C[P];
  int i,j,x,y,m,M,r,co,index;
  unsigned int tempFirst, tempSec, temp;
	/* tempFirst, tempSec used to construct first, second child */
	/* co is crossover point */
  double q[P]; 		/* cumulative probabilities */
  double total, z;

  relativefitness(A,p,pfitness); /* calculate relative fitness of the
				    individuals */

  /* determine cumulative probabilities 
     q[i] = sum_{j<=i} p_j
  */
  q[0]=p[0];
  for (i=1;i<P;i++) q[i] = q[i-1] + p[i]; 
 
  /* amplify fit individuals placing into B, using ROULETTE Wheel */ 
  for (i=0;i<P;i++) {
    z = (double) rand()/ (RAND_MAX+1.0);  /* random z in (0,1( */
    j=0;
    while (q[j] < z) j++;	/* find j s.t. q[j-1] < z <= q[j] */ 
    B[i] = A[j];
  }
 
  /* using crossover probability cprob, place certain individuals
     from B into C, place a "bar" m above those individuals copied into
     C and then copy the remaining individuals not chosen for crossover
     from B into the "top" part of C. The individuals chosen for crossover
     thus appear in the "bottom" part of C, below index m. */ 
  m=0; M=P-1;		/* bottom and top indices of C */
  for(i=0; i<P; i++) {
    z = (double) rand()/(RAND_MAX+1.0);
    if(z<cprob) {
      C[m] = B[i];
      m++;
    } else {
      C[M] = B[i];
      M--;
    }
  }

  /* need an even number of individuals for crossover */
  if(m%2 != 0) { m++; M=m-1; }


  /* perform a "shuffle" of the bottom part of C, to ensure
     that those chosen for crossover are randomly paired. Let
     index = m, the "barrier", choose random x from 0...index-1,
     then swap C[x] and C[index-1], and decrement index. When
     index = 2, there are only 2 remaining individuals C[0],C[1],
     so these must be paired.  */
  for(index=m; index>2; index--) {
    x = rand()%index;
    temp = C[index-1];
    C[index-1] = C[x];
    C[x] = temp;
  }

  /* now perform crossover, using bit operations on the integer individuals */
  for (y=0;y<m;y=y+2) {
    co = rand()%NumBits; /* site of crossover */
	/* this should be fixed in machine indep manner to handle
           larger sizes */
    tempFirst = tempSec = 0;
    r = POW2(NumBits-1);
    /* C[y] is Mom, C[y+1] is Dad */
    /* Put left side of Mom with right side of Dad */
    for (j=NumBits-1;j>=co;j--){
      tempFirst *= 2; tempSec *= 2;
      if ((r & C[y]) > 0) tempFirst++;
      if ((r & C[y+1]) > 0) tempSec++;
      r = r/2;
    }
    /* Put left side of Dad with right side of Mom */
    for (j=co-1;j>=0;j--) {
      tempFirst *= 2; tempSec *= 2;
      if ((r & C[y+1]) > 0) tempFirst++;
      if ((r & C[y]) > 0) tempSec++;
      r = r/2;
    } 
    D[y] = tempFirst; 
    D[y+1] = tempSec;
  }
  for (y=m;y<P;y++) D[y]=C[y]; /* fill in remaining individuals, which
				  were not crossed over */
}



/* popfitness outputs (via pointers) the fitness of the overall
population, the fitness of the fittest individual of the population,
and the index of the fittest individual as an unsigned integer in F's
domain */

void popfitness(unsigned int pop[P], double *ptrToPopFit,
     double *ptrToIndFit, unsigned int *ptrToBest)
{
  double f, sum, maxFitness;
  int i;
  int bestindex;
  bestindex = 0;      /* initially assume F(pop[0]) maximal */
  sum = maxFitness = F(pop[0]);
  for(i=1;i<P;i++) {
    f = F(pop[i]);
    sum += f;
    if ( f > maxFitness ) {
      bestindex = i;
      maxFitness = f;
    }
  }
  *ptrToPopFit = sum;
  *ptrToIndFit = maxFitness;
  *ptrToBest = pop[bestindex];
}

/* relativefitness(pop,p,pfitness) computes the relative
   fitness(fitness divided by pfitness) of each individual in pop and
   returns the results in array p.  */
void relativefitness(unsigned int pop[P], double p[P], double pfitness)
{
  int i;
  for(i=0;i<P;i++) p[i] = F(pop[i])/pfitness;
}

void mutate(unsigned int pop[P])
{
  int i,j, bit;
  double rreal;
  unsigned int temp, r;
  for(i=0; i<P; i++) {
    temp = 0;
    r=POW2(NumBits-1);
    for (j=0; j<NumBits; j++) {
      temp *=2;
      bit = pop[i]&r;
      rreal = (double) rand()/ (RAND_MAX+1.0);
      if (rreal < mprob) {
        if (bit==0) temp++;
      }
      else
        if (bit>0) temp++;
      r = r/2; 
      /* toggle bit with probability mprob */
    }
    pop[i]=temp;
  }
}

/* refresh copies population npop to pop
 */
void refresh(unsigned int pop[MAX], unsigned int npop[MAX])
{
  int i;
  for(i=0; i<P; i++)
    pop[i] = npop[i];
}