src/lc.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 <algorithm>
  24 #include <map>
  25 #include <set>
  26 #include <dai/lc.h>
  27 #include <dai/diffs.h>
  28 #include <dai/util.h>
  29 #include <dai/alldai.h>
  30 #include <dai/x2x.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 }
  61
  62
  63 PropertySet LC::getProperties() const {
  64     PropertySet opts;
  65     opts.Set( "tol", props.tol );
  66     opts.Set( "maxiter", props.maxiter );
  67     opts.Set( "verbose", props.verbose );
  68     opts.Set( "cavity", props.cavity );
  69     opts.Set( "updates", props.updates );
  70     opts.Set( "cavainame", props.cavainame );
  71     opts.Set( "cavaiopts", props.cavaiopts );
  72     opts.Set( "reinit", props.reinit );
  73     return opts;
  74 }
  75
  76
  77 string LC::printProperties() const {
  78     stringstream s( stringstream::out );
  79     s << "[";
  80     s << "tol=" << props.tol << ",";
  81     s << "maxiter=" << props.maxiter << ",";
  82     s << "verbose=" << props.verbose << ",";
  83     s << "cavity=" << props.cavity << ",";
  84     s << "updates=" << props.updates << ",";
  85     s << "cavainame=" << props.cavainame << ",";
  86     s << "cavaiopts=" << props.cavaiopts << ",";
  87     s << "reinit=" << props.reinit << "]";
  88     return s.str();
  89 }
  90
  91
  92 LC::LC( const FactorGraph & fg, const PropertySet &opts ) : DAIAlgFG(fg), _pancakes(), _cavitydists(), _phis(), _beliefs(), props(), maxdiff(0.0) {
  93     setProperties( opts );
  94
  95     // create pancakes
  96     _pancakes.resize(nrVars());
  97
  98     // create cavitydists
  99     for( size_t i=0; i < nrVars(); i++ )
 100         _cavitydists.push_back(Factor(delta(i)));
 101
 102     // create phis
 103     _phis.reserve( nrVars() );
 104     for( size_t i = 0; i < nrVars(); i++ ) {
 105         _phis.push_back( vector<Factor>() );
 106         _phis[i].reserve( nbV(i).size() );
 107         foreach( const Neighbor &I, nbV(i) )
 108             _phis[i].push_back( Factor( factor(I).vars() / var(i) ) );
 109     }
 110
 111     // create beliefs
 112     for( size_t i=0; i < nrVars(); i++ )
 113         _beliefs.push_back(Factor(var(i)));
 114 }
 115
 116
 117 string LC::identify() const {
 118     return string(Name) + printProperties();
 119 }
 120
 121
 122 void LC::CalcBelief (size_t i) {
 123     _beliefs[i] = _pancakes[i].marginal(var(i));
 124 }
 125
 126
 127 double LC::CalcCavityDist (size_t i, const std::string &name, const PropertySet &opts) {
 128     Factor Bi;
 129     double maxdiff = 0;
 130
 131     if( props.verbose >= 2 )
 132         cout << "Initing cavity " << var(i) << "(" << delta(i).size() << " vars, " << delta(i).states() << " states)" << endl;
 133
 134     if( props.cavity == Properties::CavityType::UNIFORM )
 135         Bi = Factor(delta(i));
 136     else {
 137         InfAlg *cav = newInfAlg( name, *this, opts );
 138         cav->makeCavity( i );
 139
 140         if( props.cavity == Properties::CavityType::FULL )
 141             Bi = calcMarginal( *cav, cav->fg().delta(i), props.reinit );
 142         else if( props.cavity == Properties::CavityType::PAIR )
 143             Bi = calcMarginal2ndO( *cav, cav->fg().delta(i), props.reinit );
 144         else if( props.cavity == Properties::CavityType::PAIR2 ) {
 145             vector<Factor> pairbeliefs = calcPairBeliefsNew( *cav, cav->fg().delta(i), props.reinit );
 146             for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
 147                 Bi *= pairbeliefs[ij];
 148         } else if( props.cavity == Properties::CavityType::PAIRINT ) {
 149             Bi = calcMarginal( *cav, cav->fg().delta(i), props.reinit );
 150
 151             // Set interactions of order > 2 to zero
 152             size_t N = delta(i).size();
 153             Real *p = &(*Bi.p().p().begin());
 154             x2x::p2logp (N, p);
 155             x2x::logp2w (N, p);
 156             x2x::fill (N, p, 2, 0.0);
 157             x2x::w2logp (N, p);
 158 //            x2x::logpnorm (N, p);
 159             x2x::logp2p (N, p);
 160         } else if( props.cavity == Properties::CavityType::PAIRCUM ) {
 161             Bi = calcMarginal( *cav, cav->fg().delta(i), props.reinit );
 162
 163             // Set cumulants of order > 2 to zero
 164             size_t N = delta(i).size();
 165             Real *p = &(*Bi.p().p().begin());
 166             x2x::p2m (N, p);
 167             x2x::m2c (N, p, N);
 168             x2x::fill (N, p, 2, 0.0);
 169             x2x::c2m (N, p, N);
 170             x2x::m2p (N, p);
 171         }
 172         maxdiff = cav->maxDiff();
 173         delete cav;
 174     }
 175     Bi.normalize( Prob::NORMPROB );
 176     _cavitydists[i] = Bi;
 177
 178     return maxdiff;
 179 }
 180
 181
 182 double LC::InitCavityDists( const std::string &name, const PropertySet &opts ) {
 183     double tic = toc();
 184
 185     if( props.verbose >= 1 ) {
 186         cout << "LC::InitCavityDists:  ";
 187         if( props.cavity == Properties::CavityType::UNIFORM )
 188             cout << "Using uniform initial cavity distributions" << endl;
 189         else if( props.cavity == Properties::CavityType::FULL )
 190             cout << "Using full " << name << opts << "...";
 191         else if( props.cavity == Properties::CavityType::PAIR )
 192             cout << "Using pairwise " << name << opts << "...";
 193         else if( props.cavity == Properties::CavityType::PAIR2 )
 194             cout << "Using pairwise(new) " << name << opts << "...";
 195     }
 196
 197     double maxdiff = 0.0;
 198     for( size_t i = 0; i < nrVars(); i++ ) {
 199         double md = CalcCavityDist(i, name, opts);
 200         if( md > maxdiff )
 201             maxdiff = md;
 202     }
 203     init();
 204
 205     if( props.verbose >= 1 ) {
 206         cout << "used " << toc() - tic << " clocks." << endl;
 207     }
 208
 209     return maxdiff;
 210 }
 211
 212
 213 long LC::SetCavityDists( std::vector<Factor> &Q ) {
 214     if( props.verbose >= 1 )
 215         cout << "LC::SetCavityDists:  Setting initial cavity distributions" << endl;
 216     if( Q.size() != nrVars() )
 217         return -1;
 218     for( size_t i = 0; i < nrVars(); i++ ) {
 219         if( _cavitydists[i].vars() != Q[i].vars() ) {
 220             return i+1;
 221         } else
 222             _cavitydists[i] = Q[i];
 223     }
 224     init();
 225     return 0;
 226 }
 227
 228
 229 void LC::init() {
 230     for( size_t i = 0; i < nrVars(); ++i )
 231         foreach( const Neighbor &I, nbV(i) )
 232             if( props.updates == Properties::UpdateType::SEQRND )
 233                 _phis[i][I.iter].randomize();
 234             else
 235                 _phis[i][I.iter].fill(1.0);
 236     for( size_t i = 0; i < nrVars(); i++ ) {
 237         _pancakes[i] = _cavitydists[i];
 238
 239         foreach( const Neighbor &I, nbV(i) ) {
 240             _pancakes[i] *= factor(I);
 241             if( props.updates == Properties::UpdateType::SEQRND )
 242               _pancakes[i] *= _phis[i][I.iter];
 243         }
 244
 245         _pancakes[i].normalize( Prob::NORMPROB );
 246
 247         CalcBelief(i);
 248     }
 249 }
 250
 251
 252 Factor LC::NewPancake (size_t i, size_t _I, bool & hasNaNs) {
 253     size_t I = nbV(i)[_I];
 254     Factor piet = _pancakes[i];
 255
 256     // recalculate _pancake[i]
 257     VarSet Ivars = factor(I).vars();
 258     Factor A_I;
 259     for( VarSet::const_iterator k = Ivars.begin(); k != Ivars.end(); k++ )
 260         if( var(i) != *k )
 261             A_I *= (_pancakes[findVar(*k)] * factor(I).inverse()).partSum( Ivars / var(i) );
 262     if( Ivars.size() > 1 )
 263         A_I ^= (1.0 / (Ivars.size() - 1));
 264     Factor A_Ii = (_pancakes[i] * factor(I).inverse() * _phis[i][_I].inverse()).partSum( Ivars / var(i) );
 265     Factor quot = A_I.divided_by(A_Ii);
 266
 267     piet *= quot.divided_by( _phis[i][_I] ).normalized( Prob::NORMPROB );
 268     _phis[i][_I] = quot.normalized( Prob::NORMPROB );
 269
 270     piet.normalize( Prob::NORMPROB );
 271
 272     if( piet.hasNaNs() ) {
 273         cout << "LC::NewPancake(" << i << ", " << _I << "):  has NaNs!" << endl;
 274         hasNaNs = true;
 275     }
 276
 277     return piet;
 278 }
 279
 280
 281 double LC::run() {
 282     if( props.verbose >= 1 )
 283         cout << "Starting " << identify() << "...";
 284     if( props.verbose >= 2 )
 285         cout << endl;
 286
 287     double tic = toc();
 288     Diffs diffs(nrVars(), 1.0);
 289
 290     double md = InitCavityDists( props.cavainame, props.cavaiopts );
 291     if( md > maxdiff )
 292         maxdiff = md;
 293
 294     vector<Factor> old_beliefs;
 295     for(size_t i=0; i < nrVars(); i++ )
 296         old_beliefs.push_back(belief(i));
 297
 298     bool hasNaNs = false;
 299     for( size_t i=0; i < nrVars(); i++ )
 300         if( _pancakes[i].hasNaNs() ) {
 301             hasNaNs = true;
 302             break;
 303         }
 304     if( hasNaNs ) {
 305         cout << "LC::run:  initial _pancakes has NaNs!" << endl;
 306         return -1.0;
 307     }
 308
 309     size_t nredges = nrEdges();
 310     vector<Edge> update_seq;
 311     update_seq.reserve( nredges );
 312     for( size_t i = 0; i < nrVars(); ++i )
 313         foreach( const Neighbor &I, nbV(i) )
 314             update_seq.push_back( Edge( i, I.iter ) );
 315
 316     size_t iter = 0;
 317
 318     // do several passes over the network until maximum number of iterations has
 319     // been reached or until the maximum belief difference is smaller than tolerance
 320     for( iter=0; iter < props.maxiter && diffs.maxDiff() > props.tol; iter++ ) {
 321         // Sequential updates
 322         if( props.updates == Properties::UpdateType::SEQRND )
 323             random_shuffle( update_seq.begin(), update_seq.end() );
 324
 325         for( size_t t=0; t < nredges; t++ ) {
 326             size_t i = update_seq[t].first;
 327             size_t _I = update_seq[t].second;
 328             _pancakes[i] = NewPancake( i, _I, hasNaNs);
 329             if( hasNaNs )
 330                 return -1.0;
 331             CalcBelief( i );
 332         }
 333
 334         // compare new beliefs with old ones
 335         for(size_t i=0; i < nrVars(); i++ ) {
 336             diffs.push( dist( belief(i), old_beliefs[i], Prob::DISTLINF ) );
 337             old_beliefs[i] = belief(i);
 338         }
 339
 340         if( props.verbose >= 3 )
 341             cout << "LC::run:  maxdiff " << diffs.maxDiff() << " after " << iter+1 << " passes" << endl;
 342     }
 343
 344     if( diffs.maxDiff() > maxdiff )
 345         maxdiff = diffs.maxDiff();
 346
 347     if( props.verbose >= 1 ) {
 348         if( diffs.maxDiff() > props.tol ) {
 349             if( props.verbose == 1 )
 350                 cout << endl;
 351                 cout << "LC::run:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.maxDiff() << endl;
 352         } else {
 353             if( props.verbose >= 2 )
 354                 cout << "LC::run:  ";
 355                 cout << "converged in " << iter << " passes (" << toc() - tic << " clocks)." << endl;
 356         }
 357     }
 358
 359     return diffs.maxDiff();
 360 }
 361
 362
 363 } // end of namespace dai