/*_________________________________________________________

Test functions
log normal adaptive mutation
Selection operator: Tournament
__________________________________________________________*/

\User declarations :
#define SIZE 10
		#define X_MIN -1.
#define X_MAX 1.
#define ITER 50      
#define Abs(x) ((x) < 0 ? -(x) : (x))
#define MAX(x,y) ((x)>(y)?(x):(y))
#define MIN(x,y) ((x)<(y)?(x):(y))
#define SIGMA  1.                     /*  mutation parameter */
#define PI 3.141592654

 
float pMutPerGene=0.1;
#define N_LIM 10
int n = N_LIM;

double (*Fitness)(double *, int);   // pointeur sur la fonction de fitness designee dans argv[1]
double Sphere(double *, int);
double AckleyPath(double *, int);
double Easom(double *, int);
double Griewangk(double *, int);
double Rastrigin(double *, int);
float Rosenbrock(float *, int);
float Schwefel(float *, int);
float Weierstrass(float *, int);

struct gpuOptions{
};

\end

\User functions:
//fitness function
#include <math.h>

inline float Sphere(float x[SIZE], int n) // Ex-DeJong function 1
{
   float ret = 0;
    	for (int i = 0;i<n; i++) ret += x[ i ] * x[ i ];
      	return ret;
}

inline float AckleyPath(float x[SIZE], int n)
{
    //Parameters
    float a = 20.; 
    float b = 0.2; 
    float c = 2*PI;
    float Scale = 32.768; // pour matcher avec les limites [-1,1]
 
    //Function computation
    float sum_x2 = 0.0;
    float sum_cos_cx = 0.0;
    for (int i=0; i<n; i++) {
		sum_x2 += Scale*Scale*x[i]*x[i];
    		sum_cos_cx += cos(c*Scale*x[i]);
		}
    float res = -a * exp( -b * sqrt( 1/(float)n * sum_x2) ) - exp( 1/(float)n * sum_cos_cx ) + a + exp(1);
    return res;
}
inline float Easom(float x[SIZE], int n) 
{
    float res = 1.;
    //function computation
    for (int i=0; i<n; i++) 
	res *= cos(100.*x[i])*exp(-(100.*x[i]-PI)*(100.*x[i]-PI)); 
    return (1.-res); // pour avoir un minimum a 0.
}
inline float Griewangk(float x[SIZE], int n) 
{
   float res, sum = 0., prod =1.;
   float Scale = 600.; // pour matcher avec les limites [-1,1]
 
	for (int i=1; i<=n; i++) {
		prod *= cos( Scale*x[i]/sqrt( (float)i ) );
    		sum += (Scale*Scale*x[i]*x[i]/4000.);
		}
  	res = sum - prod + 1.;
    	return res;
}

inline float Rastrigin(float x[SIZE], int n)
{
   float res = 0.;
   float Scale = 5.12; // pour matcher avec les limites [-1,1]

   for (int i = 0;i<n; i++) 
	res += ((Scale*Scale*x[i]*x[i])-10*cos( 2*PI*Scale*x[i]));
   return (10.*n + res);
}

inline float Rosenbrock(float x[SIZE], int n)
{
   float res = 0.;
   float Scale = 2.048; // pour matcher avec les limites [-1,1]

   for (int i = 0;i<n; i++) 
    	res += 100.*((Scale*x[i+1] - Scale*Scale*x[i]*x[i])*(Scale*x[i+1] 
               - Scale*Scale*x[i]*x[i]))+(1-Scale*x[i])*(1-Scale*x[i]);
   return (res);
}

inline float Schwefel(float x[SIZE], int n)
{
   float res = 0.;
   float Scale = 500.; // pour matcher avec les limites [-1,1]

   for (int i = 0;i<n; i++) 
    	res += (-Scale*x[i]*sin( sqrt(Abs(Scale*x[i]))));
   return ((float)n*418.9829 + res);
}

__device__ __host__ inline float Weierstrass(float x[SIZE], int n)  // Weierstrass multimidmensionnel h = 0.25
{
   float res = 0.;
   float val[SIZE];
   float b=2.;
   float h = 0.25;

   for (int i = 0;i<n; i++) {
	val[i] = 0.;
    	for (int k=0;k<ITER;k++)
		val[i] += pow(b,-(float)k*h) * sin(pow(b,(float)k)*x[i]);
	res += Abs(val[i]);
	}
   return (res);
} 

