
/*
 *  The source code of "Enhanced Opposition-Based Differential Evolution for Solving 
 *  High-Dimensional Continuous Optimization Problems (GODE)"
 *
 * 	The GODE code is used for Special Issue of Soft Computing: 
 *	A Fusion of Foundations, Methodologies and Applications on 
 *	Scalability of Evolutionary Algorithms and other Metaheuristics 
 *	for Large Scale Continuous Optimization Problems.
 *
 *	If you have any questions about it, pleas contact 
 *	Dr. Hui Wang <wanghui_cug@yahoo.com.cn>.
 *
*/

#include "global.h"


// define the individual class
class Indiv 
{
public:
	Indiv();
	virtual ~Indiv();
	void CopyObj(Indiv ob);  // copy function
	double x[D];    // position
	long double fit;   //fitness
	long double Func(int nf);  // evaluation function
 	double low[NF+1];     // lower boundary
 	double up[NF+1];    // upper boundary
};

//---------global  variables-------------------------
long	rnd_uni_init;                 /* serves as a seed for rnd_uni()   */
int		eval;  // the number of function evaluations (FEs)
int		seed;   // random seed
int		r1, r2, r3;					/*  three random indexes for DE    */

// boundaries definition for GOBL
double a[D];
double b[D];
//------------------------- average results------------
//long double avg_fit[Run][200];
int avg_eval=0;
int last_eval = 0;
long double bestvalue[Run];
//----------------------------------
double k = 1.0;       // a parameter for GOBL, and it is randomly generated in the evolution
double p_o = 0.05;    // the rate of opposition

double scale = 0.5;  // F
double CR = 0.9;     // CR

Indiv pop[ps];   // current population
Indiv gopop[ps];   // the opposite population  
Indiv Best;        // the best individual found so far
//--------------------------------save results into an array
long double n_results[Run];  // save the best fitness value for each run
long double avg_n_results[Run][100]; // every 50*D evals, we save the current best result.
double avg_curve[100];
long double best, median, worst, mean, sd;
int n_run = 0;

//------------------------------------- output datas to files
ofstream out; // output the results
//----------------------------------function declare
void Initialize_Seeds();  // initialize the seeds
double rnd_uni(long *idum);   //  retrun a random value between 0 and 1
double rand_double(double l, double u); //  retrun a random value between l and u
void InitPop(int nf);              // population initialization
int IncreaseEval();     // FEs = FEs + 1
void SaveResult();
int FindBest();
void DE_Generation(int nf);  // DE
void SelectR(int candidate);   // Randomly select three different individuals
void Print(int nf);
void Rand1Bin(int candidate, int nf);      // DE/rand/1/bin
void Rand1Exp(int candidate, int nf);    // DE/rand/1/exp
void GOBL(int nf);               // generalized opposition-based learning
void UpdateAB();              //update the boundaries of current population
void Select_n_Best();             // select n best individuals
void UpdateBest();                // update the best individuals
void Run_All();           // Run GODE on all test functions over 25 trials.
void Run_One_Function();    // Run GODE on one function
void Statistics();        // calculate mean, maximum, minimum, median, etc.
void Show_Best_Solution();   // print the best solution
void GODE_Generation(int nf);   // the framework of GODE
//---------------------------------------------



void main()
{
	Initialize_Seeds(); 	// initialize the seeds
//	Run_One_Function();  // Run GODE on one function, please modify the function ID in global.h
	Run_All();   // Run GODE on 19 functions
}

void Initialize_Seeds(){
	srand((unsigned)time(NULL));
	do{seed=rand();}while(seed==0);
	rnd_uni_init = -(long)seed;
}

void Run_One_Function(){
   InitPop(F);
   GODE_Generation(F);
}

