use rank mu selector ; bugfix estimator's linear algebra : mu is useless in estimator ; arx = pop^T ; store D as a diagonal ; cwise prod for covar recomposition ; more asserts

edo/application/cmaes/main.cpp

@@ -57,7 +57,7 @@ int main(int ac, char** av)
     eoState state;
 
     // Instantiate all needed parameters for EDA algorithm
-    double selection_rate = parser.createParam((double)0.5, "selection_rate", "Selection Rate", 'R', section).value(); // R
+    //double selection_rate = parser.createParam((double)0.5, "selection_rate", "Selection Rate", 'R', section).value(); // R
 
     unsigned long max_eval = parser.getORcreateParam((unsigned long)0, "maxEval", "Maximum number of evaluations (0 = none)", 'E', "Stopping criterion").value(); // E
 
@@ -69,10 +69,10 @@ int main(int ac, char** av)
 
     edoNormalAdaptive<RealVec> distribution(dim);
 
-    eoSelect< RealVec >* selector = new eoDetSelect< RealVec >( selection_rate );
+    eoSelect< RealVec >* selector = new eoRankMuSelect< RealVec >( mu );
     state.storeFunctor(selector);
 
-    edoEstimator< Distrib >* estimator = new edoEstimatorNormalAdaptive<RealVec>( distribution, mu );
+    edoEstimator< Distrib >* estimator = new edoEstimatorNormalAdaptive<RealVec>( distribution );
     state.storeFunctor(estimator);
 
     eoEvalFunc< RealVec >* plainEval = new Rosenbrock< RealVec >();

edo/src/edoEstimatorNormalAdaptive.h

