/* $Revision: 1.3 $ */

/*WBL 13 Apr 2016 for twins and 6-mux */
/*use -D on gcc command line for selection and popsize */
/*WBL 31 Mar 2010 for 6-mux */

/*
 * Program:      tiny_gp_b.c
 *
 * Author:       Riccardo Poli (email: R.Poli@cs.bham.ac.uk)
 *               School of Computer Science
 *               The University of Birmingham
 *               UK
 *
 * Creation date: August 2000
 *
 * Description: A tiny C implementation of GP which works with Boolean problems.
 */

/* This file is Copyright (c) 2000 Riccardo Poli.  All rights
 * reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met: 1. Redistributions of source code must retain the above
 * copyright notice, this list of conditions and the following
 * disclaimer.  2.  Redistributions in binary form must reproduce the
 * above copyright notice, this list of conditions and the following
 * disclaimer in the documentation and/or other materials provided
 * with the distribution.  3.  Neither the name of the Author nor the
 * names of its contributors may be used to endorse or promote
 * products derived from this software without specific prior written
 * permission.
 *
 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESS OR IMPLIED
 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
 * OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
 * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
 * DAMAGE.
 */

#include <stdio.h>
#include <stdlib.h>
#include <sys/types.h>
#include <time.h>
#include <string.h>
#include <assert.h>

enum {X1, X2, X3, X4, X5, X6, /*X7, X8, X9, X10,*/ AND, OR, NOR, NAND, XOR, EQ /* ,NOT */  };

#define FSET_START AND
#define FSET_END NAND

int x1, x2, x3, x4, x5, x6, x7, x8, x9, x10;

char *program;

float fbest = 0.0;
float avg_len;
long seed;


#ifndef SHORT_OUTPUT
#define SHORT_OUTPUT 0
#endif

/* #define DIVERSITY */
#define MAX_LEN 64

#ifndef NTH_KILL
#define NTH_KILL 2
#endif

#ifndef CAN_KILL_BEST
#define CAN_KILL_BEST 1
#endif

#ifndef POPSIZE
#define POPSIZE 1000
#endif

#ifndef DEPTH
#define DEPTH   6
#endif

#ifndef PMUT
#define PMUT    0.05
#endif

#ifndef GENERATIONS
#define GENERATIONS 100
#endif

#ifndef TSIZE
#define TSIZE 2
#endif

#define	INT_MAX 2147483647

unsigned long run() /* Interpreter */
{
  switch ( *program++ )
    {
    case X1  : return( x1 );
    case X2  : return( x2 );
    case X3  : return( x3 );
    case X4  : return( x4 );
    case X5  : return( x5 );
    case X6  : return( x6 );
      /*
    case X7  : return( x7 );
    case X8  : return( x8 );
    case X9  : return( x9 );
    case X10 : return( x10 );
      */
      /* case NOT : return 1 & ( ~run() ); */
    case AND : return 1 & ( run() & run() );
    case OR  : return 1 & ( run() | run() );
    case XOR : return 1 & ( run() ^ run() );
    case NAND : return 1 & ( ~ (run() & run()) );
    case NOR  : return 1 & ( ~ (run() | run()) );
    case EQ  : return  1 & ( ~( run() ^ run()) );
    }
  return 0; // This shouldn't happen! (making eclipse happy)
}

char *traverse( char *buffer )
{
  switch(*buffer)
    {
    case X1:
    case X2:
    case X3:
    case X4:
    case X5:
    case X6:
      /*
    case X7:
    case X8:
    case X9:
    case X10:
      */
      return( ++buffer );

/*    case NOT: return( traverse( ++buffer ) );
*/
    case AND:
    case OR:
    case XOR:
    case NAND:
    case NOR:
    case EQ:
      return( traverse( traverse( ++buffer ) ) );

    }
  return ""; // This shouldn't happen! (making eclipse happy)
}

unsigned long mux6(/*bool test[6]*/) { 
  const int a = 
    x1 /*test[0]*/+
    2*x2; //test[1]+
    //2;
  //return test[a];
  switch(a) {
  case 0 : return x3;
  case 1 : return x4;
  case 2 : return x5;
  case 3 : return x6;
  }
}

unsigned long e10parity() /* Even-10 parity function */
{
  return 1 & (~(x1^x2^x3^x4^x5^x6^x7^x8^x9^x10));
}

unsigned long targets[1024];

void setup_fitness()
{
  int i = 0;

  for( x1 = 0; x1 < 2; x1 ++ )
  for( x2 = 0; x2 < 2; x2 ++ )
  for( x3 = 0; x3 < 2; x3 ++ )
  for( x4 = 0; x4 < 2; x4 ++ )
  for( x5 = 0; x5 < 2; x5 ++ )
  for( x6 = 0; x6 < 2; x6 ++ )
    /*
  for( x7 = 0; x7 < 1; x7 ++ )
  for( x8 = 0; x8 < 1; x8 ++ )
  for( x9 = 0; x9 < 1; x9 ++ )
  for( x10 = 0; x10 < 1; x10 ++ )
    */
    targets[i++] = mux6(); //e10parity();
/*		  printf("%lx\n",targets[i-1]); */
}


