#include <algorithm> //std::max(..) std::min(..)
#include <iostream> // cout
#include <vector>
#include <map>
#include <math.h>
#include <time.h>


using namespace std;
#define TrainSetSize 20
#define PI 3.141592653589793238463
#define epoch 1000000

#define POPULATIONSIZE 50 

#define DBG_ASSERT(f) { if(!(f)){ __debugbreak(); }; }
double sigmoidF(double x) { return (1.0f / (1.0f + std::exp(-x))); }
//double random(double min, double max) { return min + rand() / (RAND_MAX / (max - min + 1.0) + 1.0); }
double random(double min, double max) { return min + (max - min) * ((double)rand() / (double)RAND_MAX); }

double clamp(double n, double lower, double upper) { return std::max(lower, std::min(n, upper)); }

struct neuron {
    neuron() { id = rand();  bias = random(-1, 1); outputF = -1; }
    int                        id;

    std::map<int, neuron*>    incomingtargets;

    double                      bias;
    std::map<int, double >    incomingweights;
    double outputF; // f(x)  forward

    void connect(neuron* n, double weight = random(-1, 1)) {
        n->incomingtargets[this->id] = this;
        n->incomingweights[this->id] = weight;
    }
    double activate(const double* input) {
        if (input != NULL) {
            //this->outputB = 1;      // f'(x)
            this->outputF = *input; // f(x)
        }
        else {
            double sum = 0; // sum (x * w) + b
            map<int, neuron*>::iterator it;
            for (it = this->incomingtargets.begin(); it != this->incomingtargets.end(); it++)
            {
                sum += (this->incomingtargets[it->first]->outputF * this->incomingweights[it->first]);
            }
            sum += this->bias;
            this->outputF = sigmoidF(sum);// f(x)
        }
        return this->outputF;
    }
};

struct network {
    vector<neuron*> inputs;//  Input Layer w/ 1 neurons
    vector<neuron*> hiddens;// Hiddent Layer w/ 3 neurons
    vector<neuron*> outputs;// Output Layer w/ 1 neuron
                            // Connect Input Layer to Hidden Layer
    double error;

    network()
    {
        error = 0;
        inputs.push_back(new neuron()); // Input Layer w/ 1 neurons
        hiddens.push_back(new neuron());  hiddens.push_back(new neuron()); hiddens.push_back(new neuron()); // Hiddent Layer w/ 3 neurons
        outputs.push_back(new neuron());  // Output Layer w/ 1 neuron

        inputs[0]->connect(hiddens[0]);        inputs[0]->connect(hiddens[1]);        inputs[0]->connect(hiddens[2]);
        // Connect Hidden Layer to Output Layer
        hiddens[0]->connect(outputs[0]);    hiddens[1]->connect(outputs[0]);    hiddens[2]->connect(outputs[0]);
    }
    ~network()
    {
        delete inputs[0];
        delete hiddens[0];
        delete hiddens[1];
        delete hiddens[2];
        delete outputs[0];
    }

    double activate(const double* input, const double* output) {
        for (int i = 0; i < 1; i++) inputs[i]->activate(input);
        for (int i = 0; i < 3; i++) hiddens[i]->activate(NULL);
        for (int i = 0; i < 1; i++) outputs[i]->activate(NULL);
        if (output != NULL) error += abs(*output - outputs[0]->outputF);
        return outputs[0]->outputF;
    }

    double GetFitness() const { return error; }

    int getDimension() const
    {
        return 5 + 6;

    } // 5 neurons (5xbias), 6 connections (6xweights)

    double& getData(int dimensionIndex)
    {
        // bias and incomingweights
        neuron* allneurons[] = { inputs[0], hiddens[0], hiddens[1], hiddens[2], outputs[0] };
        int count = 0;
        for (int i = 0; i<5; i++)
        {
            if (count == dimensionIndex)
            {
                return allneurons[i]->bias;
            }
            count++;
            for (int k = 0; k<(int)allneurons[i]->incomingweights.size(); k++)
            {
                if (count == dimensionIndex)
                {
                    return allneurons[i]->incomingweights[k];
                }
                count++;
            }
        }
        DBG_ASSERT(false);
        return allneurons[0]->bias; // never get here!
    }

    void setData(int dimensionIndex, double val)
    {
        // bias and incomingweights
        neuron* allneurons[] = { inputs[0], hiddens[0], hiddens[1], hiddens[2], outputs[0] };
        int count = 0;
        for (int i = 0; i<5; i++)
        {
            if (count == dimensionIndex)
            {
                allneurons[i]->bias = val;
                return;
            }
            count++;
            for (int k = 0; k<(int)allneurons[i]->incomingweights.size(); k++)
            {
                if (count == dimensionIndex)
                {
                    allneurons[i]->incomingweights[k] = val;
                    return;
                }
                count++;
            }
        }
        DBG_ASSERT(false);
        //return allneurons[0]->bias; // never get here!
    }

};


