Merge branch 'master' of ssh://eodev.git.sourceforge.net/gitroot/eodev/eodev

edo/src/edo

@@ -59,6 +59,7 @@ Authors:
 #include "edoRepairer.h"
 #include "edoRepairerDispatcher.h"
 #include "edoRepairerRound.h"
+#include "edoRepairerModulo.h"
 #include "edoBounder.h"
 #include "edoBounderNo.h"
 #include "edoBounderBound.h"

edo/src/edoBounderUniform.h

@@ -49,6 +49,10 @@ public:
         assert( this->max().size() > 0 );
 
         assert( sol.size() > 0);
+        assert( sol.size() == this->min().size() );
+
+        eo::log << eo::debug << "BounderUniform: from sol = " << sol;
+        eo::log.flush();
 
         unsigned int size = sol.size();
         for (unsigned int d = 0; d < size; ++d) {
@@ -58,6 +62,8 @@ public:
                 sol[d] = rng.uniform( this->min()[d], this->max()[d] );
             }
         } // for d in size
+        
+        eo::log << eo::debug << "\tto sol = " << sol << std::endl;
     }
 };
 

edo/src/edoRepairerApply.h

@@ -0,0 +1,100 @@
+/*
+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) 2011 Thales group
+*/
+/*
+Authors:
+    Johann Dréo <johann.dreo@thalesgroup.com>
+*/
+
+#ifndef _edoRepairerApply_h
+#define _edoRepairerApply_h
+
+#include <algorithm>
+
+#include "edoRepairer.h"
+
+
+template < typename EOT, typename F = typename EOT::AtomType(typename EOT::AtomType) >
+class edoRepairerApply : public edoRepairer<EOT>
+{
+public:
+    edoRepairerApply( F function ) : _function(function) {}
+
+protected:
+    F * _function;
+};
+
+
+/** Apply an arbitrary unary function as a repairer on each item of the solution
+ *
+ * By default, the signature of the expected function is "EOT::AtomType(EOT::AtomType)"
+ *
+ * @ingroup Repairers
+ */
+template < typename EOT, typename F = typename EOT::AtomType(typename EOT::AtomType)>
+class edoRepairerApplyUnary : public edoRepairerApply<EOT,F>
+{
+public:
+    edoRepairerApplyUnary( F function ) : edoRepairerApply<EOT,F>(function) {}
+
+    virtual void operator()( EOT& sol )
+    {
+        std::transform( sol.begin(), sol.end(), sol.begin(), *(this->_function) );
+        sol.invalidate();
+    }
+};
+
+
+/** Apply an arbitrary binary function as a repairer on each item of the solution,
+ * the second argument of the function being fixed and given at instanciation.
+ *
+ * @see edoRepairerApplyUnary
+ *
+ * @ingroup Repairers
+ */
+template < typename EOT, typename F = typename EOT::AtomType(typename EOT::AtomType, typename EOT::AtomType)>
+class edoRepairerApplyBinary : public edoRepairerApply<EOT,F>
+{
+public:
+    typedef typename EOT::AtomType ArgType;
+
+    edoRepairerApplyBinary(   
+            F function, 
+            ArgType arg 
+        ) : edoRepairerApply<EOT,F>(function), _arg(arg) {}
+
+    virtual void operator()( EOT& sol )
+    {
+        // call the binary function on each item
+        // TODO find a way to use std::transform here? Or would it be too bloated?
+        for(typename EOT::iterator it = sol.begin(); it != sol.end(); ++it ) {
+            *it = (*(this->_function))( *it, _arg );
+        }
+        sol.invalidate();
+    }
+
+protected:
+    ArgType _arg;
+};
+
+
+#endif // !_edoRepairerApply_h
+

edo/src/edoRepairerDispatcher.h

@@ -38,23 +38,69 @@ Authors:
  * of indexes).
  *
  * Only work on EOT that implements the "push_back( EOT::AtomType )" and 
- * "operator[](uint)" and "at(uint)" methods.
+ * "operator[](uint)" and "at(uint)" methods (i.e. random access containers).
  *
  * Expects _addresses_ of the repairer operators.
  *
+ * Use the second template type if you want a different container to store
+ * indexes. You can use any iterable. For example, you may want to use a set if
+ * you need to be sure that indexes are use only once:
+ *     edoRepairerDispatcher<EOT, std::set<unsigned int> > rpd; 
+ *     std::set<unsigned int> idx(1,1);
+ *     idx.insert(2);
+ *     rpd.add( idx, &repairer );
+ *
+ * A diagram trying to visually explain how it works:
+ \ditaa
+
+       |
+ /-\   |   /------------\
+ | +---|---+ Dispatcher |
+ | |   v   |            |
+ | |+-----+| --------------------------------+
+ | || x_0 ||  +-+-+-+   |   +------------\   |   /-\
+ | |+-----+|  |2|3|5+*----*-* Repairer A +---|---+ |
+ | || x_1 ||  +-+-+-+   | | |            |   v   | |
+ | |+-----+|            | | |            |+-----+| |
+ | || x_2 ||            | | |            || x_2 || |
+ | |+-----+|            | | |            |+-----+| |
+ | || x_3 ||            | | |            || x_3 || |
+ | |+-----+|            | | |            |+-----+| |
+ | || x_4 ||            | | |            || x_5 || |           
+ | |+-----+|            | | |            |+-----+| |
+ | || x_5 ||            | | |            |   |   | |
+ | |+-----+|            | | |            +---|---+ |
+ | || x_6 ||            | | \------------/   |   \-/
+ | |+-----+| <-------------------------------+
+ | || x_7 ||            | |
+ | |+-----+|  +-+-+     | |
+ | || x_8 ||  |2|3+*------+
+ | |+-----+|  +-+-+     |
+ | || x_9 ||            |
+ | |+-----+|  +-+-+     |   +------------\       /-\
+ | |   |   |  |1|5+*--------* Repairer B +-------+ |
+ | |   |   |  +-+-+     |   |            |       | |
+ | |   |   |            |   |            |       | |
+ | |   |   |            |   |            +-------+ |
+ | +---|---+            |   \------------/       \-/
+ \-/   |   \------------/
+       v
+
+ \endditaa
+
  * @example t-dispatcher-round.cpp
  *
  * @ingroup Repairers
  */
