/*Generalized GP, for many bins.  From gp9.c, on 20 August 1997.*/

#include <stdio.h>

#define MAX_BIN 16
  /* max. no. of possible bins */
#define MAX_PRINTS 200
  /* max. no. of times fitness values will be printed */
#define MAX_NUMBER_OF_NODES_POSSIBLE 200
#define MAX_POPULATION_SIZE 1000
#define aa 16807
#define mm 2147483647
#define qq 127773
#define rr 2836  
long int jran = 2147483647; long int fpseed = 2; /*all these help in random number generation*/
extern double drand48();
int nodecount=30;
int bincount=2; /*2 for bipartitioning */
int single_bin_count;

int op_choice=1;
char outfilename[80];
int mut_prob = 5; /*EXPRESSED AS A PERCENTAGE (5%) */
float crossover_prob = 100; /* EXPRESSED AS A PERCENTAGE (100%) */
float pop_sum_fitness=0.0;
float  pop_sum_cost = 0.0;
int max_num_of_generations, pop_size, location_of_best_individual=0; 
double cost[MAX_POPULATION_SIZE], w[MAX_NUMBER_OF_NODES_POSSIBLE][MAX_NUMBER_OF_NODES_POSSIBLE], fitness[MAX_POPULATION_SIZE];

int current_population[MAX_POPULATION_SIZE][MAX_NUMBER_OF_NODES_POSSIBLE], bestsol[MAX_NUMBER_OF_NODES_POSSIBLE];
FILE *fp; 

