Fixed a bug in JTree::findMaximum() (reported by zhengyun84 and Dhruv Batra):
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Wed, 11 Aug 2010 07:16:10 +0000 (09:16 +0200)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Wed, 11 Aug 2010 07:16:10 +0000 (09:16 +0200)
if one or more variables had a MAP belief with more than one maximum, an
incorrect MAP state could result.

12 files changed:
ChangeLog
Makefile
README
examples/example.cpp
include/dai/bp.h
include/dai/daialg.h
include/dai/doc.h
include/dai/jtree.h
src/bp.cpp
src/daialg.cpp
src/jtree.cpp
tests/jtreemapbug.fg [new file with mode: 0644]

index a9aad67..1e24851 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,11 @@
+git master
+----------
+
+* Fixed a bug in JTree::findMaximum() (reported by zhengyun84 and Dhruv Batra):
+  if one or more variables had a MAP belief with more than one maximum, an 
+  incorrect MAP state could result.
+
+
 libDAI-0.2.6 (2010-08-05)
 -------------------------
 
index df6d0fb..50fbf0f 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -15,8 +15,8 @@ include Makefile.ALL
 include Makefile.conf
 
 # Set version and date
-DAI_VERSION="0.2.6"
-DAI_DATE="August 5, 2010"
+DAI_VERSION="git master"
+DAI_DATE="August 11, 2010 (or later)"
 
 # Directories of libDAI sources
 # Location of libDAI headers
diff --git a/README b/README
index e02ea38..c3c674c 100644 (file)
--- a/README
+++ b/README
@@ -2,8 +2,8 @@ libDAI  -  A free/open source C++ library for Discrete Approximate Inference
 
 -------------------------------------------------------------------------------
 
-Version:  0.2.6
-Date:     August 5, 2010
+Version:  git master
+Date:     August 11, 2010 (or later)
 See also: http://www.libdai.org
 
 -------------------------------------------------------------------------------
