Fixed bug in findMaximum(): inconsistent max-beliefs are now detected,

instead of returning a MAP state with zero joint probability
(reported by Hynek Urban)

diff --git a/ChangeLog b/ChangeLog
index 757aec7..ce5a257 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,6 @@
+* Fixed bug in findMaximum(): inconsistent max-beliefs are now detected,
+  instead of returning a MAP state with zero joint probability
+  (reported by Hynek Urban)
 * Fixed a Boost-related bug in src/util.cpp (reported by Avneesh Saluja)
   (the random seed needs to be an unsigned int on some platforms)
 * Fixed two bugs in examples/example_sprinkler_gibbs (reported by Priya)

diff --git a/examples/example.cpp b/examples/example.cpp
index 86cf473..ac6efd0 100644 (file)
--- a/examples/example.cpp
+++ b/examples/example.cpp
@@ -12,6 +12,9 @@
 #include <iostream>
 #include <map>
 #include <dai/alldai.h>  // Include main libDAI header file
+#include <dai/jtree.h>
+#include <dai/bp.h>
+#include <dai/decmap.h>
 
 
 using namespace dai;

diff --git a/include/dai/bp.h b/include/dai/bp.h
index 8d35bf6..8e8d38d 100644 (file)
--- a/include/dai/bp.h
+++ b/include/dai/bp.h
@@ -6,6 +6,7 @@
  *
  *  Copyright (C) 2006-2010  Joris Mooij  [joris dot mooij at libdai dot org]
  *  Copyright (C) 2006-2007  Radboud University Nijmegen, The Netherlands
+ *  Copyright (C) 2008-2009  Giuseppe Passino
  */
 
 

diff --git a/include/dai/factorgraph.h b/include/dai/factorgraph.h
index 47195ca..f15e4d4 100644 (file)
--- a/include/dai/factorgraph.h
+++ b/include/dai/factorgraph.h
@@ -231,7 +231,7 @@ class FactorGraph {
         std::vector<VarSet> maximalFactorDomains() const;
 
         /// Evaluates the log score (i.e., minus the energy) of the joint configuration \a statevec
-        Real logScore( const std::vector<size_t>& statevec );
+        Real logScore( const std::vector<size_t>& statevec ) const;
     //@}
 
     /// \name Backup/restore mechanism for factors

diff --git a/include/dai/jtree.h b/include/dai/jtree.h
index ad2dec3..86b898e 100644 (file)
--- a/include/dai/jtree.h
+++ b/include/dai/jtree.h
@@ -127,7 +127,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 { return dai::findMaximum( *this ); }
+        std::vector<std::size_t> findMaximum() const;
         virtual void init() {}
         virtual void init( const VarSet &/*ns*/ ) {}
         virtual Real run();

diff --git a/include/dai/prob.h b/include/dai/prob.h
index 633cfa4..41e33d5 100644 (file)
--- a/include/dai/prob.h
+++ b/include/dai/prob.h
@@ -6,6 +6,7 @@
  *
  *  Copyright (C) 2006-2009  Joris Mooij  [joris dot mooij at libdai dot org]
  *  Copyright (C) 2006-2007  Radboud University Nijmegen, The Netherlands
+ *  Copyright (C) 2008-2009  Giuseppe Passino
  */
 
 

diff --git a/src/bp.cpp b/src/bp.cpp
index 5944d62..a9a8c20 100644 (file)
--- a/src/bp.cpp
+++ b/src/bp.cpp
@@ -6,6 +6,8 @@
  *
  *  Copyright (C) 2006-2010  Joris Mooij  [joris dot mooij at libdai dot org]
  *  Copyright (C) 2006-2007  Radboud University Nijmegen, The Netherlands
+ *  Copyright (C) 2008-2009  Giuseppe Passino
+ *  Copyright (C) 2008       Claudio Lima
  */
 
 

diff --git a/src/daialg.cpp b/src/daialg.cpp
index 01ec3d9..1dc10b4 100644 (file)
--- a/src/daialg.cpp
+++ b/src/daialg.cpp
@@ -6,6 +6,7 @@
  *
  *  Copyright (C) 2006-2009  Joris Mooij  [joris dot mooij at libdai dot org]
  *  Copyright (C) 2006-2007  Radboud University Nijmegen, The Netherlands
+ *  Copyright (C) 2008-2009  Giuseppe Passino
  */
 
 
@@ -216,73 +217,57 @@ std::vector<size_t> findMaximum( const InfAlg& obj ) {
     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] )
+    scheduledFactors.push( 0 );
+    while( !scheduledFactors.empty() ) {
+        size_t I = scheduledFactors.top();
+        scheduledFactors.pop();
+        if( visitedFactors[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;
+        visitedFactors[I] = true;
+
+        // 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() );
+        size_t maxcount = 0;
+        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] ) {
-                    allDetermined = false;
+                if( visitedVars[j.node] && maximum[j.node] != s(obj.fg().var(j.node)) ) {
+                    allowedState = 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 ) {
+            // If it is consistent, check if its probability is larger than what we have seen so far
+            if( allowedState ) {
+                if( probF[s] > maxProb ) {
                     maxState = s;
                     maxProb = probF[s];
-                }
+                    maxcount = 1;
+                } else
+                    maxcount++;
             }
-
-            // 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);
-                }
+        }
+        DAI_ASSERT( maxProb != 0.0 );
+        DAI_ASSERT( obj.fg().factor(I).p()[maxState] != 0.0 );
+
+        // 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);
             }
         }
     }

