/*************************************************************************
* 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 "MovingPeak.h"


MovingPeak *MovingPeak::msp_MPs=0;

void  MovingPeak::freeMemory(){
	int i;

	for (i=0; i< m_numPeaks; i++){
        delete [] mpp_prevMovement[i];
	}
	delete [] mpp_prevMovement;
	delete []  mp_shift ;
	delete []  mp_coveredPeaks ;
	delete [] mp_isTracked;
	delete [] mp_heightOrder;
	delete [] mp_found;
	mpp_prevMovement=0;
	mp_shift=0 ;
	mp_coveredPeaks=0 ;
	mp_isTracked=0;
	mp_heightOrder=0;
	mp_found=0;
}
void MovingPeak::allocateMemory(const int rDimNum, const int rPeaks){
    mp_shift = new double[rDimNum];
	mp_coveredPeaks = new int[rPeaks];
	mpp_prevMovement = new double*[rPeaks];
	mp_isTracked= new int [rPeaks];
	mp_heightOrder=new int [rPeaks];
	mp_found=new bool[rPeaks];
	for (int i=0; i< rPeaks; i++){
        mpp_prevMovement[i] =new double[rDimNum];
	}

}


bool MovingPeak::readData(){

	/***************************************
		//		m_F		Evaluation Function
		//		1		constant_basis_func()
		//		2		five_peak_basis_func()
		//		3		peak_function1()
		//		4		peak_function_cone()
		//		5		peak_function_hilly()
		//		6		peak_function_twin()
     **************************************
		in>>temp>>m_vlength; // distance by which the peaks are moved, severity
		 lambda determines whether there is a direction of the movement, or whether
		they are totally random. For lambda = 1.0 each move has the same direction,
		while for lambda = 0.0, each move has a random direction
		//in>>temp>>lambda;
		//in>>temp>>m_useBasisFunction;  if set to 1, a static landscape (basis_function) is included in the fitness evaluation
	}*/

	m_F=4;
	m_lambda=0;
	m_useBasisFunction=0;
	m_vlength=1.0;
	m_calculateRightPeak = 1; /* saves computation time if not needed and set to 0 */
	m_minHeight = 30.0;
	m_maxHeight = 70.0;
	m_standardHeight = 50.0;
	m_minWidth = 1.0 ;
	m_maxWidth = 12.0;
	m_standardWidth = 0.0;
	m_heightSeverity=7.0; /* severity of height changes, larger numbers  mean larger severity */
	m_widthSeverity = 1.0; /* severity of width changes, larger numbers mean larger severity */
	
	return true;

}
MovingPeak::MovingPeak(const int rId, const int rDimNumber, const Encoding rEncoding, const int rNumPeaks,float const rChangingRatio,const bool rFlagDimChange,const bool rFlagNumPeakChange):DynamicContinuous(rId,rDimNumber,rEncoding,rNumPeaks){
    m_flagDimensionChange=rFlagDimChange;
		setNumPeaksChange(rFlagNumPeakChange);
   // m_flagNumPeaksChange=rFlagNumPeakChange;

	if(!readData()) exit(0);
	m_peaksFound=0;

	allocateMemory(m_dimNumber,m_numPeaks);
	strcpy(ma_name,"Moving Peak Benchmark");
	initialize();
	setNumberofChanges((int)(rChangingRatio*m_numPeaks));

	m_proPar<<"Vlength:"<<m_vlength<<"; ";


}
void MovingPeak::initialize(){
    int i=0,j=0;
    setSearchRange<double>(0,100);
    setProblemType(MAX_OPT);

	for ( i=0; i< m_numPeaks; i++)
	for ( j=0; j< m_dimNumber; j++){
	  mpp_peak[i][j] = 100.0*Global::gp_uniformPro->Next();
	  mpp_prevMovement[i][j] = Global::gp_uniformPro->Next()-0.5;
	}

	if (m_standardHeight <= 0.0){
        for ( i=0; i< m_numPeaks; i++) mp_height[i]=(m_maxHeight-m_minHeight)*Global::gp_uniformPro->Next()+m_minHeight;
	}else{
        for (i=0; i< m_numPeaks; i++) mp_height[i]= m_standardHeight;
	}

	if (m_standardWidth <= 0.0){
        for (i=0; i< m_numPeaks; i++)
            mp_width[i]= (m_maxWidth-m_minWidth)*Global::gp_uniformPro->Next()+m_minWidth;
	}else{
	for (i=0; i< m_numPeaks; i++)
	  mp_width[i]= m_standardWidth;
	}

	calculateGlobalOptima();
	for (i=0; i< m_numPeaks; i++) {
		mp_heightOrder[i]=i;
		mp_found[i]=false;
	}
	gQuickSort(mp_height,mp_heightOrder,0,m_numPeaks-1);

	for ( i=0; i< m_numPeaks; i++) mp_isTracked[i]=0;
	for (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);
}