-
+template < typename EOT, typename ICT = std::vector<unsigned int> >
-template < typename EOT >
 class edoRepairerDispatcher 
     : public edoRepairer<EOT>, 
              std::vector< 
-                  std::pair< std::vector< unsigned int >, edoRepairer< EOT >* > 
+                  std::pair< ICT, edoRepairer< EOT >* > 
              >
 {
 public:
+
     //! Empty constructor
     edoRepairerDispatcher() : 
         std::vector< 
@@ -63,7 +109,7 @@ public:
     {}
 
     //! Constructor with a single index set and repairer operator
-    edoRepairerDispatcher( std::vector<unsigned int> idx, edoRepairer<EOT>* op ) :
+    edoRepairerDispatcher( ICT idx, edoRepairer<EOT>* op ) :
         std::vector< 
             std::pair< std::vector< unsigned int >, edoRepairer< EOT >* > 
         >() 
@@ -72,7 +118,7 @@ public:
     }
 
     //! Add more indexes set and their corresponding repairer operator address to the list
-    void add( std::vector<unsigned int> idx, edoRepairer<EOT>* op )
+    void add( ICT idx, edoRepairer<EOT>* op )
     {
         assert( idx.size() > 0 );
         assert( op != NULL );
@@ -83,15 +129,27 @@ public:
     //! Repair a solution by calling several repair operator on subset of indexes
     virtual void operator()( EOT& sol )
     {
-        // ipair is an iterator that points on a pair
+//        std::cout << "in dispatcher, sol = " << sol << std::endl;
+
+        // ipair is an iterator that points on a pair of <indexes,repairer>
         for( typename edoRepairerDispatcher<EOT>::iterator ipair = this->begin(); ipair != this->end(); ++ipair ) {
+
+            assert( ipair->first.size() <= sol.size() ); // assert there is less indexes than items in the whole solution
+
             // a partial copy of the sol
             EOT partsol;
 
+//            std::cout << "\tusing indexes = ";
             // j is an iterator that points on an uint
             for( std::vector< unsigned int >::iterator j = ipair->first.begin(); j != ipair->first.end(); ++j ) {
+
+//                std::cout << *j << " ";
+//                std::cout.flush();
+
                 partsol.push_back( sol.at(*j) );
             } // for j
+//            std::cout << std::endl;
+//            std::cout << "\tpartial sol = " << partsol << std::endl;
 
             assert( partsol.size() > 0 );
 
@@ -108,6 +166,8 @@ public:
                 } // for j
             } // context for k
         } // for ipair
+        
+        sol.invalidate();
     }
 };
 

edo/src/edoRepairerModulo.h

@@ -0,0 +1,47 @@
+/*
+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) 2011 Thales group
+*/
+/*
+Authors:
+    Johann Dréo <johann.dreo@thalesgroup.com>
+*/
+
+#ifndef _edoRepairerModulo_h
+#define _edoRepairerModulo_h
+
+#include <cmath>
+
+#include "edoRepairerApply.h"
+
+/**
+ *
+ * @ingroup Repairers
+ */
+template < typename EOT >
+class edoRepairerModulo: public edoRepairerApplyBinary<EOT>
+{
+public:
+    edoRepairerModulo<EOT>( double denominator ) : edoRepairerApplyBinary<EOT>( std::fmod, denominator ) {}
+};
+
+
+#endif // !_edoRepairerModulo_h
+

edo/src/edoRepairerRound.h

@@ -30,39 +30,81 @@ Authors:
 
 #include <cmath>
 
-#include "edoRepairer.h"
+#include "edoRepairerApply.h"
 
-/**
+
+/** A repairer that calls "floor" on each items of a solution
+ *
+ * Just a proxy to "edoRepairerApplyUnary<EOT, EOT::AtomType(EOT::AtomType)> rep( std::floor);"
  *
  * @ingroup Repairers
  */
 template < typename EOT >
-class edoRepairerFloor : public edoRepairer<EOT>
+class edoRepairerFloor : public edoRepairerApplyUnary<EOT>
 {
 public:
-    virtual void operator()( EOT& sol )
+    edoRepairerFloor() : edoRepairerApplyUnary<EOT>( std::floor ) {}
-    {
-        for( unsigned int i=0; i < sol.size(); ++i ) {
-            sol[i] = floor( sol[i] );
-        }
-    }
 };
 
-/**
+
+/** A repairer that calls "ceil" on each items of a solution
+ *
+ * @see edoRepairerFloor
  *
  * @ingroup Repairers
  */
 template < typename EOT >
