Implemented various heuristics for choosing a variable elimination sequence in JTree
[libdai.git] / src / hak.cpp
index 7232615..81698dc 100644 (file)
@@ -24,6 +24,25 @@ using namespace std;
 const char *HAK::Name = "HAK";
 
 
 const char *HAK::Name = "HAK";
 
 
+/// Sets factor entries that lie between 0 and \a epsilon to \a epsilon
+template <class T>
+TFactor<T>& makePositive( TFactor<T> &f, T epsilon ) {
+    for( size_t t = 0; t < f.states(); t++ )
+        if( (0 < f[t]) && (f[t] < epsilon) )
+            f[t] = epsilon;
+    return f;
+}
+
+/// Sets factor entries that are smaller (in absolute value) than \a epsilon to 0
+template <class T>
+TFactor<T>& makeZero( TFactor<T> &f, T epsilon ) {
+    for( size_t t = 0; t < f.states(); t++ )
+        if( f[t] < epsilon && f[t] > -epsilon )
+            f[t] = 0;
+    return f;
+}
+
+
 void HAK::setProperties( const PropertySet &opts ) {
     DAI_ASSERT( opts.hasKey("tol") );
     DAI_ASSERT( opts.hasKey("maxiter") );
 void HAK::setProperties( const PropertySet &opts ) {
     DAI_ASSERT( opts.hasKey("tol") );
     DAI_ASSERT( opts.hasKey("maxiter") );
@@ -31,7 +50,7 @@ void HAK::setProperties( const PropertySet &opts ) {
     DAI_ASSERT( opts.hasKey("doubleloop") );
     DAI_ASSERT( opts.hasKey("clusters") );
 
     DAI_ASSERT( opts.hasKey("doubleloop") );
     DAI_ASSERT( opts.hasKey("clusters") );
 
-    props.tol = opts.getStringAs<double>("tol");
+    props.tol = opts.getStringAs<Real>("tol");
     props.maxiter = opts.getStringAs<size_t>("maxiter");
     props.verbose = opts.getStringAs<size_t>("verbose");
     props.doubleloop = opts.getStringAs<bool>("doubleloop");
     props.maxiter = opts.getStringAs<size_t>("maxiter");
     props.verbose = opts.getStringAs<size_t>("verbose");
     props.doubleloop = opts.getStringAs<bool>("doubleloop");
@@ -42,7 +61,7 @@ void HAK::setProperties( const PropertySet &opts ) {
     else
         DAI_ASSERT( props.clusters != Properties::ClustersType::LOOP );
     if( opts.hasKey("damping") )
     else
         DAI_ASSERT( props.clusters != Properties::ClustersType::LOOP );
     if( opts.hasKey("damping") )
-        props.damping = opts.getStringAs<double>("damping");
+        props.damping = opts.getStringAs<Real>("damping");
     else
         props.damping = 0.0;
     if( opts.hasKey("init") )
     else
         props.damping = 0.0;
     if( opts.hasKey("init") )
@@ -81,12 +100,12 @@ string HAK::printProperties() const {
 }
 
 
 }
 
 
-void HAK::constructMessages() {
+void HAK::construct() {
     // Create outer beliefs
     _Qa.clear();
     _Qa.reserve(nrORs());
     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
     // Create outer beliefs
     _Qa.clear();
     _Qa.reserve(nrORs());
     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
-        _Qa.push_back( Factor( OR(alpha).vars() ) );
+        _Qa.push_back( Factor( OR(alpha) ) );
 
     // Create inner beliefs
     _Qb.clear();
 
     // Create inner beliefs
     _Qb.clear();
@@ -115,16 +134,15 @@ void HAK::constructMessages() {
 HAK::HAK( const RegionGraph &rg, const PropertySet &opts ) : DAIAlgRG(rg), _Qa(), _Qb(), _muab(), _muba(), _maxdiff(0.0), _iters(0U), props() {
     setProperties( opts );
 
 HAK::HAK( const RegionGraph &rg, const PropertySet &opts ) : DAIAlgRG(rg), _Qa(), _Qb(), _muab(), _muba(), _maxdiff(0.0), _iters(0U), props() {
     setProperties( opts );
 
-    constructMessages();
+    construct();
 }
 
 
 void HAK::findLoopClusters( const FactorGraph & fg, std::set<VarSet> &allcl, VarSet newcl, const Var & root, size_t length, VarSet vars ) {
     for( VarSet::const_iterator in = vars.begin(); in != vars.end(); in++ ) {
         VarSet ind = fg.delta( fg.findVar( *in ) );
 }
 
 
 void HAK::findLoopClusters( const FactorGraph & fg, std::set<VarSet> &allcl, VarSet newcl, const Var & root, size_t length, VarSet vars ) {
     for( VarSet::const_iterator in = vars.begin(); in != vars.end(); in++ ) {
         VarSet ind = fg.delta( fg.findVar( *in ) );
-        if( (newcl.size()) >= 2 && ind.contains( root ) ) {
+        if( (newcl.size()) >= 2 && ind.contains( root ) )
             allcl.insert( newcl | *in );
             allcl.insert( newcl | *in );
-        }
         else if( length > 1 )
             findLoopClusters( fg, allcl, newcl | *in, root, length - 1, ind / newcl );
     }
         else if( length > 1 )
             findLoopClusters( fg, allcl, newcl | *in, root, length - 1, ind / newcl );
     }
@@ -137,9 +155,12 @@ HAK::HAK(const FactorGraph & fg, const PropertySet &opts) : DAIAlgRG(), _Qa(), _
     vector<VarSet> cl;
     if( props.clusters == Properties::ClustersType::MIN ) {
         cl = fg.Cliques();
     vector<VarSet> cl;
     if( props.clusters == Properties::ClustersType::MIN ) {
         cl = fg.Cliques();
+        constructCVM( fg, cl );
     } else if( props.clusters == Properties::ClustersType::DELTA ) {
     } else if( props.clusters == Properties::ClustersType::DELTA ) {
+        cl.reserve( fg.nrVars() );
         for( size_t i = 0; i < fg.nrVars(); i++ )
         for( size_t i = 0; i < fg.nrVars(); i++ )
-            cl.push_back(fg.Delta(i));
+            cl.push_back( fg.Delta(i) );
+        constructCVM( fg, cl );
     } else if( props.clusters == Properties::ClustersType::LOOP ) {
         cl = fg.Cliques();
         set<VarSet> scl;
     } else if( props.clusters == Properties::ClustersType::LOOP ) {
         cl = fg.Cliques();
         set<VarSet> scl;
@@ -155,12 +176,37 @@ HAK::HAK(const FactorGraph & fg, const PropertySet &opts) : DAIAlgRG(), _Qa(), _
             for( vector<VarSet>::const_iterator cli = cl.begin(); cli != cl.end(); cli++ )
                 cerr << *cli << endl;
         }
             for( vector<VarSet>::const_iterator cli = cl.begin(); cli != cl.end(); cli++ )
                 cerr << *cli << endl;
         }
+        constructCVM( fg, cl );
+    } else if( props.clusters == Properties::ClustersType::BETHE ) {
+        // build outer regions (the cliques)
+        cl = fg.Cliques();
+        size_t nrEdges = 0;
+        for( size_t c = 0; c < cl.size(); c++ )
+            nrEdges += cl[c].size();
+
+        // build inner regions (single variables)
+        vector<Region> irs;
+        irs.reserve( fg.nrVars() );
+        for( size_t i = 0; i < fg.nrVars(); i++ )
+            irs.push_back( Region( fg.var(i), 1.0 ) );
+
+        // build edges (an outer and inner region are connected if the outer region contains the inner one)
+        // and calculate counting number for inner regions
+        vector<std::pair<size_t, size_t> > edges;
+        edges.reserve( nrEdges );
+        for( size_t c = 0; c < cl.size(); c++ )
+            for( size_t i = 0; i < irs.size(); i++ )
+                if( cl[c] >> irs[i] ) {
+                    edges.push_back( make_pair( c, i ) );
+                    irs[i].c() -= 1.0;
+                }
+
+        // build region graph
+        RegionGraph::construct( fg, cl, irs, edges );
     } else
         DAI_THROW(UNKNOWN_ENUM_VALUE);
 
     } else
         DAI_THROW(UNKNOWN_ENUM_VALUE);
 
-    RegionGraph rg(fg,cl);
-    RegionGraph::operator=(rg);
-    constructMessages();
+    construct();
 
     if( props.verbose >= 3 )
         cerr << Name << " regiongraph: " << *this << endl;
 
     if( props.verbose >= 3 )
         cerr << Name << " regiongraph: " << *this << endl;
@@ -173,12 +219,14 @@ string HAK::identify() const {
 
 
 void HAK::init( const VarSet &ns ) {
 
 
 void HAK::init( const VarSet &ns ) {
-    for( vector<Factor>::iterator alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
-        if( alpha->vars().intersects( ns ) ) {
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ )
+        if( _Qa[alpha].vars().intersects( ns ) ) {
             if( props.init == Properties::InitType::UNIFORM )
             if( props.init == Properties::InitType::UNIFORM )
-                alpha->fill( 1.0 / alpha->states() );
+                _Qa[alpha].setUniform();
             else
             else
-                alpha->randomize();
+                _Qa[alpha].randomize();
+            _Qa[alpha] *= OR(alpha);
+            _Qa[alpha].normalize();
         }
 
     for( size_t beta = 0; beta < nrIRs(); beta++ )
         }
 
     for( size_t beta = 0; beta < nrIRs(); beta++ )
@@ -202,24 +250,27 @@ void HAK::init( const VarSet &ns ) {
 
 
 void HAK::init() {
 
 
 void HAK::init() {
-    for( vector<Factor>::iterator alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
         if( props.init == Properties::InitType::UNIFORM )
         if( props.init == Properties::InitType::UNIFORM )
-            alpha->fill( 1.0 / alpha->states() );
+            _Qa[alpha].setUniform();
         else
         else
-            alpha->randomize();
+            _Qa[alpha].randomize();
+        _Qa[alpha] *= OR(alpha);
+        _Qa[alpha].normalize();
+    }
 
 
-    for( vector<Factor>::iterator beta = _Qb.begin(); beta != _Qb.end(); beta++ )
+    for( size_t beta = 0; beta < nrIRs(); beta++ )
         if( props.init == Properties::InitType::UNIFORM )
         if( props.init == Properties::InitType::UNIFORM )
-            beta->fill( 1.0 / beta->states() );
+            _Qb[beta].setUniform();
         else
         else
-            beta->randomize();
+            _Qb[beta].randomize();
 
     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
         foreach( const Neighbor &beta, nbOR(alpha) ) {
             size_t _beta = beta.iter;
             if( props.init == Properties::InitType::UNIFORM ) {
 
     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
         foreach( const Neighbor &beta, nbOR(alpha) ) {
             size_t _beta = beta.iter;
             if( props.init == Properties::InitType::UNIFORM ) {
-                muab( alpha, _beta ).fill( 1.0 / muab( alpha, _beta ).states() );
-                muba( alpha, _beta ).fill( 1.0 / muab( alpha, _beta ).states() );
+                muab( alpha, _beta ).setUniform();
+                muba( alpha, _beta ).setUniform();
             } else {
                 muab( alpha, _beta ).randomize();
                 muba( alpha, _beta ).randomize();
             } else {
                 muab( alpha, _beta ).randomize();
                 muba( alpha, _beta ).randomize();
@@ -228,7 +279,7 @@ void HAK::init() {
 }
 
 
 }
 
 
-double HAK::doGBP() {
+Real HAK::doGBP() {
     if( props.verbose >= 1 )
         cerr << "Starting " << identify() << "...";
     if( props.verbose >= 3)
     if( props.verbose >= 1 )
         cerr << "Starting " << identify() << "...";
     if( props.verbose >= 3)
@@ -241,17 +292,19 @@ double HAK::doGBP() {
         DAI_ASSERT( nbIR(beta).size() + IR(beta).c() != 0.0 );
 
     // Keep old beliefs to check convergence
         DAI_ASSERT( nbIR(beta).size() + IR(beta).c() != 0.0 );
 
     // Keep old beliefs to check convergence
-    vector<Factor> old_beliefs;
-    old_beliefs.reserve( nrVars() );
+    vector<Factor> oldBeliefsV;
+    oldBeliefsV.reserve( nrVars() );
     for( size_t i = 0; i < nrVars(); i++ )
     for( size_t i = 0; i < nrVars(); i++ )
-        old_beliefs.push_back( belief( var(i) ) );
-
-    // Differences in single node beliefs
-    Diffs diffs(nrVars(), 1.0);
+        oldBeliefsV.push_back( beliefV(i) );
+    vector<Factor> oldBeliefsF;
+    oldBeliefsF.reserve( nrFactors() );
+    for( size_t I = 0; I < nrFactors(); I++ )
+        oldBeliefsF.push_back( beliefF(I) );
 
     // do several passes over the network until maximum number of iterations has
     // been reached or until the maximum belief difference is smaller than tolerance
 
     // do several passes over the network until maximum number of iterations has
     // been reached or until the maximum belief difference is smaller than tolerance
-    for( _iters = 0; _iters < props.maxiter && diffs.maxDiff() > props.tol; _iters++ ) {
+    Real maxDiff = INFINITY;
+    for( _iters = 0; _iters < props.maxiter && maxDiff > props.tol; _iters++ ) {
         for( size_t beta = 0; beta < nrIRs(); beta++ ) {
             foreach( const Neighbor &alpha, nbIR(beta) ) {
                 size_t _beta = alpha.dual;
         for( size_t beta = 0; beta < nrIRs(); beta++ ) {
             foreach( const Neighbor &alpha, nbIR(beta) ) {
                 size_t _beta = alpha.dual;
@@ -324,24 +377,30 @@ double HAK::doGBP() {
         }
 
         // Calculate new single variable beliefs and compare with old ones
         }
 
         // Calculate new single variable beliefs and compare with old ones
-        for( size_t i = 0; i < nrVars(); i++ ) {
-            Factor new_belief = belief( var( i ) );
-            diffs.push( dist( new_belief, old_beliefs[i], Prob::DISTLINF ) );
-            old_beliefs[i] = new_belief;
+        maxDiff = -INFINITY;
+        for( size_t i = 0; i < nrVars(); ++i ) {
+            Factor b = beliefV(i);
+            maxDiff = std::max( maxDiff, dist( b, oldBeliefsV[i], Prob::DISTLINF ) );
+            oldBeliefsV[i] = b;
+        }
+        for( size_t I = 0; I < nrFactors(); ++I ) {
+            Factor b = beliefF(I);
+            maxDiff = std::max( maxDiff, dist( b, oldBeliefsF[I], Prob::DISTLINF ) );
+            oldBeliefsF[I] = b;
         }
 
         if( props.verbose >= 3 )
         }
 
         if( props.verbose >= 3 )
-            cerr << Name << "::doGBP:  maxdiff " << diffs.maxDiff() << " after " << _iters+1 << " passes" << endl;
+            cerr << Name << "::doGBP:  maxdiff " << maxDiff << " after " << _iters+1 << " passes" << endl;
     }
 
     }
 
-    if( diffs.maxDiff() > _maxdiff )
-        _maxdiff = diffs.maxDiff();
+    if( maxDiff > _maxdiff )
+        _maxdiff = maxDiff;
 
     if( props.verbose >= 1 ) {
 
     if( props.verbose >= 1 ) {
-        if( diffs.maxDiff() > props.tol ) {
+        if( maxDiff > props.tol ) {
             if( props.verbose == 1 )
                 cerr << endl;
             if( props.verbose == 1 )
                 cerr << endl;
-            cerr << Name << "::doGBP:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << diffs.maxDiff() << endl;
+            cerr << Name << "::doGBP:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << maxDiff << endl;
         } else {
             if( props.verbose >= 2 )
                 cerr << Name << "::doGBP:  ";
         } else {
             if( props.verbose >= 2 )
                 cerr << Name << "::doGBP:  ";
@@ -349,11 +408,11 @@ double HAK::doGBP() {
         }
     }
 
         }
     }
 
-    return diffs.maxDiff();
+    return maxDiff;
 }
 
 
 }
 
 
-double HAK::doDoubleLoop() {
+Real HAK::doDoubleLoop() {
     if( props.verbose >= 1 )
         cerr << "Starting " << identify() << "...";
     if( props.verbose >= 3)
     if( props.verbose >= 1 )
         cerr << "Starting " << identify() << "...";
     if( props.verbose >= 3)
@@ -365,7 +424,7 @@ double HAK::doDoubleLoop() {
     vector<FRegion> org_ORs = ORs;
 
     // Save original inner counting numbers and set negative counting numbers to zero
     vector<FRegion> org_ORs = ORs;
 
     // Save original inner counting numbers and set negative counting numbers to zero
-    vector<double> org_IR_cs( nrIRs(), 0.0 );
+    vector<Real> org_IR_cs( nrIRs(), 0.0 );
     for( size_t beta = 0; beta < nrIRs(); beta++ ) {
         org_IR_cs[beta] = IR(beta).c();
         if( IR(beta).c() < 0.0 )
     for( size_t beta = 0; beta < nrIRs(); beta++ ) {
         org_IR_cs[beta] = IR(beta).c();
         if( IR(beta).c() < 0.0 )
@@ -373,18 +432,19 @@ double HAK::doDoubleLoop() {
     }
 
     // Keep old beliefs to check convergence
     }
 
     // Keep old beliefs to check convergence
-    vector<Factor> old_beliefs;
-    old_beliefs.reserve( nrVars() );
+    vector<Factor> oldBeliefsV;
+    oldBeliefsV.reserve( nrVars() );
     for( size_t i = 0; i < nrVars(); i++ )
     for( size_t i = 0; i < nrVars(); i++ )
-        old_beliefs.push_back( belief( var(i) ) );
-
-    // Differences in single node beliefs
-    Diffs diffs(nrVars(), 1.0);
+        oldBeliefsV.push_back( beliefV(i) );
+    vector<Factor> oldBeliefsF;
+    oldBeliefsF.reserve( nrFactors() );
+    for( size_t I = 0; I < nrFactors(); I++ )
+        oldBeliefsF.push_back( beliefF(I) );
 
     size_t outer_maxiter   = props.maxiter;
 
     size_t outer_maxiter   = props.maxiter;
-    double outer_tol       = props.tol;
+    Real   outer_tol       = props.tol;
     size_t outer_verbose   = props.verbose;
     size_t outer_verbose   = props.verbose;
-    double org_maxdiff     = _maxdiff;
+    Real   org_maxdiff     = _maxdiff;
 
     // Set parameters for inner loop
     props.maxiter = 5;
 
     // Set parameters for inner loop
     props.maxiter = 5;
@@ -392,7 +452,8 @@ double HAK::doDoubleLoop() {
 
     size_t outer_iter = 0;
     size_t total_iter = 0;
 
     size_t outer_iter = 0;
     size_t total_iter = 0;
-    for( outer_iter = 0; outer_iter < outer_maxiter && diffs.maxDiff() > outer_tol; outer_iter++ ) {
+    Real maxDiff = INFINITY;
+    for( outer_iter = 0; outer_iter < outer_maxiter && maxDiff > outer_tol; outer_iter++ ) {
         // Calculate new outer regions
         for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
             OR(alpha) = org_ORs[alpha];
         // Calculate new outer regions
         for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
             OR(alpha) = org_ORs[alpha];
@@ -405,16 +466,22 @@ double HAK::doDoubleLoop() {
             return 1.0;
 
         // Calculate new single variable beliefs and compare with old ones
             return 1.0;
 
         // Calculate new single variable beliefs and compare with old ones
+        maxDiff = -INFINITY;
         for( size_t i = 0; i < nrVars(); ++i ) {
         for( size_t i = 0; i < nrVars(); ++i ) {
-            Factor new_belief = belief( var( i ) );
-            diffs.push( dist( new_belief, old_beliefs[i], Prob::DISTLINF ) );
-            old_beliefs[i] = new_belief;
+            Factor b = beliefV(i);
+            maxDiff = std::max( maxDiff, dist( b, oldBeliefsV[i], Prob::DISTLINF ) );
+            oldBeliefsV[i] = b;
+        }
+        for( size_t I = 0; I < nrFactors(); ++I ) {
+            Factor b = beliefF(I);
+            maxDiff = std::max( maxDiff, dist( b, oldBeliefsF[I], Prob::DISTLINF ) );
+            oldBeliefsF[I] = b;
         }
 
         total_iter += Iterations();
 
         if( props.verbose >= 3 )
         }
 
         total_iter += Iterations();
 
         if( props.verbose >= 3 )
-            cerr << Name << "::doDoubleLoop:  maxdiff " << diffs.maxDiff() << " after " << total_iter << " passes" << endl;
+            cerr << Name << "::doDoubleLoop:  maxdiff " << maxDiff << " after " << total_iter << " passes" << endl;
     }
 
     // restore _maxiter, _verbose and _maxdiff
     }
 
     // restore _maxiter, _verbose and _maxdiff
@@ -423,8 +490,8 @@ double HAK::doDoubleLoop() {
     _maxdiff = org_maxdiff;
 
     _iters = total_iter;
     _maxdiff = org_maxdiff;
 
     _iters = total_iter;
-    if( diffs.maxDiff() > _maxdiff )
-        _maxdiff = diffs.maxDiff();
+    if( maxDiff > _maxdiff )
+        _maxdiff = maxDiff;
 
     // Restore original outer regions
     ORs = org_ORs;
 
     // Restore original outer regions
     ORs = org_ORs;
@@ -434,10 +501,10 @@ double HAK::doDoubleLoop() {
         IR(beta).c() = org_IR_cs[beta];
 
     if( props.verbose >= 1 ) {
         IR(beta).c() = org_IR_cs[beta];
 
     if( props.verbose >= 1 ) {
-        if( diffs.maxDiff() > props.tol ) {
+        if( maxDiff > props.tol ) {
             if( props.verbose == 1 )
                 cerr << endl;
             if( props.verbose == 1 )
                 cerr << endl;
-                cerr << Name << "::doDoubleLoop:  WARNING: not converged within " << outer_maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << diffs.maxDiff() << endl;
+                cerr << Name << "::doDoubleLoop:  WARNING: not converged within " << outer_maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << maxDiff << endl;
             } else {
                 if( props.verbose >= 3 )
                     cerr << Name << "::doDoubleLoop:  ";
             } else {
                 if( props.verbose >= 3 )
                     cerr << Name << "::doDoubleLoop:  ";
@@ -445,11 +512,11 @@ double HAK::doDoubleLoop() {
             }
         }
 
             }
         }
 
-    return diffs.maxDiff();
+    return maxDiff;
 }
 
 
 }
 
 
-double HAK::run() {
+Real HAK::run() {
     if( props.doubleloop )
         return doDoubleLoop();
     else
     if( props.doubleloop )
         return doDoubleLoop();
     else
@@ -469,17 +536,13 @@ Factor HAK::belief( const VarSet &ns ) const {
         for( alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
             if( alpha->vars() >> ns )
                 break;
         for( alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
             if( alpha->vars() >> ns )
                 break;
-        DAI_ASSERT( alpha != _Qa.end() );
+        if( alpha == _Qa.end() )
+            DAI_THROW(BELIEF_NOT_AVAILABLE);
         return( alpha->marginal(ns) );
     }
 }
 
 
         return( alpha->marginal(ns) );
     }
 }
 
 
-Factor HAK::belief( const Var &n ) const {
-    return belief( (VarSet)n );
-}
-
-
 vector<Factor> HAK::beliefs() const {
     vector<Factor> result;
     for( size_t beta = 0; beta < nrIRs(); beta++ )
 vector<Factor> HAK::beliefs() const {
     vector<Factor> result;
     for( size_t beta = 0; beta < nrIRs(); beta++ )