//
//                neural_ga.cpp
//This code uses a genetic algorithm to train a
//single-hidden-layer feedfoward artificial neural network.
//
//The user specifices a data file containing training vectors of the form
//(input1, input2, ... , inputn, output1, output2, ... , outputn).  The file
//must be ASCII text format, with spaces separating the components of each
//vector, and hard returns separating the vectors.
//The user defines the number of nodes in each layer of the network.
//
//The genome for each individual is a series of floating point numbers
//representing weights on the output nodes and then on the hidden nodes.  The
//genetic algorithm (after Montana and Davis) uses fitness-proportionate
//selection, multipoint crossover which keeps weights on each node together,
//and mutatation.  The mutation operator adds
//a small random number to each weight selected for mutation.
//
//The training vectors are split (80/20) into a training set and a test set.
//The training set data are used to evolve the population; at the end of
//training (based only on completing the specified number of generations)
//the most fit individual is selected using the test set, and the weight values
//this individual are written into a file entitled "trainw.dat".
//
//This code is written as a Windows console application.  Removing the conio.h
//include and the getch() statements will make it generic.
//
//
//Copyright 1998 Michael J. Wax
//You may use this code as is, or incorporate it in another program,
//as long as proper attribution is given.  However, I make no warranties,
//express or implied, regarding the performance of this code, its freedom from
//error, or its suitability for use in a given application.
//
//

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

//use 80% of training vectors for training, 20% for testing

class node {
  public:
  double *weight, output;
  //constructor
  node () {}
  double calcout (double &activation){
	 return 1/(1+exp(-activation));
	 }
  };

class hiddennode : public node {
  };

class outputnode : public node {
  };

class datavector {
  public:
  double *vinput, *voutput;
  };

class chromosome {
  public:
  double fitness, *gene;
  //constructor initializes gene
  chromosome () {}
  //mutation function - adds random number between -1 and 1 to gene
  //probability of mutation on each gene is 0.01
  void mutate (int gene_number) {
    if (double (rand()) / double (RAND_MAX) < 0.01) {
      gene[gene_number] += 2 * (double (rand()) / double (RAND_MAX) - 0.5);
      }
    }
  };

void swap (chromosome list[], int a, int b) {
  chromosome v;
  v = list[a];
  list[a] = list[b];
  list[b] = v;
  }

void quicksort (chromosome list[], int left, int right) {
  int i, j;
  chromosome v;
  if (right > left) {
    v = list[right];
    i = left - 1;
    j = right;
    for (;;) {
      while (list[++i].fitness < v.fitness);
      while (list[--j].fitness > v.fitness);
      if (i >= j) break;
      swap (list, i, j);
      }
    swap (list, i, right);
    quicksort (list, left, i-1);
    quicksort (list, i+1, right);
    }
  }

