From ccd77037ec2c202592ecda2e03648633da35f9a7 Mon Sep 17 00:00:00 2001 From: nojhan Date: Tue, 29 Nov 2011 15:06:19 +0100 Subject: [PATCH 01/10] bugfix: IQR should not be negative (oups :-) --- eo/src/utils/eoStat.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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; } } From e77b27257937c559c223c9bbd1c9b51f006335db Mon Sep 17 00:00:00 2001 From: nojhan Date: Thu, 1 Dec 2011 10:20:17 +0100 Subject: [PATCH 02/10] bad doc string for eoSIGContinue --- eo/src/eoSIGContinue.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From 5115deb7dd255d4a7443732aad0e8a6dc209acc4 Mon Sep 17 00:00:00 2001 From: nojhan Date: Mon, 12 Dec 2011 22:22:04 +0100 Subject: [PATCH 03/10] an eoStat that simply count generations --- eo/src/utils/checkpointing | 1 + eo/src/utils/eoGenCounter.h | 46 +++++++++++++++++++++++++++++++++++++ 2 files changed, 47 insertions(+) create mode 100644 eo/src/utils/eoGenCounter.h 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 From 3b66f04fd6f896053a6817788a96545953dd7dec Mon Sep 17 00:00:00 2001 From: nojhan Date: Wed, 14 Dec 2011 15:49:27 +0100 Subject: [PATCH 04/10] refactoring Cholesky decomposition in several classes --- edo/src/utils/edoCholesky.h | 280 ++++++++++++++++++++++++++++++++++++ 1 file changed, 280 insertions(+) create mode 100644 edo/src/utils/edoCholesky.h diff --git a/edo/src/utils/edoCholesky.h b/edo/src/utils/edoCholesky.h new file mode 100644 index 00000000..ea0c2bca --- /dev/null +++ b/edo/src/utils/edoCholesky.h @@ -0,0 +1,280 @@ +/* +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() + * */ + Cholesky( 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. + */ + Cholesky(const CovarMat& V) : + _L(ublas::zero_matrix(V.size1(),V.size2())) + { + (*this)( V ); + } + + /** Compute the factorization and cache the result */ + virtual void operator()( const CovarMat& V ) = 0; + + /** Compute the factorization and return the result */ + virtual const FactorMat& operator()( const CovarMat& V ) + { + (*this)( 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 +{ + virtual void operator()( 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) = 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 += _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[ + } + + inline virtual L_i_i( const 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 +{ + inline virtual L_i_i( const 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 +{ + inline virtual L_i_i( const 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 +{ + virtual void operator()( 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 Date: Wed, 14 Dec 2011 15:53:57 +0100 Subject: [PATCH 05/10] beautifying code for Cholesky classes --- edo/src/utils/edoCholesky.h | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/edo/src/utils/edoCholesky.h b/edo/src/utils/edoCholesky.h index ea0c2bca..be17f5fc 100644 --- a/edo/src/utils/edoCholesky.h +++ b/edo/src/utils/edoCholesky.h @@ -129,6 +129,7 @@ protected: template< typename T > class CholeskyLLT : public CholeskyBase { +public: virtual void operator()( const CovarMat& V ) { unsigned int N = assert_properties( V ); @@ -164,6 +165,8 @@ class CholeskyLLT : public CholeskyBase } // for i in [1,N[ } + + /** The step of the standard LLT algorithm where round off errors may appear */ inline virtual L_i_i( const CovarMat& V, const unsigned int& i, const double& sum ) const { // round-off errors may appear here @@ -173,7 +176,6 @@ class CholeskyLLT : public CholeskyBase }; - /** 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 @@ -185,6 +187,7 @@ class CholeskyLLT : public CholeskyBase template< typename T > class CholeskyLLTabs : public CholeskyLLT { +public: inline virtual L_i_i( const CovarMat& V, const unsigned int& i, const double& sum ) const { /***** ugly hack *****/ @@ -203,6 +206,7 @@ class CholeskyLLTabs : public CholeskyLLT template< typename T > class CholeskyLLTzero : public CholeskyLLT { +public: inline virtual L_i_i( const CovarMat& V, const unsigned int& i, const double& sum ) const { T Lii; @@ -226,6 +230,7 @@ class CholeskyLLTzero : public CholeskyLLT template< typename T > class CholeskyLDLT : public CholeskyBase { +public: virtual void operator()( const CovarMat& V ) { // use "int" everywhere, because of the "j-1" operation @@ -256,24 +261,24 @@ class CholeskyLDLT : public CholeskyBase } // for i in rows } // for j in columns - _L = compute_L( L, D ); + _L = root( L, D ); } - inline FactorMat compute_L( FactorMat& L, FactorMat& D ) + inline FactorMat root( FactorMat& L, FactorMat& D ) { // 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; + FactorMat sqrt_D = D; for(int i=0; i Date: Fri, 23 Dec 2011 16:05:05 +0100 Subject: [PATCH 06/10] cholesky classes are now in separate files --- edo/src/edoSamplerNormalMulti.h | 310 -------------------------------- 1 file changed, 310 deletions(-) diff --git a/edo/src/edoSamplerNormalMulti.h b/edo/src/edoSamplerNormalMulti.h index 6099ed95..a017d493 100644 --- a/edo/src/edoSamplerNormalMulti.h +++ b/edo/src/edoSamplerNormalMulti.h @@ -48,316 +48,6 @@ 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; - } - - 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(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; - } - - - /** 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 ); - } - } - - - /** 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) = 0; - } - - 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) From 4805f72a806a48b79fe115586ab4e2c51f454c6f Mon Sep 17 00:00:00 2001 From: nojhan Date: Fri, 23 Dec 2011 16:31:35 +0100 Subject: [PATCH 07/10] use the new cholesky api in edoSamplerNormalMulti --- edo/src/edoSamplerNormalMulti.h | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/edo/src/edoSamplerNormalMulti.h b/edo/src/edoSamplerNormalMulti.h index e0516685..b7e73c75 100644 --- a/edo/src/edoSamplerNormalMulti.h +++ b/edo/src/edoSamplerNormalMulti.h @@ -32,6 +32,7 @@ Authors: #include #include +#include #include #include @@ -48,8 +49,8 @@ class edoSamplerNormalMulti : public edoSampler< EOD > public: typedef typename EOT::AtomType AtomType; - edoSamplerNormalMulti( edoRepairer & repairer, typename Cholesky::Method use = Cholesky::absolute ) - : edoSampler< EOD >( repairer), _cholesky(use) + edoSamplerNormalMulti( edoRepairer & repairer ) + : edoSampler< EOD >( repairer) {} @@ -59,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 ); @@ -80,7 +81,7 @@ public: } protected: - Cholesky _cholesky; + cholesky::CholeskyLLT _cholesky; }; #endif // !_edoSamplerNormalMulti_h From ebf33d81772c384c5db93a0bbf4df7fedd1af53c Mon Sep 17 00:00:00 2001 From: nojhan Date: Fri, 23 Dec 2011 16:50:51 +0100 Subject: [PATCH 08/10] bugfixes concerning the mew api --- edo/src/utils/edoCholesky.h | 46 ++++++++++++++++++------------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/edo/src/utils/edoCholesky.h b/edo/src/utils/edoCholesky.h index be17f5fc..e0437d00 100644 --- a/edo/src/utils/edoCholesky.h +++ b/edo/src/utils/edoCholesky.h @@ -48,26 +48,26 @@ public: /** Instanciate without computing anything, you are responsible of * calling the algorithm and getting the result with operator() * */ - Cholesky( size_t s1 = 1, size_t s2 = 1 ) : + 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. */ - Cholesky(const CovarMat& V) : + CholeskyBase(const CovarMat& V) : _L(ublas::zero_matrix(V.size1(),V.size2())) { (*this)( V ); } /** Compute the factorization and cache the result */ - virtual void operator()( const CovarMat& V ) = 0; + virtual void factorize( const CovarMat& V ) = 0; /** Compute the factorization and return the result */ virtual const FactorMat& operator()( const CovarMat& V ) { - (*this)( V ); + this->factorize( V ); return decomposition(); } @@ -130,16 +130,16 @@ template< typename T > class CholeskyLLT : public CholeskyBase { public: - virtual void operator()( const CovarMat& V ) + virtual void factorize( const typename CholeskyBase::CovarMat& V ) { unsigned int N = assert_properties( V ); unsigned int i=0, j=0, k; - _L(0, 0) = sqrt( V(0, 0) ); + this->_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); + this->_L(j, 0) = V(0, j) / this->_L(0, 0); } // end of the matrix @@ -147,19 +147,19 @@ public: // diagonal double sum = 0.0; for ( k = 0; k < i; ++k) { - sum += _L(i, k) * _L(i, k); + sum += this->_L(i, k) * this->_L(i, k); } - _L(i,i) = L_i_i( V, i, sum ); + 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 += _L(j, k) * _L(i, k); + sum += this->_L(j, k) * this->_L(i, k); } - _L(j, i) = (V(j, i) - sum) / _L(i, i); + this->_L(j, i) = (V(j, i) - sum) / this->_L(i, i); } // for j in ]i,N[ } // for i in [1,N[ @@ -167,7 +167,7 @@ public: /** The step of the standard LLT algorithm where round off errors may appear */ - inline virtual L_i_i( const CovarMat& V, const unsigned int& i, const double& sum ) const + 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 ); @@ -188,7 +188,7 @@ template< typename T > class CholeskyLLTabs : public CholeskyLLT { public: - inline virtual L_i_i( const CovarMat& V, const unsigned int& i, const double& sum ) const + 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) ); @@ -207,7 +207,7 @@ template< typename T > class CholeskyLLTzero : public CholeskyLLT { public: - inline virtual L_i_i( const CovarMat& V, const unsigned int& i, const double& sum ) const + 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 ) { @@ -231,15 +231,15 @@ template< typename T > class CholeskyLDLT : public CholeskyBase { public: - virtual void operator()( const CovarMat& V ) + 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 ); - FactorMat L = ublas::zero_matrix(N,N); - FactorMat D = ublas::zero_matrix(N,N); + 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 FactorMat root( FactorMat& L, FactorMat& D ) + inline typename CholeskyBase::FactorMat root( typename CholeskyBase::FactorMat& L, typename CholeskyBase::FactorMat& D ) { - // now compute the final symetric matrix: _L = L D^1/2 + // 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 - FactorMat sqrt_D = D; - for(int i=0; i::FactorMat sqrt_D = D; + for( int i=0; i < D.size1(); ++i) { sqrt_D(i,i) = sqrt(D(i,i)); } - // the factorization is thus _L*D^1/2 + // the factorization is thus this->_L*D^1/2 return ublas::prod( L, sqrt_D ); } }; From d10325d1ad38fb324f073b1ae7c94485362a92d9 Mon Sep 17 00:00:00 2001 From: nojhan Date: Thu, 12 Jan 2012 16:49:50 +0100 Subject: [PATCH 09/10] bypass repairer calls if the index list is empty in the dispatcher --- edo/src/edoRepairerDispatcher.h | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/edo/src/edoRepairerDispatcher.h b/edo/src/edoRepairerDispatcher.h index 117e8dfb..1bf9cb39 100644 --- a/edo/src/edoRepairerDispatcher.h +++ b/edo/src/edoRepairerDispatcher.h @@ -131,6 +131,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 +141,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 +154,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 From 10a5712594cc74d5089feb3b3907613642508b8f Mon Sep 17 00:00:00 2001 From: nojhan Date: Thu, 12 Jan 2012 17:47:55 +0100 Subject: [PATCH 10/10] replace the assert on empty index in the dispatcher by a warning --- edo/src/edoRepairerDispatcher.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/edo/src/edoRepairerDispatcher.h b/edo/src/edoRepairerDispatcher.h index 1bf9cb39..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) );