/*WBL 30 March 2018 */
#define version "$Revision: 1.48 $"

//compile gcc -DNDEBUG -O2 px_nk.c -lm

/*Modification:
WBL 23 Apr 2018 clean for release, incorporate pch.h so single file, enable NDEBUG
                allow more K values
WBL 22 Apr 2018 try inspiration from frchicanog/EfficientHillClimbers-master/src/main/java/neo/landscape/theory/apps/pseudoboolean/px/PartitionCrossover.java  Jul 16  2017
WBL 16 Apr 2018 revert r1.34 (again 2) call freduced once
WBL  9 Apr 2018 r1.18 still not ok, try GA2 to do tintos_2015_foga algorithm 3
WBL  2 Apr 2018 Add lsfi like tintos_2015_foga algorithm 2
WBL  1 Apr 2018 Add pc and other changes to GA to make like tintos_2015_foga algorithm 2
   Add px
   Results in Darrell Whitley's PPSN-2016 tutorial in tintos_2015_foga.pdf table 5
WBL 31 Mar 2018 ensure each locus is independent
*/

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include <unistd.h>
#include <string.h>
#include <float.h>

//pick up Park_miller random numbers gp/solexa_3623/source gpquick_seed
int gpquick_seed;
//#include "pch.h"
// pch.h   Adapted from GPQUICK-2.1
// W. Langdon cs.ucl.ac.uk 5 May 1994 (21 Feb 95 make it inline)
// Use Park-Miller random numbers

// QPQUICK
// Standard header files for compilers supporting Pre-Compiled Headers

//WBL 31 Mar 18  for gcc 4.8.5 NB C not C++
//WBL 19 Oct 03  Port to Alpha gcc 3.3
//WBL 21 Dec 02  Port to Microsoft Visual C++

#include <limits.h>

//////////  Random number functions cause portability problems.  Resolve here.
// Define	rand_0to1()  (random float from 0 to 1)
//			rnd(x) (random integer from 0 to x-1)
//			rndomize() (initialize with a time based seed) 			
// park-miller.cc   Park-Miller random number generator
// W. Langdon cs.ucl.ac.uk 5 May 1994

#ifndef LONG_MAX
error LONG_MAX missing
#endif

inline int intrnd (int* seed) // 1<=seed<=m
{
#if LONG_MAX > (16807*2147483647)
long int const seed_= *seed;
long int const a    = 16807;      //ie 7**5
int const m    = 2147483647; //ie 2**31-1
	*seed = (seed_ * a)%m;
	return *seed;
#else
double const a    = 16807;      //ie 7**5
double const m    = 2147483647; //ie 2**31-1

	double temp = *seed * a;
	*seed = (int) (temp - m * floor ( temp / m ));
	//cout<<"seed "<<seed<<endl;
	return *seed;
#endif
}

	extern int gpquick_seed;

        inline float rand_0to1()
                { return (float)(intrnd(&gpquick_seed)) / 2147483647.0; }

        inline double drand_0to1()
                { return (double)(intrnd(&gpquick_seed)) / 2147483647.0; }

	inline int rnd(int __num)
		{ return (intrnd(&gpquick_seed) % __num); }
	inline void rndomize(void)
		{ do {
			gpquick_seed = (unsigned) time(NULL);
		     } while (gpquick_seed <= 0);
                }
//end pch.h
inline int myrnd(const int __num, int* prng_NK2) {
  const int save = gpquick_seed;
  gpquick_seed = *prng_NK2;
  const float ans = rnd(__num);
  *prng_NK2 = gpquick_seed;
  gpquick_seed = save;
  return ans;
}

inline float myrand_0to1(int* prng_NK1) {
  const int save = gpquick_seed;
  gpquick_seed = *prng_NK1;
  const float ans = rand_0to1();
  *prng_NK1 = gpquick_seed;
  gpquick_seed = save;
  return ans;
}


//set once globals
int K;
int NK_rand;
//since K is small we might as well initialise f before we start
int NK_fit_stride;
float* NK_fit;
int* neighbour;

