diff --git a/edo/src/edoEstimatorNormalMulti.h b/edo/src/edoEstimatorNormalMulti.h index 48aa4e96..76d535da 100644 --- a/edo/src/edoEstimatorNormalMulti.h +++ b/edo/src/edoEstimatorNormalMulti.h @@ -34,8 +34,8 @@ Authors: #include "edoNormalMulti.h" #ifdef WITH_BOOST -//! edoEstimatorNormalMulti< EOT > +//! edoEstimatorNormalMulti< EOT > template < typename EOT > class edoEstimatorNormalMulti : public edoEstimator< edoNormalMulti< EOT > > { @@ -43,95 +43,95 @@ public: class CovMatrix { public: - typedef typename EOT::AtomType AtomType; + typedef typename EOT::AtomType AtomType; - CovMatrix( const eoPop< EOT >& pop ) - { - //------------------------------------------------------------- - // Some checks before starting to estimate covar - //------------------------------------------------------------- + CovMatrix( const eoPop< EOT >& pop ) + { + //------------------------------------------------------------- + // Some checks before starting to estimate covar + //------------------------------------------------------------- - unsigned int p_size = pop.size(); // population size - assert(p_size > 0); + unsigned int p_size = pop.size(); // population size + assert(p_size > 0); - unsigned int s_size = pop[0].size(); // solution size - assert(s_size > 0); + unsigned int s_size = pop[0].size(); // solution size + assert(s_size > 0); - //------------------------------------------------------------- + //------------------------------------------------------------- - //------------------------------------------------------------- - // Copy the population to an ublas matrix - //------------------------------------------------------------- + //------------------------------------------------------------- + // Copy the population to an ublas matrix + //------------------------------------------------------------- - ublas::matrix< AtomType > sample( p_size, s_size ); + ublas::matrix< AtomType > sample( p_size, s_size ); - for (unsigned int i = 0; i < p_size; ++i) - { - for (unsigned int j = 0; j < s_size; ++j) - { - sample(i, j) = pop[i][j]; - } - } + for (unsigned int i = 0; i < p_size; ++i) + { + for (unsigned int j = 0; j < s_size; ++j) + { + sample(i, j) = pop[i][j]; + } + } - //------------------------------------------------------------- + //------------------------------------------------------------- - _varcovar.resize(s_size); + _varcovar.resize(s_size); - //------------------------------------------------------------- - // variance-covariance matrix are symmetric (and semi-definite - // positive), thus a triangular storage is sufficient - // - // variance-covariance matrix computation : transpose(A) * A - //------------------------------------------------------------- + //------------------------------------------------------------- + // variance-covariance matrix are symmetric (and semi-definite + // positive), thus a triangular storage is sufficient + // + // variance-covariance matrix computation : transpose(A) * A + //------------------------------------------------------------- - ublas::symmetric_matrix< AtomType, ublas::lower > var = ublas::prod( ublas::trans( sample ), sample ); + ublas::symmetric_matrix< AtomType, ublas::lower > var = ublas::prod( ublas::trans( sample ), sample ); - // Be sure that the symmetric matrix got the good size + // Be sure that the symmetric matrix got the good size - assert(var.size1() == s_size); - assert(var.size2() == s_size); - assert(var.size1() == _varcovar.size1()); - assert(var.size2() == _varcovar.size2()); + assert(var.size1() == s_size); + assert(var.size2() == s_size); + assert(var.size1() == _varcovar.size1()); + assert(var.size2() == _varcovar.size2()); - //------------------------------------------------------------- + //------------------------------------------------------------- - // TODO: to remove the comment below + // TODO: to remove the comment below - // for (unsigned int i = 0; i < s_size; ++i) - // { - // // triangular LOWER matrix, thus j is not going further than i - // for (unsigned int j = 0; j <= i; ++j) - // { - // // we want a reducted covariance matrix - // _varcovar(i, j) = var(i, j) / p_size; - // } - // } + // for (unsigned int i = 0; i < s_size; ++i) + // { + // // triangular LOWER matrix, thus j is not going further than i + // for (unsigned int j = 0; j <= i; ++j) + // { + // // we want a reducted covariance matrix + // _varcovar(i, j) = var(i, j) / p_size; + // } + // } - _varcovar = var / p_size; + _varcovar = var / p_size; - _mean.resize(s_size); // FIXME: check if it is really used because of the assignation below + _mean.resize(s_size); // FIXME: check if it is really used because of the assignation below - // unit vector - ublas::scalar_vector< AtomType > u( p_size, 1 ); + // unit vector + ublas::scalar_vector< AtomType > u( p_size, 1 ); - // sum over columns - _mean = ublas::prod( ublas::trans( sample ), u ); + // sum over columns + _mean = ublas::prod( ublas::trans( sample ), u ); - // division by n - _mean /= p_size; - } + // division by n + _mean /= p_size; + } - const ublas::symmetric_matrix< AtomType, ublas::lower >& get_varcovar() const {return _varcovar;} + const ublas::symmetric_matrix< AtomType, ublas::lower >& get_varcovar() const {return _varcovar;} - const ublas::vector< AtomType >& get_mean() const {return _mean;} + const ublas::vector< AtomType >& get_mean() const {return _mean;} private: - ublas::symmetric_matrix< AtomType, ublas::lower > _varcovar; - ublas::vector< AtomType > _mean; + ublas::symmetric_matrix< AtomType, ublas::lower > _varcovar; + ublas::vector< AtomType > _mean; }; public: @@ -139,21 +139,109 @@ public: edoNormalMulti< EOT > operator()(eoPop& pop) { - unsigned int popsize = pop.size(); - assert(popsize > 0); + unsigned int popsize = pop.size(); + assert(popsize > 0); - unsigned int dimsize = pop[0].size(); - assert(dimsize > 0); + unsigned int dimsize = pop[0].size(); + assert(dimsize > 0); - CovMatrix cov( pop ); + CovMatrix cov( pop ); - return edoNormalMulti< EOT >( cov.get_mean(), cov.get_varcovar() ); + return edoNormalMulti< EOT >( cov.get_mean(), cov.get_varcovar() ); } }; #else #ifdef WITH_EIGEN +//! edoEstimatorNormalMulti< EOT > +template < typename EOT > +class edoEstimatorNormalMulti : public edoEstimator< edoNormalMulti< EOT > > +{ +public: + class CovMatrix + { + public: + typedef typename EOT::AtomType AtomType; + typedef Eigen::Matrix< AtomType, Eigen::Dynamic, 1> Vector; + typedef Eigen::Matrix< AtomType, Eigen::Dynamic, Eigen::Dynamic> Matrix; + + CovMatrix( const eoPop< EOT >& pop ) + { + // Some checks before starting to estimate covar + unsigned int p_size = pop.size(); // population size + assert(p_size > 0); + unsigned int s_size = pop[0].size(); // solution size + assert(s_size > 0); + + // Copy the population to an ublas matrix + //ublas::matrix< AtomType > sample( p_size, s_size ); + Matrix sample( p_size, s_size ); + + for (unsigned int i = 0; i < p_size; ++i) { + for (unsigned int j = 0; j < s_size; ++j) { + sample(i, j) = pop[i][j]; + } + } + + // _varcovar.resize(s_size); + + // variance-covariance matrix are symmetric, thus a triangular storage is sufficient + // variance-covariance matrix computation : transpose(A) * A + //ublas::symmetric_matrix< AtomType, ublas::lower > var = ublas::prod( ublas::trans( sample ), sample ); + Matrix var = sample.transpose() * sample; + + // Be sure that the symmetric matrix got the good size + assert(var.innerSize() == s_size); + assert(var.outerSize() == s_size); + assert(var.innerSize() == _varcovar.innerSize()); + assert(var.outerSize() == _varcovar.outerSize()); + + _varcovar = var / p_size; + + // _mean.resize(s_size); // FIXME: check if it is really used because of the assignation below + + // unit vector + // ublas::scalar_vector< AtomType > u( p_size, 1 ); + Vector u( p_size, 1); + + // sum over columns + // _mean = ublas::prod( ublas::trans( sample ), u ); + _mean = sample.transpose() * u; + + // division by n + _mean /= p_size; + } + + // const ublas::symmetric_matrix< AtomType, ublas::lower >& get_varcovar() const {return _varcovar;} + const Matrix& get_varcovar() const {return _varcovar;} + + // const ublas::vector< AtomType >& get_mean() const {return _mean;} + const Vector& get_mean() const {return _mean;} + + private: + // ublas::symmetric_matrix< AtomType, ublas::lower > _varcovar; + Matrix _varcovar; + // ublas::vector< AtomType > _mean; + Vector _mean; + }; + +public: + typedef typename EOT::AtomType AtomType; + + edoNormalMulti< EOT > operator()(eoPop& pop) + { + unsigned int popsize = pop.size(); + assert(popsize > 0); + + unsigned int dimsize = pop[0].size(); + assert(dimsize > 0); + + CovMatrix cov( pop ); + + return edoNormalMulti< EOT >( cov.get_mean(), cov.get_varcovar() ); + } +}; #endif // WITH_EIGEN #endif // WITH_BOOST diff --git a/edo/src/edoNormalMulti.h b/edo/src/edoNormalMulti.h index 5ddd6c8e..d5f76871 100644 --- a/edo/src/edoNormalMulti.h +++ b/edo/src/edoNormalMulti.h @@ -21,8 +21,8 @@ Copyright (C) 2010 Thales group */ /* Authors: - Johann Dreo - Caner Candan + Johann Dreo + Caner Candan */ #ifndef _edoNormalMulti_h @@ -39,7 +39,6 @@ Copyright (C) 2010 Thales group namespace ublas = boost::numeric::ublas; //! edoNormalMulti< EOT > - template < typename EOT > class edoNormalMulti : public edoDistrib< EOT > { @@ -51,18 +50,18 @@ public: const ublas::vector< AtomType >& mean, const ublas::symmetric_matrix< AtomType, ublas::lower >& varcovar ) - : _mean(mean), _varcovar(varcovar) + : _mean(mean), _varcovar(varcovar) { - assert(_mean.size() > 0); - assert(_mean.size() == _varcovar.size1()); - assert(_mean.size() == _varcovar.size2()); + assert(_mean.size() > 0); + assert(_mean.size() == _varcovar.size1()); + assert(_mean.size() == _varcovar.size2()); } unsigned int size() { - assert(_mean.size() == _varcovar.size1()); - assert(_mean.size() == _varcovar.size2()); - return _mean.size(); + assert(_mean.size() == _varcovar.size1()); + assert(_mean.size() == _varcovar.size2()); + return _mean.size(); } ublas::vector< AtomType > mean() const {return _mean;} @@ -77,6 +76,44 @@ private: #else #ifdef WITH_EIGEN +#include + +template < typename EOT > +class edoNormalMulti : public edoDistrib< EOT > +{ +public: + typedef typename EOT::AtomType AtomType; + typedef Eigen::Matrix< AtomType, Eigen::Dynamic, 1> Vector; // Note: by default, Eigen is column-major + typedef Eigen::Matrix< AtomType, Eigen::Dynamic, Eigen::Dynamic> Matrix; + + edoNormalMulti( + const Vector & mean, + const Matrix & varcovar + ) + : _mean(mean), _varcovar(varcovar) + { + assert(_mean.innerSize() > 0); + assert(_mean.innerSize() == _varcovar.innerSize()); + assert(_mean.innerSize() == _varcovar.outerSize()); + } + + unsigned int size() + { + assert(_mean.innerSize() == _varcovar.innerSize()); + assert(_mean.innerSize() == _varcovar.outerSize()); + return _mean.innerSize(); + } + + Vector mean() const {return _mean;} + Matrix varcovar() const {return _varcovar;} + +private: + Vector _mean; + Matrix _varcovar; +}; + + + #endif // WITH_EIGEN #endif // WITH_BOOST diff --git a/edo/src/edoNormalMultiCenter.h b/edo/src/edoNormalMultiCenter.h index 548a7d06..8dac8461 100644 --- a/edo/src/edoNormalMultiCenter.h +++ b/edo/src/edoNormalMultiCenter.h @@ -31,7 +31,6 @@ Authors: #include "edoModifierMass.h" #include "edoNormalMulti.h" - #ifdef WITH_BOOST //! edoNormalMultiCenter< EOT > @@ -44,15 +43,29 @@ public: void operator() ( edoNormalMulti< EOT >& distrib, EOT& mass ) { - ublas::vector< AtomType > mean( distrib.size() ); - std::copy( mass.begin(), mass.end(), mean.begin() ); - distrib.mean() = mean; + ublas::vector< AtomType > mean( distrib.size() ); + std::copy( mass.begin(), mass.end(), mean.begin() ); + distrib.mean() = mean; } }; #else #ifdef WITH_EIGEN +template < typename EOT > +class edoNormalMultiCenter : public edoModifierMass< edoNormalMulti< EOT > > +{ +public: + typedef typename EOT::AtomType AtomType; + + void operator() ( edoNormalMulti< EOT >& distrib, EOT& mass ) + { + assert( distrib.size() == mass.innerSize() ); + Eigen::Matrix< AtomType, Eigen::Dynamic, 1 > mean( mass ); + distrib.mean() = mean; + } +}; + #endif // WITH_EIGEN #endif // WITH_BOOST diff --git a/edo/src/edoSamplerNormalMulti.h b/edo/src/edoSamplerNormalMulti.h index 5a2bbdf8..5b067a14 100644 --- a/edo/src/edoSamplerNormalMulti.h +++ b/edo/src/edoSamplerNormalMulti.h @@ -33,12 +33,6 @@ Authors: #include -#ifdef WITH_BOOST - -#include -#include -#include - /** 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: @@ -46,6 +40,13 @@ Authors: * - compute the Cholesky decomposition L of V (i.e. such as V=LL*) * - return X = M + LT */ + +#ifdef WITH_BOOST + +#include +#include +#include + template< class EOT, typename EOD = edoNormalMulti< EOT > > class edoSamplerNormalMulti : public edoSampler< EOD > { @@ -90,6 +91,60 @@ protected: #else #ifdef WITH_EIGEN +template< class EOT, typename EOD = edoNormalMulti< EOT > > +class edoSamplerNormalMulti : public edoSampler< EOD > +{ +public: + typedef typename EOT::AtomType AtomType; + typedef Eigen::Matrix< AtomType, Eigen::Dynamic, 1> Vector; + typedef Eigen::Matrix< AtomType, Eigen::Dynamic, Eigen::Dynamic> Matrix; + + edoSamplerNormalMulti( edoRepairer & repairer ) + : edoSampler< EOD >( repairer) + {} + + + EOT sample( EOD& distrib ) + { + unsigned int size = distrib.size(); + assert(size > 0); + + // L = cholesky decomposition of varcovar + + // Computes L and D such as V = L D L^T + Eigen::LDLT cholesky( distrib.varcovar() ); + Matrix L0 = cholesky.matrixL(); + Eigen::Diagonal D = cholesky.vectorD(); + + // 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 + Eigen::Diagonal sqrtD = D.cwiseSqrt(); + + Matrix L = L0 * D; + + + // T = vector of size elements drawn in N(0,1) + Vector T( size ); + for ( unsigned int i = 0; i < size; ++i ) { + T( i ) = rng.normal(); + } + + // LT = L * T + Vector LT = L * T; + + // solution = means + LT + Vector mean = distrib.mean(); + Vector typed_solution = mean + LT; + EOT solution( size ); + for( unsigned int i = 0; i < mean.innerSize(); i++ ) { + solution.push_back( typed_solution(i) ); + } + + return solution; + } +}; #endif // WITH_EIGEN #endif // WITH_BOOST diff --git a/edo/src/utils/edoStatNormalMulti.h b/edo/src/utils/edoStatNormalMulti.h index 32d975ed..de917f6c 100644 --- a/edo/src/utils/edoStatNormalMulti.h +++ b/edo/src/utils/edoStatNormalMulti.h @@ -28,16 +28,24 @@ Authors: #ifndef _edoStatNormalMulti_h #define _edoStatNormalMulti_h -#include +#include #include "edoStat.h" #include "edoNormalMulti.h" - #ifdef WITH_BOOST -//! edoStatNormalMulti< EOT > +#include +#else +#ifdef WITH_EIGEN + + // include nothing + +#endif // WITH_EIGEN +#endif // WITH_BOOST + +//! edoStatNormalMulti< EOT > template < typename EOT > class edoStatNormalMulti : public edoDistribStat< edoNormalMulti< EOT > > { @@ -47,34 +55,28 @@ public: using edoDistribStat< edoNormalMulti< EOT > >::value; edoStatNormalMulti( std::string desc = "" ) - : edoDistribStat< edoNormalMulti< EOT > >( desc ) + : edoDistribStat< edoNormalMulti< EOT > >( desc ) {} void operator()( const edoNormalMulti< EOT >& distrib ) { - value() = "\n# ====== multi normal distribution dump =====\n"; + value() = "\n# ====== multi normal distribution dump =====\n"; - std::ostringstream os; + std::ostringstream os; - os << distrib.mean() << " " << distrib.varcovar() << std::endl; + os << distrib.mean() << " " << distrib.varcovar() << std::endl; - // ublas::vector< AtomType > mean = distrib.mean(); - // std::copy(mean.begin(), mean.end(), std::ostream_iterator< std::string >( os, " " )); + // ublas::vector< AtomType > mean = distrib.mean(); + // std::copy(mean.begin(), mean.end(), std::ostream_iterator< std::string >( os, " " )); - // ublas::symmetric_matrix< AtomType, ublas::lower > varcovar = distrib.varcovar(); - // std::copy(varcovar.begin(), varcovar.end(), std::ostream_iterator< std::string >( os, " " )); + // ublas::symmetric_matrix< AtomType, ublas::lower > varcovar = distrib.varcovar(); + // std::copy(varcovar.begin(), varcovar.end(), std::ostream_iterator< std::string >( os, " " )); - // os << std::endl; + // os << std::endl; - value() += os.str(); + value() += os.str(); } }; -#else -#ifdef WITH_EIGEN - -#endif // WITH_EIGEN -#endif // WITH_BOOST - #endif // !_edoStatNormalMulti_h