src/lc.cpp

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