add the Eigen library implementations of normal distributions computations

This commit is contained in:
nojhan 2012-07-09 18:47:35 +02:00
commit f3e1562a14
5 changed files with 301 additions and 106 deletions

View file

@ -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<EOT>& 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<EOT>& 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

View file

@ -21,8 +21,8 @@ Copyright (C) 2010 Thales group
*/
/*
Authors:
Johann Dreo <johann.dreo@thalesgroup.com>
Caner Candan <caner.candan@thalesgroup.com>
Johann Dreo <johann.dreo@thalesgroup.com>
Caner Candan <caner.candan@thalesgroup.com>
*/
#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 <Eigen/Dense>
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

View file

@ -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

View file

@ -33,12 +33,6 @@ Authors:
#include <edoSampler.h>
#ifdef WITH_BOOST
#include <utils/edoCholesky.h>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
/** 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 <utils/edoCholesky.h>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
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<EOT> & 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<Matrix> cholesky( distrib.varcovar() );
Matrix L0 = cholesky.matrixL();
Eigen::Diagonal<const Matrix> 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<const Matrix> 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

View file

@ -28,16 +28,24 @@ Authors:
#ifndef _edoStatNormalMulti_h
#define _edoStatNormalMulti_h
#include <boost/numeric/ublas/io.hpp>
#include<sstream>
#include "edoStat.h"
#include "edoNormalMulti.h"
#ifdef WITH_BOOST
//! edoStatNormalMulti< EOT >
#include <boost/numeric/ublas/io.hpp>
#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