unsigned long nodes_evaluated = 0;

int length( char *Prog ) 
{
  return traverse( Prog ) - Prog;
}

int fitness_function( char *Prog ) /* Even-10 parity fitness function */
{
  int fit = 0, i = 0, len;
  unsigned result;

  len = traverse( Prog ) - Prog;

  for( x1 = 0; x1 < 2; x1 ++ )
  for( x2 = 0; x2 < 2; x2 ++ )
  for( x3 = 0; x3 < 2; x3 ++ )
  for( x4 = 0; x4 < 2; x4 ++ )
  for( x5 = 0; x5 < 2; x5 ++ )
  for( x6 = 0; x6 < 2; x6 ++ )
    /*
  for( x7 = 0; x7 < 1; x7 ++ )
  for( x8 = 0; x8 < 1; x8 ++ )
  for( x9 = 0; x9 < 1; x9 ++ )
  for( x10 = 0; x10 < 1; x10 ++ )
    */
    {
      program = Prog;
      result = run();
      nodes_evaluated += len;

      if ( result == targets[i++] ) fit ++;
    }

  /*  if ( ( fit < fbest || CAN_KILL_BEST) && len > avg_len )
    if ( lrand48() % NTH_KILL == 0 )
      return( 0 );
  */
  return( fit );
}


int grow( char *buffer, int pos, int max, int depth )
{
  char prim = lrand48() % 2;

  if ( pos >= max )
    return( -1 );

  if ( pos == 0 )
    prim = 1;

  if ( prim == 0 || depth == 0 )
    {
      prim = lrand48() % (FSET_START - X1) + X1;
      buffer[pos] = prim;
      return(pos+1);
    }
  else
    {
      prim = lrand48() % (FSET_END - FSET_START + 1) + FSET_START;
      switch(prim)
	{
	  /*
	case NOT: buffer[pos] = prim;
	  return( grow( buffer, pos+1, max, depth-1) );
	  */
	case AND:
	case OR:
	case XOR:
	case NAND:
	case NOR:
	case EQ:

	  buffer[pos] = prim;
	  return( grow( buffer, grow( buffer, pos+1, max,depth-1), max,depth-1 ) );
	}
    }
  return -1; // This shouldn't happen! (making eclipse happy)
}

char *print_indiv( char *buffer )
{
  switch(*buffer)
    {
    case X1:
    case X2:
    case X3:
    case X4:
    case X5:
    case X6:
      /*
    case X7:
    case X8:
    case X9:
    case X10:
      */
      printf( "X%d ", *buffer + 1);
      return( ++buffer );

/*    case NOT: printf( "NOT ");
      return( print_indiv( ++buffer ) );
*/
    case AND: printf( "AND(");
      buffer = print_indiv( ++buffer );
      break;
    case OR: printf( "OR(");
      buffer = print_indiv( ++buffer );
      break;
    case XOR: printf( "XOR(");
      buffer = print_indiv( ++buffer );
      break;
    case NAND: printf( "NAND(");
      buffer = print_indiv( ++buffer );
      break;
    case NOR: printf( "NOR(");
      buffer = print_indiv( ++buffer );
      break;
    case EQ: printf( "EQ(");
      buffer = print_indiv( ++buffer );
      break;
    }
  printf( ",");
  buffer = print_indiv( buffer );
  printf( ")");
  return buffer;
  //return ""; // This shouldn't happen! (making eclipse happy)
}


char *create_random_indiv( int depth )
{
  char buffer[MAX_LEN];
  char *ind;
  int len;

  len = grow( buffer, 0, MAX_LEN, depth );


  while (len < 0 )
    len = grow( buffer, 0, MAX_LEN, depth );

  ind = (char *) malloc( len * sizeof(char));

  memcpy( ind, buffer, len );
  return( ind );
}

char **create_random_pop(int n, int depth, int *fitness )
{
  char **pop = (char **) malloc( n * sizeof(char *));
  int i;

  for ( i = 0; i < n; i ++ )
    {
      pop[i] = create_random_indiv( depth );
      fitness[i] = fitness_function( pop[i] );
#ifdef DEBUG
  print_indiv( pop[i]);
  printf( "\n");
#endif
    }
  return( pop );
}


int stats( int *fitness, char **pop )
{
  int i, best = lrand48() % POPSIZE, node_count = 0;
  fbest = fitness[best];

  for ( i = 0; i < POPSIZE; i ++ )
    {
      node_count +=  traverse( pop[i] ) - pop[i];

      if ( fitness[i] > fbest )
		  {
			  best = i;
			  fbest = fitness[i];
		  }
    }
  avg_len = (float) node_count / POPSIZE;
  return ( best );
}

