diff --git a/edo/application/cmaes/main.cpp b/edo/application/cmaes/main.cpp index 18c3f093f..4da2fdbed 100644 --- a/edo/application/cmaes/main.cpp +++ b/edo/application/cmaes/main.cpp @@ -41,8 +41,8 @@ Authors: #include "Sphere.h" -typedef eoReal EOT; -typedef edoNormalMulti< EOT > Distrib; +typedef eoReal RealVec; +typedef edoNormalAdaptive< RealVec > Distrib; int main(int ac, char** av) @@ -57,26 +57,34 @@ int main(int ac, char** av) eoState state; // Instantiate all needed parameters for EDA algorithm - double selection_rate = parser.createParam((double)0.5, "selection_rate", "Selection Rate", 'R', section).value(); // R - - eoSelect< EOT >* selector = new eoDetSelect< EOT >( selection_rate ); - state.storeFunctor(selector); - - edoEstimator< Distrib >* estimator = new edoEstimatorNormalMulti< EOT >(); - state.storeFunctor(estimator); - - eoEvalFunc< EOT >* plainEval = new Rosenbrock< EOT >(); - state.storeFunctor(plainEval); + //double selection_rate = parser.createParam((double)0.5, "selection_rate", "Selection Rate", 'R', section).value(); // R unsigned long max_eval = parser.getORcreateParam((unsigned long)0, "maxEval", "Maximum number of evaluations (0 = none)", 'E', "Stopping criterion").value(); // E - eoEvalFuncCounterBounder< EOT > eval(*plainEval, max_eval); + + unsigned int dim = parser.createParam((unsigned int)10, "dimension-size", "Dimension size", 'd', section).value(); // d + + + double mu = dim / 2; + + + edoNormalAdaptive distribution(dim); + + eoSelect< RealVec >* selector = new eoRankMuSelect< RealVec >( mu ); + state.storeFunctor(selector); + + edoEstimator< Distrib >* estimator = new edoEstimatorNormalAdaptive( distribution ); + state.storeFunctor(estimator); + + eoEvalFunc< RealVec >* plainEval = new Rosenbrock< RealVec >(); + state.storeFunctor(plainEval); + + eoEvalFuncCounterBounder< RealVec > eval(*plainEval, max_eval); eoRndGenerator< double >* gen = new eoUniformGenerator< double >(-5, 5); state.storeFunctor(gen); - unsigned int dimension_size = parser.createParam((unsigned int)10, "dimension-size", "Dimension size", 'd', section).value(); // d - eoInitFixedLength< EOT >* init = new eoInitFixedLength< EOT >( dimension_size, *gen ); + eoInitFixedLength< RealVec >* init = new eoInitFixedLength< RealVec >( dim, *gen ); state.storeFunctor(init); @@ -84,28 +92,28 @@ int main(int ac, char** av) // Generation of population from do_make_pop (creates parameters, manages persistance and so on...) // ... and creates the parameters: L P r S // this first sampler creates a uniform distribution independently from our distribution (it does not use edoUniform). - eoPop< EOT >& pop = do_make_pop(parser, state, *init); + eoPop< RealVec >& pop = do_make_pop(parser, state, *init); // (2) First evaluation before starting the research algorithm apply(eval, pop); // Prepare bounder class to set bounds of sampling. // This is used by edoSampler. - edoBounder< EOT >* bounder = - new edoBounderRng< EOT >( EOT(pop[0].size(), -5), EOT(pop[0].size(), 5), *gen); // FIXME do not use hard-coded bounds + edoBounder< RealVec >* bounder = + new edoBounderRng< RealVec >( RealVec(dim, -5), RealVec(dim, 5), *gen); // FIXME do not use hard-coded bounds state.storeFunctor(bounder); // Prepare sampler class with a specific distribution - edoSampler< Distrib >* sampler = new edoSamplerNormalMulti< EOT >( *bounder ); + edoSampler< Distrib >* sampler = new edoSamplerNormalAdaptive< RealVec >( *bounder ); state.storeFunctor(sampler); - + // stopping criteria // ... and creates the parameter letters: C E g G s T - eoContinue< EOT >& eo_continue = do_make_continue(parser, state, eval); - + eoContinue< RealVec >& eo_continue = do_make_continue(parser, state, eval); + // population output - eoCheckPoint< EOT >& pop_continue = do_make_checkpoint(parser, state, eval, eo_continue); - + eoCheckPoint< RealVec >& pop_continue = do_make_checkpoint(parser, state, eval, eo_continue); + // distribution output edoDummyContinue< Distrib >* dummy_continue = new edoDummyContinue< Distrib >(); state.storeFunctor(dummy_continue); @@ -115,9 +123,9 @@ int main(int ac, char** av) // eoEPRemplacement causes the using of the current and previous // sample for sampling. - eoReplacement< EOT >* replacor = new eoEPReplacement< EOT >(pop.size()); + eoReplacement< RealVec >* replacor = new eoEPReplacement< RealVec >(pop.size()); state.storeFunctor(replacor); - + // Some stuff to display helper when we are using -h option if (parser.userNeedsHelp()) { @@ -129,40 +137,11 @@ int main(int ac, char** av) make_verbose(parser); make_help(parser); - // population output (after helper) - // - // FIXME: theses objects are instanciated there in order to avoid a folder - // removing as edoFileSnapshot does within ctor. - edoPopStat< EOT >* popStat = new edoPopStat; - state.storeFunctor(popStat); - pop_continue.add(*popStat); - - edoFileSnapshot* fileSnapshot = new edoFileSnapshot("EDA_ResPop"); - state.storeFunctor(fileSnapshot); - fileSnapshot->add(*popStat); - pop_continue.add(*fileSnapshot); - - // distribution output (after helper) - edoDistribStat< Distrib >* distrib_stat = new edoStatNormalMulti< EOT >(); - state.storeFunctor(distrib_stat); - - distribution_continue->add( *distrib_stat ); - - // eoMonitor* stdout_monitor = new eoStdoutMonitor(); - // state.storeFunctor(stdout_monitor); - // stdout_monitor->add(*distrib_stat); - // distribution_continue->add( *stdout_monitor ); - - eoFileMonitor* file_monitor = new eoFileMonitor("eda_distribution_bounds.txt"); - state.storeFunctor(file_monitor); - file_monitor->add(*distrib_stat); - distribution_continue->add( *file_monitor ); - - eoPopLoopEval popEval( eval ); + eoPopLoopEval popEval( eval ); // EDA algorithm configuration - edoAlgo< Distrib >* algo = new edoEDA< Distrib > - (popEval, *selector, *estimator, *sampler, *replacor, + edoAlgo< Distrib >* algo = new edoAlgoAdaptive< Distrib > + (distribution, popEval, *selector, *estimator, *sampler, *replacor, pop_continue, *distribution_continue ); diff --git a/edo/build_gcc_linux_release b/edo/build_gcc_linux_release index 389d9b745..a1a0f53b9 100755 --- a/edo/build_gcc_linux_release +++ b/edo/build_gcc_linux_release @@ -4,5 +4,6 @@ mkdir -p release cd release cmake -DWITH_EIGEN=1 .. #cmake -DWITH_BOOST=1 .. + make cd .. diff --git a/edo/src/edoAlgoAdaptive.h b/edo/src/edoAlgoAdaptive.h new file mode 100644 index 000000000..dbb8046fc --- /dev/null +++ b/edo/src/edoAlgoAdaptive.h @@ -0,0 +1,204 @@ +/* +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 _edoAlgoAdaptive_h +#define _edoAlgoAdaptive_h + +#include + +#include + +#include "edoAlgo.h" +#include "edoEstimator.h" +#include "edoModifierMass.h" +#include "edoSampler.h" +#include "edoContinue.h" + +//! edoEDA< D > + +/** A generic stochastic search template for algorithms that need a distribution parameter. + */ +template < typename EOD > +class edoAlgoAdaptive : public edoAlgo< EOD > +{ +public: + //! Alias for the type EOT + typedef typename EOD::EOType EOType; + + //! Alias for the atom type + typedef typename EOType::AtomType AtomType; + + //! Alias for the fitness + typedef typename EOType::Fitness Fitness; + +public: + + /*! + Takes algo operators, all are mandatory + + \param distrib A distribution to use, if you want to update this parameter (e.gMA-ES) instead of replacing it (e.g. an EDA) + \param evaluation Evaluate a population + \param selector Selection of the best candidate solutions in the population + \param estimator Estimation of the distribution parameters + \param sampler Generate feasible solutions using the distribution + \param replacor Replace old solutions by new ones + \param pop_continuator Stopping criterion based on the population features + \param distribution_continuator Stopping criterion based on the distribution features + */ + edoAlgoAdaptive( + EOD & distrib, + eoPopEvalFunc < EOType > & evaluator, + eoSelect< EOType > & selector, + edoEstimator< EOD > & estimator, + edoSampler< EOD > & sampler, + eoReplacement< EOType > & replacor, + eoContinue< EOType > & pop_continuator, + edoContinue< EOD > & distribution_continuator + ) : + _distrib(distrib), + _evaluator(evaluator), + _selector(selector), + _estimator(estimator), + _sampler(sampler), + _replacor(replacor), + _pop_continuator(pop_continuator), + _dummy_continue(), + _distribution_continuator(distribution_continuator) + {} + + + //! constructor without an edoContinue + /*! + Takes algo operators, all are mandatory + + \param distrib A distribution to use, if you want to update this parameter (e.gMA-ES) instead of replacing it (e.g. an EDA) + \param evaluation Evaluate a population + \param selector Selection of the best candidate solutions in the population + \param estimator Estimation of the distribution parameters + \param sampler Generate feasible solutions using the distribution + \param replacor Replace old solutions by new ones + \param pop_continuator Stopping criterion based on the population features + */ + edoAlgoAdaptive ( + EOD & distrib, + eoPopEvalFunc < EOType > & evaluator, + eoSelect< EOType > & selector, + edoEstimator< EOD > & estimator, + edoSampler< EOD > & sampler, + eoReplacement< EOType > & replacor, + eoContinue< EOType > & pop_continuator + ) : + _distrib( distrib ), + _evaluator(evaluator), + _selector(selector), + _estimator(estimator), + _sampler(sampler), + _replacor(replacor), + _pop_continuator(pop_continuator), + _dummy_continue(), + _distribution_continuator( _dummy_continue ) + {} + + /** Call the algorithm + * + * \param pop the population of candidate solutions + * \return void + */ + void operator ()(eoPop< EOType > & pop) + { + assert(pop.size() > 0); + + eoPop< EOType > current_pop; + eoPop< EOType > selected_pop; + + // update the extern distribution passed to the estimator (cf. CMA-ES) + // OR replace the dummy distribution for estimators that do not need extern distributions (cf. EDA) + _distrib = _estimator(pop); + + // Evaluating a first time the candidate solutions + // The first pop is not supposed to be evaluated (@see eoPopLoopEval). + // _evaluator( current_pop, pop ); + + do { + // (1) Selection of the best points in the population + _selector(pop, selected_pop); + assert( selected_pop.size() > 0 ); + + // (2) Estimation of the distribution parameters + _distrib = _estimator(selected_pop); + + // (3) sampling + // The sampler produces feasible solutions (@see edoSampler that + // encapsulate an edoBounder) + current_pop.clear(); + for( unsigned int i = 0; i < pop.size(); ++i ) { + current_pop.push_back( _sampler(_distrib) ); + } + + // (4) Evaluate new solutions + _evaluator( pop, current_pop ); + + // (5) Replace old solutions by new ones + _replacor(pop, current_pop); // e.g. copy current_pop in pop + + } while( _distribution_continuator( _distrib ) && _pop_continuator( pop ) ); + } // operator() + + +protected: + + //! The distribution that you want to update + EOD & _distrib; + + //! A full evaluation function. + eoPopEvalFunc & _evaluator; + + //! A EOType selector + eoSelect & _selector; + + //! A EOType estimator. It is going to estimate distribution parameters. + edoEstimator & _estimator; + + //! A D sampler + edoSampler & _sampler; + + //! A EOType replacor + eoReplacement & _replacor; + + //! A EOType population continuator + eoContinue & _pop_continuator; + + //! A D continuator that always return true + edoDummyContinue _dummy_continue; + + //! A D continuator + edoContinue & _distribution_continuator; + +}; + +#endif // !_edoAlgoAdaptive_h + diff --git a/edo/src/edoAlgoStateless.h b/edo/src/edoAlgoStateless.h new file mode 100644 index 000000000..45bcea06f --- /dev/null +++ b/edo/src/edoAlgoStateless.h @@ -0,0 +1,109 @@ +/* +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 _edoAlgoStateless_h +#define _edoAlgoStateless_h + +#include "edoAlgoAdaptive.h" + +/** A generic stochastic search template for algorithms that need a distribution parameter but replace it rather than update it + * + * This use a default dummy distribution, for algorithms willing to replace it instead of updating + * Thus we can instanciate _distrib on this and replace it at the first iteration with an estimator. + * This is why an edoDistrib must have an empty constructor. + */ +template < typename EOD > +class edoAlgoStateless : public edoAlgoAdaptive< EOD > +{ +public: + //! Alias for the type EOT + typedef typename EOD::EOType EOType; + + //! Alias for the atom type + typedef typename EOType::AtomType AtomType; + + //! Alias for the fitness + typedef typename EOType::Fitness Fitness; + +public: + + /** Full constructor + \param evaluation Evaluate a population + \param selector Selection of the best candidate solutions in the population + \param estimator Estimation of the distribution parameters + \param sampler Generate feasible solutions using the distribution + \param replacor Replace old solutions by new ones + \param pop_continuator Stopping criterion based on the population features + \param distribution_continuator Stopping criterion based on the distribution features + + You are not supposed to override the tmp_distrib default initalization, or else use edoAlgoAdaptive + */ + edoAlgoStateless( + eoPopEvalFunc < EOType > & evaluator, + eoSelect< EOType > & selector, + edoEstimator< EOD > & estimator, + edoSampler< EOD > & sampler, + eoReplacement< EOType > & replacor, + eoContinue< EOType > & pop_continuator, + edoContinue< EOD > & distribution_continuator, + EOD* tmp_distrib = (new EOD()) + ) : + edoAlgoAdaptive( *tmp_distrib, evaluator, selector, estimator, sampler, replacor, pop_continuator, distribution_continuator) + {} + + /** Constructor without an edoContinue + + \param evaluation Evaluate a population + \param selector Selection of the best candidate solutions in the population + \param estimator Estimation of the distribution parameters + \param sampler Generate feasible solutions using the distribution + \param replacor Replace old solutions by new ones + \param pop_continuator Stopping criterion based on the population features + + You are not supposed to override the tmp_distrib default initalization, or else use edoAlgoAdaptive + */ + edoAlgoStateless ( + eoPopEvalFunc < EOType > & evaluator, + eoSelect< EOType > & selector, + edoEstimator< EOD > & estimator, + edoSampler< EOD > & sampler, + eoReplacement< EOType > & replacor, + eoContinue< EOType > & pop_continuator, + EOD* tmp_distrib = (new EOD()) + ) : + edoAlgoAdaptive( *tmp_distrib, evaluator, selector, estimator, sampler, replacor, pop_continuator) + {} + + ~edoAlgoStateless() + { + // delete the temporary distrib allocated in constructors + delete &(this->_distrib); + } +}; + +#endif // !_edoAlgoStateless_h + diff --git a/edo/src/edoEstimatorAdaptive.h b/edo/src/edoEstimatorAdaptive.h new file mode 100644 index 000000000..a07f47d1f --- /dev/null +++ b/edo/src/edoEstimatorAdaptive.h @@ -0,0 +1,55 @@ +/* +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 _edoEstimatorAdaptive_h +#define _edoEstimatorAdaptive_h + +#include +#include + +#include "edoEstimator.h" + +/** An interface that explicits the needs for a permanent distribution + * that will be updated by operators. + */ +template < typename EOD > +class edoEstimatorAdaptive : public edoEstimator +{ +public: + typedef typename EOD::EOType EOType; + + edoEstimatorAdaptive( EOD& distrib ) : _distrib(distrib) {} + + // virtual D operator() ( eoPop< EOT >& )=0 (provided by eoUF< A1, R >) + + EOD & distribution() const { return _distrib; } + +protected: + EOD & _distrib; +}; + +#endif // !_edoEstimatorAdaptive_h diff --git a/edo/src/edoEstimatorNormalAdaptive.h b/edo/src/edoEstimatorNormalAdaptive.h new file mode 100644 index 000000000..a2566772b --- /dev/null +++ b/edo/src/edoEstimatorNormalAdaptive.h @@ -0,0 +1,249 @@ +/* +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; // column vectors @see edoNormalAdaptive + typedef typename EOD::Matrix Matrix; + + edoEstimatorNormalAdaptive( EOD& distrib ) : + edoEstimatorAdaptive( distrib ), + _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 < pop_size; ++i ) { + double w_i = log( pop_size + 0.5 ) - log( i + 1 ); + weights(i) = w_i; + sum_w += w_i; + } + // normalization of weights + weights /= sum_w; + + assert( weights.size() == pop_size); + 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( N, lambda ); + + // copy the pop (most probably a vector of vectors) in a Eigen3 matrix + for( unsigned int d = 0; d < N; ++d ) { + for( unsigned int i = 0; i < lambda; ++i ) { + arx(d,i) = pop[i][d]; // NOTE: pop = arx.transpose() + } // dimensions + } // individuals + + // muXone array for weighted recombination + Eigen::VectorXd weights = edoCMAESweights( lambda ); + assert( weights.size() == lambda ); + + // 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(); + assert( invsqrtC.innerSize() == d.coord_sys().innerSize() ); + assert( invsqrtC.outerSize() == d.coord_sys().outerSize() ); + + // 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(); + assert( xold.size() == N ); + // xmean ( N, 1 ) = arx( N, lambda ) * weights( lambda, 1 ) + Vector xmean = arx * weights; + assert( xmean.size() == N ); + d.mean( xmean ); + + + /********************************************************************** + * 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 xmu( N, lambda); + xmu = xold.rowwise().replicate(lambda); + assert( xmu.innerSize() == N ); + assert( xmu.outerSize() == lambda ); + Matrix artmp = (1.0/d.sigma()) * (arx - xmu); + // Matrix artmp = (1.0/d.sigma()) * arx - xold.colwise().replicate(lambda); + assert( artmp.innerSize() == N && artmp.outerSize() == lambda ); + + + /********************************************************************** + * 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(); + assert( D.innerSize() == N && D.outerSize() == N ); + + // from variance to standard deviations + D.cwiseSqrt(); + d.scaling( D.diagonal() ); + } + + return d; + } // operator() + +protected: + + 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 new file mode 100644 index 000000000..7c6e80a6c --- /dev/null +++ b/edo/src/edoNormalAdaptive.h @@ -0,0 +1,121 @@ + +/* +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 Dreo + Pierre Savéant +*/ + +#ifndef _edoNormalAdaptive_h +#define _edoNormalAdaptive_h + +#include "edoDistrib.h" + +#ifdef WITH_EIGEN + +#include + +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; // column vectors ( n lines, 1 column) + typedef Eigen::Matrix< AtomType, Eigen::Dynamic, Eigen::Dynamic> Matrix; + + edoNormalAdaptive( unsigned int dim = 1 ) : + _dim(dim), + _mean( Vector::Zero(dim) ), + _C( Matrix::Identity(dim,dim) ), + _B( Matrix::Identity(dim,dim) ), + _D( Vector::Constant( dim, 1) ), + _sigma(1.0), + _p_c( Vector::Zero(dim) ), + _p_s( Vector::Zero(dim) ) + { + 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 _C;} + Matrix coord_sys() const {return _B;} + Vector scaling() const {return _D;} + double sigma() const {return _sigma;} + Vector path_covar() const {return _p_c;} + Vector path_sigma() const {return _p_s;} + + void mean( Vector m ) { _mean = m; assert( m.size() == _dim ); } + void covar( Matrix c ) { _C = c; assert( c.innerSize() == _dim && c.outerSize() == _dim ); } + void coord_sys( Matrix b ) { _B = b; assert( b.innerSize() == _dim && b.outerSize() == _dim ); } + void scaling( Vector d ) { _D = d; assert( d.size() == _dim ); } + void sigma( double s ) { _sigma = s; assert( s != 0.0 );} + void path_covar( Vector p ) { _p_c = p; assert( p.size() == _dim ); } + void path_sigma( Vector p ) { _p_s = p; assert( p.size() == _dim ); } + +private: + unsigned int _dim; + Vector _mean; // + Matrix _C; // covariance matrix + Matrix _B; // eigen vectors / coordinates system + Vector _D; // eigen values / scaling + double _sigma; // + Vector _p_c; // evolution path for C + Vector _p_s; // evolution path for sigma +}; + +#endif // WITH_EIGEN + +#endif // !_edoNormalAdaptive_h diff --git a/edo/src/edoSamplerNormalAdaptive.h b/edo/src/edoSamplerNormalAdaptive.h new file mode 100644 index 000000000..e8188105a --- /dev/null +++ b/edo/src/edoSamplerNormalAdaptive.h @@ -0,0 +1,91 @@ +/* +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 _edoSamplerNormalAdaptive_h +#define _edoSamplerNormalAdaptive_h + +#include +#include + +#include + +/** Sample points in a multi-normal law defined by a mean vector and a covariance matrix. + * + * Given M the mean vector and V the covariance matrix, of order n: + * - draw a vector T in N(0,I) (i.e. each value is drawn in a normal law with mean=0 an stddev=1) + * - compute the Cholesky decomposition L of V (i.e. such as V=LL*) + * - return X = M + LT + */ +#ifdef WITH_EIGEN + +template< class EOT, typename EOD = edoNormalAdaptive< EOT > > +class edoSamplerNormalAdaptive : public edoSampler< EOD > +{ +public: + typedef typename EOT::AtomType AtomType; + + typedef typename EOD::Vector Vector; + typedef typename EOD::Matrix Matrix; + + edoSamplerNormalAdaptive( edoRepairer & repairer ) + : edoSampler< EOD >( repairer) + {} + + + EOT sample( EOD& distrib ) + { + unsigned int N = distrib.size(); + assert( N > 0); + + // T = vector of size elements drawn in N(0,1) + Vector T( N ); + for ( unsigned int i = 0; i < N; ++i ) { + T( i ) = rng.normal(); + } + assert(T.innerSize() == N ); + assert(T.outerSize() == 1); + + // mean(N,1) + sigma * B(N,N) * ( D(N,1) .* T(N,1) ) + Vector sol = distrib.mean() + + distrib.sigma() + * distrib.coord_sys() * (distrib.scaling().cwiseProduct(T) ); // C * T = B * (D .* T) + assert( sol.size() == N ); + /*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( N ); + for( unsigned int i = 0; i < N; i++ ) { + solution[i]= sol(i); + } + + return solution; + } +}; +#endif // WITH_EIGEN + +#endif // !_edoSamplerNormalAdaptive_h diff --git a/eo/src/eoRankMuSelect.h b/eo/src/eoRankMuSelect.h new file mode 100644 index 000000000..6e0a737ea --- /dev/null +++ b/eo/src/eoRankMuSelect.h @@ -0,0 +1,54 @@ +/* +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) 2012 Thales group +*/ +/* +Authors: + Johann Dréo +*/ + + +#ifndef _eoRankMuSelect_h +#define _eoRankMuSelect_h + +#include "eoDetSelect.h" + +/** Selects the "Mu" bests individuals. + * + * Note: sorts the population before trucating it. + * + * @ingroup Selectors +*/ +template +class eoRankMuSelect : public eoDetSelect +{ +public : + // false, because mu is not a rate + eoRankMuSelect( unsigned int mu ) : eoDetSelect( mu, false ) {} + + void operator()(const eoPop& source, eoPop& dest) + { + eoPop tmp( source ); + tmp.sort(); + eoDetSelect::operator()( tmp, dest ); + } +}; + +#endif // !_eoRankMuselect_h