/*---------------------------------------------------------*/
float myrand() /* returns a float >= 0 and <1 */
{  long int n;
   jran = (1283 + (jran * 106)) % 6075; /* 0 <= jran < 6075 */
   return ((float) (jran) / 6075 );
}/***********************************************************************/
int preflip(prob) /* To help in random number generation -----------------*/
float prob;
{  long int lo, hi, test;  float u;

  hi = fpseed / qq;
  lo = fpseed % (long int) qq;
  test = aa*lo - rr*hi;

  if (test > 0) 
    fpseed = test;
  else
    fpseed = test + mm;

  return ( (srand(fpseed) / mm) >= prob);
}/***********************************************************************/
int rnd(low, high) /* Generates an integer between low and high (inclusive) */
int low, high;
{ int n;
  if (low == high) return low;
  else
    {  if (low >high) {n=low; low=high; high=n;}
       n = (int)( drand48() *(high-low+1)) + low;
       return( (n>high)? high: n);
   }
}/* ============================================================ */
void printind(xx) int xx; /* print current_population[xx]  ----------------*/
{int j; for (j=0; j<nodecount; j++)   fprintf(fp, "%1d", current_population[xx][j]);
}/*------------------------------------------------------------*/
void printinds(indiv_vector) int *indiv_vector; /*print vector */
{int j; for (j=0; j<nodecount; j++) fprintf(fp, "%1d", indiv_vector[j]);
}/* ------------------------------------------------------------ */
void printpop() /*print entire current population*/
{int i,j;
 for (i=0; i<pop_size; i++) {
    fprintf(fp, "\n%2d ", i);
    for (j=0; j<nodecount; j++)      fprintf(fp, "%1d", current_population[i][j]);
    fprintf(fp, " cost %3f\n", cost[i]);
  } 
    fprintf(fp, "\n");  
}/* ------------------------------------------------------------*/
double  calc_cost(xp) int *xp;
{int j, k;  double kost = 0.0;
 for (j = 0; j < nodecount; j++) 
    for (k = 0; k < j; k++) 
      if (*(xp+j) != *(xp+k)) kost += w[j][k];
 if (kost <0)
   { printf("ERROR calc_cost for ");
     for (j = 0; j < nodecount; j++) printf("%d",*(xp+j));
     printf("\n");
   }
 return kost;
}/*------------------------------------------------------------*/
double calc_pop_cost() /* Calculate costs for current pop. members, returns average */
{ int i, zmin=9999999, zmax=-999999; float avg_cost;

  location_of_best_individual= pop_sum_fitness= pop_sum_cost = 0;

  for (i=0; i<pop_size; i++)
   { cost[i] = calc_cost(current_population[i]);
     pop_sum_cost += cost[i];
     pop_sum_fitness += 1/(cost[i]*cost[i]);
     if (cost[i] > zmax) zmax = cost[i];
     if (cost[i] < zmin) {
	zmin = cost[i];
	location_of_best_individual = i;
      }
    }
  avg_cost = pop_sum_cost/pop_size;
  for (i=0; i<nodecount; i++) bestsol[i] = current_population[location_of_best_individual][i];
  for (i=0; i<pop_size; i++) 
    fitness[i] = (1.0/(cost[i]*cost[i] * pop_sum_fitness));
  return (avg_cost); 
}/* ------------------------------------------------------------ */
void generate_graph()
{  int 		i,j;  double u;
  for(i=0; i< nodecount; i++){
    w[i][i]=0;
    for(j=i; j < nodecount; j++){
      if (i!=j)
	  w[i][j]=(int)(21.0* drand48()); /*each weight is between 0 and 20*/
      w[i][i] += w[i][j];
    }
    for(j=0; j< i; j++){
      w[i][j]=w[j][i];
      w[i][i] += w[i][j];
    }
  }
}/*============================================================*/
void generate_initial_set()
{ int i, j, temp, count[MAX_BIN];
  for (i=0; i< pop_size; i++) 
    {for (j=0; j < bincount; j++) count[j]=0;
     for (j=0; j < nodecount; j++) 
       {do 
	  temp = (int) (myrand() * bincount) ;
        while (count[temp] == single_bin_count);
	current_population[i][j] = temp;
	count[temp] = count[temp] +1;
      }
   }
}/*============================================================*/
void mutate(xp) int *xp;
{int i;
  for (i=0; i<nodecount; i++)
    if (mut_prob > rnd(0,100)) xp[i]= (int)(drand48()*bincount); 
}/* ------------------------------------------------------------ */
void adjust_elt(xp) int *xp; /* To balance all bins */
{ int i,big,temp, count[MAX_BIN], bin,bigval, smallval; float epsilon=0.0001;
  for (temp=0; temp<bincount; temp++)    count[temp]=0;
  for (temp=0; temp<nodecount; temp++)  
    {bin = *(xp + temp); count[bin] = count[bin] +1;}

  big = (int)( drand48() * (nodecount-epsilon));
  for (i=0;(i<nodecount) && (count[xp[big]] <= single_bin_count);i++)
     big = ((big >= nodecount-1) ?0 : big+1);

  while (i<nodecount) /*implies that some bin occurs more often than single_bin_count*/
   {
    smallval = (int)( drand48() * bincount);
    for (i=0;(i<bincount) && (count[smallval] >= single_bin_count);i++)
       smallval = ((smallval >= bincount-1) ?0 : smallval+1);
    if (count[smallval] >= single_bin_count)
       printf("Didn't find small val!\n"); /*This shouldn't really happen*/
    bigval = *(xp+big); 
    *(xp+big)=smallval; 
    (count[bigval])--;     (count[smallval])++;

    big = (int)( drand48() * (nodecount-epsilon));
    for (i=0;(i<nodecount) && (count[xp[big]] <= single_bin_count);i++)
      big = ((big >= nodecount-1) ?0 : big+1);
  }
}/* ------------------------------------------------------------ */