@@ -45,12 +45,11 @@ class edoEstimatorNormalAdaptive : public edoEstimatorAdaptive< EOD >
 {
 public:
     typedef typename EOT::AtomType AtomType;
-    typedef typename EOD::Vector Vector;
+    typedef typename EOD::Vector Vector; // column vectors @see edoNormalAdaptive
     typedef typename EOD::Matrix Matrix;
 
-    edoEstimatorNormalAdaptive( EOD& distrib, unsigned int mu ) :
+    edoEstimatorNormalAdaptive( EOD& distrib ) :
         edoEstimatorAdaptive<EOD>( distrib ),
-        _mu(mu),
         _calls(0),
         _eigeneval(0)
     {}
@@ -61,14 +60,15 @@ private:
         // compute recombination weights
         Eigen::VectorXd weights( pop_size );
         double sum_w = 0;
-        for( unsigned int i = 0; i < _mu; ++i ) {
+        for( unsigned int i = 0; i < pop_size; ++i ) {
-            double w_i = log( _mu + 0.5 ) - log( i + 1 );
+            double w_i = log( pop_size + 0.5 ) - log( i + 1 );
             weights(i) = w_i;
             sum_w += w_i;
         }
         // normalization of weights
         weights /= sum_w;
 
+        assert( weights.size() == pop_size);
         return weights;
     }
 
@@ -97,17 +97,18 @@ public:
         // Here, if we are in canonical CMA-ES,
         // pop is supposed to be the mu ranked better solutions,
         // as the rank mu selection is supposed to have occured.
-        Matrix arx( pop.size(), N );
+        Matrix arx( N, lambda );
 
         // copy the pop (most probably a vector of vectors) in a Eigen3 matrix
-        for( unsigned int i = 0; i < lambda; ++i ) {
+        for( unsigned int d = 0; d < N; ++d ) {
-            for( unsigned int d = 0; d < N; ++d ) {
+            for( unsigned int i = 0; i < lambda; ++i ) {
-                arx(i,d) = pop[i][d];
+                arx(d,i) = pop[i][d]; // NOTE: pop = arx.transpose()
             } // dimensions
         } // individuals
 
         // muXone array for weighted recombination
-        Eigen::VectorXd weights = edoCMAESweights( N );
+        Eigen::VectorXd weights = edoCMAESweights( lambda );
+        assert( weights.size() == lambda );
 
         // FIXME exposer les constantes dans l'interface
 
@@ -137,6 +138,8 @@ public:
         Matrix invsqrtC =
             d.coord_sys() * d.scaling().asDiagonal().inverse()
             * d.coord_sys().transpose();
+        assert( invsqrtC.innerSize() == d.coord_sys().innerSize() );
+        assert( invsqrtC.outerSize() == d.coord_sys().outerSize() );
 
         // expectation of ||N(0,I)|| == norm(randn(N,1))
         double chiN = sqrt(N)*(1-1/(4*N)+1/(21*pow(N,2)));
@@ -148,8 +151,11 @@ public:
 
         // compute weighted mean into xmean
         Vector xold = d.mean();
-        d.mean( arx * weights );
+        assert( xold.size() == N );
-        Vector xmean = d.mean();
+        //  xmean ( N, 1 ) = arx( N, lambda ) * weights( lambda, 1 )
+        Vector xmean = arx * weights;
+        assert( xmean.size() == N );
+        d.mean( xmean );
 
 
         /**********************************************************************
@@ -157,9 +163,9 @@ public:
          *********************************************************************/
 
         // cumulation for sigma
-        d.path_sigma( 
+        d.path_sigma(
-                (1.0-cs)*d.path_sigma() + sqrt(cs*(2.0-cs)*mueff)*invsqrtC*(xmean-xold)/d.sigma()
+            (1.0-cs)*d.path_sigma() + sqrt(cs*(2.0-cs)*mueff)*invsqrtC*(xmean-xold)/d.sigma()
-                );
+        );
 
         // sign of h
         double hsig;
@@ -173,10 +179,16 @@ public:
 
         // cumulation for the covariance matrix
         d.path_covar(
-                (1.0-cc)*d.path_covar() + hsig*sqrt(cc*(2.0-cc)*mueff)*(xmean-xold) / d.sigma()
+            (1.0-cc)*d.path_covar() + hsig*sqrt(cc*(2.0-cc)*mueff)*(xmean-xold) / d.sigma()
-                );
+        );
 
-        Matrix artmp = (1.0/d.sigma()) * arx - xold.rowwise().replicate(_mu);
+        Matrix xmu( N, lambda);
+        xmu = xold.rowwise().replicate(lambda);
+        assert( xmu.innerSize() == N );
+        assert( xmu.outerSize() == lambda );
+        Matrix artmp = (1.0/d.sigma()) * (arx - xmu);
+        // Matrix artmp = (1.0/d.sigma()) * arx - xold.colwise().replicate(lambda);
+        assert( artmp.innerSize() == N && artmp.outerSize() == lambda );
 
 
         /**********************************************************************
@@ -214,10 +226,11 @@ public:
             Eigen::SelfAdjointEigenSolver<Matrix> eigensolver( d.covar() ); // FIXME use JacobiSVD?
             d.coord_sys( eigensolver.eigenvectors() );
             Matrix D = eigensolver.eigenvalues().asDiagonal();
+            assert( D.innerSize() == N && D.outerSize() == N );
 
             // from variance to standard deviations
             D.cwiseSqrt();
-            d.scaling( D );
+            d.scaling( D.diagonal() );
         }
 
         return d;
@@ -225,7 +238,6 @@ public:
 
 protected:
 
-    unsigned int _mu;
     unsigned int _calls;
     unsigned int _eigeneval;
 

edo/src/edoNormalAdaptive.h

@@ -41,10 +41,11 @@ class edoNormalAdaptive : public edoDistrib< EOT >
 public:
     //typedef EOT EOType;
     typedef typename EOT::AtomType AtomType;
-    typedef Eigen::Matrix< AtomType, Eigen::Dynamic, 1> Vector;
+    typedef Eigen::Matrix< AtomType, Eigen::Dynamic, 1> Vector; // column vectors ( n lines, 1 column)
     typedef Eigen::Matrix< AtomType, Eigen::Dynamic, Eigen::Dynamic> Matrix;
 
     edoNormalAdaptive( unsigned int dim ) :
+        _dim(dim),
         _mean( Vector::Zero(dim) ),
         _C( Matrix::Identity(dim,dim) ),
         _B( Matrix::Identity(dim,dim) ),
@@ -53,7 +54,7 @@ public:
         _p_c( Vector::Zero(dim) ),
         _p_s( Vector::Zero(dim) )
     {
-        assert( dim > 0);
+        assert( _dim > 0);
     }
 
     edoNormalAdaptive( unsigned int dim, 
@@ -88,23 +89,24 @@ public:
         return _mean.innerSize();
     }
 
-    Vector mean() const {return _mean;}
+    Vector mean()       const {return _mean;}
-    Matrix covar() const {return _C;}
+    Matrix covar()      const {return _C;}
-    Matrix coord_sys() const {return _B;}
+    Matrix coord_sys()  const {return _B;}
-    Vector scaling() const {return _D;}
+    Vector scaling()    const {return _D;}
-    double sigma() const {return _sigma;}
+    double sigma()      const {return _sigma;}
     Vector path_covar() const {return _p_c;}
     Vector path_sigma() const {return _p_s;}
 
-    void mean(       Vector m ) { _mean = m; }
+    void mean(       Vector m ) { _mean = m;  assert( m.size() == _dim ); }
-    void covar(      Matrix c ) { _C = c; }
+    void covar(      Matrix c ) { _C = c;     assert( c.innerSize() == _dim && c.outerSize() == _dim ); }
-    void coord_sys(  Matrix b ) { _B = b; }
+    void coord_sys(  Matrix b ) { _B = b;     assert( b.innerSize() == _dim && b.outerSize() == _dim ); }
-    void scaling(    Vector d ) { _D = d; }
+    void scaling(    Vector d ) { _D = d;     assert( d.size() == _dim ); }
-    void sigma(      double s ) { _sigma = s; }
+    void sigma(      double s ) { _sigma = s; assert( s != 0.0 );}
-    void path_covar( Vector p ) { _p_c = p; }
+    void path_covar( Vector p ) { _p_c = p;   assert( p.size() == _dim ); }
-    void path_sigma( Vector p ) { _p_s = p; }
+    void path_sigma( Vector p ) { _p_s = p;   assert( p.size() == _dim ); }
 
 private:
+    unsigned int _dim;
     Vector _mean; // 
     Matrix _C; // covariance matrix
     Matrix _B; // eigen vectors / coordinates system

edo/src/edoSamplerNormalAdaptive.h

@@ -58,24 +58,28 @@ public:
 
     EOT sample( EOD& distrib )
     {
-        unsigned int size = distrib.size();
+        unsigned int N = distrib.size();
-        assert(size > 0);
+        assert( N > 0);
 
         // T = vector of size elements drawn in N(0,1)
-        Vector T( size );
+        Vector T( N );
-        for ( unsigned int i = 0; i < size; ++i ) {
+        for ( unsigned int i = 0; i < N; ++i ) {
             T( i ) = rng.normal();
         }
-        assert(T.innerSize() == size);
+        assert(T.innerSize() == N );
         assert(T.outerSize() == 1);
 
-        Vector sol = distrib.mean() + distrib.sigma() * distrib.coord_sys() * (distrib.scaling().dot(T) );
+        // 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)
+        assert( sol.size() == N );
         /*Vector sol = distrib.mean() + distrib.sigma()
             * distrib.coord_sys().dot( distrib.scaling().dot( T ) );*/
 
         // copy in the EOT structure (more probably a vector)
-        EOT solution( size );
+        EOT solution( N );
-        for( unsigned int i = 0; i < size; i++ ) {
+        for( unsigned int i = 0; i < N; i++ ) {
             solution[i]= sol(i);
         }