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