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-2010  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         constructCVM( fg, cl );
 159     } else if( props.clusters == Properties::ClustersType::DELTA ) {
 160         cl.reserve( fg.nrVars() );
 161         for( size_t i = 0; i < fg.nrVars(); i++ )
 162             cl.push_back( fg.Delta(i) );
 163         constructCVM( fg, cl );
 164     } else if( props.clusters == Properties::ClustersType::LOOP ) {
 165         cl = fg.Cliques();
 166         set<VarSet> scl;
 167         for( size_t i0 = 0; i0 < fg.nrVars(); i0++ ) {
 168             VarSet i0d = fg.delta(i0);
 169             if( props.loopdepth > 1 )
 170                 findLoopClusters( fg, scl, fg.var(i0), fg.var(i0), props.loopdepth - 1, fg.delta(i0) );
 171         }
 172         for( set<VarSet>::const_iterator c = scl.begin(); c != scl.end(); c++ )
 173             cl.push_back(*c);
 174         if( props.verbose >= 3 ) {
 175             cerr << Name << " uses the following clusters: " << endl;
 176             for( vector<VarSet>::const_iterator cli = cl.begin(); cli != cl.end(); cli++ )
 177                 cerr << *cli << endl;
 178         }
 179         constructCVM( fg, cl );
 180     } else if( props.clusters == Properties::ClustersType::BETHE ) {
 181         // build outer regions (the cliques)
 182         cl = fg.Cliques();
 183         size_t nrEdges = 0;
 184         for( size_t c = 0; c < cl.size(); c++ )
 185             nrEdges += cl[c].size();
 186
 187         // build inner regions (single variables)
 188         vector<Region> irs;
 189         irs.reserve( fg.nrVars() );
 190         for( size_t i = 0; i < fg.nrVars(); i++ )
 191             irs.push_back( Region( fg.var(i), 1.0 ) );
 192
 193         // build edges (an outer and inner region are connected if the outer region contains the inner one)
 194         // and calculate counting number for inner regions
 195         vector<std::pair<size_t, size_t> > edges;
 196         edges.reserve( nrEdges );
 197         for( size_t c = 0; c < cl.size(); c++ )
 198             for( size_t i = 0; i < irs.size(); i++ )
 199                 if( cl[c] >> irs[i] ) {
 200                     edges.push_back( make_pair( c, i ) );
 201                     irs[i].c() -= 1.0;
 202                 }
 203
 204         // build region graph
 205         RegionGraph::construct( fg, cl, irs, edges );
 206     } else
 207         DAI_THROW(UNKNOWN_ENUM_VALUE);
 208
 209     construct();
 210
 211     if( props.verbose >= 3 )
 212         cerr << Name << " regiongraph: " << *this << endl;
 213 }
 214
 215
 216 string HAK::identify() const {
 217     return string(Name) + printProperties();
 218 }
 219
 220
 221 void HAK::init( const VarSet &ns ) {
 222     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
 223         if( _Qa[alpha].vars().intersects( ns ) ) {
 224             if( props.init == Properties::InitType::UNIFORM )
 225                 _Qa[alpha].setUniform();
 226             else
 227                 _Qa[alpha].randomize();
 228             _Qa[alpha] *= OR(alpha);
 229             _Qa[alpha].normalize();
 230         }
 231
 232     for( size_t beta = 0; beta < nrIRs(); beta++ )
 233         if( IR(beta).intersects( ns ) ) {
 234             if( props.init == Properties::InitType::UNIFORM )
 235                 _Qb[beta].fill( 1.0 );
 236             else
 237                 _Qb[beta].randomize();
 238             foreach( const Neighbor &alpha, nbIR(beta) ) {
 239                 size_t _beta = alpha.dual;
 240                 if( props.init == Properties::InitType::UNIFORM ) {
 241                     muab( alpha, _beta ).fill( 1.0 );
 242                     muba( alpha, _beta ).fill( 1.0 );
 243                 } else {
 244                     muab( alpha, _beta ).randomize();
 245                     muba( alpha, _beta ).randomize();
 246                 }
 247             }
 248         }
 249 }
 250
 251
 252 void HAK::init() {
 253     for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
 254         if( props.init == Properties::InitType::UNIFORM )
 255             _Qa[alpha].setUniform();
 256         else
 257             _Qa[alpha].randomize();
 258         _Qa[alpha] *= OR(alpha);
 259         _Qa[alpha].normalize();
 260     }
 261
 262     for( size_t beta = 0; beta < nrIRs(); beta++ )
 263         if( props.init == Properties::InitType::UNIFORM )
 264             _Qb[beta].setUniform();
 265         else
 266             _Qb[beta].randomize();
 267
 268     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
 269         foreach( const Neighbor &beta, nbOR(alpha) ) {
 270             size_t _beta = beta.iter;
 271             if( props.init == Properties::InitType::UNIFORM ) {
 272                 muab( alpha, _beta ).setUniform();
 273                 muba( alpha, _beta ).setUniform();
 274             } else {
 275                 muab( alpha, _beta ).randomize();
 276                 muba( alpha, _beta ).randomize();
 277             }
 278         }
 279 }
 280
 281
 282 Real HAK::doGBP() {
 283     if( props.verbose >= 1 )
 284         cerr << "Starting " << identify() << "...";
 285     if( props.verbose >= 3)
 286         cerr << endl;
 287
 288     double tic = toc();
 289
 290     // Check whether counting numbers won't lead to problems
 291     for( size_t beta = 0; beta < nrIRs(); beta++ )
 292         DAI_ASSERT( nbIR(beta).size() + IR(beta).c() != 0.0 );
 293
 294     // Keep old beliefs to check convergence
 295     vector<Factor> oldBeliefsV;
 296     oldBeliefsV.reserve( nrVars() );
 297     for( size_t i = 0; i < nrVars(); i++ )
 298         oldBeliefsV.push_back( beliefV(i) );
 299     vector<Factor> oldBeliefsF;
 300     oldBeliefsF.reserve( nrFactors() );
 301     for( size_t I = 0; I < nrFactors(); I++ )
 302         oldBeliefsF.push_back( beliefF(I) );
 303
 304     // do several passes over the network until maximum number of iterations has
 305     // been reached or until the maximum belief difference is smaller than tolerance
 306     Real maxDiff = INFINITY;
 307     for( _iters = 0; _iters < props.maxiter && maxDiff > props.tol; _iters++ ) {
 308         for( size_t beta = 0; beta < nrIRs(); beta++ ) {
 309             foreach( const Neighbor &alpha, nbIR(beta) ) {
 310                 size_t _beta = alpha.dual;
 311                 muab( alpha, _beta ) = _Qa[alpha].marginal(IR(beta)) / muba(alpha,_beta);
 312                 /* TODO: INVESTIGATE THIS PROBLEM
 313                  *
 314                  * In some cases, the muab's can have very large entries because the muba's have very
 315                  * small entries. This may cause NANs later on (e.g., multiplying large quantities may
 316                  * result in +inf; normalization then tries to calculate inf / inf which is NAN).
 317                  * A fix of this problem would consist in normalizing the messages muab.
 318                  * However, it is not obvious whether this is a real solution, because it has a
 319                  * negative performance impact and the NAN's seem to be a symptom of a fundamental
 320                  * numerical unstability.
 321                  */
 322                  muab(alpha,_beta).normalize();
 323             }
 324
 325             Factor Qb_new;
 326             foreach( const Neighbor &alpha, nbIR(beta) ) {
 327                 size_t _beta = alpha.dual;
 328                 Qb_new *= muab(alpha,_beta) ^ (1 / (nbIR(beta).size() + IR(beta).c()));
 329             }
 330
 331             Qb_new.normalize();
 332             if( Qb_new.hasNaNs() ) {
 333                 // TODO: WHAT TO DO IN THIS CASE?
 334                 cerr << Name << "::doGBP:  Qb_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                 _Qb[beta] = Qb_new; // no damping for double loop
 344             else
 345                 _Qb[beta] = (Qb_new^(1.0 - props.damping)) * (_Qb[beta]^props.damping);
 346
 347             foreach( const Neighbor &alpha, nbIR(beta) ) {
 348                 size_t _beta = alpha.dual;
 349                 muba(alpha,_beta) = _Qb[beta] / muab(alpha,_beta);
 350
 351                 /* TODO: INVESTIGATE WHETHER THIS HACK (INVENTED BY KEES) TO PREVENT NANS MAKES SENSE
 352                  *
 353                  *   muba(beta,*alpha).makePositive(1e-100);
 354                  *
 355                  */
 356
 357                 Factor Qa_new = OR(alpha);
 358                 foreach( const Neighbor &gamma, nbOR(alpha) )
 359                     Qa_new *= muba(alpha,gamma.iter);
 360                 Qa_new ^= (1.0 / OR(alpha).c());
 361                 Qa_new.normalize();
 362                 if( Qa_new.hasNaNs() ) {
 363                     cerr << Name << "::doGBP:  Qa_new has NaNs!" << endl;
 364                     return 1.0;
 365                 }
 366                 /* TODO: WHAT IS THE PURPOSE OF THE FOLLOWING CODE?
 367                  *
 368                  *   _Qb[beta] = Qb_new.makeZero(1e-100);
 369                  */
 370
 371                 if( props.doubleloop || props.damping == 0.0 )
 372                     _Qa[alpha] = Qa_new; // no damping for double loop
 373                 else
 374                     // FIXME: GEOMETRIC DAMPING IS SLOW!
 375                 _Qa[alpha] = (Qa_new^(1.0 - props.damping)) * (_Qa[alpha]^props.damping);
 376             }
 377         }
 378
 379         // Calculate new single variable beliefs and compare with old ones
 380         maxDiff = -INFINITY;
 381         for( size_t i = 0; i < nrVars(); ++i ) {
 382             Factor b = beliefV(i);
 383             maxDiff = std::max( maxDiff, dist( b, oldBeliefsV[i], Prob::DISTLINF ) );
 384             oldBeliefsV[i] = b;
 385         }
 386         for( size_t I = 0; I < nrFactors(); ++I ) {
 387             Factor b = beliefF(I);
 388             maxDiff = std::max( maxDiff, dist( b, oldBeliefsF[I], Prob::DISTLINF ) );
 389             oldBeliefsF[I] = b;
 390         }
 391
 392         if( props.verbose >= 3 )
 393             cerr << Name << "::doGBP:  maxdiff " << maxDiff << " after " << _iters+1 << " passes" << endl;
 394     }
 395
 396     if( maxDiff > _maxdiff )
 397         _maxdiff = maxDiff;
 398
 399     if( props.verbose >= 1 ) {
 400         if( maxDiff > props.tol ) {
 401             if( props.verbose == 1 )
 402                 cerr << endl;
 403             cerr << Name << "::doGBP:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << maxDiff << endl;
 404         } else {
 405             if( props.verbose >= 2 )
 406                 cerr << Name << "::doGBP:  ";
 407             cerr << "converged in " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
 408         }
 409     }
 410
 411     return maxDiff;
 412 }
 413
 414
 415 Real HAK::doDoubleLoop() {
 416     if( props.verbose >= 1 )
 417         cerr << "Starting " << identify() << "...";
 418     if( props.verbose >= 3)
 419         cerr << endl;
 420
 421     double tic = toc();
 422
 423     // Save original outer regions
 424     vector<FRegion> org_ORs = ORs;
 425
 426     // Save original inner counting numbers and set negative counting numbers to zero
 427     vector<Real> org_IR_cs( nrIRs(), 0.0 );
 428     for( size_t beta = 0; beta < nrIRs(); beta++ ) {
 429         org_IR_cs[beta] = IR(beta).c();
 430         if( IR(beta).c() < 0.0 )
 431             IR(beta).c() = 0.0;
 432     }
 433
 434     // Keep old beliefs to check convergence
 435     vector<Factor> oldBeliefsV;
 436     oldBeliefsV.reserve( nrVars() );
 437     for( size_t i = 0; i < nrVars(); i++ )
 438         oldBeliefsV.push_back( beliefV(i) );
 439     vector<Factor> oldBeliefsF;
 440     oldBeliefsF.reserve( nrFactors() );
 441     for( size_t I = 0; I < nrFactors(); I++ )
 442         oldBeliefsF.push_back( beliefF(I) );
 443
 444     size_t outer_maxiter   = props.maxiter;
 445     Real   outer_tol       = props.tol;
 446     size_t outer_verbose   = props.verbose;
 447     Real   org_maxdiff     = _maxdiff;
 448
 449     // Set parameters for inner loop
 450     props.maxiter = 5;
 451     props.verbose = outer_verbose ? outer_verbose - 1 : 0;
 452
 453     size_t outer_iter = 0;
 454     size_t total_iter = 0;
 455     Real maxDiff = INFINITY;
 456     for( outer_iter = 0; outer_iter < outer_maxiter && maxDiff > outer_tol; outer_iter++ ) {
 457         // Calculate new outer regions
 458         for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
 459             OR(alpha) = org_ORs[alpha];
 460             foreach( const Neighbor &beta, nbOR(alpha) )
 461                 OR(alpha) *= _Qb[beta] ^ ((IR(beta).c() - org_IR_cs[beta]) / nbIR(beta).size());
 462         }
 463
 464         // Inner loop
 465         if( isnan( doGBP() ) )
 466             return 1.0;
 467
 468         // Calculate new single variable beliefs and compare with old ones
 469         maxDiff = -INFINITY;
 470         for( size_t i = 0; i < nrVars(); ++i ) {
 471             Factor b = beliefV(i);
 472             maxDiff = std::max( maxDiff, dist( b, oldBeliefsV[i], Prob::DISTLINF ) );
 473             oldBeliefsV[i] = b;
 474         }
 475         for( size_t I = 0; I < nrFactors(); ++I ) {
 476             Factor b = beliefF(I);
 477             maxDiff = std::max( maxDiff, dist( b, oldBeliefsF[I], Prob::DISTLINF ) );
 478             oldBeliefsF[I] = b;
 479         }
 480
 481         total_iter += Iterations();
 482
 483         if( props.verbose >= 3 )
 484             cerr << Name << "::doDoubleLoop:  maxdiff " << maxDiff << " after " << total_iter << " passes" << endl;
 485     }
 486
 487     // restore _maxiter, _verbose and _maxdiff
 488     props.maxiter = outer_maxiter;
 489     props.verbose = outer_verbose;
 490     _maxdiff = org_maxdiff;
 491
 492     _iters = total_iter;
 493     if( maxDiff > _maxdiff )
 494         _maxdiff = maxDiff;
 495
 496     // Restore original outer regions
 497     ORs = org_ORs;
 498
 499     // Restore original inner counting numbers
 500     for( size_t beta = 0; beta < nrIRs(); ++beta )
 501         IR(beta).c() = org_IR_cs[beta];
 502
 503     if( props.verbose >= 1 ) {
 504         if( maxDiff > props.tol ) {
 505             if( props.verbose == 1 )
 506                 cerr << endl;
 507                 cerr << Name << "::doDoubleLoop:  WARNING: not converged within " << outer_maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << maxDiff << endl;
 508             } else {
 509                 if( props.verbose >= 3 )
 510                     cerr << Name << "::doDoubleLoop:  ";
 511                 cerr << "converged in " << total_iter << " passes (" << toc() - tic << " seconds)." << endl;
 512             }
 513         }
 514
 515     return maxDiff;
 516 }
 517
 518
 519 Real HAK::run() {
 520     if( props.doubleloop )
 521         return doDoubleLoop();
 522     else
 523         return doGBP();
 524 }
 525
 526
 527 Factor HAK::belief( const VarSet &ns ) const {
 528     vector<Factor>::const_iterator beta;
 529     for( beta = _Qb.begin(); beta != _Qb.end(); beta++ )
 530         if( beta->vars() >> ns )
 531             break;
 532     if( beta != _Qb.end() )
 533         return( beta->marginal(ns) );
 534     else {
 535         vector<Factor>::const_iterator alpha;
 536         for( alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
 537             if( alpha->vars() >> ns )
 538                 break;
 539         if( alpha == _Qa.end() )
 540             DAI_THROW(BELIEF_NOT_AVAILABLE);
 541         return( alpha->marginal(ns) );
 542     }
 543 }
 544
 545
 546 vector<Factor> HAK::beliefs() const {
 547     vector<Factor> result;
 548     for( size_t beta = 0; beta < nrIRs(); beta++ )
 549         result.push_back( Qb(beta) );
 550     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
 551         result.push_back( Qa(alpha) );
 552     return result;
 553 }
 554
 555
 556 Real HAK::logZ() const {
 557     Real s = 0.0;
 558     for( size_t beta = 0; beta < nrIRs(); beta++ )
 559         s += IR(beta).c() * Qb(beta).entropy();
 560     for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
 561         s += OR(alpha).c() * Qa(alpha).entropy();
 562         s += (OR(alpha).log(true) * Qa(alpha)).sum();
 563     }
 564     return s;
 565 }
 566
 567
 568 } // end of namespace dai