add FastGA foundry and eoStandardBitMutation variants

eo/src/eo

@@ -162,6 +162,7 @@
 #include "eoForge.h"
 #include "eoAlgoFoundry.h"
 #include "eoAlgoFoundryEA.h"
+#include "eoAlgoFoundryFastGA.h"
 #include "eoEvalFoundryEA.h"
 
 //-----------------------------------------------------------------------------

eo/src/eoAlgoFoundry.h

@@ -32,7 +32,7 @@
  * and hold the link to the encoding. @see eoAlgoFoundryEA
  *
  * As with eoForgeVector, adding a managed operator
- *is done through public member variable's `add` method,
+ * is done through public member variable's `add` method,
  * which takes the class name as template and its constructor's parameters
  * as arguments. For example:
  * @code
@@ -74,7 +74,7 @@ class eoOperatorFoundry : public eoForgeVector< Itf >
  * @endcode
  *
  * In a second step, the operators to be used should be selected
- * by indicating their index, just like if the foundry was a array:
+ * by indicating their index, just like if the foundry was an array:
  * @code
  * foundry = {0, 1, 2};
  * //         ^  ^  ^

eo/src/eoAlgoFoundryFastGA.h

@@ -0,0 +1,251 @@
+
+/*
+   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
+
+   © 2020 Thales group
+
+    Authors:
+        Johann Dreo <johann.dreo@thalesgroup.com>
+*/
+
+#ifndef _eoAlgoFoundryFastGA_H_
+#define _eoAlgoFoundryFastGA_H_
+
+#include <array>
+#include <tuple>
+#include <limits>
+
+/** A class that assemble an eoEasyEA on the fly, given a combination of available operators.
+ *
+ * The foundry should first be set up with sets of operators
+ * for the main modules of an EA:
+ * continuators, crossovers, mutations, selection and replacement operators.
+ *
+ * This is done through public member variable's `add` method,
+ * which takes the class name as template and its constructor's parameters
+ * as arguments. For example:
+ * @code
+ * foundry.selectors.add< eoStochTournamentSelect<EOT> >( 0.5 );
+ * @endcode
+ *
+ * @warning If the constructor takes a reference YOU SHOULD ABSOLUTELY wrap it
+ * in a `std::ref`, or it will silently be passed as a copy,
+ * which would effectively disable any link between operators.
+ *
+ * In a second step, the operators to be used should be selected
+ * by indicating their index, passing an array of eight elements:
+ * @code
+ * foundry.select({0, 1, 2, 3, 4, 5, 6, 7});
+ * @endcode
+ *
+ * @note: by default, the firsts of the eight operators are selected.
+ *
+ * If you don't (want to) recall the order of the operators in the encoding,
+ * you can use the `index()` member, for example:
+ * @code
+ * foundry.at(foundry.continuators.index()) = 2; // select the third continuator
+ * @endcode
+ *
+ * Now, you can call the fourdry just like any eoAlgo, by passing it an eoPop:
+ * @code
+ * foundry(pop);
+ * @encode
+ * It will instantiate the needed operators (only) and the algorithm itself on-the-fly,
+ * and then run it.
+ *
+ * @note: Thanks to the underlying eoOperatorFoundry, not all the added operators are instantiated.
+ * Every instantiation is deferred upon actual use. That way, you can still reconfigure them
+ * at any time with `eoForgeOperator::setup`, for example:
+ * @code
+ * foundry.selector.at(0).setup(0.5); // using constructor's arguments
+ * @endcode
+ *
+ * @ingroup Foundry
+ * @ingroup Algorithms
+ */
+template<class EOT>
+class eoAlgoFoundryFastGA : public eoAlgoFoundry<EOT>
+{
+    public:
+        /** The constructon only take an eval, because all other operators
+         * are stored in the public containers.
+         */
+        eoAlgoFoundryFastGA( eoInit<EOT> & init, eoEvalFunc<EOT>& eval, size_t max_evals = 10000, size_t max_restarts = std::numeric_limits<size_t>::max() ) :
+            eoAlgoFoundry<EOT>(8),
+            continuators(0, true), // Always re-instantiate continuators, because they hold a state.
+            crossover_rates(1, false),
+            crossovers(2, false),
+            mutation_rates(3, false),
+            mutations(4, false),
+            selectors(5, false),
+            pop_sizes(6, false),
+            replacements(7, false),
+            _eval(eval),
+            _init(init),
+            _max_evals(max_evals),
+            _max_restarts(max_restarts)
+        { }
+
+    public:
+
+        /* Operators containers @{ */
+        eoOperatorFoundry< eoContinue<EOT>    > continuators;
+        eoOperatorFoundry< double             > crossover_rates;
+        eoOperatorFoundry< eoQuadOp<EOT>      > crossovers;
+        eoOperatorFoundry< double             > mutation_rates;
+        eoOperatorFoundry< eoMonOp<EOT>       > mutations;
+        eoOperatorFoundry< eoSelectOne<EOT>   > selectors;
+        eoOperatorFoundry< size_t             > pop_sizes;
+        eoOperatorFoundry< eoReplacement<EOT> > replacements;
+        /* @} */
+
+        /** instantiate and call the pre-selected algorithm.
+         */
+        void operator()(eoPop<EOT>& pop)
+        {
+            assert(continuators.size() > 0); assert(this->at(continuators.index()) < continuators.size());
+            assert(  crossover_rates.size() > 0); assert(this->at(  crossover_rates.index()) <   crossover_rates.size());
+            assert(  crossovers.size() > 0); assert(this->at(  crossovers.index()) <   crossovers.size());
+            assert(   mutation_rates.size() > 0); assert(this->at(   mutation_rates.index()) <    mutation_rates.size());
+            assert(   mutations.size() > 0); assert(this->at(   mutations.index()) <    mutations.size());
+            assert(   selectors.size() > 0); assert(this->at(   selectors.index()) <    selectors.size());
+            assert(   pop_sizes.size() > 0); assert(this->at(   pop_sizes.index()) <    pop_sizes.size());
+            assert(replacements.size() > 0); assert(this->at(replacements.index()) < replacements.size());
+
+            // Crossover or clone
+            double cross_rate = this->crossover_rate();
+            eoProportionalOp<EOT> cross;
+            // Cross-over that produce only one offspring,
+            // made by wrapping the quad op (which produce 2 offsprings)
+            // in a bin op (which ignore the second offspring).
+            eoQuad2BinOp<EOT> single_cross(this->crossover());
+            cross.add(single_cross, cross_rate);
+            eoBinCloneOp<EOT> cross_clone;
+            cross.add(cross_clone, 1 - cross_rate); // Clone
+
+            // Mutation or clone
+            double mut_rate = this->mutation_rate();
+            eoProportionalOp<EOT> mut;
+            mut.add(this->mutation(), mut_rate);
+            eoMonCloneOp<EOT> mut_clone;
+            mut.add(mut_clone, 1 - mut_rate); // FIXME TBC
+
+            // Apply mutation after cross-over.
+            eoSequentialOp<EOT> variator;
+            variator.add(cross,1.0);
+            variator.add(mut,1.0);
+
+            // All variatiors
+            double lambda = this->pop_size();
+            eoGeneralBreeder<EOT> breeder(this->selector(), variator, lambda, /*as rate*/false);
+
+            // Objective function calls counter
+            eoEvalCounterThrowException<EOT> eval(_eval, _max_evals);
+            eoPopLoopEval<EOT> pop_eval(eval);
+
+            // Algorithm itself
+            eoEasyEA<EOT> algo = eoEasyEA<EOT>(this->continuator(), pop_eval, breeder, this->replacement());
+
+            // Restart wrapper
+            eoAlgoPopReset<EOT> reset_pop(_init, pop_eval);
+            eoGenContinue<EOT> restart_cont(_max_restarts);
+            eoAlgoRestart<EOT> restart(eval, algo, restart_cont, reset_pop);
+
+            try {
+                restart(pop);
+            } catch(eoMaxEvalException e) {
+                // In case some solutions were not evaluated when max eval occured.
+                eoPopLoopEval<EOT> pop_last_eval(_eval);
+                pop_last_eval(pop,pop);
+            }
+        }
+
+        /** Return an approximate name of the selected algorithm.
+         *
+         * @note: does not take into account parameters of the operators,
+         * only show class names.
+         */
+        std::string name()
+        {
+            std::ostringstream name;
+            name << this->at(continuators.index()) << " (" << this->continuator().className() << ") + ";
+            name << this->at(crossover_rates.index())   << " (" << this->crossover_rate().className()   << ") + ";
+            name << this->at(crossovers.index())   << " (" << this->crossover().className()   << ") + ";
+            name << this->at(mutation_rates.index())    << " (" << this->mutation_rate().className()    << ") + ";
+            name << this->at(mutations.index())    << " (" << this->mutation().className()    << ") + ";
+            name << this->at(selectors.index())    << " (" << this->selector().className()    << ") + ";
+            name << this->at(pop_sizes.index()) << " (" << this->pop_size().className() << ")";
+            name << this->at(replacements.index()) << " (" << this->replacement().className() << ")";
+            return name.str();
+        }
+
+    protected:
+        eoEvalFunc<EOT>& _eval;
+        eoInit<EOT>& _init;
+        const size_t _max_evals;
+        const size_t _max_restarts;
+
+    public:
+        eoContinue<EOT>& continuator()
+        {
+            assert(this->at(continuators.index()) < continuators.size());
+            return continuators.instantiate(this->at(continuators.index()));
+        }
+
+        double& crossover_rate()
+        {
+            assert(this->at(crossover_rates.index()) < crossover_rates.size());
+            return crossover_rates.instantiate(this->at(crossover_rates.index()));
+        }
+
+        eoQuadOp<EOT>& crossover()
+        {
+            assert(this->at(crossovers.index()) < crossovers.size());
+            return crossovers.instantiate(this->at(crossovers.index()));
+        }
+
+        double& mutation_rate()
+        {
+            assert(this->at(mutation_rates.index()) < mutation_rates.size());
+            return mutation_rates.instantiate(this->at(mutation_rates.index()));
+        }
+
+        eoMonOp<EOT>& mutation()
+        {
+            assert(this->at(mutations.index()) < mutations.size());
+            return mutations.instantiate(this->at(mutations.index()));
+        }
+
+        eoSelectOne<EOT>& selector()
+        {
+            assert(this->at(selectors.index()) < selectors.size());
+            return selectors.instantiate(this->at(selectors.index()));
+        }
+
+        size_t& pop_size()
+        {
+            assert(this->at(pop_sizes.index()) < pop_sizes.size());
+            return pop_sizes.instantiate(this->at(pop_sizes.index()));
+        }
+
+        eoReplacement<EOT>& replacement()
+        {
+            assert(this->at(replacements.index()) < replacements.size());
+            return replacements.instantiate(this->at(replacements.index()));
+        }
+
+};
+
+#endif // _eoAlgoFoundryFastGA_H_