-class edoRepairerCeil : public edoRepairer<EOT>
+class edoRepairerCeil : public edoRepairerApplyUnary<EOT>
 {
 public:
-    virtual void operator()( EOT& sol )
+    edoRepairerCeil() : edoRepairerApplyUnary<EOT>( std::ceil ) {}
+};
+
+
+// FIXME find a way to put this function as a member of edoRepairerRoundDecimals
+template< typename ArgType >
+ArgType edoRound( ArgType val, ArgType prec = 1.0 )
+{ 
+    return (val > 0.0) ? 
+        floor(val * prec + 0.5) / prec : 
+         ceil(val * prec - 0.5) / prec ; 
+}
+
+/** A repairer that round values at a given a precision. 
+ *
+ * e.g. if prec=0.1, 8.06 will be rounded to 8.1
+ *
+ * @see edoRepairerFloor
+ * @see edoRepairerCeil
+ *
+ * @ingroup Repairers
+ */
+template < typename EOT >
+class edoRepairerRoundDecimals : public edoRepairerApplyBinary<EOT>
+{
+public:
+    typedef typename EOT::AtomType ArgType;
+
+    //! Generally speaking, we expect decimals being <= 1, but it can work for higher values
+    edoRepairerRoundDecimals( ArgType decimals ) : edoRepairerApplyBinary<EOT>( edoRound<ArgType>, 1 / decimals ) 
     {
-        for( unsigned int i=0; i < sol.size(); ++i ) {
+        assert( decimals <= 1.0 );
-            sol[i] = ceil( sol[i] );
+        assert( 1/decimals >= 1.0 );
-        }
     }
 };
 
 
+/** A repairer that do a rounding around val+0.5
+ *
+ * @see edoRepairerRoundDecimals
+ *
+ * @ingroup Repairers
+ */
+template < typename EOT >
+class edoRepairerRound : public edoRepairerRoundDecimals<EOT>
+{
+public:
+    edoRepairerRound() : edoRepairerRoundDecimals<EOT>( 1.0 ) {}
+};
+
 #endif // !_edoRepairerRound_h

edo/src/edoSamplerNormalMono.h

@@ -28,6 +28,8 @@ Authors:
 #ifndef _edoSamplerNormalMono_h
 #define _edoSamplerNormalMono_h
 
+#include <cmath>
+
 #include <utils/eoRNG.h>
 
 #include "edoSampler.h"
@@ -47,27 +49,25 @@ public:
 
     edoSamplerNormalMono( edoRepairer<EOT> & repairer ) : edoSampler< D >( repairer) {}
 
-    EOT sample( edoNormalMono< EOT >& distrib )
+    EOT sample( edoNormalMono<EOT>& distrib )
     {
         unsigned int size = distrib.size();
         assert(size > 0);
 
-        // Point we want to sample to get higher a set of points
+        // The point we want to draw
         // (coordinates in n dimension)
         // x = {x1, x2, ..., xn}
         EOT solution;
         
         // Sampling all dimensions
-        for (unsigned int i = 0; i < size; ++i)
+        for (unsigned int i = 0; i < size; ++i) {
-            {
             AtomType mean = distrib.mean()[i];
             AtomType variance = distrib.variance()[i];
-            AtomType random = rng.normal(mean, variance);
+            // should use the standard deviation, which have the same scale than the mean
-
+            AtomType random = rng.normal(mean, sqrt(variance) );
-            assert(variance >= 0);
 
             solution.push_back(random);
-            }
+        }
 
         return solution;
     }

edo/src/edoSamplerNormalMulti.h

@@ -28,151 +28,391 @@ Authors:
 #ifndef _edoSamplerNormalMulti_h
 #define _edoSamplerNormalMulti_h
 
+#include <cmath>
+#include <limits>
+
 #include <edoSampler.h>
 #include <boost/numeric/ublas/lu.hpp>
 #include <boost/numeric/ublas/symmetric.hpp>
 