diff --git a/src/factorgraph.cpp b/src/factorgraph.cpp
index 416ecf8..7238656 100644 (file)
--- a/src/factorgraph.cpp
+++ b/src/factorgraph.cpp
@@ -6,6 +6,7 @@
  *
  *  Copyright (C) 2006-2009  Joris Mooij  [joris dot mooij at libdai dot org]
  *  Copyright (C) 2006-2007  Radboud University Nijmegen, The Netherlands
+ *  Copyright (C) 2008       Christian Wojek
  */
 
 
@@ -166,6 +167,8 @@ std::istream& operator>> ( std::istream& is, FactorGraph &fg ) {
             Ivars.push_back( Var(labels[mi], dims[mi]) );
         }
         facs.push_back( Factor( VarSet( Ivars.begin(), Ivars.end(), Ivars.size() ), (Real)0 ) );
+        if( verbose >= 3 )
+            cerr << "  vardims: " << vardims << endl;
 
         // calculate permutation object
         Permute permindex( Ivars );
@@ -337,7 +340,7 @@ vector<VarSet> FactorGraph::maximalFactorDomains() const {
 }
 
 
-Real FactorGraph::logScore( const std::vector<size_t>& statevec ) {
+Real FactorGraph::logScore( const std::vector<size_t>& statevec ) const {
     // Construct a State object that represents statevec
     // This decouples the representation of the joint state in statevec from the factor graph
     map<Var, size_t> statemap;

diff --git a/src/jtree.cpp b/src/jtree.cpp
index d2253bb..14f2e63 100644 (file)
--- a/src/jtree.cpp
+++ b/src/jtree.cpp
@@ -10,6 +10,7 @@
 
 
 #include <iostream>
+#include <stack>
 #include <dai/jtree.h>
 
 
@@ -580,4 +581,71 @@ std::pair<size_t,long double> boundTreewidth( const FactorGraph &fg, greedyVaria
 }
 
 
+std::vector<size_t> JTree::findMaximum() const {
+    vector<size_t> maximum( nrVars() );
+    vector<bool> visitedVars( nrVars(), false );
+    vector<bool> visitedORs( nrORs(), false );
+    stack<size_t> scheduledORs;
+    scheduledORs.push( 0 );
+    while( !scheduledORs.empty() ) {
+        size_t alpha = scheduledORs.top();
+        scheduledORs.pop();
+        if( visitedORs[alpha] )
+            continue;
+        visitedORs[alpha] = true;
+
+        // Get marginal of outer region alpha 
+        Prob probF = Qa[alpha].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( OR(alpha).vars() );
+        size_t maxcount = 0;
+        for( State s( OR(alpha).vars() ); s.valid(); ++s ) {
+            // First, calculate whether this state is consistent with variables that
+            // have been assigned already
+            bool allowedState = true;
+            foreach( const Var& j, OR(alpha).vars() ) {
+                size_t j_index = findVar(j);
+                if( visitedVars[j_index] && maximum[j_index] != s(j) ) {
+                    allowedState = false;
+                    break;
+                }
+            }
+            // If it is consistent, check if its probability is larger than what we have seen so far
+            if( allowedState ) {
+                if( probF[s] > maxProb ) {
+                    maxState = s;
+                    maxProb = probF[s];
+                    maxcount = 1;
+                } else
+                    maxcount++;
+            }
+        }
+        DAI_ASSERT( maxProb != 0.0 );
+        DAI_ASSERT( Qa[alpha][maxState] != 0.0 );
+
+        // Decode the argmax
+        foreach( const Var& j, OR(alpha).vars() ) {
+            size_t j_index = findVar(j);
+            if( visitedVars[j_index] ) {
+                // We have already visited j earlier - hopefully our state is consistent
+                if( maximum[j_index] != maxState( j ) )
+                    DAI_THROWE(RUNTIME_ERROR,"MAP state inconsistent due to loops");
+            } else {
+                // We found a consistent state for variable j
+                visitedVars[j_index] = true;
+                maximum[j_index] = maxState( j );
+                foreach( const Neighbor &beta, nbOR(alpha) )
+                    foreach( const Neighbor &alpha2, nbIR(beta) )
+                        if( !visitedORs[alpha2] )
+                            scheduledORs.push(alpha2);
+            }
+        }
+    }
+    return maximum;
+}
+
+
 } // end of namespace dai

diff --git a/src/matlab/dai.cpp b/src/matlab/dai.cpp
index 0289992..258310c 100644 (file)
--- a/src/matlab/dai.cpp
+++ b/src/matlab/dai.cpp
@@ -6,6 +6,7 @@
  *
  *  Copyright (C) 2006-2010  Joris Mooij  [joris dot mooij at libdai dot org]
  *  Copyright (C) 2006-2007  Radboud University Nijmegen, The Netherlands
+ *  Copyright (C) 2009       Sebastian Nowozin
  */