diff --git a/edo/src/edoRepairerDispatcher.h b/edo/src/edoRepairerDispatcher.h index 117e8dfb..b2c0777b 100644 --- a/edo/src/edoRepairerDispatcher.h +++ b/edo/src/edoRepairerDispatcher.h @@ -120,7 +120,10 @@ public: //! Add more indexes set and their corresponding repairer operator address to the list void add( ICT idx, edoRepairer* op ) { - assert( idx.size() > 0 ); + //assert( idx.size() > 0 ); +#ifndef NDEBUG + eo::log << eo::warnings << "A repairer is added to the dispatcher while having an empty index list, nothing will be repaired" << std::endl; +#endif assert( op != NULL ); this->push_back( std::make_pair(idx, op) ); @@ -131,6 +134,7 @@ public: { // std::cout << "in dispatcher, sol = " << sol << std::endl; + // iterate over { indexe, repairer } // ipair is an iterator that points on a pair of for( typename edoRepairerDispatcher::iterator ipair = this->begin(); ipair != this->end(); ++ipair ) { @@ -140,6 +144,8 @@ public: EOT partsol; // std::cout << "\tusing indexes = "; +// + // iterate over 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 ) { @@ -151,6 +157,9 @@ public: // std::cout << std::endl; // std::cout << "\tpartial sol = " << partsol << std::endl; + if( partsol.size() == 0 ) { + continue; + } assert( partsol.size() > 0 ); // apply the repairer on the partial copy diff --git a/edo/src/edoSamplerNormalMulti.h b/edo/src/edoSamplerNormalMulti.h index dc3519eb..b7e73c75 100644 --- a/edo/src/edoSamplerNormalMulti.h +++ b/edo/src/edoSamplerNormalMulti.h @@ -32,6 +32,7 @@ Authors: #include #include +#include #include #include @@ -48,340 +49,8 @@ class edoSamplerNormalMulti : public edoSampler< EOD > public: typedef typename EOT::AtomType AtomType; - - /** 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: - //! The covariance-matrix is symetric - typedef ublas::symmetric_matrix< AtomType, ublas::lower > CovarMat; - - //! 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(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 - std::string format(const MT& mat ) - { - std::ostringstream out; - - for( unsigned int i=0; i(Vl,Vl); - - eo::log << eo::debug << std::endl << "Covariance matrix:" << format( V ) << std::endl; - -#ifndef NDEBUG - assert(Vl > 0); - - unsigned int Vc = V.size2(); // number of columns - assert(Vc > 0); - assert( Vl == Vc ); - - // 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 ); - } - - /* 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 - - return Vl; - } - - - /** 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 < 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); - } - - // 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 < 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, 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 ); - - 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); - } - - if( V(i,i) - sum >= 0 ) { - _L(i,i) = sqrt( V(i,i) - sum); - } else { - _L(i,i) = std::numeric_limits::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, - * 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(N,N); - FactorMat D = ublas::zero_matrix(N,N); - D(0,0) = V(0,0); - - for( int j=0; j & repairer, typename Cholesky::Method use = Cholesky::absolute ) - : edoSampler< EOD >( repairer), _cholesky(use) + edoSamplerNormalMulti( edoRepairer & repairer ) + : edoSampler< EOD >( repairer) {} @@ -391,7 +60,7 @@ public: assert(size > 0); // L = cholesky decomposition of varcovar - const typename Cholesky::FactorMat& L = _cholesky( distrib.varcovar() ); + const typename cholesky::CholeskyBase::FactorMat& L = _cholesky( distrib.varcovar() ); // T = vector of size elements drawn in N(0,1) ublas::vector< AtomType > T( size ); @@ -412,7 +81,7 @@ public: } protected: - Cholesky _cholesky; + cholesky::CholeskyLLT _cholesky; }; #endif // !_edoSamplerNormalMulti_h diff --git a/edo/src/utils/edoCholesky.h b/edo/src/utils/edoCholesky.h new file mode 100644 index 00000000..e0437d00 --- /dev/null +++ b/edo/src/utils/edoCholesky.h @@ -0,0 +1,285 @@ +/* +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 + Caner Candan +*/ + +namespace cholesky { + +/** 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. + */ +template< typename T > +class CholeskyBase +{ +public: + //! The covariance-matrix is symetric + typedef ublas::symmetric_matrix< T, ublas::lower > CovarMat; + + //! The factorization matrix is triangular + // FIXME check if triangular types behaviour is like having 0 + typedef ublas::matrix< T > FactorMat; + + /** Instanciate without computing anything, you are responsible of + * calling the algorithm and getting the result with operator() + * */ + CholeskyBase( size_t s1 = 1, size_t s2 = 1 ) : + _L(ublas::zero_matrix(s1,s2)) + {} + + /** Computation is made at instanciation and then cached in a member variable, + * use decomposition() to get the result. + */ + CholeskyBase(const CovarMat& V) : + _L(ublas::zero_matrix(V.size1(),V.size2())) + { + (*this)( V ); + } + + /** Compute the factorization and cache the result */ + virtual void factorize( const CovarMat& V ) = 0; + + /** Compute the factorization and return the result */ + virtual const FactorMat& operator()( const CovarMat& V ) + { + this->factorize( V ); + return decomposition(); + } + + //! The decomposition of the covariance matrix + const FactorMat & decomposition() const + { + return _L; + } + +protected: + + /** 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(Vl,Vl); + +#ifndef NDEBUG + assert(Vl > 0); + + unsigned int Vc = V.size2(); // number of columns + assert(Vc > 0); + assert( Vl == Vc ); + + // 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 ); + } + + /* 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 + + return Vl; + } + + //! The decomposition is a (lower) symetric matrix, just like the covariance matrix + FactorMat _L; +}; + + +/** 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. + */ +template< typename T > +class CholeskyLLT : public CholeskyBase +{ +public: + virtual void factorize( const typename CholeskyBase::CovarMat& V ) + { + unsigned int N = assert_properties( V ); + + unsigned int i=0, j=0, k; + this->_L(0, 0) = sqrt( V(0, 0) ); + + // end of the column + for ( j = 1; j < N; ++j ) { + this->_L(j, 0) = V(0, j) / this->_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 += this->_L(i, k) * this->_L(i, k); + } + + this->_L(i,i) = L_i_i( V, i, sum ); + + for ( j = i + 1; j < N; ++j ) { // rows + // one element + sum = 0.0; + for ( k = 0; k < i; ++k ) { + sum += this->_L(j, k) * this->_L(i, k); + } + + this->_L(j, i) = (V(j, i) - sum) / this->_L(i, i); + + } // for j in ]i,N[ + } // for i in [1,N[ + } + + + /** The step of the standard LLT algorithm where round off errors may appear */ + inline virtual T L_i_i( const typename CholeskyBase::CovarMat& V, const unsigned int& i, const double& sum ) const + { + // round-off errors may appear here + assert( V(i,i) - sum >= 0 ); + return sqrt( V(i,i) - sum ); + } +}; + + +/** 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. + */ +template< typename T > +class CholeskyLLTabs : public CholeskyLLT +{ +public: + inline virtual T L_i_i( const typename CholeskyBase::CovarMat& V, const unsigned int& i, const double& sum ) const + { + /***** ugly hack *****/ + return sqrt( fabs( V(i,i) - sum) ); + } +}; + + +/** This standard algorithm makes use of square root but do not fail + * 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. + */ +template< typename T > +class CholeskyLLTzero : public CholeskyLLT +{ +public: + inline virtual T L_i_i( const typename CholeskyBase::CovarMat& V, const unsigned int& i, const double& sum ) const + { + T Lii; + if( V(i,i) - sum >= 0 ) { + Lii = sqrt( V(i,i) - sum); + } else { + /***** ugly hack *****/ + Lii = 0; + } + return Lii; + } +}; + + +/** This alternative algorithm do not use square root in an inner loop, + * 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 + */ +template< typename T > +class CholeskyLDLT : public CholeskyBase +{ +public: + virtual void factorize( const typename CholeskyBase::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 ); + + typename CholeskyBase::FactorMat L = ublas::zero_matrix(N,N); + typename CholeskyBase::FactorMat D = ublas::zero_matrix(N,N); + D(0,0) = V(0,0); + + for( int j=0; j_L = root( L, D ); + } + + + inline typename CholeskyBase::FactorMat root( typename CholeskyBase::FactorMat& L, typename CholeskyBase::FactorMat& D ) + { + // now compute the final symetric matrix: this->_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 + typename CholeskyBase::FactorMat sqrt_D = D; + for( int i=0; i < D.size1(); ++i) { + sqrt_D(i,i) = sqrt(D(i,i)); + } + + // the factorization is thus this->_L*D^1/2 + return ublas::prod( L, sqrt_D ); + } +}; + +} // namespace cholesky diff --git a/eo/CMakeLists.txt b/eo/CMakeLists.txt index 3cab6563..c10ac3c5 100644 --- a/eo/CMakeLists.txt +++ b/eo/CMakeLists.txt @@ -16,7 +16,7 @@ INCLUDE(eo-conf.cmake OPTIONAL) PROJECT(EO) # CMake > 2.8 is needed, because of the FindOpenMP feature -cmake_minimum_required(VERSION 2.8) +#cmake_minimum_required(VERSION 2.8) #SET(PROJECT_VERSION_MAJOR 1) #SET(PROJECT_VERSION_MINOR 1) diff --git a/eo/src/eoSIGContinue.h b/eo/src/eoSIGContinue.h index d7db55d3..2199b49e 100644 --- a/eo/src/eoSIGContinue.h +++ b/eo/src/eoSIGContinue.h @@ -48,7 +48,7 @@ void set_bool(int) } /** - Ctrl C handling: this eoContinue tells whether the user pressed Ctrl C + A continuator that stops if a given signal is received during the execution */ template< class EOT> class eoSIGContinue: public eoContinue diff --git a/eo/src/utils/checkpointing b/eo/src/utils/checkpointing index dea1432c..5d49e589 100644 --- a/eo/src/utils/checkpointing +++ b/eo/src/utils/checkpointing @@ -17,6 +17,7 @@ #include #include #include +#include // and make_help - any better suggestion to include it? void make_help(eoParser & _parser); diff --git a/eo/src/utils/eoGenCounter.h b/eo/src/utils/eoGenCounter.h new file mode 100644 index 00000000..a06949c2 --- /dev/null +++ b/eo/src/utils/eoGenCounter.h @@ -0,0 +1,46 @@ +/* + 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 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + (c) Thales group 2011 + + Author: johann.dreo@thalesgroup.com + */ + +#ifndef _eoGenCounter_h +#define _eoGenCounter_h + +#include +#include + +/** + An eoStat that simply gives the current generation index + + @ingroup Stats +*/ +class eoGenCounter : public eoUpdater, public eoValueParam +{ +public: + eoGenCounter( unsigned int start = 0, std::string label = "Gen" ) : eoValueParam(start, label), _nb(start) {} + + virtual void operator()() + { + value() = _nb++; + } + +private: + unsigned int _nb; +}; + +#endif diff --git a/eo/src/utils/eoStat.h b/eo/src/utils/eoStat.h index 916bcf43..32c3bf9f 100644 --- a/eo/src/utils/eoStat.h +++ b/eo/src/utils/eoStat.h @@ -448,7 +448,7 @@ public: std::nth_element( pop.begin(), pop.begin()+quartile*3, pop.end() ); typename EOT::Fitness Q3 = pop[quartile*3].fitness(); - value() = Q1 - Q3; + value() = Q3 - Q1; } }