src/lc.cpp

   1 /*  This file is part of libDAI - http://www.libdai.org/
   2  *
   3  *  Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
   4  *
   5  *  Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
   6  */
   7
   8
   9 #include <dai/dai_config.h>
  10 #ifdef DAI_WITH_LC
  11
  12
  13 #include <iostream>
  14 #include <algorithm>
  15 #include <map>
  16 #include <set>
  17 #include <dai/lc.h>
  18 #include <dai/util.h>
  19 #include <dai/alldai.h>
  20
  21
  22 namespace dai {
  23
  24
  25 using namespace std;
  26
  27
  28 void LC::setProperties( const PropertySet &opts ) {
  29     DAI_ASSERT( opts.hasKey("tol") );
  30     DAI_ASSERT( opts.hasKey("maxiter") );
  31     DAI_ASSERT( opts.hasKey("cavity") );
  32     DAI_ASSERT( opts.hasKey("updates") );
  33
  34     props.tol = opts.getStringAs<Real>("tol");
  35     props.maxiter = opts.getStringAs<size_t>("maxiter");
  36     props.cavity = opts.getStringAs<Properties::CavityType>("cavity");
  37     props.updates = opts.getStringAs<Properties::UpdateType>("updates");
  38     if( opts.hasKey("verbose") )
  39         props.verbose = opts.getStringAs<size_t>("verbose");
  40     else
  41         props.verbose = 0;
  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         bforeach( 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 void LC::CalcBelief (size_t i) {
 113     _beliefs[i] = _pancakes[i].marginal(var(i));
 114 }
 115
 116
 117 Factor LC::belief (const VarSet &ns) const {
 118     if( ns.size() == 0 )
 119         return Factor();
 120     else if( ns.size() == 1 )
 121         return beliefV( findVar( *(ns.begin()) ) );
 122     else {
 123         DAI_THROW(BELIEF_NOT_AVAILABLE);
 124         return Factor();
 125     }
 126 }
 127
 128
 129 Real LC::CalcCavityDist (size_t i, const std::string &name, const PropertySet &opts) {
 130     Factor Bi;
 131     Real maxdiff = 0;
 132
 133     if( props.verbose >= 2 )
 134         cerr << "Initing cavity " << var(i) << "(" << delta(i).size() << " vars, " << delta(i).nrStates() << " states)" << endl;
 135
 136     if( props.cavity == Properties::CavityType::UNIFORM )
 137         Bi = Factor(delta(i));
 138     else {
 139         InfAlg *cav = newInfAlg( name, *this, opts );
 140         cav->makeCavity( i );
 141
 142         if( props.cavity == Properties::CavityType::FULL )
 143             Bi = calcMarginal( *cav, cav->fg().delta(i), props.reinit );
 144         else if( props.cavity == Properties::CavityType::PAIR ) {
 145             vector<Factor> pairbeliefs = calcPairBeliefs( *cav, cav->fg().delta(i), props.reinit, false );
 146             for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
 147                 Bi *= pairbeliefs[ij];
 148         } else if( props.cavity == Properties::CavityType::PAIR2 ) {
 149             vector<Factor> pairbeliefs = calcPairBeliefs( *cav, cav->fg().delta(i), props.reinit, true );
 150             for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
 151                 Bi *= pairbeliefs[ij];
 152         }
 153         maxdiff = cav->maxDiff();
 154         delete cav;
 155     }
 156     Bi.normalize();
 157     _cavitydists[i] = Bi;
 158
 159     return maxdiff;
 160 }
 161
 162
 163 Real LC::InitCavityDists( const std::string &name, const PropertySet &opts ) {
 164     double tic = toc();
 165
 166     if( props.verbose >= 1 ) {
 167         cerr << this->name() << "::InitCavityDists:  ";
 168         if( props.cavity == Properties::CavityType::UNIFORM )
 169             cerr << "Using uniform initial cavity distributions" << endl;
 170         else if( props.cavity == Properties::CavityType::FULL )
 171             cerr << "Using full " << name << opts << "...";
 172         else if( props.cavity == Properties::CavityType::PAIR )
 173             cerr << "Using pairwise " << name << opts << "...";
 174         else if( props.cavity == Properties::CavityType::PAIR2 )
 175             cerr << "Using pairwise(new) " << name << opts << "...";
 176     }
 177
 178     Real maxdiff = 0.0;
 179     for( size_t i = 0; i < nrVars(); i++ ) {
 180         Real md = CalcCavityDist(i, name, opts);
 181         if( md > maxdiff )
 182             maxdiff = md;
 183     }
 184
 185     if( props.verbose >= 1 ) {
 186         cerr << this->name() << "::InitCavityDists used " << toc() - tic << " seconds." << endl;
 187     }
 188
 189     return maxdiff;
 190 }
 191
 192
 193 long LC::SetCavityDists( std::vector<Factor> &Q ) {
 194     if( props.verbose >= 1 )
 195         cerr << name() << "::SetCavityDists:  Setting initial cavity distributions" << endl;
 196     if( Q.size() != nrVars() )
 197         return -1;
 198     for( size_t i = 0; i < nrVars(); i++ ) {
 199         if( _cavitydists[i].vars() != Q[i].vars() ) {
 200             return i+1;
 201         } else
 202             _cavitydists[i] = Q[i];
 203     }
 204     return 0;
 205 }
 206
 207
 208 void LC::init() {
 209     for( size_t i = 0; i < nrVars(); ++i )
 210         bforeach( const Neighbor &I, nbV(i) )
 211             if( props.updates == Properties::UpdateType::SEQRND )
 212                 _phis[i][I.iter].randomize();
 213             else
 214                 _phis[i][I.iter].fill(1.0);
 215 }
 216
 217
 218 Factor LC::NewPancake (size_t i, size_t _I, bool & hasNaNs) {
 219     size_t I = nbV(i)[_I];
 220     Factor piet = _pancakes[i];
 221
 222     // recalculate _pancake[i]
 223     VarSet Ivars = factor(I).vars();
 224     Factor A_I;
 225     for( VarSet::const_iterator k = Ivars.begin(); k != Ivars.end(); k++ )
 226         if( var(i) != *k )
 227             A_I *= (_pancakes[findVar(*k)] * factor(I).inverse()).marginal( Ivars / var(i), false );
 228     if( Ivars.size() > 1 )
 229         A_I ^= (1.0 / (Ivars.size() - 1));
 230     Factor A_Ii = (_pancakes[i] * factor(I).inverse() * _phis[i][_I].inverse()).marginal( Ivars / var(i), false );
 231     Factor quot = A_I / A_Ii;
 232     if( props.damping != 0.0 )
 233         quot = (quot^(1.0 - props.damping)) * (_phis[i][_I]^props.damping);
 234
 235     piet *= quot / _phis[i][_I].normalized();
 236     _phis[i][_I] = quot.normalized();
 237
 238     piet.normalize();
 239
 240     if( piet.hasNaNs() ) {
 241         cerr << name() << "::NewPancake(" << i << ", " << _I << "):  has NaNs!" << endl;
 242         hasNaNs = true;
 243     }
 244
 245     return piet;
 246 }
 247
 248
 249 Real LC::run() {
 250     if( props.verbose >= 1 )
 251         cerr << "Starting " << identify() << "...";
 252     if( props.verbose >= 2 )
 253         cerr << endl;
 254
 255     double tic = toc();
 256
 257     Real md = InitCavityDists( props.cavainame, props.cavaiopts );
 258     if( md > _maxdiff )
 259         _maxdiff = md;
 260
 261     for( size_t i = 0; i < nrVars(); i++ ) {
 262         _pancakes[i] = _cavitydists[i];
 263
 264         bforeach( const Neighbor &I, nbV(i) ) {
 265             _pancakes[i] *= factor(I);
 266             if( props.updates == Properties::UpdateType::SEQRND )
 267               _pancakes[i] *= _phis[i][I.iter];
 268         }
 269
 270         _pancakes[i].normalize();
 271
 272         CalcBelief(i);
 273     }
 274
 275     vector<Factor> oldBeliefsV;
 276     for( size_t i = 0; i < nrVars(); i++ )
 277         oldBeliefsV.push_back( beliefV(i) );
 278
 279     bool hasNaNs = false;
 280     for( size_t i=0; i < nrVars(); i++ )
 281         if( _pancakes[i].hasNaNs() ) {
 282             hasNaNs = true;
 283             break;
 284         }
 285     if( hasNaNs ) {
 286         cerr << name() << "::run:  initial _pancakes has NaNs!" << endl;
 287         return 1.0;
 288     }
 289
 290     size_t nredges = nrEdges();
 291     vector<Edge> update_seq;
 292     update_seq.reserve( nredges );
 293     for( size_t i = 0; i < nrVars(); ++i )
 294         bforeach( const Neighbor &I, nbV(i) )
 295             update_seq.push_back( Edge( i, I.iter ) );
 296
 297     // do several passes over the network until maximum number of iterations has
 298     // been reached or until the maximum belief difference is smaller than tolerance
 299     Real maxDiff = INFINITY;
 300     for( _iters = 0; _iters < props.maxiter && maxDiff > props.tol; _iters++ ) {
 301         // Sequential updates
 302         if( props.updates == Properties::UpdateType::SEQRND )
 303             random_shuffle( update_seq.begin(), update_seq.end(), rnd );
 304
 305         for( size_t t=0; t < nredges; t++ ) {
 306             size_t i = update_seq[t].first;
 307             size_t _I = update_seq[t].second;
 308             _pancakes[i] = NewPancake( i, _I, hasNaNs);
 309             if( hasNaNs )
 310                 return 1.0;
 311             CalcBelief( i );
 312         }
 313
 314         // compare new beliefs with old ones
 315         maxDiff = -INFINITY;
 316         for( size_t i = 0; i < nrVars(); i++ ) {
 317             maxDiff = std::max( maxDiff, dist( beliefV(i), oldBeliefsV[i], DISTLINF ) );
 318             oldBeliefsV[i] = beliefV(i);
 319         }
 320
 321         if( props.verbose >= 3 )
 322             cerr << name() << "::run:  maxdiff " << maxDiff << " after " << _iters+1 << " passes" << endl;
 323     }
 324
 325     if( maxDiff > _maxdiff )
 326         _maxdiff = maxDiff;
 327
 328     if( props.verbose >= 1 ) {
 329         if( maxDiff > props.tol ) {
 330             if( props.verbose == 1 )
 331                 cerr << endl;
 332                 cerr << name() << "::run:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << maxDiff << endl;
 333         } else {
 334             if( props.verbose >= 2 )
 335                 cerr << name() << "::run:  ";
 336                 cerr << "converged in " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
 337         }
 338     }
 339
 340     return maxDiff;
 341 }
 342
 343
 344 } // end of namespace dai
 345
 346
 347 #endif