Merge branch 'master' of https://github.com/nojhan/eodev

edo/CMakeLists.txt

@@ -125,7 +125,7 @@ SET(SAMPLE_SRCS)
 
 ADD_SUBDIRECTORY(src)
 ADD_SUBDIRECTORY(application)
-#ADD_SUBDIRECTORY(test)
+ADD_SUBDIRECTORY(test)
 ADD_SUBDIRECTORY(doc)
 
 ######################################################################################

edo/src/edoEstimatorNormalMono.h

@@ -39,64 +39,84 @@ Authors:
 template < typename EOT >
 class edoEstimatorNormalMono : public edoEstimator< edoNormalMono< EOT > >
 {
-public:
-    typedef typename EOT::AtomType AtomType;
-
-    class Variance
-    {
     public:
-        Variance() : _sumvar(0){}
+        typedef typename EOT::AtomType AtomType;
 
-        void update(AtomType v)
+        //! Knuth's algorithm, online variance, numericably stable
+        class Variance
         {
-            _n++;
+            public:
+                Variance() : _n(0), _mean(0), _M2(0) {}
 
-            AtomType d = v - _mean;
-
-            _mean += 1 / _n * d;
-            _sumvar += (_n - 1) / _n * d * d;
-        }
-
-        AtomType get_mean() const {return _mean;}
-        AtomType get_var() const {return _sumvar / (_n - 1);}
-        AtomType get_std() const {return sqrt( get_var() );}
-
-    private:
-        AtomType _n;
-        AtomType _mean;
-        AtomType _sumvar;
-    };
-
-public:
-    edoNormalMono< EOT > operator()(eoPop<EOT>& pop)
-    {
-        unsigned int popsize = pop.size();
-        assert(popsize > 0);
-
-        unsigned int dimsize = pop[0].size();
-        assert(dimsize > 0);
-
-        std::vector< Variance > var( dimsize );
-
-        for (unsigned int i = 0; i < popsize; ++i)
-            {
-            for (unsigned int d = 0; d < dimsize; ++d)
+                void update(AtomType x)
                 {
-                var[d].update( pop[i][d] );
+                    _n++;
+
+                    AtomType delta = x - _mean;
+
+                    _mean += delta / _n;
+                    _M2 += delta * ( x - _mean );
+                }
+
+                AtomType mean() const {return _mean;}
+
+                //! Population variance
+                AtomType var_n() const {
+                    assert( _n > 0 );
+                    return _M2 / _n;
+                }
+
+                /** Sample variance (using Bessel's correction)
+                 * is an unbiased estimate of the population variance,
+                 * but it has uniformly higher mean squared error
+                 */
+                AtomType var() const {
+                    assert( _n > 1 );
+                    return _M2 / (_n - 1);
+                }
+
+                //! Population standard deviation
+                AtomType std_n() const {return sqrt( var_n() );}
+
+                //! Sample standard deviation, is a biased estimate of the population standard deviation
+                AtomType std() const {return sqrt( var() );}
+
+            private:
+                AtomType _n;
+                AtomType _mean;
+                AtomType _M2;
+        };
+
+    public:
+        edoNormalMono< EOT > operator()(eoPop<EOT>& pop)
+        {
+            unsigned int popsize = pop.size();
+            assert(popsize > 0);
+
+            unsigned int dimsize = pop[0].size();
+            assert(dimsize > 0);
+
+            std::vector< Variance > var( dimsize );
+
+            for (unsigned int i = 0; i < popsize; ++i)
+            {
+                for (unsigned int d = 0; d < dimsize; ++d)
+                {
+                    var[d].update( pop[i][d] );
                 }
             }
 
-        EOT mean( dimsize );
-        EOT variance( dimsize );
+            EOT mean( dimsize );
+            EOT variance( dimsize );
 
-        for (unsigned int d = 0; d < dimsize; ++d)
+            for (unsigned int d = 0; d < dimsize; ++d)
             {
-            mean[d] = var[d].get_mean();
-            variance[d] = var[d].get_var();
+                mean[d] = var[d].mean();
+                variance[d] = var[d].var_n();
             }
 
-        return edoNormalMono< EOT >( mean, variance );
-    }
+            return edoNormalMono< EOT >( mean, variance );
+        }
 };
 
 #endif // !_edoEstimatorNormalMono_h

