256 lines
7.8 KiB
C++
256 lines
7.8 KiB
C++
/*
|
|
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 <johann.dreo@thalesgroup.com>
|
|
Caner Candan <caner.candan@thalesgroup.com>
|
|
*/
|
|
|
|
|
|
#ifndef _edoEstimatorNormalMulti_h
|
|
#define _edoEstimatorNormalMulti_h
|
|
|
|
|
|
#include "edoEstimator.h"
|
|
#include "edoNormalMulti.h"
|
|
|
|
#ifdef WITH_BOOST
|
|
#include <boost/numeric/ublas/symmetric.hpp>
|
|
#include <boost/numeric/ublas/lu.hpp>
|
|
namespace ublas = boost::numeric::ublas;
|
|
#else
|
|
#ifdef WITH_EIGEN
|
|
#include <Eigen/Dense>
|
|
#endif // WITH_EIGEN
|
|
#endif // WITH_BOOST
|
|
|
|
|
|
/** An estimator for edoNormalMulti
|
|
*
|
|
* Exists in two implementations, using either
|
|
* <a href="http://www.boost.org/doc/libs/1_50_0/libs/numeric/ublas/doc/index.htm">Boost::uBLAS</a> (if compiled WITH_BOOST)
|
|
* or <a href="http://eigen.tuxfamily.org">Eigen3</a> (WITH_EIGEN).
|
|
*
|
|
* @ingroup Estimators
|
|
* @ingroup EMNA
|
|
* @ingroup Multinormal
|
|
*/
|
|
template < typename EOT, typename D=edoNormalMulti<EOT> >
|
|
class edoEstimatorNormalMulti : public edoEstimator<D>
|
|
{
|
|
|
|
#ifdef WITH_BOOST
|
|
public:
|
|
class CovMatrix
|
|
{
|
|
public:
|
|
typedef typename EOT::AtomType AtomType;
|
|
|
|
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 );
|
|
|
|
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 (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 );
|
|
|
|
// 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());
|
|
|
|
//-------------------------------------------------------------
|
|
|
|
|
|
// 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;
|
|
// }
|
|
// }
|
|
|
|
_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 );
|
|
|
|
// sum over columns
|
|
_mean = ublas::prod( ublas::trans( sample ), u );
|
|
|
|
// division by n
|
|
_mean /= p_size;
|
|
}
|
|
|
|
const ublas::symmetric_matrix< AtomType, ublas::lower >& get_varcovar() const {return _varcovar;}
|
|
|
|
const ublas::vector< AtomType >& get_mean() const {return _mean;}
|
|
|
|
private:
|
|
ublas::symmetric_matrix< AtomType, ublas::lower > _varcovar;
|
|
ublas::vector< AtomType > _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() );
|
|
}
|
|
};
|
|
|
|
#else
|
|
#ifdef WITH_EIGEN
|
|
|
|
public:
|
|
class CovMatrix
|
|
{
|
|
public:
|
|
typedef typename EOT::AtomType AtomType;
|
|
typedef typename D::Vector Vector;
|
|
typedef typename D::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
|
|
}; // class edoNormalMulti
|
|
|
|
#endif // !_edoEstimatorNormalMulti_h
|