evoga.cpp

 /********************************************************************************
 *  FARSA Experiments Library                                                   *
 *  Copyright (C) 2007-2012                                                     *
 *  Stefano Nolfi <stefano.nolfi@istc.cnr.it>                                   *
 *  Onofrio Gigliotta <onofrio.gigliotta@istc.cnr.it>                           *
 *  Gianluca Massera <emmegian@yahoo.it>                                        *
 *  Tomassino Ferrauto <tomassino.ferrauto@istc.cnr.it>                         *
 *                                                                              *
 *  This program is free software; you can redistribute it and/or modify        *
 *  it under the terms of the GNU General Public License as published by        *
 *  the Free Software Foundation; either version 2 of the License, or           *
 *  (at your option) any later version.                                         *
 *                                                                              *
 *  This program is distributed in the hope that it will be useful,             *
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of              *
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               *
 *  GNU General Public License for more details.                                *
 *                                                                              *
 *  You should have received a copy of the GNU General Public License           *
 *  along with this program; if not, write to the Free Software                 *
 *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA  *
 ********************************************************************************/

#include "evoga.h"
#include "evodataviewer.h"
#include "logger.h"
#include "randomgenerator.h"
#include <configurationhelper.h>
#include <QVector>
#include <QThreadPool>
#include <QtConcurrentMap>
#include <QtAlgorithms>
#include <QTime>
#include <QFile>

#include <cmath>

#include <Eigen/Core>
#include <Eigen/Dense>
#include <unsupported/Eigen/MatrixFunctions>

#include <exception>

// All the suff below is to avoid warnings on Windows about the use of unsafe
// functions. This should be only a temporary workaround, the solution is stop
// using C string and file functions...
#if defined(_MSC_VER)
    #pragma warning(push)
    #pragma warning(disable:4996)
#endif

namespace farsa {

class EvaluatorThreadForEvoga
{
public:
    EvaluatorThreadForEvoga(Evoga *ga, EvoRobotExperiment *exp) :
        m_ga(ga),
        m_exp(exp)
    {
    }

    ~EvaluatorThreadForEvoga()
    {
        delete m_exp;
    }

    void setGenotype(int id)
    {
        m_id = id;
        m_exp->setNetParameters(m_ga->getGenes(id));
    }

    int getGenotypeId() const
    {
        return m_id;
    }

    void run()
    {
        m_exp->doAllTrialsForIndividual(m_id);
        m_fitness = m_exp->getFitness();
    }

    double getFitness() const
    {
        return m_fitness;
    }

    EvoRobotExperiment* getExperiment()
    {
        return m_exp;
    }

private:
    Evoga *const m_ga;

    EvoRobotExperiment *const m_exp;

    int m_id;

    double m_fitness;
};

void runEvaluatorThreadForEvoga(EvaluatorThreadForEvoga* e)
{
    e->run();
}

class EvaluatorThreadForXnes
{
public:
    EvaluatorThreadForXnes(Evoga *ga, EvoRobotExperiment *exp) :
        m_ga(ga),
        m_exp(exp)
    {
    }

    ~EvaluatorThreadForXnes()
    {
        delete m_exp;
    }

    void setGenotype(int id, float* genes)
    {
        m_id = id;
        m_exp->setNetParameters(genes);
    }

    int getGenotypeId() const
    {
        return m_id;
    }

    void run()
    {
        m_exp->doAllTrialsForIndividual(m_id);
        m_fitness = m_exp->getFitness();
        m_error = m_exp->getError();
    }

    double getFitness() const
    {
        return m_fitness;
    }

    double getError() const
    {
        return m_error;
    }

    EvoRobotExperiment* getExperiment()
    {
        return m_exp;
    }

private:
    Evoga *const m_ga;

    EvoRobotExperiment *const m_exp;

    int m_id;

    double m_fitness;

    double m_error;
};

void runEvaluatorThreadForXnes(EvaluatorThreadForXnes* e)
{
    e->run();
}

struct FitnessAndId
{
    float fitness;

