Refactor edoBinomialMulti to allow more complex data structures

Refactor distribution, sampler and estimator related to the multi-binomial distribution.
This introduce tomic methods which may be overloaded for data structures more complex than eoReal of vector of bool (the
default implentation).

edo/src/edoBinomialMulti.h

@@ -49,6 +49,13 @@ public:
     edoBinomialMulti( T initial_probas )
         : T(initial_probas) {}
 
+    /** Initialize all the probabilities to a constant
+     *
+     * 0.5 by default
+     */
+    edoBinomialMulti( unsigned int rows, unsigned int cols, double proba=0.5 )
+        : T::Constant(rows,cols,proba) {}
+
     /** Constructor without any assumption.
      */
     edoBinomialMulti() {}

edo/src/edoEstimatorBinomialMulti.h

@@ -41,55 +41,72 @@ Authors:
 template< class EOT, class D = edoBinomialMulti<EOT> >
 class edoEstimatorBinomialMulti : public edoEstimator<D>
 {
-    protected:
-        D eot2d( EOT from, unsigned int rows, unsigned int cols ) // FIXME maybe more elegant with Eigen::Map?
-        {
-            assert( rows > 0 );
-            assert( from.size() == rows );
-            assert( cols > 0 );
+protected:
+    /** Decide whether a given element of the distribution is true or false.
+     *
+     * The default implementation is to set the item to the value of the atom itself
+     * (which is a boolean in the basic version).
+     * If you have a more complex data structure, you can just overload this.
+     */
+    virtual void make( D & to, unsigned int i, unsigned int j, typename EOT::AtomType::const_iterator iatom )
+    {
+        to(i,j) = *iatom;
+    }
 
-            D to( Eigen::MatrixXd(rows, cols) );
-            for( unsigned int i=0; i < rows; ++i ) {
-                assert( from[i].size() == cols );
-                for( unsigned int j=0; j < cols; ++j ) {
-                    to(i,j) = from[i][j];
-                }
+    /** Transliterate a EOT in a boolean matrix
+     */
+    D eot2d( const EOT & from, unsigned int rows, unsigned int cols ) // FIXME maybe more elegant with Eigen::Map?
+    {
+        assert( rows > 0 );
+        assert( from.size() == rows );
+        assert( cols > 0 );
+
+        D to( Eigen::MatrixXd(rows, cols) );
+        unsigned int i=0;
+        for( typename EOT::const_iterator irow = from.begin(), end=from.end(); irow != end; ++irow ) {
+            assert( irow->size() == cols );
+            unsigned int j=0;
+            for( typename EOT::AtomType::const_iterator icol = irow->begin(), end=irow->end(); icol != end; ++icol ) {
+                make( to, i, j, icol );
+                j++;
             }
-
-            return to;
+            i++;
         }
 
-    public:
-        /** The expected EOT interface is the same as an Eigen3::MatrixXd.
-         */
-        D operator()( eoPop<EOT>& pop )
-        {
-            unsigned int popsize = pop.size();
-            assert(popsize > 0);
+        return to;
+    }
 
-            unsigned int rows = pop[0].size();
-            assert( rows > 0 );
-            unsigned int cols = pop[0][0].size();
-            assert( cols > 0 );
+public:
+    /** The expected EOT interface is the same as an Eigen3::MatrixXd.
+     */
+    D operator()( eoPop<EOT>& pop )
+    {
+        unsigned int popsize = pop.size();
+        assert(popsize > 0);
 
-            D probas( D::Zero(rows, cols) );
+        unsigned int rows = pop.begin()->size();
+        assert( rows > 0 );
+        unsigned int cols = pop.begin()->begin()->size();
+        assert( cols > 0 );
 
-            // We still need a loop over pop, because it is an eoVector
-            for (unsigned int i = 0; i < popsize; ++i) {
-                D indiv = eot2d( pop[i], rows, cols );
-                assert( indiv.rows() == rows && indiv.cols() == cols );
+        D probas( D::Zero(rows, cols) );
 
-                // the EOT matrix should be filled with 1 or 0 only
-                assert( indiv.sum() <= popsize );
+        // We still need a loop over pop, because it is an eoVector
+        for (unsigned int i = 0; i < popsize; ++i) {
+            D indiv = eot2d( pop[i], rows, cols );
+            assert( indiv.rows() == rows && indiv.cols() == cols );
 
-                probas += indiv / popsize;
+            // the EOT matrix should be filled with 1 or 0 only
+            assert( indiv.sum() <= popsize );
 
-                // sum and scalar product, no size pb expected
-                assert( probas.rows() == rows && probas.cols() == cols );
-            }
+            probas += indiv / popsize;
 
-            return probas;
+            // sum and scalar product, no size pb expected
+            assert( probas.rows() == rows && probas.cols() == cols );
         }
+
+        return probas;
+    }
 };
 
 #endif // WITH_EIGEN

edo/src/edoSamplerBinomialMulti.h

@@ -44,6 +44,24 @@ template< class EOT, class D = edoBinomialMulti<EOT> >
 class edoSamplerBinomialMulti : public edoSampler<D>
 {
 public:
+    typedef typename EOT::AtomType AtomType;
+
+    /** Called if the sampler draw the item at (i,j)
+     *
+     * The default implementation is to push back a true boolean.
+     * If you have a more complex data structure, you can just overload this.
+     */
+    virtual void make_true( AtomType & atom, unsigned int i, unsigned int j )
+    {
+        atom.push_back( 1 );
+    }
+
+    /** @see make_true */
+    virtual void make_false( AtomType & atom, unsigned int i, unsigned int j )
+    {
+        atom.push_back( 0 );
+    }
+
     EOT sample( D& distrib )
     {
         unsigned int rows = distrib.rows();
@@ -52,17 +70,22 @@ public:
         assert(cols > 0);
 
         // The point we want to draw
-        // x = {x1, x2, ..., xn}
+        // X = {x1, x2, ..., xn}
+        // with xn a container of booleans
         EOT solution;
 
         // Sampling all dimensions
         for( unsigned int i = 0; i < rows; ++i ) {
-            typename EOT::AtomType vec;
+            AtomType atom;
             for( unsigned int j = 0; j < cols; ++j ) {
                 // Toss a coin, biased by the proba of being 1.
-                vec.push_back( rng.flip( distrib(i,j) ) );
+                if( rng.flip( distrib(i,j) ) ) {
+                    make_true( atom, i, j );
+                } else {
+                    make_false( atom, i, j );
+                }
             }
-            solution.push_back( vec );
+            solution.push_back( atom );
         }
 
         return solution;