-//! edoSamplerNormalMulti< EOT >
+/** Sample points in a multi-normal law defined by a mean vector and a covariance matrix.
-
+ *
-template< class EOT, typename D = edoNormalMulti< EOT > >
+ * Given M the mean vector and V the covariance matrix, of order n:
-class edoSamplerNormalMulti : public edoSampler< D >
+ *   - 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
+ */
+template< class EOT, typename EOD = edoNormalMulti< EOT > >
+class edoSamplerNormalMulti : public edoSampler< EOD >
 {
 public:
     typedef typename EOT::AtomType AtomType;
 
-    edoSamplerNormalMulti( edoRepairer<EOT> & repairer ) : edoSampler< D >( repairer) {}
 
+    /** Cholesky decomposition, given a matrix V, return a matrix L
+     * such as V = L L^T (L^T being the transposed of L).
+     *
+     * Need a symmetric and positive definite matrix as an input, which 
+     * should be the case of a non-ill-conditionned covariance matrix.
+     * Thus, expect a (lower) triangular matrix.
+     */
     class Cholesky
     {
     public:
-        Cholesky( const ublas::symmetric_matrix< AtomType, ublas::lower >& V)
+        //! The covariance-matrix is symetric
-        {
+        typedef ublas::symmetric_matrix< AtomType, ublas::lower > CovarMat;
-            unsigned int Vl = V.size1();
 
+        //! The factorization matrix is triangular
+        // FIXME check if triangular types behaviour is like having 0
+        typedef ublas::matrix< AtomType > FactorMat;
+
+        enum Method { 
+            //! use the standard algorithm, with square root @see factorize_LLT
+            standard,
+            //! use the algorithm using absolute value within the square root @see factorize_LLT_abs
+            absolute, 
+            //! use the method that set negative square roots to zero @see factorize_LLT_zero
+            zeroing,
+            //! use the robust algorithm, without square root @see factorize_LDLT
+            robust 
+        };
+        Method _use;
+
+
+        /** Instanciate without computing anything, you are responsible of 
+         * calling the algorithm and getting the result with operator()
+         * */
+        Cholesky( Cholesky::Method use = standard ) : _use(use) {}
+
+
+        /** Computation is made at instanciation and then cached in a member variable, 
+         * use decomposition() to get the result.
+         *
+         * Use the standard unstable method by default.
+         */
+        Cholesky(const CovarMat& V, Cholesky::Method use = standard ) :
+            _use(use), _L(ublas::zero_matrix<AtomType>(V.size1(),V.size2()))
+        {
+            factorize( V );
+        }
+
+
+        /** Compute the factorization and return the result 
+         */
+        const FactorMat& operator()( const CovarMat& V )
+        {
+            factorize( V );
+            return decomposition();
+        }
+
+        //! The decomposition of the covariance matrix
+        const FactorMat & decomposition() const 
+        {
+            return _L;
+        }
+
+        template<typename MT>
+        std::string format(const MT& mat )
+        {
+            std::ostringstream out;
+
+            for( unsigned int i=0; i<mat.size1(); ++i) {
+                out << std::endl;
+                for( unsigned int j=0; j<mat.size2(); ++j) {
+                    out << mat(i,j) << "\t";
+                } // columns
+            } // rows
+
+            return out.str();
+        }
+
+
+    protected:
+
+        //! The decomposition is a (lower) symetric matrix, just like the covariance matrix
+        FactorMat _L;
+       
+        /** Assert that the covariance matrix have the required properties and returns its dimension.
+         *
+         * Note: if compiled with NDEBUG, will not assert anything and just return the dimension.
+         */
+        unsigned assert_properties( const CovarMat& V )
+        {
+            unsigned int Vl = V.size1(); // number of lines
+
+            // the result goes in _L
+            _L = ublas::zero_matrix<AtomType>(Vl,Vl);
+
+            eo::log << eo::debug << std::endl << "Covariance matrix:" << format( V ) << std::endl;
+
+#ifndef NDEBUG
             assert(Vl > 0);
 
-            unsigned int Vc = V.size2();
+            unsigned int Vc = V.size2(); // number of columns
-
             assert(Vc > 0);
-
             assert( Vl == Vc );
 
-            _L.resize(Vl);
+            // partial assert that V is semi-positive definite
+            // assert that all diagonal elements are positives
+            for( unsigned int i=0; i < Vl; ++i ) {
+                assert( V(i,i) > 0 );
+            }
 
-            unsigned int i,j,k;
+            /* FIXME what is the more efficient way to check semi-positive definite? Candidates are:
+                 * perform the cholesky factorization
+                 * check if all eigenvalues are positives
+                 * check if all of the leading principal minors are positive
+                 */
+#endif
 
-            // first column
+            return Vl;
-            i=0;
+        }
 
-            // diagonal
+
-            j=0;
+        /** Actually performs the factorization with the method given at 
+         * instanciation. Results are cached in _L.
+         */
+        void factorize( const CovarMat& V )
+        {
+            if( _use == standard ) {
+                factorize_LLT( V );
+            } else if( _use == absolute ) {
+                factorize_LLT_abs( V );
+            } else if( _use == zeroing ) {
+                factorize_LLT_zero( V );
+            } else if( _use == robust ) {
+                factorize_LDLT( V );
+            }
+            
+            eo::log << eo::debug << std::endl << "Decomposed matrix:" << format( _L ) << std::endl;
+        }
+
+
+        /** This standard algorithm makes use of square root and is thus subject
+          * to round-off errors if the covariance matrix is very ill-conditioned.
+          *
+          * Compute L such that V = L L^T
+          *
+          * When compiled in debug mode and called on ill-conditionned matrix,
+          * will raise an assert before calling the square root on a negative number.
+          */
+        void factorize_LLT( const CovarMat& V)
+        {
+            unsigned int N = assert_properties( V );
+
+            unsigned int i=0, j=0, k;
             _L(0, 0) = sqrt( V(0, 0) );
 
             // end of the column
-            for ( j = 1; j < Vc; ++j )
+            for ( j = 1; j < N; ++j ) {
-            {
                 _L(j, 0) = V(0, j) / _L(0, 0);
             }
 
             // end of the matrix
-            for ( i = 1; i < Vl; ++i ) // each column
+            for ( i = 1; i < N; ++i ) { // each column
-            {
-
                 // diagonal
                 double sum = 0.0;
+                for ( k = 0; k < i; ++k) {
+                    sum += _L(i, k) * _L(i, k);
+                }
 
-                for ( k = 0; k < i; ++k)
+                // round-off errors may appear here
-                {
+                assert( V(i,i) - sum >= 0 );
+                _L(i,i) = sqrt( V(i,i) - sum );
+
+                for ( j = i + 1; j < N; ++j ) { // rows
+                    // one element
+                    sum = 0.0;
+                    for ( k = 0; k < i; ++k ) {
+                        sum += _L(j, k) * _L(i, k);
+                    }
+
+                    _L(j, i) = (V(j, i) - sum) / _L(i, i);
+
+                } // for j in ]i,N[
+            } // for i in [1,N[
+        }
+
+
+        /** This standard algorithm makes use of square root but do not fail
+          * if the covariance matrix is very ill-conditioned.
+          * Here, we propagate the error by using the absolute value before
+          * computing the square root.
+          *
+          * Be aware that this increase round-off errors, this is just a ugly
+          * hack to avoid crash.
+          */
+        void factorize_LLT_abs( const CovarMat & V)
+        {
+            unsigned int N = assert_properties( V );
+
+            unsigned int i=0, j=0, k;
+
+            _L(0, 0) = sqrt( V(0, 0) );
+
+            // end of the column
+            for ( j = 1; j < N; ++j ) {
+                _L(j, 0) = V(0, j) / _L(0, 0);
+            }
+
+            // end of the matrix
+            for ( i = 1; i < N; ++i ) { // each column
+                // diagonal
+                double sum = 0.0;
+                for ( k = 0; k < i; ++k) {
                     sum += _L(i, k) * _L(i, k);
                 }
 
                 _L(i,i) = sqrt( fabs( V(i,i) - sum) );
 
-                for ( j = i + 1; j < Vl; ++j ) // rows
+                for ( j = i + 1; j < N; ++j ) { // rows
-                {
                     // one element
                     sum = 0.0;
-
+                    for ( k = 0; k < i; ++k ) {
-                    for ( k = 0; k < i; ++k )
-                    {
                         sum += _L(j, k) * _L(i, k);
                     }
 
                     _L(j, i) = (V(j, i) - sum) / _L(i, i);
-                }
+
-            }
+                } // for j in ]i,N[
+            } // for i in [1,N[
         }
 
-        const ublas::symmetric_matrix< AtomType, ublas::lower >& get_L() const {return _L;}
 
-    private:
+        /** This standard algorithm makes use of square root but do not fail
-        ublas::symmetric_matrix< AtomType, ublas::lower > _L;
+          * if the covariance matrix is very ill-conditioned.
-    };
+          * Here, if the diagonal difference ir negative, we set it to zero.
+          *
+          * Be aware that this increase round-off errors, this is just a ugly
+          * hack to avoid crash.
+          */
+        void factorize_LLT_zero( const CovarMat & V)
+        {
+            unsigned int N = assert_properties( V );
 
-    edoSamplerNormalMulti( edoBounder< EOT > & bounder )
+            unsigned int i=0, j=0, k;
-        : edoSampler< edoNormalMulti< EOT > >( bounder )
-    {}
 
-    EOT sample( edoNormalMulti< EOT >& distrib )
+            _L(0, 0) = sqrt( V(0, 0) );
-    {
-        unsigned int size = distrib.size();
 
-        assert(size > 0);
+            // end of the column
-
+            for ( j = 1; j < N; ++j ) {
-
+                _L(j, 0) = V(0, j) / _L(0, 0);
-        //-------------------------------------------------------------
-        // Cholesky factorisation gererating matrix L from covariance
-        // matrix V.
-        // We must use cholesky.get_L() to get the resulting matrix.
-        //
-        // L = cholesky decomposition of varcovar
-        //-------------------------------------------------------------
-
-        Cholesky cholesky( distrib.varcovar() );
-        ublas::symmetric_matrix< AtomType, ublas::lower > L = cholesky.get_L();
-
-        //-------------------------------------------------------------
-
-
-        //-------------------------------------------------------------
-        // T = vector of size elements drawn in N(0,1) rng.normal(1.0)
-        //-------------------------------------------------------------
-
-        ublas::vector< AtomType > T( size );
-
-        for ( unsigned int i = 0; i < size; ++i )
-            {
-            T( i ) = rng.normal( 1.0 );
             }
 
-        //-------------------------------------------------------------
+            // end of the matrix
+            for ( i = 1; i < N; ++i ) { // each column
+                // diagonal
+                double sum = 0.0;
+                for ( k = 0; k < i; ++k) {
+                    sum += _L(i, k) * _L(i, k);
+                }
+
+                if(  V(i,i) - sum >= 0 ) {
+                    _L(i,i) = sqrt( V(i,i) - sum);
+                } else {
+                    _L(i,i) = std::numeric_limits<double>::epsilon();
+                }
+
+                for ( j = i + 1; j < N; ++j ) { // rows
+                    // one element
+                    sum = 0.0;
+                    for ( k = 0; k < i; ++k ) {
+                        sum += _L(j, k) * _L(i, k);
+                    }
+
+                    _L(j, i) = (V(j, i) - sum) / _L(i, i);
+
+                } // for j in ]i,N[
+            } // for i in [1,N[
+        }
 
 
-        //-------------------------------------------------------------
+        /** This alternative algorithm do not use square root in an inner loop,
-        // LT = prod( L, T )
+         * but only for some diagonal elements of the matrix D.
-        //-------------------------------------------------------------
+         *
+         * Computes L and D such as V = L D L^T. 
+         * The factorized matrix is (L D^1/2), because V = (L D^1/2) (L D^1/2)^T
+         */
+        void factorize_LDLT( const CovarMat& V)
+        {
+            // use "int" everywhere, because of the "j-1" operation
+            int N = assert_properties( V );
+            // example of an invertible matrix whose decomposition is undefined
+            assert( V(0,0) != 0 ); 
 
+            FactorMat L = ublas::zero_matrix<AtomType>(N,N);
+            FactorMat D = ublas::zero_matrix<AtomType>(N,N);
+            D(0,0) = V(0,0);
+            
+            for( int j=0; j<N; ++j ) { // each columns
+                L(j, j) = 1;
+
+                D(j,j) = V(j,j);
+                for( int k=0; k<=j-1; ++k) { // sum
+                    D(j,j) -= L(j,k) * L(j,k) * D(k,k);
+                }
+
+                for( int i=j+1; i<N; ++i ) { // remaining rows
+
+                    L(i,j) = V(i,j);
+                    for( int k=0; k<=j-1; ++k) { // sum
+                        L(i,j) -= L(i,k)*L(j,k) * D(k,k);
+                    }
+                    L(i,j) /= D(j,j);
+
+                } // for i in rows
+            } // for j in columns
+
+            // now compute the final symetric matrix: _L = L D^1/2
+            // remember that V = ( L D^1/2) ( L D^1/2)^T
+
+            // fortunately, the square root of a diagonal matrix is the square 
+            // root of all its elements
+            FactorMat D12 = D;
+            for(int i=0; i<N; ++i) {
+                D12(i,i) = sqrt(D(i,i));
+            }
+
+            // the factorization is thus _L*D^1/2
+            _L = ublas::prod( L, D12);
+        }
+        
+    }; // class Cholesky
+
+
+    edoSamplerNormalMulti( edoRepairer<EOT> & repairer, typename Cholesky::Method use = Cholesky::absolute ) 
+        : edoSampler< EOD >( repairer), _cholesky(use)
+    {}
+
+
+    EOT sample( EOD& distrib )
+    {
+        unsigned int size = distrib.size();
+        assert(size > 0);
+
+        // L = cholesky decomposition of varcovar
+        const typename Cholesky::FactorMat& L = _cholesky( distrib.varcovar() );
+
+        // T = vector of size elements drawn in N(0,1)
+        ublas::vector< AtomType > T( size );
+        for ( unsigned int i = 0; i < size; ++i ) {
+            T( i ) = rng.normal();
+        }
+
+        // LT = L * T
         ublas::vector< AtomType > LT = ublas::prod( L, T );
 
-        //-------------------------------------------------------------
-
-
-        //-------------------------------------------------------------
         // solution = means + LT
-        //-------------------------------------------------------------
-
         ublas::vector< AtomType > mean = distrib.mean();
-
         ublas::vector< AtomType > ublas_solution = mean + LT;
-
         EOT solution( size );
-
         std::copy( ublas_solution.begin(), ublas_solution.end(), solution.begin() );
 
-        //-------------------------------------------------------------
-
         return solution;
     }
+
+protected:
+    Cholesky _cholesky;
 };
 
 #endif // !_edoSamplerNormalMulti_h

edo/test/CMakeLists.txt

@@ -33,12 +33,14 @@ LINK_DIRECTORIES(${Boost_LIBRARY_DIRS})
 INCLUDE_DIRECTORIES(${CMAKE_SOURCE_DIR}/application/common)
 
 SET(SOURCES
+  t-cholesky
   t-edoEstimatorNormalMulti
   t-mean-distance
   t-bounderno
   t-uniform
   t-continue
   t-dispatcher-round
+  t-repairer-modulo
   )
 
 FOREACH(current ${SOURCES})

edo/test/t-cholesky.cpp

@@ -0,0 +1,250 @@
+
+/*
+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 <johann.dreo@thalesgroup.com>
+*/
+
+//#include <vector>
+#include <cstdlib>
+#include <iostream>
+#include <sstream>
+#include <limits>
+#include <iomanip>
+#include <ctime>
+
+#include <eo>
+#include <es.h>
+#include <edo>
+
+typedef eoReal< eoMinimizingFitness > EOT;
+typedef edoNormalMulti<EOT> EOD;
+
+
+void setformat( std::ostream& out )
+{
+    out << std::right;
+    out << std::setfill(' ');
+    out << std::setw( 5 + std::numeric_limits<double>::digits10);
+    out << std::setprecision(std::numeric_limits<double>::digits10);
+    out << std::setiosflags(std::ios_base::showpoint);
+}
+
+
+template<typename MT>
+std::string format(const MT& mat )
+{
+    std::ostringstream out;
+    setformat(out);
+
+    for( unsigned int i=0; i<mat.size1(); ++i) {
+        for( unsigned int j=0; j<mat.size2(); ++j) {
+            out << mat(i,j) << "\t";
+        } // columns
+        out << std::endl;
+    } // rows
+
+    return out.str();
+}
+
+
+template< typename T >
+T round( T val, T prec = 1.0 )
+{ 
+    return (val > 0.0) ? 
+        floor(val * prec + 0.5) / prec : 
+         ceil(val * prec - 0.5) / prec ; 
+}
+
+
+template<typename MT>
+bool equal( const MT& M1, const MT& M2, double prec /* = 1/std::numeric_limits<double>::digits10 ???*/ )
+{
+    if( M1.size1() != M2.size1() || M1.size2() != M2.size2() ) {
+        return false;
+    }
+
+    for( unsigned int i=0; i<M1.size1(); ++i ) {
+        for( unsigned int j=0; j<M1.size2(); ++j ) {
+            if( round(M1(i,j),prec) != round(M2(i,j),prec) ) {
+                std::cout << "round(M(" << i << "," << j << "," << prec << ") == " 
+                    << round(M1(i,j),prec) << " != " << round(M2(i,j),prec) << std::endl;
+                return false;
+            }
+        }
+    }
+
+    return true;
+}
+
+
+template<typename MT >
+MT error( const MT& M1, const MT& M2 )
+{
+    assert( M1.size1() == M2.size1() && M1.size1() == M2.size2() );
+
+    MT Err = ublas::zero_matrix<double>(M1.size1(),M1.size2());
+
+    for( unsigned int i=0; i<M1.size1(); ++i ) {
+        for( unsigned int j=0; j<M1.size2(); ++j ) {
+            Err(i,j) = M1(i,j) - M2(i,j);
+        }
+    }
+
+    return Err;
+}
+
+
+template<typename MT >
+double trigsum( const MT& M )
+{
+    double sum;
+    for( unsigned int i=0; i<M.size1(); ++i ) {
+        for( unsigned int j=i; j<M.size2(); ++j ) { // triangular browsing
+            sum += fabs( M(i,j) ); // absolute deviation
+        }
+    }
+    return sum;
+}
+
+
+template<typename T>
+double sum( const T& c )
+{
+     return std::accumulate(c.begin(), c.end(), 0);
+}
+
+
+int main(int argc, char** argv)
+{
+    srand(time(0));
+
+    unsigned int N = 4; // size of matrix
+    unsigned int R = 1000; // nb of repetitions
+
+    if( argc >= 2 ) {
+        N = std::atoi(argv[1]);
+    }
+    if( argc >= 3 ) {
+        R = std::atoi(argv[2]);
+    }
+
+    std::cout << "Usage: t-cholesky [matrix size] [repetitions]" << std::endl;
+    std::cout << "matrix size = " << N << std::endl;
+    std::cout << "repetitions = " << R << std::endl;
+
+    typedef edoSamplerNormalMulti<EOT,EOD>::Cholesky::CovarMat CovarMat;
+    typedef edoSamplerNormalMulti<EOT,EOD>::Cholesky::FactorMat FactorMat;
+
+    edoSamplerNormalMulti<EOT,EOD>::Cholesky LLT( edoSamplerNormalMulti<EOT,EOD>::Cholesky::standard );
+    edoSamplerNormalMulti<EOT,EOD>::Cholesky LLTa( edoSamplerNormalMulti<EOT,EOD>::Cholesky::absolute );
+    edoSamplerNormalMulti<EOT,EOD>::Cholesky LLTz( edoSamplerNormalMulti<EOT,EOD>::Cholesky::zeroing );
+    edoSamplerNormalMulti<EOT,EOD>::Cholesky LDLT( edoSamplerNormalMulti<EOT,EOD>::Cholesky::robust );
+
+    std::vector<double> s0,s1,s2,s3;
+    for( unsigned int n=0; n<R; ++n ) {
+
+        // a variance-covariance matrix of size N*N
+        CovarMat V(N,N);
+
+        // random covariance matrix
+        for( unsigned int i=0; i<N; ++i) {
+            V(i,i) = std::pow(rand(),2); // variance should be >= 0
+            for( unsigned int j=i+1; j<N; ++j) {
+                V(i,j) = rand();
+            }
+        }
+
+        FactorMat L0 = LLT(V);
+        CovarMat V0 = ublas::prod( L0, ublas::trans(L0) );
+        s0.push_back( trigsum(error(V,V0)) );
+
+        FactorMat L1 = LLTa(V);
+        CovarMat V1 = ublas::prod( L1, ublas::trans(L1) );
+        s1.push_back( trigsum(error(V,V1)) );
+
+        FactorMat L2 = LLTz(V);
+        CovarMat V2 = ublas::prod( L2, ublas::trans(L2) );
+        s2.push_back( trigsum(error(V,V2)) );
+
+        FactorMat L3 = LDLT(V);
+        CovarMat V3 = ublas::prod( L3, ublas::trans(L3) );
+        s3.push_back( trigsum(error(V,V3)) );
+    }
+
+    std::cout << "Average error:" << std::endl;
+    std::cout << "\tLLT:  " << sum(s0)/R << std::endl;
+    std::cout << "\tLLTa: " << sum(s1)/R << std::endl;
+    std::cout << "\tLLTz: " << sum(s2)/R << std::endl;
+    std::cout << "\tLDLT: " << sum(s3)/R << std::endl;
+
+//    double precision = 1e-15;
+//    if( argc >= 4 ) {
+//        precision = std::atof(argv[3]);
+//    }
+//    std::cout << "precision = " << precision << std::endl;
+//    std::cout << "usage: t-cholesky [N] [precision]" << std::endl;
+//    std::cout << "N = " << N << std::endl;
+//    std::cout << "precision = " << precision << std::endl;
+//    std::string linesep = "--------------------------------------------------------------------------------------------";
+//    std::cout << linesep << std::endl;
+//
+//    setformat(std::cout);
+//
+//    std::cout << "Covariance matrix" << std::endl << format(V) << std::endl;
+//    std::cout << linesep << std::endl;
+//
+//    edoSamplerNormalMulti<EOT,EOD>::Cholesky LLT( edoSamplerNormalMulti<EOT,EOD>::Cholesky::standard );
+//    FactorMat L0 = LLT(V);
+//    CovarMat V0 = ublas::prod( L0, ublas::trans(L0) );
+//    CovarMat E0 = error(V,V0);
+//    std::cout << "LLT" << std::endl << format(E0) << std::endl;
+//    std::cout << trigsum(E0) << std::endl;
+//    std::cout << "LLT" << std::endl << format(L0) << std::endl;
+//    std::cout << "LLT covar" << std::endl << format(V0) << std::endl;
+//    assert( equal(V0,V,precision) );
+//    std::cout << linesep << std::endl;
+//
+//    edoSamplerNormalMulti<EOT,EOD>::Cholesky LLTa( edoSamplerNormalMulti<EOT,EOD>::Cholesky::absolute );
+//    FactorMat L1 = LLTa(V);
+//    CovarMat V1 = ublas::prod( L1, ublas::trans(L1) );
+//    CovarMat E1 = error(V,V1);
+//    std::cout << "LLT abs" << std::endl << format(E1) << std::endl;
+//    std::cout << trigsum(E1) << std::endl;
+//    std::cout << "LLT abs" << std::endl << format(L1) << std::endl;
+//    std::cout << "LLT covar" << std::endl << format(V1) << std::endl;
+//    assert( equal(V1,V,precision) );
+//    std::cout << linesep << std::endl;
+//    
+//    edoSamplerNormalMulti<EOT,EOD>::Cholesky LDLT( edoSamplerNormalMulti<EOT,EOD>::Cholesky::robust );
+//    FactorMat L2 = LDLT(V);
+//    CovarMat V2 = ublas::prod( L2, ublas::trans(L2) );
+//    CovarMat E2 = error(V,V2);
+//    std::cout << "LDLT" << std::endl << format(E2) << std::endl;
+//    std::cout << trigsum(E2) << std::endl;
+//    std::cout << "LDLT" << std::endl << format(L2) << std::endl;
+//    std::cout << "LDLT covar" << std::endl << format(V2) << std::endl;
+//    assert( equal(V2,V,precision) );
+//    std::cout << linesep << std::endl;
+    
+}

edo/test/t-dispatcher-round.cpp

@@ -38,10 +38,18 @@ int main(void)
     sol.push_back(1.1);
     sol.push_back(3.9);
     sol.push_back(3.9);
-    // we expect {1,2,3,4}
+    sol.push_back(5.4);
+    sol.push_back(5.6);
+    sol.push_back(7.011);
+    sol.push_back(8.09);
+    sol.push_back(8.21);
+    
+    std::cout << "expect: INVALID  9 1 2 3 4 5 6 7 8.1 8.2" << std::endl;
 
     edoRepairer<EOT>* rep1 = new edoRepairerFloor<EOT>();
     edoRepairer<EOT>* rep2 = new edoRepairerCeil<EOT>();
+    edoRepairer<EOT>* rep3 = new edoRepairerRound<EOT>();
+    edoRepairer<EOT>* rep4 = new edoRepairerRoundDecimals<EOT>( 10 );
 
     std::vector<unsigned int> indexes1;
     indexes1.push_back(0);
@@ -51,8 +59,19 @@ int main(void)
     indexes2.push_back(1);
     indexes2.push_back(3);
 
+    std::vector<unsigned int> indexes3;
+    indexes3.push_back(4);
+    indexes3.push_back(5);
+
+    std::vector<unsigned int> indexes4;
+    indexes4.push_back(6);
+    indexes4.push_back(7);
+    indexes4.push_back(8);
+
     edoRepairerDispatcher<EOT> repare( indexes1, rep1 );
     repare.add( indexes2, rep2 );
+    repare.add( indexes3, rep3 );
+    repare.add( indexes4, rep4 );
 
     repare(sol);
 

edo/test/t-repairer-modulo.cpp

@@ -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 <johann.dreo@thalesgroup.com>
+*/
+
+#define _USE_MATH_DEFINES
+#include <math.h>
+
+#include <eo>
+#include <edo>
+#include <es.h>
+
+typedef eoReal< eoMinimizingFitness > EOT;
+
+int main(void)
+{
+    EOT sol;
+    sol.push_back( M_PI * 1 );
+    sol.push_back( M_PI * 2 );
+    sol.push_back( M_PI * 3 );
+    sol.push_back( M_PI * 4 );
+    sol.push_back( M_PI * 4 + M_PI / 2 );
+    sol.push_back( M_PI * 5 + M_PI / 2 );
+    // we expect {pi,0,pi,0,pi/2,pi+pi/2}
+    std::cout << "expect: INVALID  4 3.14159 0 3.14159 0 1.5708 4.71239" << std::endl;
+
+    edoRepairer<EOT>* repare = new edoRepairerModulo<EOT>( 2 * M_PI ); // modulo 2pi
+
+    (*repare)(sol);
+
+    std::cout << sol << std::endl;
+
+    return 0;
+}

eo/src/utils/eoRNG.h

@@ -216,6 +216,7 @@ public :
     /** Gaussian deviate
 
     Zero mean Gaussian deviate with standard deviation 1.
+    Note: Use the Marsaglia polar method.
 
     @return Random Gaussian deviate
     */