/***************************************************************************
 *  File name : longruns.c (modifed from tourn.c)
 * 10000 runs, 300 pop_size
 *  First written: 30 July 1997 by Chilukuri K. Mohan at Syracuse University 
 ***************************************************************************/    
#include <stdio.h>
#include <malloc.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <sys/time.h>
struct    timeval *tp;
struct    timezone *tzp;


int *listofprimes; /*  primes >= 13 */
int max_prime_count =  40;
time_t t=0.0;
int max_trials = 10;
int maxpos =0;
float mut_rate = 0.033333333;
long int jran = 2147483647;
float	epsilon   = 0.01;
float   pop_mean_fitness;
float   pop_max_fitness;

int	max_num_of_generation = 10000;
int	indiv_length = 300;
int sizeby3 = 100;
int	print_freq = 250;
int	pop_size   = 300;
int     sel_choice = 1;
int     op_choice = 0;

float **gradient;
float *onegrad;
float    *sum_maxf;
int	**population;
int	**offspring_set;
float	*newfit;
float   *oldfit;
int     no_of_mutation = 0;
/*********** function prototypes ***********************************/

double when();
void initialization();
float mean_pop_fitness();
float max_pop_fitness();
void init_population();

void execute_one_generation();

int fitness_proportionate();
void change_contents();
void d_tournament();
void nd_tournament();
void tournament();
void replacement();
void elitist();

float fitness_modified();
void crossover();
void mutation();
int  ux();
void one_ptx();
void two_ptx();
int  pux();
void  hux();
float weight();
float variation();
float fitness();
float max();
float min();
float triple();
float ugly3();
float easy3();

float myrand();

float getaprime();

extern struct timeval;
extern struct timezone;


double when(){           /* time function to record current time*/
 struct timeval tp;
 struct timezone tzp;
 gettimeofday(&tp,&tzp);
 return ((double) tp.tv_sec + (double) tp.tv_usec*1e-6);
}

/***********************************************************************
 * main
 **********************************************************************/
main ()
{
  int  argc,	i,j,k,	no_of_prints, trial_num,	no_gen=0;
  double     start_time = 0.0;

  initialization();   
  no_of_prints = max_num_of_generation/print_freq;
  for (i=0;i < no_of_prints; i++)     sum_maxf[i] = 0.0;
     start_time = when();
     for (trial_num=0; trial_num < max_trials; trial_num++)
     {
       init_population(); 
       for (i=0; i < no_of_prints; i++) {
	 for (j=0; j<print_freq; j++) 		execute_one_generation();
	 sum_maxf[i] = 	 sum_maxf[i] + pop_max_fitness;
       }
     }
     gettimeofday(tp,tzp);
     no_gen = print_freq;
     for (i=0; i < no_of_prints; i++) {
       if (no_gen<100) 	 fprintf(stderr,"  ");
	 else 	 if (no_gen<1000) 	 fprintf(stderr," ");
       fprintf(stderr,"%d gens., max. fitness %7.2f \n",
	       no_gen,sum_maxf[i]/max_trials);
       no_gen = no_gen + print_freq;
   }
      free	(population);
      free	(offspring_set);            
fprintf(stderr,"\n Average Elapsed Time = %f \n ", (when() - start_time)/max_trials);
}