inline float NKi(const int i, const int N, const char chrome[N]){ //K etc global
  assert(i>=0 && i<N);
  const int L0 = chrome[i];
  assert(L0==0||L0==1);
  unsigned int L = L0;
  int k;
  for(k=0;k<K;k++) {
    const int n = neighbour[i*K+k];
    assert(n>=0 && n<N);
    const unsigned int l = chrome[n];
    assert(l==0||l==1);
    L = L + (l<<(k+1));
  }
  assert(L<NK_fit_stride);
  return NK_fit[NK_fit_stride*i+L];
}

int* neighbour2_forw;
int* neighbour2_head;
int* neighbour2_loc;

int NKn_skip = 0;

inline float NKn(const int loc, const int N, const char chrome[N], char done[N]){
  float ans = 0;
  int sanity=0;
  int i;
  for(i=neighbour2_head[loc]; i != -1; i=neighbour2_forw[i]) {
    assert(i>=0 && i<N*(K+1));
    const int f = neighbour2_loc[i];
    if(done==NULL || done[f]==0) { //only do each f once
      ans += NKi(f,N,chrome);
      if(done!=NULL) done[f] = 1;
    } else NKn_skip++;
    assert(sanity++<=N);
  }
  assert(sanity>0);
  return ans;
}

inline float NK(const int N, const char chrome[N]){
  int i;
  float ans = 0;
  for(i=0;i<N;i++) { ans += NKi(i,N,chrome); }
//printf("NK(%d,chrome) %f\n",N,ans);

return ans;
}

int nstat2;
double sum_nimp;
double sum2_nimp;