index a7dae54..86cf473 100644 (file)
@@ -133,12 +133,12 @@ int main( int argc, char *argv[] ) {
         // Report exact MAP factor marginals
         cout << "Exact MAP factor marginals:" << endl;
         for( size_t I = 0; I < fg.nrFactors(); I++ )
-            cout << jtmap.belief(fg.factor(I).vars()) << "=" << jtmap.beliefF(I) << endl;
+            cout << jtmap.belief(fg.factor(I).vars()) << " == " << jtmap.beliefF(I) << endl;
 
         // Report max-product factor marginals
         cout << "Approximate (max-product) MAP factor marginals:" << endl;
         for( size_t I = 0; I < fg.nrFactors(); I++ )
-            cout << mp.belief(fg.factor(I).vars()) << "=" << mp.beliefF(I) << endl;
+            cout << mp.belief(fg.factor(I).vars()) << " == " << mp.beliefF(I) << endl;
 
         // Report exact MAP joint state
         cout << "Exact MAP state (log score = " << fg.logScore( jtmapstate ) << "):" << endl;
index deb4fee..1400ca3 100644 (file)
@@ -199,7 +199,7 @@ class BP : public DAIAlgFG {
         virtual Real logZ() const;
         /** \pre Assumes that run() has been called and that \a props.inference == \c MAXPROD
          */
-        std::vector<std::size_t> findMaximum() const;
+        std::vector<std::size_t> findMaximum() const { return dai::findMaximum( *this ); }
         virtual void init();
         virtual void init( const VarSet &ns );
         virtual Real run();
index df88c1f..3a28712 100644 (file)
@@ -254,6 +254,7 @@ typedef DAIAlg<RegionGraph> DAIAlgRG;
  */
 Factor calcMarginal( const InfAlg& obj, const VarSet& vs, bool reInit );
 
+
 /// Calculates beliefs for all pairs of variables in \a vs using inference algorithm \a obj.
 /** calcPairBeliefs() works by 
  *  - clamping single variables in \a vs and calculating the partition sum and the single variable beliefs for each clamped state, if \a accurate == \c false;
@@ -269,6 +270,12 @@ Factor calcMarginal( const InfAlg& obj, const VarSet& vs, bool reInit );
 std::vector<Factor> calcPairBeliefs( const InfAlg& obj, const VarSet& vs, bool reInit, bool accurate=false );
 
 
+/// Calculates the joint state of all variables that has maximum probability, according to the inference algorithm \a obj
+/** \note Before this method is called, obj.run() should have been called.
+ */
+std::vector<size_t> findMaximum( const InfAlg& obj );
+
+
 } // end of namespace dai
 
 
index 064163b..7ce2486 100644 (file)
@@ -35,8 +35,8 @@
 
 /** \mainpage Reference manual for libDAI - A free/open source C++ library for Discrete Approximate Inference methods
  *  \author Joris Mooij (with contributions of Frederik Eaton)
- *  \version 0.2.6
- *  \date August 5, 2010
+ *  \version git master
+ *  \date August 11, 2010 (or later)
  *
  *  <hr size="1">
  *  \section about About libDAI
index 5ec7f9c..00364b7 100644 (file)
@@ -130,7 +130,7 @@ class JTree : public DAIAlgRG {
         virtual Real logZ() const;
         /** \pre Assumes that run() has been called and that \a props.inference == \c MAXPROD
          */
-        std::vector<std::size_t> findMaximum() const;
+        std::vector<std::size_t> findMaximum() const { return dai::findMaximum( *this ); }
         virtual void init() {}
         virtual void init( const VarSet &/*ns*/ ) {}
         virtual Real run();
index da8f343..4b165f5 100644 (file)
@@ -14,7 +14,6 @@
 #include <map>
 #include <set>
 #include <algorithm>
-#include <stack>
 #include <dai/bp.h>
 #include <dai/util.h>
 #include <dai/properties.h>
@@ -481,85 +480,4 @@ void BP::updateResidual( size_t i, size_t _I, Real r ) {
 }
 
 
-std::vector<size_t> BP::findMaximum() const {
-    vector<size_t> maximum( nrVars() );
-    vector<bool> visitedVars( nrVars(), false );
-    vector<bool> visitedFactors( nrFactors(), false );
-    stack<size_t> scheduledFactors;
-    for( size_t i = 0; i < nrVars(); ++i ) {
-        if( visitedVars[i] )
-            continue;
-        visitedVars[i] = true;
-
-        // Maximise with respect to variable i
-        Prob prod;
-        calcBeliefV( i, prod );
-        maximum[i] = prod.argmax().first;
-
-        foreach( const Neighbor &I, nbV(i) )
-            if( !visitedFactors[I] )
-                scheduledFactors.push(I);
-
-        while( !scheduledFactors.empty() ){
-            size_t I = scheduledFactors.top();
-            scheduledFactors.pop();
-            if( visitedFactors[I] )
-                continue;
-            visitedFactors[I] = true;
-
-            // Evaluate if some neighboring variables still need to be fixed; if not, we're done
-            bool allDetermined = true;
-            foreach( const Neighbor &j, nbF(I) )
-                if( !visitedVars[j.node] ) {
-                    allDetermined = false;
-                    break;
-                }
-            if( allDetermined )
-                continue;
-
-            // Calculate product of incoming messages on factor I
-            Prob prod2;
-            calcBeliefF( I, prod2 );
-
-            // The allowed configuration is restrained according to the variables assigned so far:
-            // pick the argmax amongst the allowed states
-            Real maxProb = -numeric_limits<Real>::max();
-            State maxState( factor(I).vars() );
-            for( State s( factor(I).vars() ); s.valid(); ++s ){
-                // First, calculate whether this state is consistent with variables that
-                // have been assigned already
-                bool allowedState = true;
-                foreach( const Neighbor &j, nbF(I) )
-                    if( visitedVars[j.node] && maximum[j.node] != s(var(j.node)) ) {
-                        allowedState = false;
-                        break;
-                    }
-                // If it is consistent, check if its probability is larger than what we have seen so far
-                if( allowedState && prod2[s] > maxProb ) {
-                    maxState = s;
-                    maxProb = prod2[s];
-                }
-            }
-
-            // Decode the argmax
-            foreach( const Neighbor &j, nbF(I) ) {
-                if( visitedVars[j.node] ) {
-                    // We have already visited j earlier - hopefully our state is consistent
-                    if( maximum[j.node] != maxState(var(j.node)) && props.verbose >= 1 )
-                        cerr << "BP::findMaximum - warning: maximum not consistent due to loops." << endl;
-                } else {
-                    // We found a consistent state for variable j
-                    visitedVars[j.node] = true;
-                    maximum[j.node] = maxState( var(j.node) );
-                    foreach( const Neighbor &J, nbV(j) )
-                        if( !visitedFactors[J] )
-                            scheduledFactors.push(J);
-                }
-            }
-        }
-    }
-    return maximum;
-}
-
-
 } // end of namespace dai
index 0e6b218..01ec3d9 100644 (file)
@@ -10,6 +10,7 @@
 
 
 #include <vector>
+#include <stack>
 #include <dai/daialg.h>
 
 
@@ -210,4 +211,83 @@ vector<Factor> calcPairBeliefs( const InfAlg & obj, const VarSet& vs, bool reIni
 }
 
 
