src/bp.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 <sstream>
  14 #include <map>
  15 #include <set>
  16 #include <algorithm>
  17 #include <stack>
  18 #include <dai/bp.h>
  19 #include <dai/util.h>
  20 #include <dai/properties.h>
  21
  22
  23 namespace dai {
  24
  25
  26 using namespace std;
  27
  28
  29 const char *BP::Name = "BP";
  30
  31
  32 #define DAI_BP_FAST 1
  33
  34
  35 void BP::setProperties( const PropertySet &opts ) {
  36     DAI_ASSERT( opts.hasKey("tol") );
  37     DAI_ASSERT( opts.hasKey("maxiter") );
  38     DAI_ASSERT( opts.hasKey("logdomain") );
  39     DAI_ASSERT( opts.hasKey("updates") );
  40
  41     props.tol = opts.getStringAs<double>("tol");
  42     props.maxiter = opts.getStringAs<size_t>("maxiter");
  43     props.logdomain = opts.getStringAs<bool>("logdomain");
  44     props.updates = opts.getStringAs<Properties::UpdateType>("updates");
  45
  46     if( opts.hasKey("verbose") )
  47         props.verbose = opts.getStringAs<size_t>("verbose");
  48     else
  49         props.verbose = 0;
  50     if( opts.hasKey("damping") )
  51         props.damping = opts.getStringAs<double>("damping");
  52     else
  53         props.damping = 0.0;
  54     if( opts.hasKey("inference") )
  55         props.inference = opts.getStringAs<Properties::InfType>("inference");
  56     else
  57         props.inference = Properties::InfType::SUMPROD;
  58 }
  59
  60
  61 PropertySet BP::getProperties() const {
  62     PropertySet opts;
  63     opts.Set( "tol", props.tol );
  64     opts.Set( "maxiter", props.maxiter );
  65     opts.Set( "verbose", props.verbose );
  66     opts.Set( "logdomain", props.logdomain );
  67     opts.Set( "updates", props.updates );
  68     opts.Set( "damping", props.damping );
  69     opts.Set( "inference", props.inference );
  70     return opts;
  71 }
  72
  73
  74 string BP::printProperties() const {
  75     stringstream s( stringstream::out );
  76     s << "[";
  77     s << "tol=" << props.tol << ",";
  78     s << "maxiter=" << props.maxiter << ",";
  79     s << "verbose=" << props.verbose << ",";
  80     s << "logdomain=" << props.logdomain << ",";
  81     s << "updates=" << props.updates << ",";
  82     s << "damping=" << props.damping << ",";
  83     s << "inference=" << props.inference << "]";
  84     return s.str();
  85 }
  86
  87
  88 void BP::construct() {
  89     // create edge properties
  90     _edges.clear();
  91     _edges.reserve( nrVars() );
  92     _edge2lut.clear();
  93     if( props.updates == Properties::UpdateType::SEQMAX )
  94         _edge2lut.reserve( nrVars() );
  95     for( size_t i = 0; i < nrVars(); ++i ) {
  96         _edges.push_back( vector<EdgeProp>() );
  97         _edges[i].reserve( nbV(i).size() );
  98         if( props.updates == Properties::UpdateType::SEQMAX ) {
  99             _edge2lut.push_back( vector<LutType::iterator>() );
 100             _edge2lut[i].reserve( nbV(i).size() );
 101         }
 102         foreach( const Neighbor &I, nbV(i) ) {
 103             EdgeProp newEP;
 104             newEP.message = Prob( var(i).states() );
 105             newEP.newMessage = Prob( var(i).states() );
 106
 107             if( DAI_BP_FAST ) {
 108                 newEP.index.reserve( factor(I).states() );
 109                 for( IndexFor k( var(i), factor(I).vars() ); k.valid(); ++k )
 110                     newEP.index.push_back( k );
 111             }
 112
 113             newEP.residual = 0.0;
 114             _edges[i].push_back( newEP );
 115             if( props.updates == Properties::UpdateType::SEQMAX )
 116                 _edge2lut[i].push_back( _lut.insert( make_pair( newEP.residual, make_pair( i, _edges[i].size() - 1 ))) );
 117         }
 118     }
 119 }
 120
 121
 122 void BP::init() {
 123     double c = props.logdomain ? 0.0 : 1.0;
 124     for( size_t i = 0; i < nrVars(); ++i ) {
 125         foreach( const Neighbor &I, nbV(i) ) {
 126             message( i, I.iter ).fill( c );
 127             newMessage( i, I.iter ).fill( c );
 128             if( props.updates == Properties::UpdateType::SEQMAX )
 129                 updateResidual( i, I.iter, 0.0 );
 130         }
 131     }
 132 }
 133
 134
 135 void BP::findMaxResidual( size_t &i, size_t &_I ) {
 136     DAI_ASSERT( !_lut.empty() );
 137     LutType::const_iterator largestEl = _lut.end();
 138     --largestEl;
 139     i  = largestEl->second.first;
 140     _I = largestEl->second.second;
 141 }
 142
 143
 144 void BP::calcNewMessage( size_t i, size_t _I ) {
 145     // calculate updated message I->i
 146     size_t I = nbV(i,_I);
 147
 148     Factor Fprod( factor(I) );
 149     Prob &prod = Fprod.p();
 150     if( props.logdomain )
 151         prod.takeLog();
 152
 153     // Calculate product of incoming messages and factor I
 154     foreach( const Neighbor &j, nbF(I) )
 155         if( j != i ) { // for all j in I \ i
 156             // prod_j will be the product of messages coming into j
 157             Prob prod_j( var(j).states(), props.logdomain ? 0.0 : 1.0 );
 158             foreach( const Neighbor &J, nbV(j) )
 159                 if( J != I ) { // for all J in nb(j) \ I
 160                     if( props.logdomain )
 161                         prod_j += message( j, J.iter );
 162                     else
 163                         prod_j *= message( j, J.iter );
 164                 }
 165
 166             // multiply prod with prod_j
 167             if( !DAI_BP_FAST ) {
 168                 /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
 169                 if( props.logdomain )
 170                     Fprod += Factor( var(j), prod_j );
 171                 else
 172                     Fprod *= Factor( var(j), prod_j );
 173             } else {
 174                 /* OPTIMIZED VERSION */
 175                 size_t _I = j.dual;
 176                 // ind is the precalculated IndexFor(j,I) i.e. to x_I == k corresponds x_j == ind[k]
 177                 const ind_t &ind = index(j, _I);
 178                 for( size_t r = 0; r < prod.size(); ++r )
 179                     if( props.logdomain )
 180                         prod[r] += prod_j[ind[r]];
 181                     else
 182                         prod[r] *= prod_j[ind[r]];
 183             }
 184         }
 185
 186     if( props.logdomain ) {
 187         prod -= prod.max();
 188         prod.takeExp();
 189     }
 190
 191     // Marginalize onto i
 192     Prob marg;
 193     if( !DAI_BP_FAST ) {
 194         /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
 195         if( props.inference == Properties::InfType::SUMPROD )
 196             marg = Fprod.marginal( var(i) ).p();
 197         else
 198             marg = Fprod.maxMarginal( var(i) ).p();
 199     } else {
 200         /* OPTIMIZED VERSION */
 201         marg = Prob( var(i).states(), 0.0 );
 202         // ind is the precalculated IndexFor(i,I) i.e. to x_I == k corresponds x_i == ind[k]
 203         const ind_t ind = index(i,_I);
 204         if( props.inference == Properties::InfType::SUMPROD )
 205             for( size_t r = 0; r < prod.size(); ++r )
 206                 marg[ind[r]] += prod[r];
 207         else
 208             for( size_t r = 0; r < prod.size(); ++r )
 209                 if( prod[r] > marg[ind[r]] )
 210                     marg[ind[r]] = prod[r];
 211         marg.normalize();
 212     }
 213
 214     // Store result
 215     if( props.logdomain )
 216         newMessage(i,_I) = marg.log();
 217     else
 218         newMessage(i,_I) = marg;
 219
 220     // Update the residual if necessary
 221     if( props.updates == Properties::UpdateType::SEQMAX )
 222         updateResidual( i, _I , dist( newMessage( i, _I ), message( i, _I ), Prob::DISTLINF ) );
 223 }
 224
 225
 226 // BP::run does not check for NANs for performance reasons
 227 // Somehow NaNs do not often occur in BP...
 228 double BP::run() {
 229     if( props.verbose >= 1 )
 230         cerr << "Starting " << identify() << "...";
 231     if( props.verbose >= 3)
 232         cerr << endl;
 233
 234     double tic = toc();
 235     Diffs diffs(nrVars(), 1.0);
 236
 237     vector<Edge> update_seq;
 238
 239     vector<Factor> old_beliefs;
 240     old_beliefs.reserve( nrVars() );
 241     for( size_t i = 0; i < nrVars(); ++i )
 242         old_beliefs.push_back( beliefV(i) );
 243
 244     size_t nredges = nrEdges();
 245
 246     if( props.updates == Properties::UpdateType::SEQMAX ) {
 247         // do the first pass
 248         for( size_t i = 0; i < nrVars(); ++i )
 249             foreach( const Neighbor &I, nbV(i) ) {
 250                 calcNewMessage( i, I.iter );
 251             }
 252     } else {
 253         update_seq.reserve( nredges );
 254         /// \todo Investigate whether performance increases by switching the order of following two loops:
 255         for( size_t i = 0; i < nrVars(); ++i )
 256             foreach( const Neighbor &I, nbV(i) )
 257                 update_seq.push_back( Edge( i, I.iter ) );
 258     }
 259
 260     // do several passes over the network until maximum number of iterations has
 261     // been reached or until the maximum belief difference is smaller than tolerance
 262     for( _iters=0; _iters < props.maxiter && diffs.maxDiff() > props.tol; ++_iters ) {
 263         if( props.updates == Properties::UpdateType::SEQMAX ) {
 264             // Residuals-BP by Koller et al.
 265             for( size_t t = 0; t < nredges; ++t ) {
 266                 // update the message with the largest residual
 267                 size_t i, _I;
 268                 findMaxResidual( i, _I );
 269                 updateMessage( i, _I );
 270
 271                 // I->i has been updated, which means that residuals for all
 272                 // J->j with J in nb[i]\I and j in nb[J]\i have to be updated
 273                 foreach( const Neighbor &J, nbV(i) ) {
 274                     if( J.iter != _I ) {
 275                         foreach( const Neighbor &j, nbF(J) ) {
 276                             size_t _J = j.dual;
 277                             if( j != i )
 278                                 calcNewMessage( j, _J );
 279                         }
 280                     }
 281                 }
 282             }
 283         } else if( props.updates == Properties::UpdateType::PARALL ) {
 284             // Parallel updates
 285             for( size_t i = 0; i < nrVars(); ++i )
 286                 foreach( const Neighbor &I, nbV(i) )
 287                     calcNewMessage( i, I.iter );
 288
 289             for( size_t i = 0; i < nrVars(); ++i )
 290                 foreach( const Neighbor &I, nbV(i) )
 291                     updateMessage( i, I.iter );
 292         } else {
 293             // Sequential updates
 294             if( props.updates == Properties::UpdateType::SEQRND )
 295                 random_shuffle( update_seq.begin(), update_seq.end() );
 296
 297             foreach( const Edge &e, update_seq ) {
 298                 calcNewMessage( e.first, e.second );
 299                 updateMessage( e.first, e.second );
 300             }
 301         }
 302
 303         // calculate new beliefs and compare with old ones
 304         for( size_t i = 0; i < nrVars(); ++i ) {
 305             Factor nb( beliefV(i) );
 306             diffs.push( dist( nb, old_beliefs[i], Prob::DISTLINF ) );
 307             old_beliefs[i] = nb;
 308         }
 309
 310         if( props.verbose >= 3 )
 311             cerr << Name << "::run:  maxdiff " << diffs.maxDiff() << " after " << _iters+1 << " passes" << endl;
 312     }
 313
 314     if( diffs.maxDiff() > _maxdiff )
 315         _maxdiff = diffs.maxDiff();
 316
 317     if( props.verbose >= 1 ) {
 318         if( diffs.maxDiff() > props.tol ) {
 319             if( props.verbose == 1 )
 320                 cerr << endl;
 321                 cerr << Name << "::run:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << diffs.maxDiff() << endl;
 322         } else {
 323             if( props.verbose >= 3 )
 324                 cerr << Name << "::run:  ";
 325                 cerr << "converged in " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
 326         }
 327     }
 328
 329     return diffs.maxDiff();
 330 }
 331
 332
 333 void BP::calcBeliefV( size_t i, Prob &p ) const {
 334     p = Prob( var(i).states(), props.logdomain ? 0.0 : 1.0 );
 335     foreach( const Neighbor &I, nbV(i) )
 336         if( props.logdomain )
 337             p += newMessage( i, I.iter );
 338         else
 339             p *= newMessage( i, I.iter );
 340 }
 341
 342
 343 void BP::calcBeliefF( size_t I, Prob &p ) const {
 344     Factor Fprod( factor( I ) );
 345     Prob &prod = Fprod.p();
 346
 347     if( props.logdomain )
 348         prod.takeLog();
 349
 350     foreach( const Neighbor &j, nbF(I) ) {
 351         // prod_j will be the product of messages coming into j
 352         Prob prod_j( var(j).states(), props.logdomain ? 0.0 : 1.0 );
 353         foreach( const Neighbor &J, nbV(j) )
 354             if( J != I ) { // for all J in nb(j) \ I
 355                 if( props.logdomain )
 356                     prod_j += newMessage( j, J.iter );
 357                 else
 358                     prod_j *= newMessage( j, J.iter );
 359             }
 360
 361         // multiply prod with prod_j
 362         if( !DAI_BP_FAST ) {
 363             /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
 364             if( props.logdomain )
 365                 Fprod += Factor( var(j), prod_j );
 366             else
 367                 Fprod *= Factor( var(j), prod_j );
 368         } else {
 369             /* OPTIMIZED VERSION */
 370             size_t _I = j.dual;
 371             // ind is the precalculated IndexFor(j,I) i.e. to x_I == k corresponds x_j == ind[k]
 372             const ind_t & ind = index(j, _I);
 373
 374             for( size_t r = 0; r < prod.size(); ++r ) {
 375                 if( props.logdomain )
 376                     prod[r] += prod_j[ind[r]];
 377                 else
 378                     prod[r] *= prod_j[ind[r]];
 379             }
 380         }
 381     }
 382
 383     p = prod;
 384 }
 385
 386
 387 Factor BP::beliefV( size_t i ) const {
 388     Prob p;
 389     calcBeliefV( i, p );
 390
 391     if( props.logdomain ) {
 392         p -= p.max();
 393         p.takeExp();
 394     }
 395     p.normalize();
 396
 397     return( Factor( var(i), p ) );
 398 }
 399
 400
 401 Factor BP::beliefF( size_t I ) const {
 402     Prob p;
 403     calcBeliefF( I, p );
 404
 405     if( props.logdomain ) {
 406         p -= p.max();
 407         p.takeExp();
 408     }
 409     p.normalize();
 410
 411     return( Factor( factor(I).vars(), p ) );
 412 }
 413
 414
 415 Factor BP::belief( const Var &n ) const {
 416     return( beliefV( findVar( n ) ) );
 417 }
 418
 419
 420 vector<Factor> BP::beliefs() const {
 421     vector<Factor> result;
 422     for( size_t i = 0; i < nrVars(); ++i )
 423         result.push_back( beliefV(i) );
 424     for( size_t I = 0; I < nrFactors(); ++I )
 425         result.push_back( beliefF(I) );
 426     return result;
 427 }
 428
 429
 430 Factor BP::belief( const VarSet &ns ) const {
 431     if( ns.size() == 1 )
 432         return belief( *(ns.begin()) );
 433     else {
 434         size_t I;
 435         for( I = 0; I < nrFactors(); I++ )
 436             if( factor(I).vars() >> ns )
 437                 break;
 438         DAI_ASSERT( I != nrFactors() );
 439         return beliefF(I).marginal(ns);
 440     }
 441 }
 442
 443
 444 Real BP::logZ() const {
 445     Real sum = 0.0;
 446     for(size_t i = 0; i < nrVars(); ++i )
 447         sum += (1.0 - nbV(i).size()) * beliefV(i).entropy();
 448     for( size_t I = 0; I < nrFactors(); ++I )
 449         sum -= dist( beliefF(I), factor(I), Prob::DISTKL );
 450     return sum;
 451 }
 452
 453
 454 string BP::identify() const {
 455     return string(Name) + printProperties();
 456 }
 457
 458
 459 void BP::init( const VarSet &ns ) {
 460     for( VarSet::const_iterator n = ns.begin(); n != ns.end(); ++n ) {
 461         size_t ni = findVar( *n );
 462         foreach( const Neighbor &I, nbV( ni ) ) {
 463             double val = props.logdomain ? 0.0 : 1.0;
 464             message( ni, I.iter ).fill( val );
 465             newMessage( ni, I.iter ).fill( val );
 466             if( props.updates == Properties::UpdateType::SEQMAX )
 467                 updateResidual( ni, I.iter, 0.0 );
 468         }
 469     }
 470 }
 471
 472
 473 void BP::updateMessage( size_t i, size_t _I ) {
 474     if( recordSentMessages )
 475         _sentMessages.push_back(make_pair(i,_I));
 476     if( props.damping == 0.0 ) {
 477         message(i,_I) = newMessage(i,_I);
 478         if( props.updates == Properties::UpdateType::SEQMAX )
 479             updateResidual( i, _I, 0.0 );
 480     } else {
 481         message(i,_I) = (message(i,_I) ^ props.damping) * (newMessage(i,_I) ^ (1.0 - props.damping));
 482         if( props.updates == Properties::UpdateType::SEQMAX )
 483             updateResidual( i, _I, dist( newMessage(i,_I), message(i,_I), Prob::DISTLINF ) );
 484     }
 485 }
 486
 487
 488 void BP::updateResidual( size_t i, size_t _I, double r ) {
 489     EdgeProp* pEdge = &_edges[i][_I];
 490     pEdge->residual = r;
 491
 492     // rearrange look-up table (delete and reinsert new key)
 493     _lut.erase( _edge2lut[i][_I] );
 494     _edge2lut[i][_I] = _lut.insert( make_pair( r, make_pair(i, _I) ) );
 495 }
 496
 497
 498 std::vector<size_t> BP::findMaximum() const {
 499     vector<size_t> maximum( nrVars() );
 500     vector<bool> visitedVars( nrVars(), false );
 501     vector<bool> visitedFactors( nrFactors(), false );
 502     stack<size_t> scheduledFactors;
 503     for( size_t i = 0; i < nrVars(); ++i ) {
 504         if( visitedVars[i] )
 505             continue;
 506         visitedVars[i] = true;
 507
 508         // Maximise with respect to variable i
 509         Prob prod;
 510         calcBeliefV( i, prod );
 511         maximum[i] = max_element( prod.begin(), prod.end() ) - prod.begin();
 512
 513         foreach( const Neighbor &I, nbV(i) )
 514             if( !visitedFactors[I] )
 515                 scheduledFactors.push(I);
 516
 517         while( !scheduledFactors.empty() ){
 518             size_t I = scheduledFactors.top();
 519             scheduledFactors.pop();
 520             if( visitedFactors[I] )
 521                 continue;
 522             visitedFactors[I] = true;
 523
 524             // Evaluate if some neighboring variables still need to be fixed; if not, we're done
 525             bool allDetermined = true;
 526             foreach( const Neighbor &j, nbF(I) )
 527                 if( !visitedVars[j.node] ) {
 528                     allDetermined = false;
 529                     break;
 530                 }
 531             if( allDetermined )
 532                 continue;
 533
 534             // Calculate product of incoming messages on factor I
 535             Prob prod2;
 536             calcBeliefF( I, prod2 );
 537
 538             // The allowed configuration is restrained according to the variables assigned so far:
 539             // pick the argmax amongst the allowed states
 540             Real maxProb = numeric_limits<Real>::min();
 541             State maxState( factor(I).vars() );
 542             for( State s( factor(I).vars() ); s.valid(); ++s ){
 543                 // First, calculate whether this state is consistent with variables that
 544                 // have been assigned already
 545                 bool allowedState = true;
 546                 foreach( const Neighbor &j, nbF(I) )
 547                     if( visitedVars[j.node] && maximum[j.node] != s(var(j.node)) ) {
 548                         allowedState = false;
 549                         break;
 550                     }
 551                 // If it is consistent, check if its probability is larger than what we have seen so far
 552                 if( allowedState && prod2[s] > maxProb ) {
 553                     maxState = s;
 554                     maxProb = prod2[s];
 555                 }
 556             }
 557
 558             // Decode the argmax
 559             foreach( const Neighbor &j, nbF(I) ) {
 560                 if( visitedVars[j.node] ) {
 561                     // We have already visited j earlier - hopefully our state is consistent
 562                     if( maximum[j.node] != maxState(var(j.node)) && props.verbose >= 1 )
 563                         cerr << "BP::findMaximum - warning: maximum not consistent due to loops." << endl;
 564                 } else {
 565                     // We found a consistent state for variable j
 566                     visitedVars[j.node] = true;
 567                     maximum[j.node] = maxState( var(j.node) );
 568                     foreach( const Neighbor &J, nbV(j) )
 569                         if( !visitedFactors[J] )
 570                             scheduledFactors.push(J);
 571                 }
 572             }
 573         }
 574     }
 575     return maximum;
 576 }
 577
 578
 579 } // end of namespace dai