src/hak.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 <map>
  24 #include <dai/hak.h>
  25 #include <dai/util.h>
  26 #include <dai/diffs.h>
  27 #include <dai/exceptions.h>
  28
  29
  30 namespace dai {
  31
  32
  33 using namespace std;
  34
  35
  36 const char *HAK::Name = "HAK";
  37
  38
  39 void HAK::setProperties( const PropertySet &opts ) {
  40     assert( opts.hasKey("tol") );
  41     assert( opts.hasKey("maxiter") );
  42     assert( opts.hasKey("verbose") );
  43     assert( opts.hasKey("doubleloop") );
  44     assert( opts.hasKey("clusters") );
  45
  46     props.tol = opts.getStringAs<double>("tol");
  47     props.maxiter = opts.getStringAs<size_t>("maxiter");
  48     props.verbose = opts.getStringAs<size_t>("verbose");
  49     props.doubleloop = opts.getStringAs<bool>("doubleloop");
  50     props.clusters = opts.getStringAs<Properties::ClustersType>("clusters");
  51
  52     if( opts.hasKey("loopdepth") )
  53         props.loopdepth = opts.getStringAs<size_t>("loopdepth");
  54     else
  55         assert( props.clusters != Properties::ClustersType::LOOP );
  56     if( opts.hasKey("damping") )
  57         props.damping = opts.getStringAs<double>("damping");
  58     else
  59         props.damping = 0.0;
  60 }
  61
  62
  63 PropertySet HAK::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( "doubleloop", props.doubleloop );
  69     opts.Set( "clusters", props.clusters );
  70     opts.Set( "loopdepth", props.loopdepth );
  71     opts.Set( "damping", props.damping );
  72     return opts;
  73 }
  74
  75
  76 string HAK::printProperties() const {
  77     stringstream s( stringstream::out );
  78     s << "[";
  79     s << "tol=" << props.tol << ",";
  80     s << "maxiter=" << props.maxiter << ",";
  81     s << "verbose=" << props.verbose << ",";
  82     s << "doubleloop=" << props.doubleloop << ",";
  83     s << "clusters=" << props.clusters << ",";
  84     s << "loopdepth=" << props.loopdepth << ",";
  85     s << "damping=" << props.damping << "]";
  86     return s.str();
  87 }
  88
  89
  90 void HAK::constructMessages() {
  91     // Create outer beliefs
  92     _Qa.clear();
  93     _Qa.reserve(nrORs());
  94     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
  95         _Qa.push_back( Factor( OR(alpha).vars() ) );
  96
  97     // Create inner beliefs
  98     _Qb.clear();
  99     _Qb.reserve(nrIRs());
 100     for( size_t beta = 0; beta < nrIRs(); beta++ )
 101         _Qb.push_back( Factor( IR(beta) ) );
 102
 103     // Create messages
 104     _muab.clear();
 105     _muab.reserve( nrORs() );
 106     _muba.clear();
 107     _muba.reserve( nrORs() );
 108     for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
 109         _muab.push_back( vector<Factor>() );
 110         _muba.push_back( vector<Factor>() );
 111         _muab[alpha].reserve( nbOR(alpha).size() );
 112         _muba[alpha].reserve( nbOR(alpha).size() );
 113         foreach( const Neighbor &beta, nbOR(alpha) ) {
 114             _muab[alpha].push_back( Factor( IR(beta) ) );
 115             _muba[alpha].push_back( Factor( IR(beta) ) );
 116         }
 117     }
 118 }
 119
 120
 121 HAK::HAK( const RegionGraph &rg, const PropertySet &opts ) : DAIAlgRG(rg), _Qa(), _Qb(), _muab(), _muba(), _maxdiff(0.0), _iters(0U), props() {
 122     setProperties( opts );
 123
 124     constructMessages();
 125 }
 126
 127
 128 void HAK::findLoopClusters( const FactorGraph & fg, std::set<VarSet> &allcl, VarSet newcl, const Var & root, size_t length, VarSet vars ) {
 129     for( VarSet::const_iterator in = vars.begin(); in != vars.end(); in++ ) {
 130         VarSet ind = fg.delta( fg.findVar( *in ) );
 131         if( (newcl.size()) >= 2 && ind.contains( root ) ) {
 132             allcl.insert( newcl | *in );
 133         }
 134         else if( length > 1 )
 135             findLoopClusters( fg, allcl, newcl | *in, root, length - 1, ind / newcl );
 136     }
 137 }
 138
 139
 140 HAK::HAK(const FactorGraph & fg, const PropertySet &opts) : DAIAlgRG(), _Qa(), _Qb(), _muab(), _muba(), _maxdiff(0.0), _iters(0U), props() {
 141     setProperties( opts );
 142
 143     vector<VarSet> cl;
 144     if( props.clusters == Properties::ClustersType::MIN ) {
 145         cl = fg.Cliques();
 146     } else if( props.clusters == Properties::ClustersType::DELTA ) {
 147         for( size_t i = 0; i < fg.nrVars(); i++ )
 148             cl.push_back(fg.Delta(i));
 149     } else if( props.clusters == Properties::ClustersType::LOOP ) {
 150         cl = fg.Cliques();
 151         set<VarSet> scl;
 152         for( size_t i0 = 0; i0 < fg.nrVars(); i0++ ) {
 153             VarSet i0d = fg.delta(i0);
 154             if( props.loopdepth > 1 )
 155                 findLoopClusters( fg, scl, fg.var(i0), fg.var(i0), props.loopdepth - 1, fg.delta(i0) );
 156         }
 157         for( set<VarSet>::const_iterator c = scl.begin(); c != scl.end(); c++ )
 158             cl.push_back(*c);
 159         if( props.verbose >= 3 ) {
 160             cout << Name << " uses the following clusters: " << endl;
 161             for( vector<VarSet>::const_iterator cli = cl.begin(); cli != cl.end(); cli++ )
 162                 cout << *cli << endl;
 163         }
 164     } else
 165         DAI_THROW(INTERNAL_ERROR);
 166
 167     RegionGraph rg(fg,cl);
 168     RegionGraph::operator=(rg);
 169     constructMessages();
 170
 171     if( props.verbose >= 3 )
 172         cout << Name << " regiongraph: " << *this << endl;
 173 }
 174
 175
 176 string HAK::identify() const {
 177     return string(Name) + printProperties();
 178 }
 179
 180
 181 void HAK::init( const VarSet &ns ) {
 182     for( vector<Factor>::iterator alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
 183         if( alpha->vars().intersects( ns ) )
 184             alpha->fill( 1.0 / alpha->states() );
 185
 186     for( size_t beta = 0; beta < nrIRs(); beta++ )
 187         if( IR(beta).intersects( ns ) ) {
 188             _Qb[beta].fill( 1.0 );
 189             foreach( const Neighbor &alpha, nbIR(beta) ) {
 190                 size_t _beta = alpha.dual;
 191                 muab( alpha, _beta ).fill( 1.0 );
 192                 muba( alpha, _beta ).fill( 1.0 );
 193             }
 194         }
 195 }
 196
 197
 198 void HAK::init() {
 199     for( vector<Factor>::iterator alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
 200         alpha->fill( 1.0 / alpha->states() );
 201
 202     for( vector<Factor>::iterator beta = _Qb.begin(); beta != _Qb.end(); beta++ )
 203         beta->fill( 1.0 / beta->states() );
 204
 205     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
 206         foreach( const Neighbor &beta, nbOR(alpha) ) {
 207             size_t _beta = beta.iter;
 208             muab( alpha, _beta ).fill( 1.0 / muab( alpha, _beta ).states() );
 209             muba( alpha, _beta ).fill( 1.0 / muab( alpha, _beta ).states() );
 210         }
 211 }
 212
 213
 214 double HAK::doGBP() {
 215     if( props.verbose >= 1 )
 216         cout << "Starting " << identify() << "...";
 217     if( props.verbose >= 3)
 218         cout << endl;
 219
 220     double tic = toc();
 221
 222     // Check whether counting numbers won't lead to problems
 223     for( size_t beta = 0; beta < nrIRs(); beta++ )
 224         assert( nbIR(beta).size() + IR(beta).c() != 0.0 );
 225
 226     // Keep old beliefs to check convergence
 227     vector<Factor> old_beliefs;
 228     old_beliefs.reserve( nrVars() );
 229     for( size_t i = 0; i < nrVars(); i++ )
 230         old_beliefs.push_back( belief( var(i) ) );
 231
 232     // Differences in single node beliefs
 233     Diffs diffs(nrVars(), 1.0);
 234
 235     // do several passes over the network until maximum number of iterations has
 236     // been reached or until the maximum belief difference is smaller than tolerance
 237     for( _iters = 0; _iters < props.maxiter && diffs.maxDiff() > props.tol; _iters++ ) {
 238         for( size_t beta = 0; beta < nrIRs(); beta++ ) {
 239             foreach( const Neighbor &alpha, nbIR(beta) ) {
 240                 size_t _beta = alpha.dual;
 241                 muab( alpha, _beta ) = _Qa[alpha].marginal(IR(beta)).divided_by( muba(alpha,_beta) );
 242                 /* TODO: INVESTIGATE THIS PROBLEM
 243                  *
 244                  * In some cases, the muab's can have very large entries because the muba's have very
 245                  * small entries. This may cause NANs later on (e.g., multiplying large quantities may
 246                  * result in +inf; normalization then tries to calculate inf / inf which is NAN).
 247                  * A fix of this problem would consist in normalizing the messages muab.
 248                  * However, it is not obvious whether this is a real solution, because it has a
 249                  * negative performance impact and the NAN's seem to be a symptom of a fundamental
 250                  * numerical unstability.
 251                  */
 252                  muab(alpha,_beta).normalize();
 253             }
 254
 255             Factor Qb_new;
 256             foreach( const Neighbor &alpha, nbIR(beta) ) {
 257                 size_t _beta = alpha.dual;
 258                 Qb_new *= muab(alpha,_beta) ^ (1 / (nbIR(beta).size() + IR(beta).c()));
 259             }
 260
 261             Qb_new.normalize();
 262             if( Qb_new.hasNaNs() ) {
 263                 // TODO: WHAT TO DO IN THIS CASE?
 264                 cout << Name << "::doGBP:  Qb_new has NaNs!" << endl;
 265                 return 1.0;
 266             }
 267             /* TODO: WHAT IS THE PURPOSE OF THE FOLLOWING CODE?
 268              *
 269              *   _Qb[beta] = Qb_new.makeZero(1e-100);
 270              */
 271
 272             if( props.doubleloop || props.damping == 0.0 )
 273                 _Qb[beta] = Qb_new; // no damping for double loop
 274             else
 275                 _Qb[beta] = (Qb_new^(1.0 - props.damping)) * (_Qb[beta]^props.damping);
 276
 277             foreach( const Neighbor &alpha, nbIR(beta) ) {
 278                 size_t _beta = alpha.dual;
 279                 muba(alpha,_beta) = _Qb[beta].divided_by( muab(alpha,_beta) );
 280
 281                 /* TODO: INVESTIGATE WHETHER THIS HACK (INVENTED BY KEES) TO PREVENT NANS MAKES SENSE
 282                  *
 283                  *   muba(beta,*alpha).makePositive(1e-100);
 284                  *
 285                  */
 286
 287                 Factor Qa_new = OR(alpha);
 288                 foreach( const Neighbor &gamma, nbOR(alpha) )
 289                     Qa_new *= muba(alpha,gamma.iter);
 290                 Qa_new ^= (1.0 / OR(alpha).c());
 291                 Qa_new.normalize();
 292                 if( Qa_new.hasNaNs() ) {
 293                     cout << Name << "::doGBP:  Qa_new has NaNs!" << endl;
 294                     return 1.0;
 295                 }
 296                 /* TODO: WHAT IS THE PURPOSE OF THE FOLLOWING CODE?
 297                  *
 298                  *   _Qb[beta] = Qb_new.makeZero(1e-100);
 299                  */
 300
 301                 if( props.doubleloop || props.damping == 0.0 )
 302                     _Qa[alpha] = Qa_new; // no damping for double loop
 303                 else
 304                     // FIXME: GEOMETRIC DAMPING IS SLOW!
 305                 _Qa[alpha] = (Qa_new^(1.0 - props.damping)) * (_Qa[alpha]^props.damping);
 306             }
 307         }
 308
 309         // Calculate new single variable beliefs and compare with old ones
 310         for( size_t i = 0; i < nrVars(); i++ ) {
 311             Factor new_belief = belief( var( i ) );
 312             diffs.push( dist( new_belief, old_beliefs[i], Prob::DISTLINF ) );
 313             old_beliefs[i] = new_belief;
 314         }
 315
 316         if( props.verbose >= 3 )
 317             cout << Name << "::doGBP:  maxdiff " << diffs.maxDiff() << " after " << _iters+1 << " passes" << endl;
 318     }
 319
 320     if( diffs.maxDiff() > _maxdiff )
 321         _maxdiff = diffs.maxDiff();
 322
 323     if( props.verbose >= 1 ) {
 324         if( diffs.maxDiff() > props.tol ) {
 325             if( props.verbose == 1 )
 326                 cout << endl;
 327             cout << Name << "::doGBP:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << diffs.maxDiff() << endl;
 328         } else {
 329             if( props.verbose >= 2 )
 330                 cout << Name << "::doGBP:  ";
 331             cout << "converged in " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
 332         }
 333     }
 334
 335     return diffs.maxDiff();
 336 }
 337
 338
 339 double HAK::doDoubleLoop() {
 340     if( props.verbose >= 1 )
 341         cout << "Starting " << identify() << "...";
 342     if( props.verbose >= 3)
 343         cout << endl;
 344
 345     double tic = toc();
 346
 347     // Save original outer regions
 348     vector<FRegion> org_ORs = ORs;
 349
 350     // Save original inner counting numbers and set negative counting numbers to zero
 351     vector<double> org_IR_cs( nrIRs(), 0.0 );
 352     for( size_t beta = 0; beta < nrIRs(); beta++ ) {
 353         org_IR_cs[beta] = IR(beta).c();
 354         if( IR(beta).c() < 0.0 )
 355             IR(beta).c() = 0.0;
 356     }
 357
 358     // Keep old beliefs to check convergence
 359     vector<Factor> old_beliefs;
 360     old_beliefs.reserve( nrVars() );
 361     for( size_t i = 0; i < nrVars(); i++ )
 362         old_beliefs.push_back( belief( var(i) ) );
 363
 364     // Differences in single node beliefs
 365     Diffs diffs(nrVars(), 1.0);
 366
 367     size_t outer_maxiter   = props.maxiter;
 368     double outer_tol       = props.tol;
 369     size_t outer_verbose   = props.verbose;
 370     double org_maxdiff     = _maxdiff;
 371
 372     // Set parameters for inner loop
 373     props.maxiter = 5;
 374     props.verbose = outer_verbose ? outer_verbose - 1 : 0;
 375
 376     size_t outer_iter = 0;
 377     size_t total_iter = 0;
 378     for( outer_iter = 0; outer_iter < outer_maxiter && diffs.maxDiff() > outer_tol; outer_iter++ ) {
 379         // Calculate new outer regions
 380         for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
 381             OR(alpha) = org_ORs[alpha];
 382             foreach( const Neighbor &beta, nbOR(alpha) )
 383                 OR(alpha) *= _Qb[beta] ^ ((IR(beta).c() - org_IR_cs[beta]) / nbIR(beta).size());
 384         }
 385
 386         // Inner loop
 387         if( isnan( doGBP() ) )
 388             return 1.0;
 389
 390         // Calculate new single variable beliefs and compare with old ones
 391         for( size_t i = 0; i < nrVars(); ++i ) {
 392             Factor new_belief = belief( var( i ) );
 393             diffs.push( dist( new_belief, old_beliefs[i], Prob::DISTLINF ) );
 394             old_beliefs[i] = new_belief;
 395         }
 396
 397         total_iter += Iterations();
 398
 399         if( props.verbose >= 3 )
 400             cout << Name << "::doDoubleLoop:  maxdiff " << diffs.maxDiff() << " after " << total_iter << " passes" << endl;
 401     }
 402
 403     // restore _maxiter, _verbose and _maxdiff
 404     props.maxiter = outer_maxiter;
 405     props.verbose = outer_verbose;
 406     _maxdiff = org_maxdiff;
 407
 408     _iters = total_iter;
 409     if( diffs.maxDiff() > _maxdiff )
 410         _maxdiff = diffs.maxDiff();
 411
 412     // Restore original outer regions
 413     ORs = org_ORs;
 414
 415     // Restore original inner counting numbers
 416     for( size_t beta = 0; beta < nrIRs(); ++beta )
 417         IR(beta).c() = org_IR_cs[beta];
 418
 419     if( props.verbose >= 1 ) {
 420         if( diffs.maxDiff() > props.tol ) {
 421             if( props.verbose == 1 )
 422                 cout << endl;
 423                 cout << Name << "::doDoubleLoop:  WARNING: not converged within " << outer_maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << diffs.maxDiff() << endl;
 424             } else {
 425                 if( props.verbose >= 3 )
 426                     cout << Name << "::doDoubleLoop:  ";
 427                 cout << "converged in " << total_iter << " passes (" << toc() - tic << " seconds)." << endl;
 428             }
 429         }
 430
 431     return diffs.maxDiff();
 432 }
 433
 434
 435 double HAK::run() {
 436     if( props.doubleloop )
 437         return doDoubleLoop();
 438     else
 439         return doGBP();
 440 }
 441
 442
 443 Factor HAK::belief( const VarSet &ns ) const {
 444     vector<Factor>::const_iterator beta;
 445     for( beta = _Qb.begin(); beta != _Qb.end(); beta++ )
 446         if( beta->vars() >> ns )
 447             break;
 448     if( beta != _Qb.end() )
 449         return( beta->marginal(ns) );
 450     else {
 451         vector<Factor>::const_iterator alpha;
 452         for( alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
 453             if( alpha->vars() >> ns )
 454                 break;
 455         assert( alpha != _Qa.end() );
 456         return( alpha->marginal(ns) );
 457     }
 458 }
 459
 460
 461 Factor HAK::belief( const Var &n ) const {
 462     return belief( (VarSet)n );
 463 }
 464
 465
 466 vector<Factor> HAK::beliefs() const {
 467     vector<Factor> result;
 468     for( size_t beta = 0; beta < nrIRs(); beta++ )
 469         result.push_back( Qb(beta) );
 470     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
 471         result.push_back( Qa(alpha) );
 472     return result;
 473 }
 474
 475
 476 Real HAK::logZ() const {
 477     Real sum = 0.0;
 478     for( size_t beta = 0; beta < nrIRs(); beta++ )
 479         sum += IR(beta).c() * Qb(beta).entropy();
 480     for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
 481         sum += OR(alpha).c() * Qa(alpha).entropy();
 482         sum += (OR(alpha).log0() * Qa(alpha)).totalSum();
 483     }
 484     return sum;
 485 }
 486
 487
 488 } // end of namespace dai