// population initialization
void InitPop(int nf){
	int i,j;
	eval = 0;
	for(i=0; i<ps; i++){
		for(j=0; j<D; j++)
			pop[i].x[j] = rand_double(pop[i].low[nf], pop[i].up[nf]);
		//evaluate
		pop[i].Func(nf);
		IncreaseEval(); // FEs = FEs + 1
	}
	//Init gopop
	UpdateAB();
	k = rnd_uni(&rnd_uni_init);
	for(i=0; i<ps; i++){
		for(int j=0; j<D; j++){
			gopop[i].x[j] =  k*(a[j]+b[j]) - pop[i].x[j];		
			//	= k*(pop[i].low[nf]+pop[i].up[nf])-pop[i].x[j];  // opposition
			if(gopop[i].x[j] < pop[i].low[nf] || gopop[i].x[j] > pop[i].up[nf])
				gopop[i].x[j] = rand_double(pop[i].low[nf], pop[i].up[nf]);
		}
		gopop[i].Func(nf);
		IncreaseEval(); // FEs = FEs + 1
	}
//	Select_n_Best();  // select n fittest indivs
	int best_id = FindBest();
	Best.CopyObj(pop[best_id]);
//	cout<<"best_id = "<<best_id<<"	best_fit = "<<pop[best_id].fit<<"	"<<Best.fit<<endl;
}
// DE 
void DE_Generation(int nf){
	while(eval < MAX_FEs){
		int candidate;
		for(candidate=0; candidate<ps; candidate++){
		//	Rand1Bin(candidate, nf);
			Rand1Exp(candidate, nf);
		}
		if(eval % 10000 == 0)
  		 	Print(nf);
	}
}
//  rand/1/exp
void Rand1Exp(int candidate, int nf){
	int n, i;
	SelectR(candidate);
	//Create trial
	Indiv *trial = new Indiv();
	trial->CopyObj(pop[candidate]);
	n = (int)(rnd_uni(&rnd_uni_init)*D);
	for(i=0; (rnd_uni(&rnd_uni_init) < CR) && (i < D); i++){
		trial->x[n] = pop[r1].x[n] + scale*(pop[r2].x[n] - pop[r3].x[n]); 
		//check boundary
		if(trial->x[n] < trial->low[nf])
			trial->x[n] = trial->low[nf];
		if(trial->x[n] > trial->up[nf])
			trial->x[n] = trial->up[nf];
		n = (n+1)%D;
	}
	trial->Func(nf);   // evaluate the trail vector
	IncreaseEval();   // FEs = FEs + 1
	//compare the trial vector with the current vector
	if(trial->fit < pop[candidate].fit){
		pop[candidate].CopyObj(*trial);
		if(trial->fit< Best.fit){
			Best.fit = trial->fit;
			Best.CopyObj(*trial);  //update the best
		}
	}
	delete trial;
}
// rand/1/bin
void Rand1Bin(int candidate, int nf){
	int L, n;
	//randomly generate 3 integers, r1, r2 and r3
	SelectR(candidate);
	//Create trial vector
	Indiv *trial = new Indiv();
	trial->CopyObj(pop[candidate]);
	n = (int)(rnd_uni(&rnd_uni_init)*D);
	for(L=0; L<D; L++){
		if((rnd_uni(&rnd_uni_init) < CR) || (L == (D-1))){
			trial->x[n] = pop[r1].x[n] + scale*(pop[r2].x[n] - pop[r3].x[n]);
			if(trial->x[n] < trial->low[nf] || trial->x[n] > trial->up[nf]) 		//check boundary
				trial->x[n] = pop[candidate].x[n];			
		}			
		n = (n+1)%D;		
	}
	trial->Func(nf); // evaluate the trail vector
	IncreaseEval();  // FEs = FEs + 1
	//compare the trial vector with the current vector
	if(trial->fit < pop[candidate].fit){
		pop[candidate].CopyObj(*trial);
		if(trial->fit< Best.fit){
			Best.fit = trial->fit;
			Best.CopyObj(*trial);  //update the best
		}
	}
	delete trial;
}
// generalized opposition-based learning
void GOBL(int nf){
	int i,j;
	k = rnd_uni(&rnd_uni_init);
	for(i=0; i<ps; i++){
		for(j=0; j<D; j++){
			gopop[i].x[j] = k*(a[j]+b[j]) - pop[i].x[j];
			if(gopop[i].x[j] < pop[i].low[nf] || gopop[i].x[j] > pop[i].up[nf])
				gopop[i].x[j] = rand_double(pop[i].low[nf], pop[i].up[nf]);
		}
		gopop[i].Func(nf);
		IncreaseEval(); // FEs = FEs + 1
	}
}
// The framework of GODE
void GODE_Generation(int nf){
	int candidate;
	while(eval < MAX_FEs){
		//generalized oppostion-based learning
		if(rnd_uni(&rnd_uni_init) < p_o){
		//update interval of the definition
			UpdateAB();
			GOBL(nf);
			//select n best indev from gopop and pop
			Select_n_Best();
		}
		else{
			// DE with rand/1/exp strategy
			for(candidate=0; candidate<ps; candidate++)
				Rand1Exp(candidate, nf);
		}
		UpdateBest();
 		if(eval % (100*D) == 0)
  	   		Print(nf);
		if(Best.fit <=1e-14) 
			Best.fit = 0.0;
	} //end-while
	Print(nf);
}

