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

  Project:  calculate optimum solutions for the hash definition problem

  Module:   HashDef.c (main program)

  Purpose:  Complete program

  Date:     1st November 2001

  Author:   D Cornforth

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

/* The author retains all rights to this program. You may copy and edit this program
   for any non-commercial purpose without permission from the author, as long as you
   acknowledge the source. Suggested acknowledgement:
   
   Cornforth, D. 2002 Evolution in the orange box - A new approach to the sphere-
   packing problem in CMAC-based neural networks. Proceedings of the 15th Australian
   Joint Conference on Artificial Intelligence.
*/

#pragma hdrstop
#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>

#pragma argsused

#define  MAX_DIMS   50
#define  MAX_LINE   100

enum Boolean {False, True};

// data with module scope
int Alleles = 0;
int *Coprime = NULL;
Boolean Best = False;
double *BestFitness = NULL;
int **BestVector = NULL;
int Dims = 0;
int EliteCount;
double *Fitness = NULL;
int Generations;
int **Genotype = NULL;
int MutationRange = 200;   // prob 1% of being -1, 1% of being +1
int Population = 0;
Boolean Prime = True;   // use coprime array
int Runs = 0;
Boolean Verbose = False;

// functions with module scope
void BestResult (int layers);
Boolean checkfile (ifstream* inp);
void CreateCoprimes (int layers);
void Crossover (int* parent1, int* parent2, int newpop);
void Evaluate (char* QuantFile);
double GetFitness (int layers, int* vector);
void InitialisePop (void);
void main (int argc, char *argv[]);
void MedianResult (int layers);
void Mutation (int layers, int start, int stop);
Boolean Reproduce (int layers);
int Roulette (double *cumulative);
void ShowBest (int layers, int run);
void SortIndex (double* item, int* index, int length);

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

  Function:   main

  Purpose:

  Arguments:  none

  Returns:    none

------------------------------------------------------------------------------*/
void main (int argc, char *argv[])
{
  int gen, pop, layers, min_layers, max_layers, run;

  time_t t;

  srand((unsigned) time(&t));

  if (argc < 2)
  {
    cerr << "Basic operation:\n";
    cerr << "     hashdef dims pop generations runs -> will produce table for 2 - 50 layers\n";
    cerr << "     hashdef dims pop generations runs layers -> will do one layer\n";
    cerr << "     hashdef filename -> will evaluate the grid spacings found in the file\n";
    cerr << "Additional flags (must come last):\n";
    cerr << "     -c -> don't restrict to coprimes.\n";
    cerr << "     -d or -v -> turn on debug (verbose) mode.\n";
    cerr << "     -b -> show best result from multiple runs (default: median result).\n";
    exit (0);
  }

  if (argc < 5)   // assume evaluation mode
  {
    if (argc == 3)
    {
      if ((strcmp (argv[2], "-d") == 0) || (strcmp (argv[2], "-v") == 0))
      {
        Verbose = True;   // show debug info
        cerr << "Verbose mode\n";
      }
    }
    Prime = False;
    Evaluate (argv[1]);
    exit (0);
  }

  Dims = atoi (argv[1]);
  if ((Dims < 2) || (Dims > MAX_DIMS))
  {
    cerr << "Invalid number of dims.\n";
    exit (0);
  }
  Population = atoi (argv[2]);
  if (Population < 2)
  {
    cerr << "Invalid population size.\n";
    exit (0);
  }
  Generations = atoi (argv[3]);
  if (Generations < 2)
  {
    cerr << "Invalid number of generations.\n";
    exit (0);
  }
  Runs = atoi (argv[4]);
  if (Runs < 1)
  {
    cerr << "Invalid number of runs.\n";
    exit (0);
  }

  min_layers = 2;
  max_layers = 50;
  int next_arg = 5;

  if (argc > 5)
  {
    if (*argv[next_arg] != '-')   // if this is not a switch
    {
      min_layers = atoi (argv[next_arg]);
      next_arg++;
      if (min_layers > 1)
      {
        max_layers = min_layers;
      }
      else
      {
        cerr << "Invalid number of layers.\n";
        exit (0);
      }
    }
  }

  Prime = True;   // default - use coprimes
  Verbose = False;
  Best = False;
  while (argc > next_arg)   // all remaining args are switches
  {
    if (strcmp (argv[next_arg], "-b") == 0)
    {
      Best = True;   // report best result, not median
    }
    else if (strcmp (argv[next_arg], "-c") == 0)
    {
      Prime = False;   // don't restrict to coprimes
    }
    else if (strcmp (argv[next_arg], "-d") == 0)
    {
      Verbose = True;   // show debug info
    }
    else if (strcmp (argv[next_arg], "-v") == 0)
    {
      Verbose = True;   // show debug info
    }
    next_arg++;
  }

  EliteCount = Population / 10;   // number of individuals passed to next generation
  if (EliteCount < 2)
  {
    EliteCount = 2;
  }
  if ((Population - EliteCount) % 2)   // if remainder is odd
  {
    EliteCount++;
  }

  Fitness = new double[Population];
  Genotype = new int*[Population];
  for (pop=0; pop<Population; pop++)
  {
    Genotype[pop] = new int[Dims];
  }

  for (layers=min_layers; layers<=max_layers; layers++)
  {
    cerr << "Layers: " << layers << endl;
    if (Prime)
    {
      Coprime = new int[layers];
      CreateCoprimes (layers);
    }
    BestVector = new int*[Runs];
    BestFitness = new double[Runs];
    for (run=0; run<Runs; run++)
    {
      InitialisePop ();
      for (gen=0; gen<Generations; gen++)
      {
        if (Reproduce (layers))   // if converged
        {
          cerr << "Population has converged.\n";
          break;
        }
      }
      BestVector[run] = new int[Dims];
      ShowBest (layers, run);
    }
    if (Best)
    {
      BestResult (layers);   // calc best vector from all the runs
    }
    else
    {
      MedianResult (layers);
    }
    for (run=0; run<Runs; run++)
    {
      delete BestVector[run];
    }
    delete[] BestVector;
    delete[] BestFitness;
    if (Prime)
    {
      delete[] Coprime;
    }
  }

  delete[] Fitness;
  for (pop=0; pop<Population; pop++)
  {
    delete Genotype[pop];
  }
  delete[] Genotype;
}

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

  Function:     BestResult

  Purpose:      Find the vector with the highest fitness from all runs

  Arguments:    ifstream

  Returns:

