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

Merge branch 'master' of ssh://eodev.git.sourceforge.net/gitroot/eodev/eodev

Browse source
This commit is contained in:
Caner Candan 2012-01-19 02:35:18 +01:00
commit 55cbeb0ca1
8 changed files with 350 additions and 340 deletions
Download patch file Download diff file Expand all files Collapse all files

11
edo/src/edoRepairerDispatcher.h
View file

@ -120,7 +120,10 @@ public:
//! Add more indexes set and their corresponding repairer operator address to the list
void add( ICT idx, edoRepairer<EOT>* op )
{
assert( idx.size() > 0 );
//assert( idx.size() > 0 );
#ifndef NDEBUG
eo::log << eo::warnings << "A repairer is added to the dispatcher while having an empty index list, nothing will be repaired" << std::endl;
#endif
assert( op != NULL );
this->push_back( std::make_pair(idx, op) );
@ -131,6 +134,7 @@ public:
{
// std::cout << "in dispatcher, sol = " << sol << std::endl;
// iterate over { indexe, repairer }
// ipair is an iterator that points on a pair of <indexes,repairer>
for( typename edoRepairerDispatcher<EOT>::iterator ipair = this->begin(); ipair != this->end(); ++ipair ) {
@ -140,6 +144,8 @@ public:
EOT partsol;
// std::cout << "\tusing indexes = ";
//
// iterate over 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 ) {
@ -151,6 +157,9 @@ public:
// std::cout << std::endl;
// std::cout << "\tpartial sol = " << partsol << std::endl;
if( partsol.size() == 0 ) {
continue;
}
assert( partsol.size() > 0 );
// apply the repairer on the partial copy

341
edo/src/edoSamplerNormalMulti.h
View file

@ -32,6 +32,7 @@ Authors:
#include <limits>
#include <edoSampler.h>
#include <utils/edoCholesky.h>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
@ -48,340 +49,8 @@ class edoSamplerNormalMulti : public edoSampler< EOD >
public:
typedef typename EOT::AtomType AtomType;
/** Cholesky decomposition, given a matrix V, return a matrix 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.
* Thus, expect a (lower) triangular matrix.
*/
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 {
//! use the standard algorithm, with square root @see factorize_LLT
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
};
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 CovarMat& V, Cholesky::Method use = standard ) :
_use(use), _L(ublas::zero_matrix<AtomType>(V.size1(),V.size2()))
{
factorize( V );
}
/** Compute the factorization and return the result
*/
const FactorMat& operator()( const CovarMat& V )
{
factorize( V );
return decomposition();
}
//! The decomposition of the covariance matrix
const FactorMat & decomposition() const
{
return _L;
}
template<typename MT>
std::string format(const MT& mat )
{
std::ostringstream out;
for( unsigned int i=0; i<mat.size1(); ++i) {
out << std::endl;
for( unsigned int j=0; j<mat.size2(); ++j) {
out << mat(i,j) << "\t";
} // columns
} // rows
return out.str();
}
protected:
//! The decomposition is a (lower) symetric matrix, just like the covariance matrix
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 CovarMat& V )
{
unsigned int Vl = V.size1(); // number of lines
// the result goes in _L
_L = ublas::zero_matrix<AtomType>(Vl,Vl);
eo::log << eo::debug << std::endl << "Covariance matrix:" << format( V ) << std::endl;
#ifndef NDEBUG
assert(Vl > 0);
unsigned int Vc = V.size2(); // number of columns
assert(Vc > 0);
assert( Vl == Vc );
// 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 );
}
/* 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 CovarMat& V )
{
if( _use == standard ) {
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 );
}
eo::log << eo::debug << std::endl << "Decomposed matrix:" << format( _L ) << std::endl;
}
/** 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 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);
}
// round-off errors may appear here
assert( V(i,i) - sum >= 0 );
_L(i,i) = sqrt( 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 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_LLT_abs( 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);
}
_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 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) = std::numeric_limits<double>::epsilon();
}
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.
*
* 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 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 );
FactorMat L = ublas::zero_matrix<AtomType>(N,N);
FactorMat D = ublas::zero_matrix<AtomType>(N,N);
D(0,0) = V(0,0);
for( int j=0; j<N; ++j ) { // each columns
L(j, j) = 1;
D(j,j) = V(j,j);
for( int k=0; k<=j-1; ++k) { // sum
D(j,j) -= L(j,k) * L(j,k) * D(k,k);
}
for( int i=j+1; i<N; ++i ) { // remaining rows
L(i,j) = V(i,j);
for( int k=0; k<=j-1; ++k) { // sum
L(i,j) -= L(i,k)*L(j,k) * D(k,k);
}
L(i,j) /= D(j,j);
} // for i in rows
} // for j in columns
// now compute the final symetric matrix: _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
FactorMat D12 = D;
for(int i=0; i<N; ++i) {
D12(i,i) = sqrt(D(i,i));
}
// the factorization is thus _L*D^1/2
_L = ublas::prod( L, D12);
}
}; // class Cholesky
edoSamplerNormalMulti( edoRepairer<EOT> & repairer, typename Cholesky::Method use = Cholesky::absolute )
: edoSampler< EOD >( repairer), _cholesky(use)
edoSamplerNormalMulti( edoRepairer<EOT> & repairer )
: edoSampler< EOD >( repairer)
{}
@ -391,7 +60,7 @@ public:
assert(size > 0);
// L = cholesky decomposition of varcovar
const typename Cholesky::FactorMat& L = _cholesky( distrib.varcovar() );
const typename cholesky::CholeskyBase<AtomType>::FactorMat& L = _cholesky( distrib.varcovar() );
// T = vector of size elements drawn in N(0,1)
ublas::vector< AtomType > T( size );
@ -412,7 +81,7 @@ public:
}
protected:
Cholesky _cholesky;
cholesky::CholeskyLLT<AtomType> _cholesky;
};
#endif // !_edoSamplerNormalMulti_h

285
edo/src/utils/edoCholesky.h Normal file
View file

@ -0,0 +1,285 @@
/*
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>
*/
namespace cholesky {
/** Cholesky decomposition, given a matrix V, return a matrix 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.
* Thus, expect a (lower) triangular matrix.
*/
template< typename T >
class CholeskyBase
{
public:
//! The covariance-matrix is symetric
typedef ublas::symmetric_matrix< T, ublas::lower > CovarMat;
//! The factorization matrix is triangular
// FIXME check if triangular types behaviour is like having 0
typedef ublas::matrix< T > FactorMat;
/** Instanciate without computing anything, you are responsible of
* calling the algorithm and getting the result with operator()
* */
CholeskyBase( size_t s1 = 1, size_t s2 = 1 ) :
_L(ublas::zero_matrix<T>(s1,s2))
{}
/** Computation is made at instanciation and then cached in a member variable,
* use decomposition() to get the result.
*/
CholeskyBase(const CovarMat& V) :
_L(ublas::zero_matrix<T>(V.size1(),V.size2()))
{
(*this)( V );
}
/** Compute the factorization and cache the result */
virtual void factorize( const CovarMat& V ) = 0;
/** Compute the factorization and return the result */
virtual const FactorMat& operator()( const CovarMat& V )
{
this->factorize( V );
return decomposition();
}
//! The decomposition of the covariance matrix
const FactorMat & decomposition() const
{
return _L;
}
protected:
/** 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 CovarMat& V )
{
unsigned int Vl = V.size1(); // number of lines
// the result goes in _L
_L = ublas::zero_matrix<T>(Vl,Vl);
#ifndef NDEBUG
assert(Vl > 0);
unsigned int Vc = V.size2(); // number of columns
assert(Vc > 0);
assert( Vl == Vc );
// 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 );
}
/* 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;
}
//! The decomposition is a (lower) symetric matrix, just like the covariance matrix
FactorMat _L;
};
/** 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.
*/
template< typename T >
class CholeskyLLT : public CholeskyBase<T>
{
public:
virtual void factorize( const typename CholeskyBase<T>::CovarMat& V )
{
unsigned int N = assert_properties( V );
unsigned int i=0, j=0, k;
this->_L(0, 0) = sqrt( V(0, 0) );
// end of the column
for ( j = 1; j < N; ++j ) {
this->_L(j, 0) = V(0, j) / this->_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 += this->_L(i, k) * this->_L(i, k);
}
this->_L(i,i) = L_i_i( V, i, sum );
for ( j = i + 1; j < N; ++j ) { // rows
// one element
sum = 0.0;
for ( k = 0; k < i; ++k ) {
sum += this->_L(j, k) * this->_L(i, k);
}
this->_L(j, i) = (V(j, i) - sum) / this->_L(i, i);
} // for j in ]i,N[
} // for i in [1,N[
}
/** The step of the standard LLT algorithm where round off errors may appear */
inline virtual T L_i_i( const typename CholeskyBase<T>::CovarMat& V, const unsigned int& i, const double& sum ) const
{
// round-off errors may appear here
assert( V(i,i) - sum >= 0 );
return sqrt( V(i,i) - sum );
}
};
/** 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.
*/
template< typename T >
class CholeskyLLTabs : public CholeskyLLT<T>
{
public:
inline virtual T L_i_i( const typename CholeskyBase<T>::CovarMat& V, const unsigned int& i, const double& sum ) const
{
/***** ugly hack *****/
return sqrt( fabs( V(i,i) - sum) );
}
};
/** 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.
*/
template< typename T >
class CholeskyLLTzero : public CholeskyLLT<T>
{
public:
inline virtual T L_i_i( const typename CholeskyBase<T>::CovarMat& V, const unsigned int& i, const double& sum ) const
{
T Lii;
if( V(i,i) - sum >= 0 ) {
Lii = sqrt( V(i,i) - sum);
} else {
/***** ugly hack *****/
Lii = 0;
}
return Lii;
}
};
/** 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 L^T.
* The factorized matrix is (L D^1/2), because V = (L D^1/2) (L D^1/2)^T
*/
template< typename T >
class CholeskyLDLT : public CholeskyBase<T>
{
public:
virtual void factorize( const typename CholeskyBase<T>::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 );
typename CholeskyBase<T>::FactorMat L = ublas::zero_matrix<T>(N,N);
typename CholeskyBase<T>::FactorMat D = ublas::zero_matrix<T>(N,N);
D(0,0) = V(0,0);
for( int j=0; j<N; ++j ) { // each columns
L(j, j) = 1;
D(j,j) = V(j,j);
for( int k=0; k<=j-1; ++k) { // sum
D(j,j) -= L(j,k) * L(j,k) * D(k,k);
}
for( int i=j+1; i<N; ++i ) { // remaining rows
L(i,j) = V(i,j);
for( int k=0; k<=j-1; ++k) { // sum
L(i,j) -= L(i,k)*L(j,k) * D(k,k);
}
L(i,j) /= D(j,j);
} // for i in rows
} // for j in columns
this->_L = root( L, D );
}
inline typename CholeskyBase<T>::FactorMat root( typename CholeskyBase<T>::FactorMat& L, typename CholeskyBase<T>::FactorMat& D )
{
// 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
typename CholeskyBase<T>::FactorMat sqrt_D = D;
for( int i=0; i < D.size1(); ++i) {
sqrt_D(i,i) = sqrt(D(i,i));
}
// the factorization is thus this->_L*D^1/2
return ublas::prod( L, sqrt_D );
}
};
} // namespace cholesky

