/*
 Chad Brewbaker
5-24-2002
seedMixer.c

   
This program searches the space of golomb seeds by making a greedy construction with a seed, and then using random elements out of that construction as new seeds


*/
#include <stdlib.h>
#include <stdio.h>
#include <time.h>

//Play with these
#define RULER_SIZE 8 //Length of the completed ruler
#define SEED_SIZE 4 //Size of the seed
#define GENERATIONS 1000 //Number of random starts
#define LOOPS 100 //How long each random seed is evolved


//Mess with these at your own peril. Quick hack aritifacts.
#define MAX_NUM 90//Maximum size of a seed number
#define RULER_NUM 100//population size
#define PXOVER .80
#define PMUTATE .1
#define LUT_SIZE 500//look up table size
#define PENALTY_MULT 200


int seeds[RULER_NUM][SEED_SIZE];//population storage
int newSeeds[RULER_NUM][SEED_SIZE];//nextPopulation storage
double seedFit[RULER_NUM];//a seed's fitness
double rSeedFit[RULER_NUM];//relative fitness
double cSeedFit[RULER_NUM];//cumulative fitness
int worstFit;//worst score for each generation
//stores the superseed
int elete[SEED_SIZE];
int superScore;

int usedTable[MAX_NUM*MAX_NUM];//array used for sorting/searching
int rulerMap[RULER_SIZE];//a copy of the completed ruler

//circular buffer
//-New element:kill off the oldest and replace with new
//-Old element:swap oldest with the current and incrament buffer
//
//
int lookUpTable[LUT_SIZE][SEED_SIZE];
double lookUpTableFit[LUT_SIZE];
int oldest=0;//oldest member of the lookup table






void exchange(int i, int j)
{ 
  int temp;
  temp=usedTable[i];
  usedTable[i]=usedTable[j];
  usedTable[j]=temp;
}

int partition(int p,int r)
{
  int i,j;
  int x;
  x=usedTable[r];
  i=p-1;
  for(j=p; j<=r-1; j++)
    {
      if(usedTable[j]<=x)
        {
          i+=1;
          exchange(i,j);
        }
    }
  exchange(i+1,r);
  return i+1;
}

void quicksort(int p, int r)
{
  int q;
  if(p<r)
    {
      q=partition(p,r);
      quicksort(p,q-1);
      quicksort(q+1,r);
    }
}

int ABS(int a, int b)
{
  int temp = a-b;
  if(temp<0)
    temp = temp*-1;
  return temp;
}

int compare (const void * a, const void * b)
{
  return ( *(int*)a - *(int*)b );
}




int checkList(int size)
{
  int i,j,k;
  int score=0;
  //printf("\nSIZE:%d",size);
  for(i=0; i<size-1;i++)
    {
      if(usedTable[i]==usedTable[i+1])
	score++; 
    }
  return score;
}

int badSeed(int seedNum)
{
  int i,j,k;
  int isBad=0;
  int tableSize=0;
  //see if seed is an acceptable Golomb Ruler
  for(i=0;i<SEED_SIZE;i++)
    {
      usedTable[tableSize]=seeds[seedNum][i];
      tableSize++;
    }
  for(i=0;i<SEED_SIZE-1;i++)
    {
      for(j=i+1;j<SEED_SIZE;j++)
	{
	  usedTable[tableSize]=ABS(seeds[seedNum][i],seeds[seedNum][j]);
	    tableSize++;
	}
    }
  quicksort(0,tableSize-1);
  isBad=checkList(tableSize);
  return isBad;
}


void init(int seedNum)
{
  int i,j,k;
  int repeats;
  
  //initialize a seeds
  
  for(i=0;i<SEED_SIZE; i++)
    {
      do
	{
	  repeats=0;
	  seeds[seedNum][i]=(rand()%MAX_NUM)+1;
	  for(j=0;j<i;j++)
	    {
	      //make sure the seeds don't repeat
	      if(seeds[seedNum][j]==seeds[seedNum][i])
		{
		  repeats++;
		  break;
		}
	    }
	}
      while(repeats);
    }
}


void reInit(int seedNum)
{
  int i,j,k;
  int repeats;
  
  //re-initialize a seed
  
  for(i=0;i<SEED_SIZE; i++)
    {
      do
	{
	  repeats=0;
	  seeds[seedNum][i]=rulerMap[(rand()%(RULER_SIZE-1))+1];
	  for(j=0;j<i;j++)
	    {
	      //make sure the seeds don't repeat
	      if(seeds[seedNum][j]==seeds[seedNum][i])
		{
		  repeats++;
		  break;
		}
	    }
	}
      while(repeats);
    }
}

int min(int mapSize)
{
  int i;
  int best;
  best =9999999;
  for(i=0;i<mapSize;i++)
    {
      if(rulerMap[i]<best)
	best=rulerMap[i];
    }
  return best;
}

int max(int mapSize)
{
  int i;
  int best;
  best =-999999;
  for(i=0;i<mapSize;i++)
    {
      if(rulerMap[i]>best)
	best=rulerMap[i];
    }
  return best;
}