void UpdateBest(){
	int id = FindBest();
	if(pop[id].fit < Best.fit)
		Best.CopyObj(pop[id]);
}
// select N_p best individuals from current population and the opposite pop
void Select_n_Best(){
	int i,j, index;
	long double val=pop[0].fit;
	Indiv *temp1 = new Indiv[2*ps];
	Indiv *temp2 = new Indiv();
	for(i=0; i<ps; i++){
		temp1[i].CopyObj(pop[i]);
		temp1[i+ps].CopyObj(gopop[i]);
	}
	for(i=0; i<ps; i++){
		val = temp1[i].fit;
		index = i;
		//find the best individual from temp1[i...2*ps]
		for(j=i+1; j<2*ps; j++){
			if(temp1[j].fit < val){
				val = temp1[j].fit;
				index = j;
			}
		}
		//temp1[index] is the best individual
		// change the position between temp1[i] and temp[index]
		temp2->CopyObj(temp1[i]);
		temp1[i].CopyObj(temp1[index]);
		temp1[index].CopyObj(*temp2);
		pop[i].CopyObj(temp1[i]);
	}
	delete []temp1;
	delete temp2;
}
// update the boundaries in the current population
void UpdateAB(){
	int i,j;
	double min1=1e10;
	double max1=-1e10;
	for(i=0; i<D; i++){
		min1=pop[0].x[i];
		max1=pop[0].x[i];
		for(j=1; j<ps; j++){
			if(pop[j].x[i]>max1)
				max1=pop[j].x[i];
			if(pop[j].x[i]<min1)
				min1=pop[j].x[i];
		}
		a[i]=min1;
		b[i]=max1;
	}
}
// print the results
void Print(int nf){
	cout<<"F = "<<nf<<" D = "<<D<<" Run = "<<n_run+1<<" eval = "<<eval<<" p_o = "<<p_o<<" BestFitness = "<<Best.fit<<endl;	
}

// Randomly select three different individuals 
void SelectR(int candidate){
	do{
		r1 = (int)(rnd_uni(&rnd_uni_init)*ps);
	}while(r1==candidate);
	do{
		r2 = (int)(rnd_uni(&rnd_uni_init)*ps);
	}while(r2== r1 || r2==candidate);
	do{
		r3 = (int)(rnd_uni(&rnd_uni_init)*ps);
	}while(r3==r2 || r3==r2 || r3==candidate);
}

