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