diff --git a/edo/src/edo b/edo/src/edo index a6a0861b5..e6a03a1ed 100644 --- a/edo/src/edo +++ b/edo/src/edo @@ -79,7 +79,7 @@ Authors: #include "edoContinue.h" #include "edoCombinedContinue.h" -#include "edoContAdaptiveIllCond.h" +#include "edoContAdaptiveIllCovar.h" #include "edoContAdaptiveFinite.h" #include "utils/edoCheckPoint.h" diff --git a/edo/src/edoContAdaptiveFinite.h b/edo/src/edoContAdaptiveFinite.h index a01775759..c92e6d4ff 100644 --- a/edo/src/edoContAdaptiveFinite.h +++ b/edo/src/edoContAdaptiveFinite.h @@ -47,24 +47,43 @@ public: bool operator()(const D& d) { - // Try to finite_check in most probably ill-conditioned order. - return finite_check(d.covar()) - and finite_check(d.path_covar()) - and finite_check(d.coord_sys()) - and finite_check(d.scaling()) - and finite_check(d.path_sigma()) - and finite_check(d.sigma()) - ; + bool fin_sigma = is_finite(d.sigma() ); + bool fin_path_sigma = is_finite(d.path_sigma()); + bool fin_scaling = is_finite(d.scaling() ); + bool fin_coord_sys = is_finite(d.coord_sys() ); + bool fin_path_covar = is_finite(d.path_covar()); + bool fin_covar = is_finite(d.covar() ); + + bool all_finite = fin_covar + and fin_path_covar + and fin_coord_sys + and fin_scaling + and fin_path_sigma + and fin_sigma; + + if( not all_finite ) { + eo::log << eo::progress << "STOP because parameters are not finite: "; + if( not fin_covar ) { eo::log << eo::errors << "covar, "; } + if( not fin_path_covar ) { eo::log << eo::errors << "path_covar, "; } + if( not fin_coord_sys ) { eo::log << eo::errors << "coord_sys, "; } + if( not fin_scaling ) { eo::log << eo::errors << "scaling, "; } + if( not fin_path_sigma ) { eo::log << eo::errors << "path_sigma, "; } + if( not fin_sigma ) { eo::log << eo::errors << "sigma"; } + eo::log << eo::errors << std::endl; + } + return all_finite; } virtual std::string className() const { return "edoContAdaptiveFinite"; } protected: - bool finite_check(const Matrix& mat) const + bool is_finite(const Matrix& mat) const { for(long i=0; i -*/ - -#ifndef _edoContAdaptiveIllCond_h -#define _edoContAdaptiveIllCond_h - -#ifdef WITH_EIGEN - -#include - -#include "edoContinue.h" - -/** A continuator that check if any matrix among the parameters - * of an edoNormalAdaptive distribution are ill-conditioned. - * - * If the condition number of the covariance matrix - * or the coordinate system matrix are strictly greater - * than the threshold given at construction, it will ask for a stop. - * - * @ingroup Continuators - */ -template -class edoContAdaptiveIllCond : public edoContinue -{ -public: - using EOType = typename D::EOType; - using Matrix = typename D::Matrix; - using Vector = typename D::Vector; - - edoContAdaptiveIllCond( double threshold = 1e6) : - _threshold(threshold) - { } - - bool operator()(const D& d) - { - if( condition(d.covar()) > _threshold - or condition(d.coord_sys()) > _threshold ) { - return false; - } else { - return true; - } - } - - virtual std::string className() const { return "edoContAdaptiveIllCond"; } - -public: - // Public function in case someone would want to dimensionate the condition threshold. - //! Returns the condition number - bool condition(const Matrix& mat) const - { - Eigen::JacobiSVD svd(mat); - return svd.singularValues()(0) / svd.singularValues()(svd.singularValues().size()-1); - } - - const double _threshold; -}; - -#endif // WITH_EIGEN - -#endif diff --git a/edo/src/edoContAdaptiveIllCovar.h b/edo/src/edoContAdaptiveIllCovar.h new file mode 100644 index 000000000..409d6db33 --- /dev/null +++ b/edo/src/edoContAdaptiveIllCovar.h @@ -0,0 +1,112 @@ +/* +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) 2020 Thales group +*/ +/* +Authors: + Johann Dréo +*/ + +#ifndef _edoContAdaptiveIllCovar_h +#define _edoContAdaptiveIllCovar_h + +#ifdef WITH_EIGEN + +#include + +#include "edoContinue.h" + +/** A continuator that check if the covariance matrix + * of an edoNormalAdaptive distribution is ill-conditioned. + * + * If the condition number of the covariance matrix + * is strictly greater than the threshold given at construction, + * it will ask for a stop. + * + * @ingroup Continuators + */ +template +class edoContAdaptiveIllCovar : public edoContinue +{ +public: + using EOType = typename D::EOType; + using Matrix = typename D::Matrix; + using Vector = typename D::Vector; + + edoContAdaptiveIllCovar( double threshold = 1e6) : + _threshold(threshold) + { } + + bool operator()(const D& d) + { + Eigen::SelfAdjointEigenSolver eigensolver( d.covar() ); + + auto info = eigensolver.info(); + if(info == Eigen::ComputationInfo::NumericalIssue) { + eo::log << eo::warnings << "WARNING: the eigen decomposition of the covariance matrix" + << " did not satisfy the prerequisites." << std::endl; + } else if(info == Eigen::ComputationInfo::NoConvergence) { + eo::log << eo::warnings << "WARNING: the eigen decomposition of the covariance matrix" + << " did not converged." << std::endl; + } else if(info == Eigen::ComputationInfo::InvalidInput) { + eo::log << eo::warnings << "WARNING: the eigen decomposition of the covariance matrix" + << " had invalid inputs." << std::endl; + } + if(info != Eigen::ComputationInfo::Success) { + eo::log << eo::progress << "STOP because the covariance matrix" + << " cannot be decomposed" << std::endl; +#ifndef NDEBUG + eo::log << eo::xdebug + << "mean:\n" << d.mean() << std::endl + << "sigma:" << d.sigma() << std::endl + << "coord_sys:\n" << d.coord_sys() << std::endl + << "scaling:\n" << d.scaling() << std::endl; +#endif + return false; + + }else { + Matrix EV = eigensolver.eigenvalues(); + double condition = EV.maxCoeff() / EV.minCoeff(); + + if( not std::isfinite(condition) ) { + eo::log << eo::progress << "STOP because the covariance matrix" + << " condition is not finite." << std::endl; + return false; + + } else if( condition >= _threshold ) { + eo::log << eo::progress << "STOP because the covariance matrix" + << " is ill-conditionned (condition number: " << condition << ")" << std::endl; + return false; + + } else { + return true; + } + } + } + + virtual std::string className() const { return "edoContAdaptiveIllCovar"; } + +protected: + const double _threshold; +}; + +#endif // WITH_EIGEN + +#endif diff --git a/edo/src/edoEstimatorNormalAdaptive.h b/edo/src/edoEstimatorNormalAdaptive.h index f1adf723f..be8c88d15 100644 --- a/edo/src/edoEstimatorNormalAdaptive.h +++ b/edo/src/edoEstimatorNormalAdaptive.h @@ -227,14 +227,15 @@ public: // Matrix CS = C.triangularView() + C.triangularView().transpose(); d.covar( C ); - Eigen::SelfAdjointEigenSolver eigensolver( d.covar() ); // FIXME use JacobiSVD? - d.coord_sys( eigensolver.eigenvectors() ); + Eigen::SelfAdjointEigenSolver eigensolver( d.covar() ); + Matrix mD = eigensolver.eigenvalues().asDiagonal(); - assert( mD.innerSize() == N && mD.outerSize() == N ); // from variance to standard deviations mD.cwiseSqrt(); d.scaling( mD.diagonal() ); + + d.coord_sys( eigensolver.eigenvectors() ); } return d; diff --git a/edo/src/edoSamplerNormalAdaptive.h b/edo/src/edoSamplerNormalAdaptive.h index 051e14835..f1b33d10f 100644 --- a/edo/src/edoSamplerNormalAdaptive.h +++ b/edo/src/edoSamplerNormalAdaptive.h @@ -72,11 +72,32 @@ public: // mean(N,1) + sigma * B(N,N) * ( D(N,1) .* T(N,1) ) Vector sol = distrib.mean() + distrib.sigma() - * distrib.coord_sys() * (distrib.scaling().cwiseProduct(T) ); // C * T = B * (D .* T) + * distrib.coord_sys() + * (distrib.scaling().cwiseProduct(T) ); // C * T = B * (D .* T) assert( sol.size() == N ); /*Vector sol = distrib.mean() + distrib.sigma() * distrib.coord_sys().dot( distrib.scaling().dot( T ) );*/ +#ifndef NDEBUG + bool is_finite = true; + for(long i=0; i