src/mf.cpp

   1 /*  Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]
   2     Radboud University Nijmegen, The Netherlands
   3
   4     This file is part of libDAI.
   5
   6     libDAI is free software; you can redistribute it and/or modify
   7     it under the terms of the GNU General Public License as published by
   8     the Free Software Foundation; either version 2 of the License, or
   9     (at your option) any later version.
  10
  11     libDAI is distributed in the hope that it will be useful,
  12     but WITHOUT ANY WARRANTY; without even the implied warranty of
  13     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14     GNU General Public License for more details.
  15
  16     You should have received a copy of the GNU General Public License
  17     along with libDAI; if not, write to the Free Software
  18     Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  19 */
  20
  21
  22 #include <iostream>
  23 #include <sstream>
  24 #include <map>
  25 #include <set>
  26 #include <dai/mf.h>
  27 #include <dai/diffs.h>
  28 #include <dai/util.h>
  29
  30
  31 namespace dai {
  32
  33
  34 using namespace std;
  35
  36
  37 const char *MF::Name = "MF";
  38
  39
  40 void MF::setProperties( const PropertySet &opts ) {
  41     assert( opts.hasKey("tol") );
  42     assert( opts.hasKey("maxiter") );
  43
  44     props.tol = opts.getStringAs<double>("tol");
  45     props.maxiter = opts.getStringAs<size_t>("maxiter");
  46     if( opts.hasKey("verbose") )
  47         props.verbose = opts.getStringAs<size_t>("verbose");
  48     else
  49         props.verbose = 0U;
  50     if( opts.hasKey("damping") )
  51         props.damping = opts.getStringAs<double>("damping");
  52     else
  53         props.damping = 0.0;
  54 }
  55
  56
  57 PropertySet MF::getProperties() const {
  58     PropertySet opts;
  59     opts.Set( "tol", props.tol );
  60     opts.Set( "maxiter", props.maxiter );
  61     opts.Set( "verbose", props.verbose );
  62     opts.Set( "damping", props.damping );
  63     return opts;
  64 }
  65
  66
  67 string MF::printProperties() const {
  68     stringstream s( stringstream::out );
  69     s << "[";
  70     s << "tol=" << props.tol << ",";
  71     s << "maxiter=" << props.maxiter << ",";
  72     s << "verbose=" << props.verbose << ",";
  73     s << "damping=" << props.damping << "]";
  74     return s.str();
  75 }
  76
  77
  78 void MF::construct() {
  79     // create beliefs
  80     _beliefs.clear();
  81     _beliefs.reserve( nrVars() );
  82     for( size_t i = 0; i < nrVars(); ++i )
  83         _beliefs.push_back( Factor( var(i) ) );
  84 }
  85
  86
  87 string MF::identify() const {
  88     return string(Name) + printProperties();
  89 }
  90
  91
  92 void MF::init() {
  93     for( vector<Factor>::iterator qi = _beliefs.begin(); qi != _beliefs.end(); qi++ )
  94         qi->fill(1.0);
  95 }
  96
  97
  98 double MF::run() {
  99     double tic = toc();
 100
 101     if( props.verbose >= 1 )
 102         cout << "Starting " << identify() << "...";
 103
 104     size_t pass_size = _beliefs.size();
 105     Diffs diffs(pass_size * 3, 1.0);
 106
 107     size_t t=0;
 108     for( t=0; t < (props.maxiter*pass_size) && diffs.maxDiff() > props.tol; t++ ) {
 109         // choose random Var i
 110         size_t i = (size_t) (nrVars() * rnd_uniform());
 111
 112         Factor jan;
 113         Factor piet;
 114         foreach( const Neighbor &I, nbV(i) ) {
 115             Factor henk;
 116             foreach( const Neighbor &j, nbF(I) ) // for all j in I \ i
 117                 if( j != i )
 118                     henk *= _beliefs[j];
 119             piet = factor(I).log0();
 120             piet *= henk;
 121             piet = piet.partSum(var(i));
 122             piet = piet.exp();
 123             jan *= piet;
 124         }
 125
 126         jan.normalize();
 127
 128         if( jan.hasNaNs() ) {
 129             cout << Name << "::run():  ERROR: jan has NaNs!" << endl;
 130             return 1.0;
 131         }
 132
 133         if( props.damping != 0.0 )
 134             jan = (jan^(1.0 - props.damping)) * (_beliefs[i]^props.damping);
 135         diffs.push( dist( jan, _beliefs[i], Prob::DISTLINF ) );
 136
 137         _beliefs[i] = jan;
 138     }
 139
 140     _iters = t / pass_size;
 141     if( diffs.maxDiff() > _maxdiff )
 142         _maxdiff = diffs.maxDiff();
 143
 144     if( props.verbose >= 1 ) {
 145         if( diffs.maxDiff() > props.tol ) {
 146             if( props.verbose == 1 )
 147                 cout << endl;
 148             cout << Name << "::run:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << diffs.maxDiff() << endl;
 149         } else {
 150             if( props.verbose >= 2 )
 151                 cout << Name << "::run:  ";
 152             cout << "converged in " << t / pass_size << " passes (" << toc() - tic << " seconds)." << endl;
 153         }
 154     }
 155
 156     return diffs.maxDiff();
 157 }
 158
 159
 160 Factor MF::beliefV( size_t i ) const {
 161     Factor piet;
 162     piet = _beliefs[i];
 163     piet.normalize();
 164     return(piet);
 165 }
 166
 167
 168 Factor MF::belief (const VarSet &ns) const {
 169     if( ns.size() == 1 )
 170         return belief( *(ns.begin()) );
 171     else {
 172         assert( ns.size() == 1 );
 173         return Factor();
 174     }
 175 }
 176
 177
 178 Factor MF::belief (const Var &n) const {
 179     return( beliefV( findVar( n ) ) );
 180 }
 181
 182
 183 vector<Factor> MF::beliefs() const {
 184     vector<Factor> result;
 185     for( size_t i = 0; i < nrVars(); i++ )
 186         result.push_back( beliefV(i) );
 187     return result;
 188 }
 189
 190
 191 Real MF::logZ() const {
 192     Real sum = 0.0;
 193
 194     for(size_t i=0; i < nrVars(); i++ )
 195         sum -= beliefV(i).entropy();
 196     for(size_t I=0; I < nrFactors(); I++ ) {
 197         Factor henk;
 198         foreach( const Neighbor &j, nbF(I) )  // for all j in I
 199             henk *= _beliefs[j];
 200         henk.normalize();
 201         Factor piet;
 202         piet = factor(I).log0();
 203         piet *= henk;
 204         sum -= piet.totalSum();
 205     }
 206
 207     return -sum;
 208 }
 209
 210
 211 void MF::init( const VarSet &ns ) {
 212     for( size_t i = 0; i < nrVars(); i++ ) {
 213         if( ns.contains(var(i) ) )
 214             _beliefs[i].fill( 1.0 );
 215     }
 216 }
 217
 218
 219 } // end of namespace dai