float gauss()
/* Generates a normally distributed random value with variance 1 and 0 mean.
    Algorithm based on "gasdev" from Numerical recipes' pg. 203. */
{
  int iset = 0;
  float gset = 0.0;
  float v1 = 0.0, v2 = 0.0, r = 0.0;
  float factor = 0.0;

  if (iset) {
        iset = 0;
        return gset;
      	}
  else {    
        do {
            v1 = (float)globalRandomGenerator->randFloat(0.,1.) * 2.0 - 1.0;
            v2 = (float)globalRandomGenerator->randFloat(0.,1.) * 2.0 - 1.0;
            r = v1 * v1 + v2 * v2;
	                }
        while (r > 1.0);
        factor = sqrt (-2.0 * log (r) / r);
        gset = v1 * factor;
        iset = 1;
        return (v2 * factor);
    	}
}
\end

\GenomeClass::display:
	 for( size_t i=0 ; i<SIZE ; i++){
	      //     cout << Genome.x[i] << ":" << Genome.sigma[i] << "|";
	      printf("%.02f:%.02f|",Genome.x[i],Genome.sigma[i]);
	 }	      
\end
\Initialisation function:
  std::cout<<"*************  n: "<<n<<std::endl;
\end

\Finalization function:
cout << "finalization function called" << endl;
\end

\User classes :
GenomeClass { 
  float x[SIZE];
  float sigma[SIZE]; // auto-adaptative mutation parameter
}
\end

\GenomeClass::initialiser : // "initializer" is also accepted
  for(int i=0; i<n; i++ ) {
     	Genome.x[i] = (float)random(X_MIN,X_MAX);
	Genome.sigma[i]=(float)random(0.,0.5);
	}
\end

\GenomeClass::crossover :
  for (int i=0; i<n; i++)
  {
     float alpha = (float)globalRandomGenerator->getRandomIntMax(1.); // barycentric crossover
     child1.x[i] = alpha*parent1.x[i] + (1.-alpha)*parent2.x[i];
     //if (&child2) child2.x[i] = alpha*parent2.x[i] + (1.-alpha)*parent1.x[i];
  }
\end

\GenomeClass::mutator : // Must return the number of mutations
  int NbMut=0;
  float pond = 1./sqrt((float)n);

    for (int i=0; i<n; i++)
    if (tossCoin(pMutPerGene)){
    	NbMut++;
       	Genome.sigma[i] = Genome.sigma[i] * exp(SIGMA*pond*(float)gauss());
       	Genome.sigma[i] = MIN(0.5,Genome.sigma[0]);              
       	Genome.sigma[i] = MAX(0.,Genome.sigma[0]);
       	Genome.x[i] += Genome.sigma[i]*(float)gauss();
       	Genome.x[i] = MIN(X_MAX,Genome.x[i]);              // pour eviter les depassements
       	Genome.x[i] = MAX(X_MIN,Genome.x[i]);
    	}
return NbMut;
\end

\GenomeClass::evaluator : // Returns the score
  float Score= 0.0;
  float Point[SIZE];
  for (int i=0; i<N_LIM; i++) Point[i] = Genome.x[i];
  Score= Weierstrass(Point, N_LIM);         
  return Score;
\end

\User Makefile options: 
CPPFLAGS+=
\end



\Default run parameters :        // Please let the parameters appear in this order
 Number of generations : 100     // NB_GEN
  Mutation probability : 1       // MUT_PROB
  Crossover probability : 1      // XOVER_PROB
  Population size : 20     // POP_SIZE
  Offspring size : 20 // 40%  
  Genitors selector: Tournament 2
  Final reduce: Tournament 3

    Reduce parents: Tournament 2
    Surviving parents : 5

    Reduce offspring: Tournament 2
    Surviving offspring: 17

		    	   //  Selection operator : Tournament // RouletteWheel, Deterministic, Ranking, Random

		   //Replacement strategy : Plus    // Comma, SteadyState, Generational
		   //  Discarding operator : Worst    // Best, Tournament, Parent, Random
  Evaluator goal : Minimise      // Maximise
  Elitism: Strong
  Elite: 1
\end