src/hak.cpp

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