Commit bdc3da2a authored by Konstantinos Chatzilygeroudis's avatar Konstantinos Chatzilygeroudis
Browse files

Removed old cmaes

parent bd92a748
This code is from: https://www.lri.fr/~hansen/cmaesintro.html
#include <stddef.h>
#include <stdlib.h>
#include <math.h>
#include <float.h> /* DBL_MAX */
#include <stdio.h>
#include "boundary_transformation.h"
static int do_assertions = 1;
static unsigned long _index(boundary_transformation_t* t, unsigned long i);
static void _FatalError(const char* s);
static double default_lower[1];
static double default_upper[1];
void boundary_transformation_init(boundary_transformation_t* t,
double const* lower_bounds,
double const* upper_bounds,
unsigned long len_of_bounds)
{
unsigned i;
double const* ub, *lb;
default_lower[0] = DBL_MAX / -1e2;
default_upper[0] = DBL_MAX / 1e2;
t->lower_bounds = lower_bounds;
t->upper_bounds = upper_bounds;
t->len_of_bounds = len_of_bounds;
if (lower_bounds == NULL && len_of_bounds <= 1) /* convenience default */
t->lower_bounds = default_lower;
if (upper_bounds == NULL && len_of_bounds <= 1)
t->upper_bounds = default_upper;
if (len_of_bounds == 0) {
t->lower_bounds = default_lower;
t->upper_bounds = default_upper;
t->len_of_bounds = 1;
}
if (t->lower_bounds == NULL || t->upper_bounds == NULL)
_FatalError("init: input upper_bounds or lower_bounds was NULL and "
"len_of_bounds > 1");
/* compute boundaries in pre-image space, al and au */
t->al = calloc(t->len_of_bounds, sizeof(double));
t->au = calloc(t->len_of_bounds, sizeof(double));
if (!t->al || !t->au)
_FatalError(" in _init(): could not allocate memory");
lb = t->lower_bounds;
ub = t->upper_bounds;
for (i = 0; i < t->len_of_bounds; ++i) {
if (lb[i] == ub[i])
_FatalError("in _init: lower and upper bounds must be different in all "
"variables");
/* between lb+al and ub-au transformation is the identity */
t->al[i] = fmin((ub[i] - lb[i]) / 2., (1. + fabs(lb[i])) / 20.);
t->au[i] = fmin((ub[i] - lb[i]) / 2., (1. + fabs(ub[i])) / 20.);
}
}
void boundary_transformation_exit(boundary_transformation_t* t)
{
if (t->al)
free(t->al);
if (t->au)
free(t->au);
}
void boundary_transformation(boundary_transformation_t* t, double const* x,
double* y, unsigned long len)
{
double lb, ub, al, au;
unsigned long i;
boundary_transformation_shift_into_feasible_preimage(t, x, y, len);
for (i = 0; i < len; ++i) {
lb = t->lower_bounds[_index(t, i)];
ub = t->upper_bounds[_index(t, i)];
al = t->al[_index(t, i)];
au = t->au[_index(t, i)];
if (y[i] < lb + al)
y[i] = lb + (y[i] - (lb - al)) * (y[i] - (lb - al)) / 4. / al;
else if (y[i] > ub - au)
y[i] = ub - (y[i] - (ub + au)) * (y[i] - (ub + au)) / 4. / au;
}
}
void boundary_transformation_discrete(boundary_transformation_t* t,
double const* x, double* y,
unsigned long len)
{
double lb, ub, al, au;
unsigned long i;
// boundary_transformation_shift_into_feasible_preimage(t, x, y, len);
for (i = 0; i < len; ++i) {
lb = t->lower_bounds[_index(t, i)];
ub = t->upper_bounds[_index(t, i)];
al = t->al[_index(t, i)];
au = t->au[_index(t, i)];
y[i] = x[i];
if (y[i] < lb)
y[i] = lb;
else if (y[i] > ub)
y[i] = ub;
}
}
void boundary_transformation_shift_into_feasible_preimage(
boundary_transformation_t* t, double const* x, double* y,
unsigned long len)
{
double lb, ub, al, au, r, xlow, xup;
unsigned long i;
for (i = 0; i < len; ++i) {
lb = t->lower_bounds[_index(t, i)];
ub = t->upper_bounds[_index(t, i)];
al = t->al[_index(t, i)];
au = t->al[_index(t, i)];
xlow = lb - 2 * al - (ub - lb) / 2.0;
xup = ub + 2 * au + (ub - lb) / 2.0;
r = 2 * (ub - lb + al + au); /* == xup - xlow == period of the transformation */
y[i] = x[i];
if (y[i] < xlow) { /* shift up */
y[i] += r * (1 + (int)((xlow - y[i]) / r));
}
if (y[i] > xup) { /* shift down */
y[i] -= r * (1 + (int)((y[i] - xup) / r));
/* printf(" \n%f\n", fmod(y[i] - ub - au, r)); */
}
if (y[i] < lb - al) /* mirror */
y[i] += 2 * (lb - al - y[i]);
if (y[i] > ub + au)
y[i] -= 2 * (y[i] - ub - au);
if ((y[i] < lb - al - 1e-15) || (y[i] > ub + au + 1e-15)) {
printf("BUG in boundary_transformation_shift_into_feasible_preimage: "
"lb=%f, ub=%f, al=%f au=%f, y=%f\n",
lb, ub, al, au, y[i]);
_FatalError("BUG");
}
}
}
void boundary_transformation_inverse(boundary_transformation_t* t,
double const* x, double* y,
unsigned long len)
{
double lb, ub, al, au;
unsigned long i;
for (i = 0; i < len; ++i) {
lb = t->lower_bounds[_index(t, i)];
ub = t->upper_bounds[_index(t, i)];
al = t->al[_index(t, i)];
au = t->al[_index(t, i)];
y[i] = x[i];
if (y[i] < lb + al)
y[i] = (lb - al) + 2 * sqrt(al * (y[i] - lb));
else if (y[i] > ub - au)
y[i] = (ub + au) - 2 * sqrt(au * (ub - y[i]));
}
if (11 < 3 || do_assertions) {
double* z = calloc(len, sizeof(double));
for (i = 0; i < len; ++i)
z[i] = y[i];
boundary_transformation(t, z, y, len);
for (i = 0; i < len; ++i)
if (fabs(y[i] - x[i]) > 1e-14)
printf(" difference for index %ld should be zero, is %f ", i,
y[i] - x[i]);
for (i = 0; i < len; ++i)
y[i] = z[i];
free(z);
}
}
static unsigned long _index(boundary_transformation_t* t, unsigned long i)
{
return i < t->len_of_bounds ? i : t->len_of_bounds - 1;
}
static void _FatalError(const char* s)
{
printf("Fatal error in boundary_transformation: %s\n", s);
exit(1);
}
#if defined(__cplusplus)
extern "C" {
#endif
/*
* INTERFACE: the relevant functions are
*
* * boundary_transformation_init(this, l_bound, u_bound, len)
* * boundary_transformation_exit(this)
* * boundary_transformation(this, x, y, len)
*
* implements a smooth mapping double *x -> double *y that guarantees
* elements of y to be within specified boundaries. The mapping is piecewise
* either linear or quadratic and can achieve values arbitrarily close to and
* on the boundaries. The middle of the domain l_bound + (u_bound-l_bound) / 2.0
* always maps to itself. Typically, feasible values not close to the boundaries
* are mapped to itself, preventing any numerical subtleties. Specifically,
* al, au > 0 are internally chosen offsets. The mapping
* [l_bound - al, u_bound + au] -> [l_bound, u_bound] is monotonous, bijective
* and invertible and furthermore values in [l_bound + al, u_bound - au] map to
* itself.
*
* The method is robust against very small/large boundary values, say
* -1e99 and/or 1e99, to emulated unbounded variables. In this case values
* between -1e98 and 1e98 are never changed, i.e. mapped to itself.
*
*/
typedef struct
{
double const* lower_bounds; /* array of size len_of_bounds */
double const* upper_bounds; /* array of size len_of_bounds */
unsigned long len_of_bounds; /* in case, last value is recycled */
double* al; /* "add"-on to lower boundary preimage, same length as bounds */
double* au; /* add-on to upper boundary preimage, same length as bounds */
} boundary_transformation_t;
/* set lower and upper bounds, the values lower_bounds[len_of_bounds - 1] and
* upper_bounds[len_of_bounds - 1] are recycled for any element >=
* len_of_bounds.
* If len_of_bounds == 0, no bounds are assumed. If len_of_bounds == 1, the
* zero pointer is allowed for lower_bounds or upper_bounds and indicates no
* respective bounds. "no bounds" is "emulated" using the very small/large value
* of DBL_MAX / -1e2 and DBL_MAX / 1e2, respectively. */
void boundary_transformation_init(boundary_transformation_t*,
double const* lower_bounds,
double const* upper_bounds,
unsigned long len_of_bounds);
/* release memory */
void boundary_transformation_exit(boundary_transformation_t*);
/* on return, y is guaranteed to have all values within the boundaries.
* The caller inputs x and is responsible for having allocated y in that
* y[len-1] = x[len-1] is a valid operation. x==y is valid input, but
* will fail together with cmaes when x is an element of the population
* returned by cmaes_SamplePopulation (these elements are of type
* double const * for a reason).
* */
void boundary_transformation(boundary_transformation_t*, double const* x,
double* y,
unsigned long len); /* new value into y */
void boundary_transformation_discrete(boundary_transformation_t*,
double const* x, double* y,
unsigned long len); /* new value into y */
/* after
* boundary_transformation(b,x,y,l) ;
* the two consecutive calls
* boundary_transformation_inverse(b,y,x,l) ; boundary_transformation(b,x,y,l)
* ;
* have no effect on y anymore (but they might change x!).
* */
void boundary_transformation_inverse(boundary_transformation_t* t,
double const* y, double* x,
unsigned long len); /* new value into x */
/* used by function boundary_transformation. After applying the shift,
* boundary_transformation_shift_into_feasible_preimage(b,x,x,l)
* the two consecutive calls
* boundary_transformation(b,x,y,l) ; boundary_transformation_inverse(b,y,x,l)
* ;
* have no effect on x anymore */
void boundary_transformation_shift_into_feasible_preimage(
boundary_transformation_t* t, double const* x, double* x_shifted,
unsigned long len); /* new value into x_shifted */
#if defined(__cplusplus)
}
#endif
This diff is collapsed.
/* --------------------------------------------------------- */
/* --- File: cmaes.h ----------- Author: Nikolaus Hansen --- */
/* ---------------------- last modified: IX 2010 --- */
/* --------------------------------- by: Nikolaus Hansen --- */
/* --------------------------------------------------------- */
/*
CMA-ES for non-linear function minimization.
Copyright (C) 1996, 2003-2010 Nikolaus Hansen.
e-mail: nikolaus.hansen (you know what) inria.fr
License: see file cmaes.c
*/
#ifndef NH_cmaes_h /* only include ones */
#define NH_cmaes_h
#include <time.h>
typedef struct
/* random_t
* sets up a pseudo random number generator instance
*/
{
/* Variables for Uniform() */
long int startseed;
long int aktseed;
long int aktrand;
long int* rgrand;
/* Variables for Gauss() */
short flgstored;
double hold;
} random_t;
typedef struct
/* timings_t
* time measurement, used to time eigendecomposition
*/
{
/* for outside use */
double totaltime; /* zeroed by calling re-calling timings_start */
double totaltotaltime;
double tictoctime;
double lasttictoctime;
/* local fields */
clock_t lastclock;
time_t lasttime;
clock_t ticclock;
time_t tictime;
short istic;
short isstarted;
double lastdiff;
double tictoczwischensumme;
} timings_t;
typedef struct
/* readpara_t
* collects all parameters, in particular those that are read from
* a file before to start. This should split in future?
*/
{
char*
filename; /* keep record of the file that was taken to read parameters */
/* input parameters */
int N; /* problem dimension, must stay constant, should be unsigned or long?
*/
unsigned int seed;
double* xstart;
double* typicalX;
int typicalXcase;
double* rgInitialStds;
double* rgDiffMinChange;
/* termination parameters */
double stopMaxFunEvals;
double facmaxeval;
double stopMaxIter;
struct
{
int flg;
double val;
} stStopFitness;
double stopTolFun;
double stopTolFunHist;
double stopTolX;
double stopTolUpXFactor;
/* internal evolution strategy parameters */
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 */
double diagonalCov; /* number of initial iterations */
struct
{
int flgalways;
double modulo;
double maxtime;
} updateCmode;
double facupdateCmode;
/* supplementary variables */
char* weigkey;
char resumefile[99];
const char** rgsformat;
void** rgpadr;
const char** rgskeyar;
double*** rgp2adr;
int n1para, n1outpara;
int n2para;
} readpara_t;
typedef struct
/* cmaes_t
* CMA-ES "object"
*/
{
const char* version;
/* char *signalsFilename; */
readpara_t sp;
random_t rand; /* random number generator */
double sigma; /* step size */
double* rgxmean; /* mean x vector, "parent" */
double* rgxbestever;
double** rgrgx; /* range of x-vectors, lambda offspring */
int* index; /* sorting index of sample pop. */
double* arFuncValueHist;
short flgIniphase; /* not really in use anymore */
short flgStop;
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 */
double* rgFuncValue;
double* publicFitness; /* returned by cmaes_init() */
double gen; /* Generation number */
double countevals;
double state; /* 1 == sampled, 2 == not in use anymore, 3 == updated */
double maxdiagC; /* repeatedly used for output */
double mindiagC;
double maxEW;
double minEW;
char sOutString[330]; /* 4x80 */
short flgEigensysIsUptodate;
short flgCheckEigen; /* control via cmaes_signals.par */
double genOfEigensysUpdate;
timings_t eigenTimings;
double dMaxSignifKond;
double dLastMinEWgroesserNull;
short flgresumedone;
time_t printtime;
time_t writetime; /* ideally should keep track for each output file */
time_t firstwritetime;
time_t firstprinttime;
} cmaes_t;
#endif
#ifndef CMAES_INTERFACE_H
#define CMAES_INTERFACE_H
#ifdef __cplusplus
extern "C" {
#endif
/* --------------------------------------------------------- */
/* --- File: cmaes_interface.h - Author: Nikolaus Hansen --- */
/* ---------------------- last modified: IV 2007 --- */
/* --------------------------------- by: Nikolaus Hansen --- */
/* --------------------------------------------------------- */
/*
CMA-ES for non-linear function minimization.
Copyright (C) 1996, 2003, 2007 Nikolaus Hansen.
e-mail: hansen AT lri.fr
License: see file cmaes.c
*/
#include "cmaes.h"
/* --------------------------------------------------------- */
/* ------------------ Interface ---------------------------- */
/* --------------------------------------------------------- */
/* --- initialization, constructors, destructors --- */
double* cmaes_init(cmaes_t*, int dimension, double* xstart, double* stddev,
long seed, int lambda, const char* input_parameter_filename);
void cmaes_resume_distribution(cmaes_t* evo_ptr, char* filename);
void cmaes_exit(cmaes_t*);
/* --- core functions --- */
double* const* cmaes_SamplePopulation(cmaes_t*);
double* cmaes_UpdateDistribution(cmaes_t*, const double* rgFitnessValues);
const char* cmaes_TestForTermination(cmaes_t*);
/* --- additional functions --- */
double* const* cmaes_ReSampleSingle(cmaes_t* t, int index);
double const* cmaes_ReSampleSingle_old(cmaes_t*, double* rgx);
double* cmaes_SampleSingleInto(cmaes_t* t, double* rgx);
void cmaes_UpdateEigensystem(cmaes_t*, int flgforce);
/* --- getter functions --- */
double cmaes_Get(cmaes_t*, char const* keyword);
const double* cmaes_GetPtr(cmaes_t*,
char const* keyword); /* e.g. "xbestever" */
double* cmaes_GetNew(cmaes_t* t,
char const* keyword); /* user is responsible to free */
double* cmaes_GetInto(
cmaes_t* t, char const* keyword,
double* mem); /* allocs if mem==NULL, user is responsible to free */
/* --- online control and output --- */
void cmaes_ReadSignals(cmaes_t*, char const* filename);
void cmaes_WriteToFile(cmaes_t*, const char* szKeyWord,
const char* output_filename);
char* cmaes_SayHello(cmaes_t*);
/* --- misc --- */
double* cmaes_NewDouble(int n); /* user is responsible to free */
void cmaes_FATAL(char const* s1, char const* s2, char const* s3,
char const* s4);
#ifdef __cplusplus
} /* Extern c */
#endif
#endif
......@@ -70,6 +70,6 @@ def create(bld):
source='/combinations/combinations_' + str(i) + '.cpp',
includes='. .. ../../',
target='/combinations/combinations_' + str(i),
uselib='BOOST EIGEN TBB SFERES',
uselib='BOOST EIGEN TBB SFERES LIBCMAES',
use='limbo')
i = i + 1
......@@ -8,8 +8,6 @@ def build(bld):
bld.recurse('benchmarks')
bld.stlib(source=' \
cmaes/cmaes.c \
cmaes/boundary_transformation.c \
ehvi/ehvi_calculations.cc \
ehvi/ehvi_montecarlo.cc \
ehvi/ehvi_sliceupdate.cc \
......@@ -21,8 +19,6 @@ def build(bld):
def build_extensive_tests(bld):
bld.stlib(source=' \
cmaes/cmaes.c \
cmaes/boundary_transformation.c \
ehvi/ehvi_calculations.cc \
ehvi/ehvi_montecarlo.cc \
ehvi/ehvi_sliceupdate.cc \
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment