cleaner numerical errors management for EDO adaptive algos

- Change the ill-condition continuator to use eigen decomposition of the
covariance matrix, just like in the adaptive estimator.
- Add a warning message in adaptive sampler.

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"
 

edo/src/edoContAdaptiveFinite.h

@@ -47,24 +47,43 @@ public:
 
     bool operator()(const D& d)
     {
-        // Try to finite_check in most probably ill-conditioned order.
+        bool fin_sigma      = is_finite(d.sigma()     );
-        return finite_check(d.covar())
+        bool fin_path_sigma = is_finite(d.path_sigma());
-           and finite_check(d.path_covar())
+        bool fin_scaling    = is_finite(d.scaling()   );
-           and finite_check(d.coord_sys())
+        bool fin_coord_sys  = is_finite(d.coord_sys() );
-           and finite_check(d.scaling())
+        bool fin_path_covar = is_finite(d.path_covar());
-           and finite_check(d.path_sigma())
+        bool fin_covar      = is_finite(d.covar()     );
-           and finite_check(d.sigma())
+
-           ;
+        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<mat.rows(); ++i) {
             for(long j=0; j<mat.cols(); ++j) {
-                if(not finite_check(mat(i,j))) {
+                // Double negation because one want to escape
+                // as soon as one element is not finite.
+                if(not is_finite(mat(i,j))) {
                     return false;
                 }
             }
@@ -72,22 +91,22 @@ protected:
         return true;
     }
 
-    bool finite_check(const Vector& vec) const
+    bool is_finite(const Vector& vec) const
     {
         for(long i=0; i<vec.size(); ++i) {
-            if(not finite_check(vec[i])) {
+            if(not is_finite(vec[i])) {
                 return false;
             }
         }
         return true;
     }
 
-    bool finite_check(const typename EOType::AtomType& x) const
+    bool is_finite(const typename EOType::AtomType& x) const
     {
-        if(not std::isfinite(x)) {
+        if(std::isfinite(x)) {
-            return false;
-        } else {
             return true;
+        } else {
+            return false;
         }
     }
 };

edo/src/edoContAdaptiveIllCond.h

@@ -1,83 +0,0 @@
-/*
-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 <johann.dreo@thalesgroup.com>
-*/
-
-#ifndef _edoContAdaptiveIllCond_h
-#define _edoContAdaptiveIllCond_h
-
-#ifdef WITH_EIGEN
-
-#include<Eigen/Dense>
-
-#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 D>
-class edoContAdaptiveIllCond : public edoContinue<D>
-{
-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<Matrix> svd(mat);
-        return svd.singularValues()(0) / svd.singularValues()(svd.singularValues().size()-1);
-    }
-
-    const double _threshold;
-};
-
-#endif // WITH_EIGEN
-
-#endif

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 <johann.dreo@thalesgroup.com>
+*/
+
+#ifndef _edoContAdaptiveIllCovar_h
+#define _edoContAdaptiveIllCovar_h
+
+#ifdef WITH_EIGEN
+
+#include<Eigen/Dense>
+
+#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 D>
+class edoContAdaptiveIllCovar : public edoContinue<D>
+{
+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<Matrix> 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

edo/src/edoEstimatorNormalAdaptive.h

@@ -227,14 +227,15 @@ public:
             // Matrix CS = C.triangularView<Eigen::Upper>() + C.triangularView<Eigen::StrictlyUpper>().transpose();
             d.covar( C );
 
-            Eigen::SelfAdjointEigenSolver<Matrix> eigensolver( d.covar() ); // FIXME use JacobiSVD?
+            Eigen::SelfAdjointEigenSolver<Matrix> eigensolver( d.covar() );
-            d.coord_sys( eigensolver.eigenvectors() );
+
             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;

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<sol.size(); ++i) {
+                    if(not std::isfinite(sol(i))) {
+                        is_finite = false;
+                    }
+                }
+                if(not is_finite) {
+                    eo::log << eo::warnings << "WARNING: sampled solution is not finite"
+                       << " (the search should stop after this warning)" << std::endl;
+                    eo::log << eo::debug << sol << std::endl;
+                    eo::log << eo::xdebug
+                        << "mean:\n" << distrib.mean() << std::endl
+                        << "sigma:" << distrib.sigma() << std::endl
+                        << "coord_sys:\n" << distrib.coord_sys() << std::endl
+                        << "scaling:\n" << distrib.scaling() << std::endl;
+                }
+                // assert(is_finite);
+#endif
+
         // copy in the EOT structure (more probably a vector)
         EOT solution( N );
         for( unsigned int i = 0; i < N; i++ ) {