+std::vector<size_t> findMaximum( const InfAlg& obj ) {
+    vector<size_t> maximum( obj.fg().nrVars() );
+    vector<bool> visitedVars( obj.fg().nrVars(), false );
+    vector<bool> visitedFactors( obj.fg().nrFactors(), false );
+    stack<size_t> scheduledFactors;
+    for( size_t i = 0; i < obj.fg().nrVars(); ++i ) {
+        if( visitedVars[i] )
+            continue;
+        visitedVars[i] = true;
+
+        // Maximise with respect to variable i
+        Prob prod = obj.beliefV(i).p();
+        maximum[i] = prod.argmax().first;
+
+        foreach( const Neighbor &I, obj.fg().nbV(i) )
+            if( !visitedFactors[I] )
+                scheduledFactors.push(I);
+
+        while( !scheduledFactors.empty() ){
+            size_t I = scheduledFactors.top();
+            scheduledFactors.pop();
+            if( visitedFactors[I] )
+                continue;
+            visitedFactors[I] = true;
+
+            // Evaluate if some neighboring variables still need to be fixed; if not, we're done
+            bool allDetermined = true;
+            foreach( const Neighbor &j, obj.fg().nbF(I) )
+                if( !visitedVars[j.node] ) {
+                    allDetermined = false;
+                    break;
+                }
+            if( allDetermined )
+                continue;
+
+            // Get marginal of factor I
+            Prob probF = obj.beliefF(I).p();
+
+            // The allowed configuration is restrained according to the variables assigned so far:
+            // pick the argmax amongst the allowed states
+            Real maxProb = -numeric_limits<Real>::max();
+            State maxState( obj.fg().factor(I).vars() );
+            for( State s( obj.fg().factor(I).vars() ); s.valid(); ++s ) {
+                // First, calculate whether this state is consistent with variables that
+                // have been assigned already
+                bool allowedState = true;
+                foreach( const Neighbor &j, obj.fg().nbF(I) )
+                    if( visitedVars[j.node] && maximum[j.node] != s(obj.fg().var(j.node)) ) {
+                        allowedState = false;
+                        break;
+                    }
+                // If it is consistent, check if its probability is larger than what we have seen so far
+                if( allowedState && probF[s] > maxProb ) {
+                    maxState = s;
+                    maxProb = probF[s];
+                }
+            }
+
+            // Decode the argmax
+            foreach( const Neighbor &j, obj.fg().nbF(I) ) {
+                if( visitedVars[j.node] ) {
+                    // We have already visited j earlier - hopefully our state is consistent
+                    if( maximum[j.node] != maxState( obj.fg().var(j.node) ) )
+                        DAI_THROWE(RUNTIME_ERROR,"MAP state inconsistent due to loops");
+                } else {
+                    // We found a consistent state for variable j
+                    visitedVars[j.node] = true;
+                    maximum[j.node] = maxState( obj.fg().var(j.node) );
+                    foreach( const Neighbor &J, obj.fg().nbV(j) )
+                        if( !visitedFactors[J] )
+                            scheduledFactors.push(J);
+                }
+            }
+        }
+    }
+    return maximum;
+}
+
+
 } // end of namespace dai
index f0e8a7d..b71e62a 100644 (file)
@@ -10,7 +10,6 @@
 
 
 #include <iostream>
-#include <stack>
 #include <dai/jtree.h>
 
 
@@ -589,41 +588,4 @@ std::pair<size_t,long double> boundTreewidth( const FactorGraph &fg, greedyVaria
 }
 
 
-std::vector<size_t> JTree::findMaximum() const {
-#ifdef DAI_DEBUG
-    // check consistency of variables and factors
-    for( size_t I = 0; I < nrFactors(); I++ ) {
-        size_t linearState = beliefF(I).p().argmax().first;
-        map<Var,size_t> state = calcState( factor(I).vars(), linearState );
-        foreach( const Neighbor& i, nbF(I) ) {
-            if( state[var(i)] != beliefV(i).p().argmax().first ) {
-                cerr << "ERROR: inconsistency between MAP beliefs: factor belief " << beliefF(I) << ": " << state << " is inconsistent with " << beliefV(i) << endl;
-                DAI_THROW(RUNTIME_ERROR);
-            }
-        }
-    }
-    // check consistency of 
-    for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
-        size_t linearStateQa = Qa[alpha].p().argmax().first;
-        map<Var,size_t> stateQa = calcState( Qa[alpha].vars(), linearStateQa );
-        foreach( const Neighbor& beta, nbOR(alpha) ) {
-            size_t linearStateQb = Qb[beta].p().argmax().first;
-            map<Var,size_t> stateQb = calcState( Qb[beta].vars(), linearStateQb );
-            for( map<Var,size_t>::const_iterator it = stateQb.begin(); it != stateQb.end(); it++ )
-                if( stateQa[it->first] != it->second ) {
-                    cerr << "ERROR: inconsistency between MAP beliefs: OR belief " << Qa[alpha] << " (with MAP " << stateQa << ") is inconsistent with IR belief " << Qb[beta] << " (with MAP " << stateQb << ")" << endl;
-                    DAI_THROW(RUNTIME_ERROR);
-                }
-        }
-    }
-#endif
-
-    vector<size_t> maximum;
-    maximum.reserve( nrVars() );
-    for( size_t i = 0; i < nrVars(); i++ )
-        maximum.push_back( beliefV(i).p().argmax().first );
-    return maximum;
-}
-
-
 } // end of namespace dai
diff --git a/tests/jtreemapbug.fg b/tests/jtreemapbug.fg
new file mode 100644 (file)
index 0000000..a64313c
--- /dev/null
@@ -0,0 +1,10 @@
+1
+
+2
+0 1
+2 2
+4
+0      1
+1      2
+2      2
+3      1