Multiple changes: changes in build system, one workaround and one bug fix
[libdai.git] / src / fbp.cpp
index 0e67cc5..c1adbe2 100644 (file)
@@ -1,14 +1,15 @@
 /*  This file is part of libDAI - http://www.libdai.org/
  *
- *  libDAI is licensed under the terms of the GNU General Public License version
- *  2, or (at your option) any later version. libDAI is distributed without any
- *  warranty. See the file COPYING for more details.
+ *  Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
  *
- *  Copyright (C) 2009  Frederik Eaton
- *  Copyright (C) 2009  Joris Mooij  [joris dot mooij at libdai dot org]
+ *  Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
  */
 
 
+#include <dai/dai_config.h>
+#ifdef DAI_WITH_FBP
+
+
 #include <dai/fbp.h>
 
 
@@ -21,86 +22,91 @@ namespace dai {
 using namespace std;
 
 
-const char *FBP::Name = "FBP";
-
-
-string FBP::identify() const {
-    return string(Name) + printProperties();
-}
-
-
-/* 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++ ) {
         sum += (beliefF(I) * factor(I).log(true)).sum();  // FBP
-        sum += scaleF(I) * beliefF(I).entropy();  // FBP
+        sum += Weight(I) * beliefF(I).entropy();  // FBP
     }
     for( size_t i = 0; i < nrVars(); ++i ) {
         Real c_i = 0.0;
-        foreach( const Neighbor &I, nbV(i) )
-            c_i += scaleF(I);
-        sum += (1.0 - c_i) * beliefV(i).entropy();  // FBP
+        bforeach( const Neighbor &I, nbV(i) )
+            c_i += Weight(I);
+        if( c_i != 1.0 )
+            sum += (1.0 - c_i) * beliefV(i).entropy();  // FBP
     }
     return sum;
 }
 
 
-/* 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 scale = scaleF(I); // FBP: c_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 *= (1/scale); // FBP
+        prod /= c_I; // FBP
     } else
-        prod ^= (1/scale); // FBP
-    
+        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
+    bforeach( const Neighbor &j, nbF(I) )
+        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) )
+            bforeach( 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 );
+                } else if( c_I != 1.0 ) {
+                    // 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);
                 }
 
-
-            size_t _I = j.dual;
-            // FBP: now multiply by m_Ij^(1-1/c_I)
-            if(props.logdomain)
-                prod_j += message( j, _I)*(1-1/scale);
-            else
-                prod_j *= message( j, _I)^(1-1/scale);
-
             // 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]];
+                        prod.set( r, prod[r] + prod_j[ind[r]] );
                     else
-                        prod[r] *= prod_j[ind[r]];
+                        prod.set( r, 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();
@@ -110,28 +116,28 @@ 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);
         if( props.inference == Properties::InfType::SUMPROD )
             for( size_t r = 0; r < prod.size(); ++r )
-                marg[ind[r]] += prod[r];
+                marg.set( ind[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.set( ind[r], prod[r] );
         marg.normalize();
     }
 
     // FBP
-    marg ^= scale;
+    marg ^= c_I;
 
     // Store result
     if( props.logdomain )
@@ -141,73 +147,17 @@ void FBP::calcNewMessage( size_t i, size_t _I ) {
 
     // Update the residual if necessary
     if( props.updates == Properties::UpdateType::SEQMAX )
-        updateResidual( i, _I , dist( newMessage( i, _I ), message( i, _I ), Prob::DISTLINF ) );
+        updateResidual( i, _I , dist( newMessage( i, _I ), message( i, _I ), DISTLINF ) );
 }
 
 
-/* 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 scale = scaleF(I); // FBP: c_I
-
-    Factor Fprod( factor(I) );
-    Prob &prod = Fprod.p();
-    
-    if( props.logdomain ) {
-        prod.takeLog();
-        prod /= scale; // FBP
-    } else
-        prod ^= (1/scale); // 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 );
-            }
-
-        size_t _I = j.dual;
-
-        // FBP: now multiply by m_Ij^(1-1/c_I)
-        if( props.logdomain )
-            prod_j += newMessage( j, _I)*(1-1/scale);
-        else
-            prod_j *= newMessage( j, _I)^(1-1/scale);
-
-        // 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();
-    _scale_factor.resize( nrFactors(), 1.0 );
+    _weight.resize( nrFactors(), 1.0 );
 }
 
 
 } // end of namespace dai
+
+
+#endif