void MovingPeak::reset(){
    int i=0,j=0;

    for ( i=0; i< m_numPeaks; i++)
	for ( j=0; j< m_dimNumber; j++){
	  mpp_peak[i][j] = 100.0*Global::gp_uniformPro->Next();
	  mpp_prevMovement[i][j] = Global::gp_uniformPro->Next()-0.5;
	}

	if (m_standardHeight <= 0.0){
        for ( i=0; i< m_numPeaks; i++) mp_height[i]=(m_maxHeight-m_minHeight)*Global::gp_uniformPro->Next()+m_minHeight;
	}else{
        for (i=0; i< m_numPeaks; i++) mp_height[i]= m_standardHeight;
	}

	if (m_standardWidth <= 0.0){
        for (i=0; i< m_numPeaks; i++)
            mp_width[i]= (m_maxWidth-m_minWidth)*Global::gp_uniformPro->Next()+m_minWidth;
	}else{
	for (i=0; i< m_numPeaks; i++)
	  mp_width[i]= m_standardWidth;
	}

	calculateGlobalOptima();
    m_changeCounter=0;
	//for ( i=0; i< m_numPeaks; i++) mp_isTracked[i]=0;  no need to be reset to 0 for different runs, accumulate the values of all runs

	for (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);
}

/* free disc space at end of program */
MovingPeak::~MovingPeak()
{
    freeMemory();
    MovingPeak::msp_MPs=0;
}

/* current_peak_calc determines the peak of the current best individual */
void MovingPeak::currentPeakCalc (double *gen)
{
  int i;
  double maximum = -100000.0, dummy;


  m_currentPeak = 0;
  maximum = functionSelection(gen, 0);
  for(i=1; i<m_numPeaks; i++){
    dummy = functionSelection(gen, i);
    if (dummy > maximum){
      maximum = dummy;
      m_currentPeak = i;
    }
  }
}


/* evaluation function */
double MovingPeak::evaluate (double const *gen, bool rFlag)
{

    int i;
    double maximum = -100000.0, dummy;

    for(i=0; i<m_numPeaks; i++)    {
        dummy = functionSelection(const_cast<double *>(gen), i);
        if (dummy > maximum)      maximum = dummy;
    }

    if (m_useBasisFunction) {

    dummy = functionSelection(const_cast<double *>(gen),-1);
    /* If value of basis function is higher return it */
    if (maximum < dummy)
      maximum = dummy;
    }
	if(rFlag) isTracked(gen);
    if(rFlag)    m_evals++;
    if(rFlag&&m_evals%m_changeFre==0&&m_evals<Global::g_tEvals){
		DynamicProblem::change();
		for (i=0; i< m_numPeaks; i++) {
			mp_heightOrder[i]=i;
			mp_found[i]=false;
		}
		gQuickSort(mp_height,mp_heightOrder,0,m_numPeaks-1);
	}
    return(maximum);
}

/* dummy evaluation function allows to evaluate without being counted */
double MovingPeak::dummyEval (double *gen)
{
  int i;
  double maximum = -100000.0, dummy;

  for(i=0; i<m_numPeaks; i++)
    {
    dummy = functionSelection(gen, i);
    if (dummy > maximum)
      maximum = dummy;
    }

  if (m_useBasisFunction) {

    dummy = functionSelection(gen,-1);
    /* If value of basis function is higher return it */
    if (maximum < dummy)
      maximum = dummy;
  }
  return(maximum);
}