void onept_xover(dadp, momp, xdp, xmp) /* combined with mutation and adjustment */
int dadp, momp, *xdp, *xmp; 
{ int i, startp, lnth, temp, count[MAX_BIN],bin;
  for (i=0; i<nodecount; i++)
    { *(xdp+i) = current_population[dadp][i];      *(xmp+i) = current_population[momp][i];
    }
  if (crossover_prob >= rnd(0,100))
        {  startp = rnd(0, nodecount-1);
	   for (i=startp; i<nodecount; i++)
	    { temp = *(xdp+i); *(xdp+i) = *(xmp+i); *(xmp+i) = temp;
	    }
	 }
  mutate(xdp);   mutate(xmp);
  adjust_elt(xdp);  
  adjust_elt(xmp);
}       /* ------------------------------------------------------------ */
void twopt_xover(dadp, momp, xdp, xmp) /* combined with mutation and adjustment */
int dadp, momp, *xdp, *xmp; 
{ int i, startp, endp, lnth, temp; 

  for (i=0; i<nodecount; i++)
    {
      *(xdp+i) = current_population[dadp][i];
      *(xmp+i) = current_population[momp][i];
    }
  
  if ((dadp != momp) &&  (crossover_prob >= rnd(0,100)))
        {  startp = rnd(0, nodecount-1);  endp = rnd(0, nodecount-1);
	   if (startp>endp) {temp=startp; startp=endp; endp=temp;}
	   for (i=startp; i<endp; i++)
	    {temp = *(xdp+i); *(xdp+i) = *(xmp+i); *(xmp+i) = temp;}
	 }
  mutate(xdp);   mutate(xmp);
  adjust_elt(xdp);
  adjust_elt(xmp);
}/* ------------------------------------------------------------ */
void uniform_xover(dadp, momp, xdp, xmp) /* combined with mutation and adjustment */
int dadp, momp, *xdp, *xmp; 
{ int i;
  for (i=0; i<nodecount; i++)
    if (drand48()>0.5)
      {*(xdp+i) = current_population[dadp][i];      *(xmp+i) = current_population[momp][i];    }
    else
      {*(xmp+i) = current_population[dadp][i];      *(xdp+i) = current_population[momp][i];    }
  mutate(xdp);   mutate(xmp);
  adjust_elt(xdp);
  adjust_elt(xmp);
}/* ------------------------------------------------------------ */
int get_parent() /*roulette wheel*/
{float current= fitness[0]; int i=0;     float randval = drand48();
 while (randval > current)      { i++; current += fitness[i];}
 return((i<pop_size) ? i : pop_size-1);
}/*----------------------------------------------------*/
void reproduction()
{ int i, j, ngene[MAX_POPULATION_SIZE][MAX_NUMBER_OF_NODES_POSSIBLE], xstrd[MAX_NUMBER_OF_NODES_POSSIBLE], xstrm[MAX_NUMBER_OF_NODES_POSSIBLE], pa, ma;
  for (i=0; i<pop_size; i += 2)
    {switch (op_choice) 
	{	
	 case 2: pa = get_parent(); while (pa == (ma=get_parent()));
	  twopt_xover(pa,ma, xstrd, xstrm);break;

	 case 3: pa = get_parent(); while (pa == (ma=get_parent()));
	  uniform_xover(pa,ma, xstrd, xstrm);break;

	 default: pa = get_parent(); while (pa == (ma=get_parent()));
	  onept_xover(pa,ma, xstrd, xstrm);
	 }
      for (j=0; j<nodecount; j++)
	{
	  ngene[i][j] = xstrd[j];
	  ngene[i+1][j] = xstrm[j];
	}
   }
  for (i=0; i<pop_size; i++) 
    if ((i != location_of_best_individual)|| (calc_cost(ngene[i]) < cost[location_of_best_individual]))
      for (j=0; j<nodecount; j++) current_population[i][j] = ngene[i][j];
 }