//Find the best vector in the current population
int FindBest(){
	long double temp=pop[0].fit;
	int index=0;
	for(int i=1;i<ps;i++){
		if(pop[i].fit<temp){
			temp=pop[i].fit;
			index=i;
		}
	}
	return index;
}

//FEs = FEs + 1;
int IncreaseEval(){
	eval = eval+1;  //FEs = FEs + 1;
	if(eval % (50*D) == 0)
		SaveResult();
	return eval;
}
//save the results
void SaveResult(){
	//out<<Best.fit<<endl;
	avg_n_results[n_run][eval/(50*D)-1] = Best.fit; // eval = k*(50*D) && k>=1
}

// return a random value within [0,1]
double rnd_uni(long *idum){
  long j;
  long k;
  static long idum2=123456789;
  static long iy=0;
  static long iv[NTAB];
  double temp;

  if (*idum <= 0)
  {
    if (-(*idum) < 1) *idum=1;
    else *idum = -(*idum);
    idum2=(*idum);
    for (j=NTAB+7;j>=0;j--)
    {
      k=(*idum)/IQ1;
      *idum=IA1*(*idum-k*IQ1)-k*IR1;
      if (*idum < 0) *idum += IM1;
      if (j < NTAB) iv[j] = *idum;
    }
    iy=iv[0];
  }
  k=(*idum)/IQ1;
  *idum=IA1*(*idum-k*IQ1)-k*IR1;
  if (*idum < 0) *idum += IM1;
  k=idum2/IQ2;
  idum2=IA2*(idum2-k*IQ2)-k*IR2;
  if (idum2 < 0) idum2 += IM2;
  j=iy/NDIV;
  iy=iv[j]-idum2;
  iv[j] = *idum;
  if (iy < 1) iy += IMM1;
  if ((temp=AM*iy) > RNMX) return RNMX;
  else return temp;
}

