Merge branch 'master' of git.tuebingen.mpg.de:libdai
[libdai.git] / src / mf.cpp
index 526c776..42ede4a 100644 (file)
@@ -1,22 +1,12 @@
-/*  Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]
-    Radboud University Nijmegen, The Netherlands
-    
-    This file is part of libDAI.
-
-    libDAI is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    libDAI is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with libDAI; if not, write to the Free Software
-    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
-*/
+/*  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-2010  Joris Mooij  [joris dot mooij at libdai dot org]
+ *  Copyright (C) 2006-2007  Radboud University Nijmegen, The Netherlands
+ */
 
 
 #include <iostream>
@@ -24,7 +14,6 @@
 #include <map>
 #include <set>
 #include <dai/mf.h>
-#include <dai/diffs.h>
 #include <dai/util.h>
 
 
@@ -34,25 +23,39 @@ namespace dai {
 using namespace std;
 
 
-const char *MF::Name = "MF";
-
-
 void MF::setProperties( const PropertySet &opts ) {
-    assert( opts.hasKey("tol") );
-    assert( opts.hasKey("maxiter") );
-    assert( opts.hasKey("verbose") );
+    DAI_ASSERT( opts.hasKey("tol") );
+    DAI_ASSERT( opts.hasKey("maxiter") );
 
-    props.tol = opts.getStringAs<double>("tol");
+    props.tol = opts.getStringAs<Real>("tol");
     props.maxiter = opts.getStringAs<size_t>("maxiter");
-    props.verbose = opts.getStringAs<size_t>("verbose");
+    if( opts.hasKey("verbose") )
+        props.verbose = opts.getStringAs<size_t>("verbose");
+    else
+        props.verbose = 0U;
+    if( opts.hasKey("damping") )
+        props.damping = opts.getStringAs<Real>("damping");
+    else
+        props.damping = 0.0;
+    if( opts.hasKey("init") )
+        props.init = opts.getStringAs<Properties::InitType>("init");
+    else
+        props.init = Properties::InitType::UNIFORM;
+    if( opts.hasKey("updates") )
+        props.updates = opts.getStringAs<Properties::UpdateType>("updates");
+    else
+        props.updates = Properties::UpdateType::NAIVE;
 }
 
 
 PropertySet MF::getProperties() const {
     PropertySet opts;
-    opts.Set( "tol", props.tol );
-    opts.Set( "maxiter", props.maxiter );
-    opts.Set( "verbose", props.verbose );
+    opts.set( "tol", props.tol );
+    opts.set( "maxiter", props.maxiter );
+    opts.set( "verbose", props.verbose );
+    opts.set( "damping", props.damping );
+    opts.set( "init", props.init );
+    opts.set( "updates", props.updates );
     return opts;
 }
 
@@ -62,115 +65,127 @@ string MF::printProperties() const {
     s << "[";
     s << "tol=" << props.tol << ",";
     s << "maxiter=" << props.maxiter << ",";
-    s << "verbose=" << props.verbose << "]";
+    s << "verbose=" << props.verbose << ",";
+    s << "init=" << props.init << ",";
+    s << "updates=" << props.updates << ",";
+    s << "damping=" << props.damping << "]";
     return s.str();
 }
 
 
-void MF::create() {
-    // clear beliefs
+void MF::construct() {
+    // create beliefs
     _beliefs.clear();
     _beliefs.reserve( nrVars() );
-
-    // create beliefs
     for( size_t i = 0; i < nrVars(); ++i )
-        _beliefs.push_back(Factor(var(i)));
+        _beliefs.push_back( Factor( var(i) ) );
 }
 
 
-string MF::identify() const { 
-    return string(Name) + printProperties();
+void MF::init() {
+    if( props.init == Properties::InitType::UNIFORM )
+        for( size_t i = 0; i < nrVars(); i++ )
+            _beliefs[i].fill( 1.0 );
+    else
+        for( size_t i = 0; i < nrVars(); i++ )
+            _beliefs[i].randomize();
 }
 
 
-void MF::init() {
-    for( vector<Factor>::iterator qi = _beliefs.begin(); qi != _beliefs.end(); qi++ )
-        qi->fill(1.0);
+Factor MF::calcNewBelief( size_t i ) {
+    Factor result;
+    foreach( const Neighbor &I, nbV(i) ) {
+        Factor belief_I_minus_i;
+        foreach( const Neighbor &j, nbF(I) ) // for all j in I \ i
+            if( j != i )
+                belief_I_minus_i *= _beliefs[j];
+        Factor f_I = factor(I);
+        if( props.updates == Properties::UpdateType::NAIVE )
+            f_I.takeLog();
+        Factor msg_I_i = (belief_I_minus_i * f_I).marginal( var(i), false );
+        if( props.updates == Properties::UpdateType::NAIVE )
+            result *= msg_I_i.exp();
+        else
+            result *= msg_I_i;
+    }
+    result.normalize();
+    return result;
 }
 
 
-double MF::run() {
+Real MF::run() {
+    if( props.verbose >= 1 )
+        cerr << "Starting " << identify() << "...";
+
     double tic = toc();
 
-    if( props.verbose >= 1 )
-        cout << "Starting " << identify() << "...";
+    vector<size_t> update_seq;
+    update_seq.reserve( nrVars() );
+    for( size_t i = 0; i < nrVars(); i++ )
+        update_seq.push_back( i );
 
-    size_t pass_size = _beliefs.size();
-    Diffs diffs(pass_size * 3, 1.0);
+    // do several passes over the network until maximum number of iterations has
+    // been reached or until the maximum belief difference is smaller than tolerance
+    Real maxDiff = INFINITY;
+    for( _iters = 0; _iters < props.maxiter && maxDiff > props.tol; _iters++ ) {
+        random_shuffle( update_seq.begin(), update_seq.end() );
 
-    size_t t=0;
-    for( t=0; t < (props.maxiter*pass_size) && diffs.maxDiff() > props.tol; t++ ) {
-        // choose random Var i
-        size_t i = (size_t) (nrVars() * rnd_uniform());
+        maxDiff = -INFINITY;
+        foreach( const size_t &i, update_seq ) {
+            Factor nb = calcNewBelief( i );
 
-        Factor jan;
-        Factor piet;
-        foreach( const Neighbor &I, nbV(i) ) {
-            Factor henk;
-            foreach( const Neighbor &j, nbF(I) ) // for all j in I \ i
-                if( j != i )
-                    henk *= _beliefs[j];
-            piet = factor(I).log0();
-            piet *= henk;
-            piet = piet.partSum(var(i));
-            piet = piet.exp();
-            jan *= piet; 
-        }
+            if( nb.hasNaNs() ) {
+                cerr << name() << "::run():  ERROR: new belief of variable " << var(i) << " has NaNs!" << endl;
+                return 1.0;
+            }
 
-        jan.normalize( Prob::NORMPROB );
+            if( props.damping != 0.0 )
+                nb = (nb^(1.0 - props.damping)) * (_beliefs[i]^props.damping);
 
-        if( jan.hasNaNs() ) {
-            cout << "MF::run():  ERROR: jan has NaNs!" << endl;
-            return 1.0;
+            maxDiff = std::max( maxDiff, dist( nb, _beliefs[i], DISTLINF ) );
+            _beliefs[i] = nb;
         }
 
-        diffs.push( dist( jan, _beliefs[i], Prob::DISTLINF ) );
-
-        _beliefs[i] = jan;
+        if( props.verbose >= 3 )
+            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 )
-                cout << endl;
-            cout << "MF::run:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.maxDiff() << endl;
+                cerr << endl;
+            cerr << name() << "::run:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << maxDiff << endl;
         } else {
-            if( props.verbose >= 2 )
-                cout << "MF::run:  ";
-            cout << "converged in " << t / pass_size << " passes (" << toc() - tic << " clocks)." << endl;
+            if( props.verbose >= 3 )
+                cerr << name() << "::run:  ";
+            cerr << "converged in " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
         }
     }
 
-    return diffs.maxDiff();
+    return maxDiff;
 }
 
 
-Factor MF::beliefV (size_t i) const {
-    Factor piet;
-    piet = _beliefs[i];
-    piet.normalize( Prob::NORMPROB );
-    return(piet);
+Factor MF::beliefV( size_t i ) const {
+    return _beliefs[i].normalized();
 }
 
 
 Factor MF::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 {
-        assert( ns.size() == 1 );
+        DAI_THROW(BELIEF_NOT_AVAILABLE);
         return Factor();
     }
 }
 
 
-Factor MF::belief (const Var &n) const {
-    return( beliefV( findVar( n ) ) );
-}
-
-
 vector<Factor> MF::beliefs() const {
     vector<Factor> result;
     for( size_t i = 0; i < nrVars(); i++ )
@@ -180,30 +195,33 @@ vector<Factor> MF::beliefs() const {
 
 
 Real MF::logZ() const {
-    Real sum = 0.0;
-    
-    for(size_t i=0; i < nrVars(); i++ )
-        sum -= beliefV(i).entropy();
-    for(size_t I=0; I < nrFactors(); I++ ) {
+    Real s = 0.0;
+
+    for( size_t i = 0; i < nrVars(); i++ )
+        s -= beliefV(i).entropy();
+    for( size_t I = 0; I < nrFactors(); I++ ) {
         Factor henk;
         foreach( const Neighbor &j, nbF(I) )  // for all j in I
             henk *= _beliefs[j];
-        henk.normalize( Prob::NORMPROB );
+        henk.normalize();
         Factor piet;
-        piet = factor(I).log0();
+        piet = factor(I).log(true);
         piet *= henk;
-        sum -= piet.totalSum();
+        s -= piet.sum();
     }
 
-    return -sum;
+    return -s;
 }
 
 
 void MF::init( const VarSet &ns ) {
-    for( size_t i = 0; i < nrVars(); i++ ) {
-        if( ns.contains(var(i) ) )
-            _beliefs[i].fill( 1.0 );
-    }
+    for( size_t i = 0; i < nrVars(); i++ )
+        if( ns.contains(var(i) ) ) {
+            if( props.init == Properties::InitType::UNIFORM )
+                _beliefs[i].fill( 1.0 );
+            else
+                _beliefs[i].randomize();
+        }
 }