/* ------------------------------------------------------------ */
main(argc, argv)
int argc;
char *argv[];
{
  int i, j, worsti, count0, count1;
  double av_cost, sum_init_cost=0, sum_best_cost=0, zmax, AvCost[MAX_PRINTS], BestCost[MAX_PRINTS];
  int rseed = 10; 
  int print_freq, print_num, print_count, trial_num = 0;
  long int clock;  float  avg_time;

if (argc==1)
    {
      printf("Usage:\n");
      printf("gp op_choice bin_count num_nodes[1-200] population[1-1000] epochs[1-10000] mutation_percent[1-100] print_frequency outputfile\n");
      printf("op_choice: 1 for 1PTX, 2 for 2PTX, 3 for UX \n");
      return;
    }
else    {
  sscanf(argv[1],"%d", &op_choice);
  sscanf(argv[2],"%d", &bincount);
  sscanf(argv[3],"%d", &nodecount);
  sscanf(argv[4], "%d", &pop_size);
  sscanf(argv[5], "%d", &max_num_of_generations);
  sscanf(argv[6], "%d", &mut_prob);
  sscanf(argv[7], "%d", &print_freq);

  print_count = max_num_of_generations/print_freq;
  if (print_count > MAX_PRINTS) 
  { print_count = MAX_PRINTS; print_freq = max_num_of_generations/print_count;}

  if ((fp = fopen(argv[8], "w+t")) == NULL)
    {
      printf("Could not open output file; please retype filename:\n");
      scanf("%d",outfilename);
      if ((fp = fopen(outfilename, "w+t")) == NULL)
	{
	  printf("Could not open file again; please check directory permissions.  Bye.\n");
      exit(0);
	}
    }
  rewind(fp);
}/* closing else */

 single_bin_count= nodecount/bincount;
  
 generate_graph();

fprintf(fp, "Nodes = %d, Population Size = %d, Max. no. of generations = %d, Operator = ",
	  nodecount, pop_size, max_num_of_generations);
  switch (op_choice)
{ case 2:   fprintf(fp, "2PTX"); break;
  case 3:    fprintf(fp, "UX"); break;
  default:    fprintf(fp, "1PTX"); 
}
clock = time(0);

for (i=0; i< print_count; i++)   BestCost[i] = AvCost[i] = 0; 

  generate_initial_set();  av_cost = calc_pop_cost();
    fprintf(fp,"\nINITIAL avg.=  %10.2f, best = %10.2f,",
	       av_cost, 
	       cost[location_of_best_individual]);
    fprintf(fp, " individual= "); printind(location_of_best_individual); fprintf(fp, "\n");

  sum_init_cost =   sum_init_cost + av_cost; 
  sum_best_cost =   sum_best_cost + cost[location_of_best_individual];

  for (i=0; i< max_num_of_generations/print_freq; i++)
     { for (j=0; j<print_freq; j++) 
	 {reproduction();    av_cost = calc_pop_cost(); }
       AvCost[i] =  AvCost[i] + av_cost; BestCost[i] += cost[location_of_best_individual];
     }
    fprintf(fp," FINAL avg.=  %10.2f, best = %10.2f",
	       av_cost, 
	       cost[location_of_best_individual]);
    fprintf(fp, ", individual= "); printind(location_of_best_individual); fprintf(fp, "\n");

 avg_time =  ((float)(time(0) - clock))/max_num_of_generations;

fprintf(fp,"\nGenerations, Average Cost, Best Cost\n");
printf("\nGenerations, Average Cost, Best Cost\n");
fprintf(fp,"------------------------------------\n");
printf("------------------------------------\n");

    fprintf(fp,"   0       %10.2f       %10.2f\n", (float)(sum_init_cost), (float)(sum_best_cost) );
    printf("0       %10.2f       %10.2f\n", (float)(sum_init_cost), (float)(sum_best_cost));
   
for (i=0; i<max_num_of_generations/print_freq; i++)
{ print_num= (i+1)*print_freq;
  if (print_num <10)    fprintf(fp, "   ");
  else if (print_num <100)  fprintf(fp, "  ");
  else if (print_num <1000)  fprintf(fp, " ");

    fprintf(fp,"%d       %10.2f       %10.2f\n", print_num, AvCost[i], BestCost[i]);
     printf("%d       %10.2f       %10.2f\n", print_num, AvCost[i], BestCost[i]);
   }
  fprintf(fp, "Average time per generation = %8.4f\n", avg_time);
  fclose(fp);
}/*------------------------------------------------------------*/