------------------------------------------------------------------------------*/
void BestResult (int layers)
{
  int run, best_ndx, dim;
  double most_fit = 0;

  for (run=0; run<Runs; run++)
  {
    if (BestFitness[run] > most_fit)
    {
      most_fit = BestFitness[run];
      best_ndx = run;
    }
  }

  cout << layers;
  for (dim=0; dim<Dims; dim++)
  {
    cout << " " << BestVector[best_ndx][dim];
  }
  cout << " " << most_fit << endl;
}

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

  Function:     checkfile

  Purpose:      check if a file was opened for reading

  Arguments:    ifstream

  Returns:

------------------------------------------------------------------------------*/
Boolean checkfile (ifstream* inp)
{
  Boolean success = True;
  filebuf* fb;

  fb = inp->rdbuf ();
  if (!fb->is_open ())
  {
    success = False;
  }

  return (success);
}

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

  Function:   CreateCoprimes

  Purpose:    create array containing suitable candidates for displacement vector
              these are coprimes of the number of layers, less than half no of layers

  Arguments:  none

  Returns:    none

------------------------------------------------------------------------------*/
void CreateCoprimes (int layers)
{
  int layer, ndx = 0;
  Boolean *noprime;

  Alleles = layers / 2 + 1;
  noprime = new Boolean[Alleles];

  // make a list of all numbers which are not coprimes
  for (layer=0; layer<Alleles; layer++)   // set the rest to True for now
  {
    noprime[layer] = True;   // assume this can be used
  }

  // eliminate non primes
  for (layer=2; layer<Alleles; layer++)   // for each layer number
  {
    if (noprime[layer] && ((double)layers / layer) == (double)(layers / layer))
    {   // if layers is a multiple of this number
      ndx = layer;
      while (ndx < Alleles)   // check all multiples too
      {
        noprime[ndx] = False;   // add to list of non primes
        ndx += layer;   // get next multiple of the non prime
      }
    }
  }

  Coprime[0] = 1;
  ndx = 1;
  for (layer=2; layer<Alleles; layer++)   // for each layer number
  {
    if (noprime[layer])
    {
      Coprime[ndx] = layer;
      ndx++;
    }
  }
  Alleles = ndx;

  if (Verbose)
  {
    cout << "Coprimes of " << layers << ":";
    for (ndx=0; ndx<Alleles; ndx++)   // for each layer number
    {
      cout << " " << Coprime[ndx];
    }
    cout << endl;
  }

  delete[] noprime;
}

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

  Function:   Crossover

  Purpose:

  Arguments:  oldpop - pointer to temporary copy of parents
              parent1, parent2 - indexes to the two parents in the array of old population
              newpop - index to position to put children

  Returns:    none