int min(const int x, const int y) {return (x<=y)? x : y;}
int max(const int x, const int y) {return (x>=y)? x : y;}
int twin_fitness(const int *fitness, const int competitor)
{
#ifdef RAND
  return 0;
#endif
#ifdef KIN
  return fitness[competitor ^ 1]; //kin selection
#endif
#ifdef MIN
  return min(fitness[competitor & (~0)], fitness[competitor & (~1)]);
#endif
#ifdef MAX
  return max(fitness[competitor & (~0)], fitness[competitor & (~1)]);
#endif
#ifdef MEAN
  return fitness[competitor & (~0)] + fitness[competitor & (~1)];
#endif
  return fitness[competitor];
}

int tournament( const int *fitness, const int tsize, const int pend_replace)
{
  int best, fbest = -1, i, competitor;

  for ( i = 0; i < tsize; i ++ )
    {
      int ntries=0;
      do {
	competitor = lrand48() % POPSIZE;
	assert(++ntries<1000);
      } while(competitor==pend_replace);

      if ( twin_fitness(fitness,competitor) > fbest )
	{
	  fbest = twin_fitness(fitness,competitor);
	  best = competitor;
	}
    }
  return( best );
}

int negative_tournament( int *fitness, int tsize )
{
  int worst, fworst = INT_MAX, i, competitor;

  for ( i = 0; i < tsize; i ++ )
    {
      competitor = lrand48() % POPSIZE;
      if ( twin_fitness(fitness,competitor) < fworst )
	{
	  fworst = twin_fitness(fitness,competitor);
	  worst = competitor;
	}
    }
  return( worst );
}

int replace(int *fitness)
{
  return negative_tournament(fitness, TSIZE);
}

char *crossover1( char *parent1,  char *parent2,
		  char *xo1start, char *xo2start)
{
  char *xo1end, *xo2end, *offspring;
  int len1 = traverse( parent1 ) - parent1;
  int len2 = traverse( parent2 ) - parent2;
  int lenoff;

  xo1end = traverse( xo1start );
  xo2end = traverse( xo2start );

  lenoff = (xo1start - parent1) + (xo2end - xo2start) + (parent1+len1-xo1end);

  offspring = (char *) malloc( lenoff * sizeof(char));

  memcpy( offspring, parent1, (xo1start - parent1) );
  memcpy( offspring + (xo1start - parent1), xo2start,  (xo2end - xo2start) );
  memcpy( offspring + (xo1start - parent1) + (xo2end - xo2start), xo1end,
	  (parent1+len1-xo1end) );

  return( offspring );
}

char *xopoint( char *parent)
{
  int len = traverse( parent ) - parent;
  return parent + lrand48() % len;
}

char *crossover( char *parent1, char *parent2 )
{
  char *xo1start, *xo2start;
  xo1start =  xopoint(parent1);
  xo2start =  xopoint(parent2);

  return crossover1(parent1, parent2, xo1start, xo2start);
}

char *mutation( char *parent, const double pmut )
{
  char *mutsite;

  const int len = traverse( parent ) - parent; int i;

  //printf("mutation( ");  print_indiv( parent );

  for (i = 0; i < len; i ++ )
    {
      if ( drand48() < pmut )
	{
	  mutsite =  parent + i;
	  switch(*mutsite)
	    {
	    case X1:
	    case X2:
	    case X3:
	    case X4:
	    case X5:
	    case X6:
	      /*
	    case X7:
	    case X8:
	    case X9:
	    case X10:
	      */
	      *mutsite = lrand48() % (FSET_START - X1 ) + X1;
	      break;

	    case AND:
	    case OR:
	    case XOR:
	    case NAND:
	    case NOR:
	    case EQ:
	      *mutsite = lrand48() % (FSET_END - FSET_START + 1) + FSET_START;
	    }
	}
    }
  //printf(" --> ");  print_indiv( parent );  printf(") mutation\n");

  assert(len == traverse( parent ) - parent);
  return( parent );
}


void print_parms()
{
  printf("\n\n----------------------------------\n6 Multiplexor problem twin $Revision: 1.3 $\n----------------------------------\n");
#ifdef RAND
  printf("RAND");
#endif
#ifdef KIN
  printf("KIN");
#endif
#ifdef MIN
  printf("MIN");
#endif
#ifdef MAX
  printf("MAX");
#endif
#ifdef MEAN
  printf("MEAN");
#endif
  printf("\n");
  printf("SEED=%ld\nMAX_LEN=%d\nNTH_KILL=%d\nCAN_KILL_BEST=%d\nPOPSIZE=%d\nDEPTH=%d\nPMUT=%lg\nGENERATIONS=%d\nTSIZE=%d\n----------------------------------\n",
 seed, MAX_LEN, NTH_KILL, CAN_KILL_BEST, POPSIZE, DEPTH, PMUT, GENERATIONS, TSIZE );
}