2
eo/CMakeLists.txt
View file

@ -16,7 +16,7 @@ INCLUDE(eo-conf.cmake OPTIONAL)
PROJECT(EO)
# CMake > 2.8 is needed, because of the FindOpenMP feature
cmake_minimum_required(VERSION 2.8)
#cmake_minimum_required(VERSION 2.8)
#SET(PROJECT_VERSION_MAJOR 1)
#SET(PROJECT_VERSION_MINOR 1)

2
eo/src/eoSIGContinue.h
View file

@ -48,7 +48,7 @@ void set_bool(int)
}
/**
Ctrl C handling: this eoContinue tells whether the user pressed Ctrl C
A continuator that stops if a given signal is received during the execution
*/
template< class EOT>
class eoSIGContinue: public eoContinue<EOT>

1
eo/src/utils/checkpointing
View file

@ -17,6 +17,7 @@
#include <utils/eoMOFitnessStat.h>
#include <utils/eoPopStat.h>
#include <utils/eoTimeCounter.h>
#include <utils/eoGenCounter.h>
// and make_help - any better suggestion to include it?
void make_help(eoParser & _parser);

46
eo/src/utils/eoGenCounter.h Normal file
View file

@ -0,0 +1,46 @@
/*
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 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
(c) Thales group 2011
Author: johann.dreo@thalesgroup.com
*/
#ifndef _eoGenCounter_h
#define _eoGenCounter_h
#include <string>
#include <utils/eoStat.h>
/**
An eoStat that simply gives the current generation index
@ingroup Stats
*/
class eoGenCounter : public eoUpdater, public eoValueParam<unsigned int>
{
public:
eoGenCounter( unsigned int start = 0, std::string label = "Gen" ) : eoValueParam<unsigned int>(start, label), _nb(start) {}
virtual void operator()()
{
value() = _nb++;
}
private:
unsigned int _nb;
};
#endif

2
eo/src/utils/eoStat.h
View file

@ -448,7 +448,7 @@ public:
std::nth_element( pop.begin(), pop.begin()+quartile*3, pop.end() );
typename EOT::Fitness Q3 = pop[quartile*3].fitness();
value() = Q1 - Q3;
value() = Q3 - Q1;
}
}