using System;
using System.Collections.Generic;
using RandomNumbers;

namespace GeneticAlgorithm
{
    public class GA
    {
        private int nsteps, population, nbits;
        private double beta, energy_LB, temp; //min_energy is LB on objective
        private double[] J, fitness, probability, min_energy;
        private double[,] spins, newspins, energy; //the cuurent and future design vector(s). The energy vector is time history of objective
        RandomGenerator randU, randBit, randG;
        #region Constructor & Initialization
        public GA(int nsteps, int population, double beta)
        {
            this.nsteps = nsteps;
            this.population = population;
            this.beta = beta;
            initialize();
        }
        private void initialize()
        {
            // Required data fields
            nbits = 100;
            J = new double[nbits];
            probability = new double[population];
            fitness = new double[population];
            spins = new double[population, nbits];
            newspins = new double[population, nbits];
            energy = new double[population, nsteps];
            min_energy = new double[nsteps];
            // Random Number Generators
            randG = new RandomGaussian(1, 0, 1); //(seed, mean, variance)
            randU = new RandomUniform(1, 0, 1); //(seed, upper_bound, lower_bound)
            randBit = new RandomUniform(1, 0, nbits-1);
            // Set the J's and find the lower bound on the energy
            energy_LB = 0;
            for (int i = 0; i < nbits; i++)
            {
                J[i] = randG.Next();
                energy_LB += -Math.Abs(J[i]);                
            }
            //Create the initial population
            for (int i = 0; i < population; i++)
            {
                for (int j = 0; j < nbits; j++)
                {
                    if (randU.Next() > 0.5) { spins[i, j] = 1; } 
                    else { spins[i, j] = -1; }
                }
            }
        }
        public double LB { get { return energy_LB; } }
        #endregion Constructor & Initialization

        public double[] Run()
        {
            double sum = 0;
            double sum_prob = 0;
            int index1, index2, crossover;
            for (int t = 0; t < nsteps; t++)
            {
                //1. Evaluate the energy (objective function) for each member of the population, and assign a probability of reproduction
                // This is known as 'Roulette Wheel Selection'
                for (int i = 0; i < population; i++)
                {
                    for (int j = 0; j < nbits - 1; j++)
                    {
                        energy[i, t] += spins[i, j] * spins[i, j + 1] * J[j];
                    }
                    energy[i, t] += spins[i, nbits - 1] * spins[i, 0] * J[nbits - 1];//periodic BC's
                    fitness[i] = Math.Exp(-beta * (energy[i, t] - energy_LB)); 
                    sum += fitness[i]; //sum all fitnesses
                }
                for (int i = 0; i < population; i++)
                {
                    //scale by sum of fitnesses so that probabilities of selection are 0 <= p <= 1. 
                    probability[i] = sum_prob + fitness[i]/sum;
                }
                for (int i = 1; i < population; i++)
                {
                    //calculate the cumulative sum of the probabilities
                    probability[i] += probability[i - 1];
                }
                probability[population - 1] = 1; //correct rounding errors
                //2. Reproduce from the old population
                for (int i = 0; i < population; i++)
                {
                    //(a) Draw from distribution
                    index1 = 0; index2 = 0;
                    double r1 = randU.Next();
                    double r2 = randU.Next();
                    for (int k = 0; k < population; k++)
                    {
                        if (r1 > probability[k]) { index1 += 1; } // population members which contribute little to the cdf are likely to be passed over                         
                    }
                    for (int k = 0; k < population; k++)
                    {
                        if (r2 > probability[k]) { index2 += 1; }
                    }
                    //(b) Crossover
                    crossover = (int)randBit.Next(); //choose a crossover location
                    for (int j = 0; j < crossover; j++)
                    {
                        newspins[i, j] = spins[index1, j];
                    }
                    for (int j = crossover; j < nbits; j++)
                    {
                        newspins[i, j] = spins[index2, j];
                    }
                    //(c) Mutate a random bit
                    newspins[i, (int)randBit.Next()] *= -1;
                }
                //3. Find the best solution in the population
                temp = 1e10;
                for (int i = 0; i < population; i++)
                {
                    if (energy[i, t] < temp) { temp = energy[i, t]; }
                }
                min_energy[t] = temp;
                spins = newspins;
            }
            return min_energy;
        }
    }
}