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

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.

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


/*-----------------------------------
   Genetic
   P. Clote, Feb 1999

This class implements a genetic algorithm, in order
to determine the maximum of a function
   f : {0,1,2,...,2^n-1}  ---> Reals
where n = NumBits, a constant. 
-------------------------------------*/



/*-----------------------------------------------
This class is used as a "record" for values population 
fitness, index in array of individuals to the fittest 
individual, and the fitness value of the fittest individual.
------------------------------------------------*/
class Fitness {
   double PopFit;
   double IndFit;
   int Best;
}



class Genetic {
   static final int MAX_STEPS = 10000;  
             // Arbitrary upper bound for convergence 
   static final int NumBits = 7;
   static final int MAX = (int) Math.pow(2,NumBits);
             // equals 2^NumBits 
    static final int P = 2*NumBits;
             // population size is 2*log(MAX)
   static final double cprob = 0.25;
             // crossover prob 
   static final double mprob = 0.001;
             // mutation prob 
   static final double change = 0.00001;
             // termination condition 


   /*--------------------------------
     F is the function, whose maximum we want
     to compute. This function can be changed.
   --------------------------------*/
   static double F(int x) { 
       return (Math.sin(Math.pow(Math.sqrt(x),1.5)) + 2); // (-Math.pow(x,4)*Math.sin(x)+10*Math.pow(x,3)+20);
   }   

   static void printFunctionValues() {
      System.out.println("\n\nFunction Values \n");
      for (int i=0;i<MAX;i++)
         System.out.println(F(i));
   }



   public static void main ( String argv[] ) {

   
     /*-------------  variable declarations -----------*/
     int pop[] = new int[P];
     int npop[] = new int[P];
     int best;
       // npop is new population, 
       // best is individual with largest F(x) 
     double pfitness, nfitness;
       // pfitness is the previous generation's pop fitness
       // pfitness is the next generation's population fitness
 
     int i, step;
     double p[] = new double[P];	
       // p[i] is relative fitness of i-th individual, 
       // i.e. the fitness of i divided by the 
       // total population fitness  
     boolean stop;
       // condition to terminate main loop
   
     Fitness fitness = new Fitness();
       // define and create space for object fitness





     
     /*--------------------------------------
     The method main uses the following methods,
     defined later:
      
       (1) static void initialize(int pop[]) 
       (2) static int POW2(int x) 
       (3) static void crossover(int A[], int D[], 
                double p[], double pfitness) 
       (4) static Fitness popfitness(int pop[]) 
       (5) static void relativefitness(int pop[], 
                  double p[], double pfitness) 
       (6) static void mutate(int pop[]) 
       (7) static void refresh(int pop[], int npop[]) 
  
     --------------------------------------*/
  



     
   
     /*-------- initialize the population --------*/
     initialize(pop);
     fitness = popfitness(pop);
         // fitness is the initial population fitness
     relativefitness(pop,p,fitness.PopFit); 
       // calculate relative fitness array p[]
   
   
   
     /*-------- the main program loop -----------------*/
     step = 0;
     best = fitness.Best;
          // best is best fitness value of an individual
          // in the initial population 
     do {
       step++;
       pfitness = fitness.PopFit;
          // pfitness is "previous" population fitness
          // before crossover
       crossover(pop,npop,p,pfitness); 	
          // perform crossover 
       mutate(npop); 	
          //  perform pointwise mutations 
       fitness = popfitness(npop);
          //  determine fitness of new population 
       nfitness = fitness.PopFit;
          // pfitness is new population fitness
       if(fitness.IndFit > F(best)) 
          best = fitness.Best;
             // best is best fitness value of an individual
             // in the current population 
       refresh(pop,npop); 
       stop = ((step >= MAX_STEPS) ||
           (Math.abs(pfitness-nfitness) < change)); }

     while ( ! stop );
  

 
     /*-------------- output ------------*/
     System.out.println("Last population fitness values\n");
     System.out.println("  INDIVIDUAL\t FITNESS\n");
     System.out.println("--------------------" 
          + "---------------------------");
     for(i=0;i<P;i++) 
        System.out.println(i + "\t" + pop[i] 
          + "\t" + F(pop[i]));
     System.out.println("\n\n");
     System.out.println("Average fitness in last population: " 
          + nfitness/P);
     System.out.println("Number of iteration steps: " 
          + step);
     System.out.println("Difference in pop fitness: " 
          + Math.abs(pfitness-nfitness));
          // absolute value of difference
     System.out.println("F achieves max of "
          + F(best) + " at " + best);


     printFunctionValues(); 
          // can plot F(x) using gnuplot
   } // end of main
   