eo/src/eoEasyEA.h

@@ -272,6 +272,8 @@ template<class EOT> class eoEasyEA: public eoAlgo<EOT>
 
               replace(_pop, offspring); // after replace, the new pop. is in _pop
 
+              std::cout << _pop << std::endl;
+
               if (pSize > _pop.size())
                 throw eoException("Population shrinking!");
               else if (pSize < _pop.size())

eo/src/eoMergeReduce.h

@@ -97,7 +97,8 @@ class eoCommaReplacement : public eoMergeReduce<EOT>
         virtual void operator()(eoPop<EOT>& _parents, eoPop<EOT>& _offspring)
         {
             // There must be more offsprings than parents, or else an exception will be raised
-            assert( _offspring.size() >= _parents.size() );
+            // Removed this assertion, which do not hold in some special case of singe indivudal algorithms.
+            //assert( _offspring.size() >= _parents.size() );
 
             eoMergeReduce<EOT>::operator()( _parents, _offspring );
         }

eo/src/eoOpContainer.h

@@ -106,6 +106,7 @@ public:
         for (size_t i = 0; i < rates.size(); ++i) {
             _pop.seekp(pos);
             do {
+            // std::clog << "Before:" << _pop.offspring().size() << " offsprings" << std::endl;
                 if (eo::rng.flip(rates[i])) {
                     //            try
                     //            {
@@ -126,6 +127,11 @@ public:
 
                 if (!_pop.exhausted())
                     ++_pop;
+
+            // std::clog << "After:" << _pop.offspring().size() << " offsprings" << std::endl;
+            // std::clog << _pop.offspring() << std::endl;
+            // std::clog << std::endl;
+
             }
             while (!_pop.exhausted());
         }
@@ -155,6 +161,7 @@ public:
 
       try
       {
+          // std::clog << "\t" << ops[i]->className() << std::endl;
         (*ops[i])(_pop);
         ++_pop;
       }

eo/src/ga.h

@@ -34,6 +34,7 @@
 
 // the operators
 #include "ga/eoBitOp.h"
+#include "ga/eoStandardBitMutation.h"
 
 // #include <ga/eoBitOpFactory.h> to be corrected - thanks someone!
 

eo/src/ga/eoStandardBitMutation.h

@@ -0,0 +1,213 @@
+
+#ifndef _eoStandardBitMutation_h_
+#define _eoStandardBitMutation_h_
+
+#include "../utils/eoRNG.h"
+
+/** Standard bit mutation with mutation rate p:
+ * choose k from the binomial distribution Bin(n,p) and apply flip_k(x).
+ *
+ * @ingroup Bitstrings
+ * @ingroup Variators
+ */
+template<class EOT>
+class eoStandardBitMutation : public eoMonOp<EOT>
+{
+    public:
+        eoStandardBitMutation(double rate = 0.5) :
+            _rate(rate),
+            _nb(1),
+            _bitflip(_nb)
+        {}
+
+        virtual bool operator()(EOT& chrom)
+        {
+            _nb = eo::rng.binomial(chrom.size(),_rate);
+            // BitFlip operator is bound to the _nb reference,
+            // thus one don't need to re-instantiate.
+            return _bitflip(chrom);
+        }
+
+    protected:
+        double _rate;
+        unsigned _nb;
+        eoDetSingleBitFlip<EOT> _bitflip;
+};
+
+/** Uniform bit mutation with mutation rate p:
+ * choose k from the uniform distribution U(0,n) and apply flip_k(x).
+ *
+ * @ingroup Bitstrings
+ * @ingroup Variators
+ */
+template<class EOT>
+class eoUniformBitMutation : public eoMonOp<EOT>
+{
+    public:
+        eoUniformBitMutation(double rate = 0.5) :
+            _rate(rate),
+            _nb(1),
+            _bitflip(_nb)
+        {}
+
+        virtual bool operator()(EOT& chrom)
+        {
+            _nb = eo::rng.random(chrom.size());
+            // BitFlip operator is bound to the _nb reference,
+            // thus one don't need to re-instantiate.
+            return _bitflip(chrom);
+        }
+
+    protected:
+        double _rate;
+        unsigned _nb;
+        eoDetSingleBitFlip<EOT> _bitflip;
+};
+
+
+/** Conditional standard bit mutation with mutation rate p:
+ * choose k from the binomial distribution Bin(n,p) until k >0
+ * and apply flip_k(x).
+ *
+ * This is identical to sampling k from the conditional binomial
+ * distribution Bin>0(n,p) which re-assigns the probability to sample
+ * a 0 proportionally to all values i ∈ [1..n]. 
+ *
+ * @ingroup Bitstrings
+ * @ingroup Variators
+ */
+template<class EOT>
+class eoConditionalBitMutation : public eoStandardBitMutation<EOT>
+{
+    public:
+        eoConditionalBitMutation(double rate = 0.5) :
+            eoStandardBitMutation<EOT>(rate)
+        {}
+
+        virtual bool operator()(EOT& chrom)
+        {
+            assert(chrom.size()>0);
+            this->_nb = eo::rng.binomial(chrom.size()-1,this->_rate);
+            this->_nb++;
+            // BitFlip operator is bound to the _nb reference,
+            // thus one don't need to re-instantiate.
+            return this->_bitflip(chrom);
+        }
+};
+
+/** Shifted standard bit mutation with mutation rate p:
+ * choose k from the binomial distribution Bin(n,p).
+ * When k= 0, set k= 1. Apply flip_k(x).
+ *
+ * This is identical to sampling k from the conditional binomial
+ * distribution Bin0→1(n,p) which re-assigns the probability to
+ * sample a 0 to sampling k= 1.
+ *
+ * @ingroup Bitstrings
+ * @ingroup Variators
+ */
+template<class EOT>
+class eoShiftedBitMutation : public eoStandardBitMutation<EOT>
+{
+    public:
+        eoShiftedBitMutation(double rate = 0.5) :
+            eoStandardBitMutation<EOT>(rate)
+        {}
+
+        virtual bool operator()(EOT& chrom)
+        {
+            assert(chrom.size()>0);
+            this->_nb = eo::rng.binomial(chrom.size()-1,this->_rate);
+            if(this->_nb == 0) {
+                this->_nb = 1;
+            }
+            // BitFlip operator is bound to the _nb reference,
+            // thus one don't need to re-instantiate.
+            return this->_bitflip(chrom);
+        }
+};
+
+/** Mutation which size is sample in a gaussian.
+ *
+ * sample k from the normal distribution N(pn,σ^2)
+ * and apply flip_k(x).
+ *
+ * From:
+ * Furong Ye, Carola Doerr, and Thomas Back.
+ * Interpolating local and global search by controllingthe variance of standard bit mutation.
+ * In 2019 IEEE Congress on Evolutionary Computation(CEC), pages 2292–2299.
+ *
+ * In contrast to standard bit mutation, this operators allows to scale
+ * the variance of the mutation strength independently of the mean.
+ *
+ * @ingroup Bitstrings
+ * @ingroup Variators
+ */
+template<class EOT>
+class eoNormalBitMutation : public eoStandardBitMutation<EOT>
+{
+    public:
+        eoNormalBitMutation(double rate = 0.5, double variance = 1) :
+            eoStandardBitMutation<EOT>(rate),
+            _variance(variance)
+        {}
+
+        virtual bool operator()(EOT& chrom)
+        {
+            this->_nb = eo::rng.normal(this->_rate * chrom.size(), _variance);
+            if(this->_nb >= chrom.size()) {
+                this->_nb = eo::rng.random(chrom.size());
+            }
+            // BitFlip operator is bound to the _nb reference,
+            // thus one don't need to re-instantiate.
+            return this->_bitflip(chrom);
+        }
+
+    protected:
+        double _variance;
+};
+
+/** Fast mutation which size is sampled from an adaptive power law.
+ *
+ * From:
+ * Benjamin Doerr, Huu Phuoc Le, Régis Makhmara, and Ta Duy Nguyen.
+ * Fast genetic algorithms.
+ * In Proc. of Genetic and Evolutionary Computation Conference (GECCO’17), pages 777–784.ACM, 2017.
+ *
+ * @ingroup Bitstrings
+ * @ingroup Variators
+ */
+template<class EOT>
+class eoFastBitMutation : public eoStandardBitMutation<EOT>
+{
+    public:
+        eoFastBitMutation(double rate = 0.5, double beta = 1.5) :
+            eoStandardBitMutation<EOT>(rate),
+            _beta(beta)
+        {
+            assert(beta > 1);
+        }
+
+        virtual bool operator()(EOT& chrom)
+        {
+            this->_nb = powerlaw(chrom.size(),_beta);
+            // BitFlip operator is bound to the _nb reference,
+            // thus one don't need to re-instantiate.
+            return this->_bitflip(chrom);
+        }
+
+    protected:
+
+        double powerlaw(unsigned n, double beta)
+        {
+            double cnb = 0;
+            for(unsigned i=1; i<n; ++i) {
+                cnb += std::pow(i,-beta);
+            }
+            return eo::rng.powerlaw(0,n,beta) / cnb;
+        }
+
+        double _beta;
+};
+
+#endif // _eoStandardBitMutation_h_

eo/src/utils/eoRNG.h

@@ -216,7 +216,7 @@ public :
             return uniform() < bias;
         }
 
