src/lc.cpp

   1 /*  This file is part of libDAI - http://www.libdai.org/
   2  *
   3  *  libDAI is licensed under the terms of the GNU General Public License version
   4  *  2, or (at your option) any later version. libDAI is distributed without any
   5  *  warranty. See the file COPYING for more details.
   6  *
   7  *  Copyright (C) 2006-2009  Joris Mooij  [joris dot mooij at libdai dot org]
   8  *  Copyright (C) 2006-2007  Radboud University Nijmegen, The Netherlands
   9  */
  10
  11
  12 #include <iostream>
  13 #include <algorithm>
  14 #include <map>
  15 #include <set>
  16 #include <dai/lc.h>
  17 #include <dai/util.h>
  18 #include <dai/alldai.h>
  19
  20
  21 namespace dai {
  22
  23
  24 using namespace std;
  25
  26
  27 const char *LC::Name = "LC";
  28
  29
  30 void LC::setProperties( const PropertySet &opts ) {
  31     DAI_ASSERT( opts.hasKey("tol") );
  32     DAI_ASSERT( opts.hasKey("maxiter") );
  33     DAI_ASSERT( opts.hasKey("verbose") );
  34     DAI_ASSERT( opts.hasKey("cavity") );
  35     DAI_ASSERT( opts.hasKey("updates") );
  36
  37     props.tol = opts.getStringAs<Real>("tol");
  38     props.maxiter = opts.getStringAs<size_t>("maxiter");
  39     props.verbose = opts.getStringAs<size_t>("verbose");
  40     props.cavity = opts.getStringAs<Properties::CavityType>("cavity");
  41     props.updates = opts.getStringAs<Properties::UpdateType>("updates");
  42     if( opts.hasKey("cavainame") )
  43         props.cavainame = opts.getStringAs<string>("cavainame");
  44     if( opts.hasKey("cavaiopts") )
  45         props.cavaiopts = opts.getStringAs<PropertySet>("cavaiopts");
  46     if( opts.hasKey("reinit") )
  47         props.reinit = opts.getStringAs<bool>("reinit");
  48     if( opts.hasKey("damping") )
  49         props.damping = opts.getStringAs<Real>("damping");
  50     else
  51         props.damping = 0.0;
  52 }
  53
  54
  55 PropertySet LC::getProperties() const {
  56     PropertySet opts;
  57     opts.Set( "tol", props.tol );
  58     opts.Set( "maxiter", props.maxiter );
  59     opts.Set( "verbose", props.verbose );
  60     opts.Set( "cavity", props.cavity );
  61     opts.Set( "updates", props.updates );
  62     opts.Set( "cavainame", props.cavainame );
  63     opts.Set( "cavaiopts", props.cavaiopts );
  64     opts.Set( "reinit", props.reinit );
  65     opts.Set( "damping", props.damping );
  66     return opts;
  67 }
  68
  69
  70 string LC::printProperties() const {
  71     stringstream s( stringstream::out );
  72     s << "[";
  73     s << "tol=" << props.tol << ",";
  74     s << "maxiter=" << props.maxiter << ",";
  75     s << "verbose=" << props.verbose << ",";
  76     s << "cavity=" << props.cavity << ",";
  77     s << "updates=" << props.updates << ",";
  78     s << "cavainame=" << props.cavainame << ",";
  79     s << "cavaiopts=" << props.cavaiopts << ",";
  80     s << "reinit=" << props.reinit << ",";
  81     s << "damping=" << props.damping << "]";
  82     return s.str();
  83 }
  84
  85
  86 LC::LC( const FactorGraph & fg, const PropertySet &opts ) : DAIAlgFG(fg), _pancakes(), _cavitydists(), _phis(), _beliefs(), _maxdiff(0.0), _iters(0), props() {
  87     setProperties( opts );
  88
  89     // create pancakes
  90     _pancakes.resize( nrVars() );
  91
  92     // create cavitydists
  93     for( size_t i=0; i < nrVars(); i++ )
  94         _cavitydists.push_back(Factor( delta(i) ));
  95
  96     // create phis
  97     _phis.reserve( nrVars() );
  98     for( size_t i = 0; i < nrVars(); i++ ) {
  99         _phis.push_back( vector<Factor>() );
 100         _phis[i].reserve( nbV(i).size() );
 101         foreach( const Neighbor &I, nbV(i) )
 102             _phis[i].push_back( Factor( factor(I).vars() / var(i) ) );
 103     }
 104
 105     // create beliefs
 106     _beliefs.reserve( nrVars() );
 107     for( size_t i=0; i < nrVars(); i++ )
 108         _beliefs.push_back(Factor(var(i)));
 109 }
 110
 111
 112 string LC::identify() const {
 113     return string(Name) + printProperties();
 114 }
 115
 116
 117 void LC::CalcBelief (size_t i) {
 118     _beliefs[i] = _pancakes[i].marginal(var(i));
 119 }
 120
 121
 122 Real LC::CalcCavityDist (size_t i, const std::string &name, const PropertySet &opts) {
 123     Factor Bi;
 124     Real maxdiff = 0;
 125
 126     if( props.verbose >= 2 )
 127         cerr << "Initing cavity " << var(i) << "(" << delta(i).size() << " vars, " << delta(i).nrStates() << " states)" << endl;
 128
 129     if( props.cavity == Properties::CavityType::UNIFORM )
 130         Bi = Factor(delta(i));
 131     else {
 132         InfAlg *cav = newInfAlg( name, *this, opts );
 133         cav->makeCavity( i );
 134
 135         if( props.cavity == Properties::CavityType::FULL )
 136             Bi = calcMarginal( *cav, cav->fg().delta(i), props.reinit );
 137         else if( props.cavity == Properties::CavityType::PAIR ) {
 138             vector<Factor> pairbeliefs = calcPairBeliefs( *cav, cav->fg().delta(i), props.reinit, false );
 139             for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
 140                 Bi *= pairbeliefs[ij];
 141         } else if( props.cavity == Properties::CavityType::PAIR2 ) {
 142             vector<Factor> pairbeliefs = calcPairBeliefs( *cav, cav->fg().delta(i), props.reinit, true );
 143             for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
 144                 Bi *= pairbeliefs[ij];
 145         }
 146         maxdiff = cav->maxDiff();
 147         delete cav;
 148     }
 149     Bi.normalize();
 150     _cavitydists[i] = Bi;
 151
 152     return maxdiff;
 153 }
 154
 155
 156 Real LC::InitCavityDists( const std::string &name, const PropertySet &opts ) {
 157     double tic = toc();
 158
 159     if( props.verbose >= 1 ) {
 160         cerr << Name << "::InitCavityDists:  ";
 161         if( props.cavity == Properties::CavityType::UNIFORM )
 162             cerr << "Using uniform initial cavity distributions" << endl;
 163         else if( props.cavity == Properties::CavityType::FULL )
 164             cerr << "Using full " << name << opts << "...";
 165         else if( props.cavity == Properties::CavityType::PAIR )
 166             cerr << "Using pairwise " << name << opts << "...";
 167         else if( props.cavity == Properties::CavityType::PAIR2 )
 168             cerr << "Using pairwise(new) " << name << opts << "...";
 169     }
 170
 171     Real maxdiff = 0.0;
 172     for( size_t i = 0; i < nrVars(); i++ ) {
 173         Real md = CalcCavityDist(i, name, opts);
 174         if( md > maxdiff )
 175             maxdiff = md;
 176     }
 177
 178     if( props.verbose >= 1 ) {
 179         cerr << Name << "::InitCavityDists used " << toc() - tic << " seconds." << endl;
 180     }
 181
 182     return maxdiff;
 183 }
 184
 185
 186 long LC::SetCavityDists( std::vector<Factor> &Q ) {
 187     if( props.verbose >= 1 )
 188         cerr << Name << "::SetCavityDists:  Setting initial cavity distributions" << endl;
 189     if( Q.size() != nrVars() )
 190         return -1;
 191     for( size_t i = 0; i < nrVars(); i++ ) {
 192         if( _cavitydists[i].vars() != Q[i].vars() ) {
 193             return i+1;
 194         } else
 195             _cavitydists[i] = Q[i];
 196     }
 197     return 0;
 198 }
 199
 200
 201 void LC::init() {
 202     for( size_t i = 0; i < nrVars(); ++i )
 203         foreach( const Neighbor &I, nbV(i) )
 204             if( props.updates == Properties::UpdateType::SEQRND )
 205                 _phis[i][I.iter].randomize();
 206             else
 207                 _phis[i][I.iter].fill(1.0);
 208 }
 209
 210
 211 Factor LC::NewPancake (size_t i, size_t _I, bool & hasNaNs) {
 212     size_t I = nbV(i)[_I];
 213     Factor piet = _pancakes[i];
 214
 215     // recalculate _pancake[i]
 216     VarSet Ivars = factor(I).vars();
 217     Factor A_I;
 218     for( VarSet::const_iterator k = Ivars.begin(); k != Ivars.end(); k++ )
 219         if( var(i) != *k )
 220             A_I *= (_pancakes[findVar(*k)] * factor(I).inverse()).marginal( Ivars / var(i), false );
 221     if( Ivars.size() > 1 )
 222         A_I ^= (1.0 / (Ivars.size() - 1));
 223     Factor A_Ii = (_pancakes[i] * factor(I).inverse() * _phis[i][_I].inverse()).marginal( Ivars / var(i), false );
 224     Factor quot = A_I / A_Ii;
 225     if( props.damping != 0.0 )
 226         quot = (quot^(1.0 - props.damping)) * (_phis[i][_I]^props.damping);
 227
 228     piet *= quot / _phis[i][_I].normalized();
 229     _phis[i][_I] = quot.normalized();
 230
 231     piet.normalize();
 232
 233     if( piet.hasNaNs() ) {
 234         cerr << Name << "::NewPancake(" << i << ", " << _I << "):  has NaNs!" << endl;
 235         hasNaNs = true;
 236     }
 237
 238     return piet;
 239 }
 240
 241
 242 Real LC::run() {
 243     if( props.verbose >= 1 )
 244         cerr << "Starting " << identify() << "...";
 245     if( props.verbose >= 2 )
 246         cerr << endl;
 247
 248     double tic = toc();
 249     vector<Real> diffs( nrVars(), INFINITY );
 250     Real maxDiff = INFINITY;
 251
 252     Real md = InitCavityDists( props.cavainame, props.cavaiopts );
 253     if( md > _maxdiff )
 254         _maxdiff = md;
 255
 256     for( size_t i = 0; i < nrVars(); i++ ) {
 257         _pancakes[i] = _cavitydists[i];
 258
 259         foreach( const Neighbor &I, nbV(i) ) {
 260             _pancakes[i] *= factor(I);
 261             if( props.updates == Properties::UpdateType::SEQRND )
 262               _pancakes[i] *= _phis[i][I.iter];
 263         }
 264
 265         _pancakes[i].normalize();
 266
 267         CalcBelief(i);
 268     }
 269
 270     vector<Factor> old_beliefs;
 271     for(size_t i=0; i < nrVars(); i++ )
 272         old_beliefs.push_back(beliefV(i));
 273
 274     bool hasNaNs = false;
 275     for( size_t i=0; i < nrVars(); i++ )
 276         if( _pancakes[i].hasNaNs() ) {
 277             hasNaNs = true;
 278             break;
 279         }
 280     if( hasNaNs ) {
 281         cerr << Name << "::run:  initial _pancakes has NaNs!" << endl;
 282         return 1.0;
 283     }
 284
 285     size_t nredges = nrEdges();
 286     vector<Edge> update_seq;
 287     update_seq.reserve( nredges );
 288     for( size_t i = 0; i < nrVars(); ++i )
 289         foreach( const Neighbor &I, nbV(i) )
 290             update_seq.push_back( Edge( i, I.iter ) );
 291
 292     // do several passes over the network until maximum number of iterations has
 293     // been reached or until the maximum belief difference is smaller than tolerance
 294     for( _iters=0; _iters < props.maxiter && maxDiff > props.tol; _iters++ ) {
 295         // Sequential updates
 296         if( props.updates == Properties::UpdateType::SEQRND )
 297             random_shuffle( update_seq.begin(), update_seq.end() );
 298
 299         for( size_t t=0; t < nredges; t++ ) {
 300             size_t i = update_seq[t].first;
 301             size_t _I = update_seq[t].second;
 302             _pancakes[i] = NewPancake( i, _I, hasNaNs);
 303             if( hasNaNs )
 304                 return 1.0;
 305             CalcBelief( i );
 306         }
 307
 308         // compare new beliefs with old ones
 309         for(size_t i=0; i < nrVars(); i++ ) {
 310             diffs[i] = dist( beliefV(i), old_beliefs[i], Prob::DISTLINF );
 311             old_beliefs[i] = beliefV(i);
 312         }
 313         maxDiff = max( diffs );
 314
 315         if( props.verbose >= 3 )
 316             cerr << Name << "::run:  maxdiff " << maxDiff << " after " << _iters+1 << " passes" << endl;
 317     }
 318
 319     if( maxDiff > _maxdiff )
 320         _maxdiff = maxDiff;
 321
 322     if( props.verbose >= 1 ) {
 323         if( maxDiff > props.tol ) {
 324             if( props.verbose == 1 )
 325                 cerr << endl;
 326                 cerr << Name << "::run:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << maxDiff << endl;
 327         } else {
 328             if( props.verbose >= 2 )
 329                 cerr << Name << "::run:  ";
 330                 cerr << "converged in " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
 331         }
 332     }
 333
 334     return maxDiff;
 335 }
 336
 337
 338 } // end of namespace dai