From 16f97144b3b7879323ebebd341bf926acc43d751 Mon Sep 17 00:00:00 2001 From: Johann Dreo Date: Thu, 12 Jul 2012 11:27:41 +0200 Subject: [PATCH] adaptive operators that compiles (but still not work) --- edo/src/edoEstimatorNormalAdaptive.h | 237 +++++++++++++++++++++++++++ edo/src/edoNormalAdaptive.h | 34 +++- edo/src/edoSamplerNormalAdaptive.h | 6 +- 3 files changed, 271 insertions(+), 6 deletions(-) create mode 100644 edo/src/edoEstimatorNormalAdaptive.h diff --git a/edo/src/edoEstimatorNormalAdaptive.h b/edo/src/edoEstimatorNormalAdaptive.h new file mode 100644 index 000000000..0e5a7b028 --- /dev/null +++ b/edo/src/edoEstimatorNormalAdaptive.h @@ -0,0 +1,237 @@ +/* +The Evolving Distribution Objects framework (EDO) is a template-based, +ANSI-C++ evolutionary computation library which helps you to write your +own estimation of distribution algorithms. + +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. + +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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + +Copyright (C) 2010 Thales group +*/ +/* +Authors: + Johann Dréo + Pierre Savéant +*/ + + +#ifndef _edoEstimatorNormalAdaptive_h +#define _edoEstimatorNormalAdaptive_h + +#ifdef WITH_EIGEN + +#include + +#include + +#include "edoNormalAdaptive.h" +#include "edoEstimatorAdaptive.h" + + +//! edoEstimatorNormalMulti< EOT > +template< typename EOT, typename EOD = edoNormalAdaptive > +class edoEstimatorNormalAdaptive : public edoEstimatorAdaptive< EOD > +{ +public: + typedef typename EOT::AtomType AtomType; + typedef typename EOD::Vector Vector; + typedef typename EOD::Matrix Matrix; + + edoEstimatorNormalAdaptive( EOD& distrib, unsigned int mu ) : + edoEstimatorAdaptive( distrib ), + _mu(mu), + _calls(0), + _eigeneval(0) + {} + +private: + Eigen::VectorXd edoCMAESweights( unsigned int pop_size ) + { + // compute recombination weights + Eigen::VectorXd weights( pop_size ); + double sum_w = 0; + for( unsigned int i = 0; i < _mu; ++i ) { + double w_i = log( _mu + 0.5 ) - log( i + 1 ); + weights(i) = w_i; + sum_w += w_i; + } + // normalization of weights + weights /= sum_w; + + return weights; + } + +public: + void resetCalls() + { + _calls = 0; + } + + // update the distribution reference this->distribution() + edoNormalAdaptive operator()( eoPop& pop ) + { + + /********************************************************************** + * INITIALIZATION + *********************************************************************/ + + unsigned int N = pop[0].size(); // FIXME expliciter la dimension du pb ? + unsigned int lambda = pop.size(); + + // number of calls to the operator == number of generations + _calls++; + // number of "evaluations" until now + unsigned int counteval = _calls * lambda; + + // Here, if we are in canonical CMA-ES, + // pop is supposed to be the mu ranked better solutions, + // as the rank mu selection is supposed to have occured. + Matrix arx( pop.size(), N ); + + // copy the pop (most probably a vector of vectors) in a Eigen3 matrix + for( unsigned int i = 0; i < lambda; ++i ) { + for( unsigned int d = 0; d < N; ++d ) { + arx(i,d) = pop[i][d]; + } // dimensions + } // individuals + + // muXone array for weighted recombination + Eigen::VectorXd weights = edoCMAESweights( N ); + + // FIXME exposer les constantes dans l'interface + + // variance-effectiveness of sum w_i x_i + double mueff = pow(weights.sum(), 2) / (weights.array().square()).sum(); + + // time constant for cumulation for C + double cc = (4+mueff/N) / (N+4 + 2*mueff/N); + + // t-const for cumulation for sigma control + double cs = (mueff+2) / (N+mueff+5); + + // learning rate for rank-one update of C + double c1 = 2 / (pow(N+1.3,2)+mueff); + + // and for rank-mu update + double cmu = 2 * (mueff-2+1/mueff) / ( pow(N+2,2)+mueff); + + // damping for sigma + double damps = 1 + 2*std::max(0.0, sqrt((mueff-1)/(N+1))-1) + cs; + + + // shortcut to the referenced distribution + EOD& d = this->distribution(); + + // C^-1/2 + Matrix invsqrtC = + d.coord_sys() * d.scaling().asDiagonal().inverse() + * d.coord_sys().transpose(); + + // expectation of ||N(0,I)|| == norm(randn(N,1)) + double chiN = sqrt(N)*(1-1/(4*N)+1/(21*pow(N,2))); + + + /********************************************************************** + * WEIGHTED MEAN + *********************************************************************/ + + // compute weighted mean into xmean + Vector xold = d.mean(); + d.mean( arx * weights ); + Vector xmean = d.mean(); + + + /********************************************************************** + * CUMULATION: UPDATE EVOLUTION PATHS + *********************************************************************/ + + // cumulation for sigma + d.path_sigma( + (1.0-cs)*d.path_sigma() + sqrt(cs*(2.0-cs)*mueff)*invsqrtC*(xmean-xold)/d.sigma() + ); + + // sign of h + double hsig; + if( d.path_sigma().norm()/sqrt(1.0-pow((1.0-cs),(2.0*counteval/lambda)))/chiN + < 1.4 + 2.0/(N+1.0) + ) { + hsig = 1.0; + } else { + hsig = 0.0; + } + + // cumulation for the covariance matrix + d.path_covar( + (1.0-cc)*d.path_covar() + hsig*sqrt(cc*(2.0-cc)*mueff)*(xmean-xold) / d.sigma() + ); + + Matrix artmp = (1.0/d.sigma()) * arx - xold.rowwise().replicate(_mu); + + + /********************************************************************** + * COVARIANCE MATRIX ADAPTATION + *********************************************************************/ + + d.covar( + (1-c1-cmu) * d.covar() // regard old matrix + + c1 * (d.path_covar()*d.path_covar().transpose() // plus rank one update + + (1-hsig) * cc*(2-cc) * d.covar()) // minor correction if hsig==0 + + cmu * artmp * weights.asDiagonal() * artmp.transpose() // plus rank mu update + ); + + // Adapt step size sigma + d.sigma( d.sigma() * exp((cs/damps)*(d.path_sigma().norm()/chiN - 1)) ); + + + + /********************************************************************** + * DECOMPOSITION OF THE COVARIANCE MATRIX + *********************************************************************/ + + // Decomposition of C into B*diag(D.^2)*B' (diagonalization) + if( counteval - _eigeneval > lambda/(c1+cmu)/N/10 ) { // to achieve O(N^2) + _eigeneval = counteval; + + // enforce symmetry of the covariance matrix + Matrix C = d.covar(); + // FIXME edoEstimatorNormalAdaptive.h:213:44: erreur: expected primary-expression before ‘)’ token + // copy the upper part in the lower one + //C.triangularView() = C.adjoint(); + // Matrix CS = C.triangularView() + C.triangularView().transpose(); + d.covar( C ); + + Eigen::SelfAdjointEigenSolver eigensolver( d.covar() ); // FIXME use JacobiSVD? + d.coord_sys( eigensolver.eigenvectors() ); + Matrix D = eigensolver.eigenvalues().asDiagonal(); + + // from variance to standard deviations + D.cwiseSqrt(); + d.scaling( D ); + } + + return d; + } // operator() + +protected: + + unsigned int _mu; + unsigned int _calls; + unsigned int _eigeneval; + + + // EOD & distribution() inherited from edoEstimatorAdaptive +}; +#endif // WITH_EIGEN + +#endif // !_edoEstimatorNormalAdaptive_h diff --git a/edo/src/edoNormalAdaptive.h b/edo/src/edoNormalAdaptive.h index 499af91f5..9ceb6b9fc 100644 --- a/edo/src/edoNormalAdaptive.h +++ b/edo/src/edoNormalAdaptive.h @@ -39,6 +39,7 @@ template < typename EOT > class edoNormalAdaptive : public edoDistrib< EOT > { public: + //typedef EOT EOType; typedef typename EOT::AtomType AtomType; typedef Eigen::Matrix< AtomType, Eigen::Dynamic, 1> Vector; typedef Eigen::Matrix< AtomType, Eigen::Dynamic, Eigen::Dynamic> Matrix; @@ -55,13 +56,40 @@ public: assert( dim > 0); } + edoNormalAdaptive( unsigned int dim, + Vector mean, + Matrix C, + Matrix B, + Vector D, + double sigma, + Vector p_c, + Vector p_s + ) : + _mean( mean ), + _C( C ), + _B( B ), + _D( D ), + _sigma(sigma), + _p_c( p_c ), + _p_s( p_s ) + { + assert( dim > 0); + assert( _mean.innerSize() == dim ); + assert( _C.innerSize() == dim && _C.outerSize() == dim ); + assert( _B.innerSize() == dim && _B.outerSize() == dim ); + assert( _D.innerSize() == dim ); + assert( _sigma != 0.0 ); + assert( _p_c.innerSize() == dim ); + assert( _p_s.innerSize() == dim ); + } + unsigned int size() { return _mean.innerSize(); } Vector mean() const {return _mean;} - Matrix covar() const {return _covar;} + Matrix covar() const {return _C;} Matrix coord_sys() const {return _B;} Vector scaling() const {return _D;} double sigma() const {return _sigma;} @@ -73,8 +101,8 @@ public: void coord_sys( Matrix b ) { _B = b; } void scaling( Vector d ) { _D = d; } void sigma( double s ) { _sigma = s; } - void path_covar( Vector p ) { _path_covar = p; } - void path_sigma( Vector p ) { _path_sigma = p; } + void path_covar( Vector p ) { _p_c = p; } + void path_sigma( Vector p ) { _p_s = p; } private: Vector _mean; // diff --git a/edo/src/edoSamplerNormalAdaptive.h b/edo/src/edoSamplerNormalAdaptive.h index 193044f0b..edf531eff 100644 --- a/edo/src/edoSamplerNormalAdaptive.h +++ b/edo/src/edoSamplerNormalAdaptive.h @@ -69,9 +69,9 @@ public: assert(T.innerSize() == size); assert(T.outerSize() == 1); - //Vector t_sol = distrib.mean() + distrib.sigma() * distrib.coord_sys() * distrib.scaling() * T; - Vector sol = distrib.mean() + distrib.sigma() - * distrib.coord_sys().dot( distrib.scaling().dot( T ) ); + Vector sol = distrib.mean() + distrib.sigma() * distrib.coord_sys() * (distrib.scaling().dot(T) ); + /*Vector sol = distrib.mean() + distrib.sigma() + * distrib.coord_sys().dot( distrib.scaling().dot( T ) );*/ // copy in the EOT structure (more probably a vector) EOT solution( size );