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