/*************************************************************************
* Project: Library of Evolutionary Algoriths
*************************************************************************
* Author: Changhe Li & Ming Yang
* Email: changhe.lw@google.com Or yangming0702@gmail.com
* Language: C++
*************************************************************************
*  This file is part of EAlib. This library 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, or (at your option) any later version.
*************************************************************************/
// Created: 21 July 2011
// Last modified:
#include "PerformSingleObj.h"
#include "../ToolnDef/Global.h"
#include "../Problems/DOPs/DynamicContinuous.h"
#include "../Problems/FunctionOpt/BenchmarkFunction.h"
#include "../Algorithms/Algorithm.h"
#include "../Problems/DOPs/MovingPeak.h"
PerformSingleObj * PerformSingleObj::msp_perf=0;
PerformSingleObj::PerformSingleObj(bool rGOpt):m_gOptIdx(0),m_numEvals2Suc(0),m_numEvals2Converge(0),m_progrOutputFlag(true),m_convgMode(PROGR_MEAN){
    //ctor
    m_numRecords=Global::g_tEvals/Global::g_sampleFre;
    mpp_data=new double*[Global::g_numRuns];
    for(int i=0;i<Global::g_numRuns;i++) mpp_data[i]=new double[m_numRecords];

    switch(Global::gp_problem->getId()){
        case Sphere: case Sphere_Noisy: case S_Sphere: case R_Sphere: case RS_Sphere: case S_Sphere_CEC05:
        case Rastrigin: case Rastrigin_Noisy: case S_Rastrigin: case R_Rastrigin: case RS_Rastrigin:
        case S_Rastrigin_CEC05: case RS_Rastrigin_CEC05: case Weierstrass: case RS_Weierstrass: case R_Weierstrass:
        case RS_Weierstrass_CEC05: case Griewank: case R_Griewank: case RS_Griewank: case Ackley: case Ackley_Noisy:
        case S_Ackley: case R_Ackley: case RS_Ackley: case RS_Ackley_Bound_CEC05: case Step: case Quartic_Noisy:
        case Scaffer_F6: case Rosenbrock: case S_Rosenbrock: case S_Rosenbrock_CEC05: case Schwefel_2_22:
        case  Schwefel_2_22_Noisy: case S_Schwefel_2_22: case R_Schwefel_2_22: case RS_Schwefel_2_22: case Schwefel:
        case Schwefel_Noisy: case S_Schwefel: case R_Schwefel: case RS_Schwefel: case Schwefel_1_2: case Schwefel_1_2_Noisy:
        case S_Schwefel_1_2: case S_Schwefel_1_2_Noisy: case R_Schwefel_1_2: case RS_Schwefel_1_2: case  S_Schwefel_1_2_CEC05:
        case  S_Schwefel_1_2_Noisy_CEC05: case Schwefel_2_21: case Penalized_1: case Penalized_2: case Noncont_Rastrigin:
        case RS_Elliptic_CEC05: case Elliptic: case Com: case R_Com: case Com_CEC05: case H_Com_CEC05: case H_Com_Noisy_CEC05:
        case RH_Com_CEC05: case RH_Com_NarrowBasin_CEC05: case RH_Com_Bound_CEC05: case M_global1: case M_global2:
        case M_global3: case M_global4: case M_global5: case S_Sphere_CEC08: case Schwefel_2_21_CEC08: case S_Rosenbrock_CEC08:
        case S_Rastrigin_CEC08: case S_Griewank_CEC08: case S_Ackley_CEC08: case RW_Gear_Train: case RW_ParEst_FMSoundWaves:
            m_numGopts=1;
            m_recordsPerChange=Global::g_tEvals/Global::g_sampleFre;
            break;
        ///Dynamic Problems
        case CompositionDBG_DOP: case RotationDBG_DOP: case MovingPeak_DOP: case DF1_DOP: case Binary_DOP:{
             m_numGopts=Global::g_tEvals/Global::g_changeFre;
            m_recordsPerChange=Global::g_changeFre/Global::g_sampleFre;
            break;
        }
    }
    mpp_convergeTime=new int*[Global::g_numRuns];
    for(int i=0;i<Global::g_numRuns;i++) mpp_convergeTime[i]=new int [m_numGopts];

    if(rGOpt){
        mpp_gOpt=new double *[Global::g_numRuns];
        for(int i=0;i<Global::g_numRuns;i++) mpp_gOpt[i]=new double [m_numGopts];

    }else{
        mpp_gOpt=0;
    }


}
PerformSingleObj::PerformSingleObj():mpp_data(0),mpp_gOpt(0),mpp_convergeTime(0){
    //ctor

}
PerformSingleObj::~PerformSingleObj(){
    //dtor
    if(mpp_data){
        for(int i=0;i<Global::g_numRuns;i++) delete [] mpp_data[i];
        delete [] mpp_data;
        mpp_data=0;
    }
    if(mpp_convergeTime){
        for(int i=0;i<Global::g_numRuns;i++) delete [] mpp_convergeTime[i];
        delete [] mpp_convergeTime;
        mpp_convergeTime=0;
    }
     if(mpp_gOpt){
        for(int i=0;i<Global::g_numRuns;i++) delete [] mpp_gOpt[i];
        delete [] mpp_gOpt;
        mpp_gOpt=0;
    }


}

