nojhan/eodev
Archived
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Merge branch 'master' of ssh://localhost:8479/gitroot/eodev/eodev

Browse source
This commit is contained in:
nojhan 2012-07-13 15:35:02 +02:00
commit 09d08bf043
55 changed files with 1660 additions and 1064 deletions
Download patch file Download diff file Expand all files Collapse all files

222
edo/src/edoEstimatorNormalMulti.h
View file

@ -29,11 +29,13 @@ Authors:
#ifndef _edoEstimatorNormalMulti_h
#define _edoEstimatorNormalMulti_h
#include "edoEstimator.h"
#include "edoNormalMulti.h"
//! edoEstimatorNormalMulti< EOT >
#ifdef WITH_BOOST
//! edoEstimatorNormalMulti< EOT >
template < typename EOT >
class edoEstimatorNormalMulti : public edoEstimator< edoNormalMulti< EOT > >
{
@ -41,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:
@ -137,16 +139,102 @@ 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, typename EOD = edoNormalMulti<EOT> >
class edoEstimatorNormalMulti : public edoEstimator< EOD >
{
public:
class CovMatrix
{
public:
typedef typename EOT::AtomType AtomType;
typedef typename EOD::Vector Vector;
typedef typename EOD::Matrix 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
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];
}
}
// variance-covariance matrix are symmetric, thus a triangular storage is sufficient
// variance-covariance matrix computation : transpose(A) * A
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);
_varcovar = var / p_size;
// unit vector
Vector u( p_size);
u = Vector::Constant(p_size, 1);
// sum over columns
_mean = sample.transpose() * u;
// division by n
_mean /= p_size;
}
const Matrix& get_varcovar() const {return _varcovar;}
const Vector& get_mean() const {return _mean;}
private:
Matrix _varcovar;
Vector _mean;
};
public:
typedef typename EOT::AtomType AtomType;
edoNormalMulti< EOT > operator()(eoPop<EOT>& pop)
{
unsigned int p_size = pop.size();
assert(p_size > 0);
unsigned int s_size = pop[0].size();
assert(s_size > 0);
CovMatrix cov( pop );
assert( cov.get_mean().innerSize() == s_size );
assert( cov.get_mean().outerSize() == 1 );
assert( cov.get_varcovar().innerSize() == s_size );
assert( cov.get_varcovar().outerSize() == s_size );
return edoNormalMulti< EOT >( cov.get_mean(), cov.get_varcovar() );
}
};
#endif // WITH_EIGEN
#endif // WITH_BOOST
#endif // !_edoEstimatorNormalMulti_h

69
edo/src/edoNormalMulti.h
View file

@ -21,22 +21,24 @@ 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
#define _edoNormalMulti_h
#include "edoDistrib.h"
#ifdef WITH_BOOST
#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include "edoDistrib.h"
namespace ublas = boost::numeric::ublas;
//! edoNormalMulti< EOT >
template < typename EOT >
class edoNormalMulti : public edoDistrib< EOT >
{
@ -48,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;}
@ -70,4 +72,49 @@ private:
ublas::symmetric_matrix< AtomType, ublas::lower > _varcovar;
};
#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;
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
#endif // !_edoNormalMulti_h

33
edo/src/edoNormalMultiCenter.h
View file

@ -31,6 +31,8 @@ Authors:
#include "edoModifierMass.h"
#include "edoNormalMulti.h"
#ifdef WITH_BOOST
//! edoNormalMultiCenter< EOT >
template < typename EOT >
@ -41,10 +43,35 @@ 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, typename EOD = edoNormalMulti< EOT > >
class edoNormalMultiCenter : public edoModifierMass<EOD>
{
public:
typedef typename EOT::AtomType AtomType;
typedef typename EOD::Vector Vector;
void operator() ( edoNormalMulti< EOT >& distrib, EOT& mass )
{
assert( distrib.size() == mass.innerSize() );
Vector mean( distrib.size() );
for( unsigned int i=0; i < distrib.size(); i++ ) {
mean(i) = mass[i];
}
distrib.mean() = mean;
}
};
#endif // WITH_EIGEN
#endif // WITH_BOOST
#endif // !_edoNormalMultiCenter_h

92
edo/src/edoSamplerNormalMulti.h
View file

@ -32,9 +32,6 @@ Authors:
#include <limits>
#include <edoSampler.h>
#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.
*
@ -43,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 >
{
@ -84,4 +88,86 @@ protected:
cholesky::CholeskyLLT<AtomType> _cholesky;
};
#else
#ifdef WITH_EIGEN
template< class EOT, typename EOD = edoNormalMulti< EOT > >
class edoSamplerNormalMulti : public edoSampler< EOD >
{
public:
typedef typename EOT::AtomType AtomType;
typedef typename EOD::Vector Vector;
typedef typename EOD::Matrix Matrix;
edoSamplerNormalMulti( edoRepairer<EOT> & repairer )
: edoSampler< EOD >( repairer)
{}
EOT sample( EOD& distrib )
{
unsigned int size = distrib.size();
assert(size > 0);
// LsD = cholesky decomposition of varcovar
// Computes L and D such as V = L D L^T
Eigen::LDLT<Matrix> cholesky( distrib.varcovar() );
Matrix L = cholesky.matrixL();
assert(L.innerSize() == size);
assert(L.outerSize() == size);
Matrix D = cholesky.vectorD().asDiagonal();
assert(D.innerSize() == size);
assert(D.outerSize() == size);
// now compute the final symetric matrix: LsD = 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
Matrix sqrtD = D.cwiseSqrt();
assert(sqrtD.innerSize() == size);
assert(sqrtD.outerSize() == size);
Matrix LsD = L * sqrtD;
assert(LsD.innerSize() == size);
assert(LsD.outerSize() == size);
// 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();
}
assert(T.innerSize() == size);
assert(T.outerSize() == 1);
// LDT = (L D^1/2) * T
Vector LDT = LsD * T;
assert(LDT.innerSize() == size);
assert(LDT.outerSize() == 1);
// solution = means + LDT
Vector mean = distrib.mean();
assert(mean.innerSize() == size);
assert(mean.outerSize() == 1);
Vector typed_solution = mean + LDT;
assert(typed_solution.innerSize() == size);
assert(typed_solution.outerSize() == 1);
// copy in the EOT structure (more probably a vector)
EOT solution( size );
for( unsigned int i = 0; i < mean.innerSize(); i++ ) {
solution[i]= typed_solution(i);
}
assert( solution.size() == size );
return solution;
}
};
#endif // WITH_EIGEN
#endif // WITH_BOOST
#endif // !_edoSamplerNormalMulti_h

10
edo/src/utils/edoCholesky.h
View file

@ -27,6 +27,9 @@ Authors:
namespace cholesky {
#ifdef WITH_BOOST
/** Cholesky decomposition, given a matrix V, return a matrix L
* such as V = L L^T (L^T being the transposed of L).
*
@ -282,4 +285,11 @@ public:
}
};
#else
#ifdef WITH_EIGEN
#endif // WITH_EIGEN
#endif // WITH_BOOST
} // namespace cholesky

36
edo/src/utils/edoStatNormalMulti.h
View file

@ -28,13 +28,24 @@ Authors:
#ifndef _edoStatNormalMulti_h
#define _edoStatNormalMulti_h
#include <boost/numeric/ublas/io.hpp>
#include<sstream>
#include "edoStat.h"
#include "edoNormalMulti.h"
//! edoStatNormalMulti< EOT >
#ifdef WITH_BOOST
#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 > >
{
@ -44,27 +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() << std::endl << std::endl << 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();
}
};
#endif // !_edoStatNormalMulti_h