------------------------------------------------------------------------------*/
void Crossover (int* parent1, int* parent2, int newpop)
{
  int crosspoint, dim;

  crosspoint = rand () % Dims;
  for (dim=0; dim<crosspoint; dim++)
  {
    Genotype[newpop][dim] = parent1[dim];
    Genotype[newpop+1][dim] = parent2[dim];
  }
  for (dim=crosspoint; dim<Dims; dim++)
  {
    Genotype[newpop][dim] = parent2[dim];
    Genotype[newpop+1][dim] = parent1[dim];
  }
}

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

  Function:   Evaluate

  Purpose:

  Arguments:  none

  Returns:    none

------------------------------------------------------------------------------*/
void Evaluate (char* QuantFile)
{
  ifstream qfile;
  Boolean success = True;

  qfile.open (QuantFile);
  if (!checkfile (&qfile))
  {
    cerr << "Error opening file " << QuantFile << "\n";
    success = False;
  }

  if (success)   // if file was opened
  {
    cerr << "Opened " << QuantFile;
    char *c_point = strstr (QuantFile, "hash") + 4;
    //char* c_point = QuantFile + 4;
    Dims = strtol (c_point, NULL, 0);
    cerr << " Dimensions: " << Dims << endl;
    int* vector = new int[Dims];

    double fitness;
    char* in_line = new char[MAX_LINE];

    while (True)
    {
      qfile.getline (in_line, MAX_LINE);
      if (qfile.fail ())
      {
        break;
      }
      char* term = in_line;
      int layers = (int)strtol (in_line, &term, 0);
      cout << layers;
      // get quantisation offset from the line read in
      for (int dim=0; dim<Dims; dim++)
      {
        c_point = term;
        vector[dim] = (int)strtol (c_point, &term, 0);
        cout << " " << vector[dim];
      }
      fitness = GetFitness (layers, vector);
      cout << " " << fitness << endl;
    }

    delete[] vector;
  }

  qfile.close ();
}

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

  Function:   GetFitness

  Purpose:

  Arguments:  none

  Returns:    none

------------------------------------------------------------------------------*/
double GetFitness (int layers, int* genotype)
{
  int dim, ndx, temp, layer, coprimendx;
  int** vector;
  double dist, sum, min_dist, fitness;
  Boolean found = False;

  vector = new int*[layers];
  for (layer=0; layer<layers; layer++)
  {
    vector[layer] = new int[Dims];
  }

  vector[1][0] = 1;
  for (dim=1; dim<Dims; dim++)
  {
    if (Prime)
    {
      coprimendx = genotype[dim];   // the genotype contains indexes into the array of allowed values
      vector[1][dim] = Coprime[coprimendx];
    }
    else
    {
      vector[1][dim] = genotype[dim];
    }
  }

  // sort the vector in ascending order
  for (dim=1; dim<Dims-1; dim++)
  {
    for (ndx=dim+1; ndx<Dims; ndx++)
    {
      if (vector[1][dim] > vector[1][ndx])   // if not in ascending order
      {
        temp = vector[1][ndx];
        vector[1][ndx] = vector[1][dim];
        vector[1][dim] = temp;
      }
    }
  }

  // set first grid position is at all zeros
  for (dim=0; dim<Dims; dim++)
  {
    vector[0][dim] = 0;
  }

  // set the other layers by arithmetic progression modulo (Layers)
  for (layer=2; layer<layers; layer++)
  {
    for (dim=0; dim<Dims; dim++)
    {
      vector[layer][dim] = (vector[1][dim] * layer) % layers;
    }
  }

  // show the complete grid spacing vector
  if (Verbose)
  {
    for (dim=0; dim<Dims; dim++)
    {
      cout << " " << vector[1][dim];
    }
    cout << endl;
  }

  // evaluate the vector by calculating the minimum distance between successive centres checking adjacent cubes
  min_dist = layers;   // start with distance along one side, because this can be achieved
  for (layer=0; layer<layers-1; layer++)
  {
    for (ndx=layer+1; ndx<layers; ndx++)
    {
      // calculate distance from other vectors in the base hypercube
      sum = 0.0;
      found = True;
      for (dim=0; dim<Dims; dim++)
      {
        dist = fabs (vector[layer][dim] - vector[ndx][dim]);
        if (dist > min_dist)   // if dist in a single dimension is greater, Euclid cannot be less
        {
          dist = layers - dist;   // shift one of the points into an adjacent cube
          if (dist > min_dist)
          {
            found = False;
            break;   // give up, distance between these 2 points cannot be less than min_dist
          }
        }
        sum += (dist * dist);   // use squared distance to get Euclidean
      }
      if (found)
      {
        sum = sqrt (sum);
        if (sum < min_dist)
        {
          min_dist = sum;
          if (Verbose)
          {
            cout << "vectors";
            for (dim=0; dim<Dims; dim++)
            {
              cout << " " << vector[layer][dim];
            }
            cout << " compared with";
            for (dim=0; dim<Dims; dim++)
            {
              cout << " " << vector[ndx][dim];
            }
            cout << " gives distance = " << min_dist << endl;
          }
        }
      }
    }
  }

  fitness = min_dist;

  for (layer=0; layer<layers; layer++)
  {
    delete[] vector[layer];
  }
  delete[] vector;

  return (fitness);
}

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

  Function:   InitialisePop

  Purpose:    fills genotype with random numbers

  Arguments:  none

  Returns:    none

