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