// random function
double rand_double(double l, double u){
	double r;
	r = l+(u-l)*rnd_uni(&rnd_uni_init);
	return r;
}
//---------------------------------------------------------
// Run GODE on 19 functions over 25 trials
void Run_All(){
//#define		Run			25
	int i, j, h;
	char buffer[50];
	clock_t starttime, endtime;
	double totaltime, avgtime;
	ofstream out_all;
	sprintf(buffer, "D_%d_All_Results.txt", D); 
	out_all.open(buffer);
	ofstream out_mean, out_max, out_min, out_median, out_sd, out_time;
	sprintf(buffer, "D_%d_Mean.txt", D); 
	out_mean.open(buffer);
	sprintf(buffer, "D_%d_Max.txt", D); 
	out_max.open(buffer);
	sprintf(buffer, "D_%d_Min.txt", D); 
	out_min.open(buffer);
	sprintf(buffer, "D_%d_Median.txt", D); 
	out_median.open(buffer);
	sprintf(buffer, "D_%d_SD.txt", D); 
	out_sd.open(buffer);
	sprintf(buffer, "D_%d_Time.txt", D); 
	out_time.open(buffer);
	//run GODE on each function
	// i=1,2,...,19
	cout<<"------GODE for Large Scale Function Optimization------"<<endl;
	for(i = 1; i <= 19; i++){
		// Time start
		starttime = clock();
		//--------------------------------
		for(h=0; h<Run; h++)
			for(j=0; j<100; j++)
				avg_n_results[h][j] = 0.0;
	//	cout<<"------GODE for Large Scale Function Optimization------"<<endl;
		sprintf(buffer, "F%d_D_%d_GODE.txt", i, D); 
		out.open(buffer);
		cout<<"D = "<<D<<endl;
		cout<<"p_o = "<<p_o<<"  N_p = "<<ps<<endl;
		cout<<"F = "<<scale<<"  CR = "<<CR<<endl;
		// GODE on the ith function
		cout<<"--------------------------Function = "<<i<<"--------------------------"<<endl;
		// for ith function, conduct GODE over Run times
		for(n_run = 0; n_run < Run; n_run++){
			cout<<"DE"<<endl;
			// the n_run th running
			cout<<"Run = "<<n_run+1<<endl;
			InitPop(i);
			GODE_Generation(i);
			//save the best results for each run
			n_results[n_run] = Best.fit;
		//	Show_Best_Solution();
		}//end-n_run--
		// Time end
		endtime = clock();
		totaltime = (double)(endtime - starttime);
		avgtime = totaltime/((double)CLOCKS_PER_SEC*Run);
		//staticstic 
		Statistics();
		cout<<"Average Computation Time = "<<avgtime<<" s"<<endl;
		cout<<"Minimum error = "<<best<<endl;
		cout<<"Median error = "<<median<<endl;
		cout<<"Maximum error = "<<worst<<endl;
		cout<<"Mean = "<<mean<<endl;
		cout<<"Std Dev = "<<sd<<endl;
		//output the results to files
		//--------------------------------------------------
		out_mean<<mean<<endl;
		out_max<<worst<<endl;
		out_min<<best<<endl;
		out_median<<median<<endl;
		out_sd<<sd<<endl;
		out_time<<avgtime<<endl;
		//--------------------------------------------------
		// Time, best, median, worst, mean, sd, and avg_curve[]
		out<<"Function = "<<i<<endl;
		out<<"D = "<<D<<endl;
		out<<"-------------------"<<endl;
		out<<"p_o = "<<p_o<<"  N_p = "<<ps<<endl;
		out<<"F = "<<scale<<"  CR = "<<CR<<endl;
		out<<"------------------------------------------------"<<endl;
		out<<"Average Computation Time = "<<avgtime<<" s"<<endl;
		out<<"Minimum error = "<<best<<endl;
		out<<"Median error = "<<median<<endl;
		out<<"Maximum error = "<<worst<<endl;
		out<<"Mean = "<<mean<<endl;
		out<<"Std Dev = "<<sd<<endl;
		out<<"-------------------"<<endl;
		out<<best<<endl;
		out<<median<<endl;
		out<<worst<<endl;
		out<<mean<<endl;
		out<<sd<<endl;
		// output the ith function
		out_all<<"------------------------------------------------"<<endl;
		out_all<<"----------Function = "<<i<<"	D = "<<D<<endl;
		out_all<<"p_o = "<<p_o<<"  N_p = "<<ps<<endl;
		out_all<<"----------F = "<<scale<<"  CR = "<<CR<<endl;
		out_all<<"Minimum error = "<<best<<endl;
		out_all<<"Median error = "<<median<<endl;
		out_all<<"Maximum error = "<<worst<<endl;
		out_all<<"Mean = "<<mean<<endl;
		out_all<<"Std Dev = "<<sd<<endl;
		out_all<<"Average Computation Time = "<<avgtime<<" s"<<endl;
		out_all<<"-------------------"<<endl;
		out_all<<best<<endl;
		out_all<<median<<endl;
		out_all<<worst<<endl;
		out_all<<mean<<endl;
		out_all<<sd<<endl;
		out_all<<avgtime<<endl;
		out<<"------------------------Results achieved on "<<Run<<" Runs------------------------"<<endl;
		out<<endl;
		for(j=0; j<Run; j++)
			out<<n_results[j]<<endl;
		out<<endl;
		out<<"------------------------Evolutionary Curve------------------------"<<endl;
		out<<endl;
		for(j=0; j<100; j++)
			out<<avg_curve[j]<<endl;
		out.close();
	}//end - ith-function
	out_all.close();
	out_mean.close();
	out_max.close();
	out_min.close();
	out_median.close();
	out_sd.close();
	out_time.close();
}

