Cleaned up MR code
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Mon, 1 Mar 2010 20:08:11 +0000 (21:08 +0100)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Mon, 1 Mar 2010 20:08:11 +0000 (21:08 +0100)
ChangeLog
include/dai/graph.h
include/dai/mr.h
src/mr.cpp
tests/testfast.out

index 31a28b9..3d8135d 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,7 +1,10 @@
 git HEAD
 --------
 
-* Rewrote implementation of response propagation in MR
+* Cleaned up MR code:
+  - rewrote response propagation implementation with help of BBP
+  - now uses GraphAL to represent the Markov graph
+  - miscellaneous smaller cleanups
 * Fixed bug in BBPCostFunction::operator=() which prevented the desired 
   assignment from happening
 * [Stefano Pellegrini] Fixed bug in BP[logdomain=1,inference=MAXPROD]
index fe526f6..9556ca2 100644 (file)
@@ -235,6 +235,18 @@ class GraphAL {
             return sum;
         }
 
+        /// Returns the index of a given node \a n2 amongst the neighbors of \a n1
+        /** \note The time complexity is linear in the number of neighbors of \a n1
+         *  \throw OBJECT_NOT_FOUND if \a n2 is not a neighbor of \a n1
+         */
+        size_t findNb( size_t n1, size_t n2 ) {
+            for( size_t _n2 = 0; _n2 < nb(n1).size(); _n2++ )
+                if( nb( n1, _n2 ) == n2 )
+                    return _n2;
+            DAI_THROW(OBJECT_NOT_FOUND);
+            return nb(n1).size();
+        }
+
         /// Returns true if the graph is connected
         bool isConnected() const;
 
index 659df99..dff91d1 100644 (file)
@@ -5,7 +5,7 @@
  *  warranty. See the file COPYING for more details.
  *
  *  Copyright (C) 2007       Bastian Wemmenhove
- *  Copyright (C) 2007-2009  Joris Mooij  [joris dot mooij at libdai dot org]
+ *  Copyright (C) 2007-2010  Joris Mooij  [joris dot mooij at libdai dot org]
  *  Copyright (C) 2007       Radboud University Nijmegen, The Netherlands
  */
 
@@ -25,6 +25,7 @@
 #include <dai/enum.h>
 #include <dai/properties.h>
 #include <dai/exceptions.h>
+#include <dai/graph.h>
 #include <boost/dynamic_bitset.hpp>
 
 