/**********************************************************************/
void initialization()
{  int i, std_parameters;

  sel_choice =0;

  fprintf(stderr,"Operator_choice (0 newton) (1 1PTX) (2 2PTX) (3 UX):");
  scanf("%d",&op_choice);

  fprintf(stderr, "Standard parameters (0 for no, 1 for yes)? ");
  scanf("%d", &std_parameters);
  if ( std_parameters) 
    {if (!(op_choice))
       {max_num_of_generation= max_num_of_generation/10;
        print_freq = max_num_of_generation/25;
      }}
    else{
      fprintf(stderr,"Number of trial runs: ");
      scanf("%d",&max_trials);

      fprintf(stderr,"Number of generations: ");
      scanf("%d",&max_num_of_generation);

      fprintf(stderr,"Length of individuals: ");
      scanf("%d",&indiv_length);
      sizeby3 = indiv_length/3;
      mut_rate = 1.0/indiv_length;

      fprintf(stderr,
	  "Frequency (integer >=1) with which changes should be printed:");
      scanf("%d",&print_freq);

      fprintf(stderr,"Population size :");
      scanf("%d",&pop_size);

    }

  population     = (int **) calloc(pop_size,sizeof(int));
  offspring_set  = (int **) calloc(pop_size,sizeof(int));
  sum_maxf  = (float *) calloc(500,sizeof(float));
  newfit = (float *)  calloc(pop_size,sizeof(float));
  oldfit = (float *)  calloc(pop_size,sizeof(float));
  gradient = (float **)  calloc(pop_size,sizeof(float));
  onegrad  = (float *)  calloc(pop_size,sizeof(float));
  listofprimes =  (int *) calloc(max_prime_count,sizeof(int));

   listofprimes[0]=13;listofprimes[1]=17   ;listofprimes[2]=19   ;listofprimes[3]=23   ;listofprimes[4]=29   ;
listofprimes[5]=31   ;listofprimes[6]=  37 ;listofprimes[7]=41   ;listofprimes[8]=43   ;listofprimes[9]=47   ;
listofprimes[10]=53   ;listofprimes[11]=59   ;listofprimes[12]=61   ;listofprimes[13]=67   ;listofprimes[14]=71   ;
listofprimes[15]=79   ;listofprimes[16]=83   ;listofprimes[17]=89   ;listofprimes[18]=97   ;listofprimes[19]=101   ;
listofprimes[20]=103   ;listofprimes[21]=107   ;listofprimes[22]=109   ;listofprimes[23]=113   ;listofprimes[24]=119   ;
listofprimes[25]=127   ;listofprimes[26]=131   ;listofprimes[27]=133   ;listofprimes[28]=137   ;listofprimes[29]=139   ;
listofprimes[30]=149   ;listofprimes[31]=151   ;listofprimes[32]=157   ;listofprimes[33]=163   ;listofprimes[34]=167   ;
listofprimes[35]=173   ;listofprimes[36]=179   ;listofprimes[37]=181   ;listofprimes[38]=191   ;listofprimes[39]=193   ;
 }
/****************************************** function definitions */

/* * finding mean fitness value of current population  */
float mean_pop_fitness()
{ int   i;  float result = 0.0;
  for(i=0;i<pop_size;i++) result += oldfit[i];
  return(result/ (float) pop_size);
}

/* * finding max fitness value w.r.t the population; 
   "maxpos" keeps track of the position of the current best
 */
float max_pop_fitness()
{ int i;  float cur_fit, result = -99999999.0;
  for(i=0;i<pop_size;i++) {
    cur_fit = oldfit[i];
    if (cur_fit > result) {maxpos = i;	result = cur_fit;}
  }
  return(result);
}

/***********************************************************************
 * init_population
 **********************************************************************/
void init_population()
{
    int  i,j, *indv; float  *onegrad;

    for (i=0; i < pop_size;i++) {
	indv = (int *) calloc(indiv_length,sizeof(int));
	onegrad = (float *) calloc(indiv_length,sizeof(float));
		  /* allocate memory for each individual */
	if ((indv == NULL)||(onegrad == NULL)) {
	   fprintf(stderr,"insufficient memory error\n");
	   exit(1);
	}
	for(j=0;j<indiv_length;j++) 
	  indv[j] = (myrand() < 0.5); /* =  0 or 1 */
	population[i] = (offspring_set[i] = indv);
	oldfit[i] = (newfit[i] = fitness(indv));
	for(j=0;j<indiv_length;j++) {
	  indv[j] = 1- indv[j];
	  onegrad[j] = fitness(indv) - oldfit[i];
	  indv[j] = 1- indv[j];
	}
	gradient[i] = onegrad;
    }
  }
/***********************************************************************
 * execute_one_generation
 **********************************************************************/
