/*
  Last changed Time-stamp: <2006-12-05 16:57:56 alex>
  
  Simple Genetic Algorithm with an external evaluation function.

  You can find several implementations of this code on the 
  internet. As far as I know the reference therefor is
  
  Genetic Algorithms + Data Structures = Evolution Program
  Third Edition, Springer, 1996
  Author: Z. Michalewicz
  
  The original program is by Denis Cormier. 
  Modified by Alexander Donath.
  
  Uses the Mersenne Twister which can be found here:
  http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
 
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>

#include "mt19937ar.h"

/* individual */
/* change the members of the record to match your needs */
typedef struct {
  char *name;          /* rna string */
  double values[6];    /* vector of variables (k1*h, k_1, k2, k3, k4, k5) */
  double fitness;      /* GT's fitness */
  double rfitness;
  double cfitness;
} genotype_rna;

/* arrays for parent and offspring (new) population */
static genotype_rna *population = NULL;
static genotype_rna *newpopulation = NULL;

/* change those parameteres to match your needs  */
static double (*evaluateFitness)(char*, double*);
static int POPSIZE = 4;
static int NVARS = 6;
static int MAX;
static FILE *outfile;

/* 
 * returns random number from interval [from, to] 
 */
static double random_value (double from, double to)
{
  return (from + genrand_real2 ()*(to-from+1.0));
}

/* 
 * creates a random string with characters out of [ACGU] 
 */
static void random_string (char *string)
{
  int i, random;
 
  string[MAX] = '\0';
 
  for (i = 0; i < MAX; i++)
    {
      random = (int)random_value (0.0, 3.0);
        
      switch (random)
	{
	case 0: 
	  string[i] = 'A';
	  break;
	case 1: 
	  string[i] = 'C';
	  break;
	case 2: 
	  string[i] = 'G';
	  break;
	case 3: 
	  string[i] = 'U';
	  break;
	default: 
	  printf ("Something is wrong with random_string (random = %i)\n", random);
	  break; 
	}
    }
}

/* 
 * returns random nucleotide 
 */
static void random_nt (char *sequence, int pos)
{
  char abc[3];
  int random = (int)(random_value(0.0, 2.0)); 
 
  switch(sequence[pos])
    {
    case 'A': 
      strcpy(abc, "CGU");
      break;
    case 'C': 
      strcpy(abc, "AUG");
      break;
    case 'G': 
      strcpy(abc, "UCA");
      break;
    case 'U': 
      strcpy(abc, "GAC");
      break;
    default: 
      printf("Something is wrong with your sequence!\n");
      break;
    }
  sequence[pos] = abc[random];
}

/* 
 * swap two chars 
 */
static void swap(char *x, char *y)
{
  char temp;
  
  temp = *x;
  *x = *y;
  *y = temp;
}

/* 
 * allocates needed memory 
 */
static void allocate_memory (int nvars, int popsize)
{
  int i;
  
  if (popsize > 0) 
    POPSIZE = popsize;
  if (nvars > 0) 
    NVARS = nvars;

  /* initialize mersene twister */
  init_genrand (time(0));
  
  /* allocate memory */
  population = (genotype_rna *) calloc (POPSIZE+1, sizeof(genotype_rna));
  for (i = 0; i < POPSIZE+1; i++)
    {
      population[i].name = (char *) malloc((MAX+1) * sizeof(char));
      population[i].fitness = 0.0;
      population[i].rfitness = 0.0;
      population[i].cfitness = 0.0;
    }
  
  newpopulation = (genotype_rna *) calloc (POPSIZE+1, sizeof(genotype_rna));
  for (i = 0; i < POPSIZE+1; i++)
    {
      newpopulation[i].name = (char *) malloc((MAX+1) * sizeof(char));
      newpopulation[i].fitness = 0.0;
      newpopulation[i].rfitness = 0.0;
      newpopulation[i].cfitness = 0.0;
    }
}

/* 
 * releases allocated memory 
 */
static void finalize (void)
{
  int i,j;
 
  for (i = 0; i < POPSIZE+1; i++)
    {
      free (population[i].name);
      free (newpopulation[i].name);
    } 
  
  free (population);
  free (newpopulation);
  
  population = NULL;
  newpopulation = NULL;
}

/* 
 * generates a random population of rna strings 
 */
static void random_population (void)
{
  int i;
  
  for (i = 0; i < POPSIZE; i++) 
    random_string (population[i].name);
}

/* 
 * calculates fitness values for each genotype 
 */
static void evaluate (genotype_rna *pop)
{
  int i;
  
  for (i = 0; i < POPSIZE; i++) 
    pop[i].fitness = evaluateFitness (pop[i].name, pop[i].values);
}

/* 
 * keep the fittest 
 */
static void keep_the_best (genotype_rna *pop)
{
  int i, current_best = 0;

  /* locate index of fittest genotype */
  for (i = 1; i < POPSIZE; i++)
    {
      if (pop[i].fitness>pop[current_best].fitness) 
	current_best = i;
    }
  
  /* copy fittest genotype to last position of population array */
  strcpy (pop[POPSIZE].name, pop[current_best].name);
  memcpy (pop[POPSIZE].values, pop[current_best].values, (size_t)NVARS*sizeof (double));
  pop[POPSIZE].fitness = pop[current_best].fitness;
}

/* 
 * selects individuals for mating 
 */