//bit flip local search (fastest ascent hill climbing)
//Doing search in deterministic order (tintos_2015_foga does not defined the order)
void lsfi(const int N, char chrome[N]){
  int count=0;
  int sanity=0;
  int cont;
  do {
    cont=0;
    int loc;
#ifdef NDEBUG
    for(loc=0;loc<N;loc++){
      const float f0 = NKn(loc,N,chrome,NULL);
      chrome[loc] = (~chrome[loc]) & 1;
      const float f1 = NKn(loc,N,chrome,NULL);
      //unless f1 improved restore original allel
      if(f0>=f1) chrome[loc] = (~chrome[loc]) & 1;
      else {cont=1; count++;}
#else //NDEBUG
    for(loc=0;loc<N;loc++){
      const float f0 = NKn(loc,N,chrome,NULL);
      //printf("lsfi(%d,chrome) loc%d f0 NKn %f\n",N,loc,f0);
      const float f0_= NK(N,chrome);
      //printf("lsfi(%d,chrome) loc%d f0_ NK %f\n",N,loc,f0_);
      chrome[loc] = (~chrome[loc]) & 1;
      const float f1 = NKn(loc,N,chrome,NULL);
      //printf("lsfi(%d,chrome) loc%d f1 NKn %f\n",N,loc,f1);
      const float f1_= NK(N,chrome);
      //printf("lsfi(%d,chrome) loc%d f1_ NK %f\n",N,loc,f1_);
      //unless f1 improved restore original allel
      if(f0>=f1) chrome[loc] = (~chrome[loc]) & 1;
      else {cont=1; count++;}

      const float epsilon = 2*FLT_EPSILON*N;
      //printf("lsfi loc%d %f %f %f %f epsilon %f\n",loc,f0,f0_,f1,f1_,epsilon);
      if(f0 <f1) assert(f0_ <f1_+epsilon);
      if(f0==f1) assert(f0_==f1_);
      if(f0> f1) assert(f0_> f1_-epsilon);
#endif
    }
    assert(sanity++<=N*N);
  } while(cont);
  nstat2++; sum_nimp += count; sum2_nimp += count*count;
}

void twopoint(const int N, const char p1[N], const char p2[N], char child[N]){
  int i;
  const int x = rnd(N);
  const int y = rnd(N);
  for(i=0;i<N;i++) child[i] = ((i<=x && i<=y)||(i>x && i>y))? p1[i] : p2[i];
}

void uniform(const int N, const char p1[N], const char p2[N], char child[N]){
  int i;
  for(i=0;i<N;i++) child[i] = rnd(2)? p1[i] : p2[i];
}


int npartitions = 0;
int* partition_head;
int* partition_forw;
int* partition_back;

void init_partitions(const int N) {
  partition_head = malloc((N+1)*sizeof(int));
  partition_forw = malloc( N   *sizeof(int));
  partition_back = malloc( N   *sizeof(int));
}

int newpartition(const int N, const int i) {
  assert(i>=0 && i< N);
  assert(npartitions < N);
  npartitions++;
  const int p = npartitions;
  partition_back[i] = p;
  partition_forw[i] = -1;
  partition_head[p] = i;
  //printf("newpartition(%d,%d) set npartitions %d\n",N,i,npartitions);
  return p;
}
void addtopartition(const int N, const int i, const int p){
  assert(i>=0 && i< N);
  assert(p> 0 && p<=npartitions);
  partition_back[i] = p;
  partition_forw[i] = partition_head[p];
  partition_head[p] = i;
}

#ifndef NDEBUG
void sanity_chain(const int N, const int p) {
  assert(p> 0 && p<=npartitions);
  if(partition_head[p] == -1) return;
  int sanity=0;
  int loc;
  for(loc=partition_head[p]; loc != -1; loc=partition_forw[loc]) {
    assert(partition_back[loc]==p);
    assert(sanity++<=N);
  }
}
#endif
void joinpartitions(const int N, const int p, const int q){
  assert(p> 0 && p<=npartitions);
  assert(q> 0 && q<=npartitions);
  assert(p != q);
  //add to q(last) as perhaps less to update
  int I=0;
  int sanity=0;
  int loc;
  for(loc=partition_head[p]; loc != -1; loc=partition_forw[loc]) {
    assert(partition_back[loc]==p);
    partition_back[loc] = q;I=loc;
    assert(sanity++<=N);
  }
  partition_forw[I] = partition_head[q];
  partition_head[q] = partition_head[p];
  partition_head[p] = -1; //dont try and reuse chain head
#ifndef NDEBUG
  sanity_chain(N,p);
  sanity_chain(N,q);
#endif
}

//Francis uses a Java set here but the list should be short or we could use done as in NKn
void addto(const int in, int* ntodo, int todo[1 + *ntodo]) {
  //printf("addto(%d,%d ",in,*ntodo);fflush(stdout);
  int i;
  for(i=0;i < *ntodo;i++) {if(todo[i]==in) return;}
  todo[*ntodo]=in;
  *ntodo = 1 + *ntodo;
}
//chose best by calculating fitness contribution of partition
int freduced(const int N, const char p[N], const int partition){
  int todo[N];
  int ntodo = 0;
  int sanity=0;
  int loc;
  for(loc=partition; loc != -1; loc=partition_forw[loc]) {
    assert(loc>=0 && loc<N);
    int i;
    for(i=neighbour2_head[loc]; i != -1; i=neighbour2_forw[i]) {
      assert(i>=0 && i<N*(K+1));
      const int f = neighbour2_loc[i];
      addto(f, &ntodo, todo);
      assert(sanity++<=N*N);
    }
    assert(sanity>0);
  }
  char chrome[N];
  memcpy(chrome,p,N); //for simplicity ensure all loci are defined
  float a1 = 0;
  int i;
  for(i=0;i<ntodo;i++){ a1 += NKi(todo[i],N,chrome); }

  //load all of partition then do fitness evaluations
  for(loc=partition; loc != -1; loc=partition_forw[loc]) {
    chrome[loc] = (~p[loc]) & 1;
  }
  float a2 = 0;
  for(i=0;i<ntodo;i++){ a2 += NKi(todo[i],N,chrome); }
  return            (a1 > a2);
}

int nstat;
double sum_npar;
double sum2_npar;
void partition_crossover(const int N, const char p1[N], const char p2[N], char child[N]){
  //based on description in tintos_2015_foga.pdf
#ifndef NDEBUG
  memset(child,123,N);//for sanity check
#endif
  npartitions = 0;
  memset(partition_head,0,(N+1)*sizeof(int));
  memset(partition_forw,0, N   *sizeof(int));
  memset(partition_back,0, N   *sizeof(int));

  int i;
  for(i=0;i<N;i++) {
    if(p1[i] == p2[i]) child[i] = p1[i];
    else {//is loc i already in a partition?
      const int p = partition_back[i];
      const int firstp = (p>0)? p : newpartition(N,i);
      int k;
      for(k=0;k<K;k++){
	const int n = neighbour[i*K+k];
	if(p1[n] == p2[n]) continue;
	const int p = partition_back[n];
	if(p==0) addtopartition(N,n,firstp);
	else if(p != firstp) joinpartitions(N,p,firstp);
      }
    }
  }

  int count=0;
  int sanity=0;
  int q;
  for(q=1;q<=npartitions;q++) {
    const int worse = freduced(N,p1,partition_head[q]);
    int loc;
    for(loc=partition_head[q]; loc != -1; loc=partition_forw[loc]) {
#ifdef NDEBUG
      child[loc] = worse? p1[loc] : p2[loc];
#else
      assert(child[loc] == 123);
      child[loc] = worse? p1[loc]|(2+10*q) : p2[loc]|(4+10*q);
#endif
      count++;
      assert(sanity++<=N);
    }
  }
  nstat++; sum_npar += count; sum2_npar += count*count;
#ifndef NDEBUG
  /*
  printf("px i:    ");
  for(i=0;i<N;i++) {printf("%2d ",i%100);}
  printf("\n");
  printf("px p1:   ");
  for(i=0;i<N;i++) {printf("%2d ",p1[i]);}
  printf("\n");
  printf("px p2:   ");
  for(i=0;i<N;i++) {printf("%2d ",p2[i]);}
  printf("\n");
  printf("px debug:");
  */
  for(i=0;i<N;i++) {
    //printf("%2u ",(unsigned char)child[i]);
    child[i] = child[i] & 1; //clean up debug flags
  }
  //printf("\n");
  //printf("px child:");
  //for(i=0;i<N;i++) {printf("%2d ",child[i]);}
  //printf("\n");
  for(i=0;i<N;i++) assert(child[i]==0 || child[i]==1); //make sure have not missed any loci
#endif
}

void crossover(const int XOtype, const int N, const char p1[N], const char p2[N], char child[N]){
  switch(XOtype) {
  case 0: return twopoint(N,p1,p2,child);
  case 1: return uniform( N,p1,p2,child);
  case 2: return partition_crossover(N,p1,p2,child);
  default : exit(98);
  }
}

int selection(const int ts, const int pop_size, const float f[pop_size]) {
  float bestfit = -1.0f; int bestguy;
  int i;
  for(i=0;i<ts;i++) {
    const int a = rnd(pop_size);
    if(f[a]>bestfit) {bestfit = f[a]; bestguy=a;}
  }
  return bestguy;
}
//call only once
void GA(const int N, const int pop_size, const int maxgen, const int XOtype, 
	const float pc, const float pm, const int ts, const int seed_GA){
  int i,g;
  const int save = gpquick_seed;
  gpquick_seed = seed_GA;

  //try to reproduce Hybrid GA p144 FOGA 2015
  //Koza appendix F shows we only need pop_size+ts not 2*pos_size
  char* pop  = calloc(pop_size,N);
  char* pop2 = calloc(pop_size,N); 

  for(i=0;i<pop_size*N;i++) pop[i] = rnd(2);
  for(i=0;i<pop_size;i++) lsfi(N,&pop[i*N]);

  float bestfit = -1.0f; int bestguy;
  for(g=0;g<=maxgen;g++) {
    int i;
    float f[N];
    for(i=0;i<pop_size;i++) {
      f[i] = NK(N,&pop[i*N]);
      if(f[i]>bestfit) {bestfit = f[i]; bestguy=i;}
    }
    printf("Best of %d generation %3d %2d %f\n",pop_size,g,bestguy,bestfit);
    if(g==0) goto do_stats;
    memcpy(&pop2[0],&pop[bestguy],N); //elitist 1
    //every ten generations reset 90% of pop then hill climb
    //including the best solution found in the previous generation
    //since our crossover etc are not recalculating bestguy do reset on pop before
    //crossover and mutation. The effect should be as pseudo code on Q in alg 2 p144
    if(g%10==0){
      for(i=0;i<pop_size;i++) {
	if(i!=bestguy && i>=pop_size/10) {
	  int loc;
	  for(loc=0;loc<N;loc++) pop[i*N+loc] = rnd(2);
	  //printf("%d %d %d %d %d %d \n",i*N+loc,pop_size*N,i,N,loc,pop_size);
	  assert(i*N+loc<=pop_size*N);
	}
	assert(i*N<pop_size*N);
	lsfi(N,&pop[i*N]);
      }
    }

    for(i=1;i<pop_size;i++) {
      const int p1 = selection(ts,pop_size,f);
      const int p2 = selection(ts,pop_size,f);
      if(rand_0to1() < pc)
	crossover(XOtype,N,&pop[p1*N],&pop[p2*N],&pop2[i*N]);
      else {
	int loc;
	for(loc=0;loc<N;loc++)
	  pop2[i*N+loc] = (rand_0to1() < pm)? (~pop[p1*N+loc]) & 1 : pop[p1*N+loc];
      }
    }
    //next generation
    memcpy(pop,pop2,pop_size*N);

  do_stats:
    if(nstat) {
      const double mean = sum_npar/nstat;
      const double var  = (sum2_npar - mean*mean)/nstat;
      printf("Partitions generation %3d %6.2f SD %6.2f   PX %2d\n",
	     g,mean,(var>0)? sqrt(var):0.0, nstat);
      nstat=0; sum2_npar=sum_npar=0.0;
    }
    if(nstat2) {
      const double mean = sum_nimp/nstat2;
      const double var  = (sum2_nimp - mean*mean)/nstat2;
      printf("lsfi       generation %3d %6.2f SD %6.2f used %2d\n",
	     g,mean,(var>0)? sqrt(var):0.0, nstat2);
      nstat2=0; sum2_nimp=sum_nimp=0.0;
    }
    if(NKn_skip) {
      printf("NKn_skip   generation %3d %6d\n",g,NKn_skip);
      NKn_skip = 0;
    }
  }
  //gpquick_seed = save; call only once
}

int local_opt(const int N, const char in[N]){ //for stats only
  char chrome[N]; //buffer for safety
  memcpy(chrome,in,N);
  int loc;
  for(loc=0;loc<N;loc++){//based on lsfi
#ifdef NDEBUG
    const float f0 = NKn(loc,N,chrome,NULL);
    chrome[loc] = (~chrome[loc]) & 1;
    const float f1 = NKn(loc,N,chrome,NULL);
    chrome[loc] = (~chrome[loc]) & 1;
#else
    const float f0 = NKn(loc,N,chrome,NULL);
    const float f0_= NK(N,chrome);
    chrome[loc] = (~chrome[loc]) & 1;
    const float f1 = NKn(loc,N,chrome,NULL);
    const float f1_= NK(N,chrome);
    chrome[loc] = (~chrome[loc]) & 1;

    const float epsilon = 2*FLT_EPSILON*N;
    if(f0 <f1) assert(f0_ <f1_+epsilon);
    if(f0==f1) assert(f0_==f1_);
    if(f0> f1) assert(f0_> f1_-epsilon);
#endif
    if(f1>f0) return 0;
  }
  return 1; //no single bit flip improvement possible
}

//call only once
void GA2(const int N, const int pop_size, const int maxgen, const int XOtype, 
	const int seed_GA){
  int i,g;
  const int save = gpquick_seed;
  gpquick_seed = seed_GA;

  //try to reproduce Alg 3 "Recombination of Local Optima" p146 FOGA 2015
  char* pop  = calloc(pop_size,N);
  char* bestchrome = calloc(N,1); //actually only need best of pop2
  float bestfit = -1.0f;

  for(g=0;g<maxgen;g++) {
    //int nstat3;
    //double sum_nlopt;
    //double sum2_nlopt;

    int i;
    //new random population every generation
    for(i=0;i<pop_size*N;i++) pop[i] = rnd(2);
    for(i=0;i<pop_size;i++) lsfi(N,&pop[i*N]);

    {float bestfit2 = -1.0f;
    //find best of new random population and old population(if any)
    for(i=0;i<pop_size;i++) {
      const float f = NK(N,&pop[i*N]);
      if(f>bestfit) {bestfit  = f; memcpy(bestchrome,&pop[i*N],N);}
      if(f>bestfit2) bestfit2 = f; //just for stats
    }
    printf("Best of %d generation %3d %f overall %f\n",
	   pop_size,g,bestfit2,bestfit);
    }
    int count=0;
    int count2=0;
    //crossover everyone (possibly including self) with best so far
    char bestchromeQ[N];
    float bestfitQ = -1.0f;
    for(i=1;i<pop_size;i++) {//cross with best now, not best found inloop
      char temp[N];
      assert(local_opt(N,bestchrome));
      assert(local_opt(N,&pop[i*N]));
      crossover(XOtype,N,bestchrome,&pop[i*N],temp);
      const float f = NK(N,temp);
      if(f>bestfitQ) {
	const l =  (memcmp(temp,bestchrome,N)==0) ? 2 :
	          ((memcmp(temp,&pop[i*N], N)==0) ? 3 : 
	          local_opt(N,temp));
	printf("New bestxo generation %3d %d %f %f local_opt %d\n",
	       g,i,f,bestfitQ,l);
	bestfitQ = f; 
	memcpy(bestchromeQ,temp,N); //actually only used if bestQ>best ever
      }
      if(memcmp(temp,bestchrome,N)!=0 && memcmp(temp,&pop[i*N],N)!=0) {
	if(local_opt(N,temp)) count++;
	count2++;
      }
    }
    //nstat3++; sum_nlopt += count; sum2_nlopt += count*count;

    printf("Best of %d generation %3d after crossover %f\n",pop_size,g,bestfitQ);
    if(bestfitQ > bestfit){
      bestfit = bestfitQ;
      memcpy(bestchrome,bestchromeQ,N);
      lsfi(N,bestchrome);
      bestfit = NK(N,bestchrome);
      printf("Best of %d generation %3d after lsfi %f\n",pop_size,g,bestfit);
    }

    if(nstat) {
      const double mean = sum_npar/nstat;
      const double var  = (sum2_npar - mean*mean)/nstat;
      printf("Partitions generation %3d %6.2f SD %6.2f   PX %2d\n",
	     g,mean,(var>0)? sqrt(var):0.0, nstat);
      nstat=0; sum2_npar=sum_npar=0.0;
    }
    if(nstat2) {
      const double mean = sum_nimp/nstat2;
      const double var  = (sum2_nimp - mean*mean)/nstat2;
      printf("lsfi       generation %3d %6.2f SD %6.2f used %2d\n",
	     g,mean,(var>0)? sqrt(var):0.0, nstat2);
      nstat2=0; sum2_nimp=sum_nimp=0.0;
    }
    if(NKn_skip) {
      printf("NKn_skip   generation %3d %6d\n",g,NKn_skip);
      NKn_skip = 0;
    }
    //assert(nstat3==1);
    if(1){//nstat3) {
      const double mean = count; //sum_nlopt/nstat3;
      const double var  = 0;     //(sum2_nlopt - mean*mean)/nstat3;
      printf("Local opt  generation %3d %6.2f SD %6.2f ndiff %2d\n",
	     g,mean,(var>0)? sqrt(var):0.0, count2);
      //nstat3=0; sum2_nlopt=sum_nlopt=0.0;
    }
  }
  //gpquick_seed = save; call only once
}

/*--------------------------------------------------------------------------*/

//http://www.ppsn2016.org/conference/wp-content/uploads/2015/07/GrayBox-PPSN-16-Whitley.pdf

//try to reproduce table on slide 32
int main(int argc, char *argv[]){
  int i,j;
  i=1;
  const int N        = (argc>i && argv[i][0])? atoi(argv[i]) : 500;    i++;
            K        = (argc>i && argv[i][0])? atoi(argv[i]) : 1;      i++;
            NK_rand  = (argc>i && argv[i][0])? atoi(argv[i]) : 0;      i++; //Adjacent
  const int GA_XO    = (argc>i && argv[i][0])? atoi(argv[i]) : 0;      i++; //2 point
  const float pc     = (argc>i && argv[i][0])? atof(argv[i]) : 0.6;    i++;
  const float pm     = (argc>i && argv[i][0])? atof(argv[i]) : 3.0/N;  i++;
  const int ts       = (argc>i && argv[i][0])? atoi(argv[i]) : 5;      i++;
  const int seed_NK1 = (argc>i && argv[i][0])? atoi(argv[i]) : 301547; i++;
  const int seed_NK2 = (argc>i && argv[i][0])? atoi(argv[i]) : 503015; i++;
  const int seed_GA  = (argc>i && argv[i][0])? atoi(argv[i]) : 154830; i++;

  if(pc>=1.0) printf("GA2 pop20 gens20 ");
  else        printf("GA pop50 gens100 ");
  printf("on NK landscape %s ",version); fflush(stdout);
#ifdef NDEBUG
  printf("NDEBUG ");
#endif
  printf("N%d,K%d %d GA_XO=%d %s ",N,K,NK_rand,GA_XO,
	 (GA_XO==0)? "2pt" : ((GA_XO==1)? "ux" : "px"));
  if(pc>=1.0) printf("only crossover");
  else        printf("pc %g pm %g",pc,pm);
  printf(", prng seeds: %d %d %d\n",seed_NK1,seed_NK2,seed_GA);
  if(ts>=100) {
    fprintf(stderr, "Bad command line? tournament size %d stopping\n",ts);
    return 1;
  }

//since we expect K to be small we might as well initialise f before we start
  NK_fit_stride = 1<<(K+1);
  NK_fit = calloc(NK_fit_stride*N,sizeof(float));
  int prng_NK1 = seed_NK1;
  for(i=0;i<N;i++) {
  for(j=0;j<NK_fit_stride;j++) {
    const int k = NK_fit_stride*i+j;
    NK_fit[k]=myrand_0to1(&prng_NK1);
    //printf("NK_fit[%d] %f\n",k,NK_fit[k]);
  }}
  int prng_NK2 = seed_NK2;
  neighbour = calloc(N*K,sizeof(int));
  neighbour2_head = malloc(N      *sizeof(int));
  neighbour2_loc  = calloc(N*(K+1),sizeof(int));
  neighbour2_forw = calloc(N*(K+1),sizeof(int));
  for(i=0;i<N;i++) neighbour2_head[i] = i; //self
  for(i=0;i<N;i++) neighbour2_loc[i]  = i; //self
  for(i=0;i<N;i++) neighbour2_forw[i] = -1; //end
  int nskip = 0;
  int nfn = N;
  for(i=0;i<N;i++) {
  for(j=0;j<K;j++) {
    const int loc = 
    neighbour[i*K+j] = (NK_rand==0)? (i+1+j)%N : (i+1+myrnd(N-1,&prng_NK2))%N;
    //due to overlap may sometimes have less than K distinct randomly chosen neighbours
    assert(loc>=0 && loc<N);
    int found = 0; //check not already a neighbour
    int sanity=0;
    int m;
    for(m=neighbour2_head[loc]; found==0 && m != -1; m=neighbour2_forw[m]) {
      if(neighbour2_loc[m] == i) {
	assert(NK_rand!=0);
	found = 1;
      }
      assert(sanity++<=N);
    }
    //if(found) {nskip++; break;} something odd with gcc version 4.8.5 break use if-else instead
    if(found) nskip++;
    else {
    assert(nfn<N*(K+1));
    assert(neighbour2_loc[nfn] == 0);
    neighbour2_forw[nfn] = neighbour2_head[loc];
    neighbour2_head[loc] = nfn;
    neighbour2_loc[nfn]  = i;
    nfn++;
    }
  }
  if(nfn+nskip!=N+(i+1)*K) printf("i %d j %d nfn %d nskip %d\n",i,j,nfn,nskip);
  assert(nfn+nskip==N+(i+1)*K);
  }
  printf("nfn %d nskip %d\n",nfn,nskip);
  assert(nfn+nskip==N*(K+1));

  init_partitions(N);
  if(pc>=1.0)
    GA2(N,20,20,GA_XO,seed_GA); //NK landscape
  else 
  GA(N,50,100,GA_XO,pc,pm,ts,seed_GA); //NK landscape, tournament selection
  return 0;
}