// calculate average, maximum, minimum, median, the standard deviation, etc.
void Statistics(){
	int i,j;
	long double temp;
	long double sum = 0.0;
	for(i=0; i<Run; i++){
		for(j=i+1; j<Run; j++){
			if(n_results[j] < n_results[i]){
				temp = n_results[i];
				n_results[i] = n_results[j];
				n_results[j] = temp;
			}
		}// so n_results[0] saves the best results over Run runs
		sum += n_results[i];
	}
	best = n_results[0];
	median = n_results[Run/2];
	worst = n_results[Run-1];
	mean = sum/Run;
	sd = 0.0;
	for(i=0; i<Run; i++)
		sd +=(n_results[i]-mean)*(n_results[i]-mean);
	sd /= (Run-1);
	sd=sqrt(sd);	
	// average avg_n_results[Run][100]
	for(j=0; j<100; j++){
		temp = 0.0;
		for(i=0; i<Run; i++)
			temp += avg_n_results[i][j];
		temp = temp/Run;
		avg_curve[j] = temp;		
	}
}

// print the best solution
void Show_Best_Solution(){
	for(int i=0; i <D; i++){
		cout<<Best.x[i]<<" ";
		if(i%3 == 0)
			cout<<endl;
	}
}


//------------   class indiv   ------------

Indiv::Indiv(){
	low[0]=0; up[0]= 0;
	// boundarise settings (definition domains for the 19 test problems)
	//----------------------------f1-f6
	low[1]=-100; up[1]= 100;
	low[2]=-100; up[2]= 100;
	low[3]=-100; up[3]=100;
	low[4]=-5; up[4]=5;
	low[5]=-600; up[5]=600;
	low[6]=-32; up[6]=32;
	//---------------f7-f11
	low[7]=-10; up[7]=10;
	low[8]=-65.536; up[8]=65.536;
	low[9]=-100; up[9]=100;
	low[10]=-15; up[10]=15;
	low[11]=-100; up[11]=100;
	//------------------------ f12--f19
	low[12]=-100; up[12]= 100;
	low[13]=-100; up[13]= 100;
	low[14]=-5; up[14]=5;
	low[15]=-10; up[15]=10;
	low[16]=-100; up[16]=100;
	low[17]=-100; up[17]=100;
	low[18]=-5; up[18]=5;
	low[19]=-10; up[19]=10;
}

Indiv::~Indiv(){

}

// copy function
void Indiv::CopyObj(Indiv ob){
	for(int i=0; i<D; i++)
		x[i] = ob.x[i];
	fit = ob.fit;
}

long double Indiv::Func(int nf){
	int f = nf;
	switch(f){
		//-----------------f1-f6
		case 1 : fit = Shifted_Sphere(x); 
			break;
		case 2 : fit = Schwefel_Problem(x); 
			break;
		case 3 : fit = Shifted_Rosenbrock(x); 
			break;
		case 4 : fit = Shifted_Rastrigin(x);
			break;
		case 5 : fit = Shifted_Griewank(x);
			break;
		case 6 : fit = Shifted_Ackley(x);
			break;
		//-----------------f7-f11
		case 7 : fit = f_Schwefel2_22(D, x); 
			break;
		case 8 : fit = f_Schwefel1_2(D, x);
			break;
		case 9 : fit = Extended_f_10(D, x);
			break;
		case 10 : fit = f_Bohachevsky(D, x);
			break;
		case 11 : fit = f_Schaffer(D, x);
			break;
		//----------------f12-f19
		case 12 : fit = f_Hybrid_12(D, x); 
			break;
		case 13 : fit = f_Hybrid_13(D, x); 
			break;
		case 14 : fit = f_Hybrid_14(D, x); 
			break;
		case 15 : fit = f_Hybrid_15(D, x); 
			break;
		case 16 : fit = f_Hybrid_16new(D, x); 
			break;
		case 17 : fit = f_Hybrid_17new(D, x); 
			break;
		case 18 : fit = f_Hybrid_18new(D, x); 
			break;
		case 19 : fit = f_Hybrid_19new(D, x); 
			break;
		default: break;
	}
	return fit;	
}
