#include <math.h>
#include <stdlib.h>
#include <limits.h>
#include <stdio.h>
#include <ctype.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>

/**  Programme d'évaluation de la stratégie MOSES pour des fonctions tests standards f_i **/
/** Default : f_3 probleme de rastrigin **/
/** Repete l'optimisation 100 fois et reporte le nombre d'evaluations de la fonction pour un seuil donné **/


/*****  taille de la population				  ***/

#define n 20

/*****  borne inf et sup du domaine de recherche 	  ***/
/* Rastrigin */

/**/
#define inf -5.0
#define sup  5.0
/**/

/* rosenbroock
#define inf -30.0
#define sup  30.0
*/

/*****  dimension du probleme = nombre de coordonnees 	  ***/

int dim ;

/***** A = ensemble de la population - A[i] = individu i  ***/

float *A[n];
float g[n];


/***** fonction a minimiser : Rastrigin-Rosenbroock	  ***/

float f(x)
float *x;
{
float E, E1, E2,x_2,y_2;
int i;

/****** Sphere f_1 ***************
E = 0.0 ;
 for (i = 0; i < dim  ; i++)
     E += x[i]*x[i];
***************/

/****** Rastrigin f_3 ***************/
 E = 0.0 ;
 for (i = 0; i < dim  ; i++)
     E += x[i]*x[i]  - 10.0 * cos(2.0 * 3.14159 * x[i]) + 10.0 ;
/*********************************/

/****** Rosenbrook f_2 ***************
 E = 0.0 ;
 for (i = 0; i < dim -1 ; i++)
	E += (x[i] - 1.0)*(x[i] - 1.0) +
	100.0*( x[i]* x[i] - x[i+1])*( x[i]* x[i] - x[i+1]);
*********************************/

/****** f_4 ***************
 E1 = 0.0 ; E2 = 0.0;
 for (i = 0; i < dim  ; i++)
     {E1 += x[i]*x[i]; E2+= cos(2.0 * 3.14159 * x[i]);}

 E = -20.0 * exp( -0.2 * sqrt( E1/(float) dim) )
 - exp( E2/(float) dim ) + 20.0 +
 exp(1.0);

****************************/

/******** f_5***************

 x_2 = x[0]*x[0];
 y_2 = x[1]*x[1];
 E = 4 *x_2  - 2.1* x_2 *x_2  + 0.333333* x_2 *x_2 *x_2 + x[0]*x[1] - 4.0*y_2 +
 4.0 *y_2 *y_2 + 1.03162;
****/

 return E;
}

float Delta(x,y,i)
float *x,y;
int i;
{
float D1,D2;

 D1 = x[i]*x[i]  - 10.0 * cos(2.0 * 3.14159 * x[i]);
 D2 = y*y  - 10.0 * cos(2.0 * 3.14159 * y);

 return D2 -D1;
}


initab(min,max,taille)
double min, max;
int taille;
{
 int i,j;
 int aux;

 for(i=0;i<taille;i++)
   {
     A[i] = (float *) malloc(sizeof(float)*dim);
     for (j=0;j<dim;j++)
  		A[i][j] = min + (max - min)*drand48();
   }
}



int intcompare(int *i, int *j)
     {
             if (g[*i] > g[*j])
                     return (1);
             if (g[*i] < g[*j])
                     return (-1);
             return (0);
     }