@@ -33,36 +34,37 @@ namespace dai {
 
 /// Approximate inference algorithm by Montanari and Rizzo [\ref MoR05]
 /** \author Bastian Wemmenhove wrote the original implementation before it was merged into libDAI
- *  \todo Clean up code (use a BipartiteGraph-like implementation for the graph structure, and BBP for response propagation).
  */
 class MR : public DAIAlgFG {
     private:
         /// Is the underlying factor graph supported?
         bool supported;
-        /// con[i] = connectivity of spin \a i
-        std::vector<size_t>                             con;
-        /// nb[i] are the neighbours of spin \a i
-        std::vector<std::vector<size_t> >               nb;
-        /// tJ[i][_j] is the hyperbolic tangent of the interaction between spin \a i and its neighbour nb[i][_j]
+
+        /// The interaction graph (Markov graph)
+        GraphAL G;
+
+        /// Convenience shortcut
+        typedef GraphAL::Neighbor Neighbor;
+
+        /// tJ[i][_j] is the hyperbolic tangent of the interaction between spin \a i and its neighbour G.nb(i,_j)
         std::vector<std::vector<Real> >                 tJ;
         /// theta[i] is the local field on spin \a i
         std::vector<Real>                               theta;
+
         /// M[i][_j] is \f$ M^{(i)}_j \f$
         std::vector<std::vector<Real> >                 M;
-        /// The \a _j 'th neighbour of spin \a i has spin \a i as its kindex[i][_j]'th neighbour
-        std::vector<std::vector<size_t> >               kindex;
         /// Cavity correlations
         std::vector<std::vector<std::vector<Real> > >   cors;
-        /// Maximum connectivity
-        static const size_t kmax = 31;
+
         /// Type used for managing a subset of neighbors
         typedef boost::dynamic_bitset<> sub_nb;
-        /// Number of variables (spins)
-        size_t N;
+        
         /// Magnetizations
         std::vector<Real> Mag;
+        
         /// Maximum difference encountered so far
         Real _maxdiff;
+        
         /// Number of iterations needed
         size_t _iters;
 
@@ -81,9 +83,8 @@ class MR : public DAIAlgFG {
              *  - RESPPROP using response propagation ("linear response")
              *  - CLAMPING using clamping and BP
              *  - EXACT using JunctionTree
-             *  - RESPPROPOLD using response propagation ("linear response"), old implementation
              */
-            DAI_ENUM(InitType,RESPPROP,CLAMPING,EXACT,RESPPROPOLD);
+            DAI_ENUM(InitType,RESPPROP,CLAMPING,EXACT);
 
             /// Verbosity (amount of output sent to stderr)
             size_t verbose;
@@ -103,7 +104,7 @@ class MR : public DAIAlgFG {
 
     public:
         /// Default constructor
-        MR() : DAIAlgFG(), supported(), con(), nb(), tJ(), theta(), M(), kindex(), cors(), N(), Mag(), _maxdiff(), _iters(), props() {}
+        MR() : DAIAlgFG(), supported(), G(), tJ(), theta(), M(), cors(), Mag(), _maxdiff(), _iters(), props() {}
 
         /// Construct from FactorGraph \a fg and PropertySet \a opts
         /** \param opts Parameters @see Properties
@@ -133,31 +134,14 @@ class MR : public DAIAlgFG {
     //@}
 
     private:
-        /// Returns the signum of \a a
-        Real sign(Real a) { return (a >= 0) ? 1.0 : -1.0; }
-
-        /// Initialize N, con, nb, tJ, theta
-        void init(size_t Nin, Real *_w, Real *_th);
-        
-        /// Initialize kindex
-        void makekindex();
-        
         /// Initialize cors
-        Real init_cor();
+        Real calcCavityCorrelations();
         
-        /// Calculate cors using response propagation
-        Real init_cor_resp();
-        
-        /// Calculate cors using response propagation (old implementation)
-        /** \deprecated Should be removed soon
-         */
-        Real init_cor_resp_old();
-
         /// Iterate update equations for cavity fields
-        void solvemcav();
+        void propagateCavityFields();
         
         /// Calculate magnetizations
-        void solveM();
+        void calcMagnetizations();
 
         /// Calculate the product of all tJ[i][_j] for _j in A
         /** \param i variable index
index f1c6e48..0b051ac 100644 (file)
@@ -64,214 +64,8 @@ string MR::printProperties() const {
 }
 
 
-void MR::init(size_t Nin, Real *_w, Real *_th) {
-    size_t i,j;
-
-    N = Nin;
-
-    con.resize(N);
-    nb.resize(N);
-    tJ.resize(N);
-    for(i=0; i<N; i++ ) {
-        nb[i].resize(kmax);
-        tJ[i].resize(kmax);
-        con[i]=0;
-        for(j=0; j<N; j++ )
-            if( _w[i*N+j] != 0.0 ) {
-                nb[i][con[i]] = j;
-                tJ[i][con[i]] = tanh(_w[i*N+j]);
-                con[i]++;
-            }
-    }
-
-    theta.resize(N);
-    for(i=0; i<N; i++)
-      theta[i] = _th[i];
-}
-
-
-Real MR::init_cor_resp_old() {
-    Real md, maxdev;
-
-    size_t maxIters = 3000;
-    Real damping = 0.1;
-
-    vector<vector<Real> > tJ_org;
-    vector<vector<size_t> > nb_org;
-    vector<size_t> con_org;
-
-    // save original tJ, nb
-    nb_org = nb;
-    tJ_org = tJ;
-    con_org = con;
-
-    maxdev = 0.0;
-    for( size_t cavity = 0; cavity < N; cavity++ ) {    // for each spin to be removed
-        con = con_org;
-        size_t concav = con[cavity];
-
-        nb = nb_org;
-        tJ = tJ_org;
-
-        // Adapt the graph variables nb[], tJ[] and con[]
-        for(size_t i=0; i<con[cavity]; i++) {
-            size_t ij = nb[cavity][i];
-            size_t flag=0;
-            size_t j=0;
-            do{
-                if(nb[ij][j]==cavity){
-                    while(j<(con[ij]-1)){
-                        nb[ij][j]=nb[ij][j+1];
-                        tJ[ij][j] = tJ[ij][j+1];
-                        j++;
-                    }
-                flag=1;
-                }
-                j++;
-            } while(flag==0);
-        }
-        for(size_t i=0; i<con[cavity]; i++)
-            con[nb[cavity][i]]--;
-        con[cavity] = 0;
-
-        // Do everything starting from the new graph********
-
-        makekindex();
-
-        vector<Real> xfield(N*kmax,0.0);
-        vector<Real> rfield(N*kmax,0.0);
-
-        for( size_t i = 0; i < kmax*N; i++ )
-            xfield[i] = 3.0 * (2.0 * rnd_uniform() - 1.);
-
-        for( size_t _i2 = 0; _i2 < concav; _i2++ ) { // Subsequently apply a field to each cavity spin ****************
-            size_t i2 = nb[cavity][_i2];
-
-            for( size_t _i = 0; _i < con[i2]; _i++ )
-                rfield[kmax*i2 + _i] = 1.0;
-
-            size_t iters = 0;
-            do { // From here start the response and belief propagation
-                iters++;
-                md = 0.0;
-                for( size_t k = 0; k < N; k++ ){
-                    if( k != cavity ) {
-                        for( size_t _l = 0; _l < con[k]; _l++ ){
-                            Real xinter = theta[k];
-                            Real rinter = 0.;
-                            if( k == i2 ) 
-                                rinter += 1.;
-                            for( size_t _j = 0; _j < con[k]; _j++ )
-                                if( _j != _l ) {
-                                    size_t ind = kmax*nb[k][_j] + kindex[k][_j];  // index of cavity field "j \ k"
-                                    Real variab2 = tanh( xfield[ind] );
-                                    Real variab1 = tJ[k][_j] * variab2;
-                                    xinter += atanh( variab1 );
-                                    rinter += tJ[k][_j] * rfield[ind] * (1.0 - variab2*variab2) / (1.0 - variab1*variab1);
-                                }
-
-                            size_t ind = kmax * k + _l;  // index of cavity field "k \ l"
-
-                            // update xfield
-                            Real devs = xinter - xfield[ind];
-                            xfield[ind] += devs * damping;
-                            if( fabs(devs) > md )
-                                md = fabs( devs );
-
-                            // update rfield
-                            Real devs2 = rinter - rfield[ind];
-                            rfield[ind] += devs2 * damping;
-                            if( fabs(devs2) > md )
-                                md = fabs(devs2);
-                        }
-                    }
-                }
-            } while( (md > props.tol) && (iters < maxIters) ); // Precision condition reached -> BP and RP finished
-            if( iters == maxIters )
-                if( props.verbose >= 2 )
-                    cerr << "init_cor_resp_old: Convergence not reached (md=" << md << ")..." << endl;
-
-            if( md > maxdev )
-                maxdev = md;
-
-            // compute the observables (i.e. magnetizations and responses)******
-
-            for( size_t _k = 0; _k < concav; _k++ ) {
-                size_t k = nb[cavity][_k];
-                Real rinter = 0.;
-                Real xinter = theta[k];
-                if( _k != _i2 )
-                    for( size_t _j = 0; _j < con[k]; _j++ ) {
-                        size_t ind = kmax*nb[k][_j] + kindex[k][_j];  // index of cavity field "j \ k"
-                        Real variab2 = tanh( xfield[ind] );
-                        Real variab1 = tJ[k][_j] * variab2;
-                        xinter += atanh( variab1 );
-                        rinter += tJ[k][_j] * rfield[ind] * (1.0 - variab2*variab2) / (1.0 - variab1*variab1);
-                    }
-
-                if( k != i2 )
-                    cors[cavity][_i2][_k] = rinter * (1.0 - tanh(xinter) * tanh(xinter));
-                else
-                    cors[cavity][_i2][_k] = 0;
-            }
-        } // close for _i2 = 0...concav
-    }
-
-    // restore nb, tJ, con
-    tJ = tJ_org;
-    nb = nb_org;
-    con = con_org;
-
-    return maxdev;
-}
-
-
-Real MR::init_cor_resp() {
-    for( size_t i = 0; i < nrVars(); i++ ) {
-        vector<Factor> pairq;
-        BP bpcav(*this, PropertySet()("updates", string("SEQMAX"))("tol", (Real)1.0e-9)("maxiter", (size_t)10000)("verbose", (size_t)0)("logdomain", false));
-        bpcav.makeCavity( i );
-        bpcav.init();
-        bpcav.run();
-
-        BBP bbp( &bpcav, PropertySet()("verbose",(size_t)0)("tol",(Real)1.0e-9)("maxiter",(size_t)10000)("damping",(Real)0.0)("updates",string("SEQ_MAX")) );
-        for( size_t _j = 0; _j < con[i]; _j++ ) {
-            size_t j = nb[i][_j];
-
-            // Create weights for magnetization of some spin
-            Prob p( 2, 0.0 );
-            p[0] = -1.0;
-            p[1] = 1.0;
-
-            // BBP cost function would be the magnetization of spin j
-            vector<Prob> b1_adj;
-            b1_adj.reserve( nrVars() );
-            for( size_t l = 0; l < nrVars(); l++ )
-                if( l == j )
-                    b1_adj.push_back( p );
-                else
-                    b1_adj.push_back( Prob( 2, 0.0 ) );
-            bbp.init_V( b1_adj );
-
-            // run BBP to estimate adjoints
-            bbp.run();
-
-            for( size_t _k = 0; _k < con[i]; _k++ ) {
-                size_t k = nb[i][_k];
-                if( k != j )
-                    cors[i][_j][_k] = (bbp.adj_psi_V(k)[1] - bbp.adj_psi_V(k)[0]);
-                else
-                    cors[i][_j][_k] = 0.0;
-            }
-        }
-    }
-
-    return 0.0;
-}
-
-
 Real MR::T(size_t i, sub_nb A) {
-    sub_nb _nbi_min_A(con[i]);
+    sub_nb _nbi_min_A(G.nb(i).size());
     _nbi_min_A.set();
     _nbi_min_A &= ~A;
 
@@ -284,14 +78,14 @@ Real MR::T(size_t i, sub_nb A) {
 
 
 Real MR::T(size_t i, size_t _j) {
-    sub_nb j(con[i]);
+    sub_nb j(G.nb(i).size());
     j.set(_j);
     return T(i,j);
 }
 
 
 Real MR::Omega(size_t i, size_t _j, size_t _l) {
-    sub_nb jl(con[i]);
+    sub_nb jl(G.nb(i).size());
     jl.set(_j);
     jl.set(_l);
     Real Tijl = T(i,jl);
@@ -300,7 +94,7 @@ Real MR::Omega(size_t i, size_t _j, size_t _l) {
 
 
 Real MR::Gamma(size_t i, size_t _j, size_t _l1, size_t _l2) {
-    sub_nb jll(con[i]);
+    sub_nb jll(G.nb(i).size());
     jll.set(_j);
     Real Tij = T(i,jll);
     jll.set(_l1);
@@ -312,7 +106,7 @@ Real MR::Gamma(size_t i, size_t _j, size_t _l1, size_t _l2) {
 
 
 Real MR::Gamma(size_t i, size_t _l1, size_t _l2) {
-    sub_nb ll(con[i]);
+    sub_nb ll(G.nb(i).size());
     Real Ti = T(i,ll);
     ll.set(_l1);
     ll.set(_l2);
@@ -376,35 +170,34 @@ void MR::sum_subs(size_t j, sub_nb A, Real *sum_even, Real *sum_odd) {
 }
 
 
-void MR::solvemcav() {
+void MR::propagateCavityFields() {
     Real sum_even, sum_odd;
     Real maxdev;
     size_t maxruns = 1000;
 
-    makekindex();
-    for(size_t i=0; i<N; i++)
-        for(size_t _j=0; _j<con[i]; _j++)
-            M[i][_j]=0.1;
+    for( size_t i = 0; i < G.nrNodes(); i++ )
+        foreach( const Neighbor &j, G.nb(i) )
+            M[i][j.iter] = 0.1;
 
     size_t run=0;
     do {
         maxdev=0.0;
         run++;
-        for(size_t i=0; i<N; i++){ // for all i
-            for(size_t _j=0; _j<con[i]; _j++){ // for all j in N_i
-                size_t _i = kindex[i][_j];
-                size_t j = nb[i][_j];
-                DAI_ASSERT( nb[j][_i] == i );
+        for( size_t i = 0; i < G.nrNodes(); i++ ) {
+            foreach( const Neighbor &j, G.nb(i) ) {
+                size_t _j = j.iter;
+                size_t _i = G.findNb(j,i);
+                DAI_ASSERT( G.nb(j,_i) == i );
 
                 Real newM = 0.0;
                 if( props.updates == Properties::UpdateType::FULL ) {
-                    // find indices in nb[j] that do not correspond with i
-                    sub_nb _nbj_min_i(con[j]);
+                    // find indices in nb(j) that do not correspond with i
+                    sub_nb _nbj_min_i(G.nb(j).size());
                     _nbj_min_i.set();
-                    _nbj_min_i.reset(kindex[i][_j]);
+                    _nbj_min_i.reset(_i);
 
-                    // find indices in nb[i] that do not correspond with j
-                    sub_nb _nbi_min_j(con[i]);
+                    // find indices in nb(i) that do not correspond with j
+                    sub_nb _nbi_min_j(G.nb(i).size());
                     _nbi_min_j.set();
                     _nbi_min_j.reset(_j);
 
@@ -414,7 +207,7 @@ void MR::solvemcav() {
                     sum_subs(i, _nbi_min_j, &sum_even, &sum_odd);
                     Real denom = sum_even + tanh(theta[i]) * sum_odd;
                     Real numer = 0.0;
-                    for(size_t _k=0; _k<con[i]; _k++) if(_k != _j) {
+                    for(size_t _k=0; _k < G.nb(i).size(); _k++) if(_k != _j) {
                         sub_nb _nbi_min_jk(_nbi_min_j);
                         _nbi_min_jk.reset(_k);
                         sum_subs(i, _nbi_min_jk, &sum_even, &sum_odd);
@@ -423,21 +216,21 @@ void MR::solvemcav() {
                     newM -= numer / denom;
                 } else if( props.updates == Properties::UpdateType::LINEAR ) {
                     newM = T(j,_i);
-                    for(size_t _l=0; _l<con[i]; _l++) if( _l != _j )
+                    for(size_t _l=0; _l<G.nb(i).size(); _l++) if( _l != _j )
                         newM -= Omega(i,_j,_l) * tJ[i][_l] * cors[i][_j][_l];
-                    for(size_t _l1=0; _l1<con[j]; _l1++) if( _l1 != _i )
-                        for( size_t _l2=_l1+1; _l2<con[j]; _l2++) if( _l2 != _i)
+                    for(size_t _l1=0; _l1<G.nb(j).size(); _l1++) if( _l1 != _i )
+                        for( size_t _l2=_l1+1; _l2<G.nb(j).size(); _l2++) if( _l2 != _i)
                             newM += Gamma(j,_i,_l1,_l2) * tJ[j][_l1] * tJ[j][_l2] * cors[j][_l1][_l2];
                 }
 
                 Real dev = newM - M[i][_j];
 //              dev *= 0.02;
-                if( fabs(dev) >= maxdev )
-                    maxdev = fabs(dev);
+                if( abs(dev) >= maxdev )
+                    maxdev = abs(dev);
 
                 newM = M[i][_j] + dev;
-                if( fabs(newM) > 1.0 )
-                    newM = sign(newM);
+                if( abs(newM) > 1.0 )
+                    newM = (newM > 0.0) ? 1.0 : -1.0;
                 M[i][_j] = newM;
             }
         }
@@ -449,16 +242,16 @@ void MR::solvemcav() {
 
     if(run==maxruns){
         if( props.verbose >= 1 )
-            cerr << "solve_mcav: Convergence not reached (maxdev=" << maxdev << ")..." << endl;
+            cerr << "MR::propagateCavityFields: Convergence not reached (maxdev=" << maxdev << ")..." << endl;
     }
 }
 
 
-void MR::solveM() {
-    for(size_t i=0; i<N; i++) {
+void MR::calcMagnetizations() {
+    for( size_t i = 0; i < G.nrNodes(); i++ ) {
         if( props.updates == Properties::UpdateType::FULL ) {
-            // find indices in nb[i]
-            sub_nb _nbi(con[i]);
+            // find indices in nb(i)
+            sub_nb _nbi( G.nb(i).size() );
             _nbi.set();
 
             // calc numerator1 and denominator1
@@ -468,45 +261,84 @@ void MR::solveM() {
             Mag[i] = (tanh(theta[i]) * sum_even + sum_odd) / (sum_even + tanh(theta[i]) * sum_odd);
 
         } else if( props.updates == Properties::UpdateType::LINEAR ) {
-            sub_nb empty(con[i]);
+            sub_nb empty( G.nb(i).size() );
             Mag[i] = T(i,empty);
 
-            for(size_t _l1=0; _l1<con[i]; _l1++)
-                for( size_t _l2=_l1+1; _l2<con[i]; _l2++)
+            for( size_t _l1 = 0; _l1 < G.nb(i).size(); _l1++ )
+                for( size_t _l2 = _l1 + 1; _l2 < G.nb(i).size(); _l2++ )
                     Mag[i] += Gamma(i,_l1,_l2) * tJ[i][_l1] * tJ[i][_l2] * cors[i][_l1][_l2];
         }
-        if(fabs(Mag[i])>1.)
-            Mag[i] = sign(Mag[i]);
+        if( abs( Mag[i] ) > 1.0 )
+            Mag[i] = (Mag[i] > 0.0) ? 1.0 : -1.0;
     }
 }
 
 
-Real MR::init_cor() {
+Real MR::calcCavityCorrelations() {
     Real md = 0.0;
     for( size_t i = 0; i < nrVars(); i++ ) {
         vector<Factor> pairq;
-        if( props.inits == Properties::InitType::CLAMPING ) {
+        if( props.inits == Properties::InitType::EXACT ) {
+            JTree jtcav(*this, PropertySet()("updates", string("HUGIN"))("verbose", (size_t)0) );
+            jtcav.makeCavity( i );
+            pairq = calcPairBeliefs( jtcav, delta(i), false, true );
+        } else if( props.inits == Properties::InitType::CLAMPING ) {
             BP bpcav(*this, PropertySet()("updates", string("SEQMAX"))("tol", (Real)1.0e-9)("maxiter", (size_t)10000)("verbose", (size_t)0)("logdomain", false));
             bpcav.makeCavity( i );
+
             pairq = calcPairBeliefs( bpcav, delta(i), false, true );
             md = std::max( md, bpcav.maxDiff() );
-        } else if( props.inits == Properties::InitType::EXACT ) {
-            JTree jtcav(*this, PropertySet()("updates", string("HUGIN"))("verbose", (size_t)0) );
-            jtcav.makeCavity( i );
-            pairq = calcPairBeliefs( jtcav, delta(i), false, true );
-        }
-        for( size_t jk = 0; jk < pairq.size(); jk++ ) {
-            VarSet::const_iterator kit = pairq[jk].vars().begin();
-            size_t j = findVar( *(kit) );
-            size_t k = findVar( *(++kit) );
-            pairq[jk].normalize();
-            Real cor = (pairq[jk][3] - pairq[jk][2] - pairq[jk][1] + pairq[jk][0]) - (pairq[jk][3] + pairq[jk][2] - pairq[jk][1] - pairq[jk][0]) * (pairq[jk][3] - pairq[jk][2] + pairq[jk][1] - pairq[jk][0]);
-            for( size_t _j = 0; _j < con[i]; _j++ ) if( nb[i][_j] == j )
-                for( size_t _k = 0; _k < con[i]; _k++ ) if( nb[i][_k] == k ) {
-                    cors[i][_j][_k] = cor;
-                    cors[i][_k][_j] = cor;
+        } else if( props.inits == Properties::InitType::RESPPROP ) {
+            BP bpcav(*this, PropertySet()("updates", string("SEQMAX"))("tol", (Real)1.0e-9)("maxiter", (size_t)10000)("verbose", (size_t)0)("logdomain", false));
+            bpcav.makeCavity( i );
+            bpcav.makeCavity( i );
+            bpcav.init();
+            bpcav.run();
+
+            BBP bbp( &bpcav, PropertySet()("verbose",(size_t)0)("tol",(Real)1.0e-9)("maxiter",(size_t)10000)("damping",(Real)0.0)("updates",string("SEQ_MAX")) );
+            foreach( const Neighbor &j, G.nb(i) ) {
+                // Create weights for magnetization of some spin
+                Prob p( 2, 0.0 );
+                p[0] = -1.0;
+                p[1] = 1.0;
+
+                // BBP cost function would be the magnetization of spin j
+                vector<Prob> b1_adj;
+                b1_adj.reserve( nrVars() );
+                for( size_t l = 0; l < nrVars(); l++ )
+                    if( l == j )
+                        b1_adj.push_back( p );
+                    else
+                        b1_adj.push_back( Prob( 2, 0.0 ) );
+                bbp.init_V( b1_adj );
+
+                // run BBP to estimate adjoints
+                bbp.run();
+
+                foreach( const Neighbor &k, G.nb(i) ) {
+                    if( k != j )
+                        cors[i][j.iter][k.iter] = (bbp.adj_psi_V(k)[1] - bbp.adj_psi_V(k)[0]);
+                    else
+                        cors[i][j.iter][k.iter] = 0.0;
                 }
+            }
         }
+
+        if( props.inits != Properties::InitType::RESPPROP ) {
+            for( size_t jk = 0; jk < pairq.size(); jk++ ) {
+                VarSet::const_iterator kit = pairq[jk].vars().begin();
+                size_t j = findVar( *(kit) );
+                size_t k = findVar( *(++kit) );
+                pairq[jk].normalize();
+                Real cor = (pairq[jk][3] - pairq[jk][2] - pairq[jk][1] + pairq[jk][0]) - (pairq[jk][3] + pairq[jk][2] - pairq[jk][1] - pairq[jk][0]) * (pairq[jk][3] - pairq[jk][2] + pairq[jk][1] - pairq[jk][0]);
+
+                size_t _j = G.findNb(i,j);
+                size_t _k = G.findNb(i,k);
+                cors[i][_j][_k] = cor;
+                cors[i][_k][_j] = cor;
+            }
+        }
+
     }
     return md;
 }
@@ -524,44 +356,16 @@ Real MR::run() {
 
         double tic = toc();
 
-        M.resize(N);
-        for(size_t i=0; i<N; i++)
-          M[i].resize(kmax);
-
-        cors.resize(N);
-        for(size_t i=0; i<N; i++)
-          cors[i].resize(kmax);
-        for(size_t i=0; i<N; i++)
-          for(size_t j=0; j<kmax; j++)
-            cors[i][j].resize(kmax);
-
-        kindex.resize(N);
-        for(size_t i=0; i<N; i++)
-          kindex[i].resize(kmax);
-
-        if( props.inits == Properties::InitType::RESPPROP ) {
-            Real md = init_cor_resp();
-            if( md > _maxdiff )
-                _maxdiff = md;
-        } else if( props.inits == Properties::InitType::RESPPROPOLD ) {
-            Real md = init_cor_resp_old();
-            if( md > _maxdiff )
-                _maxdiff = md;
-        } else if( props.inits == Properties::InitType::EXACT ) {
-            Real md = init_cor();
-            if( md > _maxdiff )
-                _maxdiff = md;
-        }
-        else if( props.inits == Properties::InitType::CLAMPING ) {
-            Real md = init_cor();
-            if( md > _maxdiff )
-                _maxdiff = md;
-        }
+        // approximate correlations of cavity spins
+        Real md = calcCavityCorrelations();
+        if( md > _maxdiff )
+            _maxdiff = md;
 
-        solvemcav();
+        // solve messages
+        propagateCavityFields();
 
-        Mag.resize(N);
-        solveM();
+        // calculate magnetizations
+        calcMagnetizations();
 
         if( props.verbose >= 1 )
             cerr << Name << " needed " << toc() - tic << " seconds." << endl;
@@ -572,18 +376,6 @@ Real MR::run() {
 }
 
 
-void MR::makekindex() {
-    for(size_t i=0; i<N; i++)
-        for(size_t j=0; j<con[i]; j++) {
-            size_t ij = nb[i][j];       // ij is the j'th neighbour of spin i
-            size_t k=0;
-            while( nb[ij][k] != i )
-                k++;
-            kindex[i][j] = k;   // the j'th neighbour of spin i has spin i as its k'th neighbour
-        }
-}
-
-
 Factor MR::beliefV( size_t i ) const {
     if( supported ) {
         Prob x(2);
@@ -599,70 +391,90 @@ Factor MR::beliefV( size_t i ) const {
 vector<Factor> MR::beliefs() const {
     vector<Factor> result;
     for( size_t i = 0; i < nrVars(); i++ )
-        result.push_back( belief( var(i) ) );
+        result.push_back( beliefV( i ) );
     return result;
 }
 
 
-
 MR::MR( const FactorGraph &fg, const PropertySet &opts ) : DAIAlgFG(fg), supported(true), _maxdiff(0.0), _iters(0) {
     setProperties( opts );
 
+    size_t N = fg.nrVars();
+
     // check whether all vars in fg are binary
-    // check whether connectivity is <= kmax
-    for( size_t i = 0; i < fg.nrVars(); i++ )
-        if( (fg.var(i).states() > 2) || (fg.delta(i).size() > kmax) ) {
+    for( size_t i = 0; i < N; i++ )
+        if( (fg.var(i).states() > 2) ) {
             supported = false;
             break;
         }
-
     if( !supported )
-        DAI_THROWE(NOT_IMPLEMENTED,"MR only supports binary variables with low connectivity");
+        DAI_THROWE(NOT_IMPLEMENTED,"MR only supports binary variables");
 
     // check whether all interactions are pairwise or single
-    for( size_t I = 0; I < fg.nrFactors(); I++ )
-        if( fg.factor(I).vars().size() > 2 ) {
+    // and construct Markov graph
+    G = GraphAL(N);
+    for( size_t I = 0; I < fg.nrFactors(); I++ ) {
+        const Factor &psi = fg.factor(I);
+        if( psi.vars().size() > 2 ) {
             supported = false;
             break;
+        } else if( psi.vars().size() == 2 ) {
+            VarSet::const_iterator jit = psi.vars().begin();
+            size_t i = fg.findVar( *(jit) );
+            size_t j = fg.findVar( *(++jit) );
+            G.addEdge( i, j, false );
         }
-
+    }
     if( !supported )
         DAI_THROWE(NOT_IMPLEMENTED,"MR does not support higher order interactions (only single and pairwise are supported)");
 
-    // create w and th
-    size_t Nin = fg.nrVars();
-
-    Real *w = new Real[Nin*Nin];
-    Real *th = new Real[Nin];
+    // construct theta
+    theta.clear();
+    theta.resize( N, 0.0 );
 
-    for( size_t i = 0; i < Nin; i++ ) {
-        th[i] = 0.0;
-        for( size_t j = 0; j < Nin; j++ )
-            w[i*Nin+j] = 0.0;
-    }
+    // construct tJ
+    tJ.resize( N );
+    for( size_t i = 0; i < N; i++ )
+        tJ[i].resize( G.nb(i).size(), 0.0 );
 
+    // initialize theta and tJ
     for( size_t I = 0; I < fg.nrFactors(); I++ ) {
         const Factor &psi = fg.factor(I);
         if( psi.vars().size() == 1 ) {
             size_t i = fg.findVar( *(psi.vars().begin()) );
-            th[i] += 0.5 * log(psi[1] / psi[0]);
+            theta[i] += 0.5 * log(psi[1] / psi[0]);
         } else if( psi.vars().size() == 2 ) {
-            size_t i = fg.findVar( *(psi.vars().begin()) );
             VarSet::const_iterator jit = psi.vars().begin();
+            size_t i = fg.findVar( *(jit) );
             size_t j = fg.findVar( *(++jit) );
 
-            w[i*Nin+j] += 0.25 * log(psi[3] * psi[0] / (psi[2] * psi[1]));
-            w[j*Nin+i] += 0.25 * log(psi[3] * psi[0] / (psi[2] * psi[1]));
+            Real w_ij = 0.25 * log(psi[3] * psi[0] / (psi[2] * psi[1]));
+            tJ[i][G.findNb(i,j)] += w_ij;
+            tJ[j][G.findNb(j,i)] += w_ij;
 
-            th[i] += 0.25 * log(psi[3] / psi[2] * psi[1] / psi[0]);
-            th[j] += 0.25 * log(psi[3] / psi[1] * psi[2] / psi[0]);
+            theta[i] += 0.25 * log(psi[3] / psi[2] * psi[1] / psi[0]);
+            theta[j] += 0.25 * log(psi[3] / psi[1] * psi[2] / psi[0]);
         }
     }
-
-    init(Nin, w, th);
-
-    delete th;
-    delete w;
+    for( size_t i = 0; i < N; i++ )
+        foreach( const Neighbor &j, G.nb(i) )
+            tJ[i][j.iter] = tanh( tJ[i][j.iter] );
+
+    // construct M
+    M.resize( N );
+    for( size_t i = 0; i < N; i++ )
+        M[i].resize( G.nb(i).size() );
+
+    // construct cors
+    cors.resize( N );
+    for( size_t i = 0; i < N; i++ )
+        cors[i].resize( G.nb(i).size() );
+    for( size_t i = 0; i < N; i++ )
+        for( size_t _j = 0; _j < cors[i].size(); _j++ )
+            cors[i][_j].resize( G.nb(i).size() );
+
+    // construct Mag
+    Mag.resize( N );
 }
 
 
index c3f5bbf..4707b5f 100644 (file)
@@ -969,7 +969,7 @@ BBP                                         7.049e-04       2.319e-04       +1.209e-03      1.000e-09
 # ({x13}, (5.351e-01, 4.649e-01))
 # ({x14}, (6.284e-01, 3.716e-01))
 # ({x15}, (1.355e-01, 8.645e-01))
-CBP                                            3.294e-06       8.863e-07       +3.531e-06      1.000e-09       
+CBP                                            5.336e-06       2.335e-06       -6.764e-06      1.000e-09       
 # ({x0}, (3.888e-01, 6.112e-01))
 # ({x1}, (5.556e-01, 4.444e-01))
 # ({x2}, (4.587e-01, 5.413e-01))