void main() {
  double scratch, generation_fitness, activation, *wheel, pc = 0.7, sumsquarederror;
  //crossover probability is 0.7 on each node
  int counter, counter2, counter3, master, numin, numhidden, numout, num_genes, population,
   max_generations, generation = 0, numvectors, numtest, individual, individual2;
  char file_name[15];
  datavector *t_vector;
  hiddennode *hn;
  outputnode *outn;
  srand( (unsigned)time( NULL ) );

  cout << "File name containing data?  ";
  cin >> file_name;
  cout << "How many data vectors? ";
  cin >> numvectors;
  cout << "Network structure:" << endl;
  cout << " How many input nodes? ";
  cin >> numin;
  cout << " How many hidden nodes? ";
  cin >> numhidden;
  cout << " How many output nodes? ";
  cin >> numout;
  cout << "How many individuals? ";
  cin >> population;
  cout << "How many generations? ";
  cin >> max_generations;

  //set up network
  cout << "Setting up network . . ." << endl;
  hn = new hiddennode[numhidden];
  outn = new outputnode[numout];
	 for (counter = 0; counter < numhidden; ++counter) {
	   hn[counter].weight = new double[numin];
	   }
	 for (counter = 0; counter < numout; ++counter) {
	   outn[counter].weight = new double[numhidden];
		}

  //create training vectors and read in from file
  cout << "Reading in data vectors . . ." << endl;
  t_vector = new datavector[numvectors];
  ifstream training_input (file_name, ios::in);
	 for (counter = 0; counter < numvectors; ++counter) {
	 t_vector[counter].vinput = new double[numin];
	 t_vector[counter].voutput = new double[numout];
		for (counter2 = 0; counter2 < numin; ++counter2) {
		training_input >> t_vector[counter].vinput[counter2];
		}
		for (counter2 = 0; counter2 < numout; ++counter2) {
		training_input >> t_vector[counter].voutput[counter2];
		}
	 }
  training_input.close ();

  //calculate division between training and test vectors
  numtest = int (0.8*numvectors);

  cout << "Training vectors: " << numtest << "; test vectors: " << (numvectors-numtest) <<endl;
  getch();

  //create population of genotypes
  chromosome *genotype;
  chromosome *temp;
  num_genes = numin * numhidden + numhidden * numout; //= number of weights
  genotype = new chromosome[population];
  temp = new chromosome[population];      //store offspring here while creating next generation
  wheel = new double[population];

  for (individual = 0; individual < population; ++individual) {
    genotype[individual].gene = new double[num_genes];
    temp[individual].gene = new double[num_genes];
    for (counter = 0; counter < num_genes; ++counter) {
      genotype[individual].gene[counter] = double (rand()) / double (RAND_MAX) - 0.5;
      }
    }

  //
  //START TRAINING HERE
  //
  cout << "Starting first generation . . ." << endl;

  while (generation < max_generations) {
  ++generation;

    for (individual = 0, generation_fitness = 0, sumsquarederror = 0; individual < population; ++individual) {
      //for each genotype, create a network
      for (counter = 0; counter < numout; ++counter) {
        for (counter2 = 0; counter2 < numhidden; ++counter2) {
          outn[counter].weight[counter2] = genotype[individual].gene[counter*numhidden+counter2];
          }
        }
      for (counter = 0; counter < numhidden; ++counter) {
        for (counter2 = 0; counter2 < numin; ++counter2) {
          hn[counter].weight[counter2] = genotype[individual].gene[counter*numin+counter2+numout*numhidden];
          }
        }

      //evaluate error of network with respect to each training vector
      for (master = 0, scratch = 0; master < numtest; ++master) {

        //calculate outputs of hidden layer nodes
		  for (counter = 0; counter < numhidden; ++counter) {
          activation = 0;
          for (counter2 = 0; counter2 < numin; ++counter2) {
          activation += hn[counter].weight[counter2] * t_vector[master].vinput[counter2];
            }
          hn[counter].output = hn->calcout(activation);
		    }

	     //calculate outputs of output layer nodes
		  for (counter = 0; counter < numout; ++counter) {
          activation = 0;
          for (counter2 = 0; counter2 < numhidden; ++counter2) {
            activation += outn[counter].weight[counter2] * hn[counter2].output;
            }
		    outn[counter].output = outn->calcout(activation);
          }

  	     //calculate errors of output layer nodes
		  for (counter = 0; counter < numout; ++counter) {
          //add squared error of output node to total error of individual
		    scratch += (t_vector[master].voutput[counter]-outn[counter].output)*(t_vector[master].voutput[counter]-outn[counter].output);
		    }

	   } //done looking at all training vectors
    genotype[individual].fitness = 1/(scratch + 1);
    generation_fitness += genotype[individual].fitness;

      //now calculate sum-squared error over test set
      for (master = numtest; master < numvectors; ++master) {

        //calculate outputs of hidden layer nodes
		  for (counter = 0; counter < numhidden; ++counter) {
          activation = 0;
          for (counter2 = 0; counter2 < numin; ++counter2) {
          activation += hn[counter].weight[counter2] * t_vector[master].vinput[counter2];
            }
          hn[counter].output = hn->calcout(activation);
		    }

	     //calculate outputs of output layer nodes
		  for (counter = 0; counter < numout; ++counter) {
          activation = 0;
          for (counter2 = 0; counter2 < numhidden; ++counter2) {
            activation += outn[counter].weight[counter2] * hn[counter2].output;
            }
		    outn[counter].output = outn->calcout(activation);
          }

  	     //calculate errors of output layer nodes
		  for (counter = 0; counter < numout; ++counter) {
          //add squared error of output node to total error of individual
		    sumsquarederror += (t_vector[master].voutput[counter]-outn[counter].output)*(t_vector[master].voutput[counter]-outn[counter].output);
		    }

        }

    } //done all individuals in population

  //now create selection wheel
  for (individual = 0, scratch = 0; individual < population; ++individual) {
    wheel[individual] = genotype[individual].fitness / generation_fitness + scratch;
    scratch = wheel[individual];
    }

  //create new population
  for (counter = 0; counter < population; counter +=2) {
    //pick two parents randomly based on fitness
	 scratch = double (rand()) / double (RAND_MAX);
	 individual = 0;
	 while (wheel[individual] < scratch && individual < population) {++individual;}
	 //genotype[individual] is the first selected individual
	 individual2 = individual;
	 //don't select same individual twice
	 while (individual == individual2) {
		scratch = double(rand()) / double (RAND_MAX);
		individual2 = 0;
		while (wheel[individual2] < scratch && individual2 < population) {++individual2;}
		  }
    //genotype[individual2] is the second selected individual

    //multipoint crossover, keeping weights to each node together
    for (counter2 = 0; counter2 < numout; ++counter2) {
      if (double (rand())/double (RAND_MAX) < pc) {
        for (counter3 = counter2 * numhidden; counter3 < (counter2 + 1) * numhidden; ++counter3) {
          temp[counter].gene[counter3] = genotype[individual].gene[counter3];
          temp[counter+1].gene[counter3]= genotype[individual2].gene[counter3];
          }
        }
      else {
        for (counter3 = counter2 * numhidden; counter3 < (counter2+1) * numhidden; ++counter3) {
          temp[counter].gene[counter3] = genotype[individual2].gene[counter3];
          temp[counter+1].gene[counter3]= genotype[individual].gene[counter3];
          }
        }
      }
    for (counter2 = 0; counter2 < numhidden; ++counter2) {
      if (double (rand())/double (RAND_MAX) < pc) {
        for (counter3 = numout * numhidden + counter2 * numin; counter3 < numout*numhidden + (counter2+1) * numin; ++counter3) {
          temp[counter].gene[counter3] = genotype[individual].gene[counter3];
          temp[counter+1].gene[counter3]= genotype[individual2].gene[counter3];
          }
        }
      else {
        for (counter3 = numout * numhidden + counter2 * numin; counter3 < numout*numhidden + (counter2+1) * numin; ++counter3) {
          temp[counter].gene[counter3] = genotype[individual2].gene[counter3];
          temp[counter+1].gene[counter3]= genotype[individual].gene[counter3];
          }
        }
      }
    }

  //now, write new population onto old population
  for (individual = 0; individual < population; ++individual) {
	 genotype[individual] = temp[individual];
	 }

  //mutate in place
  for (individual = 0; individual < population; ++individual) {
    for (counter2 = 0; counter2 < num_genes; ++counter2) {
      genotype[individual].mutate(counter2);
      }
    }

  cout << "Generation " << generation << " fitness = " << generation_fitness
       << "/" << population
       << "; test set sum-squared error = " << sumsquarederror << endl;
  } //end while, i.e., done all generations

  //print actual and predicted output values for best individual

  for (individual = 0, generation_fitness = 0; individual < population; ++individual) {
    //for each genotype, create a network
    for (counter = 0; counter < numout; ++counter) {
      for (counter2 = 0; counter2 < numhidden; ++counter2) {
        outn[counter].weight[counter2] = genotype[individual].gene[counter*numhidden+counter2];
        }
      }
    for (counter = 0; counter < numhidden; ++counter) {
      for (counter2 = 0; counter2 < numin; ++counter2) {
        hn[counter].weight[counter2] = genotype[individual].gene[counter*numin+counter2+numout*numhidden];
        }
      }

    //evaluate error of network with respect to each training vector
    for (master = 0, scratch = 0; master < numtest; ++master) {

      //calculate outputs of hidden layer nodes
      for (counter = 0; counter < numhidden; ++counter) {
        activation = 0;
        for (counter2 = 0; counter2 < numin; ++counter2) {
          activation += hn[counter].weight[counter2] * t_vector[master].vinput[counter2];
          }
        hn[counter].output = hn->calcout(activation);
	     }

      //calculate outputs of output layer nodes
      for (counter = 0; counter < numout; ++counter) {
        activation = 0;
        for (counter2 = 0; counter2 < numhidden; ++counter2) {
          activation += outn[counter].weight[counter2] * hn[counter2].output;
          }
	     outn[counter].output = outn->calcout(activation);
        }

      //calculate errors of output layer nodes
      for (counter = 0; counter < numout; ++counter) {
        //add squared error of output node to total error of individual
  	     scratch += (t_vector[master].voutput[counter]-outn[counter].output)*(t_vector[master].voutput[counter]-outn[counter].output);
	     }

	   } //done looking at all training vectors

    genotype[individual].fitness = 1/(scratch + 1);
    generation_fitness += genotype[individual].fitness;
    }

    //sort genotypes by fitness
    quicksort(genotype, 0, population-1);

    for (individual = 0; individual < population; ++individual) {
      cout << "Fitness of individual " << individual << " = " << genotype[individual].fitness << endl;
      }

    cout << "Population fitness = " << generation_fitness << endl;
    getch();

    //print out and write to file weights of highest fitness individual
    cout << "Weights of most fit individual:" << endl;
    ofstream weights_output ("trainw.dat", ios::out);
    for (counter = 0; counter < num_genes; ++counter) {
      cout << genotype[population-1].gene[counter] << endl;
      weights_output << genotype[population-1].gene[counter] << endl;
      }
    weights_output.close();

    getch();

    //print out actual and predicted values for most fit individual
    for (counter = 0; counter < numout; ++counter) {
      for (counter2 = 0; counter2 < numhidden; ++counter2) {
        outn[counter].weight[counter2] = genotype[population-1].gene[counter*numhidden+counter2];
        }
      }
    for (counter = 0; counter < numhidden; ++counter) {
      for (counter2 = 0; counter2 < numin; ++counter2) {
        hn[counter].weight[counter2] = genotype[population-1].gene[counter*numin+counter2+numout*numhidden];
        }
      }

    //evaluate network output for best individual on test set
    for (master = numtest, sumsquarederror = 0; master < numvectors; ++master) {

      //calculate outputs of hidden layer nodes
      for (counter = 0, activation = 0; counter < numhidden; ++counter) {
        for (counter2 = 0; counter2 < numin; ++counter2) {
          activation += hn[counter].weight[counter2] * t_vector[master].vinput[counter2];
          }
        hn[counter].output = hn->calcout(activation);
	     }

      //calculate outputs of output layer nodes
      cout << "vector " << master << endl;
      for (counter = 0, activation = 0; counter < numout; ++counter) {
        for (counter2 = 0; counter2 < numhidden; ++counter2) {
          activation += outn[counter].weight[counter2] * hn[counter2].output;
          }
	     outn[counter].output = outn->calcout(activation);
        sumsquarederror += (t_vector[master].voutput[counter]-outn[counter].output) * (t_vector[master].voutput[counter]-outn[counter].output);
        cout << "(" << t_vector[master].voutput[counter] << "," << outn[counter].output << ")  ";
        }
      cout << endl;
      }
    cout << "Sum squared error of best individual on test set = " << sumsquarederror << endl;
    getch();

}
//