void execute_one_generation()
{
  int i,j;
  int *sample;

  sample = (int *) calloc(indiv_length,sizeof(int));
  if (sample == NULL) {
     fprintf(stderr,"No memory\n");exit(1);
  }

  pop_mean_fitness = mean_pop_fitness();  /* global variables */
  pop_max_fitness  = max_pop_fitness();
  for (i=0;i<pop_size;i++) {
      crossover(fitness_proportionate(),  /* was find_parent() */
		fitness_proportionate(),
		sample);
      mutation(sample);
      for(j=0;j<indiv_length;j++)	offspring_set[i][j] = sample[j];
      newfit[i] = fitness(sample);
  }
  free(sample);
  replacement();
}
/***********************************************************************/
int fitness_proportionate() /*reproduction*/
{
  float current, randval; int i=0;
  randval =  myrand() * pop_size *  pop_mean_fitness; 
           /* 0 <= randval < sum of fitness of all population members */
  current =  oldfit[0];
  while (randval > current)      { i++; current += oldfit[i];}
/*fprintf(stderr,"i=%d, randint=%d, current = %d,\n", i, randint, current);*/
  return((i<pop_size) ? i : pop_size-1); /*   return (i); */
}
/***********************************************************************/
void change_contents(index) /* UNLESS parent is the best one from prev.it. */
int index;
{ int j;
  if ((index != maxpos)      || (newfit[index] >    pop_max_fitness))
    {for (j=0;j<indiv_length;j++) 
       {population[index][j] = offspring_set[index][j];
	gradient[index][j] = -999999;
      }
     oldfit[index] = newfit[index];
   }
}
/***********************************************************************/
void d_tournament(index) /* DETERMINISTIC tournament */
int index;
{tournament(index,1);}

void nd_tournament(index) /* NON_DETERMINISTIC tournament */
int index;
{tournament(index,0);}

void tournament(index,det) /* det = 0 or 1 */
int index,det;
{int i,j,k,curbest; float curfit;
 for (i=0; i<pop_size; i++)
   { curbest = pop_size; curfit = newfit[i];
     for (j=0; j<index; j++)
       { k = (int) (myrand() * indiv_length);
	 if ((det && (oldfit[k] > curfit)) || 
	     (!(det) &&  (myrand() > curfit/(curfit + oldfit[k]))))
	   { curfit = oldfit[k]; curbest = k;}
       }
     if (curbest < pop_size)
       { for (j=0; j < indiv_length; j++)
	   offspring_set[i][j] = population[curbest][j];
	 newfit[i] = oldfit[curbest];
       } /* Later, also store gradient for reuse */
   }
 for (i=0; i<pop_size; i++) change_contents(i);
}
/***********************************************************************/
void replacement() /*  For changing contents of population */
{ int i;
  switch(sel_choice) {
    case 1: elitist(); break;
    case 2: tournament(0); break; /* non-deterministic tournament */
    case 3: tournament(1); break; /* deterministic tournament */
    default: for (i=0;i<pop_size;i++)  change_contents(i);  /* wholesale */
    }
}
/***********************************************************************/
void elitist() /* preserving better of the ith offspring and parents */
{ int i;
  for (i=0;i<pop_size;i++) 
      if (fitness(population[i]) < fitness(offspring_set[i]))
	  change_contents(i);
}
/***********************************************************************/
float fitness_modified(X, i) /*fitness when ith bit of pop[X] is changed*/
int X,  i;
{ float fXi;
  if (gradient[X][i] < -99) 
    { population[X][i] = 1-  population[X][i];
      fXi = fitness (population[X]);
      gradient[X][i] = fXi - oldfit[X];
      population[X][i] = 1-  population[X][i];
    }
  else
     fXi = oldfit[X] + gradient[X][i];
  return (fXi);
}
/*********-------------     pop[X] + pop[Y] -> New **************/
void crossover( X,  Y,  New) int  X,  Y, * New;
{ float diffXY=0;  int i, oneprime=getaprime(); int j=0;
  if (X==Y)
        for(i=0;i<indiv_length;i++)      New[i] = population[X][i];
  else{  switch (op_choice) {
            case 1: one_ptx(X,Y,New); break;
	    case 2: two_ptx(X,Y,New); break;
	    case 3: ux(X,Y,New);      break;
	    default : 
	      diffXY = oldfit[X]-oldfit[Y];
	      for(i=0;i<indiv_length;i++, j=(j+oneprime)%indiv_length)
		{
		New[j] = ((population[X][j] == population[Y][j]) 
			  ? population[X][j]
			  : ((diffXY > (fitness_modified(X,j) - fitness_modified(Y,j)))
			     ? population[X][j] 
			     : population[Y][j]));
	      }
	    }
       }
}
/*---------------------------------------------------------*/
float myrand() /* returns a float >= 0 and <1 */
{  long int n;
   jran = (1283 + (jran * 106)) % 6075; /* 0 <= jran < 6075 */
   return ((float) (jran) / 6075 );
}/***********************************************************************/
void mutation( X) int * X;
{ int i;
  for(i=0;i<indiv_length;i++)  if (myrand() < mut_rate)  X[i] = 1-X[i];
} /************* uniform crossover *****************/
int ux(x, y, Z) int x, y, * Z; {  return(pux(x,y,Z,0.5));}
/***********************************************************************/
void one_ptx( x, y, Z)/* x and y are indices; generates 1 child*/
int  x,  y, * Z;
{ int i, pick = (int) (myrand()* indiv_length);
  for(i=0;i<pick;i++)                Z[i] =  population[x][i];
  for(i=pick;i<indiv_length;i++)     Z[i] =  population[y][i];
}
/***********************************************************************/
void two_ptx(x,  y,  Z)
int x, y, * Z;
{ int i,pick2, pick1 = (int) (myrand()* indiv_length);
  pick2 = pick1+(int)(myrand()*(indiv_length-pick1));

  for(i=0;i<pick1; i++)              Z[i] =  population[x][i];
  for(i=pick1;i<pick2;i++)           Z[i] =  population[y][i];
  for(i=pick2;i<indiv_length;i++)    Z[i] =  population[x][i];
}/***********************************************************************
 * pux : parameterized uniform crossover, x and y are indices
 **********************************************************************/