edo/test/CMakeLists.txt

@@ -34,6 +34,7 @@ INCLUDE_DIRECTORIES(${CMAKE_SOURCE_DIR}/application/common)
 
 SET(SOURCES
     #t-cholesky
+  t-variance
   t-edoEstimatorNormalMulti
   t-mean-distance
   t-bounderno

edo/test/t-variance.cpp

@@ -0,0 +1,38 @@
+
+#include <iostream>
+#include <vector>
+
+#include <eo>
+#include <es.h>
+#include <edo>
+
+int main()
+{
+    typedef eoReal<eoMinimizingFitness> Vec;
+
+    eoPop<Vec> pop;
+    for( unsigned int i=1; i<7; ++i) {
+        Vec indiv(1,i);
+        pop.push_back( indiv );
+        std::clog << indiv << " ";
+    }
+    std::clog << std::endl;
+
+    edoEstimatorNormalMono<Vec> estimator;
+
+    edoNormalMono<Vec> distrib = estimator(pop);
+
+    Vec ex_mean(1,3.5);
+    Vec ex_var(1,17.5/6);
+    Vec es_mean = distrib.mean();
+    Vec es_var = distrib.variance();
+
+    std::cout << "expected  mean=" << ex_mean << " variance=" << ex_var << std::endl;
+    std::cout << "estimated mean=" << es_mean << " variance=" << es_var << std::endl;
+
+    for( unsigned int i=0; i<ex_mean.size(); ++i ) {
+        assert( es_mean[i] == ex_mean[i] );
+        assert(  es_var[i] == ex_var[i] );
+    }
+}
+

eo/src/eo

@@ -78,6 +78,7 @@
 #include <eoEvalCounterThrowException.h>
 #include <eoEvalTimeThrowException.h>
 #include <eoEvalUserTimeThrowException.h>
+#include <eoEvalKeepBest.h>
 
 // Continuators - all include eoContinue.h
 #include <eoCombinedContinue.h>

eo/src/eoEvalKeepBest.h

@@ -0,0 +1,107 @@
+
+/*
+   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;
+   version 2 of the License.
+
+   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
+
+   © 2012 Thales group
+
+    Authors:
+        Johann Dreo <johann.dreo@thalesgroup.com>
+*/
+
+#ifndef eoEvalKeepBest_H
+#define eoEvalKeepBest_H
+
+#include <eoEvalFunc.h>
+#include <utils/eoParam.h>
+
+/**
+  Evaluate with the given evaluator and keep the best individual found so far.
+
+  This is useful if you use a non-monotonic algorithm, such as CMA-ES, where the
+  population's best fitness can decrease between two generations. This is
+  sometime necessary and one can't use elitist replacors, as one do not want to
+  introduce a bias in the population.
+
+  The eoEvalBestKeep is a wrapper around a classical evaluator, that keep the
+  best individual it has found since its instanciation.
+
+  To get the best individual, you have to call best_element() on the
+  eoEvalKeepBest itself, and not on the population (or else you would get the
+  best individual found at the last generation).
+
+  Example:
+
+    MyEval true_eval;
+    eoEvalKeepBest<T> wrapped_eval( true_eval );
+
+    // as an interesting side effect, you will get the best individual since
+    // initalization.
+    eoPop<T> pop( my_init );
+    eoPopLoopEval<T> loop_eval( wrapped_eval );
+
+    loop_eval( pop );
+
+    eoEasyEA algo( …, wrapped_eval, … );
+
+    algo(pop);
+
+    // do not use pop.best_element()!
+    std::cout << wrapped_eval.best_element() << std::endl;
+
+  @ingroup Evaluation
+  */
+template<class EOT> class eoEvalKeepBest : public eoEvalFunc<EOT>, public eoValueParam<EOT>
+{
+    public :
+        eoEvalKeepBest(eoEvalFunc<EOT>& _func, std::string _name = "VeryBest. ")
+            : eoValueParam<EOT>(EOT(), _name), func(_func) {}
+
+        virtual void operator()(EOT& sol)
+        {
+            if( sol.invalid() ) {
+                func(sol); // evaluate
+
+                // if there is no best kept
+                if( this->value().invalid() ) {
+                    // take the first individual as best
+                    this->value() = sol;
+                } else {
+                    // if sol is better than the kept individual
+                    if( sol.fitness() > this->value().fitness() ) {
+                        this->value() = sol;
+                    }
+                }
+            } // if invalid
+        }
+
+        //! Return the best individual found so far.
+        EOT best_element()
+        {
+            return this->value();
+        }
+
+        /** Reset the best individual to the given one. If no individual is
+         * provided, the next evaluated one will be taken as a reference.
+         */
+        void reset( const EOT& new_best = EOT() )
+        {
+            this->value() = new_best;
+        }
+
+    protected :
+        eoEvalFunc<EOT>& func;
+};
+
+#endif

eo/test/CMakeLists.txt

@@ -34,6 +34,7 @@ ENDIF()
 ######################################################################################
 
 SET (TEST_LIST
+  t-eoEvalKeepBest
   t-eoInitVariableLength
   t-eofitness
   t-eoRandom

eo/test/t-eoEvalKeepBest.cpp

@@ -0,0 +1,78 @@
+#include <iostream>
+
+#include <es/make_real.h>
+#include <apply.h>
+#include <eoEvalKeepBest.h>
+#include "real_value.h"
+
+using namespace std;
+
+int main(int argc, char* argv[])
+{
+
+    typedef eoReal<eoMinimizingFitness> EOT;
+
+    eoParser parser(argc, argv);  // for user-parameter reading
+    eoState state;    // keeps all things allocated
+
+
+    /*********************************************
+     * problem or representation dependent stuff *
+     *********************************************/
+
+    // The evaluation fn - encapsulated into an eval counter for output
+    eoEvalFuncPtr<EOT, double, const std::vector<double>&>
+        main_eval( real_value ); // use a function defined in real_value.h
+
+    // wrap the evaluation function in a call counter
+    eoEvalFuncCounter<EOT> eval_counter(main_eval);
+
+    // the genotype - through a genotype initializer
+    eoRealInitBounded<EOT>& init = make_genotype(parser, state, EOT());
+
+    // Build the variation operator (any seq/prop construct)
+    eoGenOp<EOT>& op = make_op(parser, state, init);
+
+
+    /*********************************************
+     * Now the representation-independent things *
+     *********************************************/
+
+
+    // initialize the population - and evaluate
+    // yes, this is representation indepedent once you have an eoInit
+    eoPop<EOT>& pop   = make_pop(parser, state, init);
+
+    // stopping criteria
+    eoContinue<EOT> & term = make_continue(parser, state, eval_counter);
+
+    // things that are called at each generation
+    eoCheckPoint<EOT> & checkpoint = make_checkpoint(parser, state, eval_counter, term);
+
+    // wrap the evaluator in another one that will keep the best individual
+    // evaluated so far
+    eoEvalKeepBest<EOT> eval_keep( eval_counter );
+
+    // algorithm
+    eoAlgo<EOT>& ea = make_algo_scalar(parser, state, eval_keep, checkpoint, op);
+
+
+    /***************************************
+     * Now, call functors and DO something *
+     ***************************************/
+
+    // to be called AFTER all parameters have been read!
+    make_help(parser);
+
+    // evaluate intial population AFTER help and status in case it takes time
+    apply<EOT>(eval_keep, pop);
+
+    std::clog << "Best individual after initialization and " << eval_counter.value() << " evaluations" << std::endl;
+    std::cout << eval_keep.best_element() << std::endl;
+
+    ea(pop); // run the ea
+
+    std::cout << "Best individual after search and " << eval_counter.value() << " evaluations" << std::endl;
+    // you can also call value(), because it is an eoParam
+    std::cout << eval_keep.value() << std::endl;
+}