/*
   This is code for a stady state genetic programming system.  It is a demo
of the rtree package as well as studying a problem in simple functional 
morphing.

*/
#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
#include <math.h>

#include "rtree.h"

//global defines and variables
/* population size */
#define popsize 100
/* number of runs to do */
#define runs 20
/* number of mating events per run */
#define mevs 20000
/* initial and maximum tree size */
#define its 12
#define mts 36
/* file name for results */
#define renm "results.dat"

node *pop[popsize];     //the population

double fit[popsize];    //the fitness values (sum squared error to be minimized

node *best;             //tracking variable for best so far
double bestfit;         //and fr best value so far

void initialize();       //set all the rtree controlls to appropriate values
void initpop();          //create a random population
void delpop();           //delete a population
double findfit(node *t); //compute the fitness of a parse tree
void mate();

main(){

int mev,run;  //loop indices
fstream aus;  //output file

   aus.open(renm,ios::out);
   initialize();
   for(run=0;run<runs;run++){//do the requested number of runs
      initpop();
      for(mev=0;mev<mevs;mev++){//do the requested number of mating events
         mate();
         if(mev%1000==0)cout << bestfit << endl;
      }
      aus << bestfit << " " << best << endl;
      aus << "----------" << endl;
      delpop();
   }
   aus.close();   
}

void initialize(){ //set all the rtree controlls to appropriate values

   srand48(121461); //seed the random number generator

   TurnOff(DelOP);  //plain parse trees, no delays
   TurnOff(UA1OP);  //
   TurnOff(UA2OP);  //or
   TurnOff(BA1OP);  //
   TurnOff(BA2OP);  //ADFs

   SetNumberOfVariables(2);

}

void initpop(){          //create a random population

int i;

   best=(node *)0;  //
   bestfit=-1;
   for(i=0;i<popsize;i++){
      pop[i]=RandTree(its);
      fit[i]=findfit(pop[i]);
      if((bestfit<0)||(fit[i]<bestfit)){
         KillTree(best);
         best=CopyTree(pop[i]);
         bestfit=fit[i];
      }
   }
}

void delpop(){           //delete a population

int i;

   for(i=0;i<popsize;i++)KillTree(pop[i]);

}

double findfit(node *t){

double i,q,accu,v[2];

   accu=0;
   for(v[0]=-1;v[0]<1;v[0]+=0.05){//over the set of sample points
      v[1]=0;                     //prepare the input vector
      q=Eval(t,v)-1/(i*i+1);      //evaluate and compare to 1/(x^2+1)
      accu+=q*q;                  //sum squared error into accumulator
      v[1]=1;                     //reprepare the input vector
      q=Eval(t,v)-1/(i*i-2*i+2);  //evaluate and compare to 1/((x-1)^2+1)
      accu+=q*q;                  //sum squared error into accumulator
   }
   return(accu);  //fitness is sum of squared error
}

void mate(){

int w[4];  //tournament for mating
int sw;

   //generate four distinct population members
   w[0]=lrand48()%popsize;
   do {
      w[1]=lrand48()%popsize;
   } while(w[0]==w[1]);
   do {
      w[2]=lrand48()%popsize;
   } while((w[0]==w[2])||(w[1]==w[2]));
   do {
      w[3]=lrand48()%popsize;
   } while((w[0]==w[3])||(w[1]==w[3])||(w[2]==w[3]));

   //Hard sort so the largest (worst) fitness is in pop[w[0]], pop[w[1]] 
   if(fit[w[0]]<fit[w[3]]){sw=w[0];w[0]=w[3];w[3]=sw;}
   if(fit[w[1]]<fit[w[2]]){sw=w[1];w[1]=w[2];w[2]=sw;}
   if(fit[w[0]]<fit[w[1]]){sw=w[0];w[0]=w[1];w[1]=sw;}
   if(fit[w[2]]<fit[w[3]]){sw=w[2];w[2]=w[3];w[3]=sw;}
   if(fit[w[1]]<fit[w[2]]){sw=w[1];w[1]=w[2];w[2]=sw;}

   //deallocate worst two, copy best two over them, perform crossover
   KillTree(pop[w[0]]);
   KillTree(pop[w[1]]);
   pop[w[0]]=CopyTree(pop[w[2]]);
   pop[w[1]]=CopyTree(pop[w[3]]);
   cross(pop[w[0]],pop[w[1]]);
   chop(pop[w[0]],mts);  //ensure tree no larger than mts
   chop(pop[w[1]],mts);

   //now give each child one mutation - one subtree one node
   mutateST(pop[w[0]]);
   mutateND(pop[w[1]]);

   //find fitness and update tracking variables
   fit[w[0]]=findfit(pop[w[0]]);
   if(fit[w[0]]<bestfit){
      KillTree(best);
      best=CopyTree(pop[w[0]]);
      bestfit=fit[w[0]];
   }
   if(fit[w[1]]<bestfit){
      KillTree(best);
      best=CopyTree(pop[w[1]]);
      bestfit=fit[w[1]];
   }
}