Miscellaneous smaller improvements
[libdai.git] / src / bp.cpp
index e3043cc..5d4a79a 100644 (file)
@@ -38,7 +38,7 @@ void BP::setProperties( const PropertySet &opts ) {
     DAI_ASSERT( opts.hasKey("logdomain") );
     DAI_ASSERT( opts.hasKey("updates") );
 
-    props.tol = opts.getStringAs<double>("tol");
+    props.tol = opts.getStringAs<Real>("tol");
     props.maxiter = opts.getStringAs<size_t>("maxiter");
     props.logdomain = opts.getStringAs<bool>("logdomain");
     props.updates = opts.getStringAs<Properties::UpdateType>("updates");
@@ -48,7 +48,7 @@ void BP::setProperties( const PropertySet &opts ) {
     else
         props.verbose = 0;
     if( opts.hasKey("damping") )
-        props.damping = opts.getStringAs<double>("damping");
+        props.damping = opts.getStringAs<Real>("damping");
     else
         props.damping = 0.0;
     if( opts.hasKey("inference") )
@@ -120,7 +120,7 @@ void BP::construct() {
 
 
 void BP::init() {
-    double c = props.logdomain ? 0.0 : 1.0;
+    Real c = props.logdomain ? 0.0 : 1.0;
     for( size_t i = 0; i < nrVars(); ++i ) {
         foreach( const Neighbor &I, nbV(i) ) {
             message( i, I.iter ).fill( c );
@@ -225,41 +225,40 @@ void BP::calcNewMessage( size_t i, size_t _I ) {
 
 // BP::run does not check for NANs for performance reasons
 // Somehow NaNs do not often occur in BP...
-double BP::run() {
+Real BP::run() {
     if( props.verbose >= 1 )
         cerr << "Starting " << identify() << "...";
     if( props.verbose >= 3)
         cerr << endl;
 
     double tic = toc();
-    Diffs diffs(nrVars(), 1.0);
+    Real maxDiff = INFINITY;
 
-    vector<Edge> update_seq;
-
-    vector<Factor> old_beliefs;
-    old_beliefs.reserve( nrVars() );
+    vector<Factor> oldBeliefsV, oldBeliefsF;
+    oldBeliefsV.reserve( nrVars() );
     for( size_t i = 0; i < nrVars(); ++i )
-        old_beliefs.push_back( beliefV(i) );
+        oldBeliefsV.push_back( beliefV(i) );
+    oldBeliefsF.reserve( nrFactors() );
+    for( size_t I = 0; I < nrFactors(); ++I )
+        oldBeliefsF.push_back( beliefF(I) );
 
     size_t nredges = nrEdges();
-
+    vector<Edge> update_seq;
     if( props.updates == Properties::UpdateType::SEQMAX ) {
         // do the first pass
         for( size_t i = 0; i < nrVars(); ++i )
-            foreach( const Neighbor &I, nbV(i) ) {
+            foreach( const Neighbor &I, nbV(i) )
                 calcNewMessage( i, I.iter );
-            }
     } else {
         update_seq.reserve( nredges );
-        /// \todo Investigate whether performance increases by switching the order of following two loops:
-        for( size_t i = 0; i < nrVars(); ++i )
-            foreach( const Neighbor &I, nbV(i) )
-                update_seq.push_back( Edge( i, I.iter ) );
+        for( size_t I = 0; I < nrFactors(); I++ )
+            foreach( const Neighbor &i, nbF(I) )
+                update_seq.push_back( Edge( i, i.dual ) );
     }
 
     // do several passes over the network until maximum number of iterations has
     // been reached or until the maximum belief difference is smaller than tolerance
-    for( _iters=0; _iters < props.maxiter && diffs.maxDiff() > props.tol; ++_iters ) {
+    for( _iters=0; _iters < props.maxiter && maxDiff > props.tol; ++_iters ) {
         if( props.updates == Properties::UpdateType::SEQMAX ) {
             // Residuals-BP by Koller et al.
             for( size_t t = 0; t < nredges; ++t ) {
@@ -301,24 +300,30 @@ double BP::run() {
         }
 
         // calculate new beliefs and compare with old ones
+        maxDiff = -INFINITY;
         for( size_t i = 0; i < nrVars(); ++i ) {
-            Factor nb( beliefV(i) );
-            diffs.push( dist( nb, old_beliefs[i], Prob::DISTLINF ) );
-            old_beliefs[i] = nb;
+            Factor b( beliefV(i) );
+            maxDiff = std::max( maxDiff, dist( b, oldBeliefsV[i], Prob::DISTLINF ) );
+            oldBeliefsV[i] = b;
+        }
+        for( size_t I = 0; I < nrFactors(); ++I ) {
+            Factor b( beliefF(I) );
+            maxDiff = std::max( maxDiff, dist( b, oldBeliefsF[I], Prob::DISTLINF ) );
+            oldBeliefsF[I] = b;
         }
 
         if( props.verbose >= 3 )
-            cerr << Name << "::run:  maxdiff " << diffs.maxDiff() << " after " << _iters+1 << " passes" << endl;
+            cerr << Name << "::run:  maxdiff " << maxDiff << " after " << _iters+1 << " passes" << endl;
     }
 
-    if( diffs.maxDiff() > _maxdiff )
-        _maxdiff = diffs.maxDiff();
+    if( maxDiff > _maxdiff )
+        _maxdiff = maxDiff;
 
     if( props.verbose >= 1 ) {
-        if( diffs.maxDiff() > props.tol ) {
+        if( maxDiff > props.tol ) {
             if( props.verbose == 1 )
                 cerr << endl;
-                cerr << Name << "::run:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << diffs.maxDiff() << endl;
+                cerr << Name << "::run:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << maxDiff << endl;
         } else {
             if( props.verbose >= 3 )
                 cerr << Name << "::run:  ";
@@ -326,7 +331,7 @@ double BP::run() {
         }
     }
 
-    return diffs.maxDiff();
+    return maxDiff;
 }
 
 
@@ -412,11 +417,6 @@ Factor BP::beliefF( size_t I ) const {
 }
 
 
-Factor BP::belief( const Var &n ) const {
-    return( beliefV( findVar( n ) ) );
-}
-
-
 vector<Factor> BP::beliefs() const {
     vector<Factor> result;
     for( size_t i = 0; i < nrVars(); ++i )
@@ -428,14 +428,17 @@ vector<Factor> BP::beliefs() const {
 
 
 Factor BP::belief( const VarSet &ns ) const {
-    if( ns.size() == 1 )
-        return belief( *(ns.begin()) );
+    if( ns.size() == 0 )
+        return Factor();
+    else if( ns.size() == 1 )
+        return beliefV( findVar( *(ns.begin() ) ) );
     else {
         size_t I;
         for( I = 0; I < nrFactors(); I++ )
             if( factor(I).vars() >> ns )
                 break;
-        DAI_ASSERT( I != nrFactors() );
+        if( I == nrFactors() )
+            DAI_THROW(BELIEF_NOT_AVAILABLE);
         return beliefF(I).marginal(ns);
     }
 }
@@ -443,7 +446,7 @@ Factor BP::belief( const VarSet &ns ) const {
 
 Real BP::logZ() const {
     Real sum = 0.0;
-    for(size_t i = 0; i < nrVars(); ++i )
+    for( size_t i = 0; i < nrVars(); ++i )
         sum += (1.0 - nbV(i).size()) * beliefV(i).entropy();
     for( size_t I = 0; I < nrFactors(); ++I )
         sum -= dist( beliefF(I), factor(I), Prob::DISTKL );
@@ -460,7 +463,7 @@ void BP::init( const VarSet &ns ) {
     for( VarSet::const_iterator n = ns.begin(); n != ns.end(); ++n ) {
         size_t ni = findVar( *n );
         foreach( const Neighbor &I, nbV( ni ) ) {
-            double val = props.logdomain ? 0.0 : 1.0;
+            Real val = props.logdomain ? 0.0 : 1.0;
             message( ni, I.iter ).fill( val );
             newMessage( ni, I.iter ).fill( val );
             if( props.updates == Properties::UpdateType::SEQMAX )
@@ -478,14 +481,17 @@ void BP::updateMessage( size_t i, size_t _I ) {
         if( props.updates == Properties::UpdateType::SEQMAX )
             updateResidual( i, _I, 0.0 );
     } else {
-        message(i,_I) = (message(i,_I) ^ props.damping) * (newMessage(i,_I) ^ (1.0 - props.damping));
+        if( props.logdomain )
+            message(i,_I) = (message(i,_I) * props.damping) + (newMessage(i,_I) * (1.0 - props.damping));
+        else
+            message(i,_I) = (message(i,_I) ^ props.damping) * (newMessage(i,_I) ^ (1.0 - props.damping));
         if( props.updates == Properties::UpdateType::SEQMAX )
             updateResidual( i, _I, dist( newMessage(i,_I), message(i,_I), Prob::DISTLINF ) );
     }
 }
 
 
-void BP::updateResidual( size_t i, size_t _I, double r ) {
+void BP::updateResidual( size_t i, size_t _I, Real r ) {
     EdgeProp* pEdge = &_edges[i][_I];
     pEdge->residual = r;
 
@@ -508,7 +514,7 @@ std::vector<size_t> BP::findMaximum() const {
         // Maximise with respect to variable i
         Prob prod;
         calcBeliefV( i, prod );
-        maximum[i] = max_element( prod.begin(), prod.end() ) - prod.begin();
+        maximum[i] = prod.argmax().first;
 
         foreach( const Neighbor &I, nbV(i) )
             if( !visitedFactors[I] )