/* whenever this function is called, the peaks are changed */
void MovingPeak::randomChange(){
	int i,j;
	double sum, sum2, offset;

	for (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(i=0; i<m_numPeaks; i++){
		if(mp_whetherChange[i]==false) continue;
		/* shift peak locations */
		sum = 0.0;
		for (j=0; j<m_dimNumber; j++){
			mp_shift[j]=Global::gp_uniformPro->Next()-0.5;
			sum += mp_shift[j]*mp_shift[j];
		}
		if(sum>0.0)		sum = m_vlength/sqrt(sum);
		else  sum = 0.0;                        /* only in case of rounding errors */

		sum2=0.0;
		for (j=0; j<m_dimNumber; j++){
			mp_shift[j]=sum*(1.0-m_lambda)*mp_shift[j]+m_lambda*mpp_prevMovement[i][j];
			sum2 += mp_shift[j]*mp_shift[j];
		}
		if(sum2>0.0)sum2 = m_vlength/sqrt(sum2);
		else     sum2 = 0.0;                      /* only in case of rounding errors */

		for(j=0; j<m_dimNumber; j++){
			mp_shift[j]*=sum2;
			mpp_prevMovement[i][j]= mp_shift[j];
			if ((mpp_peak[i][j]+mpp_prevMovement[i][j]) < ((Boundary<double> *)m_searchRange[j])->lower){
				mpp_peak[i][j] = 2.0*((Boundary<double> *)m_searchRange[j])->lower-mpp_peak[i][j]-mpp_prevMovement[i][j];
				mpp_prevMovement[i][j]*=-1.0;
			}
			else if((mpp_peak[i][j]+mpp_prevMovement[i][j]) > ((Boundary<double> *)m_searchRange[j])->upper){
				mpp_peak[i][j]= 2.0*((Boundary<double> *)m_searchRange[j])->upper-mpp_peak[i][j]-mpp_prevMovement[i][j];
				mpp_prevMovement[i][j]*=-1.0;
			}else
			mpp_peak[i][j] += mpp_prevMovement[i][j];
		}
		/* change peak width */

		offset = Global::gp_normalPro->Next()*m_widthSeverity;
		if ((mp_width[i]+offset) < m_minWidth)		mp_width[i] = 2.0*m_minWidth-mp_width[i]-offset;
		else if ((mp_width[i]+offset) > m_maxWidth)	mp_width[i]= 2.0*m_maxWidth-mp_width[i]-offset;
		else	mp_width[i] += offset;
		/* change peak height */
		offset = m_heightSeverity*Global::gp_normalPro->Next();

		if ((mp_height[i]+offset) < m_minHeight)	mp_height[i] = 2.0*m_minHeight-mp_height[i]-offset;
		else if ((mp_height[i]+offset) > m_maxHeight)	mp_height[i]= 2.0*m_maxHeight-mp_height[i]-offset;
		else	mp_height[i] += offset;
	}

	calculateGlobalOptima();
	updateNumberofChanges();

}


/* Basis Functions */

/* This gives a constant value back to the eval-function that chooses the max of them */
double MovingPeak::constantBasisFunc(double *gen)
{
  return 0.0;
}

double MovingPeak::fivePeakBasisFunc(double *gen)
{
  int i,j;
  double maximum = -100000.0, dummy;
  static double basis_peak [5][7]=
{
  {8.0,  64.0,  67.0,  55.0,   4.0, 0.1, 50.0},
  {50.0,  13.0,  76.0,  15.0,   7.0, 0.1, 50.0},
  {9.0,  19.0,  27.0,  67.0, 24.0, 0.1, 50.0},
  {66.0,  87.0,  65.0,  19.0,  43.0, 0.1, 50.0},
  {76.0,  32.0,  43.0,  54.0,  65.0, 0.1, 50.0}
};
  for(i=0; i<5; i++)
    {
      dummy = (gen[0]-basis_peak[i][0])*(gen[0]-basis_peak[i][0]);
      for (j=1; j< m_dimNumber; j++)
	dummy += (gen[j]-basis_peak[i][j])*(gen[j]-basis_peak[i][j]);
      dummy = basis_peak[i][m_dimNumber+1]-(basis_peak[i][m_dimNumber]*dummy);
      if (dummy > maximum)
        maximum = dummy;
    }
  return maximum;
}

/* Peak Functions */

/* sharp peaks */
double MovingPeak::peakFunction1(double *gen, int peak_number)
{
  int j;
  double dummy;

  dummy = (gen[0]-mpp_peak[peak_number][0])*(gen[0]-mpp_peak[peak_number][0]);
  for (j=1; j< m_dimNumber; j++)
    dummy += (gen[j]-mpp_peak[peak_number][j])*(gen[j]-mpp_peak[peak_number][j]);

  return mp_height[peak_number]/(1+mp_width[peak_number]*dummy);
}

double MovingPeak::peakFunctionCone(double *gen, int peak_number)
{
  int j;
  double dummy;

  dummy =  (gen[0]-mpp_peak[peak_number][0])*(gen[0]-mpp_peak[peak_number][0]);
  for (j=1; j< m_dimNumber; j++)    dummy += (gen[j]-mpp_peak[peak_number][j])*(gen[j]-mpp_peak[peak_number][j]);

  if(dummy!=0)
  dummy =mp_height[peak_number]-mp_width[peak_number]*sqrt(dummy);


  return dummy;
}

double MovingPeak::peakFunctionHilly(double *gen, int peak_number)
{
  int j;
  double dummy;

  dummy =  (gen[0]-mpp_peak[peak_number][0])*(gen[0]-mpp_peak[peak_number][0]);
  for (j=1; j< m_dimNumber; j++)
    dummy += (gen[j]-mpp_peak[peak_number][j])*(gen[j]-mpp_peak[peak_number][j]);

  return mp_height[peak_number]- mp_width[peak_number]*dummy-0.01*sin(20.0*dummy);
}

double MovingPeak::peakFunctionTwin(double  *gen, int peak_number) /* two twin peaks moving together */
{
  int j;
  double maximum = -100000.0, dummy;
  static double twin_peak [7] = /* difference to first peak */
  {
    1.0,  1.0,  1.0,  1.0,   1.0, 0.0, 0.0,
  };

  dummy = pow(gen[0]-mpp_peak[peak_number][0],2);
  for (j=1; j< m_dimNumber; j++)
     dummy += pow(gen[j]-mpp_peak[peak_number][j],2);

  dummy = mp_height[peak_number]- mp_width[peak_number]*dummy;

  maximum = dummy;
  dummy = pow(gen[0]-(mpp_peak[peak_number][0]+twin_peak[0]),2);
  for (j=1; j< m_dimNumber; j++)
     dummy += pow(gen[j]-(mpp_peak[peak_number][j]+twin_peak[0]),2);

  dummy =  mp_height[peak_number]+twin_peak[m_dimNumber+1]-((mp_width[peak_number]+twin_peak[m_dimNumber])*dummy);
  if (dummy > maximum)
    maximum = dummy;

  return maximum;
}
double MovingPeak::functionSelection(double  *gen, int peak_number){
	double dummy;
	switch(m_F){
	case 1: {
			dummy=constantBasisFunc(gen);
			break;
		 }
	case 2: {
			dummy=fivePeakBasisFunc(gen);
			break;
		 }
		case 3: {
			dummy=peakFunction1(gen, peak_number);
			break;
		 }
	case 4: {
			dummy=peakFunctionCone(gen, peak_number);
			break;
		 }
	case 5: {
			dummy=peakFunctionHilly(gen, peak_number);
			break;
		 }
	case 6: {
			dummy=peakFunctionTwin(gen, peak_number);
			break;
		 }
	}
	return dummy;
}
/* The following procedures may be used to change the step size over time */


void MovingPeak::changeStepsizeRandom () /* assigns vlength a value from a normal distribution */
{
  m_vlength = Global::gp_normalPro->Next();
}

void MovingPeak::changeStepsizeLinear() /* sinusoidal change of the stepsize, */
{
  static int counter = 1;
  static double frequency = 3.14159/20.0;  /* returns to same value after 20 changes */

  m_vlength = 1+ sin((double)counter*frequency);
  counter ++;
}

int MovingPeak::getRightPeak()  /* returns 1 if current best individual is on highest mpp_peak, 0 otherwise */
{
    bool flag=false;

    for(int i=0;i<m_numPeaks;i++){
        if(mp_globalOptimaIdx[i]==true && m_currentPeak == i){
            flag=true;
            break;
        }
    }

    return flag;
}
void MovingPeak::setVlength(const double s){
	m_vlength=s;

    size_t start, end;
    start=m_proPar.str().find("Vlength:");
    for(int i=start;i<m_proPar.str().size();i++){
        if(m_proPar.str()[i]==';') {
            end=i;
            break;
        }
    }
    stringstream ss;
    ss<<"Vlength:"<<m_vlength<<"; ";
    string result=m_proPar.str();
    result.replace(start,end-start+1, ss.str());
    m_proPar.str(result);

}

MovingPeak * MovingPeak::getMPs(){
	if(!MovingPeak::msp_MPs){
        Throw(Logic_error("the MPB problem is not initialized"));
	}

	return MovingPeak::msp_MPs;
}
void MovingPeak::deleteMPs(){
	delete MovingPeak::msp_MPs;
	MovingPeak::msp_MPs=0;
}

bool MovingPeak::isInitialized(){
	if(MovingPeak::msp_MPs!=0) return true;
	else return false;
}

void MovingPeak::initialize(int rDim, int rPeaks,float const rChangingRatio,const bool rFlagDimChange, const bool flagNumPeakChange){
    if(MovingPeak::msp_MPs){
        Throw(Logic_error("the MPB problem has been already initialized"));
        return;
    }
    MovingPeak::msp_MPs=new MovingPeak(MovingPeak_DOP,rDim,C_DECIMAL,rPeaks,rChangingRatio,rFlagDimChange,flagNumPeakChange);
}
void MovingPeak::changeNumPeaks(){


    MovingPeak* mpb=new MovingPeak(MovingPeak_DOP,m_dimNumber,C_DECIMAL,m_numPeaksTemp,m_changePeakRatio,m_flagDimensionChange,m_flagNumPeaksChange);


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

	freeMemory();
	DynamicContinuous::freeMemory();

	allocateMemory(m_dimNumber,m_numPeaksTemp);
    DynamicContinuous::allocateMemory(m_dimNumber,m_numPeaksTemp);
  
    m_numPeaks=m_numPeaksTemp;
 
	*this=*mpb;
	delete mpb;

}

void MovingPeak::parameterSetting(DynamicProblem * rP){
	DynamicContinuous::parameterSetting(rP);

	MovingPeak *mpb=static_cast<MovingPeak *>(rP);
	int dim=m_dimNumberTemp<rP->getDimNumber()?m_dimNumberTemp:rP->getDimNumber();
	int peaks=m_numPeaks<rP->getNumberofPeak()?m_numPeaks:rP->getNumberofPeak();

	m_F=mpb->m_F;
	m_vlength=mpb->m_vlength; 
	m_lambda=mpb->m_lambda;
	m_useBasisFunction=mpb->m_useBasisFunction;
	m_calculateRightPeak=mpb->m_calculateRightPeak ; 
	m_standardHeight=mpb->m_standardHeight;
	m_standardWidth=mpb->m_standardWidth;

	m_peaksFound==mpb->m_peaksFound;
	gCopy(mp_shift,mpb->mp_shift,dim);
	gCopy(mp_coveredPeaks,mpb->mp_coveredPeaks,peaks);
	gCopy(mp_isTracked,mpb->mp_isTracked,peaks);
	gCopy(mp_heightOrder,mpb->mp_heightOrder,peaks);
	gCopy(mp_found,mpb->mp_found,peaks);


	for (int i=0; i<peaks; i++){
		gCopy(mpp_prevMovement[i],mpb->mpp_prevMovement[i],dim);
  
	}

}
MovingPeak &MovingPeak::operator=(MovingPeak &other){
	if(this==&other) return *this;

	DynamicContinuous::operator=(other);

	m_F=other.m_F;
	m_vlength=other.m_vlength; 
	m_lambda=other.m_lambda;
	m_useBasisFunction=other.m_useBasisFunction;
	m_calculateRightPeak=other.m_calculateRightPeak ; 
	m_standardHeight=other.m_standardHeight;
	m_standardWidth=other.m_standardWidth;
	m_peaksFound=other.m_peaksFound;

	gCopy(mp_shift,other.mp_shift,m_dimNumber);
	gCopy(mp_coveredPeaks,other.mp_coveredPeaks,m_numPeaks);
	gCopy(mp_isTracked,other.mp_isTracked,m_numPeaks);
	gCopy(mp_heightOrder,other.mp_heightOrder,m_numPeaks);
	gCopy(mp_found,other.mp_found,m_numPeaks);

	for (int i=0; i<m_numPeaks; i++){
		gCopy(mpp_prevMovement[i],other.mpp_prevMovement[i],m_dimNumber);
  
	}
	return *this;
}
int MovingPeak::getTrackNumber(int idex){
	return mp_isTracked[idex];
}

bool MovingPeak::isTracked(double const *gen){
	bool flag=false;
	for(int i=0;i<m_numPeaks;i++){
		double dis=0;
    	for(int j=0;j<m_dimNumber;j++) dis+=(gen[j]-mpp_peak[i][j])*(gen[j]-mpp_peak[i][j]);
		dis=sqrt(dis);
		if(dis<=0.1){//Global::g_sigma
			int j=0;
			while(mp_heightOrder[j++]!=i&&j<m_numPeaks);
			if(!mp_found[i]){
				mp_isTracked[j-1]++;
				mp_found[i]=true;
				m_peaksFound++;
			}
			flag=true;
		}
	}

	return flag;
}
int MovingPeak::getPeaksFound(){
	return m_peaksFound;
}