main()
{
  float *annexe[n] ;
  float *best;

  float g_temp[n];

  int permut[n];

  float f_arret , S, S2;
  float pmut, rayon;
  float Mut;
  float min,max, f_min;
  int compteur = 0;
  int i,j,K;
  int Imin;
  int cpt = n-1, Cpt = 0, Cpt2 = 0;
  int flag = 0;

  int Tot = 0;
  float var2 ;

 srand48(time(0));

/* saisie de la dimension du probleme */

printf("Fonction de Rastrigin\n");
printf("Entrez la dimension du problème (entier, 30 par exemple) \n");
scanf("%d",&dim);


for (i = 0; i < n ; i++) annexe[i] = (float *) malloc( sizeof(float)*dim);

best = (float *) malloc( sizeof(float)*dim);


/* test d'arret sur la valeur de la fonction.
 On arrete lorsque l'algorithme donne une evaluation inferieure a f_arret */

printf("Entrez la valeur du test d'arret : f_arret  (0.05 par exemple )  \n");
scanf("%f",&f_arret);
/*f_arret = 0.001 ;*/

/* probabilite de mutation et rayon de mutation */

printf("Entrez la probabilite de mutation (entre 0.0 et 1.0) \n");
scanf("%f",&pmut);
/*pmut = 1.0;*/

printf("Entrez le rayon de recherche local (entre 0.0 et 2.0) - stratégie uniforme \n");
scanf("%f",&rayon);
/*rayon = 0.1 ;*/

 Tot = 0; var2 = 0.0 ;

for ( K = 1 ; K <= 100 ; K++ )
{
/* initialisation de la population */

initab(inf,sup,n);
 for (i = 0 ;i < n ; i++) g[i] = f(A[i]);
f_min = g[n-1];

Cpt2 = 0;

while ( f_min  >  f_arret )
  {

/*********** Tri de la population  ***********************/



	for (i = 0; i < n ; i++) permut[i] =i ;

	qsort(permut, n, sizeof(int), intcompare);



	for (i = 0; i < n ; i++)
		{
		g_temp[i] = g[permut[i]];
		memcpy(annexe[i],A[permut[i]],sizeof(float)*dim);
 		}

         /*memcpy(best,A[permut[0]],sizeof(float)*dim);*/

 	f_min = g_temp[0];

 	for (i = 0; i < n ; i++)
		{
		g[i] = g_temp[i];
		memcpy(A[i],annexe[i],sizeof(float)*dim);
		}





/************  mutations        ****************************/

cpt = 0;
    for (i = 1 ; i < n ; i++)
      {
         if ( drand48() < pmut )
	 {
	 cpt++;
         j = rand() % dim ;
 	 Mut = A[i][j] +  rayon * ( 1.0 - 2.0 * drand48() );
	 if ( Mut < inf)  Mut = 2*inf - Mut ;
	 else if ( Mut > sup)  Mut = 2*sup - Mut ;

	 /*g[i] += Delta(A[i], Mut ,j);*/

	 A[i][j] = Mut ; g[i] = f(A[i]);
	 }
	 else
	 {
	 /*************** choix uniforme parmi les meilleurs ****************/
	 j = rand() % (1+i);
	 memmove(A[i],annexe[j],sizeof(float)*dim);
	 g[i] = g_temp[j];
 	 /********************************************************************/

	 /************** decalage *****************************************
 	 memcpy(A[i],annexe[i-1],sizeof(float)*dim); g[i] = g_temp[i-1];
	 ****************************************************************/


	 }
      }

    /*

   if (Cpt2 > 10000 )
      {
       printf(" Cpt = %d Erreur = %f \n ", Cpt,  f_min );
       Cpt2 = 0;
      }
     */

    compteur++; Cpt+=cpt; Cpt2+=cpt;
  }

 /*printf("Solution proposee : \n");
 S=0.0;S2=0.0;
 for (j=0; j < dim; j++) {S+= fabs(best[j]);S2+=best[j]*best[j];}
 printf("%f %f\n", S/(float)dim, sqrt(S2));

 printf("Evaluation finale : \n");
 printf("%f\n",f_min);
 printf("Nombre moyen d'evaluations : \n");
 printf("%f\n",  (float) Cpt/ (float) K);*/
 printf("Essai %d   Nombre d'évaluations = %d \n", K, Cpt2);	Tot += Cpt2 ; var2 += (float)Cpt2* (float)Cpt2;
 }
 printf(" \n");
 printf("Mombre moyen d'évaluations = %f\n",(float)Tot/100.0);
 /*printf("%f\n",var2);*/
 printf("s.d  = %f\n", sqrt(var2*0.01  - (float)Tot*(float)Tot*0.0001) );

}