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