Merged duplicate code (in calcBeliefF() and calcNewMessage()) in BP,FBP,TRWBP
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Wed, 3 Feb 2010 13:22:58 +0000 (14:22 +0100)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Wed, 3 Feb 2010 13:22:58 +0000 (14:22 +0100)
include/dai/bp.h
include/dai/fbp.h
include/dai/trwbp.h
src/bp.cpp
src/fbp.cpp
src/trwbp.cpp

index aa0441d..aa3e329 100644 (file)
@@ -57,8 +57,6 @@ namespace dai {
  *  \note There are two implementations, an optimized one (the default) which caches IndexFor objects,
  *  and a slower, less complicated one which is easier to maintain/understand. The slower one can be 
  *  enabled by defining DAI_BP_FAST as false in the source file.
- *
- *  \todo Merge duplicate code in calcNewMessage() and calcBeliefF()
  */
 class BP : public DAIAlgFG {
     protected:
@@ -232,6 +230,11 @@ class BP : public DAIAlgFG {
         /// Returns reference to residual for the edge between variable \a i and its \a _I 'th neighbor
         Real & residual(size_t i, size_t _I) { return _edges[i][_I].residual; }
 
+        /// Calculate the product of factor \a I and the incoming messages
+        /** If \a without_i == \c true, the message coming from variable \a i is omitted from the product
+         *  \note This function is used by calcNewMessage() and calcBeliefF()
+         */
+        virtual Prob calcIncomingMessageProduct( size_t I, bool without_i, size_t i ) const;
         /// Calculate the updated message from the \a _I 'th neighbor of variable \a i to variable \a i
         virtual void calcNewMessage( size_t i, size_t _I );
         /// Replace the "old" message from the \a _I 'th neighbor of variable \a i to variable \a i by the "new" (updated) message
@@ -243,7 +246,9 @@ class BP : public DAIAlgFG {
         /// Calculates unnormalized belief of variable \a i
         virtual void calcBeliefV( size_t i, Prob &p ) const;
         /// Calculates unnormalized belief of factor \a I
-        virtual void calcBeliefF( size_t I, Prob &p ) const;
+        virtual void calcBeliefF( size_t I, Prob &p ) const {
+            p = calcIncomingMessageProduct( I, false, 0 );
+        }
 
         /// Helper function for constructors
         virtual void construct();
index bcfa83a..85382d3 100644 (file)
@@ -52,7 +52,6 @@ namespace dai {
  *    \f[ c_i := \sum_{I \in N_i} c_I \f]
  *
  *  \todo Add nice way to set weights
- *  \todo Merge duplicate code in calcNewMessage() and calcBeliefF()
  *  \author Frederik Eaton
  */
 class FBP : public BP {
@@ -103,11 +102,19 @@ class FBP : public BP {
         void setWeights( const std::vector<Real> &c ) { _weight = c; }
 
     protected:
+        /// Calculate the product of factor \a I and the incoming messages
+        /** If \a without_i == \c true, the message coming from variable \a i is omitted from the product
+         *  \note This function is used by calcNewMessage() and calcBeliefF()
+         */
+        virtual Prob calcIncomingMessageProduct( size_t I, bool without_i, size_t i ) const;
+
         // Calculate the updated message from the \a _I 'th neighbor of variable \a i to variable \a i
         virtual void calcNewMessage( size_t i, size_t _I );
 
         // Calculates unnormalized belief of factor \a I
-        virtual void calcBeliefF( size_t I, Prob &p ) const;
+        virtual void calcBeliefF( size_t I, Prob &p ) const {
+            p = calcIncomingMessageProduct( I, false, 0 );
+        }
 
         // Helper function for constructors
         virtual void construct();
index 83f39d5..57625de 100644 (file)
@@ -45,7 +45,6 @@ namespace dai {
  *    \f[ c_i := \sum_{I \in N_i} c_I \f]
  *
  *  \note TRWBP is actually equivalent to FBP
- *  \todo Merge duplicate code in calcNewMessage() and calcBeliefF()
  *  \todo Add nice way to set weights
  *  \todo Merge code of FBP and TRWBP
  */
@@ -102,14 +101,19 @@ class TRWBP : public BP {
         void setWeights( const std::vector<Real> &c ) { _weight = c; }
 
     protected:
-        // Calculate the updated message from the \a _I 'th neighbor of variable \a i to variable \a i
-        virtual void calcNewMessage( size_t i, size_t _I );
+        /// Calculate the product of factor \a I and the incoming messages
+        /** If \a without_i == \c true, the message coming from variable \a i is omitted from the product
+         *  \note This function is used by calcNewMessage() and calcBeliefF()
+         */
+        virtual Prob calcIncomingMessageProduct( size_t I, bool without_i, size_t i ) const;
 
         /// Calculates unnormalized belief of variable \a i
         virtual void calcBeliefV( size_t i, Prob &p ) const;
 
         // Calculates unnormalized belief of factor \a I
-        virtual void calcBeliefF( size_t I, Prob &p ) const;
+        virtual void calcBeliefF( size_t I, Prob &p ) const {
+            p = calcIncomingMessageProduct( I, false, 0 );
+        }
 
         // Helper function for constructors
         virtual void construct();
index dc7dd91..9c2ea90 100644 (file)
@@ -141,6 +141,49 @@ void BP::findMaxResidual( size_t &i, size_t &_I ) {
 }
 
 
+Prob BP::calcIncomingMessageProduct( size_t I, bool without_i, size_t i ) const {
+    Factor Fprod( factor(I) );
+    Prob &prod = Fprod.p();
+    if( props.logdomain )
+        prod.takeLog();
+
+    // Calculate product of incoming messages and factor I
+    foreach( const Neighbor &j, nbF(I) )
+        if( !(without_i && (j == i)) ) {
+            // prod_j will be the product of messages coming into j
+            Prob prod_j( var(j).states(), props.logdomain ? 0.0 : 1.0 );
+            foreach( const Neighbor &J, nbV(j) )
+                if( J != I ) { // for all J in nb(j) \ I
+                    if( props.logdomain )
+                        prod_j += message( j, J.iter );
+                    else
+                        prod_j *= message( j, J.iter );
+                }
+
+            // multiply prod with prod_j
+            if( !DAI_BP_FAST ) {
+                // UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION
+                if( props.logdomain )
+                    Fprod += Factor( var(j), prod_j );
+                else
+                    Fprod *= Factor( var(j), prod_j );
+            } else {
+                // OPTIMIZED VERSION
+                size_t _I = j.dual;
+                // ind is the precalculated IndexFor(j,I) i.e. to x_I == k corresponds x_j == ind[k]
+                const ind_t &ind = index(j, _I);
+
+                for( size_t r = 0; r < prod.size(); ++r )
+                    if( props.logdomain )
+                        prod[r] += prod_j[ind[r]];
+                    else
+                        prod[r] *= prod_j[ind[r]];
+            }
+    }
+    return prod;
+}
+
+
 void BP::calcNewMessage( size_t i, size_t _I ) {
     // calculate updated message I->i
     size_t I = nbV(i,_I);
@@ -151,41 +194,7 @@ void BP::calcNewMessage( size_t i, size_t _I ) {
     else {
         Factor Fprod( factor(I) );
         Prob &prod = Fprod.p();
-        if( props.logdomain )
-            prod.takeLog();
-
-        // Calculate product of incoming messages and factor I
-        foreach( const Neighbor &j, nbF(I) )
-            if( j != i ) { // for all j in I \ i
-                // prod_j will be the product of messages coming into j
-                Prob prod_j( var(j).states(), props.logdomain ? 0.0 : 1.0 );
-                foreach( const Neighbor &J, nbV(j) )
-                    if( J != I ) { // for all J in nb(j) \ I
-                        if( props.logdomain )
-                            prod_j += message( j, J.iter );
-                        else
-                            prod_j *= message( j, J.iter );
-                    }
-
-                // multiply prod with prod_j
-                if( !DAI_BP_FAST ) {
-                    /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
-                    if( props.logdomain )
-                        Fprod += Factor( var(j), prod_j );
-                    else
-                        Fprod *= Factor( var(j), prod_j );
-                } else {
-                    /* OPTIMIZED VERSION */
-                    size_t _I = j.dual;
-                    // ind is the precalculated IndexFor(j,I) i.e. to x_I == k corresponds x_j == ind[k]
-                    const ind_t &ind = index(j, _I);
-                    for( size_t r = 0; r < prod.size(); ++r )
-                        if( props.logdomain )
-                            prod[r] += prod_j[ind[r]];
-                        else
-                            prod[r] *= prod_j[ind[r]];
-                }
-            }
+        prod = calcIncomingMessageProduct( I, true, i );
 
         if( props.logdomain ) {
             prod -= prod.max();
@@ -349,50 +358,6 @@ void BP::calcBeliefV( size_t i, Prob &p ) const {
 }
 
 
-void BP::calcBeliefF( size_t I, Prob &p ) const {
-    Factor Fprod( factor( I ) );
-    Prob &prod = Fprod.p();
-
-    if( props.logdomain )
-        prod.takeLog();
-
-    foreach( const Neighbor &j, nbF(I) ) {
-        // prod_j will be the product of messages coming into j
-        Prob prod_j( var(j).states(), props.logdomain ? 0.0 : 1.0 );
-        foreach( const Neighbor &J, nbV(j) )
-            if( J != I ) { // for all J in nb(j) \ I
-                if( props.logdomain )
-                    prod_j += newMessage( j, J.iter );
-                else
-                    prod_j *= newMessage( j, J.iter );
-            }
-
-        // multiply prod with prod_j
-        if( !DAI_BP_FAST ) {
-            /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
-            if( props.logdomain )
-                Fprod += Factor( var(j), prod_j );
-            else
-                Fprod *= Factor( var(j), prod_j );
-        } else {
-            /* OPTIMIZED VERSION */
-            size_t _I = j.dual;
-            // ind is the precalculated IndexFor(j,I) i.e. to x_I == k corresponds x_j == ind[k]
-            const ind_t & ind = index(j, _I);
-
-            for( size_t r = 0; r < prod.size(); ++r ) {
-                if( props.logdomain )
-                    prod[r] += prod_j[ind[r]];
-                else
-                    prod[r] *= prod_j[ind[r]];
-            }
-        }
-    }
-
-    p = prod;
-}
-
-
 Factor BP::beliefV( size_t i ) const {
     Prob p;
     calcBeliefV( i, p );
index 8cd51b0..b5e1a2b 100644 (file)
@@ -29,7 +29,7 @@ string FBP::identify() const {
 }
 
 
-/* This code has been copied from bp.cpp, except where comments indicate FBP-specific behaviour */
+// This code has been copied from bp.cpp, except where comments indicate FBP-specific behaviour
 Real FBP::logZ() const {
     Real sum = 0.0;
     for( size_t I = 0; I < nrFactors(); I++ ) {
@@ -46,26 +46,24 @@ Real FBP::logZ() const {
 }
 
 
-/* This code has been copied from bp.cpp, except where comments indicate FBP-specific behaviour */
-void FBP::calcNewMessage( size_t i, size_t _I ) {
-    // calculate updated message I->i
-    size_t I = nbV(i,_I);
-
+// This code has been copied from bp.cpp, except where comments indicate FBP-specific behaviour
+Prob FBP::calcIncomingMessageProduct( size_t I, bool without_i, size_t i ) const {
     Real c_I = Weight(I); // FBP: c_I
 
     Factor Fprod( factor(I) );
     Prob &prod = Fprod.p();
+
     if( props.logdomain ) {
         prod.takeLog();
         prod /= c_I; // FBP
     } else
         prod ^= (1.0 / c_I); // FBP
-    
+
     // Calculate product of incoming messages and factor I
     foreach( const Neighbor &j, nbF(I) )
-        if( j != i ) { // for all j in I \ i
-            // FBP: same as n_jI
+        if( !(without_i && (j == i)) ) {
             // prod_j will be the product of messages coming into j
+            // FBP: corresponds to messages n_jI
             Prob prod_j( var(j).states(), props.logdomain ? 0.0 : 1.0 );
             foreach( const Neighbor &J, nbV(j) )
                 if( J != I ) { // for all J in nb(j) \ I
@@ -76,30 +74,45 @@ void FBP::calcNewMessage( size_t i, size_t _I ) {
                 } else {
                     // FBP: multiply by m_Ij^(1-1/c_I)
                     if( props.logdomain )
-                        prod_j += message( j, J.iter ) * (1.0 - 1.0 / c_I);
+                        prod_j += newMessage( j, J.iter) * (1.0 - 1.0 / c_I);
                     else
-                        prod_j *= message( j, J.iter ) ^ (1.0 - 1.0 / c_I);
+                        prod_j *= newMessage( j, J.iter) ^ (1.0 - 1.0 / c_I);
                 }
 
             // multiply prod with prod_j
             if( !DAI_FBP_FAST ) {
-                /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
+                // UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION
                 if( props.logdomain )
                     Fprod += Factor( var(j), prod_j );
                 else
                     Fprod *= Factor( var(j), prod_j );
             } else {
-                /* OPTIMIZED VERSION */
+                // OPTIMIZED VERSION
                 size_t _I = j.dual;
                 // ind is the precalculated IndexFor(j,I) i.e. to x_I == k corresponds x_j == ind[k]
                 const ind_t &ind = index(j, _I);
+
                 for( size_t r = 0; r < prod.size(); ++r )
                     if( props.logdomain )
                         prod[r] += prod_j[ind[r]];
                     else
                         prod[r] *= prod_j[ind[r]];
             }
-        }
+    }
+    return prod;
+}
+
+
+// This code has been copied from bp.cpp, except where comments indicate FBP-specific behaviour
+void FBP::calcNewMessage( size_t i, size_t _I ) {
+    // calculate updated message I->i
+    size_t I = nbV(i,_I);
+
+    Real c_I = Weight(I); // FBP: c_I
+
+    Factor Fprod( factor(I) );
+    Prob &prod = Fprod.p();
+    prod = calcIncomingMessageProduct( I, true, i );
 
     if( props.logdomain ) {
         prod -= prod.max();
@@ -109,13 +122,13 @@ void FBP::calcNewMessage( size_t i, size_t _I ) {
     // Marginalize onto i
     Prob marg;
     if( !DAI_FBP_FAST ) {
-        /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
+        // UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION
         if( props.inference == Properties::InfType::SUMPROD )
             marg = Fprod.marginal( var(i) ).p();
         else
             marg = Fprod.maxMarginal( var(i) ).p();
     } else {
-        /* OPTIMIZED VERSION */
+        // OPTIMIZED VERSION
         marg = Prob( var(i).states(), 0.0 );
         // ind is the precalculated IndexFor(i,I) i.e. to x_I == k corresponds x_i == ind[k]
         const ind_t ind = index(i,_I);
@@ -144,63 +157,6 @@ void FBP::calcNewMessage( size_t i, size_t _I ) {
 }
 
 
-/* This code has been copied from bp.cpp, except where comments indicate FBP-specific behaviour */
-void FBP::calcBeliefF( size_t I, Prob &p ) const {
-    Real c_I = Weight(I); // FBP: c_I
-
-    Factor Fprod( factor(I) );
-    Prob &prod = Fprod.p();
-    
-    if( props.logdomain ) {
-        prod.takeLog();
-        prod /= c_I; // FBP
-    } else
-        prod ^= (1.0 / c_I); // FBP
-
-    foreach( const Neighbor &j, nbF(I) ) {
-        // prod_j will be the product of messages coming into j
-        // FBP: corresponds to messages n_jI
-        Prob prod_j( var(j).states(), props.logdomain ? 0.0 : 1.0 );
-        foreach( const Neighbor &J, nbV(j) )
-            if( J != I ) { // for all J in nb(j) \ I
-                if( props.logdomain )
-                    prod_j += newMessage( j, J.iter );
-                else
-                    prod_j *= newMessage( j, J.iter );
-            } else {
-                // FBP: multiply by m_Ij^(1-1/c_I)
-                if( props.logdomain )
-                    prod_j += newMessage( j, J.iter) * (1.0 - 1.0 / c_I);
-                else
-                    prod_j *= newMessage( j, J.iter) ^ (1.0 - 1.0 / c_I);
-            }
-
-        // multiply prod with prod_j
-        if( !DAI_FBP_FAST ) {
-            /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
-            if( props.logdomain )
-                Fprod += Factor( var(j), prod_j );
-            else
-                Fprod *= Factor( var(j), prod_j );
-        } else {
-            /* OPTIMIZED VERSION */
-            size_t _I = j.dual;
-            // ind is the precalculated IndexFor(j,I) i.e. to x_I == k corresponds x_j == ind[k]
-            const ind_t & ind = index(j, _I);
-
-            for( size_t r = 0; r < prod.size(); ++r ) {
-                if( props.logdomain )
-                    prod[r] += prod_j[ind[r]];
-                else
-                    prod[r] *= prod_j[ind[r]];
-            }
-        }
-    }
-
-    p = prod;
-}
-    
-
 void FBP::construct() {
     BP::construct();
     _weight.resize( nrFactors(), 1.0 );
index bda5bc7..0674e8f 100644 (file)
@@ -28,7 +28,7 @@ string TRWBP::identify() const {
 }
 
 
-/* This code has been copied from bp.cpp, except where comments indicate TRWBP-specific behaviour */
+// This code has been copied from bp.cpp, except where comments indicate TRWBP-specific behaviour
 Real TRWBP::logZ() const {
     Real sum = 0.0;
     for( size_t I = 0; I < nrFactors(); I++ ) {
@@ -45,109 +45,67 @@ Real TRWBP::logZ() const {
 }
 
 
-/* This code has been copied from bp.cpp, except where comments indicate TRWBP-specific behaviour */
-void TRWBP::calcNewMessage( size_t i, size_t _I ) {
-    // calculate updated message I->i
-    size_t I = nbV(i,_I);
-    const Var &v_i = var(i);
-    Real c_I = Weight(I); // TRWBP: c_I (\mu_I in the paper)
-
-    Prob marg;
-    if( factor(I).vars().size() == 1 ) { // optimization
-        marg = factor(I).p();
-    } else {
-        Factor Fprod( factor(I) );
-        Prob &prod = Fprod.p();
-        if( props.logdomain ) {
-            prod.takeLog();
-            prod /= c_I;         // TRWBP
-        } else
-            prod ^= (1.0 / c_I); // TRWBP
-    
-        // Calculate product of incoming messages and factor I
-        foreach( const Neighbor &j, nbF(I) )
-            if( j != i ) { // for all j in I \ i
-                const Var &v_j = var(j);
-
-                // TRWBP: corresponds to messages n_jI
-                // prod_j will be the product of messages coming into j
-                Prob prod_j( v_j.states(), props.logdomain ? 0.0 : 1.0 );
-                foreach( const Neighbor &J, nbV(j) ) {
-                    Real c_J = Weight(J);
-                    if( J != I ) { // for all J in nb(j) \ I
-                        if( props.logdomain )
-                            prod_j += message( j, J.iter ) * c_J;
-                        else
-                            prod_j *= message( j, J.iter ) ^ c_J;
-                    } else { // TRWBP: multiply by m_Ij^(c_I-1)
-                        if( props.logdomain )
-                            prod_j += message( j, J.iter ) * (c_J - 1.0);
-                        else
-                            prod_j *= message( j, J.iter ) ^ (c_J - 1.0);
-                    }
-                }
+// This code has been copied from bp.cpp, except where comments indicate TRWBP-specific behaviour
+Prob TRWBP::calcIncomingMessageProduct( size_t I, bool without_i, size_t i ) const {
+    Real c_I = Weight(I); // TRWBP: c_I
+
+    Factor Fprod( factor(I) );
+    Prob &prod = Fprod.p();
+    if( props.logdomain ) {
+        prod.takeLog();
+        prod /= c_I; // TRWBP
+    } else
+        prod ^= (1.0 / c_I); // TRWBP
 
-                // multiply prod with prod_j
-                if( !DAI_TRWBP_FAST ) {
-                    /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
+    // Calculate product of incoming messages and factor I
+    foreach( const Neighbor &j, nbF(I) )
+        if( !(without_i && (j == i)) ) {
+            const Var &v_j = var(j);
+            // prod_j will be the product of messages coming into j
+            // TRWBP: corresponds to messages n_jI
+            Prob prod_j( v_j.states(), props.logdomain ? 0.0 : 1.0 );
+            foreach( const Neighbor &J, nbV(j) ) {
+                Real c_J = Weight(J);  // TRWBP
+                if( J != I ) { // for all J in nb(j) \ I
+                    if( props.logdomain )
+                        prod_j += message( j, J.iter ) * c_J;
+                    else
+                        prod_j *= message( j, J.iter ) ^ c_J;
+                } else { // TRWBP: multiply by m_Ij^(c_I-1)
                     if( props.logdomain )
-                        Fprod += Factor( v_j, prod_j );
+                        prod_j += message( j, J.iter ) * (c_J - 1.0);
                     else
-                        Fprod *= Factor( v_j, prod_j );
-                } else {
-                    /* OPTIMIZED VERSION */
-                    size_t _I = j.dual;
-                    // ind is the precalculated IndexFor(j,I) i.e. to x_I == k corresponds x_j == ind[k]
-                    const ind_t &ind = index(j, _I);
-                    for( size_t r = 0; r < prod.size(); ++r )
-                        if( props.logdomain )
-                            prod[r] += prod_j[ind[r]];
-                        else
-                            prod[r] *= prod_j[ind[r]];
+                        prod_j *= message( j, J.iter ) ^ (c_J - 1.0);
                 }
             }
 
-        if( props.logdomain ) {
-            prod -= prod.max();
-            prod.takeExp();
-        }
-
-        // Marginalize onto i
-        if( !DAI_TRWBP_FAST ) {
-            /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
-            if( props.inference == Properties::InfType::SUMPROD )
-                marg = Fprod.marginal( v_i ).p();
-            else
-                marg = Fprod.maxMarginal( v_i ).p();
-        } else {
-            /* OPTIMIZED VERSION */
-            marg = Prob( v_i.states(), 0.0 );
-            // ind is the precalculated IndexFor(i,I) i.e. to x_I == k corresponds x_i == ind[k]
-            const ind_t ind = index(i,_I);
-            if( props.inference == Properties::InfType::SUMPROD )
-                for( size_t r = 0; r < prod.size(); ++r )
-                    marg[ind[r]] += prod[r];
-            else
-                for( size_t r = 0; r < prod.size(); ++r )
-                    if( prod[r] > marg[ind[r]] )
-                        marg[ind[r]] = prod[r];
-            marg.normalize();
+            // multiply prod with prod_j
+            if( !DAI_TRWBP_FAST ) {
+                // UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION
+                if( props.logdomain )
+                    Fprod += Factor( v_j, prod_j );
+                else
+                    Fprod *= Factor( v_j, prod_j );
+            } else {
+                // OPTIMIZED VERSION
+                size_t _I = j.dual;
+                // ind is the precalculated IndexFor(j,I) i.e. to x_I == k corresponds x_j == ind[k]
+                const ind_t &ind = index(j, _I);
+
+                for( size_t r = 0; r < prod.size(); ++r ) {
+                    if( props.logdomain )
+                        prod[r] += prod_j[ind[r]];
+                    else
+                        prod[r] *= prod_j[ind[r]];
+                }
+            }
         }
-    }
-
-    // Store result
-    if( props.logdomain )
-        newMessage(i,_I) = marg.log();
-    else
-        newMessage(i,_I) = marg;
-
-    // Update the residual if necessary
-    if( props.updates == Properties::UpdateType::SEQMAX )
-        updateResidual( i, _I , dist( newMessage( i, _I ), message( i, _I ), Prob::DISTLINF ) );
+    
+    return prod;
 }
 
 
-/* This code has been copied from bp.cpp, except where comments indicate TRWBP-specific behaviour */
+// This code has been copied from bp.cpp, except where comments indicate TRWBP-specific behaviour
 void TRWBP::calcBeliefV( size_t i, Prob &p ) const {
     p = Prob( var(i).states(), props.logdomain ? 0.0 : 1.0 );
     foreach( const Neighbor &I, nbV(i) ) {
@@ -160,67 +118,6 @@ void TRWBP::calcBeliefV( size_t i, Prob &p ) const {
 }
 
 
-/* This code has been copied from bp.cpp, except where comments indicate TRWBP-specific behaviour */
-void TRWBP::calcBeliefF( size_t I, Prob &p ) const {
-    Real c_I = Weight(I); // TRWBP: c_I
-
-    Factor Fprod( factor(I) );
-    Prob &prod = Fprod.p();
-
-    if( props.logdomain ) {
-        prod.takeLog();
-        prod /= c_I; // TRWBP
-    } else
-        prod ^= (1.0 / c_I); // TRWBP
-
-    // Calculate product of incoming messages and factor I
-    foreach( const Neighbor &j, nbF(I) ) {
-        const Var &v_j = var(j);
-
-        // TRWBP: corresponds to messages n_jI
-        // prod_j will be the product of messages coming into j
-        Prob prod_j( v_j.states(), props.logdomain ? 0.0 : 1.0 );
-        foreach( const Neighbor &J, nbV(j) ) {
-            Real c_J = Weight(J);
-            if( J != I ) { // for all J in nb(j) \ I
-                if( props.logdomain )
-                    prod_j += newMessage( j, J.iter ) * c_J;
-                else
-                    prod_j *= newMessage( j, J.iter ) ^ c_J;
-            } else { // TRWBP: multiply by m_Ij^(c_I-1)
-                if( props.logdomain )
-                    prod_j += newMessage( j, J.iter ) * (c_J - 1.0);
-                else
-                    prod_j *= newMessage( j, J.iter ) ^ (c_J - 1.0);
-            }
-        }
-
-        // multiply prod with prod_j
-        if( !DAI_TRWBP_FAST ) {
-            /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
-            if( props.logdomain )
-                Fprod += Factor( v_j, prod_j );
-            else
-                Fprod *= Factor( v_j, prod_j );
-        } else {
-            /* OPTIMIZED VERSION */
-            size_t _I = j.dual;
-            // ind is the precalculated IndexFor(j,I) i.e. to x_I == k corresponds x_j == ind[k]
-            const ind_t &ind = index(j, _I);
-
-            for( size_t r = 0; r < prod.size(); ++r ) {
-                if( props.logdomain )
-                    prod[r] += prod_j[ind[r]];
-                else
-                    prod[r] *= prod_j[ind[r]];
-            }
-        }
-    }
-    
-    p = prod;
-}
-
-
 void TRWBP::construct() {
     BP::construct();
     _weight.resize( nrFactors(), 1.0 );