diff --git a/eo/src/es/CMAParams.cpp b/eo/src/es/CMAParams.cpp new file mode 100644 index 00000000..b5e4346f --- /dev/null +++ b/eo/src/es/CMAParams.cpp @@ -0,0 +1,234 @@ + +/* + * C++ification of Nikolaus Hansen's original C-source code for the + * CMA-ES + * + * C++-ificiation performed by Maarten Keijzer (C) 2005. Licensed under + * the LGPL. Original copyright of Nikolaus Hansen can be found below + * + * + * + */ + +/* --------------------------------------------------------- */ +/* --------------------------------------------------------- */ +/* --- File: cmaes.c -------- Author: Nikolaus Hansen --- */ +/* --------------------------------------------------------- */ +/* + * CMA-ES for non-linear function minimization. + * + * Copyright (C) 1996, 2003 Nikolaus Hansen. + * e-mail: hansen@bionik.tu-berlin.de + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later + * version (see http://www.gnu.org/copyleft/lesser.html). + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * */ +/* --- 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. + * Version 2.23. + * */ + +#include +#include + +#include + +using namespace std; + +namespace eo { + +CMAParams::CMAParams(eoParser& parser) { + + string section = "CMA parameters"; + + n = parser.createParam(0, "dimensionality", "Dimensionality (N) of the problem", 'N', section, true).value(); + + maxgen = parser.createParam( + 1000, + "max-gen", + "Maximum number of generations that the system will run (needed for damping)", + 'M', + section).value(); + + + if (n == 0) { + return; + } + + defaults(n, maxgen); + + /* handle lambda */ + lambda = parser.createParam( + lambda, + "lambda", + "Number of offspring", + 'l', + section).value(); + + if (lambda < 2) { + lambda = 4+(int)(3*log((double) n)); + cerr << "Too small lambda specified, setting it to " << lambda << endl; + } + + /* handle mu */ + mu = parser.createParam( + mu, + "mu", + "Population size", + 'm', + section).value(); + + if (mu >= lambda) { + mu = lambda/2; + cerr << "Mu set larger/equal to lambda, setting it to " << mu << endl; + } + + /* handle selection weights */ + + int weight_type = parser.createParam( + 0, + "weighting", + "Weighting scheme (for 'selection'): 0 = logarithmic, 1 = equal, 2 = linear", + 'w', + section).value(); + + switch (weight_type) { + case 1: + { + for (unsigned i = 0; i < weights.size(); ++i) { + weights[i] = mu - i; + } + } + case 2: + { + weights = 1.; + } + default : + { + for (unsigned i = 0; i < weights.size(); ++i) { + weights[i] = log(mu+1.)-log(i+1.); + } + } + + } + + /* Normalize weights and set mu_eff */ + double sumw = weights.sum(); + mueff = sumw * sumw / (weights * weights).sum(); + weights /= sumw; + + + /* most of the rest depends on mu_eff, so needs to be set again */ + + /* set the others using Nikolaus logic. If you want to tweak, you can parameterize over these defaults */ + mucov = mueff; + ccumsig = (mueff + 2.) / (n + mueff + 3.); + ccumcov = 4. / (n + 4); + + double t1 = 2. / ((n+1.4142)*(n+1.4142)); + double t2 = (2.*mucov-1.) / ((n+2.)*(n+2.)+mucov); + t2 = (t2 > 1) ? 1 : t2; + t2 = (1./mucov) * t1 + (1.-1./mucov) * t2; + + ccov = t2; + + damp = 1 + std::max(0.3,(1.-(double)n/(double)maxgen)) + * (1+2*std::max(0.,sqrt((mueff-1.)/(n+1.))-1)) /* limit sigma increase */ + / ccumsig; + + vector mins(1,0.0); + mins = parser.createParam( + mins, + "min-stdev", + "Array of minimum stdevs, last one will apply for all remaining axes", + 0, + section).value(); + + if (mins.size() > n) mins.resize(n); + + if (mins.size()) { + minStdevs = mins.back(); + for (unsigned i = 0; i < mins.size(); ++i) { + minStdevs[i] = mins[i]; + } + } + + vector inits(1,0.3); + inits = parser.createParam( + inits, + "init-stdev", + "Array of initial stdevs, last one will apply for all remaining axes", + 0, + section).value(); + + if (inits.size() > n) inits.resize(n); + + if (inits.size()) { + initialStdevs = inits.back(); + for (unsigned i = 0; i < inits.size(); ++i) { + initialStdevs[i] = inits[i]; + } + } + +} + +void CMAParams::defaults(unsigned n_, unsigned maxgen_) { + n = n_; + maxgen = maxgen_; + + lambda = 4+(int)(3*log((double) n)); + mu = lambda / 2; + + weights.resize(mu); + + for (unsigned i = 0; i < weights.size(); ++i) { + weights[i] = log(mu+1.)-log(i+1.); + } + + /* Normalize weights and set mu_eff */ + double sumw = weights.sum(); + mueff = sumw * sumw / (weights * weights).sum(); + weights /= sumw; + + mucov = mueff; + ccumsig *= (mueff + 2.) / (n + mueff + 3.); + ccumcov = 4. / (n + 4); + + double t1 = 2. / ((n+1.4142)*(n+1.4142)); + double t2 = (2.*mucov-1.) / ((n+2.)*(n+2.)+mucov); + t2 = (t2 > 1) ? 1 : t2; + t2 = (1./mucov) * t1 + (1.-1./mucov) * t2; + + ccov = t2; + + damp = 1 + std::max(0.3,(1.-(double)n/maxgen)) + * (1+2*std::max(0.,sqrt((mueff-1.)/(n+1.))-1)) /* limit sigma increase */ + / ccumsig; + + minStdevs.resize(n); + minStdevs = 0.0; + + initialStdevs.resize(n); + initialStdevs = 0.3; + + +} + + +}// namespace eo diff --git a/eo/src/es/CMAParams.h b/eo/src/es/CMAParams.h new file mode 100644 index 00000000..c5b32f7b --- /dev/null +++ b/eo/src/es/CMAParams.h @@ -0,0 +1,52 @@ +/* + * C++ification of Nikolaus Hansen's original C-source code for the + * CMA-ES. + * + * Copyright (C) 1996, 2003, Nikolaus Hansen + * (C) 2005, Maarten Keijzer + * + * License: LGPL + * + */ + +#ifndef CMAPARAMS_H__ +#define CMAPARAMS_H__ + +#include + +class eoParser; +namespace eo { + +class CMAParams { + + public: + + CMAParams() { /* Call this and all values need to be set by hand */ } + CMAParams(eoParser& parser); + + void defaults(unsigned n_, unsigned maxgen_); /* apply all defaults using n and maxgen */ + + unsigned n; + unsigned maxgen; + + unsigned lambda; /* -> mu */ + unsigned mu; /* -> weights, lambda */ + + std::valarray weights; /* <- mu, -> mueff -> mucov -> ccov */ + double mueff; /* <- weights */ + + double mucov; + + double damp; /* <- ccumsig, maxeval, lambda */ + double ccumsig; /* -> damp, <- N */ + double ccumcov; + double ccov; /* <- mucov, N */ + + std::valarray minStdevs; /* Minimum standard deviations per coordinate (default = 0.0) */ + std::valarray initialStdevs; /* Initial standard deviations per coordinate (default = 0.3) */ +}; + +} // namespace eo + +#endif + diff --git a/eo/src/es/CMAState.cpp b/eo/src/es/CMAState.cpp new file mode 100644 index 00000000..f44fdf5d --- /dev/null +++ b/eo/src/es/CMAState.cpp @@ -0,0 +1,365 @@ +/* + * C++ification of Nikolaus Hansen's original C-source code for the + * CMA-ES + * + * C++-ificiation performed by Maarten Keijzer (C) 2005. Licensed under + * the LGPL. Original copyright of Nikolaus Hansen can be found below + * + * + * Some changes have been made to the original flow and logic of the + * algorithm: + * + * - Numerical issues are now treated 'before' going into the eigenvector decomposition + * (this was done out of convenience) + * - dMaxSignifiKond (smallest x such that x == x + 1.0) replaced by + * numeric_limits::epsilon() (smallest x such that 1.0 != 1.0 + x) + * + * + */ + +/* --------------------------------------------------------- */ +/* --------------------------------------------------------- */ +/* --- File: cmaes.c -------- Author: Nikolaus Hansen --- */ +/* --------------------------------------------------------- */ +/* + * CMA-ES for non-linear function minimization. + * + * Copyright (C) 1996, 2003 Nikolaus Hansen. + * e-mail: hansen@bionik.tu-berlin.de + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later + * version (see http://www.gnu.org/copyleft/lesser.html). + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * */ +/* --- 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. + * Version 2.23. + * */ + +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +using namespace std; + +namespace eo { + +struct CMAStateImpl { + + CMAParams p; + + lower_triangular_matrix C; // Covariance matrix + square_matrix B; // Eigen vectors (in columns) + valarray d; // eigen values (diagonal matrix) + valarray pc; // Evolution path + valarray ps; // Evolution path for stepsize; + + vector mean; // current mean to sample around + double sigma; // global step size + + unsigned gen; + vector fitnessHistory; + + + CMAStateImpl(const CMAParams& params_, const vector& m, double sigma_) : + p(params_), + C(p.n), B(p.n), d(p.n), pc(p.n), ps(p.n), mean(m), sigma(sigma_), + gen(0), fitnessHistory(3) + { + double trace = (p.initialStdevs * p.initialStdevs).sum(); + /* Initialize covariance structure */ + for (unsigned i = 0; i < p.n; ++i) + { + B[i][i] = 1.; + d[i] = p.initialStdevs[i] * sqrt(p.n / trace); + C[i][i] = d[i] * d[i]; + pc[i] = 0.; + ps[i] = 0.; + } + + } + + void sample(vector& v) { + unsigned n = p.n; + v.resize(n); + + vector tmp(n); + for (unsigned i = 0; i < n; ++i) + tmp[i] = d[i] * rng.normal(); + + /* add mutation (sigma * B * (D*z)) */ + for (unsigned i = 0; i < n; ++i) { + double sum = 0; + for (unsigned j = 0; j < n; ++j) { + sum += B[i][j] * tmp[j]; + } + v[i] = mean[i] + sigma * sum; + } + } + + void reestimate(const vector* >& pop, double muBest, double muWorst) { + + assert(pop.size() == p.mu); + + unsigned n = p.n; + + fitnessHistory[gen % fitnessHistory.size()] = muBest; // needed for divergence check + + vector oldmean = mean; + valarray BDz(n); + + /* calculate xmean and rgBDz~N(0,C) */ + for (unsigned i = 0; i < n; ++i) { + mean[i] = 0.; + for (unsigned j = 0; j < pop.size(); ++j) { + mean[i] += p.weights[j] * (*pop[j])[i]; + } + BDz[i] = sqrt(p.mueff)*(mean[i] - oldmean[i])/sigma; + } + + vector tmp(n); + /* calculate z := D^(-1) * B^(-1) * rgBDz into rgdTmp */ + for (unsigned i = 0; i < n; ++i) { + double sum = 0.0; + for (unsigned j = 0; j < n; ++j) { + sum += B[j][i] * BDz[j]; + } + tmp[i] = sum / d[i]; + } + + /* cumulation for sigma (ps) using B*z */ + for (unsigned i = 0; i < n; ++i) { + double sum = 0.0; + for (unsigned j = 0; j < n; ++j) + sum += B[i][j] * tmp[j]; + + ps[i] = (1. - p.ccumsig) * ps[i] + sqrt(p.ccumsig * (2. - p.ccumsig)) * sum; + } + + /* calculate norm(ps)^2 */ + double psxps = (ps * ps).sum(); + + + double chiN = sqrt((double) p.n) * (1. - 1./(4.*p.n) + 1./(21.*p.n*p.n)); + /* cumulation for covariance matrix (pc) using B*D*z~N(0,C) */ + double hsig = sqrt(psxps) / sqrt(1. - pow(1.-p.ccumsig, 2.*gen)) / chiN < 1.5 + 1./(p.n-0.5); + + pc = (1. - p.ccumcov) * pc + hsig * sqrt(p.ccumcov * (2. - p.ccumcov)) * BDz; + + /* stop initial phase (MK, this was not reachable in the org code, deleted) */ + + /* remove momentum in ps, if ps is large and fitness is getting worse */ + + if (gen >= fitnessHistory.size()) { + + // find direction from muBest and muWorst (muBest == muWorst handled seperately + double direction = muBest < muWorst? -1.0 : 1.0; + + unsigned now = gen % fitnessHistory.size(); + unsigned prev = (gen-1) % fitnessHistory.size(); + unsigned prevprev = (gen-2) % fitnessHistory.size(); + + bool fitnessWorsens = (muBest == muWorst) || // <- increase norm also when population has converged (this deviates from Hansen's scheme) + ( (direction * fitnessHistory[now] < direction * fitnessHistory[prev]) + && + (direction * fitnessHistory[now] < direction * fitnessHistory[prevprev])); + + if(psxps/p.n > 1.5 + 10.*sqrt(2./p.n) && fitnessWorsens) { + double tfac = sqrt((1 + std::max(0., log(psxps/p.n))) * p.n / psxps); + ps *= tfac; + psxps *= tfac*tfac; + } + } + + /* update of C */ + /* Adapt_C(t); not used anymore */ + if (p.ccov != 0.) { + //flgEigensysIsUptodate = 0; + + /* update covariance matrix */ + for (unsigned i = 0; i < n; ++i) { + vector::iterator c_row = C[i]; + for (unsigned j = 0; j <= i; ++j) { + c_row[j] = + (1 - p.ccov) * c_row[j] + + + p.ccov * (1./p.mucov) * pc[i] * pc[j] + + + (1-hsig) * p.ccumcov * (2. - p.ccumcov) * c_row[j]; + + /*C[i][j] = (1 - p.ccov) * C[i][j] + + sp.ccov * (1./sp.mucov) + * (rgpc[i] * rgpc[j] + + (1-hsig)*sp.ccumcov*(2.-sp.ccumcov) * C[i][j]); */ + for (unsigned k = 0; k < p.mu; ++k) { /* additional rank mu update */ + c_row[j] += p.ccov * (1-1./p.mucov) * p.weights[k] + * ( (*pop[k])[i] - oldmean[i]) + * ( (*pop[k])[j] - oldmean[j]) + / sigma / sigma; + + // * (rgrgx[index[k]][i] - rgxold[i]) + // * (rgrgx[index[k]][j] - rgxold[j]) + // / sigma / sigma; + } + } + } + } + + /* update of sigma */ + sigma *= exp(((sqrt(psxps)/chiN)-1.)/p.damp); + /* calculate eigensystem, must be done by caller */ + //cmaes_UpdateEigensystem(0); + + + /* treat minimal standard deviations and numeric problems + * Note that in contrast with the original code, some numerical issues are treated *before* we + * go into the eigenvalue calculation */ + + treatNumericalIssues(muBest, muWorst); + + gen++; // increase generation + } + + void treatNumericalIssues(double best, double worst) { + + /* treat stdevs */ + for (unsigned i = 0; i < p.n; ++i) { + if (sigma * sqrt(C[i][i]) < p.minStdevs[i]) { + // increase stdev + sigma *= exp(0.05+1./p.damp); + break; + } + } + + /* treat convergence */ + if (best == worst) { + sigma *= exp(0.2 + 1./p.damp); + } + + /* Jede Hauptachse i testen, ob x == x + 0.1 * sigma * rgD[i] * B[i] */ + /* Test if all the means are not numerically out of whack with our coordinate system*/ + for (unsigned axis = 0; axis < p.n; ++axis) { + double fac = 0.1 * sigma * d[axis]; + unsigned coord; + for (coord = 0; coord < p.n; ++coord) { + if (mean[coord] != mean[coord] + fac * B[coord][axis]) { + break; + } + } + + if (coord == p.n) { // means are way too big (little) for numerical accuraccy. Start rocking the craddle a bit more + sigma *= exp(0.2+1./p.damp); + } + + } + + /* Testen ob eine Komponente des Objektparameters festhaengt */ + /* Correct issues with scale between objective values and covariances */ + bool theresAnIssue = false; + + for (unsigned i = 0; i < p.n; ++i) { + if (mean[i] == mean[i] + 0.2 * sigma * sqrt(C[i][i])) { + C[i][i] *= (1. + p.ccov); + theresAnIssue = true; + } + } + + if (theresAnIssue) { + sigma *= exp(0.05 + 1./p.damp); + } + } + + + bool updateEigenSystem(unsigned max_tries, unsigned max_iters) { + + if (max_iters==0) max_iters = 30 * p.n; + + static double lastGoodMinimumEigenValue = 1.0; + + /* Try to get a valid calculation */ + for (unsigned tries = 0; tries < max_tries; ++tries) { + + unsigned iters = eig( p.n, C, d, B, max_iters); + if (iters < max_iters) + { // all is well + + /* find largest/smallest eigenvalues */ + double minEV = d.min(); + double maxEV = d.max(); + + /* (MK Original comment was) :Limit Condition of C to dMaxSignifKond+1 + * replaced dMaxSignifKond with 1./numeric_limits::epsilon() + * */ + if (maxEV * numeric_limits::epsilon() > minEV) { + double tmp = maxEV * numeric_limits::epsilon() - minEV; + minEV += tmp; + for (unsigned i=0;i& initial_point, const double initial_sigma) + : pimpl(new CMAStateImpl(params, initial_point, initial_sigma)) {} + +CMAState::~CMAState() { delete pimpl; } +CMAState::CMAState(const CMAState& that) : pimpl(new CMAStateImpl(*that.pimpl )) {} +CMAState& CMAState::operator=(const CMAState& that) { *pimpl = *that.pimpl; return *this; } + +void CMAState::sample(vector& v) const { pimpl->sample(v); } + +void CMAState::reestimate(const vector* >& population, double muBest, double muWorst) { pimpl->reestimate(population, muBest, muWorst); } +bool CMAState::updateEigenSystem(unsigned max_tries, unsigned max_iters) { return pimpl->updateEigenSystem(max_tries, max_iters); } + + +} // namespace eo + diff --git a/eo/src/es/CMAState.h b/eo/src/es/CMAState.h new file mode 100644 index 00000000..505a31b7 --- /dev/null +++ b/eo/src/es/CMAState.h @@ -0,0 +1,79 @@ + +/* + * C++ification of Nikolaus Hansen's original C-source code for the + * CMA-ES. + * + * Copyright (C) 1996, 2003, Nikolaus Hansen + * (C) 2005, Maarten Keijzer + * + * License: LGPL (see source file) + * + */ + +#ifndef CMASTATE_H_ +#define CMASTATE_H_ + +#include +#include + +namespace eo { + + +class CMAStateImpl; +class CMAParams; +class CMAState { + + CMAStateImpl* pimpl; /* pointer to implementation, hidden in source file */ + + public: + + CMAState(const CMAParams&, const std::vector& initial_point, const double initial_sigma = 1.0); + ~CMAState(); + CMAState(const CMAState&); + CMAState& operator=(const CMAState&); + + /** + * sample a vector from the distribution + * + * If the sample is not to your liking (i.e., not within bounds) + * you can do one of two things: + * + * a) Call sample again + * b) multiply the entire vector with a number between -1 and 1 + * + * Do not modify the sample in any other way as this will invalidate the + * internal consistency of the system. + * + * A final approach is to copy the sample and modify the copy externally (in the evaluation function) + * and possibly add a penalty depending on the size of the modification. + * + */ + void sample(std::vector& v) const; + + /** + * Reestimate covariance matrix and other internal parameters + * Does NOT update the eigen system (call that seperately) + * + * Needs a population of mu individuals, sorted on fitness, plus + * + * muBest: the best fitness in the population + * muWorst: the worst fitness in the population + */ + + void reestimate(const std::vector* >& sorted_population, double muBest, double muWorst); + + /** + * call this function after reestimate in order to update the eigen system + * It is a seperate call to allow the user to periodically skip this expensive step + * + * max_iters = 0 implies 30 * N iterations + * + * If after max_tries still no numerically sound eigen system is constructed, + * the function returns false + */ + bool updateEigenSystem(unsigned max_tries = 1, unsigned max_iters = 0); +}; + +} // namespace eo + +#endif diff --git a/eo/src/es/Makefile.am b/eo/src/es/Makefile.am index 2fc5ecb0..6a4b827d 100644 --- a/eo/src/es/Makefile.am +++ b/eo/src/es/Makefile.am @@ -1,6 +1,6 @@ ## Makefile.am for eo/src/es -lib_LIBRARIES = libes.a +lib_LIBRARIES = libes.a libcma.a libes_a_SOURCES = make_algo_scalar_es.cpp \ make_algo_scalar_real.cpp \ @@ -40,5 +40,15 @@ esinclude_HEADERS = eoEsChromInit.h \ make_op_real.h \ make_real.h +libcma_a_SOURCES = eig.cpp \ + CMAState.cpp \ + CMAParams.cpp + +cmaincludedir = $(pkgincludedir)/es + +cmainclude_HEADERS = eig.h \ + CMAState.h \ + matrices.h \ + CMAParams.h AM_CXXFLAGS = -I$(top_srcdir)/src diff --git a/eo/src/es/eig.cpp b/eo/src/es/eig.cpp new file mode 100644 index 00000000..29998bcb --- /dev/null +++ b/eo/src/es/eig.cpp @@ -0,0 +1,259 @@ + +/* + * C++ification of Nikolaus Hansen's original C-source code for the + * CMA-ES. These are the eigenvector calculations + * + * C++-ificiation performed by Maarten Keijzer (C) 2005. Licensed under + * the LGPL. Original copyright of Nikolaus Hansen can be found below + * + * This algorithm is held almost completely intact. Some other datatypes are used, + * but hardly any code has changed + * + */ + +/* --------------------------------------------------------- */ +/* --------------------------------------------------------- */ +/* --- File: cmaes.c -------- Author: Nikolaus Hansen --- */ +/* --------------------------------------------------------- */ +/* + * CMA-ES for non-linear function minimization. + * + * Copyright (C) 1996, 2003 Nikolaus Hansen. + * e-mail: hansen@bionik.tu-berlin.de + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later + * version (see http://www.gnu.org/copyleft/lesser.html). + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * */ +/* --- 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. + * Version 2.23. + * */ +#include "eig.h" + +using namespace std; + +/* ========================================================= */ +/* + Householder Transformation einer symmetrischen Matrix + auf tridiagonale Form. + -> n : Dimension + -> ma : symmetrische nxn-Matrix + <- ma : Transformationsmatrix (ist orthogonal): + Tridiag.-Matrix == <-ma * ->ma * (<-ma)^t + <- diag : Diagonale der resultierenden Tridiagonalmatrix + <- neben[0..n-1] : Nebendiagonale (==1..n-1) der res. Tridiagonalmatrix + + */ +static void +Householder( int N, square_matrix& ma, valarray& diag, double* neben) +{ + double epsilon; + int i, j, k; + double h, sum, tmp, tmp2; + + for (i = N-1; i > 0; --i) + { + h = 0.0; + if (i == 1) + neben[i] = ma[i][i-1]; + else + { + for (k = i-1, epsilon = 0.0; k >= 0; --k) + epsilon += fabs(ma[i][k]); + + if (epsilon == 0.0) + neben[i] = ma[i][i-1]; + else + { + for(k = i-1, sum = 0.0; k >= 0; --k) + { /* i-te Zeile von i-1 bis links normieren */ + ma[i][k] /= epsilon; + sum += ma[i][k]*ma[i][k]; + } + tmp = (ma[i][i-1] > 0) ? -sqrt(sum) : sqrt(sum); + neben[i] = epsilon*tmp; + h = sum - ma[i][i-1]*tmp; + ma[i][i-1] -= tmp; + for (j = 0, sum = 0.0; j < i; ++j) + { + ma[j][i] = ma[i][j]/h; + tmp = 0.0; + for (k = j; k >= 0; --k) + tmp += ma[j][k]*ma[i][k]; + for (k = j+1; k < i; ++k) + tmp += ma[k][j]*ma[i][k]; + neben[j] = tmp/h; + sum += neben[j] * ma[i][j]; + } /* for j */ + sum /= 2.*h; + for (j = 0; j < i; ++j) + { + neben[j] -= ma[i][j]*sum; + tmp = ma[i][j]; + tmp2 = neben[j]; + for (k = j; k >= 0; --k) + ma[j][k] -= (tmp*neben[k] + tmp2*ma[i][k]); + } /* for j */ + } /* else epsilon */ + } /* else i == 1 */ + diag[i] = h; + } /* for i */ + + diag[0] = 0.0; + neben[0] = 0.0; + + for (i = 0; i < N; ++i) + { + if(diag[i] != 0.0) + for (j = 0; j < i; ++j) + { + for (k = i-1, tmp = 0.0; k >= 0; --k) + tmp += ma[i][k] * ma[k][j]; + for (k = i-1; k >= 0; --k) + ma[k][j] -= tmp*ma[k][i]; + } /* for j */ + diag[i] = ma[i][i]; + ma[i][i] = 1.0; + for (k = i-1; k >= 0; --k) + ma[k][i] = ma[i][k] = 0.0; + } /* for i */ +} + +/* + QL-Algorithmus mit implizitem Shift, zur Berechnung von Eigenwerten + und -vektoren einer symmetrischen Tridiagonalmatrix. + -> n : Dimension. + -> diag : Diagonale der Tridiagonalmatrix. + -> neben[0..n-1] : Nebendiagonale (==0..n-2), n-1. Eintrag beliebig + -> mq : Matrix output von Householder() + -> maxIt : maximale Zahl der Iterationen + <- diag : Eigenwerte + <- neben : Garbage + <- mq : k-te Spalte ist normalisierter Eigenvektor zu diag[k] + + */ + +static int +QLalgo( int N, valarray& diag, square_matrix& mq, + int maxIter, double* neben) +{ + int i, j, k, kp1, l; + double tmp, diff, cneben, c1, c2, p; + int iter; + + neben[N-1] = 0.0; + for (i = 0, iter = 0; i < N && iter < maxIter; ++i) + do /* while j != i */ + { + for (j = i; j < N-1; ++j) + { + tmp = fabs(diag[j]) + fabs(diag[j+1]); + if (fabs(neben[j]) + tmp == tmp) + break; + } + if (j != i) + { + if (++iter > maxIter) return maxIter-1; + diff = (diag[i+1]-diag[i])/neben[i]/2.0; + if (diff >= 0) + diff = diag[j] - diag[i] + + neben[i] / (diff + sqrt(diff * diff + 1.0)); + else + diff = diag[j] - diag[i] + + neben[i] / (diff - sqrt(diff * diff + 1.0)); + c2 = c1 = 1.0; + p = 0.0; + for (k = j-1; k >= i; --k) + { + kp1 = k + 1; + tmp = c2 * neben[k]; + cneben = c1 * neben[k]; + if (fabs(tmp) >= fabs(diff)) + { + c1 = diff / tmp; + c2 = 1. / sqrt(c1*c1 + 1.0); + neben[kp1] = tmp / c2; + c1 *= c2; + } + else + { + c2 = tmp / diff; + c1 = 1. / sqrt(c2*c2 + 1.0); + neben[kp1] = diff / c1; + c2 *= c1; + } /* else */ + tmp = (diag[k] - diag[kp1] + p) * c2 + 2.0 * c1 * cneben; + diag[kp1] += tmp * c2 - p; + p = tmp * c2; + diff = tmp * c1 - cneben; + for (l = N-1; l >= 0; --l) /* TF-Matrix Q */ + { + tmp = mq[l][kp1]; + mq[l][kp1] = c2 * mq[l][k] + c1 * tmp; + mq[l][k] = c1 * mq[l][k] - c2 * tmp; + } /* for l */ + } /* for k */ + diag[i] -= p; + neben[i] = diff; + neben[j] = 0.0; + } /* if */ + } while (j != i); + return iter; +} /* QLalgo() */ + +/* ========================================================= */ +/* + Calculating eigenvalues and vectors. + Input: + N: dimension. + C: lower_triangular NxN-matrix. + 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. + return: number of iterations in QL-Algorithm. + */ + +namespace eo { +int +eig( int N, const lower_triangular_matrix& C, valarray& diag, square_matrix& Q, + int niter) +{ + int ret; + int i, j; + + if (niter == 0) niter = 30*N; + + for (i=0; i < N; ++i) + { + vector::const_iterator row = C[i]; + for (j = 0; j <= i; ++j) + Q[i][j] = Q[j][i] = row[j]; + } + + double* rgtmp = new double[N+1]; + Householder( N, Q, diag, rgtmp); + ret = QLalgo( N, diag, Q, niter, rgtmp+1); + delete [] rgtmp; + + return ret; +} + +} // namespace eo diff --git a/eo/src/es/eig.h b/eo/src/es/eig.h new file mode 100644 index 00000000..02aa36e0 --- /dev/null +++ b/eo/src/es/eig.h @@ -0,0 +1,25 @@ +#ifndef EIG_H__ +#define EIG_H__ + +#include +#include + +namespace eo { +/* ========================================================= */ +/* + Calculating eigenvalues and vectors. + Input: + N: dimension. + C: lower_triangular NxN-matrix. + niter: number of maximal iterations for QL-Algorithm. + Output: + diag: N eigenvalues. + Q: Columns are normalized eigenvectors. + return: number of iterations in QL-Algorithm. + */ +extern int eig( int N, const lower_triangular_matrix& C, std::valarray& diag, square_matrix& Q, + int niter = 0); + +} // namespace eo + +#endif diff --git a/eo/src/es/eoCMABreed.h b/eo/src/es/eoCMABreed.h new file mode 100644 index 00000000..3add8c9e --- /dev/null +++ b/eo/src/es/eoCMABreed.h @@ -0,0 +1,76 @@ +// -*- mode: c++; c-indent-level: 4; c++-member-init-indent: 8; comment-column: 35; fill-column: 80; -*- + +//----------------------------------------------------------------------------- +// eoCMABreed +// (c) Maarten Keijzer 2005 +/* + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + */ +//----------------------------------------------------------------------------- + + +#ifndef _EOCMABREED_H +#define _EOCMABREED_H + +#include +#include +#include + +#include + +/// TODO, handle bounds +template +class eoCMABreed : public eoBreed< eoVector > { + + eo::CMAState& state; + unsigned lambda; + + typedef eoVector EOT; + + public: + eoCMABreed(eo::CMAState& state_, unsigned lambda_) : state(state_), lambda(lambda_) {} + + void operator()(const eoPop& parents, eoPop& offspring) { + + // two temporary arrays of pointers to store the sorted population + std::vector sorted(parents.size()); + std::vector* > mu(parents.size()); + + parents.sort(sorted); + for (unsigned i = 0; i < sorted.size(); ++i) { + mu[i] = static_cast< const std::vector* >( sorted[i] ); + } + + // learn + + state.reestimate(mu, sorted[0]->fitness(), sorted.back()->fitness()); + + if (!state.updateEigenSystem(10)) { + std::cerr << "No good eigensystem found" << std::endl; + } + + // generate + offspring.resize(lambda); + + for (unsigned i = 0; i < lambda; ++i) { + state.sample( static_cast< std::vector& >( offspring[i] )); + } + + } +}; + + +#endif diff --git a/eo/src/es/eoCMAInit.h b/eo/src/es/eoCMAInit.h new file mode 100644 index 00000000..b605cd99 --- /dev/null +++ b/eo/src/es/eoCMAInit.h @@ -0,0 +1,54 @@ +// -*- mode: c++; c-indent-level: 4; c++-member-init-indent: 8; comment-column: 35; fill-column: 80; -*- + +//----------------------------------------------------------------------------- +// eoCMAInit +// (c) Maarten Keijzer 2005 +/* + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + Contact: todos@geneura.ugr.es, http://geneura.ugr.es + marc.schoenauer@polytechnique.fr + http://eeaax.cmap.polytchnique.fr/ + */ +//----------------------------------------------------------------------------- + + +#ifndef _EOCMAINIT_H +#define _EOCMAINIT_H + +#include +#include +#include + +/// TODO, handle bounds +template +class eoCMAInit : public eoInit< eoVector > { + + const eo::CMAState& state; + + typedef eoVector EOT; + + public: + eoCMAInit(const eo::CMAState& state_) : state(state_) {} + + + void operator()(EOT& v) { + state.sample(static_cast& >(v)); + v.invalidate(); + } +}; + + +#endif diff --git a/eo/src/es/make_algo_scalar_es.cpp b/eo/src/es/make_algo_scalar_es.cpp index abc52682..fa35cb36 100644 --- a/eo/src/es/make_algo_scalar_es.cpp +++ b/eo/src/es/make_algo_scalar_es.cpp @@ -53,7 +53,7 @@ /////// eoAlgo >& make_algo_scalar(eoParser& _parser, eoState& _state, eoEvalFunc >& _eval, eoContinue >& _continue, eoGenOp >& _op, eoDistance >* _dist) { - return do_make_algo_scalar(_parser, _state, _eval, _continue, _op); + return do_make_algo_scalar(_parser, _state, _eval, _continue, _op); } eoAlgo >& make_algo_scalar(eoParser& _parser, eoState& _state, eoEvalFunc >& _eval, eoContinue >& _continue, eoGenOp >& _op, eoDistance >* _dist) diff --git a/eo/src/es/matrices.h b/eo/src/es/matrices.h new file mode 100644 index 00000000..f30bd5c9 --- /dev/null +++ b/eo/src/es/matrices.h @@ -0,0 +1,42 @@ +#ifndef MATRICES_H +#define MATRICES_H + +#include + +class lower_triangular_matrix { + + unsigned n; + std::vector data; + + public: + + lower_triangular_matrix(unsigned n_ = 0) : n(n_), data(n * (n+1) / 2) {}; + + void resize(unsigned n_) { + n = n_; + data.resize(n*(n+1)/2); + } + + std::vector::iterator operator[](unsigned i) { return data.begin() + i * (i+1) / 2; } + std::vector::const_iterator operator[](unsigned i) const { return data.begin() + i*(i+1)/2; } +}; + +class square_matrix { + unsigned n; + std::vector data; + + public: + + square_matrix(unsigned n_ = 0) : n(n_), data(n * n) {}; + + void resize(unsigned n_) { + n = n_; + data.resize(n*n); + } + + std::vector::iterator operator[](unsigned i) { return data.begin() + i * n; } + std::vector::const_iterator operator[](unsigned i) const { return data.begin() + i*n; } +}; + +#endif +