PerformSingleObj& PerformSingleObj::operator=(const PerformSingleObj& rhs){
    if (this == &rhs) return *this; // handle self assignment
    //assignment operator
    return *this;
}

double PerformSingleObj::getConvergeSpeed2Best(){
    return m_speed2Best;
}
double PerformSingleObj::getConvergeSpeed2GOpt(){
    return m_speed2gOpt;

}
double PerformSingleObj::getConvergeSpeed(){

    return m_speed;

}
void PerformSingleObj::calculateConvergenceTime(){
//track back from the end of mpp_data[i], the convergence time of mpp_convergeTime[i] is the index of mpp_data[i] when the value changes.
//e.g., mpp_data has ten elements: 10 9 8 8 7 6 3 2 2 2. the convergence time is 7 as  mpp_data[7] is different from it's previous one mpp_data[6]
    m_numEvals2Converge=0;
    for(int i=0;i<Global::g_numRuns;i++) {
        for(int j=m_numGopts-1;j>=0;j--){
            for(int k=m_numRecords-(m_numGopts-j-1)*m_recordsPerChange-1;k>=m_numRecords-(m_numGopts-j)*m_recordsPerChange;k--){
                if(k==m_numRecords-(m_numGopts-j)*m_recordsPerChange){
                    mpp_convergeTime[i][j]=k;
                    m_numEvals2Converge+=(k-(m_numRecords-(m_numGopts-j)*m_recordsPerChange)+1)*Global::g_sampleFre;
                    break;
                }
                if(mpp_data[i][k]!=mpp_data[i][k-1]){
                    mpp_convergeTime[i][j]=k;
                    m_numEvals2Converge+=(k-(m_numRecords-(m_numGopts-j)*m_recordsPerChange)+1)*Global::g_sampleFre;
                    break;
                }
            }
        }
    }
    //the average number of evaluations before the algorithm converges
    m_numEvals2Converge/=(Global::g_numRuns*m_numGopts);

}
void PerformSingleObj::getStatInfor(double &mean, double &var, double &worst, double &best){
    mean=m_mean; var=m_var;
	if(mpp_gOpt){
		worst=m_largest;
        best=m_smallest;
	}else{
		if(Global::gp_problem->getProblemType()==MIN_OPT){
			worst=m_largest;
			best=m_smallest;
		}else{
			best=m_largest;
			worst=m_smallest;

		}
	}
}
double PerformSingleObj::getSucRate(){
    return m_sucRate;

}
double PerformSingleObj::getPerformance(){
    return m_performance;
}
PerformSingleObj * PerformSingleObj::getPerformSingleObj(){
     if(!PerformSingleObj::msp_perf){
            Throw(Logic_error("the performance metric method is not initialized"));
            return NULL;
        }

    return PerformSingleObj::msp_perf;
}
void PerformSingleObj::initialize(bool rFlag){

    if(PerformSingleObj::msp_perf){
        //Throw(Logic_error("the PerformSingleObj method has been already initialized"));
        return;
    }

    PerformSingleObj::msp_perf=new PerformSingleObj(rFlag);


}
void PerformSingleObj::deletePerformSingleObj(){
    delete PerformSingleObj::msp_perf;
    PerformSingleObj::msp_perf=0;

}

