Commit 893e0e9c authored by Frederic's avatar Frederic

cma-es std

parent ff7c803e
......@@ -23,6 +23,7 @@ Centre de Math
#define DREAM 3
#define CUDA 4
#define STD 5
#define CMAES 6
#define STD_FLAVOR_SO 0
#define STD_FLAVOR_MO 1
......@@ -45,7 +46,8 @@ extern char sREPLACEMENT[], sSELECTOR[], sSELECT_PRM[], sRED_PAR[], sRED_PAR_PRM
extern int nMINIMISE,nELITE;
extern bool bELITISM, bVERBOSE;
extern bool bPROP_SEQ;
extern int nPOP_SIZE, nNB_GEN, nNB_ISLANDS, nOFF_SIZE, nSURV_PAR_SIZE, nSURV_OFF_SIZE;
extern char* nGENOME_NAME;
extern int nPOP_SIZE, nNB_GEN, nNB_ISLANDS, nOFF_SIZE, nSURV_PAR_SIZE, nSURV_OFF_SIZE, nPROBLEM_DIM;
extern float fMUT_PROB, fXOVER_PROB, fREPL_PERC, fMIG_FREQ;
extern int nMIG_CLONE, nNB_MIG, nIMMIG_REPL;
extern char sMIG_SEL[], sMIGRATOR[], sIMMIG_SEL[],sMIG_TARGET_SELECTOR[];
......
......@@ -682,6 +682,8 @@ exponent ([Ee][+-]?[0-9]+)
<TEMPLATE_ANALYSIS>"\\RED_FINAL_PRM" {fprintf(fpOutputFile,"%s",sRED_FINAL_PRM);}
<TEMPLATE_ANALYSIS>"\\SURV_PAR_SIZE" {fprintf(fpOutputFile,"%d",nSURV_PAR_SIZE);}
<TEMPLATE_ANALYSIS>"\\SURV_OFF_SIZE" {fprintf(fpOutputFile,"%d",nSURV_OFF_SIZE);}
<TEMPLATE_ANALYSIS>"\\GENOME_NAME" {fprintf(fpOutputFile,"%s",nGENOME_NAME);}
<TEMPLATE_ANALYSIS>"\\PROBLEM_DIM" {fprintf(fpOutputFile,"%d",nPROBLEM_DIM);}
<TEMPLATE_ANALYSIS>"\\NB_GEN" {fprintf(fpOutputFile,"%d",nNB_GEN);}
<TEMPLATE_ANALYSIS>"\\NB_ISLANDS" {fprintf(fpOutputFile,"%d",nNB_ISLANDS);}
<TEMPLATE_ANALYSIS>"\\PROP_OR_SEQ" {fprintf(fpOutputFile,"%s",bPROP_SEQ?"Prop":"Seq");}
......@@ -1956,6 +1958,8 @@ int CEASEALexer::create(CEASEAParser* pParser, CSymbolTable* pSymTable)
if (TARGET==STD){
if(TARGET_FLAVOR == STD_FLAVOR_SO)
strcat(sTemp,"STD.tpl");
else if (TARGET_FLAVOR == CMAES)
strcat(sTemp,"CMAES.tpl");
else
strcat(sTemp,"STD_MO.tpl");
if (!(yyin = fpTemplateFile = fopen(sTemp, "r"))){
......
......@@ -45,6 +45,8 @@ int nELITE;
bool bELITISM=0;
bool bVERBOSE=0;
int nPOP_SIZE, nOFF_SIZE, nSURV_PAR_SIZE, nSURV_OFF_SIZE;
char *nGENOME_NAME;
int nPROBLEM_DIM;
int nNB_GEN;
int nNB_ISLANDS;
bool bPROP_SEQ;
......@@ -353,6 +355,7 @@ Object
if (bVERBOSE) printf(" %s pointer declared (%d bytes)\n",$2->sName,$2->nSize);
}
| Symbol '[' Expr ']' {
if((TARGET_FLAVOR==CMAES) && nPROBLEM_DIM==0 && strcmp(pCURRENT_CLASS->sName,"Genome")==0) { nGENOME_NAME=$1->sName; nPROBLEM_DIM=(int)$3;}
$1->nSize=pCURRENT_TYPE->nSize*(int)$3;
$1->pClass=pCURRENT_CLASS;
$1->pType=pCURRENT_TYPE;
......@@ -1375,6 +1378,10 @@ int main(int argc, char *argv[]){
TARGET=STD;
TARGET_FLAVOR = STD_FLAVOR_MO;
}
else if (!mystricmp(sTemp,"cmaes")) {
TARGET=STD;
TARGET_FLAVOR = CMAES;
}
else if (!mystricmp(sTemp,"v")) bVERBOSE=true;
else if (!mystricmp(sTemp,"path")) {
......
......@@ -23,4 +23,4 @@ $(LIB_NAME): $(OBJS)
$(CPPC) $(CPPFLAGS) -fpic -c $<
clean:
rm -f *.o $(LIB_NAME)
\ No newline at end of file
rm -f *.o $(LIB_NAME) $(STATIC_LIB_NAME)
\ No newline at end of file
......@@ -63,8 +63,8 @@ clean:
install:$(EXEC)
sudo cp $< /usr/bin/dev-easea
# realclean: clean
# rm -f EaseaParse.cpp EaseaParse.h EaseaLex.cpp EaseaLex.h
realclean: clean
rm -f EaseaParse.cpp EaseaParse.h EaseaLex.cpp EaseaLex.h
EaseaParse.cpp: EaseaParse.y
......
\TEMPLATE_START/**
This is program entry for CUDA template for EASEA
*/
\ANALYSE_PARAMETERS
using namespace std;
#include <iostream>
#include "EASEATools.hpp"
#include "EASEAIndividual.hpp"
#include <time.h>
RandomGenerator* globalRandomGenerator;
int main(int argc, char** argv){
parseArguments("EASEA.prm",argc,argv);
size_t parentPopulationSize = setVariable("popSize",\POP_SIZE);
size_t offspringPopulationSize = setVariable("nbOffspring",\OFF_SIZE);
float pCrossover = \XOVER_PROB;
float pMutation = \MUT_PROB;
float pMutationPerGene = 0.05;
time_t seed = setVariable("seed",time(0));
globalRandomGenerator = new RandomGenerator(seed);
std::cout << "Seed is : " << seed << std::endl;
SelectionOperator* selectionOperator = new \SELECTOR;
SelectionOperator* replacementOperator = new \RED_FINAL;
SelectionOperator* parentReductionOperator = new \RED_PAR;
SelectionOperator* offspringReductionOperator = new \RED_OFF;
float selectionPressure = \SELECT_PRM;
float replacementPressure = \RED_FINAL_PRM;
float parentReductionPressure = \RED_PAR_PRM;
float offspringReductionPressure = \RED_OFF_PRM;
string outputfile = setVariable("outputfile","");
string inputfile = setVariable("inputfile","");
EASEAInit(argc,argv);
EvolutionaryAlgorithm ea(parentPopulationSize,offspringPopulationSize,selectionPressure,replacementPressure,parentReductionPressure,offspringReductionPressure,
selectionOperator,replacementOperator,parentReductionOperator, offspringReductionOperator, pCrossover, pMutation, pMutationPerGene,
outputfile,inputfile);
StoppingCriterion* sc = new GenerationalCriterion(&ea,setVariable("nbGen",\NB_GEN));
ea.addStoppingCriterion(sc);
Population* pop = ea.getPopulation();
ea.runEvolutionaryLoop();
EASEAFinal(pop);
delete pop;
delete sc;
delete selectionOperator;
delete replacementOperator;
delete globalRandomGenerator;
return 0;
}
\START_CUDA_GENOME_CU_TPL
#include "EASEAIndividual.hpp"
#include "EASEAUserClasses.hpp"
#include <string.h>
#include <fstream>
#include <sys/time.h>
#include <math.h>
#define CMAES_TPL
extern RandomGenerator* globalRandomGenerator;
//Déclarations nécéssaires pour cma-es
#define MAX(x,y) ((x)>(y)?(x):(y))
#define MIN(x,y) ((x)<(y)?(x):(y))
typedef struct
{
/* Variables for Uniform() */
long int startseed;
long int aktseed;
long int aktalea;
long int rgalea[32];
/* Variables for Gauss() */
short flgstored;
double hold;
} aleatoire_t;
typedef struct{
//random_t rand; /* random number generator */
double sigma; /* step size */
double *rgxmean; /* mean x vector, "parent" */
double chiN;
double **C; /* lower triangular matrix: i>=j for C[i][j] */
double **B; /* matrix with normalize eigenvectors in columns */
double *rgD; /* axis lengths */
double *rgpc;
double *rgps;
double *rgxold;
double *rgout;
double *rgBDz; /* for B*D*z */
double *rgdTmp; /* temporary (random) vector used in different places */
short flgEigensysIsUptodate;
short flgCheckEigen; /* control via signals.par */
double genOfEigensysUpdate;
short flgIniphase;
int lambda; /* -> mu, <- N */
int mu; /* -> weights, (lambda) */
double mucov, mueff; /* <- weights */
double *weights; /* <- mu, -> mueff, mucov, ccov */
double damps; /* <- cs, maxeval, lambda */
double cs; /* -> damps, <- N */
double ccumcov; /* <- N */
double ccov; /* <- mucov, <- N */
int gen;
double * xstart;
double * typicalX;
int typicalXcase;
double * rgInitialStds;
double * rgDiffMinChange;
struct { int flgalways; double modulo; double maxtime; } updateCmode;
double facupdateCmode;
aleatoire_t alea; /* random number generator */
int seed;
}CMA;
CMA cma;
int nbEnfants;
\INSERT_USER_DECLARATIONS
\ANALYSE_USER_CLASSES
//Functions for cma
long alea_Start( aleatoire_t *t, long unsigned inseed)
{
long tmp;
int i;
t->flgstored = 0;
t->startseed = inseed;
if (inseed < 1)
inseed = 1;
t->aktseed = inseed;
for (i = 39; i >= 0; --i)
{
tmp = t->aktseed/127773;
t->aktseed = 16807 * (t->aktseed - tmp * 127773)
- 2836 * tmp;
if (t->aktseed < 0) t->aktseed += 2147483647;
if (i < 32)
t->rgalea[i] = t->aktseed;
}
t->aktalea = t->rgalea[0];
return inseed;
}
long alea_init(aleatoire_t *t, long unsigned inseed)
{
clock_t cloc = clock();
t->flgstored = 0;
if (inseed < 1) {
while ((long) (cloc - clock()) == 0)
; /* TODO: remove this for time critical applications? */
inseed = (long)abs(100*time(NULL)+clock());
}
return alea_Start(t, inseed);
}
double alea_Uniform( aleatoire_t *t)
{
long tmp;
tmp = t->aktseed/127773;
t->aktseed = 16807 * (t->aktseed - tmp * 127773)
- 2836 * tmp;
if (t->aktseed < 0)
t->aktseed += 2147483647;
tmp = t->aktalea / 67108865;
t->aktalea = t->rgalea[tmp];
t->rgalea[tmp] = t->aktseed;
return (double)(t->aktalea)/(2.147483647e9);
}
double alea_Gauss(aleatoire_t *t)
{
double x1, x2, rquad, fac;
if (t->flgstored)
{
t->flgstored = 0;
return t->hold;
}
do
{
x1 = 2.0 * alea_Uniform(t) - 1.0;
x2 = 2.0 * alea_Uniform(t) - 1.0;
rquad = x1*x1 + x2*x2;
} while(rquad >= 1 || rquad <= 0);
fac = sqrt(-2.0*log(rquad)/rquad);
t->flgstored = 1;
t->hold = fac * x1;
return fac * x2;
}
void CMA_init_param(CMA *t, int lambda, int mu){
double s1, s2;
double t1, t2;
int i;
clock_t cloc = clock();
t->lambda = lambda;
t->mu = mu;
/*set weights*/
t->weights = (double*)malloc(t->mu*sizeof(double));
for (i=0; i<t->mu; ++i)
t->weights[i] = log(t->mu+1.)-log(i+1.);
/* normalize weights vector and set mueff */
s1=0., s2=0.;
for (i=0; i<t->mu; ++i) {
s1 += t->weights[i];
s2 += t->weights[i]*t->weights[i];
}
t->mueff = s1*s1/s2;
for (i=0; i<t->mu; ++i)
t->weights[i] /= s1;
if(t->mu < 1 || t->mu > t->lambda || (t->mu==t->lambda && t->weights[0]==t->weights[t->mu-1])){
printf("readpara_SetWeights(): invalid setting of mu or lambda\n");
exit(0);
}
/*supplement defaults*/
t->cs = (t->mueff + 2.) / (\PROBLEM_DIM + t->mueff + 3.);
t->ccumcov = 4. / (\PROBLEM_DIM + 4);
t->mucov = t->mueff;
t1 = 2. / ((\PROBLEM_DIM+1.4142)*(\PROBLEM_DIM+1.4142));
t2 = (2.*t->mueff-1.) / ((\PROBLEM_DIM+2.)*(\PROBLEM_DIM+2.)+t->mueff);
t2 = (t2 > 1) ? 1 : t2;
t2 = (1./t->mucov) * t1 + (1.-1./t->mucov) * t2;
t->ccov = t2;
//t->stopMaxIter = ceil((double)(t->stopMaxFunEvals / t->lambda));
t->damps = 1;
t->damps = t->damps * (1 + 2*MAX(0., sqrt((t->mueff-1.)/(\PROBLEM_DIM+1.)) - 1)) * 0.3 + t->cs;
t->updateCmode.modulo = 1./t->ccov/(double)(\PROBLEM_DIM)/10.;
t->updateCmode.modulo *= t->facupdateCmode;
while ((int) (cloc - clock()) == 0)
; /* TODO: remove this for time critical applications!? */
t->seed = (unsigned int)abs(100*time(NULL)+clock());
}
static void TestMinStdDevs(CMA *t)
/* increases sigma */
{
int i;
if (t->rgDiffMinChange == NULL)
return;
for (i = 0; i < \PROBLEM_DIM; ++i)
while (t->sigma * sqrt(t->C[i][i]) < t->rgDiffMinChange[i])
t->sigma *= exp(0.05+t->cs/t->damps);
} /* cmaes_TestMinStdDevs() */
int Check_Eigen(int SIZE, double **C, double *diag, double **Q)
{
/* compute Q diag Q^T and Q Q^T to check */
int i, j, k, res = 0;
double cc, dd;
for (i=0; i < SIZE; ++i)
for (j=0; j < SIZE; ++j) {
for (cc=0.,dd=0., k=0; k < SIZE; ++k) {
cc += diag[k] * Q[i][k] * Q[j][k];
dd += Q[i][k] * Q[j][k];
}
/* check here, is the normalization the right one? */
if (fabs(cc - C[i>j?i:j][i>j?j:i])/sqrt(C[i][i]*C[j][j]) > 1e-10 && fabs(cc - C[i>j?i:j][i>j?j:i]) > 3e-14)
{
printf("cmaes_t:Eigen(): imprecise result detected \n");
++res;
}
if (fabs(dd - (i==j)) > 1e-10) {
printf("cmaes_t:Eigen(): imprecise result detected (Q not orthog.)\n");
++res;
}
}
return res;
}
double myhypot(double a, double b)
/* sqrt(a^2 + b^2) numerically stable. */
{
double r = 0;
if (fabs(a) > fabs(b)) {
r = b/a;
r = fabs(a)*sqrt(1+r*r);
} else if (b != 0) {
r = a/b;
r = fabs(b)*sqrt(1+r*r);
}
return r;
}
void Householder2(int n, double **V, double *d, double *e) {
int i,j,k;
for (j = 0; j < n; j++) {
d[j] = V[n-1][j];
}
/* Householder reduction to tridiagonal form */
for (i = n-1; i > 0; i--) {
/* Scale to avoid under/overflow */
double scale = 0.0;
double h = 0.0;
for (k = 0; k < i; k++) {
scale = scale + fabs(d[k]);
}
if (scale == 0.0) {
e[i] = d[i-1];
for (j = 0; j < i; j++) {
d[j] = V[i-1][j];
V[i][j] = 0.0;
V[j][i] = 0.0;
}
} else {
/* Generate Householder vector */
double f, g, hh;
for (k = 0; k < i; k++) {
d[k] /= scale;
h += d[k] * d[k];
}
f = d[i-1];
g = sqrt(h);
if (f > 0) {
g = -g;
}
e[i] = scale * g;
h = h - f * g;
d[i-1] = f - g;
for (j = 0; j < i; j++) {
e[j] = 0.0;
}
/* Apply similarity transformation to remaining columns */
for (j = 0; j < i; j++) {
f = d[j];
V[j][i] = f;
g = e[j] + V[j][j] * f;
for (k = j+1; k <= i-1; k++) {
g += V[k][j] * d[k];
e[k] += V[k][j] * f;
}
e[j] = g;
}
f = 0.0;
for (j = 0; j < i; j++) {
e[j] /= h;
f += e[j] * d[j];
}
hh = f / (h + h);
for (j = 0; j < i; j++) {
e[j] -= hh * d[j];
}
for (j = 0; j < i; j++) {
f = d[j];
g = e[j];
for (k = j; k <= i-1; k++) {
V[k][j] -= (f * e[k] + g * d[k]);
}
d[j] = V[i-1][j];
V[i][j] = 0.0;
}
}
d[i] = h;
}
/* Accumulate transformations */
for (i = 0; i < n-1; i++) {
double h;
V[n-1][i] = V[i][i];
V[i][i] = 1.0;
h = d[i+1];
if (h != 0.0) {
for (k = 0; k <= i; k++) {
d[k] = V[k][i+1] / h;
}
for (j = 0; j <= i; j++) {
double g = 0.0;
for (k = 0; k <= i; k++) {
g += V[k][i+1] * V[k][j];
}
for (k = 0; k <= i; k++) {
V[k][j] -= g * d[k];
}
}
}
for (k = 0; k <= i; k++) {
V[k][i+1] = 0.0;
}
}
for (j = 0; j < n; j++) {
d[j] = V[n-1][j];
V[n-1][j] = 0.0;
}
V[n-1][n-1] = 1.0;
e[0] = 0.0;
} /* Housholder() */
void QLalgo2 (int n, double *d, double *e, double **V) {
int i, k, l, m;
double f = 0.0;
double tst1 = 0.0;
double eps = 2.22e-16; /* Math.pow(2.0,-52.0); == 2.22e-16 */
/* shift input e */
for (i = 1; i < n; i++) {
e[i-1] = e[i];
}
e[n-1] = 0.0; /* never changed again */
for (l = 0; l < n; l++) {
/* Find small subdiagonal element */
if (tst1 < fabs(d[l]) + fabs(e[l]))
tst1 = fabs(d[l]) + fabs(e[l]);
m = l;
while (m < n) {
if (fabs(e[m]) <= eps*tst1) {
/* if (fabs(e[m]) + fabs(d[m]+d[m+1]) == fabs(d[m]+d[m+1])) { */
break;
}
m++;
}
/* If m == l, d[l] is an eigenvalue, */
/* otherwise, iterate. */
if (m > l) {
int iter = 0;
do { /* while (fabs(e[l]) > eps*tst1); */
double dl1, h;
double g = d[l];
double p = (d[l+1] - g) / (2.0 * e[l]);
double r = myhypot(p, 1.);
iter = iter + 1; /* Could check iteration count here */
/* Compute implicit shift */
if (p < 0) {
r = -r;
}
d[l] = e[l] / (p + r);
d[l+1] = e[l] * (p + r);
dl1 = d[l+1];
h = g - d[l];
for (i = l+2; i < n; i++) {
d[i] -= h;
}
f = f + h;
/* Implicit QL transformation. */
p = d[m];
{
double c = 1.0;
double c2 = c;
double c3 = c;
double el1 = e[l+1];
double s = 0.0;
double s2 = 0.0;
for (i = m-1; i >= l; i--) {
c3 = c2;
c2 = c;
s2 = s;
g = c * e[i];
h = c * p;
r = myhypot(p, e[i]);
e[i+1] = s * r;
s = e[i] / r;
c = p / r;
p = c * d[i] - s * g;
d[i+1] = h + s * (c * g + s * d[i]);
/* Accumulate transformation. */
for (k = 0; k < n; k++) {
h = V[k][i+1];
V[k][i+1] = s * V[k][i] + c * h;
V[k][i] = c * V[k][i] - s * h;
}
}
p = -s * s2 * c3 * el1 * e[l] / dl1;
e[l] = s * p;
d[l] = c * p;
}
/* Check for convergence. */
} while (fabs(e[l]) > eps*tst1);
}
d[l] = d[l] + f;
e[l] = 0.0;
}
/* Sort eigenvalues and corresponding vectors. */
int j;
double p;
for (i = 0; i < n-1; i++) {
k = i;
p = d[i];
for (j = i+1; j < n; j++) {
if (d[j] < p) {
k = j;
p = d[j];
}
}
if (k != i) {
d[k] = d[i];
d[i] = p;
for (j = 0; j < n; j++) {
p = V[j][i];
V[j][i] = V[j][k];
V[j][k] = p;
}
}
}
} /* QLalgo2 */
void Eigen( int SIZE, double **C, double *diag, double **Q, double *rgtmp)
{
int i, j;
if (rgtmp == NULL) /* was OK in former versions */
printf("cmaes_t:Eigen(): input parameter double *rgtmp must be non-NULL\n");
/* copy C to Q */
if (C != Q) {
for (i=0; i < SIZE; ++i)
for (j = 0; j <= i; ++j)
Q[i][j] = Q[j][i] = C[i][j];
}
Householder2( SIZE, Q, diag, rgtmp);
QLalgo2( SIZE, diag, rgtmp, Q);
}
void cmaes_UpdateEigensystem(CMA *t, int flgforce)
{
int i;
if(flgforce == 0) {
if (t->flgEigensysIsUptodate == 1)
return;
/* return on modulo generation number */
if (t->gen < t->genOfEigensysUpdate + t->updateCmode.modulo)