From 0a21943fc50eca1da70f7670052a439d24589cb2 Mon Sep 17 00:00:00 2001
From: amirabdelrahman <amira-rahman@aucegypt.edu>
Date: Mon, 27 Apr 2020 03:47:31 -0400
Subject: [PATCH] added dnn
---
CMakeLists.txt | 3 +
cpp_dnn/CMakeLists.txt | 24 +
cpp_dnn/cmaes.cpp | 3137 +++++++++++++++++++++++++++++++++++++
cpp_dnn/cmaes.h | 176 +++
cpp_dnn/cmaes_interface.h | 68 +
cpp_dnn/dnn.cpp | 483 ++++++
cpp_dnn/frep.cpp | 707 +++++++++
cpp_dnn/image.cpp | 143 ++
cpp_dnn/main.cpp | 185 +++
cpp_dnn/nelder_mead.h | 276 ++++
cpp_dnn/vector.h | 22 +
cpp_frep/CMakeLists.txt | 6 +-
cpp_frep/main.cpp | 4 +-
13 files changed, 5230 insertions(+), 4 deletions(-)
create mode 100644 cpp_dnn/CMakeLists.txt
create mode 100644 cpp_dnn/cmaes.cpp
create mode 100644 cpp_dnn/cmaes.h
create mode 100644 cpp_dnn/cmaes_interface.h
create mode 100644 cpp_dnn/dnn.cpp
create mode 100644 cpp_dnn/frep.cpp
create mode 100644 cpp_dnn/image.cpp
create mode 100644 cpp_dnn/main.cpp
create mode 100644 cpp_dnn/nelder_mead.h
create mode 100644 cpp_dnn/vector.h
diff --git a/CMakeLists.txt b/CMakeLists.txt
index a2e3b0e..665ec8e 100755
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -7,7 +7,10 @@ endif()
message(STATUS "Build type: ${CMAKE_BUILD_TYPE}")
find_package (Eigen3 3.3 REQUIRED NO_MODULE)
+find_package(Armadillo REQUIRED)
+find_package( OpenCV REQUIRED )
include(cmake/shared_settings.cmake)
add_subdirectory(cpp_frep)
+add_subdirectory(cpp_dnn)
diff --git a/cpp_dnn/CMakeLists.txt b/cpp_dnn/CMakeLists.txt
new file mode 100644
index 0000000..e7435ef
--- /dev/null
+++ b/cpp_dnn/CMakeLists.txt
@@ -0,0 +1,24 @@
+add_library(cpp_dnn_lib
+ dnn.cpp
+ vector.h
+ frep.cpp
+ nelder_mead.h
+ image.cpp
+ cmaes.cpp
+ cmaes.h
+ cmaes_interface.h
+)
+
+target_include_directories(cpp_dnn_lib PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${ARMADILLO_INCLUDE_DIRS})
+target_link_libraries(cpp_dnn_lib shared_settings Eigen3::Eigen ${ARMADILLO_LIBRARIES} ${OpenCV_LIBS})
+target_compile_features(cpp_dnn_lib PUBLIC cxx_std_17)
+
+add_executable(cpp_dnn
+ main.cpp
+)
+target_link_libraries(cpp_dnn shared_settings cpp_dnn_lib)
+
+# CMA-ES (c-cmaes)
+# From https://github.com/CMA-ES/c-cmaes.
+# add_subdirectory(cma-es)
+
diff --git a/cpp_dnn/cmaes.cpp b/cpp_dnn/cmaes.cpp
new file mode 100644
index 0000000..a38e1d4
--- /dev/null
+++ b/cpp_dnn/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/cpp_dnn/cmaes.h b/cpp_dnn/cmaes.h
new file mode 100644
index 0000000..277accb
--- /dev/null
+++ b/cpp_dnn/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/cpp_dnn/cmaes_interface.h b/cpp_dnn/cmaes_interface.h
new file mode 100644
index 0000000..be7c736
--- /dev/null
+++ b/cpp_dnn/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
+
diff --git a/cpp_dnn/dnn.cpp b/cpp_dnn/dnn.cpp
new file mode 100644
index 0000000..c0c3da6
--- /dev/null
+++ b/cpp_dnn/dnn.cpp
@@ -0,0 +1,483 @@
+#include <Eigen/Core>
+#include <armadillo>
+#include <iostream>
+#include <vector>
+
+namespace cpp_dnn {
+ using scalar=double;
+ using MatrixXI = arma::Mat<long long unsigned int>;
+ using VectorXf = Eigen::Matrix<scalar, Eigen::Dynamic, 1>;
+ using MatrixX = arma::Mat<scalar>;
+
+ //..................................................................................................
+
+ //rename to layer??
+ class Node {
+ public:
+ int _m;
+ int _n;
+ MatrixX _A;
+ bool linear=false;
+ MatrixX _W0;//(n x 1)
+ MatrixX _W;//(m x n)
+ std::string _name="Node";
+
+ Node():_m(0),_n(0){};
+
+ virtual MatrixX forward(const MatrixX& Z){
+ return Z;
+ }
+
+ virtual MatrixX backward(const MatrixX dLdA){
+ return dLdA;
+ }
+
+ virtual unsigned getWeightsNum(){
+ return this->_W.size()+this->_W0.size();
+ }
+ virtual int getM(){
+ return _m;
+ }
+ virtual int getN(){
+ return _n;
+ }
+ virtual void setW(const MatrixX& W){
+ _W=MatrixX(W);
+ }
+ virtual void setW0(const MatrixX& W0){
+ _W0=MatrixX(W0);
+ }
+
+ virtual bool isLinear(){
+ return this->linear;
+ }
+ virtual void displayWeights(){
+ std::cout<<_W0 <<'\n';
+ std::cout<<_W <<'\n';
+ }
+
+
+
+ // void sgd_step(const scalar lrate){
+ // }
+ virtual ~Node() {
+ // std::cout<<"deleted_"<<name<<" \n";
+ }
+ // private:
+ };
+ //..................................................................................................
+ class Linear: public Node {
+ public:
+ int _m;
+ int _n;
+ MatrixX _W0;//(n x 1)
+ MatrixX _W;//(m x n)
+ std::string _name="Linear";
+ MatrixX _A;
+ MatrixX _dLdW;
+ MatrixX _dLdW0;
+ bool linear=true;
+
+ Linear():_m(0),_n(0){};
+ Linear(int m,int n) {
+ _m=m;
+ _n=n;
+ // (in size, out size)
+ _W0.zeros(_n, 1);
+
+ // _W<< 1.24737338 << 0.28295388 << 0.69207227<<arma::endr
+ // << 1.58455078 << 1.32056292 <<-0.69103982<<arma::endr;
+ _W.randn(m,n);//change with tf.randomNormal(shape1,0,1.0 * m ** (-.5),'float32',0);
+ // _name="Linear("+m+","+n+")";
+
+
+ };
+
+ MatrixX forward(const MatrixX& A){
+ _A=MatrixX(A); // (m x b)
+ // std::cout << arma::size(_W.t()) << '\n';
+ // std::cout << arma::size(_A) << '\n';
+ // std::cout << _W.size()<< '\n';
+ MatrixX temp=_W.t()*_A;
+ return temp.each_col()+_W0 ; // (n x b)
+
+ }
+
+
+ MatrixX backward(const MatrixX dLdZ){ // dLdZ is (n x b)
+ _dLdW=_A*dLdZ.t();
+ _dLdW0= arma::resize(arma::sum(dLdZ,1),arma::size(_W0));
+ return _W *dLdZ;
+ }
+
+ unsigned getWeightsNum(){
+ return this->_W.size()+this->_W0.size();
+ }
+ int getM(){
+ return _m;
+ }
+ int getN(){
+ return _n;
+ }
+ void setW(const MatrixX& W){
+ _W=MatrixX(W);
+ }
+ void setW0(const MatrixX& W0){
+ _W0=MatrixX(W0);
+ }
+
+ bool isLinear(){
+ return this->linear;
+ }
+
+ void sgd_step(const scalar lrate){
+ _W= _W-(_dLdW* lrate);
+ _W0= _W0-(_dLdW0* lrate);
+ }
+
+ void displayWeights(){
+ std::cout<<_W0 <<'\n';
+ std::cout<<_W <<'\n';
+ }
+ // private:
+ };
+
+ //..................................................................................................
+
+ class Tanh: public Node {
+ public:
+ MatrixX _A;
+ std::string _name="Tanh";
+
+ Tanh() { };
+
+ MatrixX forward(const MatrixX& Z){
+ _A=arma::tanh(Z);
+ return _A;
+ }
+
+ MatrixX backward(const MatrixX dLdA){
+ return dLdA%(1.-_A%_A);
+ // (1. + _A) % (1. - _A) * 0.5;
+ }
+ unsigned getWeightsNum(){
+ return this->_W.size()+this->_W0.size();
+ }
+
+ // void sgd_step(const scalar lrate){
+ // }
+ // private:
+ };
+
+ class ReLU: public Node {
+ public:
+ MatrixX _A;
+ std::string _name="Relu";
+
+ ReLU() { };
+
+ MatrixX forward(const MatrixX& Z){
+ _A=MatrixX(Z);
+ for (unsigned j=0; j<Z.n_cols; j++)
+ for (unsigned i=0; i<Z.n_rows; i++)
+ _A(i, j) = Z(i, j) > 0.0 ? Z(i, j) : 0.0;
+ return _A;
+ }
+
+ MatrixX backward(const MatrixX dLdA){
+
+ MatrixX temp=MatrixX(dLdA);
+ for (unsigned j=0; j<temp.n_cols; j++)
+ for (unsigned i=0; i<temp.n_rows; i++)
+ temp(i, j) = _A(i, j) <= 0. ? 0.0 : dLdA(i, j);
+ return temp;
+
+ }
+ unsigned getWeightsNum(){
+ return this->_W.size()+this->_W0.size();
+ }
+
+ // void sgd_step(const scalar lrate){
+ // }
+ // private:
+ };
+
+ class Sigmoid: public Node {
+ public:
+ MatrixX _A;
+ std::string _name="Sigmoid";
+
+ Sigmoid() { };
+
+ MatrixX forward(const MatrixX& Z){
+ _A= 1.0 / (1.0 + arma::exp(-Z));
+ return _A;
+ }
+
+ MatrixX backward(const MatrixX& dLdA){
+ return dLdA % (1. - dLdA); //check wrong probably
+ }
+
+ // void sgd_step(const scalar lrate){
+ // }
+ // private:
+ };
+
+ class SoftMax: public Node {
+ public:
+ MatrixX _A;
+ std::string _name="SoftMax";
+
+ SoftMax() { };
+
+ MatrixX forward(const MatrixX& Z){
+ MatrixX temp=arma::sum(arma::exp(Z),0);
+ MatrixX temp1=arma::exp(Z);
+ _A=temp1.each_row()/temp;
+ //.each_col()
+ // tf.div(tf.exp(Z) , tf.sum(tf.exp(Z),0));
+ return _A;
+ }
+
+ MatrixX backward(const MatrixX& dLdZ){
+ return dLdZ;
+ }
+ MatrixXI predict (const MatrixX& Ypred) {
+ return arma::index_max(Ypred,0) ;
+ };
+
+ unsigned getWeightsNum(){
+ return this->_W.size()+this->_W0.size();
+ }
+
+ // void sgd_step(const scalar lrate){
+ // }
+ // private:
+ };
+
+ //..................................................................................................
+
+ class Loss: public Node {
+ public:
+ MatrixX _Ypred;
+ MatrixX _Y;
+ std::string _name;
+
+ Loss() { };
+
+ virtual MatrixX forward(const MatrixX& Ypred,const MatrixX& Y){
+
+ _Ypred =MatrixX( Ypred);
+ _Y = MatrixX(Y);
+
+ return -arma::sum( (_Y%arma::log(_Ypred)) ,0);//o?
+ //tf.sum(tf.mul(this.Y ,tf.log(this.Ypred))).neg();
+ }
+
+ virtual MatrixX backward(){
+ return _Ypred-_Y;
+ }
+ // virtual void sgd_step(const scalar lrate){
+ // }
+
+ // private:
+ };
+
+ class NLL: public Loss {
+ public:
+ MatrixX _Ypred;
+ MatrixX _Y;
+ std::string _name="NLL";
+
+ NLL() { };
+
+ MatrixX forward(const MatrixX& Ypred,const MatrixX& Y){
+
+ _Ypred =MatrixX( Ypred);
+ _Y = MatrixX(Y);
+
+ return -arma::sum( (_Y%arma::log(_Ypred)) ,0);//o?
+ //tf.sum(tf.mul(this.Y ,tf.log(this.Ypred))).neg();
+ }
+
+ MatrixX backward(){
+ return _Ypred-_Y;
+ }
+ // void sgd_step(const scalar lrate){
+ // }
+
+ // private:
+ };
+
+ class Sequential {
+ public:
+ std::vector<Node*> _modules;
+ Loss* _loss;
+ scalar _totalLoss=0;
+ // std::string _name;
+
+ Sequential(std::vector<Node*> modules,Loss* loss) {
+ _modules=modules;
+ _loss=loss;
+
+ };
+
+ MatrixX forward(const MatrixX& Xt){
+ MatrixX res=MatrixX(Xt);
+ for(unsigned i=0;i< _modules.size();i++ ){
+ res = _modules[i]->forward(res);
+ }
+ return res;
+ }
+
+ MatrixX backward(const MatrixX& delta){
+ MatrixX res=MatrixX(delta);
+ for(unsigned i=_modules.size()-1;i>=0.;i-- ){
+ res = _modules[i]->backward(res);
+ }
+ }
+ // void sgd_step(const scalar lrate){
+ // for(unsigned i=0;i< _modules.size();i++ ){
+ // _modules[i]->sgd_step(lrate);
+ // }
+ // }
+ MatrixXI step (const MatrixX& Xt,const MatrixX&Yt) { //add argument for weights
+ MatrixX Ypred = this->forward(Xt);
+ // std::cout << Ypred<< '\n';
+ // std::cout <<"Loss:"<<_loss->forward(Ypred, Yt)<< '\n';
+ // _totalLoss= arma::accu(_loss->forward(Ypred, Yt));
+ MatrixXI pred=this->predict ( Ypred);
+ MatrixXI truth=this->predict ( Yt);
+ _totalLoss=arma::accu(arma::abs(pred+(-truth)));
+ return pred;
+ };
+
+ scalar sumLoss(){
+ std::cout<<"Loss:"<<_totalLoss <<'\n';
+ return _totalLoss;
+
+ }
+
+ MatrixXI predict (const MatrixX& Ypred) {
+ return arma::index_max(Ypred,0) ;
+ };
+
+ unsigned getWeightsNum(){
+ MatrixX _W0;//(n x 1)
+ MatrixX _W;//(m x n)
+ unsigned res=0;
+ for(unsigned i=0;i< _modules.size();i++ ){
+
+ if(_modules[i]->isLinear()){
+ res += _modules[i]->getWeightsNum();
+ }
+
+ }
+ return res;
+
+ }
+
+ std::vector<scalar> slice(std::vector<scalar> const &v, int const m, int const n){
+ auto first = v.cbegin() + m;
+ auto last = v.cbegin()+ m + n ;
+
+ std::vector<scalar> vec(first, last);
+ return vec;
+ }
+
+ MatrixX vecToMatrix( std::vector<scalar> const &v,int const m, int const n){
+ MatrixX mat;
+ mat.zeros(m, n);
+ unsigned ii=0;
+ for(unsigned i=0;i<m;i++){
+ for(unsigned j=0;j<n;j++){
+ mat(i,j)=v[ii];
+ ii++;
+ }
+ }
+ return mat;
+
+ }
+
+ void displayVector(std::vector<scalar> const &v){
+
+ for(unsigned i=0;i<v.size();i++){
+ std::cout<<v[i] <<" ";
+ }
+ std::cout<<'\n';
+ }
+
+ void setWeights(std::vector<scalar> weightsVector){
+ int count=0;
+ for(unsigned i=0;i< _modules.size();i++ ){
+
+ if(_modules[i]->isLinear()){
+ int m=_modules[i]->getM();
+ int n=_modules[i]->getN();
+ int size=m*n;
+
+ //set W0
+ // displayVector(slice(weightsVector, count, n));
+ // std::cout<<vecToMatrix( slice(weightsVector, count, n),n ,1 ) <<'\n';
+ MatrixX newW0= vecToMatrix( slice(weightsVector, count, n),n ,1 );
+ _modules[i]->setW0(newW0);
+ count += n;
+
+ //set W
+ // displayVector(slice(weightsVector, count, size));
+ // std::cout<<vecToMatrix( slice(weightsVector, count, size),m ,n ) <<'\n';
+ MatrixX newW= vecToMatrix( slice(weightsVector, count, size),m ,n );
+ _modules[i]->setW(newW);
+ count += (m*n);
+ }
+ }
+ // std::cout<<count <<'\n';
+ }
+
+
+
+ // void sgd(){
+ void sgd(const MatrixX& X,const MatrixX& Y, const unsigned iters,const scalar lrate){
+ unsigned D= X.n_cols;
+ unsigned N= X.n_rows;
+ // var D=X.shape[0];
+ // var N=X.shape[1];
+ int sum_loss=0;
+ // var sum_loss=tf.tensor([0]);
+ for(unsigned i=0;i<iters;i++){
+ unsigned j=(unsigned)(rand() % N);
+ }
+
+ // for(var i=0; i < iters; i++){
+ // var j= getRandomInt(N);
+ // var Xt = X.slice([0,j], [D,1]);
+ // var Yt = Y.slice([0,j], [D,1]);
+ // var Ypred = this.forward(Xt);
+ // sum_loss.add(sum_loss, this.loss.forward(Ypred, Yt));
+ // var err = this.loss.backward();
+ // this.backward(err);
+ // this.sgd_step(lrate);
+
+ // }
+ // D, N = X.shape
+ // for it in range(iters):
+ // j = np.random.randint(N)
+ // Xt = X[:,j:j+1]; Yt = Y[:,j:j+1]
+ // print(Xt)
+ // Ypred=self.forward( Xt)
+ // l=self.loss.forward(Ypred, Yt)
+ // delta=self.loss.backward()
+ // b=self.backward( delta)
+ // s=self.sgd_step( lrate)
+ // }
+ }
+
+ // private:
+ };
+
+ //..................................................................................................
+
+
+
+
+}
diff --git a/cpp_dnn/frep.cpp b/cpp_dnn/frep.cpp
new file mode 100644
index 0000000..71393f9
--- /dev/null
+++ b/cpp_dnn/frep.cpp
@@ -0,0 +1,707 @@
+#include <iostream>
+#include <math.h>
+#include <vector>
+#include <sstream> // std::ostringstream
+
+
+
+namespace cpp_frep {
+ using scalar=float;
+
+ class Frep {
+ public:
+ scalar X;
+ scalar Y;
+ std::string name;
+ std::vector<scalar> params;
+ Frep* shape1;
+ Frep* shape2;
+ Frep* shape;
+
+ Frep():X(0.0f),Y(0.0f){};
+ Frep(std::string name_,scalar X_, scalar Y_) {
+ X=X_;
+ Y=Y_;
+ name=name_;
+ };
+
+ Frep get(){
+ return *this;
+ }
+
+ void setXY(scalar X_,scalar Y_){
+ X=X_;
+ Y=Y_;
+ }
+
+ virtual scalar eval(){
+ return 0.0f;
+ }
+ virtual void setParameters(std::vector<scalar> params_){
+ params=params_;
+ }
+ virtual void printName(){
+ std::cout<<"Frep";
+ }
+ virtual std::string getPrintedName(){
+ return "Frep";
+ }
+
+ void setShape1(Frep* shape1_){
+ shape1=shape1_;
+ }
+ void setShape2(Frep* shape2_){
+ shape2=shape2_;
+ }
+ void setShape(Frep* shape_){
+ shape=shape_;
+ }
+ virtual ~Frep() {
+ // std::cout<<"deleted_"<<name<<" \n";
+ }
+
+ // private:
+
+ };
+ //..........................SHAPES..................................//
+ class FrepPrimitive: public Frep {
+ public:
+ scalar center_x;
+ scalar center_y;
+
+ FrepPrimitive():center_x(0.0f),center_y(0.0f){};
+ FrepPrimitive(scalar center_x_, scalar center_y_) {
+ center_x=center_x_;
+ center_y=center_y_;
+ };
+
+ FrepPrimitive get(){
+ return *this;
+ }
+
+ virtual scalar eval(){
+ return 0.0f;
+ }
+
+ virtual void printName(){
+ std::cout<<name;
+ }
+
+ virtual std::string getPrintedName(){
+ return name;
+ }
+
+
+ };
+
+ class Circle: public FrepPrimitive {
+ public:
+ scalar radius;
+
+ Circle():radius(1.0f){
+ name="Circle";
+ };
+
+ Circle(scalar radius_,scalar center_x_, scalar center_y_) {
+ radius=radius_;
+ center_x=center_x_;
+ center_y=center_y_;
+ name="Circle";
+ };
+
+ scalar eval(){
+ return ((radius)*(radius)-((X-(center_x))*(X-(center_x))+(Y-(center_y))*(Y-(center_y))));
+ }
+
+ Circle get(){
+ return *this;
+ }
+
+
+ // private:
+
+ };
+
+ class Rectangle: public FrepPrimitive {
+ public:
+ scalar xmin;
+ scalar xmax;
+ scalar ymin;
+ scalar ymax;
+
+ Rectangle():xmin (-1.0f), xmax (1.0f), ymin (-1.0f), ymax (1.0f) {
+ name="Rectangle";
+ };
+ Rectangle(scalar xmin_, scalar xmax_,scalar ymin_,scalar ymax_) {
+ xmin=xmin_;
+ xmax=xmax_;
+ ymin=ymin_;
+ ymax=ymax_;
+ name="Rectangle";
+ };
+
+ scalar eval(){
+ scalar xfn = std::min(X-(xmin),(xmax)-X);
+ scalar yfn = std::min(Y-(ymin),(ymax)-Y);
+ scalar fn = std::min((xfn),(yfn));
+ return fn;
+ }
+
+ Rectangle get(){
+ return *this;
+ }
+ // private:
+
+ };
+
+ class InvoluteGear: public FrepPrimitive {
+ public:
+ scalar module =1.0f;
+ scalar angle =20.0f;
+ scalar teeth =10.0f;
+ scalar addendum =1.0f;
+ scalar dedendum =1.1f;
+ const scalar PI_F=3.14159265358979;
+
+
+ InvoluteGear():module (1.0f), addendum (1.0f), dedendum (1.1f), angle (20.0f), teeth (10.0f) {
+ name="Rectangle";
+ };
+ InvoluteGear(scalar module_, scalar addendum_,scalar dedendum_,scalar angle_,scalar teeth_) {
+ module =module_;
+ angle =angle_;
+ teeth =teeth_;
+ addendum =addendum_;
+ dedendum =dedendum_;
+ name="Gear";
+ };
+
+ scalar eval(){
+ scalar rp = module*teeth/2.0f; // pitch radius
+ scalar rb = rp*cos(angle); // base radius
+ scalar ha = addendum*module; // addendum height
+ scalar hd= dedendum*module; // dedendum height
+ scalar ai = tan(angle)-angle; // involute angle
+ scalar ap = 2.0f*PI_F/teeth; // pitch angle
+
+ angle=PI_F*angle/180.0f;
+
+ scalar fn=std::max(std::min(std::min((rp+ha)-sqrtf(X*X+Y*Y),remainderf((PI_F+atan2f(Y,X)),ap)-(sqrtf(pow(std::max(rb,sqrtf(X*X+Y*Y))/rb,2.0)-1.0)-acosf(rb/std::max(rb,sqrt(X*X+Y*Y))))),-(sqrtf(powf(std::max(rb,sqrtf(X*X+Y*Y))/rb,2.0)-1.0)-acosf(rb/std::max(rb,sqrtf(X*X+Y*Y))))-(-(ap/2+2*ai)+remainderf((PI_F+atan2f(Y,X)),ap))),(rp-hd)-sqrtf(X*X+Y*Y));
+
+ return fn;
+ }
+
+ InvoluteGear get(){
+ return *this;
+ }
+ // private:
+
+ };
+ //..................................................................//
+
+ //.......................OPERATORS..................................//
+ class FrepOperators: public Frep {
+ public:
+
+ FrepOperators(){};
+ virtual void printName(){
+ std::cout<<name<<"(";
+ shape1->printName();
+ std::cout<<",";
+ shape2->printName();
+ std::cout<<")";
+ }
+
+ virtual std::string getPrintedName(){
+ std::string frepName=name+"("+shape1->getPrintedName()+","+shape2->getPrintedName()+")" ;
+ return frepName;
+ }
+
+ };
+
+ class Add: public FrepOperators {
+ public:
+
+ Add(){
+ name="Add";
+ };
+ Add(Frep* shape1_,Frep* shape2_) {
+ shape1=shape1_;
+ shape2=shape2_;
+ name="Add";
+ };
+
+ scalar eval(){
+ //change x and y
+ shape1->setXY(X,Y);
+ shape2->setXY(X,Y);
+ // std::cout<<"addddd\n";
+ return std::max(shape1->eval(),shape2->eval());
+ }
+
+ Add get(){
+ return *this;
+ }
+
+
+ };
+
+ class Subtract: public FrepOperators {
+ public:
+ Subtract(){
+ name="Subtract";
+ };
+ Subtract(Frep* shape1_,Frep* shape2_) {
+ shape1=shape1_;
+ shape2=shape2_;
+ name="Subtract";
+ };
+
+ scalar eval(){
+ //change x and y
+ shape1->setXY(X,Y);
+ shape2->setXY(X,Y);
+ return std::min(shape1->eval(),-shape2->eval());
+ }
+
+
+ Subtract get(){
+ return *this;
+ }
+
+
+ };
+
+ class Intersect: public FrepOperators {
+ public:
+
+ Intersect(){
+ name="Intersect";
+ };
+ Intersect(Frep* shape1_,Frep* shape2_) {
+ shape1=shape1_;
+ shape2=shape2_;
+ name="Intersect";
+ };
+
+ scalar eval(){
+ //change x and y
+ shape1->setXY(X,Y);
+ shape2->setXY(X,Y);
+ return std::min(shape1->eval(),shape2->eval());
+ }
+
+
+ Intersect get(){
+ return *this;
+ }
+
+ };
+
+ class Clearance: public FrepOperators {
+ public:
+
+ scalar offset;
+
+ Clearance():offset(0.5f){
+ name="Clearance";
+ };
+ Clearance(Frep* shape1_,Frep* shape2_,scalar offset_) {
+ shape1=shape1_;
+ shape2=shape2_;
+ offset=offset_;
+ name="Clearance";
+ };
+
+ scalar eval(){
+ //change x and y
+ shape1->setXY(X,Y);
+ shape2->setXY(X,Y);
+ return std::min(shape1->eval(),-(offset+shape2->eval()));
+ }
+ void setParameters(std::vector<scalar> offset_){
+ offset=offset_[0];
+ }
+
+ Clearance get(){
+ return *this;
+ }
+ };
+
+ class Blend: public FrepOperators {
+ public:
+
+ scalar offset;
+
+ Blend():offset(0.5f){
+ name="Blend";
+ };
+ Blend(Frep* shape1_,Frep* shape2_,scalar offset_) {
+ shape1=shape1_;
+ shape2=shape2_;
+ offset=offset_;
+ name="Blend";
+ };
+
+ Blend get(){
+ return *this;
+ }
+
+ void setParameters(std::vector<scalar> offset_){
+ offset=offset_[0];
+ }
+
+
+ scalar eval(){
+ //change x and y
+ shape1->setXY(X,Y);
+ shape2->setXY(X,Y);
+ scalar fn =std::max(shape1->eval(),shape2->eval());
+ return std::max(fn,(shape1->eval())+(shape2->eval())+(offset));
+ }
+
+ };
+
+ class Morph: public FrepOperators {
+ public:
+
+ scalar weight;
+
+ Morph():weight(0.5f){
+ name="Morph";
+ };
+ Morph(Frep* shape1_,Frep* shape2_,scalar weight_) {
+ shape1=shape1_;
+ shape2=shape2_;
+ weight=weight_;
+ name="Morph";
+ };
+
+ void setParameters(std::vector<scalar> weight_){
+ weight=weight_[0];
+ }
+
+ Morph get(){
+ return *this;
+ }
+
+ scalar eval(){
+ //change x and y
+ // std::cout<<"xv cv\n";
+ shape1->setXY(X,Y);
+ shape2->setXY(X,Y);
+ // std::cout<<"moorrrphhhh\n";
+ return (weight*shape1->eval()+(1-weight)*shape2->eval());
+ }
+
+ std::string to_string_with_precision(const scalar a_value, const int n = 2)
+ {
+ std::ostringstream out;
+ out.precision(n);
+ out << std::fixed << a_value;
+ return out.str();
+ }
+
+ std::string getPrintedName(){
+ std::string frepName=name+"("+shape1->getPrintedName()+","+shape2->getPrintedName()+","+ to_string_with_precision(weight,2) +")" ;
+ return frepName;
+ }
+
+ };
+ //..................................................................//
+
+ //.......................TRANSFORMS.................................//
+ class FrepTransforms: public Frep {
+ public:
+
+ //Frep(std::string const& name): name_(std::move(name)) {}
+ FrepTransforms(){};
+ virtual void printName(){
+ std::cout<<name<<"(";
+ shape->printName();
+ std::cout<<")";
+ }
+ virtual std::string getPrintedName(){
+ std::string frepName=name+"("+shape->getPrintedName()+")" ;
+ return frepName;
+ }
+
+
+
+ };
+
+ class Translate: public FrepTransforms {
+ public:
+ scalar dx;
+ scalar dy;
+
+ Translate():dx(1.0f),dy(1.0f){
+ name="Translate";
+ };
+ Translate(Frep* shape_,scalar dx_,scalar dy_) {
+ shape= shape_;
+ dx=dx_;
+ dy=dy_;
+ name="Translate";
+ };
+
+ scalar eval(){
+ //change x and y
+ shape->setXY(X-dx,Y-dy);
+ return shape->eval();
+ }
+
+ void setParameters(std::vector<scalar> params_){
+ dx=params_[0];
+ dy=params_[1];
+ }
+
+ Translate get(){
+ return *this;
+ }
+ };
+
+ class Scale: public FrepTransforms {
+ public:
+ scalar sx;
+ scalar sy;
+ scalar ox;
+ scalar oy;
+
+ Scale():sx(0.5f),sy(0.5f),ox(0.0f),oy(0.0f){
+ name="Scale";
+ };
+ Scale(Frep* shape_,scalar ox_, scalar oy_,scalar sx_,scalar sy_) {
+ shape=shape_;
+ sx=sx_;
+ sy=sy_;
+ ox=ox_;
+ oy=oy_;
+ name="Scale";
+ };
+
+ scalar eval(){
+ //change x and y
+ shape->setXY(((ox)+(X-(ox))/(sx)),((oy)+(Y-(oy))/(sy)));
+ return shape->eval();
+ }
+
+ void setParameters(std::vector<scalar> params_){
+ sx=params_[0];
+ sy=params_[1];
+ ox=params_[2];
+ oy=params_[3];
+ }
+
+ Scale get(){
+ return *this;
+ }
+ };
+
+ class Rotate: public FrepTransforms {
+ public:
+ scalar angle;
+ scalar ox;
+ scalar oy;
+ const scalar PI_F=3.14159265358979;
+
+ Rotate():angle(45.0f),ox(0.0f),oy(0.0f){
+ name="Rotate";
+ };
+ Rotate(Frep* shape_,scalar ox_, scalar oy_, scalar angle_) {
+ shape=shape_;
+ angle=angle_;
+ ox=ox_;
+ oy=oy_;
+ name="Rotate";
+ };
+
+ Rotate get(){
+ return *this;
+ }
+
+ void setParameters(std::vector<scalar> params_){
+ angle=params_[0];
+ ox=params_[1];
+ oy=params_[2];
+ }
+
+ scalar eval(){
+ //change x and y
+
+ scalar r =angle*PI_F/180.0f;
+ scalar newX=((ox)+(X-(ox))*cos(r)+(Y-(oy))*sin(r));
+ scalar newY=((oy)-(X-(ox))*sin(r)+(Y-(oy))*cos(r));
+ shape->setXY(newX,newY);
+ // std::cout<<"rooottaattte\n";
+
+ // std::cout<< shape->name<<"\n";
+ scalar result=shape->eval();
+ // std::cout<<"dvds\n";
+ return result;
+ }
+
+
+ };
+
+ class LinearArrayX: public FrepTransforms {
+ public:
+ scalar spacing;
+ scalar number;
+ scalar size;
+
+ const scalar PI_F=3.14159265358979;
+
+ LinearArrayX():spacing(0.5),number(2.0),size(1.0){
+ name="LinearArrayX";
+ };
+ LinearArrayX(Frep* shape_,scalar spacing_,scalar number_,scalar size_) {
+ shape=shape_;
+ spacing=spacing_;
+ number=number_;
+ size=size_;
+ name="LinearArrayX";
+ };
+
+ LinearArrayX get(){
+ return *this;
+ }
+
+ void setParameters(std::vector<scalar> params_){
+ spacing=params_[0];
+ number=params_[1];
+ size=params_[2];
+ }
+
+
+ scalar eval(){
+ //change x and y
+ scalar xmin=-size;
+ scalar xmax=size;
+
+ scalar newX=remainderf( X ,spacing+xmax-xmin );
+ shape->setXY(newX,Y);
+ scalar fn=shape->eval();
+
+ fn = std::min(X-(xmin),fn);
+ fn = std::min((xmin+number*(spacing+xmax-xmin))-X,fn);
+
+ return fn;
+ }
+
+
+ };
+
+ class LinearArrayY: public FrepTransforms {
+ public:
+ scalar spacing;
+ scalar number;
+ scalar size;
+
+ const scalar PI_F=3.14159265358979;
+
+ LinearArrayY():spacing(0.5),number(2.0),size(1.0){
+ name="LinearArrayY";
+ };
+ LinearArrayY(Frep* shape_,scalar spacing_,scalar number_,scalar size_) {
+ shape=shape_;
+ spacing=spacing_;
+ number=number_;
+ size=size_;
+ name="LinearArrayY";
+ };
+
+ LinearArrayY get(){
+ return *this;
+ }
+
+ void setParameters(std::vector<scalar> params_){
+ spacing=params_[0];
+ number=params_[1];
+ size=params_[2];
+ }
+
+ scalar eval(){
+ //change x and y
+ scalar ymin=-size;
+ scalar ymax=size;
+
+ scalar newY=remainderf( Y ,spacing+ymax-ymin );
+ shape->setXY(X,newY);
+ scalar fn=shape->eval();
+
+ fn = std::min(Y-(ymin),fn);
+ fn = std::min((ymin+number*(spacing+ymax-ymin))-Y,fn);
+
+ return fn;
+ }
+
+
+ };
+
+ class PolarArray: public FrepTransforms {
+ public:
+ scalar angle;
+ scalar radius;
+ scalar ox;
+ scalar oy;
+ const scalar PI_F=3.14159265358979;
+
+ PolarArray():angle(45.0),radius(1.0),ox(0.0),oy(0.0){
+ name="PolarArray";
+ };
+ PolarArray(Frep* shape_,scalar angle_,scalar radius_) {
+ shape=shape_;
+ angle=angle_;
+ radius=radius_;
+ name="PolarArray";
+
+ };
+
+ PolarArray get(){
+ return *this;
+ }
+
+ void setParameters(std::vector<scalar> params_){
+ angle=params_[0];
+ radius=params_[1];
+ ox=params_[2];
+ oy=params_[3];
+ }
+
+ scalar eval(){
+ //change x and y
+
+ scalar r =angle*PI_F/180.0f;
+ ox=-radius;
+ oy=radius;
+ scalar tmin = atan2f(-oy,-ox);
+
+ scalar newX;
+ scalar newY;
+
+ newX=((0.0f)+(X-(0.0f))*cos(r/2.0f)+(Y-(0.0f))*sin(r/2.0f));
+ newY=((0.0f)-(X-(0.0f))*sin(r/2.0f)+(Y-(0.0f))*cos(r/2.0f));
+
+ X=newX;
+ Y=newY;
+
+ X-=radius;
+ Y+=radius;
+
+ newX=(ox)+sqrt((X-(ox))*(X-(ox))+(Y-(oy))*(Y-(oy)))*cos((tmin)+remainderf((2.0f*PI_F+atan2f(Y-(oy),X-(ox))-(tmin)),(r)));
+ newY=(oy)+sqrt((X-(ox))*(X-(ox))+(Y-(oy))*(Y-(oy)))*sin((tmin)+remainderf((2.0f*PI_F+atan2f(Y-(oy),X-(ox))-(tmin)),(r)));
+
+
+ shape->setXY(newX,newY);
+ return shape->eval();
+ }
+
+
+ };
+ //..................................................................//
+
+
+
+
+}
diff --git a/cpp_dnn/image.cpp b/cpp_dnn/image.cpp
new file mode 100644
index 0000000..0bb6dbe
--- /dev/null
+++ b/cpp_dnn/image.cpp
@@ -0,0 +1,143 @@
+#include "nelder_mead.h"
+#include "frep.cpp"
+#include <iostream>
+#include <math.h>
+#include <vector>
+#include <stdio.h>
+#include <stdlib.h>
+#include <opencv2/opencv.hpp>
+#include <armadillo>
+#include <iostream>
+
+
+
+namespace cpp_dnn {
+ using scalar=double;
+ using MatrixX = arma::Mat<scalar>;
+ using MatrixXI = arma::Mat<long long unsigned int>;
+
+ class Image {
+ public:
+ scalar res=0.01;
+ scalar bound =3.0;
+ scalar minX;
+ scalar maxX;
+ scalar minY;
+ scalar maxY;
+ scalar dx;
+ scalar dy;
+ MatrixX loc;
+ MatrixX target;
+ int height;
+ int width;
+
+ Image():res(0.01),bound(3.0){
+ minX=-bound;
+ maxX=bound;
+ minY=-bound;
+ maxY=bound;
+ dx=res;
+ dy=res;
+ setLocMatrix();
+ };
+
+ Image(scalar res, scalar bound){
+ minX=-bound;
+ maxX=bound;
+ minY=-bound;
+ maxY=bound;
+ dx=res;
+ dy=res;
+ setLocMatrix();
+
+ };
+
+ Image(scalar res, scalar boundX,scalar boundY){
+ minX=-boundX;
+ maxX=boundX;
+ minY=-boundY;
+ maxY=boundY;
+ dx=res;
+ dy=res;
+ setLocMatrix();
+ };
+
+ void setLocMatrix(){
+ height=(int)((maxX-minX)/dx);
+ width=(int)((maxY-minY)/dy);
+ this->loc.zeros(height*width,2);
+ unsigned ii=0;
+ for(scalar i=minX;i<maxX-dx;i+=dx){
+ for(scalar j=minY;j<maxY-dy;j+=dy){
+ this->loc(ii,0)=i;
+ this->loc(ii,1)=j;
+ ii++;
+ }
+ }
+ }
+
+ MatrixX getLocMatrix(){
+ return this->loc.t();
+ }
+
+ MatrixXI predict (const MatrixX& Ypred) {
+ return arma::index_max(Ypred,0) ;
+ }
+
+ MatrixX getTarget( cpp_frep::Frep* frep,std::string name){
+ target.zeros(height*width,2);
+ for (unsigned i=0; i<this->loc.n_rows; i++){
+ frep->setXY(this->loc(i,0),this->loc(i,1));
+ if(frep->eval()>0.0){
+ target(i, 0) = 0.0;
+ target(i, 1) = 1.0;
+ }else{
+ target(i, 0) = 1.0;
+ target(i, 1) = 0.0;
+ }
+ }
+ saveImage(predict(target.t()) ,"Target_"+ name);
+ return target.t();
+ }
+
+ void saveImage(const MatrixXI& Ypred,std::string name){
+ const int size=height*width*4;
+ int stride=0;
+ unsigned char pix[size];
+
+ for(unsigned i=0;i<Ypred.size();i++){
+ if(Ypred(i)>0.0){
+ pix[stride+0]=0;
+ pix[stride+1]=0;
+ pix[stride+2]=0;
+ pix[stride+3]=0;
+
+ }else{
+ pix[stride+0]=255;
+ pix[stride+1]=255;
+ pix[stride+2]=255;
+ pix[stride+3]=255;
+ }
+ stride+=4;
+ }
+
+ ////////////////create image and save it//////////////////////////
+ cv::Mat image = cv::Mat(width, height, CV_8UC4, pix);
+ std::string fileName="img/"+name+".jpg";
+ cv::imwrite( fileName, image );
+ bool preview=false;
+ if(preview){
+ cv::namedWindow("Display Image", cv::WINDOW_AUTOSIZE );
+ cv::imshow("Display Image", image);
+ cv::waitKey(0);
+ }
+
+ }
+
+ virtual ~Image() {
+ // std::cout<<"deleted_"<<name<<" \n";
+ }
+
+ };
+}
+
diff --git a/cpp_dnn/main.cpp b/cpp_dnn/main.cpp
new file mode 100644
index 0000000..9696e8a
--- /dev/null
+++ b/cpp_dnn/main.cpp
@@ -0,0 +1,185 @@
+#include "dnn.cpp"
+#include "Image.cpp"
+#include "vector.h"
+#include <iostream>
+#include <armadillo>
+
+#include <stdio.h>
+#include <stdlib.h> /* free() */
+#include <stddef.h> /* NULL */
+
+#include "cmaes_interface.h"
+
+using namespace cpp_dnn;
+using scalar=double;
+
+scalar getRandom(){//between 0 and 1
+ return ((scalar) rand() / (RAND_MAX));
+}
+
+std::vector<scalar> getRandomWeights(unsigned numWeights){
+ std::vector<scalar> v;
+
+ for(unsigned i=0;i<numWeights;i++){
+ v.push_back(getRandom()*3.0-6.0);
+ }
+ return v;
+}
+
+/* the objective (fitness) function to be minimized */
+// double fitfun(double const *x, int N) { /* function "cigtab" */
+// int i;
+// double sum = 1e4*x[0]*x[0] + 1e-4*x[1]*x[1];
+// for(i = 2; i < N; ++i)
+// sum += x[i]*x[i];
+// return sum;
+// }
+std::vector<scalar> pointerToVector(double const *x, int N){
+ std::vector<scalar> v;
+ for(int i = 0; i < N; ++i)
+ v.push_back(x[i]);
+ return v;
+}
+
+// /* the objective (fitness) function to be minimized */
+double fitfun(double const *x, int N,MatrixX X,MatrixX Y,Sequential seq) { /* function "cigtab" */
+ seq.setWeights( pointerToVector(x, N));
+ seq.step(X,Y);
+ return seq.sumLoss();
+}
+
+//--------------------------------------------------------------------------------------------------
+int main(int const, char const**) {
+
+ //Create network
+ std::vector<Node*> network;
+
+ Linear l1(2, 12);
+ Tanh a1;
+
+ Linear l2(12, 12);
+ Tanh a2;
+
+ Linear l3(12, 6);
+ Tanh a3;
+
+ Linear l4(6, 6);
+ Tanh a4;
+
+ Linear l5(6, 2);
+ SoftMax a5;
+
+ NLL l;
+ network.push_back( &l1);
+ network.push_back( &a1);
+ network.push_back( &l2);
+ network.push_back( &a2);
+ network.push_back( &l3);
+ network.push_back( &a3);
+ network.push_back( &l4);
+ network.push_back( &a4);
+ network.push_back( &l5);
+ network.push_back( &a5);
+
+ Sequential seq(network, &l);
+
+ //Create Image
+ Image Im;
+ MatrixX X=Im.getLocMatrix();
+ cpp_frep::Circle frep(1.5f,0.0f,0.0f);
+ MatrixX Y=Im.getTarget( &frep ,"Circle");
+ // Im.saveImage(seq.step(X,Y),"trial"); //test
+
+ ///////////////////////////// //Optimize neldermead
+ unsigned numWeights=seq.getWeightsNum();
+
+ // std::vector<scalar> v=getRandomWeights(numWeights);
+ // std::cout << seq.sumLoss()<< '\n';
+ // seq.setWeights(v);
+ // Im.saveImage(seq.step(X,Y),"trial1"); //test
+ // std::cout << seq.sumLoss()<< '\n';
+
+
+ // float precision = 0.001;
+ // int dimension = (int)numWeights;
+ // unsigned maxSteps=1000;
+ // NelderMeadOptimizer o(dimension, precision);
+ // std::cout << dimension<< '\n';
+
+ // request a simplex to start with
+ // Vector v(getRandomWeights(numWeights));
+ // Vector v1(getRandomWeights(numWeights));
+ // Vector v2(getRandomWeights(numWeights));
+ // o.insert(v);
+ // o.insert(v1);
+ // o.insert(v2);
+ // unsigned count=0;
+
+ // while (!o.done() && count<maxSteps) {
+ // seq.setWeights(v.get());
+ // Im.saveImage(seq.step(X,Y),"trial"+std::to_string(count));
+ // v = o.step(v, -seq.sumLoss()); //later change f(v)
+ // count++;
+
+ // }
+
+
+ ///////////////////////////// //Optimize CMAES
+ cmaes_t evo; /* an CMA-ES type struct or "object" */
+ double *arFunvals, *const*pop, *xfinal;
+ int i;
+ arFunvals = cmaes_init(&evo, 5, NULL, NULL, 0, 0, "cpp_dnn/cmaes_initials.par");
+ printf("%s\n", cmaes_SayHello(&evo));
+
+ cmaes_ReadSignals(&evo, "cmaes_signals.par"); /* write header and initial values */
+
+ /* Iterate until stop criterion holds */
+ while(!cmaes_TestForTermination(&evo))
+ {
+ /* generate lambda new search points, sample population */
+ pop = cmaes_SamplePopulation(&evo); /* do not change content of pop */
+
+ /* Here we may resample each solution point pop[i] until it
+ becomes feasible. function is_feasible(...) needs to be
+ user-defined.
+ Assumptions: the feasible domain is convex, the optimum is
+ not on (or very close to) the domain boundary, initialX is
+ feasible and initialStandardDeviations are sufficiently small
+ to prevent quasi-infinite looping. */
+ /* for (i = 0; i < cmaes_Get(&evo, "popsize"); ++i)
+ while (!is_feasible(pop[i]))
+ cmaes_ReSampleSingle(&evo, i);
+ */
+
+ /* evaluate the new search points using fitfun */
+ for (i = 0; i < cmaes_Get(&evo, "lambda"); ++i) {
+ arFunvals[i] = fitfun(pop[i], (int) cmaes_Get(&evo, "dim"),X,Y,seq);
+ }
+
+ /* update the search distribution used for cmaes_SamplePopulation() */
+ cmaes_UpdateDistribution(&evo, arFunvals);
+
+ /* read instructions for printing output or changing termination conditions */
+ cmaes_ReadSignals(&evo, "cmaes_signals.par");
+ fflush(stdout); /* useful in MinGW */
+ }
+ printf("Stop:\n%s\n", cmaes_TestForTermination(&evo)); /* print termination reason */
+ cmaes_WriteToFile(&evo, "all", "allcmaes.dat"); /* write final results */
+
+ /* get best estimator for the optimum, xmean */
+ xfinal = cmaes_GetNew(&evo, "xmean"); /* "xbestever" might be used as well */
+ cmaes_exit(&evo); /* release memory */
+
+ /* do something with final solution and finally release memory */
+ std::cout<<*xfinal<<'\n';
+
+ free(xfinal);
+
+ seq.setWeights( pointerToVector(pop[0], (int) cmaes_Get(&evo, "dim")));
+ Im.saveImage(seq.step(X,Y),"trialll");
+
+
+
+
+ return 0;
+}
diff --git a/cpp_dnn/nelder_mead.h b/cpp_dnn/nelder_mead.h
new file mode 100644
index 0000000..92ec8f5
--- /dev/null
+++ b/cpp_dnn/nelder_mead.h
@@ -0,0 +1,276 @@
+//based on https://github.com/blinry/nelder-mead-optimizer
+#include <ctime>
+#include <vector>
+#include <map>
+#include <cmath>
+#include <algorithm>
+using namespace std;
+
+// Float vector with standard operations
+class Vector {
+ public:
+ Vector() {
+ }
+ Vector(float c0) {
+ coords.push_back(c0);
+ }
+ Vector(float c0, float c1) {
+ coords.push_back(c0);
+ coords.push_back(c1);
+ }
+ Vector(float c0, float c1, float c2) {
+ coords.push_back(c0);
+ coords.push_back(c1);
+ coords.push_back(c2);
+ }
+ Vector(vector<float> coords_) {
+ coords=coords_;
+ }
+ Vector(vector<double> coords_) {
+ coords= vector<float>(coords_.begin(), coords_.end());
+ }
+ vector<double> get(){
+ return vector<double>(coords.begin(), coords.end());
+ }
+
+ // add more constructors when N gets > 3
+
+ float& operator[](int i) {
+ return coords[i];
+ }
+ float at(int i) const {
+ return coords[i];
+ }
+ int dimension() const {
+ return coords.size();
+ }
+ void prepare(int size) {
+ for (int i=0; i<size; i++) {
+ coords.push_back(0);
+ }
+ }
+ Vector operator+(Vector other) {
+ Vector result;
+ result.prepare(dimension());
+ for (int i=0; i<dimension(); i++) {
+ result[i] = coords[i] + other[i];
+ }
+ return result;
+ }
+ void operator+=(Vector other) {
+ for (int i=0; i<dimension(); i++) {
+ coords[i] += other[i];
+ }
+ }
+ Vector operator-(Vector other) {
+ Vector result;
+ result.prepare(dimension());
+ for (int i=0; i<dimension(); i++) {
+ result[i] = coords[i] - other[i];
+ }
+ return result;
+ }
+ bool operator==(Vector other) {
+ if (dimension() != other.dimension()) {
+ return false;
+ }
+ for (int i=0; i<dimension(); i++) {
+ if (other[i] != coords[i]) {
+ return false;
+ }
+ }
+ return true;
+ }
+ Vector operator*(float factor) {
+ Vector result;
+ result.prepare(dimension());
+ for (int i=0; i<dimension(); i++) {
+ result[i] = coords[i]*factor;
+ }
+ return result;
+ }
+ Vector operator/(float factor) {
+ Vector result;
+ result.prepare(dimension());
+ for (int i=0; i<dimension(); i++) {
+ result[i] = coords[i]/factor;
+ }
+ return result;
+ }
+ void operator/=(float factor) {
+ for (int i=0; i<dimension(); i++) {
+ coords[i] /= factor;
+ }
+ }
+ bool operator<(const Vector other) const {
+ for (int i=0; i<dimension(); i++) {
+ if (at(i) < other.at(i))
+ return false;
+ else if (at(i) > other.at(i))
+ return true;
+ }
+ return false;
+ }
+ float length() {
+ float sum = 0;
+ for (int i=0; i<dimension(); i++) {
+ sum += coords[i]*coords[i];
+ }
+ return pow(sum, 0.5f);
+ }
+ private:
+ vector<float> coords;
+};
+
+// This class stores known values for vectors. It throws unknown vectors.
+class ValueDB {
+ public:
+ ValueDB() {
+ }
+ float lookup(Vector vec) {
+ if (!contains(vec)) {
+ throw vec;
+ } else {
+ return values[vec];
+ }
+ }
+ void insert(Vector vec, float value) {
+ values[vec] = value;
+ }
+ void reset(){
+ values.clear();
+ }
+ private:
+ bool contains(Vector vec) {
+ map<Vector, float>::iterator it = values.find(vec); // TODO add tolerance
+ return it != values.end();
+ }
+ map<Vector, float> values;
+};
+
+class NelderMeadOptimizer {
+ public:
+ NelderMeadOptimizer(int dimension, float termination_distance=0.001) {
+ this->dimension = dimension;
+ srand(time(NULL));
+ alpha = 1;
+ gamma = 2;
+ rho = -0.5;
+ sigma = 0.5;
+ this->termination_distance = termination_distance;
+ }
+ // used in `step` to sort the vectors
+ bool operator()(const Vector& a, const Vector& b) {
+ return db.lookup(a) < db.lookup(b);
+ }
+ // termination criteria: each pair of vectors in the simplex has to
+ // have a distance of at most `termination_distance`
+ bool done() {
+ if ((int)vectors.size() < dimension) {
+ return false;
+ }
+ for (int i=0; i<dimension+1; i++) {
+ for (int j=0; j<dimension+1; j++) {
+ if (i==j) continue;
+ if ((vectors[i]-vectors[j]).length() > termination_distance) {
+ return false;
+ }
+ }
+ }
+ return true;
+ }
+
+ bool done(float score) {
+ if (score<1.0f+termination_distance && score>1.0f-termination_distance) {
+ return true;
+ }
+ return false;
+
+ }
+ void reset(){
+ vectors.clear();
+ db.reset();
+ }
+
+ void insert(Vector vec) {
+ if ((int)vectors.size() < dimension+1) {
+ vectors.push_back(vec);
+ }
+ }
+
+ Vector step(Vector vec, float score) {
+ db.insert(vec, score);
+ try {
+ if ((int)vectors.size() < dimension+1) {
+ vectors.push_back(vec);
+ }
+
+ // otherwise: optimize!
+ if ((int)vectors.size() == dimension+1) {
+ while(!done(score)) {
+ sort(vectors.begin(), vectors.end(), *this);
+ Vector cog; // center of gravity
+ cog.prepare(dimension);
+ for (int i = 1; i<=dimension; i++) {
+ cog += vectors[i];
+ }
+ cog /= dimension;
+ Vector best = vectors[dimension];
+ Vector worst = vectors[0];
+ Vector second_worst = vectors[1];
+ // reflect
+ Vector reflected = cog + (cog - worst)*alpha;
+ if (f(reflected) > f(second_worst) && f(reflected) < f(best)) {
+ vectors[0] = reflected;
+ } else if (f(reflected) > f(best)) {
+ // expand
+ Vector expanded = cog + (cog - worst)*gamma;
+ if (f(expanded) > f(reflected)) {
+ vectors[0] = expanded;
+ } else {
+ vectors[0] = reflected;
+ }
+ } else {
+ // contract
+ Vector contracted = cog + (cog - worst)*rho;
+ if (f(contracted) > f(worst)) {
+ vectors[0] = contracted;
+ } else {
+ for (int i=0; i<dimension; i++) {
+ vectors[i] = best + (vectors[i] - best)*sigma;
+ }
+ }
+ }
+ }
+
+ // algorithm is terminating, output: simplex' center of gravity
+ Vector cog;
+ for (int i = 0; i<=dimension; i++) {
+ cog += vectors[i];
+ }
+ return cog/(dimension+1);
+ } else {
+ // as long as we don't have enough vectors, request random ones,
+ // with coordinates between 0 and 1. If you want other start vectors,
+ // simply ignore these and use `step` on the vectors you want.
+ Vector result;
+ result.prepare(dimension);
+ for (int i = 0; i<dimension; ++i) {
+ result[i] = 0.001*(rand()%1000);
+ }
+ return result;
+ }
+ } catch (Vector v) {
+ return v;
+ }
+ }
+ private:
+ float f(Vector vec) {
+ return db.lookup(vec);
+ }
+ int dimension;
+ float alpha, gamma, rho, sigma;
+ float termination_distance;
+ vector<Vector> vectors;
+ ValueDB db;
+};
\ No newline at end of file
diff --git a/cpp_dnn/vector.h b/cpp_dnn/vector.h
new file mode 100644
index 0000000..5253005
--- /dev/null
+++ b/cpp_dnn/vector.h
@@ -0,0 +1,22 @@
+#ifndef CPP_DNN_VECTOR_H
+#define CPP_DNN_VECTOR_H
+
+#include <Eigen/Core>
+
+namespace cpp_dnn {
+
+//--------------------------------------------------------------------------------------------------
+template <typename T>
+using Vector2 = Eigen::Matrix<T, 2, 1>;
+
+//--------------------------------------------------------------------------------------------------
+template <typename T>
+using Vector3 = Eigen::Matrix<T, 3, 1>;
+
+//--------------------------------------------------------------------------------------------------
+template <typename T>
+using VectorX = Eigen::Matrix<T, Eigen::Dynamic, 1>;
+
+}
+
+#endif
diff --git a/cpp_frep/CMakeLists.txt b/cpp_frep/CMakeLists.txt
index 81eaa2c..d29236a 100755
--- a/cpp_frep/CMakeLists.txt
+++ b/cpp_frep/CMakeLists.txt
@@ -4,7 +4,8 @@ add_library(cpp_frep_lib
nelder_mead.h
)
find_package( OpenCV REQUIRED )
-find_package(nlohmann_json 3.2.0 REQUIRED)
+# find_package(nlohmann_json 3.2.0 REQUIRED)
+# find_package(nlohmann-json-dev REQUIRED)
target_include_directories(cpp_frep_lib PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
target_link_libraries(cpp_frep_lib shared_settings Eigen3::Eigen ${OpenCV_LIBS})
target_compile_features(cpp_frep_lib PUBLIC cxx_std_17)
@@ -12,4 +13,5 @@ target_compile_features(cpp_frep_lib PUBLIC cxx_std_17)
add_executable(cpp_frep
main.cpp
)
-target_link_libraries(cpp_frep shared_settings cpp_frep_lib nlohmann_json::nlohmann_json)
+# target_link_libraries(cpp_frep shared_settings cpp_frep_lib nlohmann_json::nlohmann_json)
+target_link_libraries(cpp_frep shared_settings cpp_frep_lib)
diff --git a/cpp_frep/main.cpp b/cpp_frep/main.cpp
index b93ae21..0ca533d 100755
--- a/cpp_frep/main.cpp
+++ b/cpp_frep/main.cpp
@@ -5,12 +5,12 @@
#include <stdio.h>
#include <stdlib.h>
#include <opencv2/opencv.hpp>
-#include <nlohmann/json.hpp>
+// #include <nlohmann/json.hpp>
using namespace cpp_frep;
-using json = nlohmann::json;
+// using json = nlohmann::json;
using scalar=float;
//--------------------------------------------------------------------------------------------------
--
GitLab