/* $Revision: 1.4 $ */

/*WBL 3 Oct 2020 generational version of Tiny GP for SIGEVOlution 13(3)
     "Multi-threaded Memory Efficient Crossover in C++ for
      Generational Genetic Programming", W. B. Langdon, 2020.
     Based on C++ code in ArXiv 2009.10460

     NB with efficient use of newpop because process new children in a
     different order pseudo random numbers are used in a different
     order hence results of crossover and mutation will be different
     in detail */

/*
 * 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 EQ

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 MAX_LEN 1000

#ifndef NTH_KILL
#define NTH_KILL 2
#endif

#ifndef CAN_KILL_BEST
#define CAN_KILL_BEST 1
#endif

#ifndef POPSIZE
#define POPSIZE 100
#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 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 < 2; x7 ++ )
  for( x8 = 0; x8 < 2; x8 ++ )
  for( x9 = 0; x9 < 2; x9 ++ )
  for( x10 = 0; x10 < 2; x10 ++ )
    targets[i++] = e10parity();
/*		  printf("%lx\n",targets[i-1]); */
}


unsigned long nodes_evaluated = 0;

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 < 2; x7 ++ )
  for( x8 = 0; x8 < 2; x8 ++ )
  for( x9 = 0; x9 < 2; x9 ++ )
  for( x10 = 0; x10 < 2; 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 ");
      return( print_indiv( print_indiv( ++buffer ) ) );

    case OR: printf( "OR ");
      return( print_indiv( print_indiv( ++buffer ) ) );

    case XOR: printf( "XOR ");
      return( print_indiv( print_indiv( ++buffer ) ) );

    case NAND: printf( "NAND ");
      return( print_indiv( print_indiv( ++buffer ) ) );

    case NOR: printf( "NOR ");
      return( print_indiv( print_indiv( ++buffer ) ) );

    case EQ: printf( "EQ ");
      return( print_indiv( print_indiv( ++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 tournament( int *fitness, int tsize )
{
  int best, fbest = -1, i, competitor;

  for ( i = 0; i < tsize; i ++ )
    {
      competitor = lrand48() % POPSIZE;
      if ( fitness[competitor] > fbest )
	{
	  fbest = fitness[competitor];
	  best = competitor;
	}
    }
  return( best );
}

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

  xo1start =  parent1 + lrand48() % len1;
  xo1end = traverse( xo1start );

  xo2start =  parent2 + lrand48() % len2;
  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 *mutation( char *parent, double pmut )
{
  char *mutsite;

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

  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;
	    }
	}
    }
  return( parent );
}


void print_parms()
{
  printf("\n\n----------------------------------\nGenerational Even 10 Parity problem  $Revision: 1.4 $\n----------------------------------\n");
  printf("SEED=%ld\nMAX_LEN=%d\nNTH_KILL=%d\nCAN_KILL_BEST=%d\nPOPSIZE=%d\nDEPTH=%d\nPMUT(not used)=%lg\nGENERATIONS=%d\nTSIZE=%d\n----------------------------------\n",
 seed, MAX_LEN, NTH_KILL, CAN_KILL_BEST, POPSIZE, DEPTH, PMUT, GENERATIONS, TSIZE );
}

/*for Memory Efficient Crossover Based on arXiv:2009.10460v1 */
int* forw; //valid 0..popsize-1
int* back; //valid 0..popsize-1, only maintained for chain2
int** children;  //popsize NULL or pointer to num_children int, -1 or 0..popsize-1
int chainhd1; //valid 0..popsize-1
int chainhd2; //valid 0..popsize-1
int* status;  // 0, 1 or 2

/*A.2.2    Functions used at the start of a new generation */
void append(const int s, int* Last, int* Chainhd) {
  const int last = *Last;
  if(last != -1) forw[last] = s;
  forw[s]    = -1;
  back[s]    = last;
  *Last      = s;
  if(*Chainhd == -1) *Chainhd = s;
}
void addchild(const int P, const int num_children, const int s) {
  //use linear search assuming num_children is small
  if(children[P]==NULL) {
    children[P] = (int*) malloc(num_children*sizeof(int));
    memset(children[P],0xff,num_children*sizeof(int));
  }
  int i=0;
  for(;i<num_children;i++){
    if(children[P][i] == -1){
      children[P][i]=s;
      break;
    }
  }
  assert(i < num_children);
  assert(children[P][i]==s);
}
/*A.2.4    Multi-threaded Functions used during a generation*/
int remchild(const int P, const int num_children, const int s, int* Last){
  /*use linear search assuming num_children is small*/
  int S = s;
  int nchild = 0;
  int ok = 0;
  int i=0;
  for(;i<num_children;i++){ /*remove only one child*/
    if(children[P][i] ==  S) {children[P][i] = -1; ok++; S=-2;}
    if(children[P][i] != -1) {*Last = children[P][i]; nchild++;}
  }
  assert(ok==1);
  if(nchild != 1) *Last -1; //be tidy
  return nchild;
}
void move21(const int active, const int s) {
  if(active == s)    return; //ignore
  //Even without threads use status in case called twice for the same guy
  if(status[s] != 2) return; //already processed

  /*remove s from chain2*/
  status[s] = 1;
  const int b = back[s];
  const int f = forw[s];
  if(chainhd2 == s) chainhd2 = f;
  forw[s] = -1;
  back[s] = -1;
  if(b != -1) forw[b] = f;
  if(f != -1) back[f] = b;

  /*Insert s at head of chain1, ok to ignore back*/
  forw[s] = chainhd1;
  chainhd1 = s;
}/*end move21*/

int main(int argc,char **argv)
{
  int gen = 0, indivs, parent1, parent2, best,  parent, len;
  int *fitness = (int *) calloc( POPSIZE, sizeof(int));
  int *newfit  = (int *) calloc( POPSIZE, sizeof(int));
  char **pop;
  char **newpop = (char **) calloc( POPSIZE, sizeof(char*));
  /*for Memory Efficient Crossover Based on arXiv:2009.10460v1 */
  int* mum       = (int *) calloc( POPSIZE, sizeof(int));
  int* dad       = (int *) calloc( POPSIZE, sizeof(int));
  int* num_children = (int *) calloc( POPSIZE, sizeof(int));
  forw      = (int *) calloc( POPSIZE, sizeof(int));
  back      = (int *) calloc( POPSIZE, sizeof(int));
  children = (int**) calloc( POPSIZE, sizeof(int*));
  status    = (int *) calloc( POPSIZE, sizeof(int)); /* 0, 1 or 2*/
  int used, max_used; /*accounting only*/

  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");
  }

  max_used = used = POPSIZE;
  for ( gen = 1; gen < GENERATIONS; gen ++ )
    {
      /*Decide parents of next generation*/
      memset(num_children,0,POPSIZE*sizeof(int));
      for ( indivs = 0; indivs < POPSIZE; indivs ++ )
	{
	  if (  lrand48() % 5 )
	    {/*crossover*/
	      parent1 = tournament( fitness, TSIZE );
	      parent2 = tournament( fitness, TSIZE );
	      num_children[parent1]++;
	      num_children[parent2]++;
	      mum[indivs] = parent1;
	      dad[indivs] = parent2;
	    }
	  else
	    {/*mutation*/
	      parent = tournament( fitness, TSIZE );
	      num_children[parent]++;
	      mum[indivs] = parent;
	      dad[indivs] = -1;
	    }
	}

      /*Classify parents by number of children and use two chains to
        order processing of the next generation so as to ensure no
        more than POPSIZE+2 trees are simultaneously needed*/
      int last1 = -1;
      int last2 = -1;
      chainhd1  = -1;
      chainhd2  = -1;
      memset(children, 0,POPSIZE*sizeof(int*));
      int s;
      for (s = 0; s < POPSIZE; s++ ) {
	const int mum_children = num_children[mum[s]];
	const int dad_children = (dad[s] != -1)? num_children[dad[s]] : 0;
	if(mum_children==1 || dad_children==1) {
	  append(s, &last1, &chainhd1);
	  status[s] = 1;
	} else {
	  append(s, &last2, &chainhd2);
	  status[s] = 2;
	}
	addchild(mum[s],mum_children,s);
	if(dad[s] != -1) addchild(dad[s],dad_children,s);
      }
      for (s = 0; s < POPSIZE; s++ ) {
	if(num_children[s]==0) {free( pop[s] ); pop[s] = NULL; used--;}
      }
      while(1) {
	/*Choose next individual from chain1 if any else chain2 else stop*/
	int i = -1;
	if(     chainhd1 != -1) {i = chainhd1;chainhd1 = forw[i];} 
	else if(chainhd2 != -1) {i = chainhd2;chainhd2 = forw[i];}
	if(i == -1) break;

	const int Mum = mum[i];
	const int Dad = dad[i];
	assert(pop[Mum]);
	assert(Dad == -1 || pop[Dad]);

	if(Dad != -1)
	  newpop[i] = crossover( pop[Mum],pop[Dad] );
	else {
	  const int parent = Mum;
	  len = traverse( pop[parent] ) - pop[parent];
	  newpop[i] = mutation( (char *) memcpy( (char *) malloc( len * sizeof(char)),
						  pop[parent], len ), PMUT );
	}
	used++;if(used>max_used){max_used = used;}

	/*Remove individual i from parents to do set.  
          If either parent has only one child left to be processed,
	  move the child (last1/last2) from chain 2 to chain 1.  
	  If either parent has no children to be processed delete it.*/
	int last1, last2;
	const int nchild1 = remchild(Mum,num_children[Mum],i,&last1);
	const int nchild2 = (Dad != -1)? remchild(Dad,num_children[Dad],i,&last2) : -1;
	if(nchild1 == 1) move21(i,last1);
	if(nchild2 == 1) move21(i,last2);
	if(nchild1 == 0) {
	  free( pop[Mum] );      pop[Mum]      = NULL; 
	  free( children[Mum] ); children[Mum] = NULL;
	  used--;
	}
	if(nchild2 == 0) {
	  free( pop[Dad] );      pop[Dad]      = NULL; 
	  free( children[Dad] ); children[Dad] = NULL;
	  used--;
	}

	newfit[i] = fitness_function( newpop[i] );
      }/*loop until both chains are empty*/

      for ( indivs = 0; indivs < POPSIZE; indivs ++ ) /*replace old pop with new*/
	{
	  assert(pop[indivs]==NULL);
	  pop[indivs]     = newpop[indivs];
	  fitness[indivs] = newfit[indivs];
	}

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

	return EXIT_SUCCESS;
}