    int id;
};

bool operator<(FitnessAndId first, FitnessAndId second)
{
    return (first.fitness < second.fitness);
}

int Evoga::mrand(int i)
{
    int r;

    r = rand();
    r = r % (int) i;
    return r;

}

Evoga::Evoga()
    : exp(NULL)
    , popSize(1)
    , glen(0)
    , elitism(false)
    , nreproducing(1)
    , noffspring(1)
    , genome()
    , bestgenome()
    , loadedIndividuals(1)
    , tfitness(NULL)
    , ntfitness(NULL)
    , statfit(NULL)
    , evolutionType()
    , nogenerations(0)
    , nreplications(1)
    , seed(1)
    , currentSeed(1)
    , mutation(0.05)
    , mutationdecay(0.01)
    , initial_mutation(0.05)
    , mutations(NULL)
    , cgen(0)
    , savebest(1)
    , fmin(0.0)
    , fmax(0.0)
    , faverage(0.0)
    , fbest(0.0)
    , fbestgen(0)
    , ccycle(0)
    , stopEvolution(false)
    , isStepByStep(false)
    , mutexStepByStep()
    , waitForNextStep()
    , numThreads(1)
    , savePopulationEachNGenerations(0)
    , averageIndividualFitnessOverGenerations(true)
    , savedConfigurationParameters()
    , savedExperimentPrefix()
{
}

Evoga::~Evoga()
{
    delete exp;
    delete[] tfitness;
    delete[] terror;
    if (statfit != NULL) {
        for(int i = 0; i < nogenerations; i++) {
            delete[] statfit[i];
        }
    }
    delete[] statfit;
    delete[] ntfitness;
    delete[] mutations;
}

void Evoga::setSeed(int s)
{
    srand(s);
    globalRNG->setSeed( s );
    currentSeed = s;
}

//return a random value between 0-1
double Evoga::drand()
{
    return (double)rand()/(double)RAND_MAX;
}

double Evoga::getNoise(double minn, double maxn)
{
    double nrange;
    if(maxn>minn)
        nrange=maxn-minn;
    else
        nrange=minn-maxn;

    return drand()*nrange+minn;
}

int Evoga::mutate(int w, double mut)
{
    int b[8];
    int val;
    int ii;

    val = w;
    for(ii=0;ii < 8;ii++) {
        b[ii] = val % 2;
        val  = val / 2;
    }
    for(ii=0;ii < 8;ii++) {
        //if (mrand(100) < percmut)
        if(drand()<mut) {
            b[ii] =(b[ii]+1)%2; // con questa modifica il bit switcha //mrand(2);
        }
    }
    w = 0;
    w += b[0] * 1;
    w += b[1] * 2;
    w += b[2] * 4;
    w += b[3] * 8;
    w += b[4] * 16;
    w += b[5] * 32;
    w += b[6] * 64;
    w += b[7] * 128;

    return(w);
}

void Evoga::putGenome(int fromgenome, int tobestgenome)
{
    if (tobestgenome < this->nreproducing) {
        for (int i = 0; i < this->glen; i++) {
            bestgenome[tobestgenome][i] = genome[fromgenome][i];
        }
    } else {
        Logger::error("putGenomeError!");
    }
}

void Evoga::getGenome(int frombestgenome, int togenome, int mut)
{
    for(int i = 0; i < this->glen; i++) {
        if (mut == 0) {
            genome[togenome][i] = bestgenome[frombestgenome][i];
        } else {
            if (mutations[i] == Evonet::DEFAULT_VALUE) {    //standard mutation
                genome[togenome][i] = mutate(bestgenome[frombestgenome][i], mutation);
            } else {                //specific mutation
                genome[togenome][i] = mutate(bestgenome[frombestgenome][i], mutations[i]);
            }
        }
    }
}

void Evoga::copyGenes(int from, int to, int mut)
{
    for(int i = 0; i < this->glen; i++) {
        if (mut == 0) {
            genome[to][i] = genome[from][i];
        } else {
            if (mutations[i] == Evonet::DEFAULT_VALUE) {    //standard mutation
                genome[to][i] = mutate(genome[from][i], mutation);
            } else {                //specific mutation
                genome[to][i] = mutate(genome[from][i], mutations[i]);
            }
        }
    }
}

//Reproduce individuals with higher ranking
void Evoga::reproduce()
{
    //to do
    //selecting best fathers
    int i;
    int bi,bx;
    double bn;

    char sbuffer[64];
    FILE *fp;

    //first of all we compute fitness stat
    this->computeFStat();


    for(bi=0;bi<this->nreproducing;bi++) {
        bn=-9999.0;
        bx=-1; //individual to be copied
        for(i=0;i< this->popSize;i++) {
            if(tfitness[i]>bn) {
                bn=tfitness[i];
                bx=i;
            }
        }

        this->putGenome(bx,bi);

        //here we save best genome
        if ((bi+1)<=this->savebest && cgen< this->nogenerations) {
            sprintf(sbuffer,"B%dS%d.gen",bi+1,this->currentSeed);
            if (cgen==0) {
                if ((fp=fopen(sbuffer, "w")) == NULL) {
                    Logger::error(QString("I cannot open file B%1S%2.gen").arg(bi+1).arg(this->currentSeed));
                }
            } else {
                if ((fp=fopen(sbuffer, "a")) == NULL) {
                    Logger::error(QString("I cannot open file B%1S%2.gen").arg(bi+1).arg(this->currentSeed));
                }
            }

            fprintf(fp,"**NET : s%d_%d.wts\n",cgen,bx);
            saveagenotype(fp,bx);
            fflush(fp);
            fclose(fp);
        }
        tfitness[bx]=-9999.0;
    }

    //reproducing best
    bx=0;
    for(bi=0;bi<this->nreproducing;bi++)
        for(i =0;i< this->noffspring;i++) {
            if(this->elitism && bi==0)
                this->getGenome(bi,bx,0);
            else
                this->getGenome(bi,bx,1);// 1 mette mutazione

            bx++;
        }

    //resetting fitness
    for (i=0;i<this->popSize;i++) tfitness[i]=0.0;
    this->saveFStat();
    cgen++;
}

void Evoga::saveBestInd()
{
    //to do
    //selecting best fathers
    int i;
    int bi;
    double bn, ffit;
    bn=-999999.0;// it was 0 possible source of bug in case of negsative fitness
    bi=-1;

    char sbuffer[64];
    FILE *fp;

    sprintf(sbuffer,"B%dS%d.gen",bi+1,this->currentSeed);
    if (cgen==0) {
        if ((fp=fopen(sbuffer, "w")) == NULL) {
            Logger::error(QString("I cannot open file B%1S%2.gen").arg(bi+1).arg(this->currentSeed));
        }
    } else {
        if ((fp=fopen(sbuffer, "a")) == NULL) {
            Logger::error(QString("I cannot open file B%1S%2.gen").arg(bi+1).arg(this->currentSeed));
        }
    }

    //finding the best simply the best, one individual
    for(i=0;i<this->popSize;i++) {
        ffit=this->tfitness[i]/this->ntfitness[i];
        if(ffit>bn) {
            bn=ffit;
            bi=i;
        }
    }

    //now saving
    fprintf(fp,"**NET : s%d_%d.wts\n",cgen,bi);
    saveagenotype(fp,bi);
    fflush(fp);
    fclose(fp);
}

void Evoga::saveBestTeam(QVector< QVector<int> > teams, QVector<double> fitness)
{
    //to do
    //selecting best fathers
    int i;
    int bi;
    double bn, ffit;
    bn=-999999.0;// it was 0 possible source of bug in case of negsative fitness
    bi=-1;
    double max;
    int indMax;

    char sbuffer[64];
    FILE *fp;

    sprintf(sbuffer,"B%dS%d.G%d.gen",bi+1,this->currentSeed,cgen);
    if ((fp=fopen(sbuffer, "w")) == NULL) {
        Logger::error(QString("I cannot open file B%1S%2.G%3.gen").arg(bi+1).arg(this->currentSeed).arg(cgen));
    }

    max = fitness[0];
    indMax = 0;
    for(int i=1;i<fitness.size();i++){
        if(fitness[i]>max){
            max = fitness[i];
            indMax = i;
        }
    }

    for(int i=0;i<numModules;i++){
        //now saving
        fprintf(fp,"**NET : s%d_%d.wts\n",cgen,teams[indMax][i]);
        saveagenotype(fp,teams[indMax][i]);
    }
    fclose(fp);
    fflush(fp);
}


//Reproduce individuals with higher ranking
void Evoga::mreproduce()
{
    //to do
    //selecting best fathers
    int i;
    int bi,bx;
    double bn;

    char sbuffer[64];
    FILE *fp;

    //first of all we compute fitness stat
    this->computeFStat();


    for(bi=0;bi<this->nreproducing;bi++) {
        bn=9999.0;
        bx=-1; //individual to be copied
        for(i=0;i< this->popSize;i++) {
            if(tfitness[i]<bn) {
                bn=tfitness[i];
                bx=i;
            }
        }

        this->putGenome(bx,bi);

        //here we save best genome
        if ((bi+1)<=this->savebest && cgen< this->nogenerations) {
            sprintf(sbuffer,"B%dS%d.gen",bi+1,this->currentSeed);
            if (cgen==0) {
                if ((fp=fopen(sbuffer, "w")) == NULL) {
                    Logger::error(QString("I cannot open file B%1S%2.gen").arg(bi+1).arg(this->currentSeed));
                }
            } else {
                if ((fp=fopen(sbuffer, "a")) == NULL) {
                    Logger::error(QString("I cannot open file B%1S%2.gen").arg(bi+1).arg(this->currentSeed));
                }
            }

            fprintf(fp,"**NET : s%d_%d.wts\n",cgen,bx);
            saveagenotype(fp,bx);
            fflush(fp);
            fclose(fp);
        }
        tfitness[bx]=9999.0;
    }

    //reproducing best
    bx=0;
    for(bi=0;bi<this->nreproducing;bi++)
        for(i =0;i< this->noffspring;i++) {
            if(this->elitism && bi==0)
                this->getGenome(bi,bx,0);
            else
                this->getGenome(bi,bx,1);

            bx++;
        }

    //resetting fitness
    for (i=0;i<this->popSize;i++) tfitness[i]=9999.0;
    this->saveFStat();
    cgen++;
}


void Evoga::printPop()
{
    for(int i = 0; i < this->popSize; i++) {
        QString output = QString("Fit %1 | ").arg(tfitness[i]);
        for(int l = 0; l < this->glen; l++) {
            output += QString("%1 ").arg(this->genome[i][l]);
        }
        Logger::info(output);
    }
}

void Evoga::printBest()
{
    for(int i = 0; i < bestgenome.size(); i++) {
        QString output = QString("Best %d | ").arg(i);
        for (int s = 0; s < this->glen; s++) {
            output += QString("%1 ").arg(this->bestgenome[i][s]);
        }
        Logger::info(output);
    }
}

void Evoga::computeFStat()
{
    int i;
    double min, max, av;

    min=max=tfitness[0];
    av=0.0;

    for(i=0;i<this->popSize;i++) {
        if(tfitness[i]<min) min=tfitness[i];
        if(tfitness[i]>max) max=tfitness[i];
        av+=tfitness[i];
    }
    this->faverage=av/(double)this->popSize;
    this->fmax=max;
    this->fmin=min;
    this->statfit[this->cgen][0]=this->fmax;
    this->statfit[this->cgen][1]=this->faverage;
    this->statfit[this->cgen][2]=this->fmin;

    if (this->fmax > this->fbest) {
        this->fbest = this->fmax;
        this->fbestgen = this->cgen;
    }
}

void Evoga::computeFStat2()
{
    int i;
    double min, max, av;

    //min=max=tfitness[0]/ntfitness[0];
    //try to fix a problem
    min=9999.00;
    max=-9999.00;
    av=0.0;

    for(i=0;i<this->popSize;i++) {
        if((tfitness[i]/ntfitness[i])<min) {
            min=tfitness[i]/ntfitness[i];
        }
        if((tfitness[i]/ntfitness[i])>max) max=tfitness[i]/ntfitness[i];
        av+=(tfitness[i]/ntfitness[i]);
    }
    this->faverage=av/(double)this->popSize;
    this->fmax=max;
    this->fmin=min;
    this->statfit[this->cgen][0]=this->fmax;
    this->statfit[this->cgen][1]=this->faverage;
    this->statfit[this->cgen][2]=this->fmin;

    if (this->fmax > this->fbest) {
        this->fbest = this->fmax;
        this->fbestgen = this->cgen;
    }
}

void Evoga::saveagenotype(FILE *fp, int ind)
{
    int j;
    fprintf(fp, "DYNAMICAL NN\n");
    for (j=0; j < this->glen; j++)
        fprintf(fp, "%d\n", this->genome[ind][j]);
    fprintf(fp, "END\n");
}

//save all current generation
void Evoga::saveallg()
{
    FILE *fp;
    char filename[64];
    int i;

    sprintf(filename,"G%dS%d.gen",cgen,currentSeed);
    if ((fp=fopen(filename, "w+")) == NULL) {
        Logger::error(QString("Cannot open file %1").arg(filename));
    } else {
        //we save
        for(i=0;i<this->popSize;i++) {
            fprintf(fp,"**NET : %d_%d_%d.wts\n",cgen,0,i);
            this->saveagenotype(fp,i);
        }
        fclose( fp );
    }
}

void Evoga::saveallgComposed(QVector< QVector<int> > composedGen)
{
    FILE *fp;
    char filename[64];
    int i;

    sprintf(filename,"G%dS%d.composed.gen",cgen,currentSeed);
    if ((fp=fopen(filename, "w+")) == NULL) {
        Logger::error(QString("Cannot open file %1").arg(filename));
    } else {
        //we save
        for(i=0;i<composedGen.size();i++) {
            fprintf(fp,"**TEAM : %d_%d_%d.wts\n",cgen,0,i);
            for(int j=0;j<numModules;j++)
                fprintf(fp,"%d ",composedGen[i][j]);

            fprintf(fp, "\nEND\n");
        }
        fclose( fp );
    }
}


void Evoga::saveFStat()
{
    FILE *fp;
    char sbuffer[128];
    sprintf(sbuffer,"statS%d.fit",currentSeed);
        if (cgen == 0)
        fp=fopen(sbuffer , "w");
    else
        fp=fopen(sbuffer , "a");

    if (fp != NULL) {
        fprintf(fp,"%.3f %.3f %.3f\n",fmax,faverage,fmin);
        fclose(fp);
    } else
        Logger::error("unable to save statistics on a file");
}

void Evoga::saveRStat(QVector<int> subsVec)
{
    FILE *fp;
    char sbuffer[128];
    sprintf(sbuffer,"statS%d.ret",currentSeed);

    if (cgen == 0)
        fp=fopen(sbuffer , "w");
    else
        fp=fopen(sbuffer , "a");

    if(fp!=NULL){
        for(int i=0;i<subsVec.size();i++){
            fprintf(fp,"%i ",subsVec[i]);
        }
        fprintf(fp,"\n");
        fclose(fp);
    } else
        Logger::error("unable to save statistics of retentions on a file");
}

void Evoga::getLastFStat( double &min, double &max, double &average ) {
    min = fmin;
    max = fmax;
    average = faverage;
}

void Evoga::loadgenotype(FILE *fp, int ind)
{
    int j;
    int v;

    fscanf(fp, "DYNAMICAL NN\n");
    for (j=0; j <this->glen; j++) {
        fscanf(fp,"%d\n",&v);//this->genome[ind][j]);
        this->genome[ind][j]=v;
    }
    fscanf(fp, "END\n");
}

int Evoga::loadallg(int gen, const char *filew)
{
    FILE *fp;
    char filename[512];//[64];
    char message[512];//[128];
    char flag[512];//[64];

    if (gen >= 0) {
        sprintf(filename, "G%dS%d.gen", gen, seed);
    } else {
        sprintf(filename, "%s", filew);//in case of B%P$S.gen
    }

    if ((fp = fopen(filename, "r")) != NULL) {
        genome.clear();
        while (true) {
            flag[0] = '\0'; //sprintf(flag,""); //flag = NULL;
            fscanf(fp, "%s : %s\n", flag, message);
            if (strcmp(flag, "**NET") == 0) {
                loadgenotype(fp, genome.addOne());
            } else {
                break;
            }
        }
        Logger::info(QString("Loaded ind: %1").arg(genome.size()));
        fclose(fp);
    } else {
        Logger::error(QString("File %1 could not be opened").arg(filename));
    }

    loadedIndividuals = genome.size();

    return genome.size();
}

void Evoga::loadallTeams(int gen, const char *filew, QVector< QVector<int> > &teams)
{
    FILE *fp;
    char filename[512];//[64];
    char message[512];//[128];
    char flag[512];//[64];

    if (gen >= 0) {
        sprintf(filename, "G%dS%d.composed.gen", gen, seed);
    } else {
        sprintf(filename, "%s", filew);//in case of B%P$S.gen
    }

    if ((fp = fopen(filename, "r")) != NULL) {
        int teamInd=0;
        while (true) {
            if(teams.size()<teamInd+1){
                teams.resize(teamInd+1);
                teams[teamInd].resize(numModules);
            }
            flag[0] = '\0'; //sprintf(flag,""); //flag = NULL;
            fscanf(fp, "%s : %s\n", flag, message);
            if (strcmp(flag, "**TEAM") == 0) {
                for(int i=0;i<numModules;i++){

                    fscanf(fp,"%i",&teams[teamInd][i]);

                    int w=0;
                }
                fscanf(fp,"\n");
                fscanf(fp, "END\n");
                teamInd++;
            } else {
                break;
            }
        }
        Logger::info(QString("Loaded teams: %1").arg(teamInd));
        fclose(fp);
    } else {
        Logger::error(QString("File %1 could not be opened").arg(filename));
    }

}

/*
 * load a .fit file (return the number loaded individuals, 0 if the file does not exists)
 */
int Evoga::loadStatistics(char *filename)
{
    FILE *fp;
    int loaded=0;
    int i=0;
    float max,av,min;
    max=min=av=-1;
    if ((fp=fopen(filename, "r")) != NULL) {
        while(fscanf(fp,"%f %f %f\n",&max,&av,&min)!=EOF) {
            this->statfit[i][0]=max;
            this->statfit[i][1]=av;
            this->statfit[i][2]=min;
            i++;
        }
        loaded=i;
        fflush(fp);
        fclose(fp);
        return(loaded);
    } else {
        return(0);
    }
}

void Evoga::randomizePop()
{
    for (int i = 0; i < genome.size(); i++) {
        for(int g = 0; g < glen; g++) {
            genome[i][g] = mrand(256);
        }
    }
}

void Evoga::setInitialPopulation(int* ge)
{
    int i,g;

    //fill genome with .phe parameters
    for(i=0; i<genome.size(); i++)
        for(g=0; g<glen; g++)
        {
            if(ge[g] == Evonet::DEFAULT_VALUE)
                genome[i][g] = mrand(256);
            else
                genome[i][g] = ge[g];
        }
}

void Evoga::setMutations(float* mut)
{
    //fill mutation vector
    for(int i=0; i<glen; i++)
        mutations[i] = mut[i];
}

int* Evoga::getGenes(int ind)
{
#if defined(__GNUC__) && defined(DEVELOPER_WARNINGS)
    #warning CONTROLLARE INDICE? ALLE VOLTE CRASHA SE NELLA GUI SI SELEZIONA UN INDICE OLTRE IL MASSIMO
#endif
    return this->genome[ind];
}

int* Evoga::getBestGenes(int ind)
{
    return this->bestgenome[ind];
}

void Evoga::resetGenerationCounter()
{
    this->cgen=0;
}

void Evoga::getPheParametersAndMutationsFromEvonet()
{
    ResourcesLocker locker(this);

    Evonet* evonet = getResource<Evonet>( "evonet" );
    if (evonet->pheFileLoaded()) {
        //temporary mutations vector
        Evonet* evonet = getResource<Evonet>( "evonet" );
        float* muts = new float[evonet->freeParameters()];
        int* pheParams = new int[evonet->freeParameters()];

        //setting the mutation vector
        evonet->getMutations(muts); //taking it from the net
        setMutations(muts);     //pushing it inside GA

        //setting initial genome
        evonet->copyPheParameters(pheParams);   //copy *phe parameters
        setInitialPopulation(pheParams);            //put *phe parameters inside genome

        delete[] muts;
        delete[] pheParams;
    }
}

void Evoga::updateGenomeFromEvonet(int ind) {

    ResourcesLocker locker(this);
    Evonet* evonet = getResource<Evonet>( "evonet" );

    float wrange = evonet->getWrange();

    for (int i = 0; i<glen; i++) {
        genome[ind][i] = -(evonet->getFreeParameter(i)-wrange)/(2*wrange)*255;
    }
    /*
    printf("%s : %d : ",__PRETTY_FUNCTION__,ind);
    for (int i = 0; i<glen; i++) {
        printf("%d ",genome[ind][i]);
    }
    printf("\n");
    */
}

QString Evoga::getEvolutionType()
{
    return evolutionType;
}

float* Evoga::xnesGetGenes(int ind)
{
    return genotypes[ind];
}

void Evoga::saveBestFitness()
{
    FILE *fp;
    char sbuffer[128];
    sprintf(sbuffer, "bestgenS%d.fit", currentSeed);
    if (cgen == 0)
    {
        fp = fopen(sbuffer, "w");
    }
    else
    {
        fp = fopen(sbuffer, "a");
    }

    if (fp != NULL)
    {
        fprintf(fp, "%d %.3f\n", fbestgen, fbest);
        fclose(fp);
    }
    else
    {
        Logger::error("unable to save best generation statistics on a file");
    }
}

void Evoga::xnesComputeStats()
{
    int i;
    double min, max, av;

    min = 9999.00;
    max = -9999.00;
    av = 0.0;

    for(i = 0; i < genotypes.size(); i++)
    {
        if((tfitness[i] / ntfitness[i]) < min)
        {
            min = tfitness[i] / ntfitness[i];
        }
        if((tfitness[i] / ntfitness[i]) > max)
        {
            max = tfitness[i] / ntfitness[i];
        }
        av += (tfitness[i] / ntfitness[i]);
    }
    this->faverage = av / (double)genotypes.size();
    this->fmax = max;
    this->fmin = min;
    this->statfit[this->cgen][0] = this->fmax;
    this->statfit[this->cgen][1] = this->faverage;
    this->statfit[this->cgen][2] = this->fmin;

    if (this->fmax > this->fbest) {
        this->fbest = this->fmax;
        this->fbestgen = this->cgen;
    }
}

void Evoga::xnesSaveStats()
{
    FILE *fp;
    char sbuffer[128];
    sprintf(sbuffer, "statS%d.fit", currentSeed);
    if (cgen == 0)
    {
        fp = fopen(sbuffer, "w");
    }
    else
    {
        fp = fopen(sbuffer, "a");
    }

    if (fp != NULL)
    {
        fprintf(fp, "%.3f %.3f %.3f\n", fmax, faverage, fmin);
        fclose(fp);
    }
    else
    {
        Logger::error("unable to save statistics on a file");
    }
}

void Evoga::xnesSaveGen(FILE *fp, float* genotype)
{
    int j;
    fprintf(fp, "DYNAMICAL NN\n");
    for (j = 0; j < this->glen; j++)
    {
        fprintf(fp, "%lf\n", genotype[j]);
    }
    fprintf(fp, "END\n");
}

void Evoga::xnesSaveAllGen()
{
    FILE *fp;
    char filename[64];
    int i;

    sprintf(filename, "G%dS%d.gen", cgen, currentSeed);
    if ((fp = fopen(filename, "w+")) == NULL)
    {
        Logger::error(QString("Cannot open file %1").arg(filename));
    }
    else
    {
        //we save
        for(i = 0; i < genotypes.size(); i++)
        {
            fprintf(fp, "**NET : %d_%d_%d.wts\n", cgen, 0, i);
            this->xnesSaveGen(fp, genotypes[i]);
        }
        fclose(fp);
    }
}

void Evoga::xnesSaveBest(int bi)
{
    char sbuffer[64];
    FILE *fp;

    sprintf(sbuffer, "B%dS%d.gen", 0, this->currentSeed);
    if (cgen==0)
    {
        if ((fp = fopen(sbuffer, "w")) == NULL)
        {
            Logger::error(QString("I cannot open file B%1S%2.gen").arg(0).arg(this->currentSeed));
        }
    }
    else
    {
        if ((fp=fopen(sbuffer, "a")) == NULL)
        {
            Logger::error(QString("I cannot open file B%1S%2.gen").arg(0).arg(this->currentSeed));
        }
    }

    //now saving
    fprintf(fp, "**NET : s%d_%d.wts\n", cgen, bi);
    xnesSaveGen(fp, genotypes[bi]);
    fflush(fp);
    fclose(fp);
}

void Evoga::xnesSaveBestPhe(int gen, int bi)
{
    char sbuffer[64];
    sprintf(sbuffer, "BestPheG%dS%d.phe", gen, this->currentSeed);
    exp->getNeuralNetwork()->save_net_blocks(sbuffer, 1);
}

void Evoga::xnesLoadGen(FILE *fp, int ind)
{
    int j;
    float v;

    fscanf(fp, "DYNAMICAL NN\n");
    for (j = 0; j < this->glen; j++)
    {
        fscanf(fp, "%f\n", &v);
        //printf("Read %lf\n", v);
        genotypes[ind][j] = v;
    }
    fscanf(fp, "END\n");
}

int Evoga::xnesLoadAllGen(int gen, const char *filew)
{
    FILE *fp;
    char filename[512];
    char message[512];
    char flag[512];
    int index;

    if (gen >= 0)
    {
        sprintf(filename, "G%dS%d.gen", gen, seed);
    }
    else
    {
        sprintf(filename, "%s", filew);//in case of B%P$S.gen
    }

    if ((fp = fopen(filename, "r")) != NULL)
    {
        genotypes.clear();
        while (true)
        {
            //printf("size of genotypes: %d\n", genotypes.size());
            flag[0] = '\0';
            fscanf(fp, "%s : %s\n", flag, message);
            if (strcmp(flag, "**NET") == 0)
            {
                genotypes.append(new float[glen]);
                index = genotypes.size() - 1;
                xnesLoadGen(fp, index);
            }
            else
            {
                break;
            }
        }
        Logger::info(QString("Loaded ind: %1").arg(index));
        fclose(fp);
    }
    else
    {
        Logger::error(QString("File %1 could not be opened").arg(filename));
    }

    loadedIndividuals = genotypes.size();

    return genotypes.size();
}

/*
 * Main function of the xNES algorithm (exponential Natural Evolution Strategies)
 */
void Evoga::xnes()
{
    // ALL VARIABLES ARE INITIALIZED TO AVOID NUMERICAL PROBLEMS

    int rp; //replication
    int gn; //generation
    const float genLen = (const float)glen; // Floating point genome length
    float offspringSize = 4.0 + floor(3.0 * log(genLen)); // Floating point number of offspring generated by the individual
    // If xnes must be combined with backpropagation in the experimental way,
    // add an element to the offspring population (the one obtained by applying backpropagation)
    if (xnesCombination)
    {
        if (xnesCombinationType > 0)
        {
            offspringSize += 1.0;
        }
    }
    const int offSize = (const int)floor(offspringSize); // Number of offspring generated by the individual
    Eigen::VectorXf ind(glen); // The individual from which the offspring will be generated
    Eigen::VectorXf dInd(glen); // The update of the individual
    Eigen::MatrixXf covMatrix(glen,glen); // Covariance Matrix
    Eigen::MatrixXf dCovM(glen,glen); // The update of the covariance matrix
    Eigen::MatrixXf expCovM(glen,glen); // Exponential matrix of CovM
    float etaInd; // Learning rate for the individual
    float etaCovM; // Learning rate for the covariance matrix
    Eigen::VectorXf utility(offSize); // Array containing the utility values
    float utilSum; // Sum of utility values
    float tmp; // Temporary variable
    Eigen::VectorXf orderedUtility(offSize); // Array containing the ordered utilities (i.e. the best individual has the highest utility value)
    float ordUtilSum; // Sum of ordered utility values (it shall be equal to 1.0)
    QVector<float> f(offSize); // Array containing the fitness of the offspring
    Eigen::VectorXf ordf(offSize); // Array containing the ordered fitness of the offspring
    QVector<int> ordoff(offSize); // Array containing the ordered indices of the offspring
    float* currOffspring = new float[glen]; // Current offspring
    Eigen::MatrixXf eye(glen,glen); // Identity matrix
    // Indices
    int i = 0;
    int j = 0;
    // Start generation index
    int startGeneration;
    // File descriptors
    char statfile[128];
    char genFile[128];
    char filename[64];
    // Default mutation range (it is used only if the <useGaussian> flag is set to false;
    // it defines the default interval for mutation values).
    const double defaultMutationRange = 0.5;
    // Use <variableStdDev> flag to se the initial value
    bool startVariableStdDev;
    double startOffspringStdDev = offspringStdDev;
    // Decay of standard deviation used to generate offspring
    // (this variable is used to reduce the fitness search space)
    const double stdDevDecay = 0.05;
    // Matrix of freeps (it is needed to compute and save the best phenotype)
    float** offspringFreep = new float*[offSize];
    for (j = 0; j < offSize; j++)
    {
        offspringFreep[j] = new float[glen];
        for (i = 0; i < glen; i++)
        {
            offspringFreep[j][i] = 0.0;
        }
    }
    // Array of backpropagation weight update
    float* backpropfreep = new float[glen];
    float* oldbackpropfreep = new float[glen];
    // File descriptors
    FILE *fp;
    FILE *backpropfp;
    FILE *offspringfp;
    Eigen::VectorXf oldInd(glen); // The individual from which the offspring will be generated
    for (i = 0; i < glen; i++)
    {
        // Initialize the old individual
        oldInd[i] = 0.0;
        // Initialize the backpropagation update and its old value
        backpropfreep[i] = 0.0;
        oldbackpropfreep[i] = 0.0;
    }
    Eigen::MatrixXf parents(numStartInd,glen); // To avoid seed dependency, try different parents and choose the best one
    float* parentFitness = new float[numStartInd];

    for (j = 0; j < numStartInd; j++)
    {
        // Initialize parents and their fitnesses
        parentFitness[j] = 0.0;
        for (i = 0; i < glen; i++)
        {
            parents(j,i) = 0.0;
        }
    }

    // Flags whether to print value for debugging
    const bool debugPrint = debugPrintInfoToFile;

    // Data for the last generation
    float* lastGenInd = new float[glen];
    for (int i = 0; i < glen; i++)
    {
        lastGenInd[i] = 0.0;
    }
    float lastGenIndFit;

    // Resize genotypes
    genotypes.resize(offSize);

    Logger::info("EVOLUTION: xNES");
    Logger::info("Number of replications: " + QString::number(nreplications));

    startGeneration = 0;

    // Initialize XNES parameters
    etaInd = 1.0;
    if (standardNes)
    {
        // Standard value (i.e. the one described in the paper)
        etaCovM = (9.0 + 3.0 * log(genLen)) / (5.0 * genLen * sqrt(genLen));
    }
    else
    {
        // Value retrieved from MATLAB implementation
        etaCovM = 0.25;
        if (etaCovM > 1.0 / genLen)
        {
            etaCovM = 1.0 / genLen;
        }
        etaCovM *= 0.5;
    }
    utilSum = 0.0;
    // Calculate the utility values before normalization
    for (j = 0; j < offSize; j++)
    {
        // u[j] is equal to the uttermost between 0 and (log(lambda / 2 + 1) - log(j))
        tmp = log(offspringSize / 2.0 + 1.0) - log((float)j + 1.0); // Remember that log(0) = -Inf, so 1 is added to index i
        if (tmp < 0.0)
        {
            tmp = 0.0;
        }
        utility[j] = tmp;
        // Update the sum of utility values
        utilSum += utility[j];
    }
    // Normalize the utility values
    for (j = 0; j < offSize; j++)
    {
        if (utilSum != 0.0)
        {
            utility[j] /= utilSum;
        }
        if (standardNes)
        {
            utility[j] -= 1.0 / genLen; // In the standard version, utility values might be negative
        }
    }

    for (i = 0; i < glen; i++)
    {
        // Define the identity matrix
        eye(i,i) = 1.0;
        for (j = 0; j < glen; j++)
        {
            if (i != j)
            {
                eye(i,j) = 0.0;
            }
            // Initialize covariance matrix, its update and its exponential
            covMatrix(i,j) = 0.0;
            dCovM(i,j) = 0.0;
            expCovM(i,j) = 0.0;
        }
        // Initialize individual and its update
        ind[i] = 0.0;
        dInd[i] = 0.0;
        // Initialize the current offspring
        currOffspring[i] = 0.0;
    }

    if (!useGaussian)
    {
        // Set mutation range to default value if it is equal to 0
        // (this parameter is used only in case the <useGaussian>
        // flag is set to false)
        if (mutationRange == 0.0)
        {
            mutationRange = defaultMutationRange;
        }
        if (offspringMutRange == 0.0)
        {
            offspringMutRange = defaultMutationRange;
        }
    }

    // Creating evaluator objects in case of a multithread simulation. Also setting the actual number of threads used
    QVector<EvaluatorThreadForXnes*> evaluators(offSize, NULL);
    if (numThreads > 1) {
        for (int i = 0; i < evaluators.size(); i++) {
            EvoRobotExperiment* newExp = savedConfigurationParameters.getObjectFromGroup<EvoRobotExperiment>(savedExperimentPrefix, false);
            newExp->setEvoga(this);
            evaluators[i] = new EvaluatorThreadForXnes(this, newExp);
        }
        QThreadPool::globalInstance()->setMaxThreadCount(numThreads);
    }

    for(rp = 0; rp < nreplications; rp++) { // replications
        setSeed(getStartingSeed() + rp);
        Logger::info(QString("Replication %1, seed: %2").arg(rp + 1).arg(getStartingSeed() + rp));
        resetGenerationCounter();
        // --- section runnable only if there is an Evonet object
        if ( hasResource( "evonet" ) ) {
            getPheParametersAndMutationsFromEvonet();
        }
        emit startingReplication( rp );

        // Set fbest to a very low value
        this->fbest = -99999.0;

        // Resetting seed in experiments
        exp->newGASeed(getCurrentSeed());
        if (numThreads > 1) {
            for (int i = 0; i < evaluators.size(); i++) {
                evaluators[i]->getExperiment()->newGASeed(getCurrentSeed());
            }
        }

        QTime evotimer;
        evotimer.start();
        // Initialization of fitnesses
        for (j = 0; j < offSize; j++) {
            tfitness[j] = 0.0;
            ntfitness[j] = 0.0;
            terror[j] = 0.0;
        }

        //code to recovery a previous evolution: Experimental
        sprintf(statfile,"statS%d.fit", getStartingSeed() + rp);
        //now check if the file exists
        DataChunk statTest(QString("stattest"),Qt::blue,2000,false);
        if (statTest.loadRawData(QString(statfile),0)) {
            startGeneration=statTest.getIndex();
            sprintf(genFile,"G%dS%d.gen", startGeneration, getStartingSeed() + rp);
            Logger::info("Recovering from startGeneration: " + QString::number(startGeneration));
            Logger::info(QString("Loading file: ") + genFile);
            loadallg(-1, genFile);
            cgen = startGeneration;
            emit recoveredInterruptedEvolution( QString(statfile) );
        } //end evolution recovery code

        double oldIndVariation = 0.0;

        // Test <numStartInd> parents and find out the best one (this is the starting "individual")

        // Initialize parents
        for (j = 0; j < numStartInd; j++)
        {
            for (i = 0; i < glen; i++)
            {
                // The initialization depends on the type of random numbers to be extracted
                if (useGaussian)
                {
                    parents(j,i) = globalRNG->getGaussian(gaussianStdDev, gaussianMean);
                }
                else
                {
                    parents(j,i) = globalRNG->getDouble(-mutationRange, mutationRange);
                }
            }
        }

        // Test parents
        for (j = 0; j < numStartInd; j++)
        {
            for (i = 0; i < glen; i++)
            {
                currOffspring[i] = parents(j,i); // To avoid the definition of an additional pointer variable
            }
            // Set the ANN parameters
            exp->setNetParameters(currOffspring);
            // Perform trials for the individual (i.e. the offspring)
            exp->doAllTrialsForIndividual(j);
            // Calculate the fitness of the individual (i.e. the offspring)
            parentFitness[j] = exp->getFitness();
        }

        // Find out the best parent
        int bestParent = 0;
        double bestParentFitness = parentFitness[0];
        for (j = 1; j < numStartInd; j++)
        {
            if (parentFitness[j] > bestParentFitness)
            {
                // Update the best parent
                bestParent = j;
                bestParentFitness = parentFitness[j];
            }
        }

        // ADD CODE TO EVOLVE BEST PARENT WITH SOME STEPS OF XNES
        // N.B.: TO DO THIS, A FUNCTION FOR EXECUTING XNES COMPUTATION
        // OUTSIDE XNES MAIN FUNCTION SHALL BE DEFINED

        // Initialize the individual for this replication
        for (i = 0; i < glen; i++)
        {
            ind[i] = parents(bestParent,i);
            // Reset the current offspring (to avoid numerical issues)
            currOffspring[i] = 0.0;
        }

        // Initialize Covariance Matrix
        for (i = 0; i < glen; i++)
        {
            for (j = 0; j < glen; j++)
            {
                covMatrix(i,j) = 0.0;
                expCovM(i,j) = 0.0;
            }
        }

        // Initialize fitnesses, ordered offspring and ordered utilities
        for (j = 0; j < offSize; j++) // Remember that the last slot is reserved for the parent
        {
            f[j] = 0.0;
            ordf[j] = 0.0;
            ordoff[j] = 0.0;
            orderedUtility[j] = 0.0;
        }

        // Initialize the matrix of offspring freeps
        for (j = 0; j < offSize; j++)
        {
            for (i = 0; i < glen; i++)
            {
                offspringFreep[j][i] = 0.0;
            }
        }

        // CODE TO PRINT SOME DEBUG INFORMATION (i.e. THE INDIVIDUAL AT THE BEGINNING OF THE EVOLUTIONARY PROCESS)
        if (debugPrint)
        {
            sprintf(filename,"freepS%d.txt", currentSeed);
            if ((fp = fopen(filename, "a")) == NULL)
            {
                Logger::error(QString("I cannot open file freep.txt"));
            }
            for (i = 0; i < glen; i++)
            {
                fprintf(fp, "%.2f\t", ind[i]);
            }
            fprintf(fp, "\n");
            fflush(fp);
        }

        // Initialise <startVariableStdDev> flag
        startVariableStdDev = variableStdDev;
        offspringStdDev = startOffspringStdDev;

        for(gn = startGeneration; gn < nogenerations; gn++) {   // generations
            evotimer.restart();
            Logger::info("Generation " + QString::number(gn));
            //Logger::info("Generation " + QString::number(gn + 1));

            // Calculate the matrix exponential
            Eigen::MatrixExponential<Eigen::MatrixXf>(covMatrix).compute(expCovM);

            // Extracting samples (i.e. generate offspring of current individual)
            Eigen::MatrixXf samples(glen,offSize);
            Eigen::MatrixXf trSamples(glen,offSize);
            Eigen::MatrixXf offspring(glen,offSize);

            for (i = 0; i < glen; i++)
            {
                for (j = 0; j < offSize; j++)
                {
                    // The initialization depends on the type of random numbers to be extracted
                    if (useGaussian)
                    {
                        samples(i,j) = globalRNG->getGaussian(offspringStdDev, gaussianMean);
                    }
                    else
                    {
                        samples(i,j) = globalRNG->getDouble(-offspringMutRange, offspringMutRange);
                    }
                }
            }

            if (xnesCombination && (xnesCombinationType > 0))
            {
                for (i = 0; i < glen; i++)
                {
                    // The last slot is reserved to backpropagation offspring,
                    // thus the current sample values are useless
                    samples(i,offSize - 1) = 0.0;
                }
            }

            exp->initGeneration(gn);

            if (xnesCombination)
            {
                // Firstly test the parent

                for (i = 0; i < glen; i++)
                {
                    currOffspring[i] = ind[i]; // This statement might sound strange (I use currOffspring to avoid defining another pointer variable)
                }

                // Set the ANN parameters
                exp->setNetParameters(currOffspring);
                // Perform trials for the individual (i.e. the offspring)
                exp->doAllTrialsForIndividual(offSize - 1);

                // If it is the last generation, save the fitness of the parent
                if ((gn + 1) == nogenerations)
                {
                    for (i = 0; i < glen; i++)
                    {
                        lastGenInd[i] = exp->getNeuralNetwork()->getFreeParameter(i);
                    }
                    lastGenIndFit = exp->getFitness();
                }

                // Getting the backpropvector for the parent
                backpropfreep = exp->getNeuralNetwork()->getCopyBackPropFreep();
                float maxVal = backpropfreep[0];
                float minVal = backpropfreep[0];
                for (i = 1; i < glen; i++)
                {
                    if (backpropfreep[i] < minVal)
                    {
                        minVal = backpropfreep[i];
                    }
                    if (backpropfreep[i] > maxVal)
                    {
                        maxVal = backpropfreep[i];
                    }
                }
                // Determine which is the normalization factor for the array
                float normFactor = 1.0;
                float factor;
                if (fabs(minVal) > fabs(maxVal))
                {
                    factor = fabs(minVal);
                }
                else
                {
                    factor = fabs(maxVal);
                }
                // Compute the normalization factor for the array
                // N.B.: if factor is outside the range, <normFactor>
                // restores backprop values in the range. If factor is
                // inside the range, <normFactor> increases backprop
                // values so that at least one entry is set to range
                normFactor = fabs(mutationLearningRate) / factor;
                // Check the amplitude of the backpropagation weight update
                for (i = 0; i < glen; i++)
                {
                    backpropfreep[i] = backpropfreep[i] * normFactor;
                }

                // CODE TO PRINT DEBUG INFORMATION ABOUT THE BACKPROPAGATION OFFSPRING
                if (debugPrint)
                {
                    sprintf(filename,"backpropfreepS%d.txt", currentSeed);
                    if ((backpropfp = fopen(filename, "a")) == NULL)
                    {
                        Logger::error(QString("I cannot open file backpropfreep.txt"));
                    }
                    double backpropoffset = 0.0;
                    for (i = 0; i < glen; i++)
                    {
                        backpropoffset += fabs(backpropfreep[i] - oldbackpropfreep[i]);
                        fprintf(backpropfp, "%.2f\t", backpropfreep[i]);
                        oldbackpropfreep[i] = backpropfreep[i];
                    }
                    fprintf(backpropfp, "\n");
                    fprintf(backpropfp, "offset %.2f \n", backpropoffset);
                    fflush(backpropfp);
                }

                if (isStopped()) { // stop evolution
                    return;
                }

                if (xnesCombinationType > 0)
                {
                    // Stores the backpropagation vector in the <samples> matrix
                    for (i = 0; i < glen; i++)
                    {
                        samples(i,offSize - 1) = backpropfreep[i];
                    }
                    // Debug information to display the offspring and the range of values
                    // (it is useless for comparison purposes)
                    if (debugPrint)
                    {
                        sprintf(filename,"offspringS%d.txt", currentSeed);
                        if ((offspringfp = fopen(filename, "a")) == NULL)
                        {
                            Logger::error(QString("I cannot open file offsrping.txt"));
                        }
                        fprintf(offspringfp, "generation %d\n", gn);
                        for (j = 0; j < offSize; j++)
                        {
                            double sampleSum = 0.0;
                            double maxVal = -99999.0;
                            double minVal = 99999.0;
                            for (i = 0; i < glen; i++)
                            {
                                fprintf(offspringfp, "%.2f\t", samples(i,j));
                                sampleSum += samples(i,j);
                                if (samples(i,j) > maxVal)
                                {
                                    maxVal = samples(i,j);
                                }
                                if (samples(i,j) < minVal)
                                {
                                    minVal = samples(i,j);
                                }
                            }
                            fprintf(offspringfp, "sum of values %.2f\t", sampleSum);
                            fprintf(offspringfp, "min value %.2f\t", minVal);
                            fprintf(offspringfp, "max value %.2f\t", maxVal);
                            fprintf(offspringfp, "\n");
                        }
                        fflush(offspringfp);
                        fclose(offspringfp);
                    }
                }
            }

            // Create the offspring
            trSamples = expCovM * samples;
            for (j = 0; j < offSize; j++)
            {
                for (i = 0; i < glen; i++)
                {
                    offspring(i,j) = ind[i] + trSamples(i,j);
                }
            }

            if (numThreads <= 1) {
                // Fill the current offspring
                for (j = 0; j < offSize; j++)
                {
                    Logger::info("Offspring " + QString::number(j));
                    for (i = 0; i < glen; i++)
                    {
                        currOffspring[i] = offspring(i,j);
                        // Save offspring genotypes
                        genotypes[j][i] = currOffspring[i];
                    }
                    // Set the ANN parameters
                    exp->setNetParameters(currOffspring);
                    // Perform trials for the individual (i.e. the offspring)
                    exp->doAllTrialsForIndividual(j);
                    // Extract the freeps for the offspring j
                    for (i = 0; i < glen; i++)
                    {
                        offspringFreep[j][i] = exp->getNeuralNetwork()->getFreeParameter(i);
                    }
                    // Calculate the fitness of the individual (i.e. the offspring)
                    f[j] = exp->getFitness();
                    // Compute the error of the individual
                    terror[j] = exp->getError();

                    if (averageIndividualFitnessOverGenerations) {
                        tfitness[j] += f[j];
                        ntfitness[j]++;
                    } else {
                        tfitness[j] = f[j];
                        ntfitness[j] = 1;
                    }
                    if (isStopped()) { // stop evolution
                        return;
                    }
                }
            }
            else
            {
                // Fill the genotypes
                for (j = 0; j < offSize; j++)
                {
                    for (i = 0; i < glen; i++)
                    {
                        genotypes[j][i] = offspring(i,j);
                    }
                }
                float* cOffspring = new float[glen];
                // We first evaluate all parents, so setting genotypes of parents (we have as many evaluators as individuals)
                for (j = 0; j < offSize; j++)
                {
                    for (i = 0; i < glen; i++)
                    {
                        cOffspring[i] = offspring(i,j);
                    }
                    evaluators[j]->setGenotype(j, cOffspring);
                }

                if (commitStep()) return; // stop the evolution process

                // Now starting parallel evaluation of parents and wating for it to finish
                QFuture<void> evaluationFuture = QtConcurrent::map(evaluators, runEvaluatorThreadForXnes);
                evaluationFuture.waitForFinished();

                if (commitStep()) return; // stop the evolution process

                // Retrieve the free parameters
                for (j = 0; j < offSize; j++)
                {
                    for (i = 0; i < glen; i++)
                    {
                        offspringFreep[j][i] = evaluators[j]->getExperiment()->getNeuralNetwork()->getFreeParameter(i);
                    }
                    f[j] = evaluators[j]->getFitness();
                    // Compute the error of the individual
                    terror[evaluators[j]->getGenotypeId()] = evaluators[j]->getError();
                }

                // We have finished evaluating parents, updating the fitness vectors
                for (j = 0; j < offSize; j++)
                {
                    if (averageIndividualFitnessOverGenerations)
                    {
                        tfitness[evaluators[j]->getGenotypeId()] += f[j];
                        ntfitness[evaluators[j]->getGenotypeId()]++;
                    }
                    else
                    {
                        tfitness[evaluators[j]->getGenotypeId()] = f[j];
                        ntfitness[evaluators[j]->getGenotypeId()] = 1;
                    }
                }
                if (commitStep()) return; // stop the evolution process
            }

            if (xnesCombination && (xnesCombinationType > 0))
            {
                // Increase the fitness of the backpropagation offspring if the
                // corresponding flag is set to true
                if (backPropOffspringFitnessIncrease)
                {
                    f[offSize - 1] += backPropOffspringFitnessIncreasePercentage * f[offSize - 1];
                }
            }

            // Write to a file the offset of the fitness of the backpropagation offspring
            // with the fitness of the best offspring (different from the backpropagation)
            if (debugPrint)
            {
                char backpropFitnessOffset[64];
                FILE* bfo;
                sprintf(backpropFitnessOffset,"backpropagationFitnessOffset.txt");
                bfo = fopen(backpropFitnessOffset, "a");
                if (bfo != NULL)
                {
                    // Look for the best fitness;
                    float bestFit = f[0];
                    for (j = 1; j < offSize - 1; j++)
                    {
                        if (f[j] > bestFit)
                        {
                            bestFit = f[j];
                        }
                    }
                    const float fitOffset = bestFit - f[offSize - 1];
                    fprintf(bfo, "Best fitness\t%lf\tBack-propagation fitness\t%lf\tOffset\t%lf\n", bestFit, f[offSize - 1], fitOffset);
                    fclose(bfo);
                }
            }

            // Now fitnesses have to be ordered in descending order
            // (in the original XNES, fitnesses are ordered in ascending order)
            int startIndex = 0;
            QVector<bool> examined(offSize);
            for (j = 0; j < offSize; j++)
            {
                examined[j] = false;
            }
            int countExamined = 0;
            while (countExamined < offSize)
            {
                bool found = false;
                i = 0;
                float minf;
                int indMinf;
                float maxf;
                int indMaxf;
                float minerr;
                while (i < offSize && !found)
                {
                    if (!examined[i])
                    {
                        if (minimization)
                        {
                            minf = f[i];
                            indMinf = i;
                        }
                        else
                        {
                            maxf = f[i];
                            indMaxf = i;
                        }
                        minerr = terror[i];
                        found = true;
                    }
                    i++;
                }
                for (j = 0; j < offSize; j++)
                {
                    if (!specialUtilityRanking)
                    {
                        // Use the standard method to compute the utility ranking
                        if (minimization)
                        {
                            if (f[j] < minf && !examined[j])
                            {
                                minf = f[j];
                                indMinf = j;
                            }
                        }
                        else
                        {
                            if (f[j] > maxf && !examined[j])
                            {
                                maxf = f[j];
                                indMaxf = j;
                            }
                        }
                    }
                    else
                    {
                        // Use the special method to compute the utility ranking
                        if ((countExamined % 2) == 0)
                        {
                            // Find the best individual depending on the fitness
                            if (minimization)
                            {
                                if (f[j] < minf && !examined[j])
                                {
                                    minf = f[j];
                                    minerr = terror[j];
                                    indMinf = j;
                                }
                            }
                            else
                            {
                                if (f[j] > maxf && !examined[j])
                                {
                                    maxf = f[j];
                                    minerr = terror[j];
                                    indMaxf = j;
                                }
                            }
                        }
                        else
                        {
                            // Find the best individual depending on the similarity error
                            if (terror[j] < minerr && !examined[j])
                            {
                                if (minimization)
                                {
                                    minf = f[j];
                                    indMinf = j;
                                }
                                else
                                {
                                    maxf = f[j];
                                    indMaxf = j;
                                }
                                minerr = terror[j];
                            }
                        }
                    }
                }
                if (xnesCombination && (xnesCombinationType > 0))
                {
                    // If the combination type is 1, the backpropagation
                    // offspring is forced to be "useful" except for the
                    // case in which the backpropagation offspring is
                    // actually the best/fittest
                    if ((countExamined == backPropOffspringUtilityRank) && !examined[offSize - 1])
                    {
                        if (minimization)
                        {
                            indMinf = offSize - 1;
                        }
                        else
                        {
                            indMaxf = offSize - 1;
                        }
                    }
                }
                if (minimization)
                {
                    ordf[startIndex] = minf;
                    ordoff[startIndex] = indMinf;
                    examined[indMinf] = true;
                }
                else
                {
                    ordf[startIndex] = maxf;
                    ordoff[startIndex] = indMaxf;
                    examined[indMaxf] = true;
                }
                countExamined++;
                startIndex++;
            }

            // Initialize the sum of ordered utility values
            ordUtilSum = 0.0;
            // Now fill the array with the utility values
            for (j = 0; j < offSize; j++)
            {
                orderedUtility[ordoff[j]] = utility[j];
                ordUtilSum += orderedUtility[ordoff[j]];
            }

            // Compute the natural gradient for both the individual and the covariance matrix

            // Create a matrix whose rows contain weights * samples
            Eigen::MatrixXf ordUtil(glen,offSize);
            Eigen::MatrixXf ordUtilSamples(glen,offSize);
            Eigen::MatrixXf Z(offSize,glen);

            for (i = 0; i < glen; i++)
            {
                for (j = 0; j < offSize; j++)
                {
                    ordUtil(i,j) = orderedUtility[j];
                }
            }

            for (i = 0; i < glen; i++)
            {
                for (j = 0; j < offSize; j++)
                {
                    ordUtilSamples(i,j) = ordUtil(i,j) * samples(i,j);
                    Z(j,i) = samples(i,j);
                }
            }

            // If you want to dissociate the effect of back-propagation on the
            // update of the covariance matrix, reset the corresponding entry
            if (dissociateCovMatrix)
            {
                for (i = 0; i < glen; i++)
                {
                    ordUtilSamples(i, offSize - 1) = 0.0;
                    Z(offSize - 1, i) = 0.0;
                }
            }

            Eigen::MatrixXf G = ordUtilSamples * Z - ordUtilSum * eye;

            // If you want to dissociate the effect of back-propagation on the
            // update of the individual, reset the corresponding entry
            if (dissociateInd)
            {
                /*for (i = 0; i < glen; i++)
                {
                    samples(i, offSize - 1) = 0.0;
                }*/
                orderedUtility[offSize - 1] = 0.0;
            }

            // Define transpose matrices for both samples and ordered utilities
            // (they are necessary to perform the calculation of gradients)
            Eigen::MatrixXf ordUtilTr(offSize,1);
            Eigen::MatrixXf v(glen,1);
            for (j = 0; j < offSize; j++)
            {
                ordUtilTr(j,0) = orderedUtility[j];
            }

            v = etaInd * expCovM * samples * ordUtilTr;
            for (i = 0; i < glen; i++)
            {
                dInd[i] = v(i,0);
            }
            dCovM = etaCovM * G;
            // Update the individual and the covariance matrix
            for (i = 0; i < glen; i++)
            {
                ind[i] += dInd[i];
                if (xnesCombination && (xnesCombinationType == 0))
                {
                    ind[i] += backpropfreep[i];
                }
            }
            covMatrix += dCovM;

            exp->endGeneration(gn);

            // Debug information to check how much the individual is changed across two generations
            if (debugPrint)
            {
                double indOffset = 0.0;
                for (i = 0; i < glen; i++)
                {
                    indOffset += fabs(ind[i] - oldInd[i]);
                    fprintf(fp, "%.2f\t", ind[i]);
                    oldInd[i] = ind[i];
                }
                fprintf(fp, "offset %.2f \n", indOffset);
                fflush(fp);
            }

            if ( commitStep() ) { return; }

            // Save the best individual (in terms of its fitness) among the offspring
            xnesSaveBest(ordoff[0]);
            // Compute and save statistics
            xnesComputeStats();
            xnesSaveStats();
            // Each <pheGen> generations, save the best phenotype (i.e. the one of the best individual)
            if (pheGen > 0)
            {
                if ((gn + 1) % pheGen == 0)
                {
                    int bestInd = ordoff[0];
                    if ((gn + 1) == nogenerations)
                    {
                        // During the last generation, the individual is also tested
                        // and it might be the best
                        if (lastGenIndFit > ordf[0])
                        {
                            exp->setNetParameters(lastGenInd);
                            xnesSaveBestPhe(gn, offSize);
                        }
                        else
                        {
                            exp->setNetParameters(offspringFreep[bestInd]);
                            xnesSaveBestPhe(gn, bestInd);
                        }
                    }
                    else
                    {
                        exp->setNetParameters(offspringFreep[bestInd]);
                        xnesSaveBestPhe(gn, bestInd);
                    }
                }
            }
            else
            {
                // If <pheGen> == 0, save the best phenotype at the end of each replication
                if ((gn + 1) == nogenerations)
                {
                    int bestInd = ordoff[0];
                    if ((gn + 1) == nogenerations)
                    {
                        // During the last generation, the individual is also tested
                        // and it might be the best
                        if (lastGenIndFit > f[bestInd])
                        {
                            exp->setNetParameters(lastGenInd);
                            xnesSaveBestPhe(gn, offSize);
                        }
                    }
                    else
                    {
                        exp->setNetParameters(offspringFreep[bestInd]);
                        xnesSaveBestPhe(gn, bestInd);
                    }
                }
            }

            emit endGeneration( cgen, fmax, faverage, fmin );
            if (commitStep()) {
                return; // stop the evolution process
            }

            cgen++;

            //always save in order to be able to resume the evolution process, but keep only the last gen file unless it has to be saved by the param savePopulationEachNGenerations
            xnesSaveAllGen();

            //remove the previous genfile unless it has to be kept because of the savePopulationEachNGenerations param
            if ((savePopulationEachNGenerations == 0) || (gn>1 && ((gn-1) % (savePopulationEachNGenerations) != 0))) {
                //EX: gn 998 = G999S1.gen --- gn 999 = G1000S1.gen --- gn = 1000 = G1001S1.gen --- gn 1001 = G1002S1.gen
                sprintf(filename,"G%dS%d.gen",gn,currentSeed);
                if(remove(filename) != 0) {
                //if(!QFile::remove(filename)) {
                    Logger::warning(QString("Error deleting temporary gen file: ") + QString::fromStdString(filename));
                }
            }

            Logger::info(QString("Generation %1 took %2 minutes").arg(gn).arg((double)evotimer.elapsed()/60000.0, 0, 'f', 2));
            fflush(stdout);

            // Reset weights and fitnesses once used
            for (j = 0; j < offSize; j++)
            {
                f[j] = 0.0;
                orderedUtility[j] = 0.0;
                tfitness[j] = 0.0;
                ntfitness[j] = 0;
                terror[j] = 0.0;
            }
            // Reset matrix exponential
            for (i = 0; i < glen; i++)
            {
                for (j = 0; j < glen; j++)
                {
                    expCovM(i,j) = 0.0;
                }
            }

            // Reset current offspring
            for (i = 0; i < glen; i++)
            {
                currOffspring[i] = 0.0;
            }

            // Reset the offspring freeps
            for (j = 0; j < offSize; j++)
            {
                for (i = 0; i < glen; i++)
                {
                    offspringFreep[j][i] = 0.0;
                }
            }

            // If the <variableStdDev> flag is set to true,
            // reduce the <offspringStdDev> until it reaches
            // the lower bound of 0.1
            if (startVariableStdDev)
            {
                double indVariation = 0.0;
                for (i = 0; i < glen; i++)
                {
                    indVariation += fabs(ind[i] - oldInd[i]);
                }
                // The variation (increase/decrease) of the offspring standard deviation
                // depends on the variation of the individual across two consecutive generations:
                // - if the variation is too big, <offspringStdDev> is halved;
                // - if the variation is too small, <offspringStdDev> is doubled;
                // - if the variation is limited, <offspringStdDev> is increased/decreased by a constant decay value.
                if (indVariation <= 1.0)
                {
                    offspringStdDev *= 2.0;
                    if (offspringStdDev > 1.0)
                    {
                        offspringStdDev = 1.0;
                    }
                }
                else if (indVariation >= 10.0)
                {
                    offspringStdDev /= 2.0;
                    if (offspringStdDev < 0.1)
                    {
                        offspringStdDev = 0.1;
                    }
                }
                else
                {
                    int multiplier = (indVariation >= oldIndVariation) ? 1 : -1;
                    offspringStdDev += multiplier * stdDevDecay;
                    if (offspringStdDev > 1.0)
                    {
                        offspringStdDev = 1.0;
                    }
                    if (offspringStdDev < 0.1)
                    {
                        offspringStdDev = 0.1;
                    }
                }
                // Update the individual variation
                oldIndVariation = indVariation;
            }
        }
        xnesSaveAllGen();
        // Save the best generation fitness statistics
        saveBestFitness();
    }

    // Deleting all evaluators
    for (int i = 0; i < evaluators.size(); i++) {
        delete evaluators[i];
    }
}

/*
 * Main function of the Genetic Algorithm (Steady State Version)
 */
void Evoga::    evolveSteadyState()
{
    int rp; //replication
    int gn; //generation
    int id; //individuals
    double fit;
    double minfit=9999;
    int    minid=-1;
    double mfit;
    float  final_mrate;
    int startGeneration=0;
    char statfile[128];
    char genFile[128];
    char filename[64];
    QVector< int > subsVec;
    double currentRetentionRate;

    subsVec.resize(popSize);
    // Resizing genome
    genome.resize(popSize * 2);

    final_mrate = mutation;
    Logger::info("EVOLUTION: steady state");
    Logger::info("Number of replications: " + QString::number(nreplications));

    // Creating evaluator objects in case of a multithread simulation. Also setting the actual number of threads used
    QVector<EvaluatorThreadForEvoga*> evaluators(popSize, NULL);
    if (numThreads > 1) {
        for (int i = 0; i < evaluators.size(); i++) {
            EvoRobotExperiment* newExp = savedConfigurationParameters.getObjectFromGroup<EvoRobotExperiment>(savedExperimentPrefix, false);
            newExp->setEvoga(this);
            evaluators[i] = new EvaluatorThreadForEvoga(this, newExp);
        }
        QThreadPool::globalInstance()->setMaxThreadCount(numThreads);
    }

    for(rp=0;rp<nreplications;rp++) {   // replications
        startGeneration = 0;
        limitationFactor = 1.0;
        mutation=initial_mutation; // initially mutation (default 50%)
        //setSeed(getSeed());
        setSeed(getStartingSeed()+rp);
        Logger::info(QString("Replication %1, seed: %2").arg(rp+1).arg(getStartingSeed()+rp));
        resetGenerationCounter();
        randomizePop();
        // --- section runnable only if there is an Evonet object
        if ( hasResource( "evonet" ) ) {
            getPheParametersAndMutationsFromEvonet();
        }

        // Set fbest to a very low value
        this->fbest = -99999.0;

        emit startingReplication( rp );

        // Resetting seed in experiments
        exp->newGASeed(getCurrentSeed());
        if (numThreads > 1) {
            for (int i = 0; i < evaluators.size(); i++) {
                evaluators[i]->getExperiment()->newGASeed(getCurrentSeed());
            }
        }

        QTime evotimer;
        evotimer.start();
        for (int i=0;i<popSize+1;i++) {
            tfitness[i]=0.0;
            ntfitness[i]=0.0;
        }
        //code to recovery a previous evolution: Experimental
        sprintf(statfile,"statS%d.fit", getStartingSeed()+rp);
        //now check if the file exists
        DataChunk statTest(QString("stattest"),Qt::blue,2000,false);
        if (statTest.loadRawData(QString(statfile),0)) {
            startGeneration=statTest.getIndex();
            sprintf(genFile,"G%dS%d.gen",startGeneration,getStartingSeed()+rp);
            Logger::info("Recovering from startGeneration: " + QString::number(startGeneration));
            Logger::info(QString("Loading file: ") + genFile);
            loadallg(-1,genFile);
            cgen=startGeneration;
            mutation=mutation-startGeneration*mutationdecay;
            if (mutation<final_mrate) mutation=final_mrate;
            // Resizing genome the loading process changed the genome size = popSize
            genome.resize(popSize * 2);
            emit recoveredInterruptedEvolution( QString(statfile) );
        } //end evolution recovery code

        for(gn=startGeneration;gn<nogenerations;gn++) { // generations
            evotimer.restart();
            currentRetentionRate = 0.0;
            Logger::info(" Generation " + QString::number(gn+1));
            // Here we do this to avoid too many complications: if we have to run no threads, we use the old
            // code, otherwise we go for the multithread code below
            if (numThreads <= 1) {
                exp->initGeneration(gn);
                if ( commitStep() ) { return; }
                // Not running with multiple threads, using the old code
                for(id=0;id<popSize;id++) { //individuals
                    fit=0.0;
                    exp->setNetParameters(getGenes(id)); // get the free parameters from the genotype
                    exp->doAllTrialsForIndividual(id);
                    fit = exp->getFitness();
                    if (averageIndividualFitnessOverGenerations) {
                        tfitness[id] += fit;
                        ntfitness[id]++;
                    } else {
                        tfitness[id] = fit;
                        ntfitness[id] = 1;
                    }
                    if (isStopped()) { // stop evolution
                        return;
                    }
                    copyGenes(id, popSize+id,1); //generate a variation by duplicating and mutating
                    tfitness[popSize+id]=0;
                    ntfitness[popSize+id]=0;

                    exp->setNetParameters(getGenes(popSize)); // get the free parameters from the genotype
                    exp->doAllTrialsForIndividual(popSize + id);
                    fit = exp->getFitness();
                    if (averageIndividualFitnessOverGenerations) {
                        tfitness[popSize+id] += fit;
                        ntfitness[popSize+id]++;
                    } else {
                        tfitness[popSize+id] = fit;
                        ntfitness[popSize+id] = 1;
                    }
                    if (isStopped()) { // stop evolution
                        return;
                    }

                }
                exp->endGeneration(gn);
                if ( commitStep() ) { return; }
            } else {
                // Multithread code

                // Calling initGeneration on all evaluator
                for (int i = 0; i < popSize; i++) {
                    evaluators[i]->getExperiment()->initGeneration(gn);
                }
                if (commitStep()) return; // stop the evolution process

                // We first evaluate all parents, so setting genotypes of parents (we have as many evaluators as individuals)
                for (int i = 0; i < popSize; i++) {
                    evaluators[i]->setGenotype(i);
                }
                if (commitStep()) return; // stop the evolution process

                // Now starting parallel evaluation of parents and wating for it to finish
                QFuture<void> evaluationFuture = QtConcurrent::map(evaluators, runEvaluatorThreadForEvoga);
                evaluationFuture.waitForFinished();
                if (commitStep()) return; // stop the evolution process

                // We have finished evaluating parents, updating the fitness vectors
                for (int i = 0; i < popSize; i++) {
                    if (averageIndividualFitnessOverGenerations) {
                        tfitness[evaluators[i]->getGenotypeId()] += evaluators[i]->getFitness();
                        ntfitness[evaluators[i]->getGenotypeId()]++;
                    } else {
                        tfitness[evaluators[i]->getGenotypeId()] = evaluators[i]->getFitness();
                        ntfitness[evaluators[i]->getGenotypeId()] = 1;
                    }
                }
                if (commitStep()) return; // stop the evolution process

                // Now we can generate all children for all individuals and set them in the evaluators
                for (int i = 0; i < popSize; i++) {
                    copyGenes(i, popSize + i, 1); //generate a variation by duplicating and mutating
                    tfitness[popSize + i] = 0;
                    ntfitness[popSize + i] = 0;
                    evaluators[i]->setGenotype(popSize + i);
                }
                if (commitStep()) return; // stop the evolution process

                // Now starting parallel evaluation of children and wating for it to finish
                evaluationFuture = QtConcurrent::map(evaluators, runEvaluatorThreadForEvoga);
                evaluationFuture.waitForFinished();
                if (commitStep()) return; // stop the evolution process

                // We have finished evaluating parents, updating the fitness vectors
                for (int i = 0; i < popSize; i++) {
                    if (averageIndividualFitnessOverGenerations) {
                        tfitness[evaluators[i]->getGenotypeId()] += evaluators[i]->getFitness();
                        ntfitness[evaluators[i]->getGenotypeId()]++;
                    } else {
                        tfitness[evaluators[i]->getGenotypeId()] = evaluators[i]->getFitness();
                        ntfitness[evaluators[i]->getGenotypeId()] = 1;
                    }
                }
                if (commitStep()) return; // stop the evolution process

                // Calling endGeneration on all evaluator
                for (int i = 0; i < popSize; i++) {
                    evaluators[i]->getExperiment()->endGeneration(gn);
                }
                if (commitStep()) return; // stop the evolution process

            }

            // ========= Selection part ========== //
            // Finally, we look for the worst individuals (parents) and substitute them with best children.
            // What we do is: we order both the parents and the children in descending order, then we take the
            // popSize best individuals. This is not the same as what is done in the sequential version of
            // the algorithm, but should be similar. We overwrite the worst parents with the best children
            QVector<FitnessAndId> parents(popSize);
            for (int i = 0; i < popSize; i++) {
                parents[i].fitness = tfitness[i] / ntfitness[i];
                parents[i].id = i;
            }
            QVector<FitnessAndId> children(popSize);
            for (int i = 0; i < popSize; i++) {
                children[i].fitness = tfitness[popSize + i] / ntfitness[popSize + i];
                children[i].id = popSize + i;
            }
            // Sorting both parents and children. They are sorted in ascending order but we need the best
            // individuals (those with the highest fitness) first
            qSort(parents);
            qSort(children);
            int p = popSize - 1;
            int c = popSize - 1;
            for (int i = 0; i < popSize; i++) {
                if(limitRetention)
                    children[c].fitness *= limitationFactor;
                if (parents[p].fitness > children[c].fitness) {
                    // No need to swap, parents are already in the population vector
                    p--;
                } else {
                    // Swapping with one of the worst parents (we know for sure that p + c = popSize)
                    copyGenes(children[c].id, parents[popSize - 1 - c].id, 0);
                    tfitness[parents[popSize - 1 - c].id] = tfitness[children[c].id];
                    ntfitness[parents[popSize - 1 - c].id] = ntfitness[children[c].id];
                    subsVec[children[c].id-popSize] = parents[popSize - 1 - c].id;
                    c--;
                    currentRetentionRate += 1.0/popSize;
                }
            }
            for (int i=0; i<=c;i++)
                subsVec[children[i].id-popSize] = -1;

            saveBestInd();
            computeFStat2();
            saveFStat();

            if(saveRetStat)
               saveRStat(subsVec);

            emit endGeneration( cgen, fmax, faverage, fmin );
            if (commitStep()) {
                return; // stop the evolution process
            }

            cgen++;
            if (mutation > final_mrate) {
                mutation -= mutationdecay;
            } else {
                mutation = final_mrate;
            }

            limitationFactor += (targetRetentionRate-currentRetentionRate)/10.0;
            if (limitationFactor > 1.0)
                limitationFactor = 1.0;

            //always save in order to be able to resume the evolution process, but keep only the last gen file unless it has to be saved by the param savePopulationEachNGenerations
            saveallg();

            //remove the previous genfile unless it has to be kept because of the savePopulationEachNGenerations param
            if ((savePopulationEachNGenerations == 0) || (gn>1 && ((gn-1) % (savePopulationEachNGenerations) != 0))) {
                //EX: gn 998 = G999S1.gen --- gn 999 = G1000S1.gen --- gn = 1000 = G1001S1.gen --- gn 1001 = G1002S1.gen
                sprintf(filename,"G%dS%d.gen",gn,currentSeed);
                if(!QFile::remove(filename)) {
                    Logger::warning(QString("Error deleting temporary gen file: ") + QString::fromStdString(filename));
                }
            }

            Logger::info(QString("Generation %1 took %2 minutes - Best fitness = %3").arg(gn+1).arg((double)evotimer.elapsed()/60000.0, 0, 'f', 2).arg(fmax));
            if(limitRetention)
                Logger::info(QString(" --- Target retention rate: %1; current rate: %2; fitness limitation factor: %3").arg(targetRetentionRate).arg(currentRetentionRate).arg(limitationFactor));
            fflush(stdout);
        }
        saveallg();

        // Save the best generation fitness statistics
        saveBestFitness();
    }

    // Deleting all evaluators
    for (int i = 0; i < evaluators.size(); i++) {
        delete evaluators[i];
    }
}

/*
 * Main function of the Genetic Algorithm (generational version with truncation selection)
 */
void Evoga::evolveGenerational()
{
    int rp; //replication
    int gn; //generation
    int id; //individuals
    double fit;
    int startGeneration=0;
    char statfile[128];
    char genFile[128];

    // Resizing genome
    genome.resize(popSize);

    for(rp=0;rp<nreplications;rp++) {// replications
        startGeneration = 0;
        setSeed(getStartingSeed()+rp);
        Logger::info(QString("Replication %1 seed: %2").arg(rp+1).arg(getStartingSeed()));
        resetGenerationCounter();
        randomizePop();
        // --- section runnable only if there is an Evonet object
        if ( hasResource( "evonet" ) ) {
            getPheParametersAndMutationsFromEvonet();
        }

        // Set fbest to a very low value
        this->fbest = -99999.0;

        emit startingReplication( rp );

        // Resetting seed in experiments
        exp->newGASeed(getCurrentSeed());

        QTime evotimer;
        evotimer.start();

        //code to recovery a previous evolution: Experimental
        sprintf(statfile,"statS%d.fit", getStartingSeed()+rp);
        //now check if the file exists
        DataChunk statTest(QString("stattest"),Qt::blue,2000,false);
        if (statTest.loadRawData(QString(statfile),0)) {
            startGeneration=statTest.getIndex();
            sprintf(genFile,"G%dS%d.gen",startGeneration,getStartingSeed()+rp);
            Logger::info("Recovering from startGeneration: " + QString::number(startGeneration));
            Logger::info(QString("Loading file: ") + genFile);
            loadallg(startGeneration,genFile);
            emit recoveredInterruptedEvolution( QString(statfile) );
        } //end evolution recovery code

        for(gn=startGeneration;gn<nogenerations;gn++) { // generations
            evotimer.restart();
            Logger::info(" Generation " + QString::number(gn+1));
            exp->initGeneration( gn );
            for(id=0;id<popSize;id++) { //individuals
                exp->setNetParameters(getGenes(id)); // get the free parameters from the genotype
                exp->doAllTrialsForIndividual(id);
                fit = exp->getFitness();
                tfitness[id]=fit;
                if (commitStep()) { // stop evolution
                    return;
                }
            }
            reproduce();

            emit endGeneration( gn, fmax, faverage, fmin );
            exp->endGeneration( gn );

            if(savePopulationEachNGenerations!=0 && gn%savePopulationEachNGenerations == 0)
                saveallg();

            Logger::info(QString("Generation %1 took %2 minutes").arg(gn+1).arg((double)evotimer.elapsed()/60000.0, 0, 'f', 2));
            fflush(stdout);
        }

        saveallg();

        // Save the best generation fitness statistics
        saveBestFitness();
    }
}


// this function generate a random seed for initialize the random number generator
unsigned int generateRandomSeed() {
    // this number is always different amongs processes because this function
    // is on the stack of the process
    unsigned long int stackMem = (unsigned long int)( generateRandomSeed );
    // time on which the process has been started
#ifdef FARSA_WIN
    unsigned long int startTime = GetTickCount();
#else
    unsigned long int startTime = time(NULL);
#endif
    // the seed is generated mixing the values above
    unsigned long int randSeed = 0;

    stackMem=stackMem-startTime;  stackMem=stackMem-randSeed;  stackMem=stackMem^(randSeed >> 13);
    startTime=startTime-randSeed;  startTime=startTime-stackMem;  startTime=startTime^(stackMem << 8);
    randSeed=randSeed-stackMem;  randSeed=randSeed-startTime;  randSeed=randSeed^(startTime >> 13);
    stackMem=stackMem-startTime;  stackMem=stackMem-randSeed;  stackMem=stackMem^(randSeed >> 12);
    startTime=startTime-randSeed;  startTime=startTime-stackMem;  startTime=startTime^(stackMem << 16);
    randSeed=randSeed-stackMem;  randSeed=randSeed-startTime;  randSeed=randSeed^(startTime >> 5);
    stackMem=stackMem-startTime;  stackMem=stackMem-randSeed;  stackMem=stackMem^(randSeed >> 3);
    startTime=startTime-randSeed;  startTime=startTime-stackMem;  startTime=startTime^(stackMem << 10);
    randSeed=randSeed-stackMem;  randSeed=randSeed-startTime;  randSeed=randSeed^(startTime >> 15);

    return randSeed%10000;
}

void Evoga::configure(ConfigurationParameters& params, QString prefix)
{
    genome.clear();
    bestgenome.clear();

    evolutionType = ConfigurationHelper::getString(params, prefix + "evolutionType", "steadyState");
    if ( evolutionType != "steadyState" && evolutionType != "generational" && evolutionType != "xnes" && evolutionType != "specializerSteadyState") {
        Logger::error( "Evoga - evolutionType has been wrongly setted. It can assume only 'steadyState', 'generational', 'xnes' or 'specializerSteadyState' values" );
        evolutionType = "steadyState";
    }
    nogenerations = ConfigurationHelper::getInt(params, prefix + "ngenerations", 100); //number of generations
    nreplications = ConfigurationHelper::getInt(params, prefix + "nreplications", 10); //number of replications
    nreproducing = ConfigurationHelper::getInt(params, prefix + "nreproducing", 20); //number of father
    noffspring = ConfigurationHelper::getInt(params, prefix + "noffspring", 5); //number of sons
    // starting seed - if not specified (or setted as 'auto') it auto generate a seed
    QString checkSeed = ConfigurationHelper::getString(params, prefix + "seed", "auto");
    if ( checkSeed == "auto" ) {
        seed = generateRandomSeed();
    } else {
        seed = checkSeed.toInt();
    }
    Logger::info( QString("Evoga - Random seed setted to ")+QString::number(seed) );
    //seed = ConfigurationHelper::getInt(params, prefix + "seed", 1234); //starting seed
    savebest = ConfigurationHelper::getInt(params, prefix + "savenbest", 1); //saving best n genotype
    elitism = ConfigurationHelper::getBool(params, prefix + "elitism", false); //activate elitism
    numThreads = ConfigurationHelper::getInt(params, prefix + "numThreads", 1); //number of concurrent threads to use
    savePopulationEachNGenerations = ConfigurationHelper::getInt(params, prefix + "savePopulationEachNGenerations", 0);
    averageIndividualFitnessOverGenerations = ConfigurationHelper::getBool(params, prefix + "averageIndividualFitnessOverGenerations", true);
    saveRetStat = ConfigurationHelper::getBool(params, prefix + "saveRetetionStatistics", false);

    limitRetention =  ConfigurationHelper::getBool(params, prefix + "limitRetention", false);
    targetRetentionRate =  ConfigurationHelper::getDouble(params, prefix + "targetRetentionRate", false);
    numModules = ConfigurationHelper::getInt(params, prefix + "numModules", 1);

    crossRate = ConfigurationHelper::getDouble(params, prefix + "crossoverRate", 0.05);;
    modulesMutationRate = ConfigurationHelper::getDouble(params, prefix + "modulesMutationRate", 0.25);
    overwriteRate = ConfigurationHelper::getDouble(params, prefix + "overwriteWrite", 0.05);
    mutateOnlyRelatives = ConfigurationHelper::getBool(params, prefix + "mutateOnlyRelatives", true);
    rankBasedProb = ConfigurationHelper::getDouble(params, prefix + "rankBasedProbability", 0.75);
    selectionType = ConfigurationHelper::getString(params, prefix + "selectionType", "rankBased");
    if ( selectionType != "rankBased" && selectionType != "species") {
        Logger::error( "Evoga - selectionType has been wrongly setted. It can assume only 'rankBased' or 'species' values. 'rankBased' will be used as default value." );
        evolutionType = "rankBased";
    }


    //mutation rate can be written both as int or as double
    mutation = ConfigurationHelper::getDouble(params, prefix + "mutation_rate", mutation); //mutation rate
    if(mutation >= 1) {
        mutation /= 100;
    }
    //speed of decay (only valid for steady-state GA)
    mutationdecay = ConfigurationHelper::getDouble(params, prefix + "mutation_decay", 0.01);
    initial_mutation = ConfigurationHelper::getDouble(params, prefix + "initial_mutation", 0.5);

    mutationRange = ConfigurationHelper::getDouble(params, prefix + "mutationRange", 0.0);
    offspringMutRange = ConfigurationHelper::getDouble(params, prefix + "offspringMutRange", 0.0);
    useGaussian = ConfigurationHelper::getBool(params, prefix + "useGaussian", false);
    gaussianMean = ConfigurationHelper::getDouble(params, prefix + "gaussianMean", 0.0);
    gaussianStdDev = ConfigurationHelper::getDouble(params, prefix + "gaussianStdDev", 1.0);
    offspringStdDev = ConfigurationHelper::getDouble(params, prefix + "offspringStdDev", 1.0);
    variableStdDev = ConfigurationHelper::getBool(params, prefix + "variableStdDev", false);
    numStartInd = ConfigurationHelper::getInt(params, prefix + "numStartInd", 1);
    mutationLearningRate = ConfigurationHelper::getDouble(params, prefix + "mutationLearningRate", 1.0);
    xnesCombination = ConfigurationHelper::getBool(params, prefix + "xnesCombination", false);
    xnesCombinationType = ConfigurationHelper::getInt(params, prefix + "xnesCombinationType", 0);
    if ( xnesCombinationType < 0 || xnesCombinationType > 1) {
        Logger::error( "Evoga - xnesCombinationType has been wrongly setted. It can assume only 0 or 1 values" );
        xnesCombinationType = 0;
    }
    backPropOffspringUtilityRank = ConfigurationHelper::getInt(params, prefix + "backPropOffspringUtilityRank", 255);
    backPropOffspringFitnessIncrease = ConfigurationHelper::getBool(params, prefix + "backPropOffspringFitnessIncrease", false);
    backPropOffspringFitnessIncreasePercentage = ConfigurationHelper::getDouble(params, prefix + "backPropOffspringFitnessIncreasePercentage", 0.1);
    if (backPropOffspringFitnessIncreasePercentage < 0.0 || backPropOffspringFitnessIncreasePercentage > 1.0)
    {
        Logger::error( "Evoga - backPropOffspringFitnessIncreasePercentage has been wrongly setted. It must be within the range [0,1] (it is a percentage!!!)" );
        backPropOffspringFitnessIncreasePercentage = 0.1;
    }
    pheGen = ConfigurationHelper::getInt(params, prefix + "pheGen", 0);
    standardNes = ConfigurationHelper::getBool(params, prefix + "standardNes", false);
    minimization = ConfigurationHelper::getBool(params, prefix + "minimization", false);
    debugPrintInfoToFile = ConfigurationHelper::getBool(params, prefix + "debugPrintInfoToFile", false);
    dissociateInd = ConfigurationHelper::getBool(params, prefix + "dissociateInd", false);
    dissociateCovMatrix = ConfigurationHelper::getBool(params, prefix + "dissociateCovMatrix", false);
    specialUtilityRanking = ConfigurationHelper::getBool(params, prefix + "specialUtilityRanking", false);

    exp = params.getObjectFromGroup<EvoRobotExperiment>(prefix + "Experiment", false);
    exp->setEvoga(this);
    Logger::info( "Created EvoRobotExperiment " + params.getValue(prefix+"Experiment/type") + " from group " + prefix + "Experiment" );

    // Copying the ConfigurationParameters object and the prefix: we will need them to create new instances of the experiment for
    // multithreading evolution
    savedConfigurationParameters = params;
    savedExperimentPrefix = prefix + "Experiment";
}

void Evoga::save(ConfigurationParameters&, QString)
{
    Logger::error("NOT IMPLEMENTED (Evoga::save)");
    abort();
}

void Evoga::describe( QString type ) {
    Descriptor d = addTypeDescription( type, "Implements the genetic algorithm developed by Stefano Nolfi" );
    d.describeEnum( "evolutionType" ).def("steadyState").values( QStringList() << "steadyState" << "generational" << "xnes" << "specializerSteadyState").props( IsMandatory ).help("Specify the type of evolution process to execute");
    d.describeInt( "ngenerations" ).def(100).limits(0,MaxInteger).help("Number of generations");
    d.describeInt( "nreplications" ).def(10).limits(1,MaxInteger).help("The number of which the evolution process will be replicated with a different random initial population");
    d.describeInt( "nreproducing" ).def(20).limits(1,MaxInteger).help("The number of individual allowed to produce offsprings; The size of populazion will be nreproducing x noffspring");
    d.describeInt( "noffspring" ).def(5).limits(1,MaxInteger).help("The number of offsprings generated by an individual; The size of populazion will be nreproducing x noffspring");
    d.describeInt( "seed" ).def(1234).limits(0,MaxInteger).help("The number used to initialize the random number generator; when a new replication will start, this value will be incremented by one to guarantee a truly different initial population for the next replica");
    d.describeInt("savenbest").def(1).limits(1,MaxInteger).help("The number of best genotypes to save each generation");
    d.describeBool("elitism").def(false).help("If use elitism or not");
    d.describeInt("numThreads").def(1).limits(1,MaxInteger).help("The number of thread used to parallelize the evaluation of individuals");
    d.describeInt("savePopulationEachNGenerations").def(0).limits(0,MaxInteger).help("If is zero only the population of the last generation are saved into a file; otherwise it saves the population each N generations done");
    d.describeReal("mutation_rate").def(0.05).limits(0,100).help("The mutation rate", "The rate at which a mutation will occur during a genotype copy; a real value below 1 (i.e. 0.12) is considered as a rate (i.e. 0.12 correspond to 12% of mutation); a value egual or above 1 is considered as a percentage of mutation (i.e. 25 correspond to 25% of mutation, or 0.25 rate of mutation)");
    d.describeReal("mutation_decay").def(0.01).limits(0,1).help("At first generation the mutation rate will be always 0.5, and at each generation done the mutation rate will be decreased by this value until it reachs the mutation_rate value");
    d.describeReal("initial_mutation").def(0.5).limits(0,1).help("The initial value of the mutation rate in case of the steadyState evolution type");
    d.describeSubgroup( "Experiment" ).props(IsMandatory).type("EvoRobotExperiment").help("The object delegated to simulate and to evaluate the fitness of an individual");
    d.describeBool("averageIndividualFitnessOverGenerations").def(true).help("Whether to average the current fitness with the previous one or not");
    d.describeBool("saveRetetionStatistics").def(false).help("Whether to save the retetions statistics or not");
    d.describeReal("mutationRange").def(0.0).help("The mutation range used to generate the individual in case of the xnes evolution type");
    d.describeReal("offspringMutRange").def(0.0).help("The mutation range used to generate the offspring in case of the xnes evolution type");
    d.describeBool("useGaussian").def(false).help("Whether to use a gaussian distribution or not in order to generate offspring in case of the xnes evolution type");
    d.describeReal("gaussianMean").def(0.0).help("The mean of the gaussian distribution used to generate both parent and offspring in case of the xnes evolution type");
    d.describeReal("gaussianStdDev").def(1.0).help("The standard deviation of the gaussian distribution used to generate parent in case of the xnes evolution type");
    d.describeReal("offspringStdDev").def(1.0).help("The standard deviation of the gaussian distribution used to generate offspring in case of the xnes evolution type");
    d.describeBool("variableStdDev").def(false).help("Whether the standard deviation of the gaussian distribution used to generate offspring is variable or constant");
    d.describeInt("numStartInd").def(1).limits(1,MaxInteger).help("The number of individuals to be tested so to find the starting individual in case of the xnes evolution type");
    d.describeReal("mutationLearningRate").def(0.0).help("The learning rate used to generate \"mutated\" offspring in case of the xnes evolution type");
    d.describeBool("xnesCombination").def(false).help("Whether to combine xnes with backpropagation algorithm or not");
    d.describeInt("xnesCombinationType").def(0).limits(0,MaxInteger).help("Type of combination between xnes and backpropagation algorithms");
    d.describeInt("backPropOffspringUtilityRank").def(255).limits(0,MaxInteger).help("The utility rank of the backpropagation algorithm in case of the second type of combination with xnes");
    d.describeBool("backPropOffspringFitnessIncrease").def(false).help("Whether to increase the fitness of the backpropagation offspring in case of the second type of combination with xnes");
    d.describeReal("backPropOffspringFitnessIncreasePercentage").def(0.1).help("Percentage of increase of the backpropagation offspring fitness value in case of the second type of combination with xnes");
    d.describeInt("pheGen").def(0).limits(1,MaxInteger).help("How often the best phenotype must be saved");
    d.describeBool("standardNes").def(false).help("Whether to use standard version of xnes or not");
    d.describeBool("minimization").def(false).help("Whether the fitness function corresponds to a minimization problem");
    d.describeBool("debugPrintInfoToFile").def(false).help("Whether to print debug information to file");
    d.describeBool("dissociateInd").def(false).help("Whether to dissociate back-propagation effect on xNES individual update");
    d.describeBool("dissociateCovMatrix").def(false).help("Whether to dissociate back-propagation effect on xNES covariance matrix update");
    d.describeBool("specialUtilityRanking").def(false).help("Whether to use a special version to compute the utility ranking");
    d.describeInt("numModules").def(1).limits(0,MaxInteger).help("The number of modules that will compose the genotype");

    d.describeReal("crossoverRate").def(0.05).help("The probability that a crossover will occur in an offspring team");
    d.describeReal("modulesMutationRate").def(0.25).help("The probability that a module will be mutated on a team");
    d.describeBool("mutateOnlyRelatives").def(true).help("Wheater a module should be replace by its own offsprings or by any offspring");
    d.describeReal("overwriteRate").def(0.05).help("The learning rate used to generate \"mutated\" offsprings in case of the xnes evolution type");
    d.describeBool("limitRetention").def(false).help("Wheater the retention should be limited or not");
    d.describeReal("targetRetentionRate").def(0.2).help("The rate in which retention should be limited");
    d.describeBool("saveRetetionStatistics").def(false).help("Whether to save the retetions statistics or not");
    d.describeEnum( "selectionType" ).def("rankBased").values( QStringList() << "rankBased" << "species").props( IsMandatory ).help("Specify the type of selection that will be used to define which individuals will survive");
    d.describeReal("rankBasedProbability").def(0.75).help("The probability that the individuals with high fitness will be selected");
}

void Evoga::postConfigureInitialization()
{
    // Sharing resources and declaring which ones we will access
    shareResourcesWith(exp);
    usableResources(QStringList() << "evonet");

    // Allocating memory
    tfitness = new double[MAXINDIVIDUALS];
    statfit = new double*[nogenerations];
    ntfitness = new double[MAXINDIVIDUALS]; //used with super ag implementation
    for(int i = 0; i < nogenerations; i++) {
        statfit[i] = new double[3];
    }
    terror = new double[MAXINDIVIDUALS];

    //allocating memory for the mutation vector
    glen = exp->getGenomeLength();
    mutations = new float[glen];
    for(int mi = 0; mi < glen; mi++) {
        mutations[mi] = Evonet::DEFAULT_VALUE;
    }

    //now we allocate memory for genome and best genome
    cgen = 0;
    popSize = noffspring * nreproducing;

    //dynamic allocation of genome
    genome.setGenomeLength(glen);
    genome.resize(popSize); // An initial set of individuals just to avoid problems...
    for (int i = 0; i < genome.size(); i++) {
        for (int r = 0; r < glen; r++) {
            genome[i][r] = i;//mrand(256);test
        }
    }

    //dynamic allocation of bestgenome
    bestgenome.setGenomeLength(glen);
    bestgenome.resize(nreproducing);

    // Only for XNES algorithm
    if (evolutionType == "xnes")
    {
        float gs = 4.0 + 3.0 * floor(log((float)glen));
        if (xnesCombination && (xnesCombinationType > 0))
        {
            gs += 1.0;
        }
        int genSize = (int)floor(gs);
        genotypes.resize(genSize);
        for (int i = 0; i < genSize; i++)
        {
            genotypes[i] = new float[glen];
            for (int j = 0; j < glen; j++)
            {
                genotypes[i][j] = 0.0;
            }
        }
    }

    if (evolutionType == "specializerSteadyState"){
        nreproducing *= numModules;
        popSize = nreproducing;
    }

    //resetting fitness
    for (int i = 0; i < MAXINDIVIDUALS; i++) {
        tfitness[i] = 0.0;//
    }

    for (int i = 0; i < nogenerations; i++) {
        statfit[i][0] = 0.0; //Average fitness of population
        statfit[i][1] = 0.0; //max fitness
        statfit[i][2] = 0.0; //min fitness
    }

    Logger::info("Evoga Configured - Number of genes: " + QString::number(glen));
}

void Evoga::evolveAllReplicas()
{
    stopEvolution = false;
    if ( evolutionType == "steadyState" ) {
        evolveSteadyState();
    } else if ( evolutionType == "generational" ) {
        evolveGenerational();
    } else if ( evolutionType == "xnes" ) {
        xnes();
    } else if ( evolutionType == "specializerSteadyState" ) {
        if(numModules==1)
            evolveSpecializerSteadyState();
        else
            evolveSpecializerSteadyState2();
    } else {
        Logger::error( QString("Evoga - request to execute a unrecognized evolution type: %1").arg(evolutionType) );
    }
}

void Evoga::stop() {
    stopEvolution = true;
    waitForNextStep.wakeAll();
}

bool Evoga::commitStep() {
    if ( isStepByStep && !stopEvolution ) {
        // will block waiting the command for going ahead
        mutexStepByStep.lock();
        waitForNextStep.wait( &mutexStepByStep );
        mutexStepByStep.unlock();
    }
    return stopEvolution;
}

bool Evoga::isStopped() {
    return stopEvolution;
}

void Evoga::resetStop() {
    stopEvolution = false;
}

void Evoga::enableStepByStep( bool enable ) {
    isStepByStep = enable;
    // if disable the step-by-step it wake any eventually blocked commitStep
    if ( !enable ) {
        waitForNextStep.wakeAll();
    }
}

bool Evoga::isEnabledStepByStep() {
    return isStepByStep;
}

void Evoga::doNextStep() {
    waitForNextStep.wakeAll();
}

EvoRobotExperiment* Evoga::getEvoRobotExperiment()
{
    return exp;
}

QVector<EvoRobotExperiment*> Evoga::getEvoRobotExperimentPool()
{
    QVector<EvoRobotExperiment*> v;
    v.append(exp);
    return v;
}

unsigned int Evoga::getCurrentGeneration()
{
    return cgen;
}

unsigned int Evoga::getStartingSeed()
{
    return seed;
}

unsigned int Evoga::getCurrentSeed()
{
    return currentSeed;
}

unsigned int Evoga::getNumReplications()
{
    return nreplications;
}

unsigned int Evoga::getNumOfGenerations() {
    return nogenerations;
}

double Evoga::getCurrentMutationRate() {
    return mutation;
}

void Evoga::setCurrentMutationRate( double mutation_rate ) {
    mutation = mutation_rate;
}

unsigned int Evoga::loadGenotypes(QString filename)
{
    if (evolutionType == "xnes")
    {
        return xnesLoadAllGen(-1, filename.toLatin1().data());
    }
    else
    {
        return loadallg(-1, filename.toLatin1().data());
    }
}

unsigned int Evoga::numLoadedGenotypes() const
{
    return loadedIndividuals;
}

int* Evoga::getGenesForIndividual(unsigned int id)
{
    return getGenes(id);
}

float* Evoga::xnesGetGenesForIndividual(unsigned int id)
{
    return xnesGetGenes(id);
}

QString Evoga::statisticsFilename(unsigned int seed)
{
    return "statS" + QString::number(seed) + QString(".fit");
}

QString Evoga::bestsFilename(unsigned int seed)
{
    return "B0S" + QString::number(seed) + QString(".gen");
}

QString Evoga::bestsFilename()
{
    return "B0S*.gen";
}

QString Evoga::generationFilename(unsigned int generation, unsigned int seed)
{
    return "G" + QString::number(generation) + "S" + QString::number(seed) + QString(".gen");
}

QString Evoga::generationFilename()
{
    return "G*S*.gen";
}

void Evoga::doNotUseMultipleThreads()
{
    numThreads = 1;
}

QString Evoga::retentionsFilename(unsigned int seed)
{
    return "retentionS" + QString::number(seed) + QString(".stat");
}


//SPECIALIZER EVOLUTION METHODS
/*
 * Main function of the Genetic Algorithm (Specializer Steady State Version)
 * --- at the moment works only for positive fitness values
 */

void Evoga::evolveSpecializerSteadyState()
{   //SEQUENTIAL VERSION
    int rp; //replication
    int gn; //generation
    int id; //individuals
    double initFit;
    double minfit=9999;
    int    minid=-1;
    double maxFit;
    int    indMax;
    float  final_mrate;
    int startGeneration=0;
    char statfile[128];
    char genFile[128];
    char filename[64];
    FILE * pFile;
    QVector< int* > actualTeam;
    QVector< int* > allIndividuals;
    QVector< double > subsFitness;
    int raffle,nEval;
    QVector< int > subsVec;
    double currentRetentionRate;

    subsVec.resize(popSize);
    actualTeam.resize(popSize);
    allIndividuals.resize(popSize*2);
    subsFitness.resize(popSize);

    final_mrate = mutation;
    Logger::info("EVOLUTION: specializer steady state (sequential)");
    Logger::info("Number of replications: " + QString::number(nreplications));

    // Creating evaluator objects in case of a multithread simulation. Also setting the actual number of threads used
    QVector<EvaluatorThreadForEvoga*> evaluators(popSize, NULL);
    if (numThreads > 1) {
        for (int i = 0; i < evaluators.size(); i++) {
            EvoRobotExperiment* newExp = savedConfigurationParameters.getObjectFromGroup<EvoRobotExperiment>(savedExperimentPrefix, false);
            newExp->setEvoga(this);
            evaluators[i] = new EvaluatorThreadForEvoga(this, newExp);
        }
        QThreadPool::globalInstance()->setMaxThreadCount(numThreads);
    }
    for(rp=0;rp<nreplications;rp++) {   // replications
        startGeneration=0;
        limitationFactor = 1.0;
        mutation=initial_mutation; // initially mutation (default 50%)
        //setSeed(getSeed());
        setSeed(getStartingSeed()+rp);
        Logger::info(QString("Replication %1, seed: %2").arg(rp+1).arg(getStartingSeed()+rp));
        resetGenerationCounter();
        randomizePop();
        // --- section runnable only if there is an Evonet object
        if ( hasResource( "evonet" ) ) {
            getPheParametersAndMutationsFromEvonet();
        }

        emit startingReplication( rp );
        QTime evotimer;
        evotimer.start();

        //code to recovery a previous evolution
        sprintf(statfile,"statS%d.fit", getStartingSeed()+rp);
        //now check if the file exists
        DataChunk statTest(QString("stattest"),Qt::blue,2000,false);
        if (statTest.loadRawData(QString(statfile),0)) {
            startGeneration=statTest.getIndex();
            sprintf(genFile,"G%dS%d.gen",startGeneration,getStartingSeed()+rp);
            Logger::info("Recovering from startGeneration: " + QString::number(startGeneration));
            Logger::info(QString("Loading file: ") + genFile);
            loadallg(startGeneration,genFile);
            cgen=startGeneration;
            mutation=mutation-startGeneration*mutationdecay;
            if (mutation<final_mrate) mutation=final_mrate;
            emit recoveredInterruptedEvolution( QString(statfile) );
        }//end evolution recovery code

        // Resizing genome
        loadedIndividuals = popSize;
        genome.resize(popSize * 2);
        for(gn=startGeneration;gn<nogenerations;gn++) { // generations
            evotimer.restart();
            Logger::info(QString("Generation %1").arg(gn+1));
            for (int i = 0; i < popSize; i++) { //produce the offsprings
                copyGenes(i, popSize + i, 1); //generate a variation by duplicating and mutating
                //evaluators[i]->setGenotype(popSize + i);
                allIndividuals[i]=actualTeam[i]=getGenesForIndividual(i);
                allIndividuals[i+popSize]=getGenesForIndividual(i+popSize);
            }
            //test the actual team of individuals - this should be done first of all so it is done in a monothread mode
            exp->initGeneration(gn);
            //copy here the individuals to be tested
            for(int i=0;i<popSize;i++)
                exp->m_weightIndividual[i] = actualTeam[i];
            if ( commitStep() ) { return; }
            initFit=0.0;
            exp->doAllTrialsForIndividual(-1); //-1 means that the team will be tested without any substitution
            initFit = exp->getFitness();
            if (isStopped()) { // stop evolution
                return;
            }
            exp->endGeneration(gn);
            if ( commitStep() ) { return; }
            //end evaluation the actual team

            for(int offspring=0;offspring<popSize;offspring++){
                if(numThreads>1){ //MULTITHREAD
                    // Creating evaluator objects and setting the actual number of threads used
                    for (int i = 0; i < evaluators.size(); i++) {
                        evaluators[i]->getExperiment()->initGeneration(gn);
                        evaluators[i]->getExperiment()->m_weightIndividual[i] = allIndividuals[offspring+popSize]; //change the individual i by the offspring being tested
                   }
                   QThreadPool::globalInstance()->setMaxThreadCount(numThreads);
                   if (commitStep()) return; // stop the evolution process

                   // Now starting parallel evaluation of parents and wating for it to finish
                   QFuture<void> evaluationFuture;
                   evaluationFuture = QtConcurrent::map(evaluators, runEvaluatorThreadForEvoga);
                   evaluationFuture.waitForFinished();
                   if (commitStep()) return; // stop the evolution process
                   //get the fitness value for each substitution
                   for (int i = 0; i < popSize; i++) {
                      subsFitness[i] = evaluators[i]->getFitness();
                      evaluators[i]->getExperiment()->m_weightIndividual[i] = actualTeam[i];
                   }

                }else{ //MONOTHREAD
                     for (int i = 0; i < popSize;i++) {
                        exp->m_weightIndividual[i] = allIndividuals[offspring+popSize]; //change the individual i by the offspring being tested
                        exp->doAllTrialsForIndividual(i); //the individual passed as parameter (i) is the substituted individual
                        subsFitness[i] = exp->getFitness();
                        exp->m_weightIndividual[i] = actualTeam[i]; //replace the individual substituted
                     }
                }

                //HERE STARTS THE SELECTION PART - each should be done one offspring at a time
                currentRetentionRate = 0.0;
                maxFit = initFit;
                indMax = -1;
                for(int i=0;i<popSize;i++){ //search which substitution produces the greatest improvement
                    if(limitRetention)
                        subsFitness[i] *= limitationFactor;
                    if(subsFitness[i]>maxFit){
                       maxFit = subsFitness[i];
                       indMax = i;
                    }
                    Logger::info(QString("Changing the policy %1 by the policy %2 produces a fitness of: %3 (original is %4)").arg(i).arg(offspring+popSize).arg(subsFitness[i]).arg(initFit));
                }
                if(indMax!=-1){ //if some substitution produces improvement, replace it with the offspring
                    copyGenes(offspring+popSize, indMax, 0);
                    actualTeam[indMax] = allIndividuals[offspring+popSize];
                    Logger::info(QString("Substituting the policy %1 by the policy %2, and increasing the fitness by %3").arg(indMax).arg(offspring+popSize).arg(maxFit-initFit));
                    initFit = maxFit;
                    subsVec[offspring] = indMax;
                    currentRetentionRate += 1.0/popSize;
                }else{
                    Logger::info(QString("Policy %1 won't be used because it presents no performance improvement.").arg(offspring+popSize));
                    subsVec[offspring] = -1;
                }
            }

            //statistics to keep compatibility -- set all values with the maxFitness of this generation
            statfit[cgen][0] = statfit[cgen][1] = statfit[cgen][2] = faverage = fmax = fmin = maxFit;
            saveFStat();

            if(saveRetStat)
               saveRStat(subsVec);

            emit endGeneration( cgen, fmax, faverage, fmin );

            if (commitStep()) {
                return; // stop the evolution process
            }

            cgen++;
            if (mutation > final_mrate) {
                mutation -= mutationdecay;
            } else {
                mutation = final_mrate;
            }

            limitationFactor += (targetRetentionRate-currentRetentionRate)/10.0;
            if (limitationFactor > 1.0)
                limitationFactor = 1.0;

            //always save in order to be able to resume the evolution process, but keep only the last gen file unless it has to be saved by the param savePopulationEachNGenerations
            saveallg();

            //remove the previous genfile unless it has to be kept because of the savePopulationEachNGenerations param
            if ((savePopulationEachNGenerations == 0) || (gn>1 && ((gn-1) % (savePopulationEachNGenerations) != 0))) {
                //EX: gn 998 = G999S1.gen --- gn 999 = G1000S1.gen --- gn = 1000 = G1001S1.gen --- gn 1001 = G1002S1.gen
                sprintf(filename,"G%dS%d.gen",gn,currentSeed);
                if(!QFile::remove(filename)) {
                    Logger::warning(QString("Error deleting temporary gen file: ") + QString::fromStdString(filename));
                }
            }

            Logger::info(QString("Generation %1 took %2 minutes").arg(gn+1).arg((double)evotimer.elapsed()/60000.0, 0, 'f', 2));
            if(limitRetention)
                Logger::info(QString(" --- Target retention rate: %1; current rate: %2; fitness limitation factor: %3").arg(targetRetentionRate).arg(currentRetentionRate).arg(limitationFactor));
            fflush(stdout);
            //---------------------------------
        }
        saveallg();
    }
    if(numThreads>1){
        // Deleting all evaluators
        for (int i = 0; i < evaluators.size(); i++) {
            delete evaluators[i];
        }
    }
}

//SPECIALIZER EVOLUTION METHODS
/*
 * Main function of the Genetic Algorithm (Specializer Steady State Version)
 *
 */

void Evoga::evolveSpecializerSteadyState2()
{   //TEAM VERSION
    int rp; //replication
    int gn; //generation

    float  final_mrate;
    int startGeneration=0;
    char statfile[128];
    char genFile[128];
    char filename[64];
    double currentRetentionRate;
    int numTeams = popSize/numModules;
    int totalTeams = numTeams + numTeams*noffspring;
    QVector< int* > allIndividuals(popSize*2);
    QVector< QVector<int> > composedGen(totalTeams);
    QVector< double > parentFitVec(numTeams);
    QVector<FitnessAndId> allTeamsFitness(totalTeams);
    QVector<FitnessAndId> selectedTeamFitness(numTeams);
    QVector<FitnessAndId> parentFitness(numTeams);
    QVector<FitnessAndId> childrenFitness(numTeams*noffspring);
    QVector< QVector<int> > savedTeams(numTeams);
    QVector<int*> indivSelected(popSize);
    QVector<int> translatedIndiv(popSize*2),subsVec(numTeams);
    QVector< QVector<int> > translatedComposedTeams(totalTeams);
    int winner,overwrited,raffle,idOffspring,numSelected;
    bool flagModificated;

    for(int i=0;i<composedGen.size();i++)
        composedGen[i].resize(numModules);

    final_mrate = mutation;
    Logger::info("EVOLUTION: specializer steady state (team representation(modules))");
    Logger::info("Number of replications: " + QString::number(nreplications));

    // Creating evaluator objects in case of a multithread simulation. Also setting the actual number of threads used
    QVector<EvaluatorThreadForEvoga*> evaluators(numTeams, NULL);
    if (numThreads > 1) {
        for (int i = 0; i < evaluators.size(); i++) {
            EvoRobotExperiment* newExp = savedConfigurationParameters.getObjectFromGroup<EvoRobotExperiment>(savedExperimentPrefix, false);
            newExp->setEvoga(this);
            evaluators[i] = new EvaluatorThreadForEvoga(this, newExp);
        }
        QThreadPool::globalInstance()->setMaxThreadCount(numThreads);
    }
    for(rp=0;rp<nreplications;rp++) {   // replications
        startGeneration=0;

        //1-INITIALIZE PARENT RANDOM TEAMS
        for(int i=0;i<numTeams;i++){ //initializing random teams
            for(int j=0;j<numModules;j++)
                composedGen[i][j] = rand()%(popSize);
        }


        limitationFactor = 1.0;
        mutation=initial_mutation; // initially mutation (default 50%)
        //setSeed(getSeed());
        setSeed(getStartingSeed()+rp);
        Logger::info(QString("Replication %1, seed: %2").arg(rp+1).arg(getStartingSeed()+rp));
        resetGenerationCounter();
        randomizePop();

        // --- section runnable only if there is an Evonet object
        if ( hasResource( "evonet" ) ) {
            getPheParametersAndMutationsFromEvonet();
        }

        emit startingReplication( rp );
        QTime evotimer;
        evotimer.start();

        //1.1-LOAD SAVED TEAMS IF THERE IS A PREVIOUS EVOLUTION
        //code to recovery a previous evolution
        sprintf(statfile,"statS%d.fit", getStartingSeed()+rp);
        //now check if the file exists
        DataChunk statTest(QString("stattest"),Qt::blue,2000,false);
        if (statTest.loadRawData(QString(statfile),0)) {
            startGeneration=statTest.getIndex();
            sprintf(genFile,"G%dS%d.gen",startGeneration,getStartingSeed()+rp);
            Logger::info("Recovering from startGeneration: " + QString::number(startGeneration));
            Logger::info(QString("Loading file: ") + genFile);
            loadallg(-1,genFile);
            sprintf(genFile,"G%dS%d.composed.gen",startGeneration,getStartingSeed()+rp);
            loadallTeams(-1,genFile,composedGen);
            cgen=startGeneration;
            mutation=mutation-startGeneration*mutationdecay;
            if (mutation<final_mrate) mutation=final_mrate;
            emit recoveredInterruptedEvolution( QString(statfile) );

        }//end evolution recovery code

        // Resizing genome
        loadedIndividuals = popSize;
        genome.resize(popSize * 2);

        for(gn=startGeneration;gn<nogenerations;gn++) { // generations
            evotimer.restart();
            Logger::info(QString("Generation %1").arg(gn+1));
            for (int i = 0; i < popSize; i++) { //produce the individuals offsprings
                copyGenes(i, popSize + i, 1); //generate a variation by duplicating and mutating
                //evaluators[i]->setGenotype(popSize + i);
                allIndividuals[i]=getGenesForIndividual(i);
                allIndividuals[i+popSize]=getGenesForIndividual(i+popSize);
            }

            //2-EVALUATE PARENTS
            if(numThreads>1){ //MULTITHREAD
               // Creating evaluator objects and setting the actual number of threads used
                for (int i = 0; i < evaluators.size(); i++) {
                    evaluators[i]->getExperiment()->initGeneration(gn);
                    evaluators[i]->getExperiment()->m_weightIndividual.resize(numModules);
                    for(int j=0;j<numModules;j++) //setting the weights of the team i to the evaluator
                        evaluators[i]->getExperiment()->m_weightIndividual[j] = allIndividuals[composedGen[i][j]];
                }
               QThreadPool::globalInstance()->setMaxThreadCount(numThreads);

               if (commitStep()) return; // stop the evolution process

               // Now starting parallel evaluation of parents and wating for it to finish
               QFuture<void> evaluationFuture;
               evaluationFuture = QtConcurrent::map(evaluators, runEvaluatorThreadForEvoga);
               evaluationFuture.waitForFinished();
               if (commitStep()) return; // stop the evolution process
               //get fitness
               for (int i = 0; i < numTeams; i++) {
                  parentFitVec[i] = evaluators[i]->getFitness();
                  parentFitness[i].id = allTeamsFitness[i].id = i;
                  parentFitness[i].fitness = allTeamsFitness[i].fitness = evaluators[i]->getFitness();
               }
            }else{ //MONOTHREAD
                exp->initGeneration(gn);
                for (int i = 0; i < numTeams;i++) { //test parents
                    exp->m_weightIndividual.resize(numModules);
                    for(int j=0;j<numModules;j++)
                        exp->m_weightIndividual[j] = allIndividuals[composedGen[i][j]];

                    exp->doAllTrialsForIndividual(i);
                    //get fitness
                    parentFitVec[i] = exp->getFitness();
                    parentFitness[i].id = allTeamsFitness[i].id = i;
                    parentFitness[i].fitness = allTeamsFitness[i].fitness = exp->getFitness();
                 }
            }
             //3-CROSSOVER, MUTATE AND DUPLICATE FOR GENERATING THE OFFSPRINGS
            for(int i=0;i<numTeams;i++){ //for each offspring to be generated

                //GENERATE ONE OFFSPRING
                for(int offspring=0;offspring<noffspring;offspring++){
                    idOffspring = i*noffspring+numTeams+offspring;
                    //copy all modules from the parent to the offspring
                    for(int j=0;j<numModules;j++)
                        composedGen[idOffspring][j] = composedGen[i][j];

                    flagModificated = false;

                    //3.1-CROSSOVER
                    raffle = rand()%100;
                    if(raffle/100.0<crossRate){
                        //select an individual according to roulette wheel
                        winner = rouletteWheel(parentFitVec);
                        //select the crossover point
                        raffle = rand()%(numModules-2);
                        int initialPoint,finalPoint;
                        if(rand()%2==0){ //define the order of the crossover: 1st parent then winner
                            //copy the individuals from the winner individual
                            initialPoint = (raffle+1);
                            finalPoint = numModules;
                        }else{//1st winner then parent
                            //copy the individuals from the original parent
                            initialPoint = 0;
                            finalPoint = numModules;
                        }
                        for(int j=initialPoint;j<finalPoint;j++)
                           composedGen[idOffspring][j] = composedGen[winner][j];
                        flagModificated = true;
                    }

                    //3.2-MUTATION
                    for(int j=0;j<numModules;j++){
                        raffle = rand()%100;
                        if(raffle/100.0<modulesMutationRate){ //check if this individual will be mutated
                            if(mutateOnlyRelatives) //if yes, check if it should be substituted by its offspring, or any offspring
                               composedGen[idOffspring][j] += popSize;
                            else
                               composedGen[idOffspring][j] = rand()%popSize+popSize;
                            flagModificated = true;
                        }
                    }

                    //3.3-DUPLICATION
                    raffle = rand()%100;
                    if(raffle/100.0<overwriteRate){
                        winner = rand()%numModules;
                        for(overwrited=rand()%numModules;overwrited==winner;overwrited=rand()%numModules);
                        composedGen[idOffspring][overwrited] = composedGen[idOffspring][winner];
                        flagModificated = true;
                    }

                    if(!flagModificated) //if the offspring is equal to the parent, repeat the process for this offspring
                        offspring--;
                }
            }

             //4-EVALUATE CHILDREN
             if(numThreads>1){ //MULTITHREAD

                 for(int offspring=0;offspring<noffspring;offspring++){
                     // Creating evaluator objects and setting the actual number of threads used
                      for (int i = 0; i < evaluators.size(); i++) {
                          idOffspring = i*noffspring+numTeams+offspring;
                          evaluators[i]->getExperiment()->m_weightIndividual.resize(composedGen[idOffspring].size());//for the case where teams can have different sizes
                          for(int j=0;j<composedGen[idOffspring].size();j++) //setting the weights of the team i to the evaluator
                              evaluators[i]->getExperiment()->m_weightIndividual[j] = allIndividuals[composedGen[idOffspring][j]];
                      }
                     QThreadPool::globalInstance()->setMaxThreadCount(numThreads);

                     if (commitStep()) return; // stop the evolution process

                     // Now starting parallel evaluation of children and wating for it to finish
                     QFuture<void> evaluationFuture;
                     evaluationFuture = QtConcurrent::map(evaluators, runEvaluatorThreadForEvoga);
                     evaluationFuture.waitForFinished();
                     if (commitStep()) return; // stop the evolution process
                     //get fitness
                     for (int i = 0; i < evaluators.size(); i++) {
                        idOffspring = i*noffspring+numTeams+offspring;
                        childrenFitness[i].id = allTeamsFitness[idOffspring].id = idOffspring;
                        childrenFitness[i].fitness = allTeamsFitness[idOffspring].fitness = evaluators[i]->getFitness();
                        if(limitRetention)
                            childrenFitness[i].fitness *= limitationFactor;
                            allTeamsFitness[idOffspring].fitness *= limitationFactor;
                     }
                 }

                 // Calling endGeneration on all evaluator
                 for (int i = 0; i < evaluators.size(); i++) {
                     evaluators[i]->getExperiment()->endGeneration(gn);
                 }
                 if (commitStep()) return; // stop the evolution process

             }else{ //MONOTHREAD
                 for (int i = 0; i < totalTeams - numTeams;i++) { //test children - testing all but the parents
                    exp->m_weightIndividual.resize(numModules);
                    for(int j=0;j<numModules;j++)
                        exp->m_weightIndividual[j] = allIndividuals[composedGen[numTeams+i][j]];
                    exp->doAllTrialsForIndividual(i+numTeams);
                    //get fitness
                    childrenFitness[i].id = allTeamsFitness[i+numTeams].id = i+numTeams;
                    childrenFitness[i].fitness = allTeamsFitness[i+numTeams].fitness = exp->getFitness();
                    if(limitRetention){
                        childrenFitness[i].fitness *= limitationFactor;
                        allTeamsFitness[i+numTeams].fitness *= limitationFactor;
                    }
                 }
                 exp->endGeneration(gn);
             }

            //5-SELECTION

            translatedIndiv.fill(-1);
            currentRetentionRate = 0.0;
            if(selectionType=="rankBased"){ //GET INDIVIDUALS WITH THE BEST FITNESS
                qSort(allTeamsFitness); //sort in ascending order
                numSelected = 0;
                for(int i=totalTeams-1;i>0;i--){ //now start with the teams with higher fitness - and proceed until numTeams were selected

                    if((rand()%100)/100.0<rankBasedProb){//select the teams with greatest fitness with a probability defined by the user
                        selectedTeamFitness[numSelected] = allTeamsFitness[i];
                        parentFitVec[numSelected] = allTeamsFitness[i].fitness;
                        numSelected++;
                    }

                    if(numSelected==numTeams)
                        break;
                }

                for(int i=0;i<totalTeams;i++)
                    Logger::info(QString("i=%3 id=%1; fitness=%2").arg(allTeamsFitness[i].id).arg(allTeamsFitness[i].fitness).arg(i));

            }else if(selectionType=="species"){ //substitute the parent by its offspring with the best fitness
                for(int i=0;i<numTeams;i++){
                    Logger::info(QString("i=%3 id=%1; fitness=%2").arg(allTeamsFitness[i].id).arg(allTeamsFitness[i].fitness).arg(i));
                    selectedTeamFitness[i] = allTeamsFitness[i]; //the parent is initially selected
                    for(int offspring=0;offspring<noffspring;offspring++){ //check each offspring team and select it if its fitness is greater than actual selected one
                        idOffspring = i*noffspring+numTeams+offspring;
                        Logger::info(QString("Offsprings =======>i=%3 id=%1; fitness=%2").arg(allTeamsFitness[idOffspring].id).arg(allTeamsFitness[idOffspring].fitness).arg(idOffspring));
                        if(selectedTeamFitness[i].fitness<allTeamsFitness[idOffspring].fitness){ //if some of its offspring has a better fitness, select it
                            selectedTeamFitness[i] = allTeamsFitness[idOffspring];
                            parentFitVec[i] = allTeamsFitness[idOffspring].fitness;
                        }
                    }
                }
            }

            //translating and completing selected individuals
            numSelected = 0;
            qSort(allTeamsFitness); //sort in ascending order
            FitnessAndId translateTeam;
            int ind = totalTeams-1;
            for(int i=0;i<totalTeams;i++){
                translatedComposedTeams[i].resize(numModules);
                if(i<numTeams){ //pick individuals from the selected teams vector
                    translateTeam = selectedTeamFitness[i];

                    if(selectedTeamFitness[i].id>=numTeams) //if the id is greater than the number of teams, it means that a parent was replaced
                        currentRetentionRate += 1.0/numTeams;

                }else{ //if there were no sufficient individuals in the selected team vector, get individuals from the teams with higher fitness that weren't selected
                    for(int i=0;i<numTeams;i++){
                        if(allTeamsFitness[ind].id==selectedTeamFitness[i].id){ //if this team is already selected pick the next one and restart the check to see if the next is already selected or not
                            ind--;
                            i=-1;
                        }
                    }
                    translateTeam = allTeamsFitness[ind];
                    ind--;
                }

                for(int j=0;j<numModules;j++){ //since the individuals that participates in the teams with the better fitness value will be selected, it will be necessary to change(translate) their id
                    if(translatedIndiv[composedGen[translateTeam.id][j]]==-1){
                       translatedIndiv[composedGen[translateTeam.id][j]] = numSelected;
                       copyGenes(composedGen[translateTeam.id][j], numSelected, 0);
                       indivSelected[numSelected] = allIndividuals[composedGen[translateTeam.id][j]];
                       numSelected++;
                    }
                    translatedComposedTeams[i][j] = translatedIndiv[composedGen[translateTeam.id][j]];
                    if(numSelected==popSize)
                        break;
                }
                if(numSelected==popSize)
                    break;
            }

           //computing statistics and saving teams for the next generation
           faverage = 0;
           fmax = fmin = parentFitVec[0];
            for(int i=0;i<numTeams;i++){
                composedGen[i] = savedTeams[i] = translatedComposedTeams[i];
                subsVec[i] = selectedTeamFitness[i].id;

                faverage += parentFitVec[i];
                if(parentFitVec[i]>fmax)
                    fmax = parentFitVec[i];
                else if(parentFitVec[i]<fmin)
                    fmin = parentFitVec[i];
            }
            faverage /= numTeams;
            //END SELECTION

            //statistics to keep compatibility -- set all values with the maxFitness of this generation
            statfit[cgen][0]=fmax;
            statfit[cgen][1]=faverage;
            statfit[cgen][2]=fmin;
            saveFStat();

            if(saveRetStat)
               saveRStat(subsVec);

            emit endGeneration( cgen, fmax, faverage, fmin );

            if (commitStep()) {
                return; // stop the evolution process
            }

            cgen++;
            if (mutation > final_mrate) {
                mutation -= mutationdecay;
            } else {
                mutation = final_mrate;
            }

            limitationFactor += (targetRetentionRate-currentRetentionRate)/10.0;
            if (limitationFactor > 1.0)
                limitationFactor = 1.0;


            //always save in order to be able to resume the evolution process, but keep only the last gen file unless it has to be saved by the param savePopulationEachNGenerations
            saveallg();
            saveallgComposed(savedTeams);
            saveBestTeam(savedTeams,parentFitVec);


            //remove the previous genfile unless it has to be kept because of the savePopulationEachNGenerations param
            if ((savePopulationEachNGenerations == 0) || (gn>1 && ((gn-1) % (savePopulationEachNGenerations) != 0))) {
                //EX: gn 998 = G999S1.gen --- gn 999 = G1000S1.gen --- gn = 1000 = G1001S1.gen --- gn 1001 = G1002S1.gen
                sprintf(filename,"G%dS%d.gen",gn,currentSeed);
                if(!QFile::remove(filename)) {
                    Logger::warning(QString("Error deleting temporary gen file: ") + QString::fromStdString(filename));
                }

                sprintf(filename,"G%dS%d.composed.gen",gn,currentSeed);
                if(!QFile::remove(filename)) {
                    Logger::warning(QString("Error deleting temporary gen file: ") + QString::fromStdString(filename));
                }

                sprintf(filename,"B%dS%d.G%d.gen",0,this->currentSeed,gn);
                if(!QFile::remove(filename)) {
                    Logger::warning(QString("Error deleting temporary gen file: ") + QString::fromStdString(filename));
                }
            }

            Logger::info(QString("Generation %1 took %2 minutes - best individual fitness %3").arg(gn+1).arg((double)evotimer.elapsed()/60000.0, 0, 'f', 2).arg(fmax));
            if(limitRetention)
                Logger::info(QString(" --- Target retention rate: %1; current rate: %2; fitness limitation factor: %3").arg(targetRetentionRate).arg(currentRetentionRate).arg(limitationFactor));
            fflush(stdout);
            //---------------------------------
        }
        if(savedTeams[0].size()!=0){
            saveallg();
            saveallgComposed(savedTeams);
            saveBestTeam(savedTeams,parentFitVec);
        }
    }
    if(numThreads>1){
        // Deleting all evaluators
        for (int i = 0; i < evaluators.size(); i++) {
            delete evaluators[i];
        }
    }
}


//END SPECIALIZER EVOLUTION METHODS

int Evoga::rouletteWheel(QVector<double> candidates){
    long double raffle, sum = 0.0;
    int i;
    for(i=0;i<candidates.size();sum += candidates[i]+1,i++);
    raffle = (rand()%((int) sum*1000))/1000.0; //generating a random double with 3 decimal places
    for(i=0;raffle>0;raffle-=(candidates[i]+1),i++);
    return i-1;
}

} // end namespace farsa

// All the suff below is to restore the warning state on Windows
#if defined(_MSC_VER)
    #pragma warning(pop)
#endif