void PerformSingleObj::record(double rObj, const bool rFlag){

    if(!rFlag) return;
    if(Global::gp_problem->getEvaluations()/Global::g_sampleFre>m_numRecords){
        Throw(Out_of_range("ERROR! the number of data to be recorded is out of the memory"));
    }
    if(Global::gp_problem->getEvaluations()%(m_recordsPerChange*Global::g_sampleFre)==1) m_bestSoFar=rObj;

    if(Global::gp_problem->getProblemType()==MIN_OPT){
        if(m_bestSoFar>rObj) m_bestSoFar=rObj;
    }else{
       if(m_bestSoFar<rObj) m_bestSoFar=rObj;
    }

    if(Global::gp_problem->getEvaluations()%Global::g_sampleFre==0){
        mpp_data[Global::g_runIdx][Global::gp_problem->getEvaluations()/Global::g_sampleFre-1]=m_bestSoFar;
    }

}
void  PerformSingleObj::calculatePerformance(){

    calculateConvergenceTime();

    //fitness decrease/increase per evaluation: m_speed=f(0)-f(covergeTime)/covergeTime
    m_speed=0;
    for(int i=0;i<Global::g_numRuns;i++){
         for(int j=m_numGopts-1;j>=0;j--)
            m_speed+=fabs(mpp_data[i][(mpp_convergeTime[i][j]/m_recordsPerChange)*m_recordsPerChange]-mpp_data[i][mpp_convergeTime[i][j]])
            /(Global::g_sampleFre*(1+mpp_convergeTime[i][j]-(mpp_convergeTime[i][j]/m_recordsPerChange)*m_recordsPerChange ));
    }
    m_speed/=(Global::g_numRuns*m_numGopts);

    //the convergence speed to the best result achieved by the algorithm
    m_speed2Best=0;
    for(int i=0;i<Global::g_numRuns;i++){
         for(int j=0;j<m_numGopts;j++){
            double t=0;
            for(int k=m_recordsPerChange*j;k<=mpp_convergeTime[i][j];k++){
                t+=fabs(mpp_data[i][k]-mpp_data[i][mpp_convergeTime[i][j]]);
            }
            m_speed2Best+=t/((mpp_convergeTime[i][j]-m_recordsPerChange*j+1)*Global::g_sampleFre);
         }

    }
    m_speed2Best/=(Global::g_numRuns*m_numGopts);

    //the convergence speed to the global optimum
    if(mpp_gOpt){
        m_speed2gOpt=0;
        for(int i=0;i<Global::g_numRuns;i++){
             for(int j=0;j<m_numGopts;j++){
                double t=0;
                for(int k=m_recordsPerChange*j;k<=mpp_convergeTime[i][j];k++){
                    t+=fabs(mpp_data[i][k]-mpp_gOpt[i][j]);
                }
                m_speed2gOpt+=t/((mpp_convergeTime[i][j]-m_recordsPerChange*j+1)*Global::g_sampleFre);
             }

        }
        m_speed2gOpt/=(Global::g_numRuns*m_numGopts);
    }

    // the mean, smallest and largest errors to the global optimum
    if(mpp_gOpt){
        m_largest=m_smallest=fabs(mpp_data[0][m_recordsPerChange-1]-mpp_gOpt[0][0]);
        m_mean=0;
        double *mean=new double[m_numGopts];

        for(int j=1;j<=m_numGopts;j++){
            double er;
             mean[j-1]=0;
            for(int i=0;i<Global::g_numRuns;i++){
                er=fabs(mpp_data[i][j*m_recordsPerChange-1]-mpp_gOpt[i][j-1]);
                if(m_largest<er) m_largest=er;
                if(m_smallest>er) m_smallest=er;
                mean[j-1]+=er;
            }
            mean[j-1]/=Global::g_numRuns;
            m_mean+=mean[j-1];
        }
        m_mean/=m_numGopts;
        m_var=0;

        for(int j=1;j<=m_numGopts;j++){
             long double var=0;
            for(int i=0;i<Global::g_numRuns;i++){
                long double x=fabs(mpp_data[i][j*m_recordsPerChange-1]-mpp_gOpt[i][j-1])-mean[j-1];
                 var+=x*x;
            }

            var=sqrt(var/(Global::g_numRuns));
            m_var+=var;
        }
        m_var/=m_numGopts;
        delete [] mean;

    }else{
        m_largest=m_smallest=mpp_data[0][m_recordsPerChange-1];
        m_mean=0;
        double *mean=new double[m_numGopts];

        for(int j=1;j<=m_numGopts;j++){
            double er;
             mean[j-1]=0;
            for(int i=0;i<Global::g_numRuns;i++){
                er=mpp_data[i][j*m_recordsPerChange-1];
                if(m_largest<er) m_largest=er;
                if(m_smallest>er) m_smallest=er;
                mean[j-1]+=er;
            }
            mean[j-1]/=Global::g_numRuns;
            m_mean+=mean[j-1];
        }
        m_mean/=m_numGopts;
        m_var=0;

        for(int j=1;j<=m_numGopts;j++){
             long double var=0;
            for(int i=0;i<Global::g_numRuns;i++){
                 long double x=mpp_data[i][j*m_recordsPerChange-1] -mean[j-1];
                 var=var+x*x;
            }

            var=sqrt(var/(Global::g_numRuns));
            m_var+=var;
        }
        m_var/=m_numGopts;
        delete [] mean;

    }

    //the success rate and the average number of evaluations to achieve the given accuracy level
    if(mpp_gOpt){
        m_sucRate=0;
        int numSuc=0;
        m_numEvals2Suc=0;
        switch(Global::gp_problem->getId()){
        case Sphere: case Sphere_Noisy: case S_Sphere: case R_Sphere: case RS_Sphere: case S_Sphere_CEC05:
        case Rastrigin: case Rastrigin_Noisy: case S_Rastrigin: case R_Rastrigin: case RS_Rastrigin:
        case S_Rastrigin_CEC05: case RS_Rastrigin_CEC05: case Weierstrass: case RS_Weierstrass: case R_Weierstrass:
        case RS_Weierstrass_CEC05: case Griewank: case R_Griewank: case RS_Griewank: case Ackley: case Ackley_Noisy:
        case S_Ackley: case R_Ackley: case RS_Ackley: case RS_Ackley_Bound_CEC05: case Step: case Quartic_Noisy:
        case Scaffer_F6: case Rosenbrock: case S_Rosenbrock: case S_Rosenbrock_CEC05: case Schwefel_2_22:
        case  Schwefel_2_22_Noisy: case S_Schwefel_2_22: case R_Schwefel_2_22: case RS_Schwefel_2_22: case Schwefel:
        case Schwefel_Noisy: case S_Schwefel: case R_Schwefel: case RS_Schwefel: case Schwefel_1_2: case Schwefel_1_2_Noisy:
        case S_Schwefel_1_2: case S_Schwefel_1_2_Noisy: case R_Schwefel_1_2: case RS_Schwefel_1_2: case  S_Schwefel_1_2_CEC05:
        case  S_Schwefel_1_2_Noisy_CEC05: case Schwefel_2_21: case Penalized_1: case Penalized_2: case Noncont_Rastrigin:
        case RS_Elliptic_CEC05: case Elliptic: case Com: case R_Com: case Com_CEC05: case H_Com_CEC05: case H_Com_Noisy_CEC05:
        case RH_Com_CEC05: case RH_Com_NarrowBasin_CEC05: case RH_Com_Bound_CEC05: case M_global1: case M_global2:
        case M_global3: case M_global4: case M_global5: case S_Sphere_CEC08: case Schwefel_2_21_CEC08: case S_Rosenbrock_CEC08:
        case S_Rastrigin_CEC08: case S_Griewank_CEC08: case S_Ackley_CEC08: case RW_Gear_Train: case RW_ParEst_FMSoundWaves:

           for(int i=0;i<Global::g_numRuns;i++){
                 for(int j=1;j<=m_numGopts;j++){
                     double accuracy=dynamic_cast <BenchmarkFunction*> (Global::gp_problem)->getAccuracy();
                    if(fabs(mpp_data[i][j*m_recordsPerChange-1]-mpp_gOpt[i][j-1])<= accuracy) {
                        m_sucRate+=1;
                        numSuc++;
                        for(int k=m_recordsPerChange*(j-1);k<=mpp_convergeTime[i][j-1];k++){
                            if(fabs(mpp_data[i][k]-mpp_gOpt[i][j-1])<accuracy){
                                m_numEvals2Suc+=(k-m_recordsPerChange*(j-1)+1);
                                break;
                            }
                        }
                    }

                 }
            }

            break;
        case CompositionDBG_DOP: case RotationDBG_DOP: case MovingPeak_DOP: case DF1_DOP:
            for(int i=0;i<Global::g_numRuns;i++){
                 for(int j=1;j<=m_numGopts;j++){
                     double accuracy=dynamic_cast<DynamicContinuous*>(Global::gp_problem)->getAccuracy();
                     if(fabs(mpp_data[i][j*m_recordsPerChange-1]-mpp_gOpt[i][j-1])<= accuracy) {
                        m_sucRate+=1;
                        numSuc++;
                        for(int k=m_recordsPerChange*(j-1);k<=mpp_convergeTime[i][j-1];k++){
                            if(fabs(mpp_data[i][k]-mpp_gOpt[i][j-1])<accuracy){
                                m_numEvals2Suc+=(k-m_recordsPerChange*(j-1)+1);
                                break;
                            }
                        }
                    }
                 }
            }

            break;
        case Binary_DOP:       break;
            /// TO DO
        }
        m_sucRate/=(Global::g_numRuns*m_numGopts);
        if(numSuc>0)        m_numEvals2Suc=m_numEvals2Suc*Global::g_sampleFre/numSuc;

    }

    // performance of the algorithm according to its convergene speed to GOPT and the best results obtained
    if(mpp_gOpt){
        m_performance=0;
        for(int i=0;i<Global::g_numRuns;i++){
             for(int j=0;j<m_numGopts;j++){
                double t=0,gap;
                if(Global::gp_problem->getProblemType()==MIN_OPT)  gap=fabs(mpp_data[i][m_recordsPerChange*(j+1)-1])+1;
                else gap=fabs(mpp_data[i][m_recordsPerChange*j])+1;

                for(int k=m_recordsPerChange*j;k<m_recordsPerChange*(j+1);k++){
                    if(Global::gp_problem->getProblemType()==MIN_OPT){
                        ///Note: preprocess the data to make sure they are greater than 0 to avoid overflow error
                        t+=1-(mpp_gOpt[i][j]+gap)/(mpp_data[i][k]+gap);
                    }else{
                        t+=1-(mpp_data[i][k]+gap)/(mpp_gOpt[i][j]+gap);
                    }

                }
                 t/=m_recordsPerChange;
                double best;
                int bst_Idx=m_recordsPerChange*(j+1)-1;

                if(Global::gp_problem->getProblemType()==MIN_OPT){
                    best=(mpp_gOpt[i][j]+gap)/(mpp_data[i][bst_Idx]+gap);
                }else{
                    best=(mpp_data[i][bst_Idx]+gap)/(mpp_gOpt[i][j]+gap);
                }
                m_performance+=best/(1+t);
             }

        }
        m_performance/=(Global::g_numRuns*m_numGopts);

    }

}
void PerformSingleObj::addGOpt(double rGOptObj){
    mpp_gOpt[Global::g_runIdx][m_gOptIdx]=rGOptObj;
    m_gOptIdx++;

}
void PerformSingleObj::resetGOptIndx(){
    m_gOptIdx=0;
}
void PerformSingleObj::setFileName(stringstream &rName){
    m_fileName.str(rName.str());
}
void PerformSingleObj::setProgrOutputFlag(bool rFlag){
    m_progrOutputFlag=rFlag;
}
void PerformSingleObj::setConvgProgMode(ConvgProgeMode rMode){
   m_convgMode=rMode;
}
void PerformSingleObj::setAlgParameter(stringstream & rPar){
    m_algPar.str(rPar.str());

}
void PerformSingleObj::outputResult(){

    if(m_progrOutputFlag){
        stringstream ss;
        ss<<"Result//"<<m_fileName.str()<<"Progr.txt";
        ofstream out(ss.str().c_str());
        double *data=new double[Global::g_numRuns];

        for(int i=0;i<m_numRecords;i++){
            for(int j=0;j<Global::g_numRuns;j++) {
				if(mpp_gOpt) data[j]=fabs(mpp_data[j][i]-mpp_gOpt[j][i/m_recordsPerChange]);
				else    data[j]=mpp_data[j][i];
            }
            gQuickSort<double>(0,Global::g_numRuns-1,data);
            double x;
            if(m_convgMode==PROGR_MEAN){
                x=0;
                for(int j=0;j<Global::g_numRuns;j++) x+=data[j];
                out<<(i+1)*Global::g_sampleFre<<" " <<x/Global::g_numRuns<<endl;
            }else if(m_convgMode==PROGR_MEDEAN){
                 out<<(i+1)*Global::g_sampleFre<<" " <<data[Global::g_numRuns/2]<<endl;
            }else if(m_convgMode==PROGR_WORST){
                if(Global::gp_problem->getProblemType()==MIN_OPT) out<<(i+1)*Global::g_sampleFre<<" " <<data[Global::g_numRuns-1]<<endl;
                else  out<<(i+1)*Global::g_sampleFre<<" " <<data[0]<<endl;
            }else if(m_convgMode==PROGR_BEST){
                if(Global::gp_problem->getProblemType()==MIN_OPT) out<<(i+1)*Global::g_sampleFre<<" " <<data[0]<<endl;
                else  out<<(i+1)*Global::g_sampleFre<<" " <<data[Global::g_numRuns-1]<<endl;
            }else{
                //...
            }
        }
        delete []data;
        out.close();
    }

    stringstream ss;
    ss<<"Result//"<<m_fileName.str()<<"Sta.txt";
    ofstream out(ss.str().c_str());
    out<<"Algorithm: "<<gAlgName[Global::g_algNumber]<<endl;
    out<<"Algorithm Parameters: "<<m_algPar.str()<<endl;
    out<<"Number of runs: "<<Global::g_numRuns<<endl;
    out<<"Problem: "<<Global::gp_problem->ma_name<<endl;
    out<<"Problem Parameters: "<<Global::gp_problem->m_proPar.str()<<endl;
    double lb,ub;
    //warning: only for the search range is the same in all dimensions
    Global::gp_problem->getSearchRange<double>(lb,ub,0);
    out<<"Lower bound: "<<lb<<" Upper bound: "<<ub<<endl<<endl;


    double mean, var, best, worst;
    getStatInfor(mean,var,worst,best);
    out<<"Mean: "<<mean<<endl;
    out<<"STD: "<<var<<endl;
    out<<"Best: "<<best<<endl;
    out<<"Worst: "<<worst<<endl;

    out<<"ConvergeSpeed: "<<m_speed<<endl;
    out<<"Number of Evals to converge: "<<m_numEvals2Converge<<endl;
    out<<"CongergeSpeed2Best: "<<m_speed2Best<<endl;
    if(mpp_gOpt){
        out<<"ConvergeSpeed2GOpt: "<<m_speed2gOpt<<endl;
        out<<"Sucess rate: "<<m_sucRate<<endl;
        out<<"Number of evals to success: "<<m_numEvals2Suc<<endl;
        out<<"Performance: "<<m_performance<<endl;
    }
	if( Global::gp_problem->getId()==MovingPeak_DOP){
		out<<"Peaks_found: "<<(double)dynamic_cast<MovingPeak*>(Global::gp_problem)->getPeaksFound()/m_numGopts/Global::g_numRuns<<endl;

		out<<"sucessly tracking rate for "<<dynamic_cast<MovingPeak*>(Global::gp_problem)->getNumberofPeak()<<" peaks"<<endl;

		for(int i=0;i<dynamic_cast<MovingPeak*>(Global::gp_problem)->getNumberofPeak();i++)
		out<<(double)dynamic_cast<MovingPeak*>(Global::gp_problem)->getTrackNumber(i)/m_numGopts/Global::g_numRuns<<endl;
	}
    out.close();

}

double PerformSingleObj::getBestSoFar(){
	return m_bestSoFar;

}