int growRuler(int rulerNum, int printMe)
{
  int i,j,k,l;  
  int temp;
  int tableSize=0;
  int oldTableSize=0;
  int good;
  int mapSize;
  int score;
  int outer;
  int maximal;
  int *pItem;
  int sameCount;
  int seedDifs;
  int properLength;
  int penalty;
  int tempTableSize;


  if(printMe)
    {
      printf("\nSEED:{");
      for(i=0;i<SEED_SIZE;i++)
	printf("%d, ",seeds[rulerNum][i]);
      printf("}");
    }

  rulerMap[0]= maximal=0;
  for(i=0;i<SEED_SIZE;i++)
    {
      rulerMap[i+1]=seeds[rulerNum][i];
      if(rulerMap[i+1]>maximal)
	maximal=rulerMap[i+1];
    }
  
  mapSize=SEED_SIZE+1;
  
  
  //add unique seed diferences to usedTable
  tableSize=0;
  for(i=0;i<mapSize-1;i++)
    {
      for(j=i+1;j<mapSize;j++)
	{
	  temp=ABS(rulerMap[i],rulerMap[j]);
	 
	  good=1;
	  for(k=0;k<tableSize;k++)
	    {
	      if(temp == usedTable[k])
		{
		  good=0;
		  break;
		}
	    }
	  if(good==1)
	    {
	      usedTable[tableSize]=temp;
	      tableSize++;
	    }
	}
    }
 

  seedDifs=tableSize;

  if(SEED_SIZE<1)
    properLength =0;
  else
    {
      i = SEED_SIZE;
      //Gauss Rocks... Sum of 1,2,...,RULER_LENGTH-1
      properLength= (i*(i+1))/2;
    }
  penalty= (properLength-tableSize)*PENALTY_MULT;
  
  //greedy search for new marks
  for(i=1;i<99999;i++)
    {

      //when a long enough ruler is made      
      if(mapSize==RULER_SIZE)
	break;
      good=1;
      
      //O(n^3)
      
      tempTableSize=tableSize;
      for(j=0;j<mapSize;j++)
	{
	  temp = ABS(i,rulerMap[j]);
	  for(k=0;k<tempTableSize;k++)
	    {
	      if(usedTable[k]==temp)
		{
		  //printf("\nBAD:%d @ %d",temp,i);
		  good=0;
		  //double break
		  k=tempTableSize+1;
		  j=mapSize+1;
		}
	    }
	  if(good)
	    {
	      usedTable[tempTableSize]=temp;
	      tempTableSize++;
	    }
	}
      
      if(good==0)
	continue; 
      //Good? Write the new mark
      tableSize=tempTableSize;
      rulerMap[mapSize]=i;
      mapSize++;
    }
  /*  
      printf("\nUT: ");
      for(i=0;i<tableSize;i++)
      printf("%d ",usedTable[i]);
  */
 //Now penalize the score for every repeated distance in the seed
  

  //find the length of the ruler
  score =rulerMap[0];
  for(i=1;i<mapSize;i++)
    {
      if(score < rulerMap[i])
	score=rulerMap[i];
    }




  if(printMe)
    {
      printf("\nRULER_SIZE:%d , MAX:%d, PENALTY:%d\n{",RULER_SIZE,score,penalty);
      for(i=0;i<RULER_SIZE;i++)
	printf("%d,",rulerMap[i]);
      printf("}\n");
    }

  //bubble sort hack to make pretty
  for(i=0;i<RULER_SIZE-1;i++)
    {
      for(j=i+1;j<RULER_SIZE;j++)
	{
	  if(rulerMap[j-1]>rulerMap[j])
	    {
	      temp = rulerMap[j];
	      rulerMap[j]=rulerMap[j-1];
	      rulerMap[j-1]=temp;
	    }
	}
    }


  
  return score+penalty;
}


void crossover(void)
{
  int i,j;
  int count=0;
  int size,position;
  int store,temp;
  double x;
  for(i=0;i<RULER_NUM;i++)
    {
      x=rand()%1000/1000.0;
      if(x<PXOVER)
	{
	  count++;
	  if(count%2==0)
	    {
	      size=rand()%SEED_SIZE;
	      position=rand()%SEED_SIZE;
	      for(j=0;j<size;j++)
		{
		  position=(position+i)%SEED_SIZE;
		  temp=seeds[i][position];
		  seeds[i][position]=seeds[store][position];
		  seeds[store][position]=temp;
		}
	    }
	  else 
	    store=i;  
	}
    }
}

void mutate(void)
{
  int i,j;
  double x;
  for(i=0;i<RULER_NUM;i++)
    {
      for(j=0;j<SEED_SIZE;j++)
	{
	  x=rand()%1000/1000.0;
	  if(x<PMUTATE) 
	    {
	      //need to add bias function???
	     seeds[i][j]=(rand()%MAX_NUM)+1;
	    }
	}
    } 
}