void main() {
    srand((unsigned int)time(NULL));

    double F  = 0.1;  //differential weight [0,2]
    double CR = 0.5f; //crossover probability [0,1]

    for (int test = 0; test < 10; test++)
    {
        F  = random(0.001, 1.9); //differential weight [0,2]
        CR = random(0.001, 0.9); //crossover probability [0,1]

        vector<network*> population;// = new network*[POPULATIONSIZE];// = new network[ POPULATIONSIZE ];
        for (int i = 0; i < POPULATIONSIZE; i++)
        {
            population.push_back(new network());
        }


        vector<pair<double, double>> trainSet;    trainSet.resize(TrainSetSize);
        for (int i = 0; i < TrainSetSize; i++) { trainSet[i] = make_pair(i / (double)TrainSetSize, sin((i / (double)TrainSetSize) * (2.0 * PI))*0.5 + 0.5); }


        //dimensionality of problem, means how many variables problem has. this case 3 (data1,data2,data3)
        const int N = population[0]->getDimension(); // 5 neurons and each neuron has connected weights and bias (should be 11 for this case)


        for (int j = 0; j < POPULATIONSIZE; j++)
        {
            for (int ts = 0; ts < TrainSetSize; ts++)
            {
                const double input = trainSet[ts].first;
                const double output = trainSet[ts].second;
                population[j]->activate(&input, &output);
            }
        }

        for (int e = 0; e < epoch; e++)
        {

            //F = random(0.001, 0.5);
            //CR = random(0.01, 0.6);

            //main loop of evolution.
            for (int j = 0; j < POPULATIONSIZE; j++)
            {


                //calculate new candidate solution

                //pick random point from population
                int x = ::floor(random(0, POPULATIONSIZE - 1));

                int a, b, c;    //pick three different random points from population
                do { a = ::floor(random(0, POPULATIONSIZE - 1)); } while (a == x);
                do { b = ::floor(random(0, POPULATIONSIZE - 1)); } while (b == x || b == a);
                do { c = ::floor(random(0, POPULATIONSIZE - 1)); } while (c == x || c == a || c == b);

                // Pick a random index [0-Dimensionality]
                int R = ::floor(random(0, N));
                //Compute the agent's new position

                network* original = population[x];
                network* candidate = new network();// = original; 

                network* individual1 = population[a];
                network* individual2 = population[b];
                network* individual3 = population[c];

                //if(i==R | i<CR)
                //candidate=a+f*(b-c)
                //else
                //candidate=x

                for (int s = 0; s < N; ++s)
                {
                    //if (s == R || random(0, 1)<CR)
                    if (s == R || random(0, 1) < CR)
                    {
                        //candidate->getData(s) = individual1->getData(s) + F*(individual2->getData(s) - individual3->getData(s));
                        double val = individual1->getData(s) + F*(individual2->getData(s) - individual3->getData(s));
                        //val = val * 0.5;
                        candidate->setData(s, val);
                        //candidate->getData(s) *= 0.5;

                        //double val = candidate->getData(s);
                        //DBG_ASSERT(val > -1000 && val < 1000);
                    }
                    else
                    {
                        //candidate->getData(s) = original->getData(s);
                        double val = original->getData(s);
                        candidate->setData(s, val);
                    }
                }// End for b

                 // calculate fitness for our new candidate
                for (int ts = 0; ts < TrainSetSize; ts++)
                {
                    const double input = trainSet[ts].first;
                    const double output = trainSet[ts].second;
                    candidate->activate(&input, &output);
                }

                double candidateFitness = candidate->GetFitness();
                double originalFitness = original->GetFitness();

                // if better replace
                if (candidateFitness < originalFitness)
                {
                    population[x] = candidate;
                    delete original;
                }
                else
                {
                    delete candidate;
                }

            }// end population

#if 1 // print out the fitness as it evolves
            {
                int    bestFitnessIndex = 0;
                double bestFitnessValue = population[bestFitnessIndex]->GetFitness();

                for (int i = 1; i < POPULATIONSIZE; ++i)
                {
                    if (population[i]->GetFitness() < bestFitnessValue)
                    {
                        bestFitnessValue = population[i]->GetFitness();
                        bestFitnessIndex = i;
                    }
                }// End for i
                cout << "epoch:" << e << " fitness:" << bestFitnessValue << "\r";
            }
#endif

        }// epoch



        int    bestFitnessIndex = 0;
        double bestFitnessValue = population[bestFitnessIndex]->GetFitness();

        for (int i = 1; i < POPULATIONSIZE; ++i)
        {
            if (population[i]->GetFitness() < bestFitnessValue)
            {
                bestFitnessValue = population[i]->GetFitness();
                bestFitnessIndex = i;
            }
        }// End for i
        DBG_ASSERT(bestFitnessIndex > -1);



#if 0
        cout << "input, idea, predicted (y =sin(x))" << "\n";
        for (int i = 0; i < 100; i++) { //Print out test results (ideal vs predicted for sin wave)
            double input = i / 100.0;
            double predicted = population[bestFitnessIndex].activate(&input, NULL);
            std::cout << (input) << ", " << sin(input * (2.0 * PI))*0.5 + 0.5 << ", " << predicted << "\n";
    } system("pause"); // wait for user to press a key before closing
#endif

#if 1
    char filename[256] = { 0 };
    sprintf(filename, "data_F=%.5f_CR=%.5f.txt", F, CR);
    FILE* fp = fopen(filename, "wt");
    fprintf(fp, "epoch %d, populationsize: %d\n", epoch, POPULATIONSIZE);
    fprintf(fp, "error %f\n", bestFitnessValue);
    fprintf(fp, "input, idea, predicted (y =sin(x))\n");
    for (int i = 0; i < 100; i++) { //Print out test results (ideal vs predicted for sin wave)
        double input = i / 100.0;
        double predicted = population[bestFitnessIndex]->activate(&input, NULL);
        double exact = sin(input * (2.0 * PI))*0.5 + 0.5;
        fprintf(fp, "%f %f %f \n", input, exact, predicted);
    }
    fflush(fp);
    fclose(fp);
    fp = NULL;

#endif

    for (int i = 0; i < POPULATIONSIZE; i++)
    {
        delete population[i];
    }
    population.resize(0);

    }// end for test

    system("pause"); // wait for user to press a key before closing 
}