#include <cstdlib>
#include <algorithm>
#include <vector>
#include <cmath>
#include "jdetevc2006.h"
#include "rnd_uni.h"

// Constructor, destructor, random
jDEtevc2006::jDEtevc2006(EvaluateFun evaluate, int D, int NP, int maxNFEVALs,
                         double *xmin, double *xmax, unsigned int RNi)
                             : evaluate(evaluate), D(D), NP(NP), maxNFEVALs(maxNFEVALs),
                             RNi(RNi),
                             idxFE(D+2), idxF(D), idxCR(D+1)  {//C after FE!
    best = trial = pop = newPop = 0; // reserved later in createRuntimeStorage!
    this->xmin = new double[D];
    this->xmax = new double[D];
    for (int i = 0; i < D; i++) {
        this->xmin[i] = xmin[i];
        this->xmax[i] = xmax[i];
    }
    NFevals = 0;
    rnd_uni_init = -(long)RNi;
    srand(RNi);

    L = D + 3; // FE, F, CR
}
jDEtevc2006::~jDEtevc2006() {
    delete [] xmin;
    delete [] xmax;
    delete [] best;
    delete [] trial;
    delete [] pop;
    delete [] newPop;
}
void jDEtevc2006::createRuntimeStorage() {
    best   = new double[L];
    trial  = new double[L];
    pop    = new double[NP * L];
    newPop = new double[NP * L];

    //nonzero init!
    const int initVal = 555; // a beast :p
    for (int i = 0; i < L; i++) {
        best[i] = initVal;
        trial[i] = initVal;
    }
    for (int i = 0; i < NP; i++)
        for (int j = 0; j < L; j++) {
            pop[i*L + j] = initVal;
            newPop[i*L + j] = initVal;
        }
}

// Utility functions
bool jDEtevc2006::isBetterFirst(const double *first, const double *second) const {
    const double &f1 = first[idxFE];
    const double &f2 = second[idxFE];

    return f1 < f2;
}
void jDEtevc2006::replaceFirstVectorIntoSecond(const double *first, double *second) const {
    for (int j = 0; j < L; j++)
        second[j] = first[j];
}
double jDEtevc2006::random() {
//    return (double)rand() / RAND_MAX;
    return rnd_uni(&rnd_uni_init);
    RNi = (69069 * RNi) % (((unsigned long long )2)<<32);
    return (double)RNi / (double)(((unsigned long long )2)<<32);
}

// Initialization & optimization
void jDEtevc2006::initPop() {
    createRuntimeStorage();
    for (int i = 0; i < NP && NFevals < maxNFEVALs; i++) {
        pop[i*L + idxF] = 0.1 + 0.9 * random(); // F
        pop[i*L + idxCR] = random(); // CR

        for (int j = 0; j < D; j++)
            pop[i*L + j] = xmin[j] + random() * (xmax[j] - xmin[j]); // uniform initialization

        evaluate(pop + i*L, pop + i*L + idxFE);
        NFevals++;
        if (i==0 || isBetterFirst(pop+i*L, best))
            replaceFirstVectorIntoSecond(pop + i*L, best);
    }
}

void jDEtevc2006::optimize(double *result, double *evaluation) {
    initPop();

    for (int g = 0; NFevals < maxNFEVALs; g++) // main generation loop
        optimizeOneGeneration(NFevals, g);

    for (int i = 0; i < D; i++) result[i] = best[i]; // return search result
    *evaluation = best[idxFE]; // return evaluation
}

void jDEtevc2006::optimizeOneGeneration(int &NFevals, int g) {
    onBeforeOptimizeOneGeneration(NFevals, g);
    for (int i = 0; i < NP && NFevals < maxNFEVALs; i++) { // iteration loop
        onBeforeOptimizeOneVector(NFevals, g, i);
 
        selfadaptControlParameters(i); //first, new control parameters

        computeMutationVector(i); // perturbs it to generate mutant vector

        computeCrossoverTrialVector(i);

        selectionTrialVsTarget(i);

        onAfterOptimizeOneVector(NFevals, g, i);
    }
    onAfterOptimizeOneGeneration(NFevals, g);
    double *tmp = newPop; newPop = pop; pop = tmp; // switch populations
}

void jDEtevc2006::onBeforeOptimizeOneGeneration(int &NFevals, int g) {
}

void jDEtevc2006::onAfterOptimizeOneGeneration(int &NFevals, int g) {
}

void jDEtevc2006::onBeforeOptimizeOneVector(int NFevals, int g, int i) {
}

void jDEtevc2006::onAfterOptimizeOneVector(int NFevals, int g, int i) {
}

void jDEtevc2006::selfadaptControlParameters(int targetIdx) {
    if (random() < 0.1) // F
        trial[idxF] = 0.1 + 0.9 * random();
    else trial[idxF] = pop[targetIdx*L + idxF];

    if (random() < 0.1) // CR
        trial[idxCR] = random();
    else trial[idxCR] = pop[targetIdx*L + idxCR];
}

void jDEtevc2006::computeMutationVector(int i) {
    int r1, r2, r3;

    do { r1 = (int)(random() * NP); } while (r1 == i);// select differing indices r1, r2, r3
    do { r2 = (int)(random() * NP); } while (r2 == i || r2 == r1);
    do { r3 = (int)(random() * NP); } while (r3 == i || r3 == r1 || r3 == r2);

    for (int j = 0; j < D; j++) // mutation
        trial[j] = pop[r1*L + j] + trial[idxF] * (pop[r2*L + j] - pop[r3*L + j]);
}

void jDEtevc2006::computeCrossoverTrialVector(int targetIdx) {
    for (int j = 0; j < D; j++) // crossover
        if (random() > trial[idxCR])
            trial[j] = pop[targetIdx*L + j];
}

bool jDEtevc2006::selectionTrialVsTarget(int i) {
    for (int j = 0; j < D; j++) {
        if (trial[j] < xmin[j]) trial[j] = xmin[j];
        if (trial[j] > xmax[j]) trial[j] = xmax[j];
    }

    evaluate(trial, trial+idxFE); // evaluation
    NFevals++;

    if (isBetterFirst(trial, pop+i*L)) { // selection
        replaceFirstVectorIntoSecond(trial, newPop + i*L);
        checkUpdateBest(best);
        return true;
    } else replaceFirstVectorIntoSecond(pop + i*L, newPop + i*L);

    return false;
}

void jDEtevc2006::checkUpdateBest(double *best) {
   if (isBetterFirst(trial, best))
        replaceFirstVectorIntoSecond(trial, best);
}