int pux(x,  y, Z, p) int  x,  y,  * Z; float p;
{ int i;   float pick; 
  for(i=0;i<indiv_length;i++) 
    Z[i] =  (myrand()  > p) ? population[y][i]: population[x][i];
}/***********************************************************************/
void hux(x,  y, Z) /* ensures half the elements come from each parent */
int  x,  y,  * Z; 
{ int i, oneprime=getaprime(); int j,xcount=0; float halflen = indiv_length/2;
  for(i=0;i<indiv_length;i++)      Z[i] = population[x][i];
  j = (int) ( myrand() * indiv_length) ;
  for (i=0; i<halflen; i++, j=(j+oneprime)%indiv_length)
	   Z[j] = population[y][j];
 }/***********************************************************************/
float weight(X)  int * X; /* sum of elements of individual X  */
{ int i, sum=0;
  for(i=0;i<indiv_length;i++)       sum += X[i];
  return((float) sum);
}/***********************************************************************/
float max( x, y) float x,y; { return ((x >= y) ? x :y);}
/**********************************************************************/
float min(x,y) float x, y;{return ((x > y) ? y :x);}
/**********************************************************************/
float fitness(X) int * X; {   return(ugly3(sizeby3,X));}
/*********** CALCULATE Hamming distance*************/
float variation(X, Y) int * X, *Y;
{ int i, sum=0;
  for(i=0;i<indiv_length;i++)           if ( Y[i] != X[i]) sum ++;
  return((float) sum);
}/***********************************************************************/
float triple(i,j,k,A) int i,j,k, *A;
{switch(A[i]+A[j]+A[k])
	 {case 0: 	     return(0.9);
	  case 1: return(0.6);
	  case 2:  return(0.3);
	  default: return(1.0);
	 }
}/*****************************************/
float ugly3(n,A) int n,*A;
{int i; float fsum=0.0;
 for (i=0; i<n; i++)
   fsum= fsum + triple(i,i+n,i+n+n,A);
 return(fsum);
}/*****************************************/
float easy3(n,A)
int n,*A;
{int i; float fsum=0.0;
 for (i=0; i<3*n; i=i+3)
   fsum= fsum + triple(i,i+1,i+2,A);
 return(fsum);
}/*****************************************/
float getaprime(){ return (listofprimes[(int)(myrand() *  max_prime_count)]);}
/* ------------------------------------------------------------ */