Merge branch 'master' of git@git.tuebingen.mpg.de:libdai
[libdai.git] / src / hak.cpp
index aaa673b..49903b8 100644 (file)
@@ -23,6 +23,7 @@
 #include <dai/hak.h>
 #include <dai/util.h>
 #include <dai/diffs.h>
 #include <dai/hak.h>
 #include <dai/util.h>
 #include <dai/diffs.h>
+#include <dai/exceptions.h>
 
 
 namespace dai {
 
 
 namespace dai {
@@ -34,30 +35,35 @@ using namespace std;
 const char *HAK::Name = "HAK";
 
 
 const char *HAK::Name = "HAK";
 
 
-bool HAK::checkProperties() {
-    if( !HasProperty("tol") )
-        return false;
-    if (!HasProperty("maxiter") )
-        return false;
-    if (!HasProperty("verbose") )
-        return false;
-    if( !HasProperty("doubleloop") )
-        return false;
-    if( !HasProperty("clusters") )
-        return false;
+void HAK::setProperties( const PropertySet &opts ) {
+    assert( opts.hasKey("tol") );
+    assert( opts.hasKey("maxiter") );
+    assert( opts.hasKey("verbose") );
+    assert( opts.hasKey("doubleloop") );
+    assert( opts.hasKey("clusters") );
     
     
-    ConvertPropertyTo<double>("tol");
-    ConvertPropertyTo<size_t>("maxiter");
-    ConvertPropertyTo<size_t>("verbose");
-    ConvertPropertyTo<bool>("doubleloop");
-    ConvertPropertyTo<ClustersType>("clusters");
-
-    if( HasProperty("loopdepth") )
-        ConvertPropertyTo<size_t>("loopdepth");
-    else if( Clusters() == ClustersType::LOOP )
-        return false;
-
-    return true;
+    props.tol = opts.getStringAs<double>("tol");
+    props.maxiter = opts.getStringAs<size_t>("maxiter");
+    props.verbose = opts.getStringAs<size_t>("verbose");
+    props.doubleloop = opts.getStringAs<bool>("doubleloop");
+    props.clusters = opts.getStringAs<Properties::ClustersType>("clusters");
+
+    if( opts.hasKey("loopdepth") )
+        props.loopdepth = opts.getStringAs<size_t>("loopdepth");
+    else
+        assert( props.clusters != Properties::ClustersType::LOOP );
+}
+
+
+PropertySet HAK::getProperties() const {
+    PropertySet opts;
+    opts.Set( "tol", props.tol );
+    opts.Set( "maxiter", props.maxiter );
+    opts.Set( "verbose", props.verbose );
+    opts.Set( "doubleloop", props.doubleloop );
+    opts.Set( "clusters", props.clusters );
+    opts.Set( "loopdepth", props.loopdepth );
+    return opts;
 }
 
 
 }
 
 
@@ -92,14 +98,14 @@ void HAK::constructMessages() {
 }
 
 
 }
 
 
-HAK::HAK(const RegionGraph & rg, const Properties &opts) : DAIAlgRG(rg, opts) {
-    assert( checkProperties() );
+HAK::HAK(const RegionGraph & rg, const PropertySet &opts ) : DAIAlgRG(rg) {
+    setProperties( opts );
 
     constructMessages();
 }
 
 
 
     constructMessages();
 }
 
 
-void HAK::findLoopClusters( const FactorGraph & fg, set<VarSet> &allcl, VarSet newcl, const Var & root, size_t length, VarSet vars ) {
+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 >> root) ) {
     for( VarSet::const_iterator in = vars.begin(); in != vars.end(); in++ ) {
         VarSet ind = fg.delta( fg.findVar( *in ) );
         if( (newcl.size()) >= 2 && (ind >> root) ) {
@@ -111,91 +117,89 @@ void HAK::findLoopClusters( const FactorGraph & fg, set<VarSet> &allcl, VarSet n
 }
 
 
 }
 
 
-HAK::HAK(const FactorGraph & fg, const Properties &opts) : DAIAlgRG(opts) {
-    assert( checkProperties() );
+HAK::HAK(const FactorGraph & fg, const PropertySet &opts) : DAIAlgRG(), props(), maxdiff(0.0) {
+    setProperties( opts );
 
     vector<VarSet> cl;
 
     vector<VarSet> cl;
-    if( Clusters() == ClustersType::MIN ) {
+    if( props.clusters == Properties::ClustersType::MIN ) {
         cl = fg.Cliques();
         cl = fg.Cliques();
-    } else if( Clusters() == ClustersType::DELTA ) {
+    } else if( props.clusters == Properties::ClustersType::DELTA ) {
         for( size_t i = 0; i < fg.nrVars(); i++ )
             cl.push_back(fg.Delta(i)); 
         for( size_t i = 0; i < fg.nrVars(); i++ )
             cl.push_back(fg.Delta(i)); 
-    } else if( Clusters() == ClustersType::LOOP ) {
+    } else if( props.clusters == Properties::ClustersType::LOOP ) {
         cl = fg.Cliques();
         set<VarSet> scl;
         for( size_t i0 = 0; i0 < fg.nrVars(); i0++ ) {
             VarSet i0d = fg.delta(i0);
         cl = fg.Cliques();
         set<VarSet> scl;
         for( size_t i0 = 0; i0 < fg.nrVars(); i0++ ) {
             VarSet i0d = fg.delta(i0);
-            if( LoopDepth() > 1 )
-                findLoopClusters( fg, scl, fg.var(i0), fg.var(i0), LoopDepth() - 1, fg.delta(i0) );
+            if( props.loopdepth > 1 )
+                findLoopClusters( fg, scl, fg.var(i0), fg.var(i0), props.loopdepth - 1, fg.delta(i0) );
         }
         for( set<VarSet>::const_iterator c = scl.begin(); c != scl.end(); c++ )
             cl.push_back(*c);
         }
         for( set<VarSet>::const_iterator c = scl.begin(); c != scl.end(); c++ )
             cl.push_back(*c);
-        if( Verbose() >= 3 ) {
+        if( props.verbose >= 3 ) {
             cout << "HAK uses the following clusters: " << endl;
             for( vector<VarSet>::const_iterator cli = cl.begin(); cli != cl.end(); cli++ )
                 cout << *cli << endl;
         }
     } else
             cout << "HAK uses the following clusters: " << endl;
             for( vector<VarSet>::const_iterator cli = cl.begin(); cli != cl.end(); cli++ )
                 cout << *cli << endl;
         }
     } else
-        throw "Invalid Clusters type";
+        DAI_THROW(INTERNAL_ERROR);
 
     RegionGraph rg(fg,cl);
     RegionGraph::operator=(rg);
     constructMessages();
 
 
     RegionGraph rg(fg,cl);
     RegionGraph::operator=(rg);
     constructMessages();
 
-    if( Verbose() >= 3 )
+    if( props.verbose >= 3 )
         cout << "HAK regiongraph: " << *this << endl;
 }
 
 
 string HAK::identify() const { 
     stringstream result (stringstream::out);
         cout << "HAK regiongraph: " << *this << endl;
 }
 
 
 string HAK::identify() const { 
     stringstream result (stringstream::out);
-    result << Name << GetProperties();
+    result << Name << getProperties();
     return result.str();
 }
 
 
 void HAK::init( const VarSet &ns ) {
     for( vector<Factor>::iterator alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
     return result.str();
 }
 
 
 void HAK::init( const VarSet &ns ) {
     for( vector<Factor>::iterator alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
-        if( alpha->vars() && ns )
-            alpha->fill( 1.0 / alpha->stateSpace() );
+        if( alpha->vars().intersects( ns ) )
+            alpha->fill( 1.0 / alpha->states() );
 
     for( size_t beta = 0; beta < nrIRs(); beta++ )
 
     for( size_t beta = 0; beta < nrIRs(); beta++ )
-        if( IR(beta) && ns ) {
-            _Qb[beta].fill( 1.0 / IR(beta).stateSpace() );
+        if( IR(beta).intersects( ns ) ) {
+            _Qb[beta].fill( 1.0 );
             foreach( const Neighbor &alpha, nbIR(beta) ) {
                 size_t _beta = alpha.dual;
             foreach( const Neighbor &alpha, nbIR(beta) ) {
                 size_t _beta = alpha.dual;
-                muab( alpha, _beta ).fill( 1.0 / IR(beta).stateSpace() );
-                muba( alpha, _beta ).fill( 1.0 / IR(beta).stateSpace() );
+                muab( alpha, _beta ).fill( 1.0 / IR(beta).states() );
+                muba( alpha, _beta ).fill( 1.0 / IR(beta).states() );
             }
         }
 }
 
 
 void HAK::init() {
             }
         }
 }
 
 
 void HAK::init() {
-    assert( checkProperties() );
-
     for( vector<Factor>::iterator alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
     for( vector<Factor>::iterator alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
-        alpha->fill( 1.0 / alpha->stateSpace() );
+        alpha->fill( 1.0 / alpha->states() );
 
     for( vector<Factor>::iterator beta = _Qb.begin(); beta != _Qb.end(); beta++ )
 
     for( vector<Factor>::iterator beta = _Qb.begin(); beta != _Qb.end(); beta++ )
-        beta->fill( 1.0 / beta->stateSpace() );
+        beta->fill( 1.0 / beta->states() );
 
     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
         foreach( const Neighbor &beta, nbOR(alpha) ) {
             size_t _beta = beta.iter;
 
     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
         foreach( const Neighbor &beta, nbOR(alpha) ) {
             size_t _beta = beta.iter;
-            muab( alpha, _beta ).fill( 1.0 / muab( alpha, _beta ).stateSpace() );
-            muba( alpha, _beta ).fill( 1.0 / muab( alpha, _beta ).stateSpace() );
+            muab( alpha, _beta ).fill( 1.0 / muab( alpha, _beta ).states() );
+            muba( alpha, _beta ).fill( 1.0 / muab( alpha, _beta ).states() );
         }
 }
 
 
 double HAK::doGBP() {
         }
 }
 
 
 double HAK::doGBP() {
-    if( Verbose() >= 1 )
+    if( props.verbose >= 1 )
         cout << "Starting " << identify() << "...";
         cout << "Starting " << identify() << "...";
-    if( Verbose() >= 3)
+    if( props.verbose >= 3)
         cout << endl;
 
         cout << endl;
 
-    clock_t tic = toc();
+    double tic = toc();
 
     // Check whether counting numbers won't lead to problems
     for( size_t beta = 0; beta < nrIRs(); beta++ )
 
     // Check whether counting numbers won't lead to problems
     for( size_t beta = 0; beta < nrIRs(); beta++ )
@@ -213,7 +217,7 @@ double HAK::doGBP() {
     size_t iter = 0;
     // do several passes over the network until maximum number of iterations has
     // been reached or until the maximum belief difference is smaller than tolerance
     size_t iter = 0;
     // 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( iter = 0; iter < MaxIter() && diffs.max() > Tol(); iter++ ) {
+    for( iter = 0; iter < props.maxiter && diffs.maxDiff() > props.tol; iter++ ) {
         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;
@@ -226,10 +230,10 @@ double HAK::doGBP() {
                 Qb_new *= muab(alpha,_beta) ^ (1 / (nbIR(beta).size() + IR(beta).c()));
             }
 
                 Qb_new *= muab(alpha,_beta) ^ (1 / (nbIR(beta).size() + IR(beta).c()));
             }
 
-            Qb_new.normalize( _normtype );
+            Qb_new.normalize( Prob::NORMPROB );
             if( Qb_new.hasNaNs() ) {
                 cout << "HAK::doGBP:  Qb_new has NaNs!" << endl;
             if( Qb_new.hasNaNs() ) {
                 cout << "HAK::doGBP:  Qb_new has NaNs!" << endl;
-                return NAN;
+                return 1.0;
             }
 //          _Qb[beta] = Qb_new.makeZero(1e-100);    // damping?
             _Qb[beta] = Qb_new;
             }
 //          _Qb[beta] = Qb_new.makeZero(1e-100);    // damping?
             _Qb[beta] = Qb_new;
@@ -243,10 +247,10 @@ double HAK::doGBP() {
                 foreach( const Neighbor &gamma, nbOR(alpha) )
                     Qa_new *= muba(alpha,gamma.iter);
                 Qa_new ^= (1.0 / OR(alpha).c());
                 foreach( const Neighbor &gamma, nbOR(alpha) )
                     Qa_new *= muba(alpha,gamma.iter);
                 Qa_new ^= (1.0 / OR(alpha).c());
-                Qa_new.normalize( _normtype );
+                Qa_new.normalize( Prob::NORMPROB );
                 if( Qa_new.hasNaNs() ) {
                     cout << "HAK::doGBP:  Qa_new has NaNs!" << endl;
                 if( Qa_new.hasNaNs() ) {
                     cout << "HAK::doGBP:  Qa_new has NaNs!" << endl;
-                    return NAN;
+                    return 1.0;
                 }
 //              _Qa[alpha] = Qa_new.makeZero(1e-100); // damping?
                 _Qa[alpha] = Qa_new;
                 }
 //              _Qa[alpha] = Qa_new.makeZero(1e-100); // damping?
                 _Qa[alpha] = Qa_new;
@@ -260,35 +264,36 @@ double HAK::doGBP() {
             old_beliefs[i] = new_belief;
         }
 
             old_beliefs[i] = new_belief;
         }
 
-        if( Verbose() >= 3 )
-            cout << "HAK::doGBP:  maxdiff " << diffs.max() << " after " << iter+1 << " passes" << endl;
+        if( props.verbose >= 3 )
+            cout << "HAK::doGBP:  maxdiff " << diffs.maxDiff() << " after " << iter+1 << " passes" << endl;
     }
 
     }
 
-    updateMaxDiff( diffs.max() );
+    if( diffs.maxDiff() > maxdiff )
+        maxdiff = diffs.maxDiff();
 
 
-    if( Verbose() >= 1 ) {
-        if( diffs.max() > Tol() ) {
-            if( Verbose() == 1 )
+    if( props.verbose >= 1 ) {
+        if( diffs.maxDiff() > props.tol ) {
+            if( props.verbose == 1 )
                 cout << endl;
                 cout << endl;
-            cout << "HAK::doGBP:  WARNING: not converged within " << MaxIter() << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.max() << endl;
+            cout << "HAK::doGBP:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.maxDiff() << endl;
         } else {
         } else {
-            if( Verbose() >= 2 )
+            if( props.verbose >= 2 )
                 cout << "HAK::doGBP:  ";
             cout << "converged in " << iter << " passes (" << toc() - tic << " clocks)." << endl;
         }
     }
 
                 cout << "HAK::doGBP:  ";
             cout << "converged in " << iter << " passes (" << toc() - tic << " clocks)." << endl;
         }
     }
 
-    return diffs.max();
+    return diffs.maxDiff();
 }
 
 
 double HAK::doDoubleLoop() {
 }
 
 
 double HAK::doDoubleLoop() {
-    if( Verbose() >= 1 )
+    if( props.verbose >= 1 )
         cout << "Starting " << identify() << "...";
         cout << "Starting " << identify() << "...";
-    if( Verbose() >= 3)
+    if( props.verbose >= 3)
         cout << endl;
 
         cout << endl;
 
-    clock_t tic = toc();
+    double tic = toc();
 
     // Save original outer regions
     vector<FRegion> org_ORs = ORs;
 
     // Save original outer regions
     vector<FRegion> org_ORs = ORs;
@@ -310,17 +315,17 @@ double HAK::doDoubleLoop() {
     // Differences in single node beliefs
     Diffs diffs(nrVars(), 1.0);
 
     // Differences in single node beliefs
     Diffs diffs(nrVars(), 1.0);
 
-    size_t  outer_maxiter   = MaxIter();
-    double  outer_tol       = Tol();
-    size_t  outer_verbose   = Verbose();
-    double  org_maxdiff     = MaxDiff();
+    size_t  outer_maxiter   = props.maxiter;
+    double  outer_tol       = props.tol;
+    size_t  outer_verbose   = props.verbose;
+    double  org_maxdiff     = maxdiff;
 
     // Set parameters for inner loop
 
     // Set parameters for inner loop
-    MaxIter( 5 );
-    Verbose( outer_verbose ? outer_verbose - 1 : 0 );
+    props.maxiter = 5;
+    props.verbose = outer_verbose ? outer_verbose - 1 : 0;
 
     size_t outer_iter = 0;
 
     size_t outer_iter = 0;
-    for( outer_iter = 0; outer_iter < outer_maxiter && diffs.max() > outer_tol; outer_iter++ ) {
+    for( outer_iter = 0; outer_iter < outer_maxiter && diffs.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];
@@ -330,7 +335,7 @@ double HAK::doDoubleLoop() {
 
         // Inner loop
         if( isnan( doGBP() ) )
 
         // Inner loop
         if( isnan( doGBP() ) )
-            return NAN;
+            return 1.0;
 
         // Calculate new single variable beliefs and compare with old ones
         for( size_t i = 0; i < nrVars(); ++i ) {
 
         // Calculate new single variable beliefs and compare with old ones
         for( size_t i = 0; i < nrVars(); ++i ) {
@@ -339,16 +344,17 @@ double HAK::doDoubleLoop() {
             old_beliefs[i] = new_belief;
         }
 
             old_beliefs[i] = new_belief;
         }
 
-        if( Verbose() >= 3 )
-            cout << "HAK::doDoubleLoop:  maxdiff " << diffs.max() << " after " << outer_iter+1 << " passes" << endl;
+        if( props.verbose >= 3 )
+            cout << "HAK::doDoubleLoop:  maxdiff " << diffs.maxDiff() << " after " << outer_iter+1 << " passes" << endl;
     }
 
     // restore _maxiter, _verbose and _maxdiff
     }
 
     // restore _maxiter, _verbose and _maxdiff
-    MaxIter( outer_maxiter );
-    Verbose( outer_verbose );
-    MaxDiff( org_maxdiff );
+    props.maxiter = outer_maxiter;
+    props.verbose = outer_verbose;
+    maxdiff = org_maxdiff;
 
 
-    updateMaxDiff( diffs.max() );
+    if( diffs.maxDiff() > maxdiff )
+        maxdiff = diffs.maxDiff();
 
     // Restore original outer regions
     ORs = org_ORs;
 
     // Restore original outer regions
     ORs = org_ORs;
@@ -357,24 +363,24 @@ double HAK::doDoubleLoop() {
     for( size_t beta = 0; beta < nrIRs(); ++beta )
         IR(beta).c() = org_IR_cs[beta];
 
     for( size_t beta = 0; beta < nrIRs(); ++beta )
         IR(beta).c() = org_IR_cs[beta];
 
-    if( Verbose() >= 1 ) {
-        if( diffs.max() > Tol() ) {
-            if( Verbose() == 1 )
+    if( props.verbose >= 1 ) {
+        if( diffs.maxDiff() > props.tol ) {
+            if( props.verbose == 1 )
                 cout << endl;
                 cout << endl;
-                cout << "HAK::doDoubleLoop:  WARNING: not converged within " << outer_maxiter << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.max() << endl;
+                cout << "HAK::doDoubleLoop:  WARNING: not converged within " << outer_maxiter << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.maxDiff() << endl;
             } else {
             } else {
-                if( Verbose() >= 3 )
+                if( props.verbose >= 3 )
                     cout << "HAK::doDoubleLoop:  ";
                 cout << "converged in " << outer_iter << " passes (" << toc() - tic << " clocks)." << endl;
             }
         }
 
                     cout << "HAK::doDoubleLoop:  ";
                 cout << "converged in " << outer_iter << " passes (" << toc() - tic << " clocks)." << endl;
             }
         }
 
-    return diffs.max();
+    return diffs.maxDiff();
 }
 
 
 double HAK::run() {
 }
 
 
 double HAK::run() {
-    if( DoubleLoop() )
+    if( props.doubleloop )
         return doDoubleLoop();
     else
         return doGBP();
         return doDoubleLoop();
     else
         return doGBP();
@@ -414,12 +420,12 @@ vector<Factor> HAK::beliefs() const {
 }
 
 
 }
 
 
-Complex HAK::logZ() const {
-    Complex sum = 0.0;
+Real HAK::logZ() const {
+    Real sum = 0.0;
     for( size_t beta = 0; beta < nrIRs(); beta++ )
     for( size_t beta = 0; beta < nrIRs(); beta++ )
-        sum += Complex(IR(beta).c()) * Qb(beta).entropy();
+        sum += IR(beta).c() * Qb(beta).entropy();
     for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
     for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
-        sum += Complex(OR(alpha).c()) * Qa(alpha).entropy();
+        sum += OR(alpha).c() * Qa(alpha).entropy();
         sum += (OR(alpha).log0() * Qa(alpha)).totalSum();
     }
     return sum;
         sum += (OR(alpha).log0() * Qa(alpha)).totalSum();
     }
     return sum;