-    /** Sample in a binomial dsitribution of size n and probability p.
+    /** Sample in a binomial distribution of size n and probability p.
 
         FIXME most naive algorithm, one should really use a rejection algorithm.
     */
@@ -229,6 +229,18 @@ public :
             return x;
         }
 
+    /** Sample in a power law distribution
+     */
+    double powerlaw(double min, double max, double gamma)
+       {
+            double x = uniform(min,max);
+            return std::pow(
+                     x * (std::pow(max,-gamma+1) - std::pow(min,-gamma+1))
+                       + std::pow(min,-gamma+1.0),
+                     1.0/(-gamma + 1.0)
+                   );
+       }
+
     /** Gaussian deviate
 
     Zero mean Gaussian deviate with standard deviation 1.

eo/test/CMakeLists.txt

@@ -75,6 +75,8 @@ set (TEST_LIST
   t-forge-algo
   t-algo-forged
   t-algo-forged-search
+  t-FastGA
+  t-eoFoundryFastGA
   )
 
 

eo/test/t-eoFoundryFastGA.cpp

@@ -0,0 +1,104 @@
+#include <iostream>
+#include <string>
+
+#include <eo>
+#include <ga.h>
+#include <utils/checkpointing>
+#include "../../problems/eval/oneMaxEval.h"
+
+using Particle = eoRealParticle<eoMaximizingFitness>;
+using Bits = eoBit<double>;
+
+// Generate a search space of 5,232,000 algorithms,
+// by enumerating candidate operators and their parameters.
+eoAlgoFoundryFastGA<Bits>& make_foundry(eoFunctorStore& store, eoInit<Bits>& init, eoEvalFunc<Bits>& eval_onemax)
+{
+    auto& foundry = store.pack< eoAlgoFoundryFastGA<Bits> >(init, eval_onemax, 20,10);
+
+    /***** Continuators ****/
+    for(size_t i=10; i < 100; i+=2 ) {
+        foundry.continuators.add< eoSteadyFitContinue<Bits> >(10,i);
+    }
+
+    for(double i=0.1; i<1.0; i+=0.1) {
+        foundry.crossover_rates.add<double>(i);
+        foundry.mutation_rates.add<double>(i);
+    }
+
+    for(size_t i=5; i<100; i+=10) {
+        foundry.pop_sizes.add<size_t>(i);
+    }
+
+    /***** Crossovers ****/
+    for(double i=0.1; i<0.9; i+=0.1) {
+        foundry.crossovers.add< eoUBitXover<Bits> >(i); // preference over 1
+    }
+    for(size_t i=1; i < 11; i+=1) {
+        foundry.crossovers.add< eoNPtsBitXover<Bits> >(i); // nb of points
+    }
+    foundry.crossovers.add< eo1PtBitXover<Bits> >();
+
+    /***** Mutations ****/
+    double p = 1.0; // Probability of flipping eath bit.
+    foundry.mutations.add< eoUniformBitMutation<Bits> >(p); // proba of flipping k bits, k drawn in uniform distrib
+    foundry.mutations.add< eoStandardBitMutation<Bits> >(p); // proba of flipping k bits, k drawn in binomial distrib
+    foundry.mutations.add< eoConditionalBitMutation<Bits> >(p); // proba of flipping k bits, k drawn in binomial distrib, minus zero
+    foundry.mutations.add< eoShiftedBitMutation<Bits> >(p); // proba of flipping k bits, k drawn in binomial distrib, changing zeros to one
+    foundry.mutations.add< eoNormalBitMutation<Bits> >(p); // proba of flipping k bits, k drawn in normal distrib
+    foundry.mutations.add< eoFastBitMutation<Bits> >(p); // proba of flipping k bits, k drawn in powerlaw distrib
+    for(size_t i=1; i < 11; i+=1) {
+        foundry.mutations.add< eoDetSingleBitFlip<Bits> >(i); // mutate k bits without duplicates
+    }
+
+    /***** Selectors *****/
+    foundry.selectors.add< eoRandomSelect<Bits> >();
+    foundry.selectors.add< eoSequentialSelect<Bits> >();
+    foundry.selectors.add< eoProportionalSelect<Bits> >();
+    for(size_t i=2; i < 10; i+=1) { // Tournament size.
+        foundry.selectors.add< eoDetTournamentSelect<Bits> >(i);
+    }
+    for(double i=0.51; i<0.91; i+=0.1) { // Tournament size as perc of pop.
+        foundry.selectors.add< eoStochTournamentSelect<Bits> >(i);
+    }
+
+    /***** Replacements ****/
+    foundry.replacements.add< eoPlusReplacement<Bits> >();
+    foundry.replacements.add< eoCommaReplacement<Bits> >();
+    foundry.replacements.add< eoSSGAWorseReplacement<Bits> >();
+    for(double i=0.51; i<0.91; i+=0.1) {
+        foundry.replacements.add< eoSSGAStochTournamentReplacement<Bits> >(i);
+    }
+    for(size_t i=2; i < 10; i+=1) {
+        foundry.replacements.add< eoSSGADetTournamentReplacement<Bits> >(i);
+    }
+
+    return foundry;
+}
+
+
+int main(int /*argc*/, char** /*argv*/)
+{
+    eo::log << eo::setlevel(eo::warnings);
+    eoFunctorStore store;
+
+    oneMaxEval<Bits> onemax_eval;
+
+    eoBooleanGenerator gen(0.5);
+    eoInitFixedLength<Bits> init(/*bitstring size=*/5, gen);
+
+    auto& foundry = make_foundry(store, init, onemax_eval);
+
+    size_t n = foundry.continuators.size() * foundry.crossovers.size() * foundry.mutations.size() * foundry.selectors.size() * foundry.replacements.size()* foundry.crossover_rates.size() * foundry.mutation_rates.size() * foundry.pop_sizes.size();
+    std::clog << n << " possible algorithms instances." << std::endl;
+
+    eoPop<Bits> pop;
+    pop.append(5,init);
+    ::apply(onemax_eval,pop);
+
+    foundry.select({0,0,0,0,0,0,0,0});
+    foundry(pop);
+
+    std::cout << "Done" << std::endl;
+    std::cout << pop << std::endl;
+    std::cout << pop.best_element() << std::endl;
+}