00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054 #include <valarray>
00055 #include <limits>
00056 #include <iostream>
00057 #include <cassert>
00058
00059 #include <utils/eoRNG.h>
00060
00061 #include <es/CMAState.h>
00062 #include <es/CMAParams.h>
00063 #include <es/matrices.h>
00064 #include <es/eig.h>
00065
00066 using namespace std;
00067
00068 namespace eo {
00069
00070 struct CMAStateImpl {
00071
00072 CMAParams p;
00073
00074 lower_triangular_matrix C;
00075 square_matrix B;
00076 valarray<double> d;
00077 valarray<double> pc;
00078 valarray<double> ps;
00079
00080 vector<double> mean;
00081 double sigma;
00082
00083 unsigned gen;
00084 vector<double> fitnessHistory;
00085
00086
00087 CMAStateImpl(const CMAParams& params_, const vector<double>& m, double sigma_) :
00088 p(params_),
00089 C(p.n), B(p.n), d(p.n), pc(p.n), ps(p.n), mean(m), sigma(sigma_),
00090 gen(0), fitnessHistory(3)
00091 {
00092 double trace = (p.initialStdevs * p.initialStdevs).sum();
00093
00094 for (unsigned i = 0; i < p.n; ++i)
00095 {
00096 B[i][i] = 1.;
00097 d[i] = p.initialStdevs[i] * sqrt(p.n / trace);
00098 C[i][i] = d[i] * d[i];
00099 pc[i] = 0.;
00100 ps[i] = 0.;
00101 }
00102
00103 }
00104
00105 void sample(vector<double>& v) {
00106 unsigned n = p.n;
00107 v.resize(n);
00108
00109 vector<double> tmp(n);
00110 for (unsigned i = 0; i < n; ++i)
00111 tmp[i] = d[i] * rng.normal();
00112
00113
00114 for (unsigned i = 0; i < n; ++i) {
00115 double sum = 0;
00116 for (unsigned j = 0; j < n; ++j) {
00117 sum += B[i][j] * tmp[j];
00118 }
00119 v[i] = mean[i] + sigma * sum;
00120 }
00121 }
00122
00123 void reestimate(const vector<const vector<double>* >& pop, double muBest, double muWorst) {
00124
00125 assert(pop.size() == p.mu);
00126
00127 unsigned n = p.n;
00128
00129 fitnessHistory[gen % fitnessHistory.size()] = muBest;
00130
00131 vector<double> oldmean = mean;
00132 valarray<double> BDz(n);
00133
00134
00135 for (unsigned i = 0; i < n; ++i) {
00136 mean[i] = 0.;
00137 for (unsigned j = 0; j < pop.size(); ++j) {
00138 mean[i] += p.weights[j] * (*pop[j])[i];
00139 }
00140 BDz[i] = sqrt(p.mueff)*(mean[i] - oldmean[i])/sigma;
00141 }
00142
00143 vector<double> tmp(n);
00144
00145 for (unsigned i = 0; i < n; ++i) {
00146 double sum = 0.0;
00147 for (unsigned j = 0; j < n; ++j) {
00148 sum += B[j][i] * BDz[j];
00149 }
00150 tmp[i] = sum / d[i];
00151 }
00152
00153
00154 for (unsigned i = 0; i < n; ++i) {
00155 double sum = 0.0;
00156 for (unsigned j = 0; j < n; ++j)
00157 sum += B[i][j] * tmp[j];
00158
00159 ps[i] = (1. - p.ccumsig) * ps[i] + sqrt(p.ccumsig * (2. - p.ccumsig)) * sum;
00160 }
00161
00162
00163 double psxps = (ps * ps).sum();
00164
00165
00166 double chiN = sqrt((double) p.n) * (1. - 1./(4.*p.n) + 1./(21.*p.n*p.n));
00167
00168 double hsig = sqrt(psxps) / sqrt(1. - pow(1.-p.ccumsig, 2.*gen)) / chiN < 1.5 + 1./(p.n-0.5);
00169
00170 pc = (1. - p.ccumcov) * pc + hsig * sqrt(p.ccumcov * (2. - p.ccumcov)) * BDz;
00171
00172
00173
00174
00175
00176 if (gen >= fitnessHistory.size()) {
00177
00178
00179 double direction = muBest < muWorst? -1.0 : 1.0;
00180
00181 unsigned now = gen % fitnessHistory.size();
00182 unsigned prev = (gen-1) % fitnessHistory.size();
00183 unsigned prevprev = (gen-2) % fitnessHistory.size();
00184
00185 bool fitnessWorsens = (muBest == muWorst) ||
00186 ( (direction * fitnessHistory[now] < direction * fitnessHistory[prev])
00187 &&
00188 (direction * fitnessHistory[now] < direction * fitnessHistory[prevprev]));
00189
00190 if(psxps/p.n > 1.5 + 10.*sqrt(2./p.n) && fitnessWorsens) {
00191 double tfac = sqrt((1 + std::max(0., log(psxps/p.n))) * p.n / psxps);
00192 ps *= tfac;
00193 psxps *= tfac*tfac;
00194 }
00195 }
00196
00197
00198
00199 if (p.ccov != 0.) {
00200
00201
00202
00203 for (unsigned i = 0; i < n; ++i) {
00204 vector<double>::iterator c_row = C[i];
00205 for (unsigned j = 0; j <= i; ++j) {
00206 c_row[j] =
00207 (1 - p.ccov) * c_row[j]
00208 +
00209 p.ccov * (1./p.mucov) * pc[i] * pc[j]
00210 +
00211 (1-hsig) * p.ccumcov * (2. - p.ccumcov) * c_row[j];
00212
00213
00214
00215
00216
00217 for (unsigned k = 0; k < p.mu; ++k) {
00218 c_row[j] += p.ccov * (1-1./p.mucov) * p.weights[k]
00219 * ( (*pop[k])[i] - oldmean[i])
00220 * ( (*pop[k])[j] - oldmean[j])
00221 / sigma / sigma;
00222
00223
00224
00225
00226 }
00227 }
00228 }
00229 }
00230
00231
00232 sigma *= exp(((sqrt(psxps)/chiN)-1.)/p.damp);
00233
00234
00235
00236
00237
00238
00239
00240
00241 treatNumericalIssues(muBest, muWorst);
00242
00243 gen++;
00244 }
00245
00246 void treatNumericalIssues(double best, double worst) {
00247
00248
00249 for (unsigned i = 0; i < p.n; ++i) {
00250 if (sigma * sqrt(C[i][i]) < p.minStdevs[i]) {
00251
00252 sigma *= exp(0.05+1./p.damp);
00253 break;
00254 }
00255 }
00256
00257
00258 if (best == worst) {
00259 sigma *= exp(0.2 + 1./p.damp);
00260 }
00261
00262
00263
00264 for (unsigned axis = 0; axis < p.n; ++axis) {
00265 double fac = 0.1 * sigma * d[axis];
00266 unsigned coord;
00267 for (coord = 0; coord < p.n; ++coord) {
00268 if (mean[coord] != mean[coord] + fac * B[coord][axis]) {
00269 break;
00270 }
00271 }
00272
00273 if (coord == p.n) {
00274 sigma *= exp(0.2+1./p.damp);
00275 }
00276
00277 }
00278
00279
00280
00281 bool theresAnIssue = false;
00282
00283 for (unsigned i = 0; i < p.n; ++i) {
00284 if (mean[i] == mean[i] + 0.2 * sigma * sqrt(C[i][i])) {
00285 C[i][i] *= (1. + p.ccov);
00286 theresAnIssue = true;
00287 }
00288 }
00289
00290 if (theresAnIssue) {
00291 sigma *= exp(0.05 + 1./p.damp);
00292 }
00293 }
00294
00295
00296 bool updateEigenSystem(unsigned max_tries, unsigned max_iters) {
00297
00298 if (max_iters==0) max_iters = 30 * p.n;
00299
00300 static double lastGoodMinimumEigenValue = 1.0;
00301
00302
00303 for (unsigned tries = 0; tries < max_tries; ++tries) {
00304
00305 unsigned iters = eig( p.n, C, d, B, max_iters);
00306 if (iters < max_iters)
00307 {
00308
00309
00310 double minEV = d.min();
00311 double maxEV = d.max();
00312
00313
00314
00315
00316 if (maxEV * numeric_limits<double>::epsilon() > minEV) {
00317 double tmp = maxEV * numeric_limits<double>::epsilon() - minEV;
00318 minEV += tmp;
00319 for (unsigned i=0;i<p.n;++i) {
00320 C[i][i] += tmp;
00321 d[i] += tmp;
00322 }
00323 }
00324 lastGoodMinimumEigenValue = minEV;
00325
00326 d = sqrt(d);
00327
00328
00329
00330
00331 return true;
00332 }
00333
00334
00335
00336
00337
00338 double summand = lastGoodMinimumEigenValue * exp((double) tries);
00339 for (unsigned i = 0; i < p.n; ++i)
00340 C[i][i] += summand;
00341
00342 }
00343
00344 return false;
00345
00346 }
00347
00348
00349 };
00350
00351 CMAState::CMAState(const CMAParams& params, const std::vector<double>& initial_point, const double initial_sigma)
00352 : pimpl(new CMAStateImpl(params, initial_point, initial_sigma)) {}
00353
00354 CMAState::~CMAState() { delete pimpl; }
00355 CMAState::CMAState(const CMAState& that) : pimpl(new CMAStateImpl(*that.pimpl )) {}
00356 CMAState& CMAState::operator=(const CMAState& that) { *pimpl = *that.pimpl; return *this; }
00357
00358 void CMAState::sample(vector<double>& v) const { pimpl->sample(v); }
00359
00360 void CMAState::reestimate(const vector<const vector<double>* >& population, double muBest, double muWorst) { pimpl->reestimate(population, muBest, muWorst); }
00361 bool CMAState::updateEigenSystem(unsigned max_tries, unsigned max_iters) { return pimpl->updateEigenSystem(max_tries, max_iters); }
00362
00363
00364 }
00365