From 905cbb7a970f54ebd0ddd49e0ed09e7daad4fb32 Mon Sep 17 00:00:00 2001
From: Erik Strand <erik.strand@cba.mit.edu>
Date: Thu, 16 Apr 2020 10:19:40 -0400
Subject: [PATCH] Copy in the reference CMA-ES implementation
---
external/CMakeLists.txt | 5 +
external/cma-es/CMakeLists.txt | 12 +
external/cma-es/cmaes.cpp | 3137 +++++++++++++++++++++++++++++
external/cma-es/cmaes.h | 176 ++
external/cma-es/cmaes_interface.h | 68 +
5 files changed, 3398 insertions(+)
create mode 100644 external/cma-es/CMakeLists.txt
create mode 100644 external/cma-es/cmaes.cpp
create mode 100644 external/cma-es/cmaes.h
create mode 100644 external/cma-es/cmaes_interface.h
diff --git a/external/CMakeLists.txt b/external/CMakeLists.txt
index f1f7d5e..651bf92 100644
--- a/external/CMakeLists.txt
+++ b/external/CMakeLists.txt
@@ -7,6 +7,11 @@ target_include_directories(catch2 INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/catch2)
add_library(clara INTERFACE)
target_include_directories(clara INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/clara)
+# CMA-ES (c-cmaes)
+# From https://github.com/CMA-ES/c-cmaes. Not using it as a submodule since their repo isn't
+# organized in the most extensible way and I'm compiling it as C++ anyway.
+add_subdirectory(cma-es)
+
# nlohmann JSON
add_library(json INTERFACE)
target_include_directories(json INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/json)
diff --git a/external/cma-es/CMakeLists.txt b/external/cma-es/CMakeLists.txt
new file mode 100644
index 0000000..bd63509
--- /dev/null
+++ b/external/cma-es/CMakeLists.txt
@@ -0,0 +1,12 @@
+add_library(cma-es
+ cmaes.cpp
+ cmaes.h
+ cmaes_interface.h
+)
+target_include_directories(cma-es PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
+target_link_libraries(cma-es shared_settings)
+# This isn't working, since shared_settings enables this and it seems to take precedence (regardless
+# of the order I put things here).
+target_compile_options(cma-es PRIVATE -Wno-zero-as-null-pointer-constant)
+# At some point I added this, doesn't seem to do anything now but I'll leave it here for a while.
+#target_compile_options(cma-es PRIVATE -Wno-stringop-overflow)
diff --git a/external/cma-es/cmaes.cpp b/external/cma-es/cmaes.cpp
new file mode 100644
index 0000000..a38e1d4
--- /dev/null
+++ b/external/cma-es/cmaes.cpp
@@ -0,0 +1,3137 @@
+/* --------------------------------------------------------- */
+/* --- File: cmaes.c -------- Author: Nikolaus Hansen --- */
+/* --------------------------------------------------------- */
+/*
+ CMA-ES for non-linear function minimization.
+
+ Copyright 1996, 2003, 2007, 2013 Nikolaus Hansen
+ e-mail: hansen .AT. lri.fr
+
+ SOURCE:
+ https://github.com/cma-es/c-cma-es
+ https://github.com/cma-es/c-cma-es/blob/master/src/cmaes.c
+
+ LICENSE: this library is free/open software and may be used
+ either under the
+
+ Apache License 2.0
+
+ or under the
+
+ GNU Lesser General Public License 2.1 or later
+
+ whichever suits best.
+
+ See also the LICENSE file
+ https://github.com/cma-es/c-cma-es/blob/master/LICENSE
+*/
+
+/* --- Changes : ---
+ 03/03/21: argument const double *rgFunVal of
+ cmaes_ReestimateDistribution() was treated incorrectly.
+ 03/03/29: restart via cmaes_resume_distribution() implemented.
+ 03/03/30: Always max std dev / largest axis is printed first.
+ 03/08/30: Damping is adjusted for large mueff.
+ 03/10/30: Damping is adjusted for large mueff always.
+ 04/04/22: Cumulation time and damping for step size adjusted.
+ No iniphase but conditional update of pc.
+ 05/03/15: in ccov-setting mucov replaced by mueff.
+ 05/10/05: revise comment on resampling in example.c
+ 05/10/13: output of "coorstddev" changed from sigma * C[i][i]
+ to correct sigma * sqrt(C[i][i]).
+ 05/11/09: Numerical problems are not anymore handled by increasing
+ sigma, but lead to satisfy a stopping criterion in
+ cmaes_Test().
+ 05/11/09: Update of eigensystem and test for numerical problems
+ moved right before sampling.
+ 06/02/24: Non-ansi array definitions replaced (thanks to Marc
+ Toussaint).
+ 06/02/25: Overflow in time measurement for runs longer than
+ 2100 seconds. This could lead to stalling the
+ covariance matrix update for long periods.
+ Time measurement completely rewritten.
+ 06/02/26: Included population size lambda as parameter to
+ cmaes_init (thanks to MT).
+ 06/02/26: Allow no initial reading/writing of parameters via
+ "non" and "writeonly" keywords for input parameter
+ filename in cmaes_init.
+ 06/02/27: Optimized code regarding time spent in updating the
+ covariance matrix in function Adapt_C2().
+ 07/08/03: clean up and implementation of an exhaustive test
+ of the eigendecomposition (via #ifdef for now)
+ 07/08/04: writing of output improved
+ 07/08/xx: termination criteria revised and more added,
+ damp replaced by damps=damp*cs, documentation improved.
+ Interface significantly changed, evaluateSample function
+ and therefore the function pointer argument removed.
+ Renaming of functions in accordance with Java code.
+ Clean up of parameter names, mainly in accordance with
+ Matlab conventions. Most termination criteria can be
+ changed online now. Many more small changes, but not in
+ the core procedure.
+ 07/10/29: ReSampleSingle() got a better interface. ReSampleSingle()
+ is now ReSampleSingle_old only for backward
+ compatibility. Also fixed incorrect documentation. The new
+ function SampleSingleInto() has an interface similar to
+ the old ReSampleSingle(), but is not really necessary.
+ 07/11/20: bug: stopMaxIter did not translate into the correct default
+ value but into -1 as default. This lead to a too large
+ damps and the termination test became true from the first
+ iteration. (Thanks to Michael Calonder)
+ 07/11/20: new default stopTolFunHist = 1e-13; (instead of zero)
+ 08/09/26: initial diagonal covariance matrix in code, but not
+ yet in interface
+ 08/09/27: diagonalCovarianceMatrix option in initials.par provided
+ 08/10/17: uncertainty handling implemented in example3.c.
+ PerturbSolutionInto() provides the optional small
+ perturbations before reevaluation.
+ 10/10/16: TestForTermination changed such that diagonalCovarianceMatrix
+ option now yields linear time behavior
+ 12/05/28: random seed > 2e9 prohibited to avoid an infinite loop on 32bit systems
+ 12/10/21: input parameter file values "no", "none" now work as "non".
+ 12/10/xx: tentative implementation of cmaes_Optimize
+ 12/10/xx: some small changes with char * mainly to prevent warnings in C++
+ 12/10/xx: added some string convenience functions isNoneStr, new_string, assign_string
+ 13/01/03: rename files example?, initials.par, signals.par
+ 14/04/29: removed bug, au = t->al[...], from the (new) boundary handling
+ code (thanks to Emmanuel Benazera for the hint)
+ 14/06/16: test of un-initialized version number removed (thanks to Paul Zimmermann)
+ 14/10/18: splitted cmaes_init() into cmaes_init_para() and cmaes_init_final(),
+ such that parameters can be set also comparatively safely within the
+ code, namely after calling cmaes_init_para(). The order of parameters
+<<<<<<< HEAD
+ of readpara_init() has changed to be the same as for cmaes_init().
+ 14/11/25 fix warnings from Microsoft C compiler (sherm)
+=======
+ of readpara_init() has changed to be the same as for cmaes_init().
+ 14/11/26: renamed exported symbols so they begin with a cmaes_prefix.
+>>>>>>> bd5aeb2a21fbc437f1b322610022f7cb1ba6a9db
+
+ Wish List
+ o make signals_filename part of cmaes_t using assign_string()
+
+ o as writing time is measure for all files at once, the display
+ cannot be independently written to a file via signals.par, while
+ this would be desirable.
+
+ o clean up sorting of eigenvalues and vectors which is done repeatedly.
+
+ o either use cmaes_Get() in cmaes_WriteToFilePtr(): revise the
+ cmaes_write that all keywords available with get and getptr are
+ recognized. Also revise the keywords, keeping backward
+ compatibility. (not only) for this it would be useful to find a
+ way how cmaes_Get() signals an unrecognized keyword. For GetPtr
+ it can return NULL.
+
+ o or break cmaes_Get() into single getter functions, being a nicer
+ interface, and compile instead of runtime error, and faster. For
+ file signals.par it does not help.
+
+ o writing data depending on timing in a smarter way, e.g. using 10%
+ of all time. First find out whether clock() is useful for measuring
+ disc writing time and then cmaes_timings_t class can be utilized.
+ For very large dimension the default of 1 seconds waiting might
+ be too small.
+
+ o allow modification of best solution depending on delivered f(xmean)
+
+ o re-write input and output procedures
+*/
+
+/* Prevent Microsoft compiler from complaining that common library functions
+ like strncpy(), ctime(), sprintf(), fopen(), fscanf(), etc. "may be unsafe".
+ This must come before the first #include.
+*/
+#ifdef _MSC_VER
+#define _CRT_SECURE_NO_WARNINGS
+#endif
+
+#include <math.h> /* sqrt() */
+#include <stddef.h> /* size_t */
+#include <stdlib.h> /* NULL, free */
+#include <string.h> /* strlen() */
+#include <stdio.h> /* sprintf(), NULL? */
+#include "cmaes_interface.h" /* <time.h> via cmaes.h */
+
+/* --------------------------------------------------------- */
+/* ------------------- Declarations ------------------------ */
+/* --------------------------------------------------------- */
+
+/* ------------------- External Visibly -------------------- */
+
+/* see cmaes_interface.h for those, not listed here */
+
+long cmaes_random_init(cmaes_random_t *, long unsigned seed /* 0==clock */);
+void cmaes_random_exit(cmaes_random_t *);
+double cmaes_random_Gauss(cmaes_random_t *); /* (0,1)-normally distributed */
+double cmaes_random_Uniform(cmaes_random_t *);
+long cmaes_random_Start(cmaes_random_t *, long unsigned seed /* 0==1 */);
+
+void cmaes_timings_init(cmaes_timings_t *timing);
+void cmaes_timings_start(cmaes_timings_t *timing); /* fields totaltime and tictoctime */
+double cmaes_timings_update(cmaes_timings_t *timing);
+void cmaes_timings_tic(cmaes_timings_t *timing);
+double cmaes_timings_toc(cmaes_timings_t *timing);
+
+void cmaes_readpara_init (cmaes_readpara_t *, int dim, const double * xstart,
+ const double * sigma, int seed, int lambda,
+ const char * filename);
+void cmaes_readpara_exit(cmaes_readpara_t *);
+void cmaes_readpara_ReadFromFile(cmaes_readpara_t *, const char *szFileName);
+void cmaes_readpara_SupplementDefaults(cmaes_readpara_t *);
+void cmaes_readpara_SetWeights(cmaes_readpara_t *, const char * mode);
+void cmaes_readpara_WriteToFile(cmaes_readpara_t *, const char *filenamedest);
+
+const double * cmaes_Optimize( cmaes_t *, double(*pFun)(double const *, int dim),
+ long iterations);
+double const * cmaes_SetMean(cmaes_t *, const double *xmean);
+double * cmaes_PerturbSolutionInto(cmaes_t *t, double *xout,
+ double const *xin, double eps);
+void cmaes_WriteToFile(cmaes_t *, const char *key, const char *name);
+void cmaes_WriteToFileAW(cmaes_t *t, const char *key, const char *name,
+ const char * append);
+void cmaes_WriteToFilePtr(cmaes_t *, const char *key, FILE *fp);
+void cmaes_ReadFromFilePtr(cmaes_t *, FILE *fp);
+void cmaes_FATAL(char const *s1, char const *s2,
+ char const *s3, char const *s4);
+
+
+/* ------------------- Locally visibly ----------------------- */
+
+static char * getTimeStr(void);
+static void TestMinStdDevs( cmaes_t *);
+/* static void WriteMaxErrorInfo( cmaes_t *); */
+
+static void Eigen( int N, double **C, double *diag, double **Q,
+ double *rgtmp);
+static int Check_Eigen( int N, double **C, double *diag, double **Q);
+static void QLalgo2 (int n, double *d, double *e, double **V);
+static void Householder2(int n, double **V, double *d, double *e);
+static void Adapt_C2(cmaes_t *t, int hsig);
+
+static void FATAL(char const *sz1, char const *s2,
+ char const *s3, char const *s4);
+static void ERRORMESSAGE(char const *sz1, char const *s2,
+ char const *s3, char const *s4);
+static int isNoneStr(const char * filename);
+static void Sorted_index( const double *rgFunVal, int *index, int n);
+static int SignOfDiff( const void *d1, const void * d2);
+static double douSquare(double);
+static double rgdouMax( const double *rgd, int len);
+static double rgdouMin( const double *rgd, int len);
+static double douMax( double d1, double d2);
+static double douMin( double d1, double d2);
+static int intMin( int i, int j);
+static int MaxIdx( const double *rgd, int len);
+static int MinIdx( const double *rgd, int len);
+static double myhypot(double a, double b);
+static double * new_double( int n);
+static void * new_void( int n, size_t size);
+static char * new_string( const char *);
+static void assign_string( char **, const char*);
+
+static const char * c_cmaes_version = "3.20.00.beta";
+
+/* --------------------------------------------------------- */
+/* ---------------- Functions: cmaes_t --------------------- */
+/* --------------------------------------------------------- */
+
+static char *
+getTimeStr(void) {
+ time_t tm = time(NULL);
+ static char s[33];
+
+ /* get time */
+ strncpy(s, ctime(&tm), 24); /* TODO: hopefully we read something useful */
+ s[24] = '\0'; /* cut the \n */
+ return s;
+}
+
+char *
+cmaes_SayHello(cmaes_t *t)
+{
+ /* write initial message */
+ sprintf(t->sOutString,
+ "(%d,%d)-CMA-ES(mu_eff=%.1f), Ver=\"%s\", dimension=%d, diagonalIterations=%ld, randomSeed=%d (%s)",
+ t->sp.mu, t->sp.lambda, t->sp.mueff, t->version, t->sp.N, (long)t->sp.diagonalCov,
+ t->sp.seed, getTimeStr());
+
+ return t->sOutString;
+}
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+void
+cmaes_init_para(cmaes_t *t, /* "this" */
+ int dimension,
+ double *inxstart,
+ double *inrgstddev, /* initial stds */
+ long int inseed,
+ int lambda,
+ const char *input_parameter_filename)
+{
+ t->version = c_cmaes_version;
+ cmaes_readpara_init(&t->sp, dimension, inxstart, inrgstddev, inseed,
+ lambda, input_parameter_filename);
+}
+
+double *
+cmaes_init_final(cmaes_t *t /* "this" */)
+/*
+ * */
+{
+ int i, j, N;
+ double dtest, trace;
+
+ if (t->version == NULL) {
+ ERRORMESSAGE("cmaes_init_final called (probably) without calling cmaes_init_para first",
+ ", which will likely lead to unexpected results",0,0);
+ printf("Error: cmaes_init_final called (probably) without calling cmaes_init_para first\n");
+ }
+ if (strcmp(c_cmaes_version, t->version) != 0) {
+ ERRORMESSAGE("cmaes_init_final called twice, which will lead to a memory leak",
+ "; use cmaes_exit() first",0,0);
+ printf("Error: cmaes_init_final called twice, which will lead to a memory leak; use cmaes_exit first\n");
+ }
+ /* assign_string(&t->signalsFilename, "cmaes_signals.par"); */
+
+ if (!t->sp.flgsupplemented) {
+ cmaes_readpara_SupplementDefaults(&t->sp);
+ if (!isNoneStr(t->sp.filename)) /* TODO: should this be done in readpara_SupplementDefaults? */
+ cmaes_readpara_WriteToFile(&t->sp, "actparcmaes.par");
+ }
+
+ t->sp.seed = cmaes_random_init( &t->rand, (long unsigned int) t->sp.seed);
+
+ N = t->sp.N; /* for convenience */
+
+ /* initialization */
+ for (i = 0, trace = 0.; i < N; ++i)
+ trace += t->sp.rgInitialStds[i]*t->sp.rgInitialStds[i];
+ t->sigma = sqrt(trace/N); /* t->sp.mueff/(0.2*t->sp.mueff+sqrt(N)) * sqrt(trace/N); */
+
+ t->chiN = sqrt((double) N) * (1. - 1./(4.*N) + 1./(21.*N*N));
+ t->flgEigensysIsUptodate = 1;
+ t->flgCheckEigen = 0;
+ t->genOfEigensysUpdate = 0;
+ cmaes_timings_init(&t->eigenTimings);
+ t->flgIniphase = 0; /* do not use iniphase, hsig does the job now */
+ t->flgresumedone = 0;
+ t->flgStop = 0;
+
+ for (dtest = 1.; dtest && dtest < 1.1 * dtest; dtest *= 2.)
+ if (dtest == dtest + 1.)
+ break;
+ t->dMaxSignifKond = dtest / 1000.; /* not sure whether this is really save, 100 does not work well enough */
+
+ t->gen = 0;
+ t->countevals = 0;
+ t->state = 0;
+ t->dLastMinEWgroesserNull = 1.0;
+ t->printtime = t->writetime = t->firstwritetime = t->firstprinttime = 0;
+
+ t->rgpc = new_double(N);
+ t->rgps = new_double(N);
+ t->rgdTmp = new_double(N+1);
+ t->rgBDz = new_double(N);
+ t->rgxmean = new_double(N+2); t->rgxmean[0] = N; ++t->rgxmean;
+ t->rgxold = new_double(N+2); t->rgxold[0] = N; ++t->rgxold;
+ t->rgxbestever = new_double(N+3); t->rgxbestever[0] = N; ++t->rgxbestever;
+ t->rgout = new_double(N+2); t->rgout[0] = N; ++t->rgout;
+ t->rgD = new_double(N);
+ t->C = (double**)new_void(N, sizeof(double*));
+ t->B = (double**)new_void(N, sizeof(double*));
+ t->publicFitness = new_double(t->sp.lambda);
+ t->rgFuncValue = new_double(t->sp.lambda+1);
+ t->rgFuncValue[0]=t->sp.lambda; ++t->rgFuncValue;
+ t->arFuncValueHist = new_double(10+(int)ceil(3.*10.*N/t->sp.lambda)+1);
+ t->arFuncValueHist[0] = (double)(10+(int)ceil(3.*10.*N/t->sp.lambda));
+ t->arFuncValueHist++;
+
+ for (i = 0; i < N; ++i) {
+ t->C[i] = new_double(i+1);
+ t->B[i] = new_double(N);
+ }
+ t->index = (int *) new_void(t->sp.lambda, sizeof(int));
+ for (i = 0; i < t->sp.lambda; ++i)
+ t->index[i] = i; /* should not be necessary */
+ t->rgrgx = (double **)new_void(t->sp.lambda, sizeof(double*));
+ for (i = 0; i < t->sp.lambda; ++i) {
+ t->rgrgx[i] = new_double(N+2);
+ t->rgrgx[i][0] = N;
+ t->rgrgx[i]++;
+ }
+
+ /* Initialize newed space */
+
+ for (i = 0; i < N; ++i)
+ for (j = 0; j < i; ++j)
+ t->C[i][j] = t->B[i][j] = t->B[j][i] = 0.;
+
+ for (i = 0; i < N; ++i)
+ {
+ t->B[i][i] = 1.;
+ t->C[i][i] = t->rgD[i] = t->sp.rgInitialStds[i] * sqrt(N / trace);
+ t->C[i][i] *= t->C[i][i];
+ t->rgpc[i] = t->rgps[i] = 0.;
+ }
+
+ t->minEW = rgdouMin(t->rgD, N); t->minEW = t->minEW * t->minEW;
+ t->maxEW = rgdouMax(t->rgD, N); t->maxEW = t->maxEW * t->maxEW;
+
+ t->maxdiagC=t->C[0][0]; for(i=1;i<N;++i) if(t->maxdiagC<t->C[i][i]) t->maxdiagC=t->C[i][i];
+ t->mindiagC=t->C[0][0]; for(i=1;i<N;++i) if(t->mindiagC>t->C[i][i]) t->mindiagC=t->C[i][i];
+
+ /* set xmean */
+ for (i = 0; i < N; ++i)
+ t->rgxmean[i] = t->rgxold[i] = t->sp.xstart[i];
+ /* use in case xstart as typicalX */
+ if (t->sp.typicalXcase)
+ for (i = 0; i < N; ++i)
+ t->rgxmean[i] += t->sigma * t->rgD[i] * cmaes_random_Gauss(&t->rand);
+
+ if (strcmp(t->sp.resumefile, "_no_") != 0)
+ cmaes_resume_distribution(t, t->sp.resumefile);
+
+ return (t->publicFitness);
+
+} /* cmaes_init_final() */
+
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+double *
+cmaes_init(cmaes_t *t, /* "this" */
+ int dimension,
+ double *inxstart,
+ double *inrgstddev, /* initial stds */
+ long int inseed,
+ int lambda,
+ const char *input_parameter_filename)
+{
+ cmaes_init_para(t, dimension, inxstart, inrgstddev, inseed,
+ lambda, input_parameter_filename);
+ return cmaes_init_final(t);
+}
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+
+#ifdef __GNUC__
+ #pragma GCC diagnostic push
+ #pragma GCC diagnostic ignored "-Wunused-result"
+#endif
+
+void
+cmaes_resume_distribution(cmaes_t *t, char *filename)
+{
+ int i, j, res, n;
+ double d;
+ FILE *fp = fopen( filename, "r");
+ if(fp == NULL) {
+ ERRORMESSAGE("cmaes_resume_distribution(): could not open '",
+ filename, "'",0);
+ return;
+ }
+ /* count number of "resume" entries */
+ i = 0; res = 0;
+ while (1) {
+ if ((res = fscanf(fp, " resume %lg", &d)) == EOF)
+ break;
+ else if (res==0)
+ fscanf(fp, " %*s");
+ else if(res > 0)
+ i += 1;
+ }
+
+ /* go to last "resume" entry */
+ n = i; i = 0; res = 0; rewind(fp);
+ while (i<n) {
+ if ((res = fscanf(fp, " resume %lg", &d)) == EOF)
+ FATAL("cmaes_resume_distribution(): Unexpected error, bug",0,0,0);
+ else if (res==0)
+ fscanf(fp, " %*s");
+ else if(res > 0)
+ ++i;
+ }
+ if (d != t->sp.N)
+ FATAL("cmaes_resume_distribution(): Dimension numbers do not match",0,0,0);
+
+ /* find next "xmean" entry */
+ while (1) {
+ if ((res = fscanf(fp, " xmean %lg", &d)) == EOF)
+ FATAL("cmaes_resume_distribution(): 'xmean' not found",0,0,0);
+ else if (res==0)
+ fscanf(fp, " %*s");
+ else if(res > 0)
+ break;
+ }
+
+ /* read xmean */
+ t->rgxmean[0] = d; res = 1;
+ for(i = 1; i < t->sp.N; ++i)
+ res += fscanf(fp, " %lg", &t->rgxmean[i]);
+ if (res != t->sp.N)
+ FATAL("cmaes_resume_distribution(): xmean: dimensions differ",0,0,0);
+
+ /* find next "path for sigma" entry */
+ while (1) {
+ if ((res = fscanf(fp, " path for sigma %lg", &d)) == EOF)
+ FATAL("cmaes_resume_distribution(): 'path for sigma' not found",0,0,0);
+ else if (res==0)
+ fscanf(fp, " %*s");
+ else if(res > 0)
+ break;
+ }
+
+ /* read ps */
+ t->rgps[0] = d; res = 1;
+ for(i = 1; i < t->sp.N; ++i)
+ res += fscanf(fp, " %lg", &t->rgps[i]);
+ if (res != t->sp.N)
+ FATAL("cmaes_resume_distribution(): ps: dimensions differ",0,0,0);
+
+ /* find next "path for C" entry */
+ while (1) {
+ if ((res = fscanf(fp, " path for C %lg", &d)) == EOF)
+ FATAL("cmaes_resume_distribution(): 'path for C' not found",0,0,0);
+ else if (res==0)
+ fscanf(fp, " %*s");
+ else if(res > 0)
+ break;
+ }
+ /* read pc */
+ t->rgpc[0] = d; res = 1;
+ for(i = 1; i < t->sp.N; ++i)
+ res += fscanf(fp, " %lg", &t->rgpc[i]);
+ if (res != t->sp.N)
+ FATAL("cmaes_resume_distribution(): pc: dimensions differ",0,0,0);
+
+ /* find next "sigma" entry */
+ while (1) {
+ if ((res = fscanf(fp, " sigma %lg", &d)) == EOF)
+ FATAL("cmaes_resume_distribution(): 'sigma' not found",0,0,0);
+ else if (res==0)
+ fscanf(fp, " %*s");
+ else if(res > 0)
+ break;
+ }
+ t->sigma = d;
+
+ /* find next entry "covariance matrix" */
+ while (1) {
+ if ((res = fscanf(fp, " covariance matrix %lg", &d)) == EOF)
+ FATAL("cmaes_resume_distribution(): 'covariance matrix' not found",0,0,0);
+ else if (res==0)
+ fscanf(fp, " %*s");
+ else if(res > 0)
+ break;
+ }
+ /* read C */
+ t->C[0][0] = d; res = 1;
+ for (i = 1; i < t->sp.N; ++i)
+ for (j = 0; j <= i; ++j)
+ res += fscanf(fp, " %lg", &t->C[i][j]);
+ if (res != (t->sp.N*t->sp.N+t->sp.N)/2)
+ FATAL("cmaes_resume_distribution(): C: dimensions differ",0,0,0);
+
+ fclose(fp);
+
+ t->flgIniphase = 0;
+ t->flgEigensysIsUptodate = 0;
+ t->flgresumedone = 1;
+ cmaes_UpdateEigensystem(t, 1);
+
+} /* cmaes_resume_distribution() */
+#ifdef __GNUC__
+ #pragma GCC diagnostic pop
+#endif
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+
+void
+cmaes_exit(cmaes_t *t)
+{
+ int i, N = t->sp.N;
+ t->version = NULL;
+ /* free(t->signals_filename) */
+ t->state = -1; /* not really useful at the moment */
+ free( t->rgpc);
+ free( t->rgps);
+ free( t->rgdTmp);
+ free( t->rgBDz);
+ free( --t->rgxmean);
+ free( --t->rgxold);
+ free( --t->rgxbestever);
+ free( --t->rgout);
+ free( t->rgD);
+ for (i = 0; i < N; ++i) {
+ free( t->C[i]);
+ free( t->B[i]);
+ }
+ for (i = 0; i < t->sp.lambda; ++i)
+ free( --t->rgrgx[i]);
+ free( t->rgrgx);
+ free( t->C);
+ free( t->B);
+ free( t->index);
+ free( t->publicFitness);
+ free( --t->rgFuncValue);
+ free( --t->arFuncValueHist);
+ cmaes_random_exit (&t->rand);
+ cmaes_readpara_exit (&t->sp);
+} /* cmaes_exit() */
+
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+double const *
+cmaes_SetMean(cmaes_t *t, const double *xmean)
+/*
+ * Distribution mean could be changed before SamplePopulation().
+ * This might lead to unexpected behaviour if done repeatedly.
+ */
+{
+ int i, N=t->sp.N;
+
+ if (t->state >= 1 && t->state < 3)
+ FATAL("cmaes_SetMean: mean cannot be set inbetween the calls of ",
+ "SamplePopulation and UpdateDistribution",0,0);
+
+ if (xmean != NULL && xmean != t->rgxmean)
+ for(i = 0; i < N; ++i)
+ t->rgxmean[i] = xmean[i];
+ else
+ xmean = t->rgxmean;
+
+ return xmean;
+}
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+double * const *
+cmaes_SamplePopulation(cmaes_t *t)
+{
+ int iNk, i, j, N=t->sp.N;
+ int flgdiag = ((t->sp.diagonalCov == 1) || (t->sp.diagonalCov >= t->gen));
+ double sum;
+ double const *xmean = t->rgxmean;
+
+ /* cmaes_SetMean(t, xmean); * xmean could be changed at this point */
+
+ /* calculate eigensystem */
+ if (!t->flgEigensysIsUptodate) {
+ if (!flgdiag)
+ cmaes_UpdateEigensystem(t, 0);
+ else {
+ for (i = 0; i < N; ++i)
+ t->rgD[i] = sqrt(t->C[i][i]);
+ t->minEW = douSquare(rgdouMin(t->rgD, N));
+ t->maxEW = douSquare(rgdouMax(t->rgD, N));
+ t->flgEigensysIsUptodate = 1;
+ cmaes_timings_start(&t->eigenTimings);
+ }
+ }
+
+ /* treat minimal standard deviations and numeric problems */
+ TestMinStdDevs(t);
+
+ for (iNk = 0; iNk < t->sp.lambda; ++iNk)
+ { /* generate scaled cmaes_random vector (D * z) */
+ for (i = 0; i < N; ++i)
+ if (flgdiag)
+ t->rgrgx[iNk][i] = xmean[i] + t->sigma * t->rgD[i] * cmaes_random_Gauss(&t->rand);
+ else
+ t->rgdTmp[i] = t->rgD[i] * cmaes_random_Gauss(&t->rand);
+ if (!flgdiag)
+ /* add mutation (sigma * B * (D*z)) */
+ for (i = 0; i < N; ++i) {
+ for (j = 0, sum = 0.; j < N; ++j)
+ sum += t->B[i][j] * t->rgdTmp[j];
+ t->rgrgx[iNk][i] = xmean[i] + t->sigma * sum;
+ }
+ }
+ if(t->state == 3 || t->gen == 0)
+ ++t->gen;
+ t->state = 1;
+
+ return(t->rgrgx);
+} /* SamplePopulation() */
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+double const *
+cmaes_ReSampleSingle_old( cmaes_t *t, double *rgx)
+{
+ int i, j, N=t->sp.N;
+ double sum;
+
+ if (rgx == NULL)
+ FATAL("cmaes_ReSampleSingle(): Missing input double *x",0,0,0);
+
+ for (i = 0; i < N; ++i)
+ t->rgdTmp[i] = t->rgD[i] * cmaes_random_Gauss(&t->rand);
+ /* add mutation (sigma * B * (D*z)) */
+ for (i = 0; i < N; ++i) {
+ for (j = 0, sum = 0.; j < N; ++j)
+ sum += t->B[i][j] * t->rgdTmp[j];
+ rgx[i] = t->rgxmean[i] + t->sigma * sum;
+ }
+ return rgx;
+}
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+double * const *
+cmaes_ReSampleSingle( cmaes_t *t, int iindex)
+{
+ int i, j, N=t->sp.N;
+ double *rgx;
+ double sum;
+ static char s[99];
+
+ if (iindex < 0 || iindex >= t->sp.lambda) {
+ sprintf(s, "index==%d must be between 0 and %d", iindex, t->sp.lambda);
+ FATAL("cmaes_ReSampleSingle(): Population member ",s,0,0);
+ }
+ rgx = t->rgrgx[iindex];
+
+ for (i = 0; i < N; ++i)
+ t->rgdTmp[i] = t->rgD[i] * cmaes_random_Gauss(&t->rand);
+ /* add mutation (sigma * B * (D*z)) */
+ for (i = 0; i < N; ++i) {
+ for (j = 0, sum = 0.; j < N; ++j)
+ sum += t->B[i][j] * t->rgdTmp[j];
+ rgx[i] = t->rgxmean[i] + t->sigma * sum;
+ }
+ return(t->rgrgx);
+}
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+double *
+cmaes_SampleSingleInto( cmaes_t *t, double *rgx)
+{
+ int i, j, N=t->sp.N;
+ double sum;
+
+ if (rgx == NULL)
+ rgx = new_double(N);
+
+ for (i = 0; i < N; ++i)
+ t->rgdTmp[i] = t->rgD[i] * cmaes_random_Gauss(&t->rand);
+ /* add mutation (sigma * B * (D*z)) */
+ for (i = 0; i < N; ++i) {
+ for (j = 0, sum = 0.; j < N; ++j)
+ sum += t->B[i][j] * t->rgdTmp[j];
+ rgx[i] = t->rgxmean[i] + t->sigma * sum;
+ }
+ return rgx;
+}
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+double *
+cmaes_PerturbSolutionInto( cmaes_t *t, double *rgx, double const *xmean, double eps)
+{
+ int i, j, N=t->sp.N;
+ double sum;
+
+ if (rgx == NULL)
+ rgx = new_double(N);
+ if (xmean == NULL)
+ FATAL("cmaes_PerturbSolutionInto(): xmean was not given",0,0,0);
+
+ for (i = 0; i < N; ++i)
+ t->rgdTmp[i] = t->rgD[i] * cmaes_random_Gauss(&t->rand);
+ /* add mutation (sigma * B * (D*z)) */
+ for (i = 0; i < N; ++i) {
+ for (j = 0, sum = 0.; j < N; ++j)
+ sum += t->B[i][j] * t->rgdTmp[j];
+ rgx[i] = xmean[i] + eps * t->sigma * sum;
+ }
+ return rgx;
+}
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+const double *
+cmaes_Optimize( cmaes_t *evo, double(*pFun)(double const *, int dim), long iterations)
+/* TODO: make signals.par another argument or, even better, part of cmaes_t */
+{
+ const char * signalsFilename = "cmaes_signals.par";
+ double *const*pop; /* sampled population */
+ const char *stop;
+ int i;
+ double startiter = evo->gen;
+
+ while(!(stop=cmaes_TestForTermination(evo)) &&
+ (evo->gen < startiter + iterations || !iterations))
+ {
+ /* Generate population of new candidate solutions */
+ pop = cmaes_SamplePopulation(evo); /* do not change content of pop */
+
+ /* Compute fitness value for each candidate solution */
+ for (i = 0; i < cmaes_Get(evo, "popsize"); ++i) {
+ evo->publicFitness[i] = (*pFun)(pop[i], evo->sp.N);
+ }
+
+ /* update search distribution */
+ cmaes_UpdateDistribution(evo, evo->publicFitness);
+
+ /* read control signals for output and termination */
+ if (signalsFilename)
+ cmaes_ReadSignals(evo, signalsFilename);
+ fflush(stdout);
+ } /* while !cmaes_TestForTermination(evo) */
+
+ /* write some data */
+ cmaes_WriteToFile(evo, "all", "allcmaes.dat");
+
+ return cmaes_GetPtr(evo, "xbestever");
+}
+
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+double *
+cmaes_UpdateDistribution( cmaes_t *t, const double *rgFunVal)
+{
+ int i, j, iNk, hsig, N=t->sp.N;
+ int flgdiag = ((t->sp.diagonalCov == 1) || (t->sp.diagonalCov >= t->gen));
+ double sum;
+ double psxps;
+
+ if(t->state == 3)
+ FATAL("cmaes_UpdateDistribution(): You need to call \n",
+ "SamplePopulation() before update can take place.",0,0);
+ if(rgFunVal == NULL)
+ FATAL("cmaes_UpdateDistribution(): ",
+ "Fitness function value array input is missing.",0,0);
+
+ if(t->state == 1) /* function values are delivered here */
+ t->countevals += t->sp.lambda;
+ else
+ ERRORMESSAGE("cmaes_UpdateDistribution(): unexpected state",0,0,0);
+
+ /* assign function values */
+ for (i=0; i < t->sp.lambda; ++i)
+ t->rgrgx[i][N] = t->rgFuncValue[i] = rgFunVal[i];
+
+
+ /* Generate index */
+ Sorted_index(rgFunVal, t->index, t->sp.lambda);
+
+ /* Test if function values are identical, escape flat fitness */
+ if (t->rgFuncValue[t->index[0]] ==
+ t->rgFuncValue[t->index[(int)t->sp.lambda/2]]) {
+ t->sigma *= exp(0.2+t->sp.cs/t->sp.damps);
+ ERRORMESSAGE("Warning: sigma increased due to equal function values\n",
+ " Reconsider the formulation of the objective function",0,0);
+ }
+
+ /* update function value history */
+ for(i = (int)*(t->arFuncValueHist-1)-1; i > 0; --i) /* for(i = t->arFuncValueHist[-1]-1; i > 0; --i) */
+ t->arFuncValueHist[i] = t->arFuncValueHist[i-1];
+ t->arFuncValueHist[0] = rgFunVal[t->index[0]];
+
+ /* update xbestever */
+ if (t->rgxbestever[N] > t->rgrgx[t->index[0]][N] || t->gen == 1)
+ for (i = 0; i <= N; ++i) {
+ t->rgxbestever[i] = t->rgrgx[t->index[0]][i];
+ t->rgxbestever[N+1] = t->countevals;
+ }
+
+ /* calculate xmean and rgBDz~N(0,C) */
+ for (i = 0; i < N; ++i) {
+ t->rgxold[i] = t->rgxmean[i];
+ t->rgxmean[i] = 0.;
+ for (iNk = 0; iNk < t->sp.mu; ++iNk)
+ t->rgxmean[i] += t->sp.weights[iNk] * t->rgrgx[t->index[iNk]][i];
+ t->rgBDz[i] = sqrt(t->sp.mueff)*(t->rgxmean[i] - t->rgxold[i])/t->sigma;
+ }
+
+ /* calculate z := D^(-1) * B^(-1) * rgBDz into rgdTmp */
+ for (i = 0; i < N; ++i) {
+ if (!flgdiag)
+ for (j = 0, sum = 0.; j < N; ++j)
+ sum += t->B[j][i] * t->rgBDz[j];
+ else
+ sum = t->rgBDz[i];
+ t->rgdTmp[i] = sum / t->rgD[i];
+ }
+
+ /* TODO?: check length of t->rgdTmp and set an upper limit, e.g. 6 stds */
+ /* in case of manipulation of arx,
+ this can prevent an increase of sigma by several orders of magnitude
+ within one step; a five-fold increase in one step can still happen.
+ */
+ /*
+ for (j = 0, sum = 0.; j < N; ++j)
+ sum += t->rgdTmp[j] * t->rgdTmp[j];
+ if (sqrt(sum) > chiN + 6. * sqrt(0.5)) {
+ rgdTmp length should be set to upper bound and hsig should become zero
+ }
+ */
+
+ /* cumulation for sigma (ps) using B*z */
+ for (i = 0; i < N; ++i) {
+ if (!flgdiag)
+ for (j = 0, sum = 0.; j < N; ++j)
+ sum += t->B[i][j] * t->rgdTmp[j];
+ else
+ sum = t->rgdTmp[i];
+ t->rgps[i] = (1. - t->sp.cs) * t->rgps[i] +
+ sqrt(t->sp.cs * (2. - t->sp.cs)) * sum;
+ }
+
+ /* calculate norm(ps)^2 */
+ for (i = 0, psxps = 0.; i < N; ++i)
+ psxps += t->rgps[i] * t->rgps[i];
+
+ /* cumulation for covariance matrix (pc) using B*D*z~N(0,C) */
+ hsig = sqrt(psxps) / sqrt(1. - pow(1.-t->sp.cs, 2*t->gen)) / t->chiN
+ < 1.4 + 2./(N+1);
+ for (i = 0; i < N; ++i) {
+ t->rgpc[i] = (1. - t->sp.ccumcov) * t->rgpc[i] +
+ hsig * sqrt(t->sp.ccumcov * (2. - t->sp.ccumcov)) * t->rgBDz[i];
+ }
+
+ /* stop initial phase */
+ if (t->flgIniphase &&
+ t->gen > douMin(1/t->sp.cs, 1+N/t->sp.mucov))
+ {
+ if (psxps / t->sp.damps / (1.-pow((1. - t->sp.cs), t->gen))
+ < N * 1.05)
+ t->flgIniphase = 0;
+ }
+
+#if 0
+ /* remove momentum in ps, if ps is large and fitness is getting worse */
+ /* This is obsolete due to hsig and harmful in a dynamic environment */
+ if(psxps/N > 1.5 + 10.*sqrt(2./N)
+ && t->arFuncValueHist[0] > t->arFuncValueHist[1]
+ && t->arFuncValueHist[0] > t->arFuncValueHist[2]) {
+ double tfac = sqrt((1 + douMax(0, log(psxps/N))) * N / psxps);
+ for (i=0; i<N; ++i)
+ t->rgps[i] *= tfac;
+ psxps *= tfac*tfac;
+ }
+#endif
+
+ /* update of C */
+
+ Adapt_C2(t, hsig);
+
+ /* Adapt_C(t); not used anymore */
+
+#if 0
+ if (t->sp.ccov != 0. && t->flgIniphase == 0) {
+ int k;
+
+ t->flgEigensysIsUptodate = 0;
+
+ /* update covariance matrix */
+ for (i = 0; i < N; ++i)
+ for (j = 0; j <=i; ++j) {
+ t->C[i][j] = (1 - t->sp.ccov) * t->C[i][j]
+ + t->sp.ccov * (1./t->sp.mucov)
+ * (t->rgpc[i] * t->rgpc[j]
+ + (1-hsig)*t->sp.ccumcov*(2.-t->sp.ccumcov) * t->C[i][j]);
+ for (k = 0; k < t->sp.mu; ++k) /* additional rank mu update */
+ t->C[i][j] += t->sp.ccov * (1-1./t->sp.mucov) * t->sp.weights[k]
+ * (t->rgrgx[t->index[k]][i] - t->rgxold[i])
+ * (t->rgrgx[t->index[k]][j] - t->rgxold[j])
+ / t->sigma / t->sigma;
+ }
+ }
+#endif
+
+
+ /* update of sigma */
+ t->sigma *= exp(((sqrt(psxps)/t->chiN)-1.)*t->sp.cs/t->sp.damps);
+
+ t->state = 3;
+
+ return (t->rgxmean);
+
+} /* cmaes_UpdateDistribution() */
+
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+static void
+Adapt_C2(cmaes_t *t, int hsig)
+{
+ int i, j, k, N=t->sp.N;
+ int flgdiag = ((t->sp.diagonalCov == 1) || (t->sp.diagonalCov >= t->gen));
+
+ if (t->sp.ccov != 0. && t->flgIniphase == 0) {
+
+ /* definitions for speeding up inner-most loop */
+ double ccov1 = douMin(t->sp.ccov * (1./t->sp.mucov) * (flgdiag ? (N+1.5) / 3. : 1.), 1.);
+ double ccovmu = douMin(t->sp.ccov * (1-1./t->sp.mucov)* (flgdiag ? (N+1.5) / 3. : 1.), 1.-ccov1);
+ double sigmasquare = t->sigma * t->sigma;
+
+ t->flgEigensysIsUptodate = 0;
+
+ /* update covariance matrix */
+ for (i = 0; i < N; ++i)
+ for (j = flgdiag ? i : 0; j <= i; ++j) {
+ t->C[i][j] = (1 - ccov1 - ccovmu) * t->C[i][j]
+ + ccov1
+ * (t->rgpc[i] * t->rgpc[j]
+ + (1-hsig)*t->sp.ccumcov*(2.-t->sp.ccumcov) * t->C[i][j]);
+ for (k = 0; k < t->sp.mu; ++k) { /* additional rank mu update */
+ t->C[i][j] += ccovmu * t->sp.weights[k]
+ * (t->rgrgx[t->index[k]][i] - t->rgxold[i])
+ * (t->rgrgx[t->index[k]][j] - t->rgxold[j])
+ / sigmasquare;
+ }
+ }
+ /* update maximal and minimal diagonal value */
+ t->maxdiagC = t->mindiagC = t->C[0][0];
+ for (i = 1; i < N; ++i) {
+ if (t->maxdiagC < t->C[i][i])
+ t->maxdiagC = t->C[i][i];
+ else if (t->mindiagC > t->C[i][i])
+ t->mindiagC = t->C[i][i];
+ }
+ } /* if ccov... */
+}
+
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+static void
+TestMinStdDevs(cmaes_t *t)
+ /* increases sigma */
+{
+ int i, N = t->sp.N;
+ if (t->sp.rgDiffMinChange == NULL)
+ return;
+
+ for (i = 0; i < N; ++i)
+ while (t->sigma * sqrt(t->C[i][i]) < t->sp.rgDiffMinChange[i])
+ t->sigma *= exp(0.05+t->sp.cs/t->sp.damps);
+
+} /* cmaes_TestMinStdDevs() */
+
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+void cmaes_WriteToFile(cmaes_t *t, const char *key, const char *name)
+{
+ cmaes_WriteToFileAW(t, key, name, "a"); /* default is append */
+}
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+void cmaes_WriteToFileAW(cmaes_t *t, const char *key, const char *name,
+ const char *appendwrite)
+{
+ const char *s = "tmpcmaes.dat";
+ FILE *fp;
+
+ if (name == NULL)
+ name = s;
+
+ fp = fopen( name, appendwrite);
+
+ if(fp == NULL) {
+ ERRORMESSAGE("cmaes_WriteToFile(): could not open '", name,
+ "' with flag ", appendwrite);
+ return;
+ }
+
+ if (appendwrite[0] == 'w') {
+ /* write a header line, very rudimentary */
+ fprintf(fp, "%% # %s (randomSeed=%d, %s)\n", key, t->sp.seed, getTimeStr());
+ } else
+ if (t->gen > 0 || strncmp(name, "outcmaesfit", 11) != 0)
+ cmaes_WriteToFilePtr(t, key, fp); /* do not write fitness for gen==0 */
+
+ fclose(fp);
+
+} /* WriteToFile */
+
+/* --------------------------------------------------------- */
+void cmaes_WriteToFilePtr(cmaes_t *t, const char *key, FILE *fp)
+
+/* this hack reads key words from input key for data to be written to
+ * a file, see file signals.par as input file. The length of the keys
+ * is mostly fixed, see key += number in the code! If the key phrase
+ * does not match the expectation the output might be strange. for
+ * cmaes_t *t == NULL it solely prints key as a header line. Input key
+ * must be zero terminated.
+ */
+{
+ int i, k, N=(t ? t->sp.N : 0);
+ char const *keyend; /* *keystart; */
+ const char *s = "few";
+ if (key == NULL)
+ key = s;
+ /* keystart = key; for debugging purpose */
+ keyend = key + strlen(key);
+
+ while (key < keyend)
+ {
+ if (strncmp(key, "axisratio", 9) == 0)
+ {
+ fprintf(fp, "%.2e", sqrt(t->maxEW/t->minEW));
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "idxminSD", 8) == 0)
+ {
+ int mini=0; for(i=N-1;i>0;--i) if(t->mindiagC==t->C[i][i]) mini=i;
+ fprintf(fp, "%d", mini+1);
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "idxmaxSD", 8) == 0)
+ {
+ int maxi=0; for(i=N-1;i>0;--i) if(t->maxdiagC==t->C[i][i]) maxi=i;
+ fprintf(fp, "%d", maxi+1);
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ /* new coordinate system == all eigenvectors */
+ if (strncmp(key, "B", 1) == 0)
+ {
+ /* int j, index[N]; */
+ int j, *iindex=(int*)(new_void(N,sizeof(int))); /* MT */
+ Sorted_index(t->rgD, iindex, N); /* should not be necessary, see end of QLalgo2 */
+ /* One eigenvector per row, sorted: largest eigenvalue first */
+ for (i = 0; i < N; ++i)
+ for (j = 0; j < N; ++j)
+ fprintf(fp, "%g%c", t->B[j][iindex[N-1-i]], (j==N-1)?'\n':'\t');
+ ++key;
+ free(iindex); /* MT */
+ }
+ /* covariance matrix */
+ if (strncmp(key, "C", 1) == 0)
+ {
+ int j;
+ for (i = 0; i < N; ++i)
+ for (j = 0; j <= i; ++j)
+ fprintf(fp, "%g%c", t->C[i][j], (j==i)?'\n':'\t');
+ ++key;
+ }
+ /* (processor) time (used) since begin of execution */
+ if (strncmp(key, "clock", 4) == 0)
+ {
+ cmaes_timings_update(&t->eigenTimings);
+ fprintf(fp, "%.1f %.1f", t->eigenTimings.totaltotaltime,
+ t->eigenTimings.tictoctime);
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ /* ratio between largest and smallest standard deviation */
+ if (strncmp(key, "stddevratio", 11) == 0) /* std dev in coordinate axes */
+ {
+ fprintf(fp, "%g", sqrt(t->maxdiagC/t->mindiagC));
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ /* standard deviations in coordinate directions (sigma*sqrt(C[i,i])) */
+ if (strncmp(key, "coorstddev", 10) == 0
+ || strncmp(key, "stddev", 6) == 0) /* std dev in coordinate axes */
+ {
+ for (i = 0; i < N; ++i)
+ fprintf(fp, "%s%g", (i==0) ? "":"\t", t->sigma*sqrt(t->C[i][i]));
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ /* diagonal of D == roots of eigenvalues, sorted */
+ if (strncmp(key, "diag(D)", 7) == 0)
+ {
+ for (i = 0; i < N; ++i)
+ t->rgdTmp[i] = t->rgD[i];
+ qsort(t->rgdTmp, (unsigned) N, sizeof(double), &SignOfDiff); /* superfluous */
+ for (i = 0; i < N; ++i)
+ fprintf(fp, "%s%g", (i==0) ? "":"\t", t->rgdTmp[i]);
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "dim", 3) == 0)
+ {
+ fprintf(fp, "%d", N);
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "eval", 4) == 0)
+ {
+ fprintf(fp, "%.0f", t->countevals);
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "few(diag(D))", 12) == 0)/* between four and six axes */
+ {
+ int add = (int)(0.5 + (N + 1.) / 5.);
+ for (i = 0; i < N; ++i)
+ t->rgdTmp[i] = t->rgD[i];
+ qsort(t->rgdTmp, (unsigned) N, sizeof(double), &SignOfDiff);
+ for (i = 0; i < N-1; i+=add) /* print always largest */
+ fprintf(fp, "%s%g", (i==0) ? "":"\t", t->rgdTmp[N-1-i]);
+ fprintf(fp, "\t%g\n", t->rgdTmp[0]); /* and smallest */
+ break; /* number of printed values is not determined */
+ }
+ if (strncmp(key, "fewinfo", 7) == 0) {
+ fprintf(fp," Iter Fevals Function Value Sigma ");
+ fprintf(fp, "MaxCoorDev MinCoorDev AxisRatio MinDii Time in eig\n");
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ }
+ if (strncmp(key, "few", 3) == 0) {
+ fprintf(fp, " %4.0f ", t->gen);
+ fprintf(fp, " %5.0f ", t->countevals);
+ fprintf(fp, "%.15e", t->rgFuncValue[t->index[0]]);
+ fprintf(fp, " %.2e %.2e %.2e", t->sigma, t->sigma*sqrt(t->maxdiagC),
+ t->sigma*sqrt(t->mindiagC));
+ fprintf(fp, " %.2e %.2e", sqrt(t->maxEW/t->minEW), sqrt(t->minEW));
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "funval", 6) == 0 || strncmp(key, "fitness", 6) == 0)
+ {
+ fprintf(fp, "%.15e", t->rgFuncValue[t->index[0]]);
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "fbestever", 9) == 0)
+ {
+ fprintf(fp, "%.15e", t->rgxbestever[N]); /* f-value */
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "fmedian", 7) == 0)
+ {
+ fprintf(fp, "%.15e", t->rgFuncValue[t->index[(int)(t->sp.lambda/2)]]);
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "fworst", 6) == 0)
+ {
+ fprintf(fp, "%.15e", t->rgFuncValue[t->index[t->sp.lambda-1]]);
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "arfunval", 8) == 0 || strncmp(key, "arfitness", 8) == 0)
+ {
+ for (i = 0; i < N; ++i)
+ fprintf(fp, "%s%.10e", (i==0) ? "" : "\t",
+ t->rgFuncValue[t->index[i]]);
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "gen", 3) == 0)
+ {
+ fprintf(fp, "%.0f", t->gen);
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "iter", 4) == 0)
+ {
+ fprintf(fp, "%.0f", t->gen);
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "sigma", 5) == 0)
+ {
+ fprintf(fp, "%.4e", t->sigma);
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "minSD", 5) == 0) /* minimal standard deviation */
+ {
+ fprintf(fp, "%.4e", sqrt(t->mindiagC));
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "maxSD", 5) == 0)
+ {
+ fprintf(fp, "%.4e", sqrt(t->maxdiagC));
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "mindii", 6) == 0)
+ {
+ fprintf(fp, "%.4e", sqrt(t->minEW));
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "0", 1) == 0)
+ {
+ fprintf(fp, "0");
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "lambda", 6) == 0 || strncmp(key, "popsi", 5) == 0 || strncmp(key, "populationsi", 12) == 0)
+ {
+ fprintf(fp, "%d", t->sp.lambda);
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "N", 1) == 0)
+ {
+ fprintf(fp, "%d", N);
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "resume", 6) == 0)
+ {
+ fprintf(fp, "\n# resume %d\n", N);
+ fprintf(fp, "xmean\n");
+ cmaes_WriteToFilePtr(t, "xmean", fp);
+ fprintf(fp, "path for sigma\n");
+ for(i=0; i<N; ++i)
+ fprintf(fp, "%g%s", t->rgps[i], (i==N-1) ? "\n":"\t");
+ fprintf(fp, "path for C\n");
+ for(i=0; i<N; ++i)
+ fprintf(fp, "%g%s", t->rgpc[i], (i==N-1) ? "\n":"\t");
+ fprintf(fp, "sigma %g\n", t->sigma);
+ /* note than B and D might not be up-to-date */
+ fprintf(fp, "covariance matrix\n");
+ cmaes_WriteToFilePtr(t, "C", fp);
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ }
+ if (strncmp(key, "xbest", 5) == 0) { /* best x in recent generation */
+ for(i=0; i<N; ++i)
+ fprintf(fp, "%s%g", (i==0) ? "":"\t", t->rgrgx[t->index[0]][i]);
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "xmean", 5) == 0) {
+ for(i=0; i<N; ++i)
+ fprintf(fp, "%s%g", (i==0) ? "":"\t", t->rgxmean[i]);
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ fprintf(fp, "%c", (*key=='+') ? '\t':'\n');
+ }
+ if (strncmp(key, "all", 3) == 0)
+ {
+ time_t ti = time(NULL);
+ fprintf(fp, "\n# --------- %s\n", asctime(localtime(&ti)));
+ fprintf(fp, " N %d\n", N);
+ fprintf(fp, " seed %d\n", t->sp.seed);
+ fprintf(fp, "function evaluations %.0f\n", t->countevals);
+ fprintf(fp, "elapsed (CPU) time [s] %.2f\n", t->eigenTimings.totaltotaltime);
+ fprintf(fp, "function value f(x)=%g\n", t->rgrgx[t->index[0]][N]);
+ fprintf(fp, "maximal standard deviation %g\n", t->sigma*sqrt(t->maxdiagC));
+ fprintf(fp, "minimal standard deviation %g\n", t->sigma*sqrt(t->mindiagC));
+ fprintf(fp, "sigma %g\n", t->sigma);
+ fprintf(fp, "axisratio %g\n", rgdouMax(t->rgD, N)/rgdouMin(t->rgD, N));
+ fprintf(fp, "xbestever found after %.0f evaluations, function value %g\n",
+ t->rgxbestever[N+1], t->rgxbestever[N]);
+ for(i=0; i<N; ++i)
+ fprintf(fp, " %12g%c", t->rgxbestever[i],
+ (i%5==4||i==N-1)?'\n':' ');
+ fprintf(fp, "xbest (of last generation, function value %g)\n",
+ t->rgrgx[t->index[0]][N]);
+ for(i=0; i<N; ++i)
+ fprintf(fp, " %12g%c", t->rgrgx[t->index[0]][i],
+ (i%5==4||i==N-1)?'\n':' ');
+ fprintf(fp, "xmean \n");
+ for(i=0; i<N; ++i)
+ fprintf(fp, " %12g%c", t->rgxmean[i],
+ (i%5==4||i==N-1)?'\n':' ');
+ fprintf(fp, "Standard deviation of coordinate axes (sigma*sqrt(diag(C)))\n");
+ for(i=0; i<N; ++i)
+ fprintf(fp, " %12g%c", t->sigma*sqrt(t->C[i][i]),
+ (i%5==4||i==N-1)?'\n':' ');
+ fprintf(fp, "Main axis lengths of mutation ellipsoid (sigma*diag(D))\n");
+ for (i = 0; i < N; ++i)
+ t->rgdTmp[i] = t->rgD[i];
+ qsort(t->rgdTmp, (unsigned) N, sizeof(double), &SignOfDiff);
+ for(i=0; i<N; ++i)
+ fprintf(fp, " %12g%c", t->sigma*t->rgdTmp[N-1-i],
+ (i%5==4||i==N-1)?'\n':' ');
+ fprintf(fp, "Longest axis (b_i where d_ii=max(diag(D))\n");
+ k = MaxIdx(t->rgD, N);
+ for(i=0; i<N; ++i)
+ fprintf(fp, " %12g%c", t->B[i][k], (i%5==4||i==N-1)?'\n':' ');
+ fprintf(fp, "Shortest axis (b_i where d_ii=max(diag(D))\n");
+ k = MinIdx(t->rgD, N);
+ for(i=0; i<N; ++i)
+ fprintf(fp, " %12g%c", t->B[i][k], (i%5==4||i==N-1)?'\n':' ');
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ } /* "all" */
+
+#if 0 /* could become generic part */
+ s0 = key;
+ d = cmaes_Get(t, key); /* TODO find way to detect whether key was found */
+ if (key == s0) /* this does not work, is always true */
+ {
+ /* write out stuff, problem: only generic format is available */
+ /* move in key until "+" or end */
+ }
+#endif
+
+ if (*key == '\0')
+ break;
+ else if (*key != '+') { /* last key was not recognized */
+ ERRORMESSAGE("cmaes_t:WriteToFilePtr(): unrecognized key '", key, "'", 0);
+ while (*key != '+' && *key != '\0' && key < keyend)
+ ++key;
+ }
+ while (*key == '+')
+ ++key;
+ } /* while key < keyend */
+
+ if (key > keyend)
+ FATAL("cmaes_t:WriteToFilePtr(): BUG regarding key sequence",0,0,0);
+
+} /* WriteToFilePtr */
+
+/* --------------------------------------------------------- */
+double
+cmaes_Get( cmaes_t *t, char const *s)
+{
+ int N=t->sp.N;
+
+ if (strncmp(s, "axisratio", 5) == 0) { /* between lengths of longest and shortest principal axis of the distribution ellipsoid */
+ return (rgdouMax(t->rgD, N)/rgdouMin(t->rgD, N));
+ }
+ else if (strncmp(s, "eval", 4) == 0) { /* number of function evaluations */
+ return (t->countevals);
+ }
+ else if (strncmp(s, "fctvalue", 6) == 0
+ || strncmp(s, "funcvalue", 6) == 0
+ || strncmp(s, "funvalue", 6) == 0
+ || strncmp(s, "fitness", 3) == 0) { /* recent best function value */
+ return(t->rgFuncValue[t->index[0]]);
+ }
+ else if (strncmp(s, "fbestever", 7) == 0) { /* ever best function value */
+ return(t->rgxbestever[N]);
+ }
+ else if (strncmp(s, "generation", 3) == 0
+ || strncmp(s, "iteration", 4) == 0) {
+ return(t->gen);
+ }
+ else if (strncmp(s, "maxeval", 4) == 0
+ || strncmp(s, "MaxFunEvals", 8) == 0
+ || strncmp(s, "stopMaxFunEvals", 12) == 0) { /* maximal number of function evaluations */
+ return(t->sp.stopMaxFunEvals);
+ }
+ else if (strncmp(s, "maxgen", 4) == 0
+ || strncmp(s, "MaxIter", 7) == 0
+ || strncmp(s, "stopMaxIter", 11) == 0) { /* maximal number of generations */
+ return(ceil(t->sp.stopMaxIter));
+ }
+ else if (strncmp(s, "maxaxislength", 5) == 0) { /* sigma * max(diag(D)) */
+ return(t->sigma * sqrt(t->maxEW));
+ }
+ else if (strncmp(s, "minaxislength", 5) == 0) { /* sigma * min(diag(D)) */
+ return(t->sigma * sqrt(t->minEW));
+ }
+ else if (strncmp(s, "maxstddev", 4) == 0) { /* sigma * sqrt(max(diag(C))) */
+ return(t->sigma * sqrt(t->maxdiagC));
+ }
+ else if (strncmp(s, "minstddev", 4) == 0) { /* sigma * sqrt(min(diag(C))) */
+ return(t->sigma * sqrt(t->mindiagC));
+ }
+ else if (strncmp(s, "N", 1) == 0 || strcmp(s, "n") == 0 ||
+ strncmp(s, "dimension", 3) == 0) {
+ return (N);
+ }
+ else if (strncmp(s, "lambda", 3) == 0
+ || strncmp(s, "samplesize", 8) == 0
+ || strncmp(s, "popsize", 7) == 0) { /* sample size, offspring population size */
+ return(t->sp.lambda);
+ }
+ else if (strncmp(s, "sigma", 3) == 0) {
+ return(t->sigma);
+ }
+ FATAL( "cmaes_Get(cmaes_t, char const * s): No match found for s='", s, "'",0);
+ return(0);
+} /* cmaes_Get() */
+
+/* --------------------------------------------------------- */
+double *
+cmaes_GetInto( cmaes_t *t, char const *s, double *res)
+{
+ int i, N = t->sp.N;
+ double const * res0 = cmaes_GetPtr(t, s);
+ if (res == NULL)
+ res = new_double(N);
+ for (i = 0; i < N; ++i)
+ res[i] = res0[i];
+ return res;
+}
+
+/* --------------------------------------------------------- */
+double *
+cmaes_GetNew( cmaes_t *t, char const *s)
+{
+ return (cmaes_GetInto(t, s, NULL));
+}
+
+/* --------------------------------------------------------- */
+const double *
+cmaes_GetPtr( cmaes_t *t, char const *s)
+{
+ int i, N=t->sp.N;
+
+ /* diagonal of covariance matrix */
+ if (strncmp(s, "diag(C)", 7) == 0) {
+ for (i = 0; i < N; ++i)
+ t->rgout[i] = t->C[i][i];
+ return(t->rgout);
+ }
+ /* diagonal of axis lengths matrix */
+ else if (strncmp(s, "diag(D)", 7) == 0) {
+ return(t->rgD);
+ }
+ /* vector of standard deviations sigma*sqrt(diag(C)) */
+ else if (strncmp(s, "stddev", 3) == 0) {
+ for (i = 0; i < N; ++i)
+ t->rgout[i] = t->sigma * sqrt(t->C[i][i]);
+ return(t->rgout);
+ }
+ /* bestever solution seen so far */
+ else if (strncmp(s, "xbestever", 7) == 0)
+ return(t->rgxbestever);
+ /* recent best solution of the recent population */
+ else if (strncmp(s, "xbest", 5) == 0)
+ return(t->rgrgx[t->index[0]]);
+ /* mean of the recent distribution */
+ else if (strncmp(s, "xmean", 1) == 0)
+ return(t->rgxmean);
+
+ return(NULL);
+}
+
+/* --------------------------------------------------------- */
+/* tests stopping criteria
+ * returns a string of satisfied stopping criterion for each line
+ * otherwise NULL
+*/
+const char *
+cmaes_TestForTermination( cmaes_t *t)
+{
+ double range, fac;
+ int iAchse, iKoo;
+ int flgdiag = ((t->sp.diagonalCov == 1) || (t->sp.diagonalCov >= t->gen));
+ static char sTestOutString[3024];
+ char * cp = sTestOutString;
+ int i, cTemp, N=t->sp.N;
+ cp[0] = '\0';
+
+ /* function value reached */
+ if ((t->gen > 1 || t->state > 1) && t->sp.stStopFitness.flg &&
+ t->rgFuncValue[t->index[0]] <= t->sp.stStopFitness.val)
+ cp += sprintf(cp, "Fitness: function value %7.2e <= stopFitness (%7.2e)\n",
+ t->rgFuncValue[t->index[0]], t->sp.stStopFitness.val);
+
+ /* TolFun */
+ range = douMax(rgdouMax(t->arFuncValueHist, (int)douMin(t->gen,*(t->arFuncValueHist-1))),
+ rgdouMax(t->rgFuncValue, t->sp.lambda)) -
+ douMin(rgdouMin(t->arFuncValueHist, (int)douMin(t->gen, *(t->arFuncValueHist-1))),
+ rgdouMin(t->rgFuncValue, t->sp.lambda));
+
+ if (t->gen > 0 && range <= t->sp.stopTolFun) {
+ cp += sprintf(cp,
+ "TolFun: function value differences %7.2e < stopTolFun=%7.2e\n",
+ range, t->sp.stopTolFun);
+ }
+
+ /* TolFunHist */
+ if (t->gen > *(t->arFuncValueHist-1)) {
+ range = rgdouMax(t->arFuncValueHist, (int)*(t->arFuncValueHist-1))
+ - rgdouMin(t->arFuncValueHist, (int)*(t->arFuncValueHist-1));
+ if (range <= t->sp.stopTolFunHist)
+ cp += sprintf(cp,
+ "TolFunHist: history of function value changes %7.2e stopTolFunHist=%7.2e",
+ range, t->sp.stopTolFunHist);
+ }
+
+ /* TolX */
+ for(i=0, cTemp=0; i<N; ++i) {
+ cTemp += (t->sigma * sqrt(t->C[i][i]) < t->sp.stopTolX) ? 1 : 0;
+ cTemp += (t->sigma * t->rgpc[i] < t->sp.stopTolX) ? 1 : 0;
+ }
+ if (cTemp == 2*N) {
+ cp += sprintf(cp,
+ "TolX: object variable changes below %7.2e \n",
+ t->sp.stopTolX);
+ }
+
+ /* TolUpX */
+ for(i=0; i<N; ++i) {
+ if (t->sigma * sqrt(t->C[i][i]) > t->sp.stopTolUpXFactor * t->sp.rgInitialStds[i])
+ break;
+ }
+ if (i < N) {
+ cp += sprintf(cp,
+ "TolUpX: standard deviation increased by more than %7.2e, larger initial standard deviation recommended \n",
+ t->sp.stopTolUpXFactor);
+ }
+
+ /* Condition of C greater than dMaxSignifKond */
+ if (t->maxEW >= t->minEW * t->dMaxSignifKond) {
+ cp += sprintf(cp,
+ "ConditionNumber: maximal condition number %7.2e reached. maxEW=%7.2e,minEW=%7.2e,maxdiagC=%7.2e,mindiagC=%7.2e\n",
+ t->dMaxSignifKond, t->maxEW, t->minEW, t->maxdiagC, t->mindiagC);
+ } /* if */
+
+ /* Principal axis i has no effect on xmean, ie.
+ x == x + 0.1 * sigma * rgD[i] * B[i] */
+ if (!flgdiag) {
+ for (iAchse = 0; iAchse < N; ++iAchse)
+ {
+ fac = 0.1 * t->sigma * t->rgD[iAchse];
+ for (iKoo = 0; iKoo < N; ++iKoo){
+ if (t->rgxmean[iKoo] != t->rgxmean[iKoo] + fac * t->B[iKoo][iAchse])
+ break;
+ }
+ if (iKoo == N)
+ {
+ /* t->sigma *= exp(0.2+t->sp.cs/t->sp.damps); */
+ cp += sprintf(cp,
+ "NoEffectAxis: standard deviation 0.1*%7.2e in principal axis %d without effect\n",
+ fac/0.1, iAchse);
+ break;
+ } /* if (iKoo == N) */
+ } /* for iAchse */
+ } /* if flgdiag */
+ /* Component of xmean is not changed anymore */
+ for (iKoo = 0; iKoo < N; ++iKoo)
+ {
+ if (t->rgxmean[iKoo] == t->rgxmean[iKoo] +
+ 0.2*t->sigma*sqrt(t->C[iKoo][iKoo]))
+ {
+ /* t->C[iKoo][iKoo] *= (1 + t->sp.ccov); */
+ /* flg = 1; */
+ cp += sprintf(cp,
+ "NoEffectCoordinate: standard deviation 0.2*%7.2e in coordinate %d without effect\n",
+ t->sigma*sqrt(t->C[iKoo][iKoo]), iKoo);
+ break;
+ }
+
+ } /* for iKoo */
+ /* if (flg) t->sigma *= exp(0.05+t->sp.cs/t->sp.damps); */
+
+ if(t->countevals >= t->sp.stopMaxFunEvals)
+ cp += sprintf(cp, "MaxFunEvals: conducted function evaluations %.0f >= %g\n",
+ t->countevals, t->sp.stopMaxFunEvals);
+ if(t->gen >= t->sp.stopMaxIter)
+ cp += sprintf(cp, "MaxIter: number of iterations %.0f >= %g\n",
+ t->gen, t->sp.stopMaxIter);
+ if(t->flgStop)
+ cp += sprintf(cp, "Manual: stop signal read\n");
+
+#if 0
+ else if (0) {
+ for(i=0, cTemp=0; i<N; ++i) {
+ cTemp += (sigma * sqrt(C[i][i]) < stopdx) ? 1 : 0;
+ cTemp += (sigma * rgpc[i] < stopdx) ? 1 : 0;
+ }
+ if (cTemp == 2*N)
+ flgStop = 1;
+ }
+#endif
+
+ if (cp - sTestOutString>320)
+ ERRORMESSAGE("Bug in cmaes_t:Test(): sTestOutString too short",0,0,0);
+
+ if (cp != sTestOutString) {
+ return sTestOutString;
+ }
+
+ return(NULL);
+
+} /* cmaes_Test() */
+
+/* --------------------------------------------------------- */
+void cmaes_ReadSignals(cmaes_t *t, char const *filename)
+{
+ const char *s = "cmaes_signals.par";
+ FILE *fp;
+ if (filename == NULL)
+ filename = s;
+/* if (filename) assign_string(&(t->signalsFilename), filename)*/
+ fp = fopen( filename, "r");
+ if(fp == NULL) {
+ return;
+ }
+ cmaes_ReadFromFilePtr( t, fp);
+ fclose(fp);
+}
+/* --------------------------------------------------------- */
+void cmaes_ReadFromFilePtr( cmaes_t *t, FILE *fp)
+/* reading commands e.g. from signals.par file
+*/
+{
+ const char *keys[15]; /* key strings for scanf */
+ char s[199], sin1[99], sin2[129], sin3[99], sin4[99];
+ int ikey, ckeys, nb;
+ double d;
+ static int flglockprint = 0;
+ static int flglockwrite = 0;
+ static long countiterlastwritten;
+ static long maxdiffitertowrite; /* to prevent long gaps at the beginning */
+ int flgprinted = 0;
+ int flgwritten = 0;
+ double deltaprinttime = (double)(time(NULL)-t->printtime); /* using clock instead might not be a good */
+ double deltawritetime = (double)(time(NULL)-t->writetime); /* idea as disc time is not CPU time? */
+ double deltaprinttimefirst = (double)(t->firstprinttime ? time(NULL)-t->firstprinttime : 0); /* time is in seconds!? */
+ double deltawritetimefirst = (double)(t->firstwritetime ? time(NULL)-t->firstwritetime : 0);
+ if (countiterlastwritten > t->gen) { /* probably restarted */
+ maxdiffitertowrite = 0;
+ countiterlastwritten = 0;
+ }
+
+ keys[0] = " stop%98s %98s"; /* s=="now" or eg "MaxIter+" %lg"-number */
+ /* works with and without space */
+ keys[1] = " print %98s %98s"; /* s==keyword for WriteFile */
+ keys[2] = " write %98s %128s %98s"; /* s1==keyword, s2==filename */
+ keys[3] = " check%98s %98s";
+ keys[4] = " maxTimeFractionForEigendecompostion %98s";
+ ckeys = 5;
+ strcpy(sin2, "tmpcmaes.dat");
+
+ if (cmaes_TestForTermination(t))
+ {
+ deltaprinttime = (double)time(NULL); /* forces printing */
+ deltawritetime = (double)time(NULL);
+ }
+ while(fgets(s, sizeof(s), fp) != NULL)
+ {
+ if (s[0] == '#' || s[0] == '%') /* skip comments */
+ continue;
+ sin1[0] = sin2[0] = sin3[0] = sin4[0] = '\0';
+ for (ikey=0; ikey < ckeys; ++ikey)
+ {
+ if((nb=sscanf(s, keys[ikey], sin1, sin2, sin3, sin4)) >= 1)
+ {
+ switch(ikey) {
+ case 0 : /* "stop", reads "stop now" or eg. stopMaxIter */
+ if (strncmp(sin1, "now", 3) == 0)
+ t->flgStop = 1;
+ else if (strncmp(sin1, "MaxFunEvals", 11) == 0) {
+ if (sscanf(sin2, " %lg", &d) == 1)
+ t->sp.stopMaxFunEvals = d;
+ }
+ else if (strncmp(sin1, "MaxIter", 4) == 0) {
+ if (sscanf(sin2, " %lg", &d) == 1)
+ t->sp.stopMaxIter = d;
+ }
+ else if (strncmp(sin1, "Fitness", 7) == 0) {
+ if (sscanf(sin2, " %lg", &d) == 1)
+ {
+ t->sp.stStopFitness.flg = 1;
+ t->sp.stStopFitness.val = d;
+ }
+ }
+ else if (strncmp(sin1, "TolFunHist", 10) == 0) {
+ if (sscanf(sin2, " %lg", &d) == 1)
+ t->sp.stopTolFunHist = d;
+ }
+ else if (strncmp(sin1, "TolFun", 6) == 0) {
+ if (sscanf(sin2, " %lg", &d) == 1)
+ t->sp.stopTolFun = d;
+ }
+ else if (strncmp(sin1, "TolX", 4) == 0) {
+ if (sscanf(sin2, " %lg", &d) == 1)
+ t->sp.stopTolX = d;
+ }
+ else if (strncmp(sin1, "TolUpXFactor", 4) == 0) {
+ if (sscanf(sin2, " %lg", &d) == 1)
+ t->sp.stopTolUpXFactor = d;
+ }
+ break;
+ case 1 : /* "print" */
+ d = 1; /* default */
+ if (sscanf(sin2, "%lg", &d) < 1 && deltaprinttimefirst < 1)
+ d = 0; /* default at first time */
+ if (deltaprinttime >= d && !flglockprint) {
+ cmaes_WriteToFilePtr(t, sin1, stdout);
+ flgprinted = 1;
+ }
+ if(d < 0)
+ flglockprint += 2;
+ break;
+ case 2 : /* "write" */
+ /* write header, before first generation */
+ if (t->countevals < t->sp.lambda && t->flgresumedone == 0)
+ cmaes_WriteToFileAW(t, sin1, sin2, "w"); /* overwrite */
+ d = 0.9; /* default is one with smooth increment of gaps */
+ if (sscanf(sin3, "%lg", &d) < 1 && deltawritetimefirst < 2)
+ d = 0; /* default is zero for the first second */
+ if(d < 0)
+ flglockwrite += 2;
+ if (!flglockwrite) {
+ if (deltawritetime >= d) {
+ cmaes_WriteToFile(t, sin1, sin2);
+ flgwritten = 1;
+ } else if (d < 1
+ && t->gen-countiterlastwritten > maxdiffitertowrite) {
+ cmaes_WriteToFile(t, sin1, sin2);
+ flgwritten = 1;
+ }
+ }
+ break;
+ case 3 : /* check, checkeigen 1 or check eigen 1 */
+ if (strncmp(sin1, "eigen", 5) == 0) {
+ if (sscanf(sin2, " %lg", &d) == 1) {
+ if (d > 0)
+ t->flgCheckEigen = 1;
+ else
+ t->flgCheckEigen = 0;
+ }
+ else
+ t->flgCheckEigen = 0;
+ }
+ break;
+ case 4 : /* maxTimeFractionForEigendecompostion */
+ if (sscanf(sin1, " %lg", &d) == 1)
+ t->sp.updateCmode.maxtime = d;
+ break;
+ default :
+ break;
+ }
+ break; /* for ikey */
+ } /* if line contains keyword */
+ } /* for each keyword */
+ } /* while not EOF of signals.par */
+ if (t->writetime == 0)
+ t->firstwritetime = time(NULL);
+ if (t->printtime == 0)
+ t->firstprinttime = time(NULL);
+
+ if (flgprinted)
+ t->printtime = time(NULL);
+ if (flgwritten) {
+ t->writetime = time(NULL);
+ if (t->gen-countiterlastwritten > maxdiffitertowrite)
+ ++maxdiffitertowrite; /* smooth prolongation of writing gaps/intervals */
+ countiterlastwritten = (long int) t->gen;
+ }
+ --flglockprint;
+ --flglockwrite;
+ flglockprint = (flglockprint > 0) ? 1 : 0;
+ flglockwrite = (flglockwrite > 0) ? 1 : 0;
+} /* cmaes_ReadFromFilePtr */
+
+/* ========================================================= */
+static int
+Check_Eigen( int N, double **C, double *diag, double **Q)
+/*
+ exhaustive test of the output of the eigendecomposition
+ needs O(n^3) operations
+
+ writes to error file
+ returns number of detected inaccuracies
+*/
+{
+ /* compute Q diag Q^T and Q Q^T to check */
+ int i, j, k, res = 0;
+ double cc, dd;
+ static char s[324];
+
+ for (i=0; i < N; ++i)
+ for (j=0; j < N; ++j) {
+ for (cc=0.,dd=0., k=0; k < N; ++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) {
+ sprintf(s, "%d %d: %.17e %.17e, %e",
+ i, j, cc, C[i>j?i:j][i>j?j:i], cc-C[i>j?i:j][i>j?j:i]);
+ ERRORMESSAGE("cmaes_t:Eigen(): imprecise result detected ",
+ s, 0, 0);
+ ++res;
+ }
+ if (fabs(dd - (i==j)) > 1e-10) {
+ sprintf(s, "%d %d %.17e ", i, j, dd);
+ ERRORMESSAGE("cmaes_t:Eigen(): imprecise result detected (Q not orthog.)",
+ s, 0, 0);
+ ++res;
+ }
+ }
+ return res;
+}
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+void
+cmaes_UpdateEigensystem(cmaes_t *t, int flgforce)
+{
+ int i, N = t->sp.N;
+
+ cmaes_timings_update(&t->eigenTimings);
+
+ if(flgforce == 0) {
+ if (t->flgEigensysIsUptodate == 1)
+ return;
+
+ /* return on modulo generation number */
+ if (t->sp.updateCmode.flgalways == 0 /* not implemented, always ==0 */
+ && t->gen < t->genOfEigensysUpdate + t->sp.updateCmode.modulo
+ )
+ return;
+
+ /* return on time percentage */
+ if (t->sp.updateCmode.maxtime < 1.00
+ && t->eigenTimings.tictoctime > t->sp.updateCmode.maxtime * t->eigenTimings.totaltime
+ && t->eigenTimings.tictoctime > 0.0002)
+ return;
+ }
+ cmaes_timings_tic(&t->eigenTimings);
+
+ Eigen( N, t->C, t->rgD, t->B, t->rgdTmp);
+
+ cmaes_timings_toc(&t->eigenTimings);
+
+ /* find largest and smallest eigenvalue, they are supposed to be sorted anyway */
+ t->minEW = rgdouMin(t->rgD, N);
+ t->maxEW = rgdouMax(t->rgD, N);
+
+ if (t->flgCheckEigen)
+ /* needs O(n^3)! writes, in case, error message in error file */
+ i = Check_Eigen( N, t->C, t->rgD, t->B);
+
+#if 0
+ /* Limit Condition of C to dMaxSignifKond+1 */
+ if (t->maxEW > t->minEW * t->dMaxSignifKond) {
+ ERRORMESSAGE("Warning: Condition number of covariance matrix at upper limit.",
+ " Consider a rescaling or redesign of the objective function. " ,"","");
+ printf("\nWarning: Condition number of covariance matrix at upper limit\n");
+ tmp = t->maxEW/t->dMaxSignifKond - t->minEW;
+ tmp = t->maxEW/t->dMaxSignifKond;
+ t->minEW += tmp;
+ for (i=0;i<N;++i) {
+ t->C[i][i] += tmp;
+ t->rgD[i] += tmp;
+ }
+ } /* if */
+ t->dLastMinEWgroesserNull = minEW;
+#endif
+
+ for (i = 0; i < N; ++i)
+ t->rgD[i] = sqrt(t->rgD[i]);
+
+ t->flgEigensysIsUptodate = 1;
+ t->genOfEigensysUpdate = t->gen;
+
+ return;
+
+} /* cmaes_UpdateEigensystem() */
+
+
+/* ========================================================= */
+static void
+Eigen( int N, double **C, double *diag, double **Q, double *rgtmp)
+/*
+ Calculating eigenvalues and vectors.
+ Input:
+ N: dimension.
+ C: symmetric (1:N)xN-matrix, solely used to copy data to Q
+ niter: number of maximal iterations for QL-Algorithm.
+ rgtmp: N+1-dimensional vector for temporal use.
+ Output:
+ diag: N eigenvalues.
+ Q: Columns are normalized eigenvectors.
+ */
+{
+ int i, j;
+
+ if (rgtmp == NULL) /* was OK in former versions */
+ FATAL("cmaes_t:Eigen(): input parameter double *rgtmp must be non-NULL", 0,0,0);
+
+ /* copy C to Q */
+ if (C != Q) {
+ for (i=0; i < N; ++i)
+ for (j = 0; j <= i; ++j)
+ Q[i][j] = Q[j][i] = C[i][j];
+ }
+
+#if 0
+ Householder( N, Q, diag, rgtmp);
+ QLalgo( N, diag, Q, 30*N, rgtmp+1);
+#else
+ Householder2( N, Q, diag, rgtmp);
+ QLalgo2( N, diag, rgtmp, Q);
+#endif
+
+}
+
+
+/* ========================================================= */
+static void
+QLalgo2 (int n, double *d, double *e, double **V) {
+ /*
+ -> n : Dimension.
+ -> d : Diagonale of tridiagonal matrix.
+ -> e[1..n-1] : off-diagonal, output from Householder
+ -> V : matrix output von Householder
+ <- d : eigenvalues
+ <- e : garbage?
+ <- V : basis of eigenvectors, according to d
+
+ Symmetric tridiagonal QL algorithm, iterative
+ Computes the eigensystem from a tridiagonal matrix in roughtly 3N^3 operations
+
+ code adapted from Java JAMA package, function tql2.
+ */
+
+ 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) { /* TODO: check the case m == n, should be rejected here!? */
+ 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. */
+#if 1
+ /* TODO: really needed here? So far not, but practical and only O(n^2) */
+ {
+ 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;
+ }
+ }
+ }
+ }
+#endif
+} /* QLalgo2 */
+
+
+/* ========================================================= */
+static void
+Householder2(int n, double **V, double *d, double *e) {
+ /*
+ Householder transformation of a symmetric matrix V into tridiagonal form.
+ -> n : dimension
+ -> V : symmetric nxn-matrix
+ <- V : orthogonal transformation matrix:
+ tridiag matrix == V * V_in * V^t
+ <- d : diagonal
+ <- e[0..n-1] : off diagonal (elements 1..n-1)
+
+ code slightly adapted from the Java JAMA package, function private tred2()
+
+ */
+
+ 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() */
+
+
+#if 0
+/* ========================================================= */
+static void
+WriteMaxErrorInfo(cmaes_t *t)
+{
+ int i,j, N=t->sp.N;
+ char *s = (char *)new_void(200+30*(N+2), sizeof(char)); s[0] = '\0';
+
+ sprintf( s+strlen(s),"\nComplete Info\n");
+ sprintf( s+strlen(s)," Gen %20.12g\n", t->gen);
+ sprintf( s+strlen(s)," Dimension %d\n", N);
+ sprintf( s+strlen(s)," sigma %e\n", t->sigma);
+ sprintf( s+strlen(s)," lastminEW %e\n",
+ t->dLastMinEWgroesserNull);
+ sprintf( s+strlen(s)," maxKond %e\n\n", t->dMaxSignifKond);
+ sprintf( s+strlen(s)," x-vector rgD Basis...\n");
+ ERRORMESSAGE( s,0,0,0);
+ s[0] = '\0';
+ for (i = 0; i < N; ++i)
+ {
+ sprintf( s+strlen(s), " %20.12e", t->rgxmean[i]);
+ sprintf( s+strlen(s), " %10.4e", t->rgD[i]);
+ for (j = 0; j < N; ++j)
+ sprintf( s+strlen(s), " %10.2e", t->B[i][j]);
+ ERRORMESSAGE( s,0,0,0);
+ s[0] = '\0';
+ }
+ ERRORMESSAGE( "\n",0,0,0);
+ free( s);
+} /* WriteMaxErrorInfo() */
+#endif
+
+/* --------------------------------------------------------- */
+/* --------------- Functions: cmaes_timings_t -------------- */
+/* --------------------------------------------------------- */
+/* cmaes_timings_t measures overall time and times between calls
+ * of tic and toc. For small time spans (up to 1000 seconds)
+ * CPU time via clock() is used. For large time spans the
+ * fall-back to elapsed time from time() is used.
+ * cmaes_timings_update() must be called often enough to prevent
+ * the fallback. */
+/* --------------------------------------------------------- */
+void
+cmaes_timings_init(cmaes_timings_t *t) {
+ t->totaltotaltime = 0;
+ cmaes_timings_start(t);
+}
+void
+cmaes_timings_start(cmaes_timings_t *t) {
+ t->totaltime = 0;
+ t->tictoctime = 0;
+ t->lasttictoctime = 0;
+ t->istic = 0;
+ t->lastclock = clock();
+ t->lasttime = time(NULL);
+ t->lastdiff = 0;
+ t->tictoczwischensumme = 0;
+ t->isstarted = 1;
+}
+
+double
+cmaes_timings_update(cmaes_timings_t *t) {
+/* returns time between last call of cmaes_timings_*() and now,
+ * should better return totaltime or tictoctime?
+ */
+ double diffc, difft;
+ clock_t lc = t->lastclock; /* measure CPU in 1e-6s */
+ time_t lt = t->lasttime; /* measure time in s */
+
+ if (t->isstarted != 1)
+ FATAL("cmaes_timings_started() must be called before using timings... functions",0,0,0);
+
+ t->lastclock = clock(); /* measures at most 2147 seconds, where 1s = 1e6 CLOCKS_PER_SEC */
+ t->lasttime = time(NULL);
+
+ diffc = (double)(t->lastclock - lc) / CLOCKS_PER_SEC; /* is presumably in [-21??, 21??] */
+ difft = difftime(t->lasttime, lt); /* is presumably an integer */
+
+ t->lastdiff = difft; /* on the "save" side */
+
+ /* use diffc clock measurement if appropriate */
+ if (diffc > 0 && difft < 1000)
+ t->lastdiff = diffc;
+
+ if (t->lastdiff < 0)
+ FATAL("BUG in time measurement", 0, 0, 0);
+
+ t->totaltime += t->lastdiff;
+ t->totaltotaltime += t->lastdiff;
+ if (t->istic) {
+ t->tictoczwischensumme += t->lastdiff;
+ t->tictoctime += t->lastdiff;
+ }
+
+ return t->lastdiff;
+}
+
+void
+cmaes_timings_tic(cmaes_timings_t *t) {
+ if (t->istic) { /* message not necessary ? */
+ ERRORMESSAGE("Warning: cmaes_timings_tic called twice without toc",0,0,0);
+ return;
+ }
+ cmaes_timings_update(t);
+ t->istic = 1;
+}
+
+double
+cmaes_timings_toc(cmaes_timings_t *t) {
+ if (!t->istic) {
+ ERRORMESSAGE("Warning: cmaes_timings_toc called without tic",0,0,0);
+ return -1;
+ }
+ cmaes_timings_update(t);
+ t->lasttictoctime = t->tictoczwischensumme;
+ t->tictoczwischensumme = 0;
+ t->istic = 0;
+ return t->lasttictoctime;
+}
+
+/* --------------------------------------------------------- */
+/* ---------------- Functions: cmaes_random_t -------------- */
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+/* X_1 exakt : 0.79788456) */
+/* chi_eins simuliert : 0.798xx (seed -3) */
+/* +-0.001 */
+/* --------------------------------------------------------- */
+/*
+ Gauss() liefert normalverteilte Zufallszahlen
+ bei vorgegebenem seed.
+*/
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+
+long
+cmaes_random_init( cmaes_random_t *t, long unsigned inseed)
+{
+ clock_t cloc = clock();
+
+ t->flgstored = 0;
+ t->rgrand = (long *) new_void(32, sizeof(long));
+ if (inseed < 1) {
+ while ((long) (cloc - clock()) == 0)
+ ; /* TODO: remove this for time critical applications? */
+ inseed = (long unsigned)labs((long)(100*time(NULL)+clock()));
+ }
+ return cmaes_random_Start(t, inseed);
+}
+
+void
+cmaes_random_exit(cmaes_random_t *t)
+{
+ free( t->rgrand);
+}
+
+/* --------------------------------------------------------- */
+long cmaes_random_Start( cmaes_random_t *t, long unsigned inseed)
+{
+ long tmp;
+ int i;
+
+ t->flgstored = 0;
+ t->startseed = inseed; /* purely for bookkeeping */
+ while (inseed > 2e9)
+ inseed /= 2; /* prevent infinite loop on 32 bit system */
+ 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->rgrand[i] = t->aktseed;
+ }
+ t->aktrand = t->rgrand[0];
+ return inseed;
+}
+
+/* --------------------------------------------------------- */
+double cmaes_random_Gauss(cmaes_random_t *t)
+{
+ double x1, x2, rquad, fac;
+
+ if (t->flgstored)
+ {
+ t->flgstored = 0;
+ return t->hold;
+ }
+ do
+ {
+ x1 = 2.0 * cmaes_random_Uniform(t) - 1.0;
+ x2 = 2.0 * cmaes_random_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;
+}
+
+/* --------------------------------------------------------- */
+double cmaes_random_Uniform( cmaes_random_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->aktrand / 67108865;
+ t->aktrand = t->rgrand[tmp];
+ t->rgrand[tmp] = t->aktseed;
+ return (double)(t->aktrand)/(2.147483647e9);
+}
+
+static char *
+szCat(const char *sz1, const char*sz2,
+ const char *sz3, const char *sz4);
+
+/* --------------------------------------------------------- */
+/* -------------- Functions: cmaes_readpara_t -------------- */
+/* --------------------------------------------------------- */
+void
+cmaes_readpara_init (cmaes_readpara_t *t,
+ int dim,
+ const double * inxstart,
+ const double * inrgsigma,
+ int inseed,
+ int lambda,
+ const char * filename)
+{
+ int i, N;
+ /* TODO: make sure cmaes_readpara_init has not been called already */
+ t->filename = NULL; /* set after successful Read */
+ t->rgsformat = (const char **) new_void(55, sizeof(char *));
+ t->rgpadr = (void **) new_void(55, sizeof(void *));
+ t->rgskeyar = (const char **) new_void(11, sizeof(char *));
+ t->rgp2adr = (double ***) new_void(11, sizeof(double **));
+ t->weigkey = (char *)new_void(7, sizeof(char));
+
+ /* All scalars: */
+ i = 0;
+ t->rgsformat[i] = " N %d"; t->rgpadr[i++] = (void *) &t->N;
+ t->rgsformat[i] = " seed %d"; t->rgpadr[i++] = (void *) &t->seed;
+ t->rgsformat[i] = " stopMaxFunEvals %lg"; t->rgpadr[i++] = (void *) &t->stopMaxFunEvals;
+ t->rgsformat[i] = " stopMaxIter %lg"; t->rgpadr[i++] = (void *) &t->stopMaxIter;
+ t->rgsformat[i] = " stopFitness %lg"; t->rgpadr[i++]=(void *) &t->stStopFitness.val;
+ t->rgsformat[i] = " stopTolFun %lg"; t->rgpadr[i++]=(void *) &t->stopTolFun;
+ t->rgsformat[i] = " stopTolFunHist %lg"; t->rgpadr[i++]=(void *) &t->stopTolFunHist;
+ t->rgsformat[i] = " stopTolX %lg"; t->rgpadr[i++]=(void *) &t->stopTolX;
+ t->rgsformat[i] = " stopTolUpXFactor %lg"; t->rgpadr[i++]=(void *) &t->stopTolUpXFactor;
+ t->rgsformat[i] = " lambda %d"; t->rgpadr[i++] = (void *) &t->lambda;
+ t->rgsformat[i] = " mu %d"; t->rgpadr[i++] = (void *) &t->mu;
+ t->rgsformat[i] = " weights %5s"; t->rgpadr[i++] = (void *) t->weigkey;
+ t->rgsformat[i] = " fac*cs %lg";t->rgpadr[i++] = (void *) &t->cs;
+ t->rgsformat[i] = " fac*damps %lg"; t->rgpadr[i++] = (void *) &t->damps;
+ t->rgsformat[i] = " ccumcov %lg"; t->rgpadr[i++] = (void *) &t->ccumcov;
+ t->rgsformat[i] = " mucov %lg"; t->rgpadr[i++] = (void *) &t->mucov;
+ t->rgsformat[i] = " fac*ccov %lg"; t->rgpadr[i++]=(void *) &t->ccov;
+ t->rgsformat[i] = " diagonalCovarianceMatrix %lg"; t->rgpadr[i++]=(void *) &t->diagonalCov;
+ t->rgsformat[i] = " updatecov %lg"; t->rgpadr[i++]=(void *) &t->updateCmode.modulo;
+ t->rgsformat[i] = " maxTimeFractionForEigendecompostion %lg"; t->rgpadr[i++]=(void *) &t->updateCmode.maxtime;
+ t->rgsformat[i] = " resume %59s"; t->rgpadr[i++] = (void *) t->resumefile;
+ t->rgsformat[i] = " fac*maxFunEvals %lg"; t->rgpadr[i++] = (void *) &t->facmaxeval;
+ t->rgsformat[i] = " fac*updatecov %lg"; t->rgpadr[i++]=(void *) &t->facupdateCmode;
+ t->n1para = i;
+ t->n1outpara = i-2; /* disregard last parameters in WriteToFile() */
+
+ /* arrays */
+ i = 0;
+ t->rgskeyar[i] = " typicalX %d"; t->rgp2adr[i++] = &t->typicalX;
+ t->rgskeyar[i] = " initialX %d"; t->rgp2adr[i++] = &t->xstart;
+ t->rgskeyar[i] = " initialStandardDeviations %d"; t->rgp2adr[i++] = &t->rgInitialStds;
+ t->rgskeyar[i] = " diffMinChange %d"; t->rgp2adr[i++] = &t->rgDiffMinChange;
+ t->n2para = i;
+
+ t->N = dim;
+ t->seed = (unsigned) inseed;
+ t->xstart = NULL;
+ t->typicalX = NULL;
+ t->typicalXcase = 0;
+ t->rgInitialStds = NULL;
+ t->rgDiffMinChange = NULL;
+ t->stopMaxFunEvals = -1;
+ t->stopMaxIter = -1;
+ t->facmaxeval = 1;
+ t->stStopFitness.flg = -1;
+ t->stopTolFun = 1e-12;
+ t->stopTolFunHist = 1e-13;
+ t->stopTolX = 0; /* 1e-11*insigma would also be reasonable */
+ t->stopTolUpXFactor = 1e3;
+
+ t->lambda = lambda;
+ t->mu = -1;
+ t->mucov = -1;
+ t->weights = NULL;
+ strcpy(t->weigkey, "log");
+
+ t->cs = -1;
+ t->ccumcov = -1;
+ t->damps = -1;
+ t->ccov = -1;
+
+ t->diagonalCov = 0; /* default is 0, but this might change in future, see below */
+
+ t->updateCmode.modulo = -1;
+ t->updateCmode.maxtime = -1;
+ t->updateCmode.flgalways = 0;
+ t->facupdateCmode = 1;
+ strcpy(t->resumefile, "_no_");
+
+ /* filename == NULL invokes default in cmaes_readpara_Read... */
+ if (!isNoneStr(filename) && (!filename || strcmp(filename, "writeonly") != 0))
+ cmaes_readpara_ReadFromFile(t, filename);
+
+ if (t->N <= 0)
+ t->N = dim;
+
+ N = t->N;
+ if (N == 0)
+ FATAL("cmaes_readpara_t(): problem dimension N undefined.\n",
+ " (no default value available).",0,0);
+ if (t->xstart == NULL && inxstart == NULL && t->typicalX == NULL) {
+ ERRORMESSAGE("Error: initialX undefined. typicalX = 0.5...0.5 used.","","","");
+ printf("\nError: initialX undefined. typicalX = 0.5...0.5 used.\n");
+ }
+ if (t->rgInitialStds == NULL && inrgsigma == NULL) {
+ /* FATAL("initialStandardDeviations undefined","","",""); */
+ ERRORMESSAGE("Error: initialStandardDeviations undefined. 0.3...0.3 used.","","","");
+ printf("\nError: initialStandardDeviations undefined. 0.3...0.3 used.\n");
+ }
+
+ if (t->xstart == NULL) {
+ t->xstart = new_double(N);
+
+ /* put inxstart into xstart */
+ if (inxstart != NULL) {
+ for (i=0; i<N; ++i)
+ t->xstart[i] = inxstart[i];
+ }
+ /* otherwise use typicalX or default */
+ else {
+ t->typicalXcase = 1;
+ for (i=0; i<N; ++i)
+ t->xstart[i] = (t->typicalX == NULL) ? 0.5 : t->typicalX[i];
+ }
+ } /* xstart == NULL */
+
+ if (t->rgInitialStds == NULL) {
+ t->rgInitialStds = new_double(N);
+ for (i=0; i<N; ++i)
+ t->rgInitialStds[i] = (inrgsigma == NULL) ? 0.3 : inrgsigma[i];
+ }
+
+ t->flgsupplemented = 0;
+
+} /* cmaes_readpara_init */
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+void cmaes_readpara_exit(cmaes_readpara_t *t)
+{
+ if (t->filename != NULL)
+ free( t->filename);
+ if (t->xstart != NULL) /* not really necessary */
+ free( t->xstart);
+ if (t->typicalX != NULL)
+ free( t->typicalX);
+ if (t->rgInitialStds != NULL)
+ free( t->rgInitialStds);
+ if (t->rgDiffMinChange != NULL)
+ free( t->rgDiffMinChange);
+ if (t->weights != NULL)
+ free( t->weights);
+
+ free((void*)t->rgsformat);
+ free(t->rgpadr);
+ free((void*)t->rgskeyar);
+ free(t->rgp2adr);
+ free(t->weigkey);
+}
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+void
+cmaes_readpara_ReadFromFile(cmaes_readpara_t *t, const char * filename)
+{
+ char s[1000];
+ const char *ss = "cmaes_initials.par";
+ int ipara, i;
+ int size;
+ FILE *fp;
+ if (filename == NULL) {
+ filename = ss;
+ }
+ t->filename = NULL; /* nothing read so far */
+ fp = fopen( filename, "r");
+ if (fp == NULL) {
+ ERRORMESSAGE("cmaes_ReadFromFile(): could not open '", filename, "'",0);
+ return;
+ }
+ for (ipara=0; ipara < t->n1para; ++ipara)
+ {
+ rewind(fp);
+ while(fgets(s, sizeof(s), fp) != NULL)
+ { /* skip comments */
+ if (s[0] == '#' || s[0] == '%')
+ continue;
+ if(sscanf(s, t->rgsformat[ipara], t->rgpadr[ipara]) == 1) {
+ if (strncmp(t->rgsformat[ipara], " stopFitness ", 13) == 0)
+ t->stStopFitness.flg = 1;
+ break;
+ }
+ }
+ } /* for */
+ if (t->N <= 0)
+ FATAL("cmaes_readpara_ReadFromFile(): No valid dimension N",0,0,0);
+ for (ipara=0; ipara < t->n2para; ++ipara)
+ {
+ rewind(fp);
+ while(fgets(s, sizeof(s), fp) != NULL) /* read one line */
+ { /* skip comments */
+ if (s[0] == '#' || s[0] == '%')
+ continue;
+ if(sscanf(s, t->rgskeyar[ipara], &size) == 1) { /* size==number of values to be read */
+ if (size > 0) {
+ *t->rgp2adr[ipara] = new_double(t->N);
+ for (i=0;i<size&&i<t->N;++i) /* start reading next line */
+ if (fscanf(fp, " %lf", &(*t->rgp2adr[ipara])[i]) != 1)
+ break;
+ if (i<size && i < t->N) {
+ ERRORMESSAGE("cmaes_readpara_ReadFromFile ", filename, ": ",0);
+ FATAL( "'", t->rgskeyar[ipara],
+ "' not enough values found.\n",
+ " Remove all comments between numbers.");
+ }
+ for (; i < t->N; ++i) /* recycle */
+ (*t->rgp2adr[ipara])[i] = (*t->rgp2adr[ipara])[i%size];
+ }
+ }
+ }
+ } /* for */
+ fclose(fp);
+ assign_string(&(t->filename), filename); /* t->filename must be freed */
+ return;
+} /* cmaes_readpara_ReadFromFile() */
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+void
+cmaes_readpara_WriteToFile(cmaes_readpara_t *t, const char *filenamedest)
+{
+ int ipara, i;
+ size_t len;
+ time_t ti = time(NULL);
+ FILE *fp = fopen( filenamedest, "a");
+ if(fp == NULL) {
+ ERRORMESSAGE("cmaes_WriteToFile(): could not open '",
+ filenamedest, "'",0);
+ return;
+ }
+ fprintf(fp, "\n# Read from %s at %s\n", t->filename ? t->filename : "",
+ asctime(localtime(&ti))); /* == ctime() */
+ for (ipara=0; ipara < 1; ++ipara) {
+ fprintf(fp, t->rgsformat[ipara], *(int *)t->rgpadr[ipara]);
+ fprintf(fp, "\n");
+ }
+ for (ipara=0; ipara < t->n2para; ++ipara) {
+ if(*t->rgp2adr[ipara] == NULL)
+ continue;
+ fprintf(fp, t->rgskeyar[ipara], t->N);
+ fprintf(fp, "\n");
+ for (i=0; i<t->N; ++i)
+ fprintf(fp, "%7.3g%c", (*t->rgp2adr[ipara])[i], (i%5==4)?'\n':' ');
+ fprintf(fp, "\n");
+ }
+ for (ipara=1; ipara < t->n1outpara; ++ipara) {
+ if (strncmp(t->rgsformat[ipara], " stopFitness ", 13) == 0)
+ if(t->stStopFitness.flg == 0) {
+ fprintf(fp, " stopFitness\n");
+ continue;
+ }
+ len = strlen(t->rgsformat[ipara]);
+ if (t->rgsformat[ipara][len-1] == 'd') /* read integer */
+ fprintf(fp, t->rgsformat[ipara], *(int *)t->rgpadr[ipara]);
+ else if (t->rgsformat[ipara][len-1] == 's') /* read string */
+ fprintf(fp, t->rgsformat[ipara], (char *)t->rgpadr[ipara]);
+ else {
+ if (strncmp(" fac*", t->rgsformat[ipara], 5) == 0) {
+ fprintf(fp, " ");
+ fprintf(fp, t->rgsformat[ipara]+5, *(double *)t->rgpadr[ipara]);
+ } else
+ fprintf(fp, t->rgsformat[ipara], *(double *)t->rgpadr[ipara]);
+ }
+ fprintf(fp, "\n");
+ } /* for */
+ fprintf(fp, "\n");
+ fclose(fp);
+} /* cmaes_readpara_WriteToFile() */
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+void
+cmaes_readpara_SupplementDefaults(cmaes_readpara_t *t)
+/* Called (only) once to finally set parameters. The settings
+ * typically depend on the current parameter values itself,
+ * where 0 or -1 may indicate to set them to a certain default
+ * value. For this reason calling `SupplementDefaults` twice
+ * might lead to unexpected results.
+*/
+{
+ double t1, t2;
+ int N = t->N;
+ clock_t cloc = clock();
+
+ if (t->flgsupplemented)
+ FATAL("cmaes_readpara_SupplementDefaults() cannot be called twice.",0,0,0);
+ if (t->seed < 1) {
+ while ((int) (cloc - clock()) == 0)
+ ; /* TODO: remove this for time critical applications!? */
+ t->seed = (unsigned int)labs((long)(100*time(NULL)+clock()));
+ }
+
+ if (t->stStopFitness.flg == -1)
+ t->stStopFitness.flg = 0;
+
+ if (t->lambda < 2)
+ t->lambda = 4+(int)(3*log((double)N));
+ if (t->mu == -1) {
+ t->mu = t->lambda/2;
+ cmaes_readpara_SetWeights(t, t->weigkey);
+ }
+ if (t->weights == NULL)
+ cmaes_readpara_SetWeights(t, t->weigkey);
+
+ if (t->cs > 0) /* factor was read */
+ t->cs *= (t->mueff + 2.) / (N + t->mueff + 3.);
+ if (t->cs <= 0 || t->cs >= 1)
+ t->cs = (t->mueff + 2.) / (N + t->mueff + 3.);
+
+ if (t->ccumcov <= 0 || t->ccumcov > 1)
+ t->ccumcov = 4. / (N + 4);
+
+ if (t->mucov < 1) {
+ t->mucov = t->mueff;
+ }
+ t1 = 2. / ((N+1.4142)*(N+1.4142));
+ t2 = (2.*t->mueff-1.) / ((N+2.)*(N+2.)+t->mueff);
+ t2 = (t2 > 1) ? 1 : t2;
+ t2 = (1./t->mucov) * t1 + (1.-1./t->mucov) * t2;
+ if (t->ccov >= 0) /* ccov holds the read factor */
+ t->ccov *= t2;
+ if (t->ccov < 0 || t->ccov > 1) /* set default in case */
+ t->ccov = t2;
+
+ if (t->diagonalCov == -1)
+ t->diagonalCov = 2 + 100. * N / sqrt((double)t->lambda);
+
+ if (t->stopMaxFunEvals == -1) /* may depend on ccov in near future */
+ t->stopMaxFunEvals = t->facmaxeval*900*(N+3)*(N+3);
+ else
+ t->stopMaxFunEvals *= t->facmaxeval;
+
+ if (t->stopMaxIter == -1)
+ t->stopMaxIter = ceil((double)(t->stopMaxFunEvals / t->lambda));
+
+ if (t->damps < 0)
+ t->damps = 1; /* otherwise a factor was read */
+ t->damps = t->damps
+ * (1 + 2*douMax(0., sqrt((t->mueff-1.)/(N+1.)) - 1)) /* basic factor */
+ * douMax(0.3, 1. - /* modify for short runs */
+ (double)N / (1e-6+douMin(t->stopMaxIter, t->stopMaxFunEvals/t->lambda)))
+ + t->cs; /* minor increment */
+
+ if (t->updateCmode.modulo < 0)
+ t->updateCmode.modulo = 1./t->ccov/(double)(N)/10.;
+ t->updateCmode.modulo *= t->facupdateCmode;
+ if (t->updateCmode.maxtime < 0)
+ t->updateCmode.maxtime = 0.20; /* maximal 20% of CPU-time */
+
+ t->flgsupplemented = 1;
+
+} /* cmaes_readpara_SupplementDefaults() */
+
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+void
+cmaes_readpara_SetWeights(cmaes_readpara_t *t, const char * mode)
+{
+ double s1, s2;
+ int i;
+
+ if(t->weights != NULL)
+ free( t->weights);
+ t->weights = new_double(t->mu);
+ if (strcmp(mode, "lin") == 0)
+ for (i=0; i<t->mu; ++i)
+ t->weights[i] = t->mu - i;
+ else if (strncmp(mode, "equal", 3) == 0)
+ for (i=0; i<t->mu; ++i)
+ t->weights[i] = 1;
+ else if (strcmp(mode, "log") == 0)
+ for (i=0; i<t->mu; ++i)
+ t->weights[i] = log(t->mu+1.)-log(i+1.);
+ else
+ for (i=0; i<t->mu; ++i)
+ t->weights[i] = log(t->mu+1.)-log(i+1.);
+
+ /* normalize weights vector and set mueff */
+ for (i=0, s1=0, s2=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]))
+ FATAL("cmaes_readpara_SetWeights(): invalid setting of mu or lambda",0,0,0);
+
+} /* cmaes_readpara_SetWeights() */
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+static int
+isNoneStr(const char * filename)
+{
+ if (filename && (strcmp(filename, "no") == 0
+ || strcmp(filename, "non") == 0
+ || strcmp(filename, "none") == 0))
+ return 1;
+
+ return 0;
+}
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+static double
+douSquare(double d)
+{
+ return d*d;
+}
+static int
+intMin( int i, int j)
+{
+ return i < j ? i : j;
+}
+static double
+douMax( double i, double j)
+{
+ return i > j ? i : j;
+}
+static double
+douMin( double i, double j)
+{
+ return i < j ? i : j;
+}
+static double
+rgdouMax( const double *rgd, int len)
+{
+ int i;
+ double max = rgd[0];
+ for (i = 1; i < len; ++i)
+ max = (max < rgd[i]) ? rgd[i] : max;
+ return max;
+}
+
+static double
+rgdouMin( const double *rgd, int len)
+{
+ int i;
+ double min = rgd[0];
+ for (i = 1; i < len; ++i)
+ min = (min > rgd[i]) ? rgd[i] : min;
+ return min;
+}
+
+static int
+MaxIdx( const double *rgd, int len)
+{
+ int i, res;
+ for(i=1, res=0; i<len; ++i)
+ if(rgd[i] > rgd[res])
+ res = i;
+ return res;
+}
+static int
+MinIdx( const double *rgd, int len)
+{
+ int i, res;
+ for(i=1, res=0; i<len; ++i)
+ if(rgd[i] < rgd[res])
+ res = i;
+ return res;
+}
+
+static 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;
+}
+
+static int SignOfDiff(const void *d1, const void * d2)
+{
+ return *((double *) d1) > *((double *) d2) ? 1 : -1;
+}
+
+#if 1
+/* dirty index sort */
+static void Sorted_index(const double *rgFunVal, int *iindex, int n)
+{
+ int i, j;
+ for (i=1, iindex[0]=0; i<n; ++i) {
+ for (j=i; j>0; --j) {
+ if (rgFunVal[iindex[j-1]] < rgFunVal[i])
+ break;
+ iindex[j] = iindex[j-1]; /* shift up */
+ }
+ iindex[j] = i; /* insert i */
+ }
+}
+#endif
+
+static void * new_void(int n, size_t size)
+{
+ static char s[70];
+ void *p = calloc((unsigned) n, size);
+ if (p == NULL) {
+ sprintf(s, "new_void(): calloc(%ld,%ld) failed",(long)n,(long)size);
+ FATAL(s,0,0,0);
+ }
+ return p;
+}
+
+double *
+cmaes_NewDouble(int n)
+{
+ return new_double(n);
+}
+
+static double * new_double(int n)
+{
+ static char s[170];
+ double *p = (double *) calloc((unsigned) n, sizeof(double));
+ if (p == NULL) {
+ sprintf(s, "new_double(): calloc(%ld,%ld) failed",
+ (long)n,(long)sizeof(double));
+ FATAL(s,0,0,0);
+ }
+ return p;
+}
+
+static char * new_string(const char *ins)
+{
+ static char s[170];
+ unsigned i;
+ char *p;
+ unsigned len = (unsigned) ((ins != NULL) ? strlen(ins) : 0);
+
+ if (len > 1000) {
+ FATAL("new_string(): input string length was larger then 1000 ",
+ "(possibly due to uninitialized char *filename)",0,0);
+ }
+
+ p = (char *) calloc( len + 1, sizeof(char));
+ if (p == NULL) {
+ sprintf(s, "new_string(): calloc(%ld,%ld) failed",
+ (long)len,(long)sizeof(char));
+ FATAL(s,0,0,0);
+ }
+ for (i = 0; i < len; ++i)
+ p[i] = ins[i];
+ return p;
+}
+static void assign_string(char ** pdests, const char *ins)
+{
+ if (*pdests)
+ free(*pdests);
+ if (ins == NULL)
+ *pdests = NULL;
+ else
+ *pdests = new_string(ins);
+}
+
+/* --------------------------------------------------------- */
+/* --------------------------------------------------------- */
+
+/* ========================================================= */
+void
+cmaes_FATAL(char const *s1, char const *s2, char const *s3,
+ char const *s4)
+{
+ time_t t = time(NULL);
+ ERRORMESSAGE( s1, s2, s3, s4);
+ ERRORMESSAGE("*** Exiting cmaes_t ***",0,0,0);
+ printf("\n -- %s %s\n", asctime(localtime(&t)),
+ s2 ? szCat(s1, s2, s3, s4) : s1);
+ printf(" *** CMA-ES ABORTED, see errcmaes.err *** \n");
+ fflush(stdout);
+ exit(1);
+}
+
+/* ========================================================= */
+static void
+FATAL(char const *s1, char const *s2, char const *s3,
+ char const *s4)
+{
+ cmaes_FATAL(s1, s2, s3, s4);
+}
+
+/* ========================================================= */
+static void ERRORMESSAGE( char const *s1, char const *s2,
+ char const *s3, char const *s4)
+{
+#if 1
+ /* static char szBuf[700]; desirable but needs additional input argument
+ sprintf(szBuf, "%f:%f", gen, gen*lambda);
+ */
+ time_t t = time(NULL);
+ FILE *fp = fopen( "errcmaes.err", "a");
+ if (!fp)
+ {
+ printf("\nFATAL ERROR: %s\n", s2 ? szCat(s1, s2, s3, s4) : s1);
+ printf("cmaes_t could not open file 'errcmaes.err'.");
+ printf("\n *** CMA-ES ABORTED *** ");
+ fflush(stdout);
+ exit(1);
+ }
+ fprintf( fp, "\n -- %s %s\n", asctime(localtime(&t)),
+ s2 ? szCat(s1, s2, s3, s4) : s1);
+ fclose (fp);
+#endif
+}
+
+/* ========================================================= */
+static char *szCat(const char *sz1, const char*sz2,
+ const char *sz3, const char *sz4)
+{
+ static char szBuf[700];
+
+ if (!sz1)
+ FATAL("szCat() : Invalid Arguments",0,0,0);
+
+ strncpy ((char *)szBuf, sz1, (unsigned)intMin( (int)strlen(sz1), 698));
+ szBuf[intMin( (int)strlen(sz1), 698)] = '\0';
+ if (sz2)
+ strncat ((char *)szBuf, sz2,
+ (unsigned)intMin((int)strlen(sz2)+1, 698 - (int)strlen((char const *)szBuf)));
+ if (sz3)
+ strncat((char *)szBuf, sz3,
+ (unsigned)intMin((int)strlen(sz3)+1, 698 - (int)strlen((char const *)szBuf)));
+ if (sz4)
+ strncat((char *)szBuf, sz4,
+ (unsigned)intMin((int)strlen(sz4)+1, 698 - (int)strlen((char const *)szBuf)));
+ return (char *) szBuf;
+}
+
+
diff --git a/external/cma-es/cmaes.h b/external/cma-es/cmaes.h
new file mode 100644
index 0000000..277accb
--- /dev/null
+++ b/external/cma-es/cmaes.h
@@ -0,0 +1,176 @@
+/* --------------------------------------------------------- */
+/* --- 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
+/* cmaes_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;
+} cmaes_random_t;
+
+typedef struct
+/* cmaes_timings_t
+ * time measurement, used to time eigendecomposition
+ */
+{
+ /* for outside use */
+ double totaltime; /* zeroed by calling re-calling cmaes_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;
+} cmaes_timings_t;
+
+typedef struct
+/* cmaes_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 */
+ short flgsupplemented;
+
+ /* 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;
+} cmaes_readpara_t;
+
+typedef struct
+/* cmaes_t
+ * CMA-ES "object"
+ */
+{
+ const char *version;
+ /* char *signalsFilename; */
+ cmaes_readpara_t sp;
+ cmaes_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;
+ cmaes_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
diff --git a/external/cma-es/cmaes_interface.h b/external/cma-es/cmaes_interface.h
new file mode 100644
index 0000000..be7c736
--- /dev/null
+++ b/external/cma-es/cmaes_interface.h
@@ -0,0 +1,68 @@
+/* --------------------------------------------------------- */
+/* --- 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
+
+ Documentation: see file docfunctions.txt
+
+ License: see file cmaes.c
+*/
+#include "cmaes.h"
+
+/* --------------------------------------------------------- */
+/* ------------------ Interface ---------------------------- */
+/* --------------------------------------------------------- */
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/* --- 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_init_para(cmaes_t *, int dimension , double *xstart,
+ double *stddev, long seed, int lambda,
+ const char *input_parameter_filename);
+double * cmaes_init_final(cmaes_t *);
+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
+} // end extern "C"
+#endif
+
--
GitLab