------------------------------------------------------------------------------*/
void InitialisePop (void)
{
  int pop, dim;

  for (pop=0; pop<Population; pop++)
  {
    for (dim=0; dim<Dims; dim++)
    {
      Genotype[pop][dim] = rand () % Alleles;
    }
  }
}

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

  Function:   MedianResult

  Purpose:    Find the vector with the median fitness from all runs of the GA

  Arguments:  none

  Returns:    none

------------------------------------------------------------------------------*/
void MedianResult (int layers)
{
  int r1, r2, dim, temp;
  int *order = new int[Runs];
  
  for (r1=0; r1<Runs; r1++)
  {
    order[r1] = r1;
  }

  // sort an index to the array of best fitness scores from each run
  for (r1=0; r1<Runs-1; r1++)
  {
    for (r2=r1+1; r2<Runs; r2++)
    {
      if (BestFitness[order[r1]] < BestFitness[order[r2]])
      {
        temp = order[r1];
        order[r1] = order[r2];
        order[r2] = temp;
      }
    }
  }

  if (Verbose)
  {
    for (r1=0; r1<Runs; r1++)
    {
      cerr << r1 << " " << BestFitness[order[r1]] << endl;
    }
  }

  cout << layers;
  for (dim=0; dim<Dims; dim++)
  {
    cout << " " << BestVector[order[Runs/2]][dim];
  }
  cout << " " << BestFitness[order[Runs/2]] << endl;
}

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

  Function:   Mutation

  Purpose:

  Arguments:  none

  Returns:    none

------------------------------------------------------------------------------*/
void Mutation (int layers, int start, int stop)
{
  int pop, dim, disturbance;

  for (pop=start; pop<stop; pop++)
  {
    for (dim=0; dim<Dims; dim++)
    {
      // valid values for disturbance are -1 and +1
      disturbance = (rand () % MutationRange) - 1;   // value from -1 to 198
      if (disturbance > 1)
      {
        disturbance = 0;   // mostly zero, sometimes can be 1 or -1
      }
      // add disturbance with wrap-around
      Genotype[pop][dim] = (Genotype[pop][dim] + disturbance) % Alleles;
      if (Genotype[pop][dim] < 0)
      {
        Genotype[pop][dim] += Alleles;   // do wrap around
      }
    }
  }
}

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

  Function:   Reproduce

  Purpose:    controls selection, new generation

  Arguments:  none

  Returns:    none

