/*************************************************************************
* 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: 11 May 2011
// Last modified:

#include "CompositionDBG.h"
CompositionDBG* CompositionDBG::msp_cdbg=0;
ComDBGFuncID CompositionDBG::ms_funID;


void CompositionDBG::freeMemory(){
    delete [] mp_stretchSeverity;
	delete [] mp_convergeSeverity;

	for(int i=0;i<msc_numComFuns;i++)
		delete [] mpp_comBoundary[i];
	delete [] mpp_comBoundary;
	delete [] mp_comFunction;
	delete []mp_fmax;
}
CompositionDBG::CompositionDBG(const int rId, const int rDimNumber, const Encoding rEncoding, const int rNumPeaks,
                               const ChangeType rT, const ComDBGFuncID rF,float const rChangingRatio,const bool rFlagDimChange,const bool rFlagNumPeakChange):RealDBG(rId,rDimNumber,rEncoding,rNumPeaks){

    allocateMemory(m_dimNumber,m_numPeaks);
    strcpy(ma_name,"Compositional DBG");
    setComBoundary();
	m_heightNormalizeSeverity=2000.;

	initialize(rT,rF,rChangingRatio,rFlagDimChange,rFlagNumPeakChange);

}

void CompositionDBG::allocateMemory(const int rDimNum, const int rPeaks){

	mpp_comBoundary=new Boundary<double>*[msc_numComFuns];
	for(int i=0;i<msc_numComFuns;i++)
		mpp_comBoundary[i]=new Boundary<double>[rDimNum];
	mp_convergeSeverity =new double[rPeaks];
	mp_stretchSeverity= new double[rPeaks];
	mp_comFunction=new ProblemTag[rPeaks];
	mp_fmax= new double[rPeaks];

}

CompositionDBG::~CompositionDBG(){

    freeMemory();
    CompositionDBG::msp_cdbg=0;

}
void CompositionDBG::setComBoundary(){
	//Sphere=0,Rastrigin,Weierstrass,Griewank,Ackley

	for(int j=0;j<m_dimNumber;j++){
			mpp_comBoundary[0][j].upper=100;
			mpp_comBoundary[0][j].lower=-100;
			mpp_comBoundary[1][j].upper=5;
			mpp_comBoundary[1][j].lower=-5;
			mpp_comBoundary[2][j].upper=0.5;
			mpp_comBoundary[2][j].lower=-0.5;
			mpp_comBoundary[3][j].upper=100;
			mpp_comBoundary[3][j].lower=-100;
			mpp_comBoundary[4][j].upper=32;
			mpp_comBoundary[4][j].lower=-32;
		}
}
void CompositionDBG::getComBoundary(const ProblemTag &f,double &l,double &u,const int rDimIdx)const{
	switch(f){
		case Sphere:
			l=mpp_comBoundary[0][rDimIdx].lower;
			u=mpp_comBoundary[0][rDimIdx].upper;
			break;
		case Rastrigin:
			l=mpp_comBoundary[1][rDimIdx].lower;
			u=mpp_comBoundary[1][rDimIdx].upper;
			break;
		case Weierstrass:
			l=mpp_comBoundary[2][rDimIdx].lower;
			u=mpp_comBoundary[2][rDimIdx].upper;
			break;
		case Griewank:
			l=mpp_comBoundary[3][rDimIdx].lower;
			u=mpp_comBoundary[3][rDimIdx].upper;
			break;
		case Ackley:
			l=mpp_comBoundary[4][rDimIdx].lower;
			u=mpp_comBoundary[4][rDimIdx].upper;
			break;
		default:
			Throw(Logic_error("No the function in the basic component fs"));
			break;
	}
}
CompositionDBG & CompositionDBG::operator=(const CompositionDBG &rCom){
	if(this==&rCom) return *this;
	RealDBG::operator =(rCom);

	setComBoundary();
	gCopy(mp_convergeSeverity,rCom.mp_convergeSeverity,m_numPeaks);
	gCopy(mp_stretchSeverity,rCom.mp_stretchSeverity,m_numPeaks);

	m_heightNormalizeSeverity=rCom.m_heightNormalizeSeverity;
	gCopy(mp_comFunction,rCom.mp_comFunction,m_numPeaks);
    gCopy(mp_fmax,rCom.mp_fmax,m_numPeaks);

	return *this;
}
void CompositionDBG::setRotationMatrix(){
	// for each basic function of dimension n(even number), R=R(l1,l2)*R(l3,l4)*....*R(ln-1,ln), 0<=li<=n
	Matrix I;

	int * d=new int[m_dimNumber];
	gInitializeRandomArray(d,m_dimNumber,false);
	for(int i=0;i<m_numPeaks;i++){
		for(int j=0;j+1<m_dimNumber;j+=2){
			double angle=2*PI*Global::gp_uniformPro->Next();				// random angle for rotation plane of d[j]-d[j+1] from d[j]th axis to d[j+1]th axis
			I.Set_Rotation(d[j],d[j+1],angle);
			if(j==0)  mp_rotationMatrix[i]=I;
			else
				mp_rotationMatrix[i]*I;
			I.Identity();
		}
	}
	delete [] d;
}
void CompositionDBG::setCovergeSevrity( const double *cs){
	gCopy(mp_convergeSeverity,cs,m_numPeaks);
}
void CompositionDBG::setStretchSeverity(){

	for(int i=0;i<m_numPeaks;i++){
		double l,u;
        getComBoundary(mp_comFunction[i],l,u);
        mp_stretchSeverity[i]=mp_convergeSeverity[i]*(((Boundary<double> *)m_searchRange[0])->upper-((Boundary<double> *)m_searchRange[0])->lower)/(u-l);

	}
}
void CompositionDBG::setBasicFunction(const ProblemTag *bf){
	gCopy(mp_comFunction,bf,m_numPeaks);
}

double CompositionDBG::evaluate(double const * x, bool rFlag){


	gCopy(mp_genes,x,m_dimNumber);


	for(int i=0;i<m_numPeaks;i++){ // calculate weight for each function
		mp_width[i]=0;
		for(int j=0;j<m_dimNumber;j++)
			mp_width[i]+=(mp_genes[j]-mpp_peak[i][j])*(mp_genes[j]-mpp_peak[i][j]);
		if(mp_width[i]!=0)	mp_width[i]=exp(-sqrt(mp_width[i]/(2*m_dimNumber*mp_convergeSeverity[i]*mp_convergeSeverity[i])));
	}

	for(int i=0;i<m_numPeaks;i++){ // calculate objective value for each function

		for(int j=0;j<m_dimNumber;j++)	// calculate the objective value of tranformation function i
			mp_genes[j]=(mp_genes[j]-mpp_peak[i][j])/mp_stretchSeverity[i];//((1+fabs(mpp_peak[i][j]/mp_searchRange[j].upper))*
		Matrix m(m_dimNumber,1);
		m.setDataRow(mp_genes,m_dimNumber);

		m*mp_rotationMatrix[i];

		gCopy(mp_genes,m.Get_Data()[0].data,m_dimNumber);
		correctSolution(mp_comFunction[i]);
		mp_fit[i]=selectFun(mp_comFunction[i]);

		mp_fit[i]=m_heightNormalizeSeverity*mp_fit[i]/fabs(mp_fmax[i]);

		gCopy(mp_genes,x,m_dimNumber);
	}
	double sumw=0,wmax;
	wmax=gExtremum(mp_width,m_numPeaks,MAX_OPT);
	for(int i=0;i<m_numPeaks;i++)
		if(mp_width[i]!=wmax)
			mp_width[i]=mp_width[i]*(1-pow(wmax,10));
	for(int i=0;i<m_numPeaks;i++)
		sumw+=mp_width[i];
	for(int i=0;i<m_numPeaks;i++)
		mp_width[i]/=sumw;
	double obj=0;
	for(int i=0;i<m_numPeaks;i++)
		obj+=mp_width[i]*(mp_fit[i]+mp_height[i]);

	if(rFlag&&m_evals%m_changeFre==0) mp_bestSoFar[0]=obj;

	if(rFlag&&mp_bestSoFar[0]>obj) mp_bestSoFar[0]=obj;

    if(rFlag)    m_evals++;

    if(rFlag&&m_evals%m_changeFre==0&&m_evals<Global::g_tEvals) DynamicProblem::change();

	return obj;
}
inline double CompositionDBG::selectFun(const ProblemTag &f){
	double value;
	switch(f){
	case Sphere:
		value=fSphere();
		break;
	case Rastrigin:
		value=fRastrigin();
		break;
	case Weierstrass:
		value=fWeierstrass();
		break;
	case Griewank:
		value=fGriewank();
		break;
	case Ackley:
		value=fAckley();
		break;
	default:
		break;
	}
	return value;
}
inline double CompositionDBG::fAckley(){
	double fitness=0;
	double s1=0,s2=0;
	for(int i=0;i<m_dimNumber;i++){
		s1+=mp_genes[i]*mp_genes[i];
		s2+=cos(2*PI*mp_genes[i]);
	}
	fitness=-20*exp(-0.2*sqrt(s1/m_dimNumber))-exp(s2/m_dimNumber)+20+E;
	return fitness;
}
inline double CompositionDBG::fGriewank(){
	double s1=0,s2=1;
	for(int i=0;i<m_dimNumber;i++){
		s1+=mp_genes[i]*mp_genes[i]/4000.;
		s2*=cos(mp_genes[i]/sqrt((double) (i+1)));
	}
	return s1-s2+1.;
}
inline double CompositionDBG::fRastrigin(){
	double fit=0;
	for(int i=0;i<m_dimNumber;i++)
		fit=fit+mp_genes[i]*mp_genes[i]-10.*cos(2*PI*mp_genes[i])+10.;
	return fit;
}
inline double CompositionDBG::fSphere(){
	double fit=0;
	for(int i=0;i<m_dimNumber;i++)
		fit+=mp_genes[i]*mp_genes[i];
	return fit;
}
inline double CompositionDBG::fWeierstrass(){
	double a=0.5,b=3;
	int kmax=20;
	double fit=0,s=0;
	for(int i=0;i<m_dimNumber;i++)
		for(int k=0;k<=kmax;k++)
			fit+=pow(a,k)*cos(2*PI*pow(b,k)*(mp_genes[i]+0.5));
	for(int k=0;k<=kmax;k++)
			s+=pow(a,k)*cos(2*PI*pow(b,k)*0.5);
	s=s*m_dimNumber;
	return fit-s;
}
void CompositionDBG::correctSolution(const ProblemTag &f){

	double l,u;
	getComBoundary(f,l,u);
	for(int j=0;j<m_dimNumber;j++){
		if(mp_genes[j]>u)  mp_genes[j]=u;
		else if(mp_genes[j]<l)  mp_genes[j]=l;
	}
}

void CompositionDBG::randomChange(){


    for (int i=0;i<m_numPeaks; i++) gCopy(mpp_prePeak[i],mpp_peak[i],m_dimNumber);
	gCopy(mp_preHeight,mp_height,m_numPeaks);
	gCopy(mp_preWidth,mp_width,m_numPeaks);

	//change the global minimum value of each function
	heightStandardChange();
	//change the position of global optimum of each function randomly
	positionStandardChange(0);

    restoreInfor();
	calculateGlobalOptima();
    updateNumberofChanges();
	m_changeType.counter++;
}
void CompositionDBG::recurrentNoisyChange(){
    for (int i=0;i<m_numPeaks; i++) gCopy(mpp_prePeak[i],mpp_peak[i],m_dimNumber);
	gCopy(mp_preHeight,mp_height,m_numPeaks);
	gCopy(mp_preWidth,mp_width,m_numPeaks);

	double initial_angle;
	double height_range=m_maxHeight-m_minHeight;

	double noisy;
	for(int i=0;i<m_numPeaks;i++){
	    if(mp_whetherChange[i]==false) continue;
		initial_angle=(double)getPeriod()*i/m_numPeaks;

		 mp_height[i]=sinValueNoisy(m_changeType.counter,m_minHeight,m_maxHeight,height_range,initial_angle,m_noisySeverity);
	}

	initial_angle=PI*(sin(2*PI*(m_changeType.counter)/m_period)+1)/12.;
	noisy=m_noisySeverity*Global::gp_normalPro->Next();
	positionStandardChange(initial_angle+noisy);

    restoreInfor();
	calculateGlobalOptima();
	 updateNumberofChanges();
	m_changeType.counter++;
}
void CompositionDBG::recurrentChange(){

    for (int i=0;i<m_numPeaks; i++) gCopy(mpp_prePeak[i],mpp_peak[i],m_dimNumber);
	gCopy(mp_preHeight,mp_height,m_numPeaks);
	gCopy(mp_preWidth,mp_width,m_numPeaks);
	double initial_angle;
	double height_range=m_maxHeight-m_minHeight;

	for(int i=0;i<m_numPeaks;i++){
	    if(mp_whetherChange[i]==false) continue;
		initial_angle=(double)m_period*i/m_numPeaks;
		mp_height[i]=m_minHeight+height_range*(sin(2*PI*(m_changeType.counter+initial_angle)/m_period)+1)/2.;
	}
	initial_angle=PI*(sin(2*PI*m_changeType.counter/m_period)+1)/12.;
	positionStandardChange(initial_angle);

    restoreInfor();
	calculateGlobalOptima();
	 updateNumberofChanges();
	m_changeType.counter++;
}
void CompositionDBG::smallStepChange(){
    for (int i=0;i<m_numPeaks; i++) gCopy(mpp_prePeak[i],mpp_peak[i],m_dimNumber);
	gCopy(mp_preHeight,mp_height,m_numPeaks);
	gCopy(mp_preWidth,mp_width,m_numPeaks);

	heightStandardChange();
	positionStandardChange(0);

    restoreInfor();
	calculateGlobalOptima();
	 updateNumberofChanges();
	m_changeType.counter++;
}
void CompositionDBG::largeStepChange(){
    for (int i=0;i<m_numPeaks; i++) gCopy(mpp_prePeak[i],mpp_peak[i],m_dimNumber);
	gCopy(mp_preHeight,mp_height,m_numPeaks);
	gCopy(mp_preWidth,mp_width,m_numPeaks);

	heightStandardChange();
	positionStandardChange(0);

	restoreInfor();
	calculateGlobalOptima();
	 updateNumberofChanges();
	m_changeType.counter++;
}
void CompositionDBG::chaoticChange(){

    for (int i=0;i<m_numPeaks; i++) gCopy(mpp_prePeak[i],mpp_peak[i],m_dimNumber);
	gCopy(mp_preHeight,mp_height,m_numPeaks);
	gCopy(mp_preWidth,mp_width,m_numPeaks);

	for(int i=0;i<m_numPeaks;i++){
	    if(mp_whetherChange[i]==false) continue;
		mp_height[i]=gChaoticValue(mp_height[i],m_minHeight,m_maxHeight);
	}

	positionStandardChange(0);

    restoreInfor();
	calculateGlobalOptima();
	 updateNumberofChanges();
	m_changeType.counter++;
}

void CompositionDBG::changeDimension(){
    /// no need to preserve  previous information, e.g., positions, heights, width....

	CompositionDBG* r_dbg=new CompositionDBG(m_id,m_dimNumberTemp,m_encoding,m_numPeaks,m_changeType.type,ms_funID,m_changePeakRatio,m_flagDimensionChange,m_flagNumPeaksChange);


	if(m_changeType.type==CT_Recurrent||m_changeType.type==CT_RecurrentNoisy){
        r_dbg->setPeriod(m_period);
	}

	r_dbg->parameterSetting(this);
	r_dbg->calculateGlobalOptima();

	freeMemory();
	RealDBG::freeMemory();
	DynamicContinuous::freeMemory();
	Problem::freeMemory();

    Problem::allocateMemory(m_dimNumberTemp);
    DynamicContinuous::allocateMemory(m_dimNumberTemp,m_numPeaks);
    RealDBG::allocateMemory(m_dimNumberTemp,m_numPeaks);
    allocateMemory(m_dimNumberTemp,m_numPeaks);

    m_dimNumber=m_dimNumberTemp;
    setPeriod(m_period);
	*this=*r_dbg;
	delete r_dbg;


}

void CompositionDBG::parameterSetting(Problem * rDP){
	RealDBG::parameterSetting(rDP);
	RealDBG *rP=dynamic_cast<RealDBG*>(rDP);

	CompositionDBG* r_dbg=static_cast<CompositionDBG*>(rP);

    int peaks=m_numPeaks<rP->getNumberofPeak()?m_numPeaks:rP->getNumberofPeak();

	gCopy(mp_comFunction,r_dbg->mp_comFunction,peaks);
//	setRotationMatrix();
	gCopy(mp_convergeSeverity,r_dbg->mp_convergeSeverity,peaks);
	gCopy(mp_stretchSeverity,r_dbg->mp_stretchSeverity,peaks);

}

CompositionDBG * CompositionDBG::getCompositionDBG(){
    if(!CompositionDBG::msp_cdbg){
        Throw(Logic_error("the CompositionDBG problem is not initialized"));
        return NULL;
    }

    return CompositionDBG::msp_cdbg;
}
void CompositionDBG::deleteCompositionDBG(){
    if(CompositionDBG::msp_cdbg){
        delete CompositionDBG::msp_cdbg;
        CompositionDBG::msp_cdbg=0;
    }
}
bool CompositionDBG::isInitialized(){
	if(CompositionDBG::msp_cdbg!=0) return true;
	else return false;
}

void CompositionDBG::initialize(int rDim, int rPeaks, const ChangeType rT, const ComDBGFuncID rF,float const rChangingRatio,const bool rFlagDimChange, const bool rFlagNumPeakChange){
    if(CompositionDBG::msp_cdbg){
        Throw(Logic_error("the CompositionDBG problem has been already initialized"));
        return;
    }
    if((unsigned int)rDim>msc_MaxDimensionNumber|| (unsigned int) rDim<msc_MinDimensionNumber){
        stringstream  oss;
        oss<<"the number of dimensions should be within ["<<msc_MinDimensionNumber<<","<<msc_MaxDimensionNumber<<"]";
        Throw(Out_of_range(oss.str().c_str()));
        return;
    }
    CompositionDBG::msp_cdbg=new CompositionDBG(CompositionDBG_DOP,rDim,C_DECIMAL,rPeaks,rT,rF,rChangingRatio,rFlagDimChange,rFlagNumPeakChange);
}

void CompositionDBG::initialize(const ChangeType rT, const ComDBGFuncID rF,float const rChangingRatio,const bool rFlagDimChange, const bool rFlagNumPeakChange){

    RealDBG::initialize(rT,rFlagDimChange,rFlagNumPeakChange);
    setProblemType(MIN_OPT);
    setNumberofChanges(rChangingRatio);
    CompositionDBG::ms_funID=rF;

    ProblemTag *basic_fun=new ProblemTag[m_numPeaks];

	switch(CompositionDBG::ms_funID){
		case COMDBG_SPHERE:
			for(int i=0;i<m_numPeaks;i++) basic_fun[i]=Sphere;
			break;
		case COMDBG_RASTRIGIN:
			for(int i=0;i<m_numPeaks;i++) basic_fun[i]=Rastrigin;
			break;
		case COMDBG_GRIEWANK:
			for(int i=0;i<m_numPeaks;i++) basic_fun[i]=Griewank;
			break;
		case COMDBG_ACKLEY:
			for(int i=0;i<m_numPeaks;i++) basic_fun[i]=Ackley;
			break;
		case COMDBG_HYBRID:
			basic_fun[0]=Sphere;		basic_fun[5]=Sphere;
			basic_fun[1]=Rastrigin;		basic_fun[6]=Rastrigin;
			basic_fun[2]=Weierstrass;	basic_fun[7]=Weierstrass;
			basic_fun[3]=Griewank;		basic_fun[8]=Griewank;
			basic_fun[4]=Ackley;		basic_fun[9]=Ackley;
			for(int i=10;i<m_numPeaks;i++) basic_fun[i]=Sphere;
			break;
	}
	setBasicFunction(basic_fun);
    delete []basic_fun;


	double *t=new double[m_numPeaks];
	for(int i=0;i<m_numPeaks;i++)t[i]=1.;
	setCovergeSevrity(t);
	setStretchSeverity();
	setRotationMatrix();

	calculateGlobalOptima();


    Matrix m(m_dimNumber,1);
	for(int i=0;i<m_numPeaks;i++){
    for(int j=0;j<m_dimNumber;j++){ // calculate the estimate max value of funciton i
        mp_genes[j]=((Boundary<double> *)m_searchRange[j])->upper;
        mp_genes[j]/=mp_stretchSeverity[i];
    }
    m.setDataRow(mp_genes,m_dimNumber);
    m*mp_rotationMatrix[i];
    gCopy(mp_genes,m.Get_Data()[0].data,m_dimNumber);
    correctSolution(mp_comFunction[i]);
    mp_fmax[i]=selectFun(mp_comFunction[i]);
    if(mp_fmax[i]==0)   Throw(Runtime_error("the estimation max value must be greater not equal to 0"));

	}
	delete [] t;

 }

void  CompositionDBG::changeNumPeaks(){

    CompositionDBG* r_dbg=new CompositionDBG(m_id,m_dimNumber,m_encoding,m_numPeaksTemp,m_changeType.type,ms_funID
                                             ,m_changePeakRatio,m_flagDimensionChange,m_flagNumPeaksChange);

	if(m_changeType.type==CT_Recurrent||m_changeType.type==CT_RecurrentNoisy)
		r_dbg->setPeriod(m_period);
	r_dbg->parameterSetting(this);


	r_dbg->calculateGlobalOptima();

	freeMemory();
	RealDBG::freeMemory();
	DynamicContinuous::freeMemory();

    DynamicContinuous::allocateMemory(m_dimNumber,m_numPeaksTemp);
    RealDBG::allocateMemory(m_dimNumber,m_numPeaksTemp);
    allocateMemory(m_dimNumber,m_numPeaksTemp);

    m_numPeaks=m_numPeaksTemp;
    setPeriod(m_period);
	*this=*r_dbg;
	delete r_dbg;
}