int selection(void)
{
  double fitTotal=0;
  double prob;
  int i,j,k;
  int **dummy;

  //make higher fitness better
  for(i=0;i<RULER_NUM;i++)
    seedFit[i]=(worstFit -seedFit[i])+1;

  for(i=0;i<RULER_NUM;i++)
    fitTotal+=seedFit[i];

  for(i=0;i<RULER_NUM;i++)
    {
      rSeedFit[i]=seedFit[i]/fitTotal;
    }
  cSeedFit[0]=rSeedFit[0];
  for(i=1; i<RULER_NUM;i++)
    {
      cSeedFit[i]=cSeedFit[i-1]+rSeedFit[i];
      
    }
  for(i=0;i<RULER_NUM;i++)
    {
      prob=rand()%1000/1000.0;
      if(prob<cSeedFit[0])
	{
	  for(j=0;j<SEED_SIZE;j++)
	    newSeeds[i][j]=seeds[0][j];
	}
      else
	{
	  for(j=0;j<RULER_NUM-1;j++)
	    {
	      if(prob>=cSeedFit[j] && prob<cSeedFit[j+1])
		{
		  for(k=0;k<SEED_SIZE;k++)
		    newSeeds[i][k]=seeds[j+1][k];
		  break; 
		}
	      
	    }
	}	  
    }
  //need to do pointer swapping
  for(i=0;i<RULER_NUM;i++)
    {
      for(j=0;j<SEED_SIZE;j++)
	seeds[i][j]=newSeeds[i][j];
    }
  return fitTotal;
}
double lookUpSeed(int seedIndex)
{
  int i,j,k;
  int sortMore=0;
  int temp;
  //might not sort the actual seed, effects the randomness??
  int seedCopy[SEED_SIZE];
  int tempSeed[SEED_SIZE];
  double tempFit;
  double returnValue;

  
  for(i=0;i<SEED_SIZE;i++)
    seedCopy[i]=seeds[seedIndex][i];


  //sort the seed
  for(i=0;i<SEED_SIZE-1;i++)
  {
    sortMore=0;
    for(j=i;j<SEED_SIZE-1;j++)
      {
	if(seedCopy[j]>seedCopy[j+1])
	  {
	    sortMore=1;
	    temp=seedCopy[j];
	    seedCopy[j]=seedCopy[j+1];
	    seedCopy[j+1]=temp;
	  }
      }
    if(sortMore==0)
      break;
  }
  
  //match to the lookup table
  temp=-1;
  for(i=0;i<LUT_SIZE;i++)
    {

      for(j=0;j<SEED_SIZE;j++)
	{ 
	  if(lookUpTable[i][j]==seedCopy[j])
	    temp =j;
	  else
	    break;
	}
      //when we have a match
      if(temp==SEED_SIZE-1)
	{
	  //swap with the oldest
	  for(j=0;j<SEED_SIZE;j++)
	    tempSeed[j]=lookUpTable[oldest][j];
	  tempFit=lookUpTableFit[oldest];
	  for(j=0;j<SEED_SIZE;j++)
	    lookUpTable[oldest][j]=lookUpTable[i][j];
	  lookUpTableFit[oldest]=lookUpTableFit[i];
	  returnValue=lookUpTableFit[i];
	  for(j=0;j<SEED_SIZE;j++)
	    lookUpTable[i][j]=tempSeed[j];
	  lookUpTableFit[i]=tempFit;
	  temp =-2;
	  break; 
	}
    
    }
  //when we have a new member write it over the oldest
  if(temp != -2)
    {
      returnValue=lookUpTableFit[oldest]=growRuler(seedIndex,0);
      for(i=0;i<SEED_SIZE;i++)
	lookUpTable[oldest][i]=seedCopy[i];
    }    
  oldest=(oldest+1)%LUT_SIZE;
  return returnValue;
}
int main()
{
  int i,j,k;
  int fitness,length;
  int bestFit=999999999;
  srand(time(NULL));
  for(k=0;k<LOOPS;k++)
    {
      //initialize one seed
      init(0); 
      //printf("\nRulers Initialized. Entering the GA.");
 
      printf("\n{");
      for(j=0;j<GENERATIONS;j++)
	{
	  //greedy construction of the seed
	  fitness=growRuler(0,0);
	  //output fitness/length
	  length=0;
	  printf("*");
	  for(i=0;i<SEED_SIZE;i++)
	    printf("%d ",seeds[0][i]);
	  printf("*");
	  for(i=0;i<RULER_SIZE;i++)
	    {
	        printf(" %d", rulerMap[i]);
	      if(rulerMap[i]>length)
		length=rulerMap[i];
	    }
	  if(length==fitness)
	    {
	      if(length<bestFit)
		bestFit=length;
	    }
	  printf("\n{%d,%d}",fitness,length);
	  if(j+1<GENERATIONS)
	    printf(",");
	  //generate a new seed from the construction
	  reInit(0);
	}
      printf("\n}%d\n",bestFit);
    }
  return 0;
}