------------------------------------------------------------------------------*/
Boolean Reproduce (int layers)
{
  int pop, ndx1, ndx2, dim, elitendx;
  int* fitorder = new int[Population];
  double *cumulative = new double[Population];
  Boolean Converge;
  int **oldpop = new int*[Population];

  // make a copy of the whole population
  for (pop=0; pop<Population; pop++)
  {
    oldpop[pop] = new int[Dims];
    for (dim=0; dim<Dims; dim++)
    {
      oldpop[pop][dim] = Genotype[pop][dim];   // copy everything
    }
    fitorder[pop] = pop;
  }

  // build array of cumulative fitness for roulette selection
  Fitness[0] = GetFitness (layers, Genotype[0]);
  cumulative[0] = Fitness[0];
  for (pop=1; pop<Population; pop++)
  {
    Fitness[pop] = GetFitness (layers, Genotype[pop]);
    cumulative[pop] = cumulative[pop-1] + Fitness[pop];
  }
  SortIndex (Fitness, fitorder, Population);   // sorts in order of Fitness

  // Elitist strategy: copy the fittest EliteCount individuals to next generation
  for (pop=0; pop<EliteCount; pop++)
  {
    elitendx = fitorder[pop];   // get index of the next fittest
    for (dim=1; dim<Dims; dim++)
    {
      Genotype[pop][dim] = oldpop[elitendx][dim];   // copy fit individual
    }
  }

  // 2 offspring from each mating, so mate half as many times as remaining population
  for (pop=EliteCount; pop<Population; pop+=2)
  {
    // select two parents
    ndx1 = Roulette (cumulative);
    ndx2 = Roulette (cumulative);
    if ((ndx1 >= Population) || (ndx2 >= Population))
    {
      if (Verbose)
      {
        cerr << "Roulette has chosen invalid individuals\n";
      }
      break;
    }
    Crossover (oldpop[ndx1], oldpop[ndx2], pop);   // produce children
  }

  // destroy pop from previous generation
  for (pop=0; pop<Population; pop++)
  {
    delete oldpop[pop];
  }
  delete[] oldpop;

  Mutation (layers, EliteCount, Population);   //  mutate all children

  // test for convergence
  Converge = True;
  for (pop=1; pop<Population; pop++)
  {
    for (dim=0; dim<Dims; dim++)
    {
      if (Genotype[pop][dim] != Genotype[0][dim])   // only have to find one difference
      {
        Converge = False;
        break;   // no point in further tests if one difference found
      }
    }
    if (!Converge)
    {
      break;
    }
  }

  delete[] fitorder;
  delete[] cumulative;

  return (Converge);
}

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

  Function:   Roulette

  Purpose:    select parents using roulette wheel

  Arguments:  none

  Returns:    none

------------------------------------------------------------------------------*/
int Roulette (double *cumulative)
{
  int pop;
  double fitness;

  fitness = cumulative[Population-1] * rand () / RAND_MAX;

  for (pop=1; pop<Population; pop++)
  {
    if (cumulative[pop] > fitness)
    {
      break;
    }
  }

  return (pop);
}

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

  Function:   ShowBest

  Purpose:    show best individual for this run

  Arguments:  none

  Returns:    none

------------------------------------------------------------------------------*/
void ShowBest (int layers, int run)
{
  int pop, max_ndx = -1, dim, ndx, temp, coprimendx;

  // find the fittest individual
  BestFitness[run] = -1.0e+30;
  for (pop=0; pop<Population; pop++)
  {
    if (Fitness[pop] > BestFitness[run])
    {
      BestFitness[run] = Fitness[pop];
      max_ndx = pop;
    }
  }

  BestVector[run][0] = 1;
  for (dim=1; dim<Dims; dim++)
  {
    if (Prime)
    {
      coprimendx = Genotype[max_ndx][dim];   // the genotype contains indexes into the array of allowed values
      BestVector[run][dim] = Coprime[coprimendx];
    }
    else
    {
      BestVector[run][dim] = Genotype[max_ndx][dim];
    }
  }

  // sort the vectors in ascending order
  for (dim=1; dim<Dims-1; dim++)
  {
    for (ndx=dim+1; ndx<Dims; ndx++)
    {
      if (BestVector[run][dim] > BestVector[run][ndx])   // if not in ascending order
      {
        temp = BestVector[run][ndx];
        BestVector[run][ndx] = BestVector[run][dim];
        BestVector[run][dim] = temp;
      }
    }
  }

  if (Verbose)
  {   // display the vector
    cout << layers;
    for (dim=0; dim<Dims; dim++)
    {
      cout << " " << BestVector[run][dim];
    }
    cout << " " << BestFitness[run] << endl;
  }
}

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

  Function:   SortIndex

  Purpose:    sort an integer array in order given by a double array

  Arguments:  item - list of numbers which define the order
              index - integer array of indices
              length - length of both arrays

  Returns:    none

------------------------------------------------------------------------------*/
void SortIndex (double* item, int* index, int length)
{
  int itemp;

  for (int ndx1=0; ndx1<length-1; ndx1++)
  {
    for (int ndx2=ndx1+1; ndx2<length; ndx2++)
    {
      if (item[index[ndx1]] < item[index[ndx2]])
      {
        itemp = index[ndx1];
        index[ndx1] = index[ndx2];   // swap index only, not actual value
        index[ndx2] = itemp;
      }
    }
  }
}