int main(int argc,char **argv)
{
  int gen = 0, indivs, best;
  int *fitness = (int *) calloc( POPSIZE, sizeof(int));
  char **pop, *newind;

  setup_fitness();

  if ( argc == 2 )
    seed = atoi(argv[1]);
  else
    seed = (long) time( (time_t *) 0 );

  print_parms();

  srand48(seed);
  pop = create_random_pop(POPSIZE, DEPTH, fitness );

  best = stats( fitness, pop );
  printf("Generation=%d Best=pop[%d] Fitness=%d Nodes Eval=%lu  Avg_size=%lg\n",gen, best,fitness[best], nodes_evaluated,avg_len);
  if (!SHORT_OUTPUT) {
    print_indiv( pop[best] );
    printf( "\n");
  }

  for ( gen = 1; gen < GENERATIONS; gen ++ )
    {
      for ( indivs = 0; indivs < POPSIZE; indivs +=2 )
	{
	  if (  lrand48() % 5 ) // if (1) //(lrand48() % 6)>1 && indivs < POPSIZE-1)  
	    {
	      const int parent1 = tournament( fitness, TSIZE, -1 );
	      const int parent2 = tournament( fitness, TSIZE, parent1 );
	      assert(parent1 != parent2);
	      int offspring1 = replace(fitness);
	      int offspring2 = offspring1 ^ 1; //lock twins together
	      assert(offspring1 != offspring2);
	      //newind = crossover( pop[parent1],pop[parent2] );
	      char *xo1start =  xopoint(pop[parent1]);
	      char *xo2start =  xopoint(pop[parent2]);
	      char *newind1 = crossover1(pop[parent1], pop[parent2], xo1start, xo2start);
	      char *newind2 = crossover1(pop[parent2], pop[parent1], xo2start, xo1start);

	      const int newfit1 = fitness_function( newind1 );
	      const int newfit2 = fitness_function( newind2 );
	      /*
	      printf("Generation=%d crossover %d(%d) F%d %d(%d) F%d => %d(%d) F%d %d(%d) F%d were L%d F%d L%d F%d\n",gen, 
		     parent1,length(pop[parent1]),fitness[parent1],
		     parent2,length(pop[parent2]),fitness[parent2],
		     offspring1,length(newind1),newfit1,
		     offspring2,length(newind2),newfit2,
		     length(pop[offspring1]),fitness[offspring1],length(pop[offspring2]),fitness[offspring2]);
	      */
	      assert(length(pop[parent1])+length(pop[parent2]) == length(newind1)+length(newind2));
	      free( pop[offspring1] );
	      free( pop[offspring2] );
	      pop[offspring1] = newind1;
	      pop[offspring2] = newind2;
	      fitness[offspring1] = newfit1;
	      fitness[offspring2] = newfit2;
	      //indivs += 2; 
	    }
	  else
	    {
	      const int parent1    = tournament( fitness, TSIZE, -1);
	      const int offspring1 = replace(fitness);

	      //printf("Generation=%d ",gen);
	      int i;
	      for(i=0;i<2;i++){
	      const int parent    = (i==0)?    parent1 : tournament( fitness, TSIZE, parent1);
	      const int offspring = (i==0)? offspring1 : offspring1 ^ 1; //lock twins together
	      assert(offspring >= 0 && offspring < POPSIZE);
	      const int len = traverse( pop[parent] ) - pop[parent];
	      newind = mutation( (char *) memcpy( (char *) malloc( len * sizeof(char)),
						  pop[parent], len ), PMUT );
	      const int newfit = fitness_function( newind );
	      /*
	      printf("mutate %d(%d) F%d => %d(%d) F%d was L%d F%d",
		     parent,   length(pop[parent]),fitness[parent],
		     offspring,length(newind),newfit,
		     length(pop[offspring]),fitness[offspring]);
	      if(i==0) printf(" "); else printf("\n");
	      */
	      free( pop[offspring] );
	      pop[offspring] = newind;
	      fitness[offspring] = newfit;
	      //indivs ++;
	      }//endfor mutate twice
	    }
	}
      best = stats( fitness, pop );
      printf("Generation=%d Best=pop[%d] Fitness=%d Nodes Eval=%lu Avg_size=%lg\n",gen, best,fitness[best], nodes_evaluated,avg_len);
      if (!SHORT_OUTPUT) {
	print_indiv( pop[best] );
	printf( "\n");
      }
      if ( fitness[best] == 64 ) return( 0 );
    }

	return EXIT_SUCCESS;
}