From 7acaa0beab5af4e5d739625c8cbe2db24adb8972 Mon Sep 17 00:00:00 2001 From: nojhan Date: Wed, 26 Oct 2011 15:35:04 +0200 Subject: [PATCH 01/24] more doc about the repairer dispatcher --- edo/src/edoRepairerDispatcher.h | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/edo/src/edoRepairerDispatcher.h b/edo/src/edoRepairerDispatcher.h index 5c4aa7d8..4a20693b 100644 --- a/edo/src/edoRepairerDispatcher.h +++ b/edo/src/edoRepairerDispatcher.h @@ -38,23 +38,31 @@ 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 > rpd; + * std::set idx(1,1); + * idx.insert(2); + * rpd.add( idx, &repairer ); + * * @example t-dispatcher-round.cpp * * @ingroup Repairers */ - -template < typename EOT > +template < typename EOT, typename ICT = std::vector > class edoRepairerDispatcher : public edoRepairer, std::vector< - std::pair< std::vector< unsigned int >, edoRepairer< EOT >* > + std::pair< ICT, edoRepairer< EOT >* > > { public: + //! Empty constructor edoRepairerDispatcher() : std::vector< @@ -63,7 +71,7 @@ public: {} //! Constructor with a single index set and repairer operator - edoRepairerDispatcher( std::vector idx, edoRepairer* op ) : + edoRepairerDispatcher( IndexContainer idx, edoRepairer* op ) : std::vector< std::pair< std::vector< unsigned int >, edoRepairer< EOT >* > >() @@ -72,7 +80,7 @@ public: } //! Add more indexes set and their corresponding repairer operator address to the list - void add( std::vector idx, edoRepairer* op ) + void add( IndexContainer idx, edoRepairer* op ) { assert( idx.size() > 0 ); assert( op != NULL ); From e6b1c63fde7e22c3eb7a8e9f4b9d3fc7a02d8b0e Mon Sep 17 00:00:00 2001 From: nojhan Date: Fri, 28 Oct 2011 16:50:15 +0200 Subject: [PATCH 02/24] assert there is less indexes than items in the whole solution given to the repairer dispatcher ; comment useless debug code --- edo/src/edoRepairerDispatcher.h | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/edo/src/edoRepairerDispatcher.h b/edo/src/edoRepairerDispatcher.h index 5c4aa7d8..5dffb3f9 100644 --- a/edo/src/edoRepairerDispatcher.h +++ b/edo/src/edoRepairerDispatcher.h @@ -83,15 +83,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 for( typename edoRepairerDispatcher::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 ); From d002dccadf1a2d97ca93702fbae7d3828506479d Mon Sep 17 00:00:00 2001 From: nojhan Date: Thu, 3 Nov 2011 10:54:07 +0100 Subject: [PATCH 03/24] cleaner debug code for uniform bounder --- edo/src/edoBounderUniform.h | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/edo/src/edoBounderUniform.h b/edo/src/edoBounderUniform.h index c90317c7..f5ea97c7 100644 --- a/edo/src/edoBounderUniform.h +++ b/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; } }; From 03b69f181919ba60fef1c22194ce729e3c47eea0 Mon Sep 17 00:00:00 2001 From: nojhan Date: Thu, 3 Nov 2011 11:13:32 +0100 Subject: [PATCH 04/24] compilfix: correct template name for dispatcher's indexes --- edo/src/edoRepairerDispatcher.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/edo/src/edoRepairerDispatcher.h b/edo/src/edoRepairerDispatcher.h index 0ee65a78..b198a14a 100644 --- a/edo/src/edoRepairerDispatcher.h +++ b/edo/src/edoRepairerDispatcher.h @@ -71,7 +71,7 @@ public: {} //! Constructor with a single index set and repairer operator - edoRepairerDispatcher( IndexContainer idx, edoRepairer* op ) : + edoRepairerDispatcher( ICT idx, edoRepairer* op ) : std::vector< std::pair< std::vector< unsigned int >, edoRepairer< EOT >* > >() @@ -80,7 +80,7 @@ public: } //! Add more indexes set and their corresponding repairer operator address to the list - void add( IndexContainer idx, edoRepairer* op ) + void add( ICT idx, edoRepairer* op ) { assert( idx.size() > 0 ); assert( op != NULL ); From 979f33b72154ff041ce0eec6535bd01acf5eccef Mon Sep 17 00:00:00 2001 From: nojhan Date: Mon, 7 Nov 2011 15:44:36 +0100 Subject: [PATCH 05/24] ditaa text diagram explaining how the repairer dispatcher works --- edo/src/edoRepairerDispatcher.h | 38 +++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/edo/src/edoRepairerDispatcher.h b/edo/src/edoRepairerDispatcher.h index b198a14a..90ee7924 100644 --- a/edo/src/edoRepairerDispatcher.h +++ b/edo/src/edoRepairerDispatcher.h @@ -50,6 +50,44 @@ Authors: * 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 From 4d45970767f0ccb6a0773b14f2f1c29fe15ba317 Mon Sep 17 00:00:00 2001 From: nojhan Date: Tue, 8 Nov 2011 17:55:01 +0100 Subject: [PATCH 06/24] refactored vectors repairers that are using a function in a more generic way, added a binary function wrapper for vectors repairers, with the modulus function as an example --- edo/src/edo | 1 + edo/src/edoRepairerApply.h | 98 ++++++++++++++++++++++++++++++++++ edo/src/edoRepairerModulo.h | 47 ++++++++++++++++ edo/src/edoRepairerRound.h | 28 ++++------ edo/test/CMakeLists.txt | 1 + edo/test/t-repairer-modulo.cpp | 55 +++++++++++++++++++ 6 files changed, 213 insertions(+), 17 deletions(-) create mode 100644 edo/src/edoRepairerApply.h create mode 100644 edo/src/edoRepairerModulo.h create mode 100644 edo/test/t-repairer-modulo.cpp diff --git a/edo/src/edo b/edo/src/edo index 69f13f11..c03acd35 100644 --- a/edo/src/edo +++ b/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" diff --git a/edo/src/edoRepairerApply.h b/edo/src/edoRepairerApply.h new file mode 100644 index 00000000..fbadacda --- /dev/null +++ b/edo/src/edoRepairerApply.h @@ -0,0 +1,98 @@ +/* +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 +*/ + +#ifndef _edoRepairerApply_h +#define _edoRepairerApply_h + +#include + +#include "edoRepairer.h" + + +template < typename EOT, typename F = typename EOT::AtomType(typename EOT::AtomType) > +class edoRepairerApply : public edoRepairer +{ +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 +{ +public: + edoRepairerApplyUnary( F function ) : edoRepairerApply(function) {} + + virtual void operator()( EOT& sol ) + { + std::transform( sol.begin(), sol.end(), sol.begin(), *(this->_function) ); + } +}; + + +/** 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 +{ +public: + typedef typename EOT::AtomType ArgType; + + edoRepairerApplyBinary( + F function, + ArgType arg + ) : edoRepairerApply(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 ); + } + } + +protected: + ArgType _arg; +}; + + +#endif // !_edoRepairerApply_h + diff --git a/edo/src/edoRepairerModulo.h b/edo/src/edoRepairerModulo.h new file mode 100644 index 00000000..21367333 --- /dev/null +++ b/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 +*/ + +#ifndef _edoRepairerModulo_h +#define _edoRepairerModulo_h + +#include + +#include "edoRepairerApply.h" + +/** + * + * @ingroup Repairers + */ +template < typename EOT > +class edoRepairerModulo: public edoRepairerApplyBinary +{ +public: + edoRepairerModulo( double denominator ) : edoRepairerApplyBinary( std::fmod, denominator ) {} +}; + + +#endif // !_edoRepairerModulo_h + diff --git a/edo/src/edoRepairerRound.h b/edo/src/edoRepairerRound.h index 03216942..21dd1cae 100644 --- a/edo/src/edoRepairerRound.h +++ b/edo/src/edoRepairerRound.h @@ -30,38 +30,32 @@ Authors: #include -#include "edoRepairer.h" +#include "edoRepairerApply.h" -/** +/** A repairer that calls "floor" on each items of a solution + * + * Just a proxy to "edoRepairerApplyUnary rep( std::floor);" * * @ingroup Repairers */ template < typename EOT > -class edoRepairerFloor : public edoRepairer +class edoRepairerFloor : public edoRepairerApplyUnary { public: - virtual void operator()( EOT& sol ) - { - for( unsigned int i=0; i < sol.size(); ++i ) { - sol[i] = floor( sol[i] ); - } - } + edoRepairerFloor() : edoRepairerApplyUnary( std::floor ) {} }; -/** +/** A repairer that calls "ceil" on each items of a solution + * + * @see edoRepairerFloor * * @ingroup Repairers */ template < typename EOT > -class edoRepairerCeil : public edoRepairer +class edoRepairerCeil : public edoRepairerApplyUnary { public: - virtual void operator()( EOT& sol ) - { - for( unsigned int i=0; i < sol.size(); ++i ) { - sol[i] = ceil( sol[i] ); - } - } + edoRepairerCeil() : edoRepairerApplyUnary( std::ceil ) {} }; diff --git a/edo/test/CMakeLists.txt b/edo/test/CMakeLists.txt index acdabc7b..61fdd68f 100644 --- a/edo/test/CMakeLists.txt +++ b/edo/test/CMakeLists.txt @@ -39,6 +39,7 @@ SET(SOURCES t-uniform t-continue t-dispatcher-round + t-repairer-modulo ) FOREACH(current ${SOURCES}) diff --git a/edo/test/t-repairer-modulo.cpp b/edo/test/t-repairer-modulo.cpp new file mode 100644 index 00000000..7a495f31 --- /dev/null +++ b/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 +*/ + +#define _USE_MATH_DEFINES +#include + +#include +#include +#include + +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* repare = new edoRepairerModulo( 2 * M_PI ); // modulo 2pi + + (*repare)(sol); + + std::cout << sol << std::endl; + + return 0; +} From 0e9d9f8c81c4c9fe2e99a61293125989c53886c8 Mon Sep 17 00:00:00 2001 From: nojhan Date: Tue, 8 Nov 2011 18:40:28 +0100 Subject: [PATCH 07/24] added a symetric round repairer around 0.5, plus a rounding at a given precision --- edo/src/edoRepairerRound.h | 43 +++++++++++++++++++++++++++++++++ edo/test/t-dispatcher-round.cpp | 21 +++++++++++++++- 2 files changed, 63 insertions(+), 1 deletion(-) diff --git a/edo/src/edoRepairerRound.h b/edo/src/edoRepairerRound.h index 21dd1cae..dc05716d 100644 --- a/edo/src/edoRepairerRound.h +++ b/edo/src/edoRepairerRound.h @@ -32,6 +32,7 @@ Authors: #include "edoRepairerApply.h" + /** A repairer that calls "floor" on each items of a solution * * Just a proxy to "edoRepairerApplyUnary rep( std::floor);" @@ -45,6 +46,7 @@ public: edoRepairerFloor() : edoRepairerApplyUnary( std::floor ) {} }; + /** A repairer that calls "ceil" on each items of a solution * * @see edoRepairerFloor @@ -59,4 +61,45 @@ public: }; +// 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=10, 8.06 will be rounded to 8.1 + * + * @see edoRepairerFloor + * @see edoRepairerCeil + * + * @ingroup Repairers + */ +template < typename EOT > +class edoRepairerRoundDecimals : public edoRepairerApplyBinary +{ +public: + typedef typename EOT::AtomType ArgType; + + edoRepairerRoundDecimals( ArgType prec ) : edoRepairerApplyBinary( edoRound, prec ) {} +}; + + +/** A repairer that do a rounding around val+0.5 + * + * @see edoRepairerRoundDecimals + * + * @ingroup Repairers + */ +template < typename EOT > +class edoRepairerRound : public edoRepairerRoundDecimals +{ +public: + edoRepairerRound() : edoRepairerRoundDecimals( 1.0 ) {} +}; + #endif // !_edoRepairerRound_h diff --git a/edo/test/t-dispatcher-round.cpp b/edo/test/t-dispatcher-round.cpp index 89fc1644..cfcd9261 100644 --- a/edo/test/t-dispatcher-round.cpp +++ b/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* rep1 = new edoRepairerFloor(); edoRepairer* rep2 = new edoRepairerCeil(); + edoRepairer* rep3 = new edoRepairerRound(); + edoRepairer* rep4 = new edoRepairerRoundDecimals( 10 ); std::vector indexes1; indexes1.push_back(0); @@ -51,8 +59,19 @@ int main(void) indexes2.push_back(1); indexes2.push_back(3); + std::vector indexes3; + indexes3.push_back(4); + indexes3.push_back(5); + + std::vector indexes4; + indexes4.push_back(6); + indexes4.push_back(7); + indexes4.push_back(8); + edoRepairerDispatcher repare( indexes1, rep1 ); repare.add( indexes2, rep2 ); + repare.add( indexes3, rep3 ); + repare.add( indexes4, rep4 ); repare(sol); From 391c9e9e221ee2febba37c51befc39517d734e62 Mon Sep 17 00:00:00 2001 From: nojhan Date: Wed, 9 Nov 2011 10:15:21 +0100 Subject: [PATCH 08/24] use decimals instead of tens --- edo/src/edoRepairerRound.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/edo/src/edoRepairerRound.h b/edo/src/edoRepairerRound.h index dc05716d..c370a03c 100644 --- a/edo/src/edoRepairerRound.h +++ b/edo/src/edoRepairerRound.h @@ -72,7 +72,7 @@ ArgType edoRound( ArgType val, ArgType prec = 1.0 ) /** A repairer that round values at a given a precision. * - * e.g. if prec=10, 8.06 will be rounded to 8.1 + * e.g. if prec=0.1, 8.06 will be rounded to 8.1 * * @see edoRepairerFloor * @see edoRepairerCeil @@ -85,7 +85,8 @@ class edoRepairerRoundDecimals : public edoRepairerApplyBinary public: typedef typename EOT::AtomType ArgType; - edoRepairerRoundDecimals( ArgType prec ) : edoRepairerApplyBinary( edoRound, prec ) {} + //! Generally speaking, we expect decimals being <= 1, but it can work for higher values + edoRepairerRoundDecimals( ArgType decimals ) : edoRepairerApplyBinary( edoRound, 1 / decimals ) {} }; From 96dc505f92919639eee28c4b855603c6b7acb993 Mon Sep 17 00:00:00 2001 From: nojhan Date: Wed, 9 Nov 2011 10:35:41 +0100 Subject: [PATCH 09/24] bugfix: invalidate the individual after having repaired it --- edo/src/edoRepairerApply.h | 2 ++ edo/src/edoRepairerDispatcher.h | 2 ++ 2 files changed, 4 insertions(+) diff --git a/edo/src/edoRepairerApply.h b/edo/src/edoRepairerApply.h index fbadacda..ca46f344 100644 --- a/edo/src/edoRepairerApply.h +++ b/edo/src/edoRepairerApply.h @@ -58,6 +58,7 @@ public: virtual void operator()( EOT& sol ) { std::transform( sol.begin(), sol.end(), sol.begin(), *(this->_function) ); + sol.invalidate(); } }; @@ -87,6 +88,7 @@ public: for(typename EOT::iterator it = sol.begin(); it != sol.end(); ++it ) { *it = (*(this->_function))( *it, _arg ); } + sol.invalidate(); } protected: diff --git a/edo/src/edoRepairerDispatcher.h b/edo/src/edoRepairerDispatcher.h index 90ee7924..117e8dfb 100644 --- a/edo/src/edoRepairerDispatcher.h +++ b/edo/src/edoRepairerDispatcher.h @@ -166,6 +166,8 @@ public: } // for j } // context for k } // for ipair + + sol.invalidate(); } }; From 1af6715ddc3e3fb48eee3846d55979a77890ecec Mon Sep 17 00:00:00 2001 From: nojhan Date: Wed, 9 Nov 2011 11:38:40 +0100 Subject: [PATCH 10/24] BUGFIX: should use the standard deviation when sampling a mono-normal law, because it have the same scale than the mean --- edo/src/edoSamplerNormalMono.h | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/edo/src/edoSamplerNormalMono.h b/edo/src/edoSamplerNormalMono.h index b3162f92..b9340b17 100644 --- a/edo/src/edoSamplerNormalMono.h +++ b/edo/src/edoSamplerNormalMono.h @@ -28,6 +28,8 @@ Authors: #ifndef _edoSamplerNormalMono_h #define _edoSamplerNormalMono_h +#include + #include #include "edoSampler.h" @@ -47,27 +49,25 @@ public: edoSamplerNormalMono( edoRepairer & repairer ) : edoSampler< D >( repairer) {} - EOT sample( edoNormalMono< EOT >& distrib ) + EOT sample( edoNormalMono& 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); - - assert(variance >= 0); + // should use the standard deviation, which have the same scale than the mean + AtomType random = rng.normal(mean, sqrt(variance) ); solution.push_back(random); - } + } return solution; } From 9fe6995df16e0c5f231b33051893102950557665 Mon Sep 17 00:00:00 2001 From: nojhan Date: Wed, 9 Nov 2011 15:48:07 +0100 Subject: [PATCH 11/24] refactoring of the cholesky decomposition --- edo/src/edoSamplerNormalMulti.h | 160 ++++++++++++++++++-------------- 1 file changed, 90 insertions(+), 70 deletions(-) diff --git a/edo/src/edoSamplerNormalMulti.h b/edo/src/edoSamplerNormalMulti.h index 7c08d22b..e5e73a01 100644 --- a/edo/src/edoSamplerNormalMulti.h +++ b/edo/src/edoSamplerNormalMulti.h @@ -28,12 +28,19 @@ Authors: #ifndef _edoSamplerNormalMulti_h #define _edoSamplerNormalMulti_h +#include + #include #include #include -//! edoSamplerNormalMulti< EOT > - +/** Sample points in a multi-normal law defined by a mean vector and a covariance matrix. + * + * Given M the mean vector and V the covariance matrix, of order n: + * - 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 D = edoNormalMulti< EOT > > class edoSamplerNormalMulti : public edoSampler< D > { @@ -42,135 +49,148 @@ public: edoSamplerNormalMulti( edoRepairer & repairer ) : edoSampler< D >( repairer) {} + /** Cholesky decomposition, given a matrix V, return a matrix L + * such as V = L Lt (Lt being the conjugate transpose 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) - { - unsigned int Vl = V.size1(); + private: + //! The decomposition is a (lower) symetric matrix, just like the covariance matrix + ublas::symmetric_matrix< AtomType, ublas::lower > _L; + public: + //! The decomposition of the covariance matrix + const ublas::symmetric_matrix< AtomType, ublas::lower >& decomposition() const {return _L;} + + /** Computation is made at instanciation and then cached in a member variable, + * use decomposition() to get the result. + */ + Cholesky( const ublas::symmetric_matrix< AtomType, ublas::lower >& V ) + { + factorize( V ); + } + + /** 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 ublas::symmetric_matrix< AtomType, ublas::lower >& V ) + { + unsigned int Vl = V.size1(); // number of lines assert(Vl > 0); - unsigned int Vc = V.size2(); - + unsigned int Vc = V.size2(); // number of columns assert(Vc > 0); - assert( Vl == Vc ); + // FIXME assert definite semi-positive + + // the result goes in _L _L.resize(Vl); - unsigned int i,j,k; + return Vl; + } - // first column - i=0; + /** This standard algorithm makes use of square root and is thus subject + * to round-off errors if the covariance matrix is very ill-conditioned. + */ + void factorize( const ublas::symmetric_matrix< AtomType, ublas::lower >& V) + { + unsigned int N = assert_properties( V ); - // diagonal - j=0; + 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) - { + for ( k = 0; k < i; ++k) { sum += _L(i, k) * _L(i, k); } - _L(i,i) = sqrt( fabs( V(i,i) - sum) ); + // round-off errors may appear here + assert( V(i,i) - sum >= 0 ); + _L(i,i) = sqrt( V(i,i) - sum ); + //_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: - ublas::symmetric_matrix< AtomType, ublas::lower > _L; - }; + /** This alternative algorithm does not use square root BUT the covariance + * matrix must be invertible. + * + * Computes L and D such as V = L D Lt + */ + /* + void factorize_robust( const ublas::symmetric_matrix< AtomType, ublas::lower >& V) + { + unsigned int N = assert_properties( V ); + + unsigned int i, j, k; + ublas::symmetric_matrix< AtomType, ublas::lower > D = ublas::zero_matrix(N); + _L(0, 0) = sqrt( V(0, 0) ); + + } + */ + + + }; // class Cholesky + edoSamplerNormalMulti( edoBounder< EOT > & bounder ) : edoSampler< edoNormalMulti< EOT > >( bounder ) {} + EOT sample( edoNormalMulti< EOT >& distrib ) { unsigned int size = distrib.size(); - assert(size > 0); - - //------------------------------------------------------------- // Cholesky factorisation gererating matrix L from covariance // matrix V. - // We must use cholesky.get_L() to get the resulting matrix. + // We must use cholesky.decomposition() to get the resulting matrix. // // L = cholesky decomposition of varcovar - //------------------------------------------------------------- - Cholesky cholesky( distrib.varcovar() ); - ublas::symmetric_matrix< AtomType, ublas::lower > L = cholesky.get_L(); + ublas::symmetric_matrix< AtomType, ublas::lower > L = cholesky.decomposition(); - //------------------------------------------------------------- - - - //------------------------------------------------------------- // 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(); + } - for ( unsigned int i = 0; i < size; ++i ) - { - T( i ) = rng.normal( 1.0 ); - } - - //------------------------------------------------------------- - - - //------------------------------------------------------------- - // LT = prod( L, T ) - //------------------------------------------------------------- - + // 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; } }; From b04a66aaea662391b090cc2d4a985bf455ccf168 Mon Sep 17 00:00:00 2001 From: nojhan Date: Wed, 9 Nov 2011 15:49:04 +0100 Subject: [PATCH 12/24] Note: normal rng use the Marsaglia polar method --- eo/src/utils/eoRNG.h | 1 + 1 file changed, 1 insertion(+) diff --git a/eo/src/utils/eoRNG.h b/eo/src/utils/eoRNG.h index 9885eff5..984e7627 100644 --- a/eo/src/utils/eoRNG.h +++ b/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 */ From fb82958253b11e40811f6aacecee216e13714e07 Mon Sep 17 00:00:00 2001 From: nojhan Date: Thu, 10 Nov 2011 11:48:14 +0100 Subject: [PATCH 13/24] more refactorization of the multi-normal sampler, now you can choose among two Cholesky factorization methods --- edo/src/edoSamplerNormalMulti.h | 140 ++++++++++++++++++++++++++++---- 1 file changed, 123 insertions(+), 17 deletions(-) diff --git a/edo/src/edoSamplerNormalMulti.h b/edo/src/edoSamplerNormalMulti.h index e5e73a01..1bff937b 100644 --- a/edo/src/edoSamplerNormalMulti.h +++ b/edo/src/edoSamplerNormalMulti.h @@ -53,52 +53,115 @@ public: * such as V = L Lt (Lt being the conjugate transpose of L). * * Need a symmetric and positive definite matrix as an input, which - * should be the case of a non-ill-conditionned covariance matrix. + * should be the case of a non-ill-conditionned covariance matrix. * Thus, expect a (lower) triangular matrix. */ class Cholesky { - private: - //! The decomposition is a (lower) symetric matrix, just like the covariance matrix - ublas::symmetric_matrix< AtomType, ublas::lower > _L; - public: - //! The decomposition of the covariance matrix - const ublas::symmetric_matrix< AtomType, ublas::lower >& decomposition() const {return _L;} + typedef ublas::symmetric_matrix< AtomType, ublas::lower > MatrixType; + + enum Method { + //! use the standard algorithm, with square root @see factorize + standard, + //! use the algorithm using absolute value within the square root @see factorize_abs + absolute, + //! use the robust algorithm, without square root + //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 ublas::symmetric_matrix< AtomType, ublas::lower >& V ) + Cholesky(const MatrixType& V, Cholesky::Method use = standard ) : _use(use) + { + factsorize( V ); + } + + + /** Compute the factorization and return the result + */ + const MatrixType& operator()( const MatrixType& V ) { factorize( V ); + return decomposition(); } + //! The decomposition of the covariance matrix + const MatrixType & decomposition() const {return _L;} + + protected: + + //! The decomposition is a (lower) symetric matrix, just like the covariance matrix + MatrixType _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 ublas::symmetric_matrix< AtomType, ublas::lower >& V ) + unsigned assert_properties( const MatrixType& V ) { unsigned int Vl = V.size1(); // number of lines + + // the result goes in _L + _L.resize(Vl); + +#ifndef NDEBUG assert(Vl > 0); unsigned int Vc = V.size2(); // number of columns assert(Vc > 0); assert( Vl == Vc ); - // FIXME assert definite semi-positive + // 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 ); + } - // the result goes in _L - _L.resize(Vl); + /* 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 MatrixType& V ) + { + if( _use == standard ) { + factorize_std( V ); + } else if( _use == absolute ) { + factorize_abs( 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. + * + * 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( const ublas::symmetric_matrix< AtomType, ublas::lower >& V) + void factorize_std( const MatrixType& V) { unsigned int N = assert_properties( V ); @@ -137,24 +200,67 @@ public: } + /** 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_abs( const MatrixType & 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 alternative algorithm does not use square root BUT the covariance * matrix must be invertible. * * Computes L and D such as V = L D Lt */ /* - void factorize_robust( const ublas::symmetric_matrix< AtomType, ublas::lower >& V) + void factorize_robust( const MatrixType& V) { unsigned int N = assert_properties( V ); unsigned int i, j, k; ublas::symmetric_matrix< AtomType, ublas::lower > D = ublas::zero_matrix(N); + _L(0, 0) = sqrt( V(0, 0) ); } */ - }; // class Cholesky @@ -173,8 +279,8 @@ public: // We must use cholesky.decomposition() to get the resulting matrix. // // L = cholesky decomposition of varcovar - Cholesky cholesky( distrib.varcovar() ); - ublas::symmetric_matrix< AtomType, ublas::lower > L = cholesky.decomposition(); + Cholesky cholesky( Cholesky::absolute ); + const typename Cholesky::MatrixType& L = cholesky( distrib.varcovar() ); // T = vector of size elements drawn in N(0,1) rng.normal(1.0) ublas::vector< AtomType > T( size ); From ce7f59d408b79c6b48eb7adbcb827cf54a528213 Mon Sep 17 00:00:00 2001 From: nojhan Date: Thu, 10 Nov 2011 14:07:58 +0100 Subject: [PATCH 14/24] pass the Cholesky method from the multi-normal sampler to the decomposition sub-class --- edo/src/edoSamplerNormalMulti.h | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/edo/src/edoSamplerNormalMulti.h b/edo/src/edoSamplerNormalMulti.h index 1bff937b..5f87152c 100644 --- a/edo/src/edoSamplerNormalMulti.h +++ b/edo/src/edoSamplerNormalMulti.h @@ -47,7 +47,6 @@ class edoSamplerNormalMulti : public edoSampler< D > public: typedef typename EOT::AtomType AtomType; - edoSamplerNormalMulti( edoRepairer & repairer ) : edoSampler< D >( repairer) {} /** Cholesky decomposition, given a matrix V, return a matrix L * such as V = L Lt (Lt being the conjugate transpose of L). @@ -264,8 +263,8 @@ public: }; // class Cholesky - edoSamplerNormalMulti( edoBounder< EOT > & bounder ) - : edoSampler< edoNormalMulti< EOT > >( bounder ) + edoSamplerNormalMulti( edoRepairer & repairer, typename Cholesky::Method use = Cholesky::absolute ) + : edoSampler< D >( repairer), _cholesky(use) {} @@ -279,8 +278,7 @@ public: // We must use cholesky.decomposition() to get the resulting matrix. // // L = cholesky decomposition of varcovar - Cholesky cholesky( Cholesky::absolute ); - const typename Cholesky::MatrixType& L = cholesky( distrib.varcovar() ); + const typename Cholesky::MatrixType& L = _cholesky( distrib.varcovar() ); // T = vector of size elements drawn in N(0,1) rng.normal(1.0) ublas::vector< AtomType > T( size ); @@ -299,6 +297,9 @@ public: return solution; } + +protected: + Cholesky _cholesky; }; #endif // !_edoSamplerNormalMulti_h From b9b3c486f95b80f27d0591f3dea4790dc85072bb Mon Sep 17 00:00:00 2001 From: nojhan Date: Thu, 10 Nov 2011 15:56:59 +0100 Subject: [PATCH 15/24] first attempt at implementing a robust cholesky --- edo/src/edoSamplerNormalMulti.h | 32 ++++++++++++++++++++++++++------ 1 file changed, 26 insertions(+), 6 deletions(-) diff --git a/edo/src/edoSamplerNormalMulti.h b/edo/src/edoSamplerNormalMulti.h index 5f87152c..5e5910fb 100644 --- a/edo/src/edoSamplerNormalMulti.h +++ b/edo/src/edoSamplerNormalMulti.h @@ -66,7 +66,7 @@ public: //! use the algorithm using absolute value within the square root @see factorize_abs absolute, //! use the robust algorithm, without square root - //robust + robust }; Method _use; @@ -150,6 +150,8 @@ public: factorize_std( V ); } else if( _use == absolute ) { factorize_abs( V ); + } else if( _use == robust ) { + factorize_robust( V ); } } @@ -183,7 +185,6 @@ public: // round-off errors may appear here assert( V(i,i) - sum >= 0 ); _L(i,i) = sqrt( V(i,i) - sum ); - //_L(i,i) = sqrt( fabs( V(i,i) - sum) ); for ( j = i + 1; j < N; ++j ) { // rows // one element @@ -247,18 +248,37 @@ public: * * Computes L and D such as V = L D Lt */ - /* void factorize_robust( const MatrixType& V) { unsigned int N = assert_properties( V ); unsigned int i, j, k; - ublas::symmetric_matrix< AtomType, ublas::lower > D = ublas::zero_matrix(N); + //MatrixType D = ublas::zero_matrix(N); + std::vector _D(N,0); - _L(0, 0) = sqrt( V(0, 0) ); + _D[0] = V(0,0); + _L(0, 0) = 1; + //_L(1,0) = 1/D[0] * V(1,0); + for( j=0; j Date: Sat, 12 Nov 2011 14:24:54 +0100 Subject: [PATCH 16/24] more pertinent names for cholesky factorization functions --- edo/src/edoSamplerNormalMulti.h | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/edo/src/edoSamplerNormalMulti.h b/edo/src/edoSamplerNormalMulti.h index 5e5910fb..c16c77ef 100644 --- a/edo/src/edoSamplerNormalMulti.h +++ b/edo/src/edoSamplerNormalMulti.h @@ -61,11 +61,11 @@ public: typedef ublas::symmetric_matrix< AtomType, ublas::lower > MatrixType; enum Method { - //! use the standard algorithm, with square root @see factorize + //! use the standard algorithm, with square root @see factorize_LLT standard, - //! use the algorithm using absolute value within the square root @see factorize_abs + //! use the algorithm using absolute value within the square root @see factorize_LLT_abs absolute, - //! use the robust algorithm, without square root + //! use the robust algorithm, without square root @see factorize_LDLT robust }; Method _use; @@ -147,11 +147,11 @@ public: void factorize( const MatrixType& V ) { if( _use == standard ) { - factorize_std( V ); + factorize_LLT( V ); } else if( _use == absolute ) { - factorize_abs( V ); + factorize_LLT_abs( V ); } else if( _use == robust ) { - factorize_robust( V ); + factorize_LDLT( V ); } } @@ -162,7 +162,7 @@ public: * 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_std( const MatrixType& V) + void factorize_LLT( const MatrixType& V) { unsigned int N = assert_properties( V ); @@ -208,7 +208,7 @@ public: * Be aware that this increase round-off errors, this is just a ugly * hack to avoid crash. */ - void factorize_abs( const MatrixType & V) + void factorize_LLT_abs( const MatrixType & V) { unsigned int N = assert_properties( V ); @@ -248,7 +248,7 @@ public: * * Computes L and D such as V = L D Lt */ - void factorize_robust( const MatrixType& V) + void factorize_LDLT( const MatrixType& V) { unsigned int N = assert_properties( V ); From b2b1a964235ea9d15e3b547660ce0a38eb88ddda Mon Sep 17 00:00:00 2001 From: nojhan Date: Sat, 12 Nov 2011 15:52:18 +0100 Subject: [PATCH 17/24] working robust cholesky factorization, with test binary --- edo/src/edoSamplerNormalMulti.h | 48 +++++++++------- edo/test/CMakeLists.txt | 1 + edo/test/t-cholesky.cpp | 97 +++++++++++++++++++++++++++++++++ 3 files changed, 126 insertions(+), 20 deletions(-) create mode 100644 edo/test/t-cholesky.cpp diff --git a/edo/src/edoSamplerNormalMulti.h b/edo/src/edoSamplerNormalMulti.h index c16c77ef..997b29c6 100644 --- a/edo/src/edoSamplerNormalMulti.h +++ b/edo/src/edoSamplerNormalMulti.h @@ -84,7 +84,7 @@ public: */ Cholesky(const MatrixType& V, Cholesky::Method use = standard ) : _use(use) { - factsorize( V ); + factorize( V ); } @@ -99,10 +99,19 @@ public: //! The decomposition of the covariance matrix const MatrixType & decomposition() const {return _L;} + /** When your using the LDLT robust decomposition (by passing the "robust" + * option to the constructor, @see factorize_LDTL), this is the diagonal + * matrix part. + */ + const MatrixType & diagonal() const {return _D;} + protected: //! The decomposition is a (lower) symetric matrix, just like the covariance matrix MatrixType _L; + + //! The diagonal matrix when using the LDLT factorization + MatrixType _D; /** Assert that the covariance matrix have the required properties and returns its dimension. @@ -243,37 +252,36 @@ public: } // for i in [1,N[ } - /** This alternative algorithm does not use square root BUT the covariance - * matrix must be invertible. + + /** This alternative algorithm do not use square root. * * Computes L and D such as V = L D Lt */ void factorize_LDLT( const MatrixType& V) { - unsigned int N = assert_properties( 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 ); - unsigned int i, j, k; - //MatrixType D = ublas::zero_matrix(N); - std::vector _D(N,0); + _D = ublas::zero_matrix(N,N); + _D(0,0) = V(0,0); + + for( int j=0; j +*/ + +#include +#include +#include + +#include +#include +#include + +typedef eoReal< eoMinimizingFitness > EOT; +typedef edoNormalMulti EOD; + +std::ostream& operator<< (std::ostream& out, const ublas::symmetric_matrix< double, ublas::lower >& mat ) +{ + for( unsigned int i=0; i::Cholesky::MatrixType MatrixType; + + // a variance-covariance matrix of size N*N + MatrixType V(N,N); + + // random covariance matrix + for( unsigned int i=0; i 0 + for( unsigned int j=i+1; j::Cholesky LLT( edoSamplerNormalMulti::Cholesky::standard ); + edoSamplerNormalMulti::Cholesky LLTa( edoSamplerNormalMulti::Cholesky::absolute ); + edoSamplerNormalMulti::Cholesky LDLT( edoSamplerNormalMulti::Cholesky::robust ); + + MatrixType L0 = LLT(V); + std::cout << "LLT" << std::endl << L0 << std::endl; + MatrixType V0 = ublas::prod( L0, ublas::trans(L0) ); + std::cout << "LLT covar" << std::endl << V0 << std::endl; + std::cout << "-----------------------------------------------------------" << std::endl; + + MatrixType L1 = LLTa(V); + std::cout << "LLT abs" << std::endl << L1 << std::endl; + MatrixType V1 = ublas::prod( L1, ublas::trans(L1) ); + std::cout << "LLT covar" << std::endl << V1 << std::endl; + std::cout << "-----------------------------------------------------------" << std::endl; + + MatrixType L2 = LDLT(V); + MatrixType D2 = LDLT.diagonal(); + std::cout << "LDLT" << std::endl << L2 << std::endl; + // ublas do not allow nested products, we should use a temporary matrix, + // thus the inline instanciation of a MatrixType + // see: http://www.crystalclearsoftware.com/cgi-bin/boost_wiki/wiki.pl?Effective_UBLAS + MatrixType V2 = ublas::prod( MatrixType(ublas::prod( L2, D2 )), ublas::trans(L2) ); + std::cout << "LDLT covar" << std::endl << V2 << std::endl; + std::cout << "-----------------------------------------------------------" << std::endl; + +} From e0f691c148672406b1f5d58dfb392065d79aab44 Mon Sep 17 00:00:00 2001 From: nojhan Date: Sat, 12 Nov 2011 20:13:18 +0100 Subject: [PATCH 18/24] the complete robust cholesky factorization -- bugfix in LDLT: return the factorization LD^1/2 instead of L, thus no needs of accessors on D --- edo/src/edoSamplerNormalMulti.h | 54 +++++++++++++++++++-------------- edo/test/t-cholesky.cpp | 8 ++--- 2 files changed, 33 insertions(+), 29 deletions(-) diff --git a/edo/src/edoSamplerNormalMulti.h b/edo/src/edoSamplerNormalMulti.h index 997b29c6..e2852447 100644 --- a/edo/src/edoSamplerNormalMulti.h +++ b/edo/src/edoSamplerNormalMulti.h @@ -41,8 +41,8 @@ Authors: * - compute the Cholesky decomposition L of V (i.e. such as V=LL*) * - return X = M + LT */ -template< class EOT, typename D = edoNormalMulti< EOT > > -class edoSamplerNormalMulti : public edoSampler< D > +template< class EOT, typename EOD = edoNormalMulti< EOT > > +class edoSamplerNormalMulti : public edoSampler< EOD > { public: typedef typename EOT::AtomType AtomType; @@ -97,23 +97,16 @@ public: } //! The decomposition of the covariance matrix - const MatrixType & decomposition() const {return _L;} - - /** When your using the LDLT robust decomposition (by passing the "robust" - * option to the constructor, @see factorize_LDTL), this is the diagonal - * matrix part. - */ - const MatrixType & diagonal() const {return _D;} + const MatrixType & decomposition() const + { + return _L; + } protected: //! The decomposition is a (lower) symetric matrix, just like the covariance matrix MatrixType _L; - //! The diagonal matrix when using the LDLT factorization - MatrixType _D; - - /** 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. @@ -222,6 +215,7 @@ public: unsigned int N = assert_properties( V ); unsigned int i=0, j=0, k; + _L(0, 0) = sqrt( V(0, 0) ); // end of the column @@ -264,27 +258,41 @@ public: // example of an invertible matrix whose decomposition is undefined assert( V(0,0) != 0 ); - _D = ublas::zero_matrix(N,N); - _D(0,0) = V(0,0); + MatrixType L(N,N); + MatrixType 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< D >( repairer), _cholesky(use) + : edoSampler< EOD >( repairer), _cholesky(use) {} - EOT sample( edoNormalMulti< EOT >& distrib ) + EOT sample( EOD& distrib ) { unsigned int size = distrib.size(); assert(size > 0); diff --git a/edo/test/t-cholesky.cpp b/edo/test/t-cholesky.cpp index 11bcaa46..c06fd726 100644 --- a/edo/test/t-cholesky.cpp +++ b/edo/test/t-cholesky.cpp @@ -85,12 +85,8 @@ int main(int argc, char** argv) std::cout << "-----------------------------------------------------------" << std::endl; MatrixType L2 = LDLT(V); - MatrixType D2 = LDLT.diagonal(); - std::cout << "LDLT" << std::endl << L2 << std::endl; - // ublas do not allow nested products, we should use a temporary matrix, - // thus the inline instanciation of a MatrixType - // see: http://www.crystalclearsoftware.com/cgi-bin/boost_wiki/wiki.pl?Effective_UBLAS - MatrixType V2 = ublas::prod( MatrixType(ublas::prod( L2, D2 )), ublas::trans(L2) ); + std::cout << "LDLT: L" << std::endl << L2 << std::endl; + MatrixType V2 = ublas::prod( L2, ublas::trans(L2) ); std::cout << "LDLT covar" << std::endl << V2 << std::endl; std::cout << "-----------------------------------------------------------" << std::endl; From fe2cebc0e871dbe35d84e696754195f81ea9384e Mon Sep 17 00:00:00 2001 From: nojhan Date: Sat, 12 Nov 2011 23:44:31 +0100 Subject: [PATCH 19/24] BUGFIX: factorized matrix are not symetric, cholesky factorization should process different types for covariance and decomposition + better format output for cholesky test --- edo/src/edoSamplerNormalMulti.h | 58 ++++++++--------- edo/test/t-cholesky.cpp | 106 ++++++++++++++++++++++++-------- 2 files changed, 110 insertions(+), 54 deletions(-) diff --git a/edo/src/edoSamplerNormalMulti.h b/edo/src/edoSamplerNormalMulti.h index e2852447..35200c79 100644 --- a/edo/src/edoSamplerNormalMulti.h +++ b/edo/src/edoSamplerNormalMulti.h @@ -49,7 +49,7 @@ public: /** Cholesky decomposition, given a matrix V, return a matrix L - * such as V = L Lt (Lt being the conjugate transpose of 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. @@ -58,7 +58,8 @@ public: class Cholesky { public: - typedef ublas::symmetric_matrix< AtomType, ublas::lower > MatrixType; + typedef ublas::symmetric_matrix< AtomType, ublas::lower > CovarMat; + typedef ublas::matrix< AtomType > FactorMat; enum Method { //! use the standard algorithm, with square root @see factorize_LLT @@ -82,7 +83,8 @@ public: * * Use the standard unstable method by default. */ - Cholesky(const MatrixType& V, Cholesky::Method use = standard ) : _use(use) + Cholesky(const CovarMat& V, Cholesky::Method use = standard ) : + _use(use), _L(ublas::zero_matrix(V.size1(),V.size2())) { factorize( V ); } @@ -90,14 +92,14 @@ public: /** Compute the factorization and return the result */ - const MatrixType& operator()( const MatrixType& V ) + const FactorMat& operator()( const CovarMat& V ) { factorize( V ); return decomposition(); } //! The decomposition of the covariance matrix - const MatrixType & decomposition() const + const FactorMat & decomposition() const { return _L; } @@ -105,18 +107,18 @@ public: protected: //! The decomposition is a (lower) symetric matrix, just like the covariance matrix - MatrixType _L; + 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 MatrixType& V ) + unsigned assert_properties( const CovarMat& V ) { unsigned int Vl = V.size1(); // number of lines // the result goes in _L - _L.resize(Vl); + _L = ublas::zero_matrix(Vl,Vl); #ifndef NDEBUG assert(Vl > 0); @@ -135,7 +137,6 @@ public: * perform the cholesky factorization * check if all eigenvalues are positives * check if all of the leading principal minors are positive - * */ #endif @@ -146,7 +147,7 @@ public: /** Actually performs the factorization with the method given at * instanciation. Results are cached in _L. */ - void factorize( const MatrixType& V ) + void factorize( const CovarMat& V ) { if( _use == standard ) { factorize_LLT( V ); @@ -161,10 +162,12 @@ public: /** 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 MatrixType& V) + void factorize_LLT( const CovarMat& V) { unsigned int N = assert_properties( V ); @@ -210,7 +213,7 @@ public: * Be aware that this increase round-off errors, this is just a ugly * hack to avoid crash. */ - void factorize_LLT_abs( const MatrixType & V) + void factorize_LLT_abs( const CovarMat & V) { unsigned int N = assert_properties( V ); @@ -247,19 +250,21 @@ public: } - /** This alternative algorithm do not use square root. + /** 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 Lt + * 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 MatrixType& V) + 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 ); - MatrixType L(N,N); - MatrixType D = ublas::zero_matrix(N,N); + 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 0); - // Cholesky factorisation gererating matrix L from covariance - // matrix V. - // We must use cholesky.decomposition() to get the resulting matrix. - // // L = cholesky decomposition of varcovar - const typename Cholesky::MatrixType& L = _cholesky( distrib.varcovar() ); + const typename Cholesky::FactorMat& L = _cholesky( distrib.varcovar() ); - // T = vector of size elements drawn in N(0,1) rng.normal(1.0) + // 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(); diff --git a/edo/test/t-cholesky.cpp b/edo/test/t-cholesky.cpp index c06fd726..5716d848 100644 --- a/edo/test/t-cholesky.cpp +++ b/edo/test/t-cholesky.cpp @@ -28,6 +28,9 @@ Authors: #include #include #include +#include +#include +#include #include #include @@ -36,26 +39,73 @@ Authors: typedef eoReal< eoMinimizingFitness > EOT; typedef edoNormalMulti EOD; -std::ostream& operator<< (std::ostream& out, const ublas::symmetric_matrix< double, ublas::lower >& mat ) + +void setformat( std::ostream& out ) { + out << std::right; + out << std::setfill(' '); + out << std::setw( 5 + std::numeric_limits::digits10); + out << std::setprecision(std::numeric_limits::digits10); + out << std::setiosflags(std::ios_base::showpoint); +} + + +template +std::string format(const MT& mat ) +{ + std::ostringstream out; + setformat(out); + for( unsigned int i=0; i +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 +bool equal( const MT& M1, const MT& M2, double prec /* = 1/std::numeric_limits::digits10 ???*/ ) +{ + if( M1.size1() != M2.size1() || M1.size2() != M2.size2() ) { + return false; + } + + for( unsigned int i=0; i::Cholesky LLT( edoSamplerNormalMulti::Cholesky::standard ); + FactorMat L0 = LLT(V); + std::cout << "LLT" << std::endl << format(L0) << std::endl; + CovarMat V0 = ublas::prod( L0, ublas::trans(L0) ); + std::cout << "LLT covar" << std::endl << format(V0) << std::endl; + assert( equal(V0,V,precision) ); + std::cout << linesep << std::endl; + edoSamplerNormalMulti::Cholesky LLTa( edoSamplerNormalMulti::Cholesky::absolute ); - edoSamplerNormalMulti::Cholesky LDLT( edoSamplerNormalMulti::Cholesky::robust ); - - MatrixType L0 = LLT(V); - std::cout << "LLT" << std::endl << L0 << std::endl; - MatrixType V0 = ublas::prod( L0, ublas::trans(L0) ); - std::cout << "LLT covar" << std::endl << V0 << std::endl; - std::cout << "-----------------------------------------------------------" << std::endl; - - MatrixType L1 = LLTa(V); - std::cout << "LLT abs" << std::endl << L1 << std::endl; - MatrixType V1 = ublas::prod( L1, ublas::trans(L1) ); - std::cout << "LLT covar" << std::endl << V1 << std::endl; - std::cout << "-----------------------------------------------------------" << std::endl; + FactorMat L1 = LLTa(V); + std::cout << "LLT abs" << std::endl << format(L1) << std::endl; + CovarMat V1 = ublas::prod( L1, ublas::trans(L1) ); + std::cout << "LLT covar" << std::endl << format(V1) << std::endl; + assert( equal(V1,V,precision) ); + std::cout << linesep << std::endl; - MatrixType L2 = LDLT(V); - std::cout << "LDLT: L" << std::endl << L2 << std::endl; - MatrixType V2 = ublas::prod( L2, ublas::trans(L2) ); - std::cout << "LDLT covar" << std::endl << V2 << std::endl; - std::cout << "-----------------------------------------------------------" << std::endl; + edoSamplerNormalMulti::Cholesky LDLT( edoSamplerNormalMulti::Cholesky::robust ); + FactorMat L2 = LDLT(V); + std::cout << "LDLT" << std::endl << format(L2) << std::endl; + CovarMat V2 = ublas::prod( L2, ublas::trans(L2) ); + std::cout << "LDLT covar" << std::endl << format(V2) << std::endl; + assert( equal(V2,V,precision) ); + std::cout << linesep << std::endl; } From 17ce8a90e3d503a55874cd09d25dd4f2b03996ca Mon Sep 17 00:00:00 2001 From: nojhan Date: Sat, 12 Nov 2011 23:53:23 +0100 Subject: [PATCH 20/24] cholesky test can take matrix size and precision arguments --- edo/test/t-cholesky.cpp | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/edo/test/t-cholesky.cpp b/edo/test/t-cholesky.cpp index 5716d848..3b9ab5f8 100644 --- a/edo/test/t-cholesky.cpp +++ b/edo/test/t-cholesky.cpp @@ -99,7 +99,16 @@ bool equal( const MT& M1, const MT& M2, double prec /* = 1/std::numeric_limits= 2 ) { + N = std::atof(argv[1]); + } + if( argc >= 3 ) { + precision = std::atof(argv[2]); + } typedef edoSamplerNormalMulti::Cholesky::CovarMat CovarMat; typedef edoSamplerNormalMulti::Cholesky::FactorMat FactorMat; @@ -115,9 +124,13 @@ int main(int argc, char** argv) } } - double precision = 1e-15; - setformat(std::cout); + 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; From 06100a6b57d8d1526c4044083fba8eaf67f62f64 Mon Sep 17 00:00:00 2001 From: nojhan Date: Sun, 13 Nov 2011 21:40:39 +0100 Subject: [PATCH 21/24] test binary now computes average errors of the differents cholesky factorization methods over random covariance matrices --- edo/src/edoSamplerNormalMulti.h | 4 + edo/test/t-cholesky.cpp | 172 +++++++++++++++++++++++--------- 2 files changed, 131 insertions(+), 45 deletions(-) diff --git a/edo/src/edoSamplerNormalMulti.h b/edo/src/edoSamplerNormalMulti.h index 35200c79..e2ce880c 100644 --- a/edo/src/edoSamplerNormalMulti.h +++ b/edo/src/edoSamplerNormalMulti.h @@ -58,7 +58,11 @@ public: 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 { diff --git a/edo/test/t-cholesky.cpp b/edo/test/t-cholesky.cpp index 3b9ab5f8..06a93aee 100644 --- a/edo/test/t-cholesky.cpp +++ b/edo/test/t-cholesky.cpp @@ -25,12 +25,13 @@ Authors: Johann Dréo */ -#include +//#include #include #include #include #include #include +#include #include #include @@ -97,66 +98,147 @@ bool equal( const MT& M1, const MT& M2, double prec /* = 1/std::numeric_limits +MT error( const MT& M1, const MT& M2 ) +{ + assert( M1.size1() == M2.size1() && M1.size1() == M2.size2() ); + + MT Err = ublas::zero_matrix(M1.size1(),M1.size2()); + + for( unsigned int i=0; i +double trigsum( const MT& M ) +{ + double sum; + for( unsigned int i=0; i +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; - double precision = 1e-15; + unsigned int N = 4; // size of matrix + unsigned int R = 1000; // nb of repetitions if( argc >= 2 ) { - N = std::atof(argv[1]); + N = std::atoi(argv[1]); } if( argc >= 3 ) { - precision = std::atof(argv[2]); + 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::Cholesky::CovarMat CovarMat; typedef edoSamplerNormalMulti::Cholesky::FactorMat FactorMat; - // a variance-covariance matrix of size N*N - CovarMat V(N,N); + edoSamplerNormalMulti::Cholesky LLT( edoSamplerNormalMulti::Cholesky::standard ); + edoSamplerNormalMulti::Cholesky LLTa( edoSamplerNormalMulti::Cholesky::absolute ); + edoSamplerNormalMulti::Cholesky LDLT( edoSamplerNormalMulti::Cholesky::robust ); - // random covariance matrix - for( unsigned int i=0; i 0 - for( unsigned int j=i+1; j s0,s1,s2; + for( unsigned int n=0; n= 0 + for( unsigned int j=i+1; j::Cholesky LLT( edoSamplerNormalMulti::Cholesky::standard ); - FactorMat L0 = LLT(V); - std::cout << "LLT" << std::endl << format(L0) << std::endl; - CovarMat V0 = ublas::prod( L0, ublas::trans(L0) ); - std::cout << "LLT covar" << std::endl << format(V0) << std::endl; - assert( equal(V0,V,precision) ); - std::cout << linesep << std::endl; - - edoSamplerNormalMulti::Cholesky LLTa( edoSamplerNormalMulti::Cholesky::absolute ); - FactorMat L1 = LLTa(V); - std::cout << "LLT abs" << std::endl << format(L1) << std::endl; - CovarMat V1 = ublas::prod( L1, ublas::trans(L1) ); - std::cout << "LLT covar" << std::endl << format(V1) << std::endl; - assert( equal(V1,V,precision) ); - std::cout << linesep << std::endl; - - edoSamplerNormalMulti::Cholesky LDLT( edoSamplerNormalMulti::Cholesky::robust ); - FactorMat L2 = LDLT(V); - std::cout << "LDLT" << std::endl << format(L2) << std::endl; - CovarMat V2 = ublas::prod( L2, ublas::trans(L2) ); - std::cout << "LDLT covar" << std::endl << format(V2) << std::endl; - assert( equal(V2,V,precision) ); - std::cout << linesep << 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::Cholesky LLT( edoSamplerNormalMulti::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::Cholesky LLTa( edoSamplerNormalMulti::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::Cholesky LDLT( edoSamplerNormalMulti::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; } From 9decda0c6a5bd09019d396b1faac9660c7d65eee Mon Sep 17 00:00:00 2001 From: nojhan Date: Tue, 15 Nov 2011 17:10:46 +0100 Subject: [PATCH 22/24] cholesky factorization with rounding to zero --- edo/src/edoSamplerNormalMulti.h | 52 +++++++++++++++++++++++++++++++++ edo/test/t-cholesky.cpp | 12 ++++++-- 2 files changed, 61 insertions(+), 3 deletions(-) diff --git a/edo/src/edoSamplerNormalMulti.h b/edo/src/edoSamplerNormalMulti.h index e2ce880c..6099ed95 100644 --- a/edo/src/edoSamplerNormalMulti.h +++ b/edo/src/edoSamplerNormalMulti.h @@ -70,6 +70,8 @@ public: 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 }; @@ -157,6 +159,8 @@ public: 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 ); } @@ -254,6 +258,54 @@ public: } + /** 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. * diff --git a/edo/test/t-cholesky.cpp b/edo/test/t-cholesky.cpp index 06a93aee..b141e6e4 100644 --- a/edo/test/t-cholesky.cpp +++ b/edo/test/t-cholesky.cpp @@ -158,9 +158,10 @@ int main(int argc, char** argv) edoSamplerNormalMulti::Cholesky LLT( edoSamplerNormalMulti::Cholesky::standard ); edoSamplerNormalMulti::Cholesky LLTa( edoSamplerNormalMulti::Cholesky::absolute ); + edoSamplerNormalMulti::Cholesky LLTz( edoSamplerNormalMulti::Cholesky::zeroing ); edoSamplerNormalMulti::Cholesky LDLT( edoSamplerNormalMulti::Cholesky::robust ); - std::vector s0,s1,s2; + std::vector s0,s1,s2,s3; for( unsigned int n=0; n= 4 ) { From e09bb5551bd6160ba0fed79a150f2da13d8e958b Mon Sep 17 00:00:00 2001 From: nojhan Date: Thu, 17 Nov 2011 13:40:26 +0100 Subject: [PATCH 23/24] assert that the decimals parameters of the round repairer is <= 1.0 --- edo/src/edoRepairerRound.h | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/edo/src/edoRepairerRound.h b/edo/src/edoRepairerRound.h index c370a03c..a882021c 100644 --- a/edo/src/edoRepairerRound.h +++ b/edo/src/edoRepairerRound.h @@ -86,7 +86,11 @@ public: typedef typename EOT::AtomType ArgType; //! Generally speaking, we expect decimals being <= 1, but it can work for higher values - edoRepairerRoundDecimals( ArgType decimals ) : edoRepairerApplyBinary( edoRound, 1 / decimals ) {} + edoRepairerRoundDecimals( ArgType decimals ) : edoRepairerApplyBinary( edoRound, 1 / decimals ) + { + assert( decimals <= 1.0 ); + assert( 1/decimals >= 1.0 ); + } }; From 164a81a10caeba5ed1964cd817b1761108d91e11 Mon Sep 17 00:00:00 2001 From: nojhan Date: Thu, 17 Nov 2011 13:41:32 +0100 Subject: [PATCH 24/24] output pretty formated matrix in debug mode in the multi-normal sampler --- edo/src/edoSamplerNormalMulti.h | 23 ++++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/edo/src/edoSamplerNormalMulti.h b/edo/src/edoSamplerNormalMulti.h index 6099ed95..dc3519eb 100644 --- a/edo/src/edoSamplerNormalMulti.h +++ b/edo/src/edoSamplerNormalMulti.h @@ -29,6 +29,7 @@ Authors: #define _edoSamplerNormalMulti_h #include +#include #include #include @@ -110,6 +111,22 @@ public: 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); @@ -164,6 +183,8 @@ public: } else if( _use == robust ) { factorize_LDLT( V ); } + + eo::log << eo::debug << std::endl << "Decomposed matrix:" << format( _L ) << std::endl; } @@ -289,7 +310,7 @@ public: if( V(i,i) - sum >= 0 ) { _L(i,i) = sqrt( V(i,i) - sum); } else { - _L(i,i) = 0; + _L(i,i) = std::numeric_limits::epsilon(); } for ( j = i + 1; j < N; ++j ) { // rows