static void selection (void)
{
  int i, j;
  double p, total_fitness = 0.0;

  /* calculate total fitness of the population */
  for (i = 0; i < POPSIZE; i++) 
    total_fitness += population[i].fitness;

  /* calculate relative fitness */
  for (i = 0; i < POPSIZE; i++) population[i].rfitness = population[i].fitness/total_fitness;

  /* calculate cumulative fitness */
  population[0].cfitness = population[0].rfitness;
  for (i = 1; i < POPSIZE; i++) 
    population[i].cfitness = population[i-1].cfitness + population[i].rfitness;
  
  /* select survivors using cumulative fitness -> save them to newpopulation */
  for (i = 0; i < POPSIZE; i++)
    {
      p = genrand_real2 ();
      for (j = 0; j < POPSIZE; j++)
	{
	  if (population[j].cfitness >= p)
	    {
	      strcpy (newpopulation[i].name, population[j].name);
	      memcpy (newpopulation[i].values, population[j].values, (size_t)NVARS*sizeof(double));
	      newpopulation[i].fitness = population[j].fitness;
	      break;
	    }
	}
    }
}

/* 
 * does the crossover 
 */
static void Xover (int one, int two)
{
  int i, point;

  /* select crossover point */
  point = (int)random_value (0, MAX-1);

  for (i = 0; i < point; i++) 
    swap (&newpopulation[one].name[i], &newpopulation[two].name[i]);
}

/* 
 * select genotypes for crossover 
 */
static void crossover (double pxover)
{
  int i, one = 0, first = 0; /* 'first' counts the number of members chosen */
  double x;
  
  for (i = 0; i < POPSIZE; i++)
    {
      x = genrand_real2();
      if (x < pxover)
	{
	  if ((++first)%2 == 0) 
	    Xover (one, i);
	  else 
	    one = i;
	}
    }
}

/* 
 * randomly mutate nucleotides 
 */
static void mutate (double pmutation)
{
  int i, j;
  double x;
  
  for (i = 0; i < POPSIZE; i++)
    {
      for (j = 0; j < MAX; j++)
	{
	  x = genrand_real2 ();
	  if (x < pmutation) 
	    random_nt (newpopulation[i].name, j);
	}
    }
}

/*
 * no comment. there's enough of this within the code.
 */
static void elitist (void)
{
  int i, j, best1, best2, worst2;

  /* elitism replacement scheme (size of elite: 1) */
  /* search fittest individual */
  best1 = 0;
  best2 = 0;
  worst2 = 0;
  /* fittest from old population */
  for (i = 1; i < POPSIZE; i++)
    {
      if (population[i].fitness > population[best1].fitness) 
	best1=i;
    }
  /* fittest and worst from new population */
  for (i = 1; i < POPSIZE; i++)
    {
      if (newpopulation[i].fitness > newpopulation[best2].fitness) 
	best2=i;
      if (newpopulation[i].fitness > newpopulation[worst2].fitness) 
	worst2=i;
    }
  
  if (newpopulation[best2].fitness > population[best1].fitness) /* fitness increase */
    {
      /* copy all individuals */
      for (i = 0; i < POPSIZE; i++)
	{
	  strcpy (population[i].name, newpopulation[i].name);
	  memcpy (population[i].values, newpopulation[i].values, (size_t)NVARS*sizeof(double));
	  population[i].fitness = newpopulation[i].fitness;
	}
    }
  else /* fitness decrease */
    {
      /* save best of old population */
      if (best1)
	{
	  strcpy (population[0].name, population[best1].name);
	  memcpy (population[0].values, population[best1].values, (size_t)NVARS*sizeof(double));
	  population[0].fitness = population[best1].fitness;
	}

      /* copy all but worst individual */
      j = 1;
      for (i = 0; i < POPSIZE; i++)
	{
	  if (i != worst2)
	    {
	      strcpy (population[j].name, newpopulation[i].name);
	      memcpy (population[j].values, newpopulation[i].values, (size_t)NVARS*sizeof(double));
	      population[j].fitness = newpopulation[i].fitness;
	      j++;
	    }
	}
    }
}

/* 
 * report to outfile 
 */
static void report (int generation)
{
  int i;
  double sum = 0.0, square_sum = 0.0, stddev, best_value;
  
  for (i = 0; i < POPSIZE; i++)
    {
      sum += population[i].fitness;
      square_sum += (population[i].fitness * population[i].fitness);
    }

  sum /= (double) POPSIZE;
  square_sum /= (double) POPSIZE;
  
  /* square_sum = sum * sum;*/
  /* stddev = sqrt((square_sum - sum_square)/(double)(POPSIZE-1)); */
  stddev = sqrt (square_sum - sum * sum);
  best_value = population[POPSIZE].fitness;

  fprintf (outfile, "%5d\t%1.6f\t%1.6f\t%1.6f\n", generation-1, best_value, sum, stddev);
  fflush (outfile);
}

/* 
 * call this function to use the GA
 * pmutation: probability of a point mutation
 * pxover: probability of a crossover
 */
void fit_parameter (int popsize, int nvars, int maxgens, int length, double pmutation, double pxover, double (*EF)(char*, double*))
{
  int i, generations = 0;
  char filename[30];
  
  sprintf(filename, "m%1.3f_x%1.3f", pmutation, pxover);
  outfile = fopen(filename,"wt");
  MAX = length;
  /* the external evaluation function */
  evaluateFitness = EF;
  allocate_memory (nvars, popsize);
  random_population ();
  
  while (generations++ < maxgens)
    {
      evaluate (population);
      keep_the_best (population);
      report (generations);
      selection ();
      crossover (pxover);
      mutate (pmutation);
      evaluate (newpopulation);
      keep_the_best (newpopulation);
      elitist ();
    }
  /* last generation */
  keep_the_best (population);
  report (generations);

  finalize();
  fclose(outfile); 
}