    /*-------- initialize the population -------*/
    static void initialize(int pop[]) {
      for(int i=0; i<P; i++)
        pop[i] = (int) (MAX * Math.random());
        // pop now consists of random integers 
        // in the domain {0,...,2^NumBits -1} 
    }
  
  
         
  
    static int POW2(int x) {
       return ((int) Math.pow(2,x));
    } // Math.pow returns a double, 
      // where we require unsigned integerr
   
  
  
 
   static void crossover(int A[], int D[], 
              double p[], double pfitness) {
     int B[] = new int[P];
     int C[] = new int[P];
     int i,j,x,y,m,M,r,co,index;
     int tempFirst, tempSec, temp;
       // tempFirst, tempSec used to construct 
       // the first and second child 
       // co is crossover point 
   
     double q[] = new double[P]; 		
       // cumulative propabilities where
       //  q[i] = sum_{j<=i} p_j 
   
     double total, z;
   
     relativefitness(A,p,pfitness); 
       // calculate relative fitness of the individuals 
   
   
   
     /*---- determine cumulative probabilities ---*/
     q[0]=p[0];
     for (i=1;i<P;i++) q[i] = p[i]+q[i-1]; 
   
   
    
     /*----------------------------------------------
       amplify fit individuals placing into B, by
       using ROULETTE Wheel 
     ----------------------------------------------*/
     for (i=0;i<P;i++) {
       z = Math.random();
           // random z in (0,1) 
       j=0;
       while (q[j] < z) j++;
           // find j such that  q[j-1] < z <= q[j] 
       B[i] = A[j];
       } 
   
   
   
     /*-----------------------------------------
     REMARK: 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 = Math.random();
       if(z<cprob) {
         C[m] = B[i];
         m++;}
       else {
         C[M] = B[i];
         M--;}}
   
   
   
   
     /*----------------------------------------
      The following is to ensure that there are
      an even number of individuals for crossover 
     ----------------------------------------*/
     if(m%2 != 0) { m++; M=m-1;}
   
   
   
   
     /*----------------------------------------------
     REMARK: 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 = (int) (Math.random() * index);
       temp = C[index-1];
       C[index-1] = C[x];
       C[x] = temp;}
   
   
   
   
   
     /*----------------------------------------------
     REMARK: Now perform crossover, using bit 
     operations on the integer individuals 
     ----------------------------------------------*/
     for (y=0;y<m;y=y+2) {
       co = (int) (Math.random() * NumBits);  
          // choose crossover point
       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, 
        // those which were not crossed over 
   }
   
   
   /*----------------------------------------------
   Method "popfitness" outputs the fitness of the overall
   population, the fitness of the fittest individual of 
   the population, and the index of the fittest individual
   -----------------------------------------------*/
   static Fitness popfitness(int pop[]) {
     Fitness aux = new Fitness();
       // auxilliary variable for returning popfitness
     double f, sum, maxFitness;  
     int bestindex;
     bestindex = 0;      
        // initially assume F(pop[0]) maximal 
     sum = maxFitness = F(pop[0]);
     for(int i=1;i<P;i++) {
       f = F(pop[i]);
       sum += f;
       if ( f > maxFitness ) {
         bestindex = i;
         maxFitness = f;} 
     }
     aux.PopFit = sum;
     aux.IndFit = maxFitness;
     aux.Best = pop[bestindex];
     return(aux);
   }
   
   
   
   static void relativefitness(int pop[], 
                double p[], double pfitness) {
     for(int i=0;i<P;i++) 
          p[i] = F(pop[i])/pfitness;
   }
   
   
   
   
   static void mutate(int pop[]) {
     int i,j, bit;
     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;
         if (Math.random() < mprob) {
           if (bit==0) temp++; }
         else
           if (bit>0) temp++;
         r = r/2; 
         // toggle bit with probability mprob 
         }
        pop[i]=temp; 
      }
   }
   
   static void refresh(int pop[], int npop[]) {
     for(int i=0; i<P; i++) pop[i] = npop[i];
   }
   
 } // end of class   
  
 


/*----------------------------------------------------
REMARK:
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  
For debugging purposes, I tried F(x) = x.
A subtle defect with this type of program in C is that for even 
simple functions, the population fitness values can overflow,
and form nonsense values (i.e. negative values, 
though all individual fitness values are positive. 
In Java I have not seen this problem, but one should
check this.


SAMPLE RUN
==========

Last population fitness values

  INDIVIDUAL	 FITNESS

-----------------------------------------------
0	111	1.4492132720584E8
1	111	1.4492132720584E8
2	119	9.133073073686197E7
3	111	1.4492132720584E8
4	119	9.133073073686197E7
5	111	1.4492132720584E8
6	119	9.133073073686197E7
7	111	1.4492132720584E8
8	111	1.4492132720584E8
9	111	1.4492132720584E8
10	111	1.4492132720584E8
11	111	1.4492132720584E8
12	111	1.4492132720584E8
13	111	1.4492132720584E8



Average fitness in last population: 1.334376279624876E8
Number of iteration steps: 7
Difference in pop fitness: 4.76837158203125E-7
F achieves max of 1.4492132720584E8 at 111


Function Values 

20.0
29.1585290151921
85.45124117078909
278.56927934715077
853.7414387988297
1869.3276716644614
2542.1224856658077
1872.5751764761874
1087.5886218306282
4606.090618328835
15460.211108893698
27970.856614108852
28426.376027657017
9989.60926119431
-10595.172176374137
849.1780920459641
.
.
.
-2.3254064966858208E8
----------------------------------------------------*/