Updated copyrights
[libdai.git] / src / treeep.cpp
index a12ff1b..e50be83 100644 (file)
@@ -1,6 +1,7 @@
-/*  Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]
-    Radboud University Nijmegen, The Netherlands
-    
+/*  Copyright (C) 2006-2008  Joris Mooij  [joris dot mooij at tuebingen dot mpg dot de]
+    Radboud University Nijmegen, The Netherlands /
+    Max Planck Institute for Biological Cybernetics, Germany
+
     This file is part of libDAI.
 
     libDAI is free software; you can redistribute it and/or modify
@@ -37,26 +38,41 @@ using namespace std;
 const char *TreeEP::Name = "TREEEP";
 
 
-bool TreeEP::checkProperties() {
-    if( !HasProperty("type") )
-        return false;
-    if( !HasProperty("tol") )
-        return false;
-    if (!HasProperty("maxiter") )
-        return false;
-    if (!HasProperty("verbose") )
-        return false;
+void TreeEP::setProperties( const PropertySet &opts ) {
+    assert( opts.hasKey("tol") );
+    assert( opts.hasKey("maxiter") );
+    assert( opts.hasKey("verbose") );
+    assert( opts.hasKey("type") );
     
-    ConvertPropertyTo<TypeType>("type");
-    ConvertPropertyTo<double>("tol");
-    ConvertPropertyTo<size_t>("maxiter");
-    ConvertPropertyTo<size_t>("verbose");
+    props.tol = opts.getStringAs<double>("tol");
+    props.maxiter = opts.getStringAs<size_t>("maxiter");
+    props.verbose = opts.getStringAs<size_t>("verbose");
+    props.type = opts.getStringAs<Properties::TypeType>("type");
+}
+
 
-    return true;
+PropertySet TreeEP::getProperties() const {
+    PropertySet opts;
+    opts.Set( "tol", props.tol );
+    opts.Set( "maxiter", props.maxiter );
+    opts.Set( "verbose", props.verbose );
+    opts.Set( "type", props.type );
+    return opts;
 }
 
 
-TreeEPSubTree::TreeEPSubTree( const DEdgeVec &subRTree, const DEdgeVec &jt_RTree, const vector<Factor> &jt_Qa, const vector<Factor> &jt_Qb, const Factor *I ) : _Qa(), _Qb(), _RTree(), _a(), _b(), _I(I), _ns(), _nsrem(), _logZ(0.0) {
+string TreeEP::printProperties() const {
+    stringstream s( stringstream::out );
+    s << "[";
+    s << "tol=" << props.tol << ",";
+    s << "maxiter=" << props.maxiter << ",";
+    s << "verbose=" << props.verbose << ",";
+    s << "type=" << props.type << "]";
+    return s.str();
+}
+
+
+TreeEP::TreeEPSubTree::TreeEPSubTree( const DEdgeVec &subRTree, const DEdgeVec &jt_RTree, const std::vector<Factor> &jt_Qa, const std::vector<Factor> &jt_Qb, const Factor *I ) : _Qa(), _Qb(), _RTree(), _a(), _b(), _I(I), _ns(), _nsrem(), _logZ(0.0) {
     _ns = _I->vars();
 
     // Make _Qa, _Qb, _a and _b corresponding to the subtree
@@ -91,10 +107,10 @@ TreeEPSubTree::TreeEPSubTree( const DEdgeVec &subRTree, const DEdgeVec &jt_RTree
 
     // Find remaining variables (which are not in the new root)
     _nsrem = _ns / _Qa[0].vars();
-};
+}
 
 
-void TreeEPSubTree::init() { 
+void TreeEP::TreeEPSubTree::init() { 
     for( size_t alpha = 0; alpha < _Qa.size(); alpha++ )
         _Qa[alpha].fill( 1.0 );
     for( size_t beta = 0; beta < _Qb.size(); beta++ )
@@ -102,7 +118,7 @@ void TreeEPSubTree::init() {
 }
 
 
-void TreeEPSubTree::InvertAndMultiply( const vector<Factor> &Qa, const vector<Factor> &Qb ) {
+void TreeEP::TreeEPSubTree::InvertAndMultiply( const std::vector<Factor> &Qa, const std::vector<Factor> &Qb ) {
     for( size_t alpha = 0; alpha < _Qa.size(); alpha++ )
         _Qa[alpha] = Qa[_a[alpha]].divided_by( _Qa[alpha] );
 
@@ -111,9 +127,7 @@ void TreeEPSubTree::InvertAndMultiply( const vector<Factor> &Qa, const vector<Fa
 }
 
 
-void TreeEPSubTree::HUGIN_with_I( vector<Factor> &Qa, vector<Factor> &Qb ) {
-    multind mi( _nsrem );
-
+void TreeEP::TreeEPSubTree::HUGIN_with_I( std::vector<Factor> &Qa, std::vector<Factor> &Qb ) {
     // Backup _Qa and _Qb
     vector<Factor> _Qa_old(_Qa);
     vector<Factor> _Qb_old(_Qb);
@@ -125,30 +139,27 @@ void TreeEPSubTree::HUGIN_with_I( vector<Factor> &Qa, vector<Factor> &Qb ) {
         Qb[_b[beta]].fill( 0.0 );
     
     // For all states of _nsrem
-    for( size_t j = 0; j < mi.max(); j++ ) {
-        vector<size_t> vi = mi.vi( j );
-        
+    for( State s(_nsrem); s.valid(); s++ ) {
         // Multiply root with slice of I
-        _Qa[0] *= _I->slice( _nsrem, j );
+        _Qa[0] *= _I->slice( _nsrem, s );
 
         // CollectEvidence
         for( size_t i = _RTree.size(); (i--) != 0; ) {
             // clamp variables in nsrem
-            size_t k = 0;
-            for( VarSet::const_iterator n = _nsrem.begin(); n != _nsrem.end(); n++, k++ )
+            for( VarSet::const_iterator n = _nsrem.begin(); n != _nsrem.end(); n++ )
                 if( _Qa[_RTree[i].n2].vars() >> *n ) {
                     Factor delta( *n, 0.0 );
-                    delta[vi[k]] = 1.0;
+                    delta[s(*n)] = 1.0;
                     _Qa[_RTree[i].n2] *= delta;
                 }
-            Factor new_Qb = _Qa[_RTree[i].n2].part_sum( _Qb[i].vars() );
+            Factor new_Qb = _Qa[_RTree[i].n2].partSum( _Qb[i].vars() );
             _Qa[_RTree[i].n1] *= new_Qb.divided_by( _Qb[i] ); 
             _Qb[i] = new_Qb;
         }
 
         // DistributeEvidence
         for( size_t i = 0; i < _RTree.size(); i++ ) {
-            Factor new_Qb = _Qa[_RTree[i].n1].part_sum( _Qb[i].vars() );
+            Factor new_Qb = _Qa[_RTree[i].n1].partSum( _Qb[i].vars() );
             _Qa[_RTree[i].n2] *= new_Qb.divided_by( _Qb[i] ); 
             _Qb[i] = new_Qb;
         }
@@ -168,16 +179,16 @@ void TreeEPSubTree::HUGIN_with_I( vector<Factor> &Qa, vector<Factor> &Qb ) {
     _logZ = 0.0;
     for( size_t alpha = 0; alpha < _Qa.size(); alpha++ ) {
         _logZ += log(Qa[_a[alpha]].totalSum());
-        Qa[_a[alpha]].normalize( Prob::NORMPROB );
+        Qa[_a[alpha]].normalize();
     }
     for( size_t beta = 0; beta < _Qb.size(); beta++ ) {
         _logZ -= log(Qb[_b[beta]].totalSum());
-        Qb[_b[beta]].normalize( Prob::NORMPROB );
+        Qb[_b[beta]].normalize();
     }
 }
 
 
-double TreeEPSubTree::logZ( const vector<Factor> &Qa, const vector<Factor> &Qb ) const {
+double TreeEP::TreeEPSubTree::logZ( const std::vector<Factor> &Qa, const std::vector<Factor> &Qb ) const {
     double sum = 0.0;
     for( size_t alpha = 0; alpha < _Qa.size(); alpha++ )
         sum += (Qa[_a[alpha]] * _Qa[alpha].log0()).totalSum();
@@ -187,73 +198,57 @@ double TreeEPSubTree::logZ( const vector<Factor> &Qa, const vector<Factor> &Qb )
 }
 
 
-TreeEP::TreeEP( const FactorGraph &fg, const Properties &opts ) : JTree(fg, opts("updates",string("HUGIN")), false) {
-    assert( checkProperties() );
+TreeEP::TreeEP( const FactorGraph &fg, const PropertySet &opts ) : JTree(fg, opts("updates",string("HUGIN")), false), _maxdiff(0.0), _iters(0), props(), _Q() {
+    setProperties( opts );
 
-    assert( fg.G.isConnected() );
+    assert( fg.isConnected() );
 
     if( opts.hasKey("tree") ) {
         ConstructRG( opts.GetAs<DEdgeVec>("tree") );
     } else {
-        if( Type() == TypeType::ORG ) {
-            // construct weighted graph with as weights a crude estimate of the
+        if( props.type == Properties::TypeType::ORG || props.type == Properties::TypeType::ALT ) {
+            // ORG: construct weighted graph with as weights a crude estimate of the
             // mutual information between the nodes
-            WeightedGraph<double> wg;
-            for( size_t i = 0; i < nrVars(); ++i ) {
-                Var v_i = var(i);
-                VarSet di = delta(i);
-                for( VarSet::const_iterator j = di.begin(); j != di.end(); j++ )
-                    if( v_i < *j ) {
-                        Factor piet;
-                        for( size_t I = 0; I < nrFactors(); I++ ) {
-                            VarSet Ivars = factor(I).vars();
-                            if( (Ivars == v_i) || (Ivars == *j) )
-                                piet *= factor(I);
-                            else if( Ivars >> (v_i | *j) )
-                                piet *= factor(I).marginal( v_i | *j );
-                        }
-                        if( piet.vars() >> (v_i | *j) ) {
-                            piet = piet.marginal( v_i | *j );
-                            Factor pietf = piet.marginal(v_i) * piet.marginal(*j);
-                            wg[UEdge(i,findVar(*j))] = real( KL_dist( piet, pietf ) );
-                        } else
-                            wg[UEdge(i,findVar(*j))] = 0;
-                    }
-            }
-
-            // find maximal spanning tree
-            ConstructRG( MaxSpanningTreePrim( wg ) );
-
-//            cout << "Constructing maximum spanning tree..." << endl;
-//            DEdgeVec MST = MaxSpanningTreePrim( wg );
-//            cout << "Maximum spanning tree:" << endl;
-//            for( DEdgeVec::const_iterator e = MST.begin(); e != MST.end(); e++ )
-//                cout << *e << endl; 
-//            ConstructRG( MST );
-        } else if( Type() == TypeType::ALT ) {
-            // construct weighted graph with as weights an upper bound on the
+            // ALT: construct weighted graph with as weights an upper bound on the
             // effective interaction strength between pairs of nodes
+            
             WeightedGraph<double> wg;
             for( size_t i = 0; i < nrVars(); ++i ) {
                 Var v_i = var(i);
                 VarSet di = delta(i);
                 for( VarSet::const_iterator j = di.begin(); j != di.end(); j++ )
                     if( v_i < *j ) {
+                        VarSet ij(v_i,*j);
                         Factor piet;
                         for( size_t I = 0; I < nrFactors(); I++ ) {
                             VarSet Ivars = factor(I).vars();
-                            if( Ivars >> (v_i | *j) )
-                                piet *= factor(I);
+                            if( props.type == Properties::TypeType::ORG ) {
+                                if( (Ivars == v_i) || (Ivars == *j) )
+                                    piet *= factor(I);
+                                else if( Ivars >> ij )
+                                    piet *= factor(I).marginal( ij );
+                            } else {
+                                if( Ivars >> ij )
+                                    piet *= factor(I);
+                            }
+                        }
+                        if( props.type == Properties::TypeType::ORG ) {
+                            if( piet.vars() >> ij ) {
+                                piet = piet.marginal( ij );
+                                Factor pietf = piet.marginal(v_i) * piet.marginal(*j);
+                                wg[UEdge(i,findVar(*j))] = KL_dist( piet, pietf );
+                            } else
+                                wg[UEdge(i,findVar(*j))] = 0;
+                        } else {
+                            wg[UEdge(i,findVar(*j))] = piet.strength(v_i, *j);
                         }
-                        wg[UEdge(i,findVar(*j))] = piet.strength(v_i, *j);
                     }
             }
 
             // find maximal spanning tree
-            ConstructRG( MaxSpanningTreePrim( wg ) );
-        } else {
-            assert( 0 == 1 );
-        }
+            ConstructRG( MaxSpanningTreePrims( wg ) );
+        } else
+            DAI_THROW(INTERNAL_ERROR);
     }
 }
 
@@ -270,11 +265,12 @@ void TreeEP::ConstructRG( const DEdgeVec &tree ) {
     for( size_t i = 0; i < Cliques.size(); i++ )
         for( size_t j = i+1; j < Cliques.size(); j++ ) {
             size_t w = (Cliques[i] & Cliques[j]).size();
-            JuncGraph[UEdge(i,j)] = w;
+            if( w )
+                JuncGraph[UEdge(i,j)] = w;
         }
     
     // Construct maximal spanning tree using Prim's algorithm
-    _RTree = MaxSpanningTreePrim( JuncGraph );
+    RTree = MaxSpanningTreePrims( JuncGraph );
 
     // Construct corresponding region graph
 
@@ -301,33 +297,32 @@ void TreeEP::ConstructRG( const DEdgeVec &tree ) {
     RecomputeORs();
 
     // Create inner regions and edges
-    IRs.reserve( _RTree.size() );
-    typedef pair<size_t,size_t> Edge;
+    IRs.reserve( RTree.size() );
     vector<Edge> edges;
-    edges.reserve( 2 * _RTree.size() );
-    for( size_t i = 0; i < _RTree.size(); i++ ) {
-        edges.push_back( Edge( _RTree[i].n1, IRs.size() ) );
-        edges.push_back( Edge( _RTree[i].n2, IRs.size() ) );
+    edges.reserve( 2 * RTree.size() );
+    for( size_t i = 0; i < RTree.size(); i++ ) {
+        edges.push_back( Edge( RTree[i].n1, IRs.size() ) );
+        edges.push_back( Edge( RTree[i].n2, IRs.size() ) );
         // inner clusters have counting number -1
-        IRs.push_back( Region( Cliques[_RTree[i].n1] & Cliques[_RTree[i].n2], -1.0 ) );
+        IRs.push_back( Region( Cliques[RTree[i].n1] & Cliques[RTree[i].n2], -1.0 ) );
     }
 
     // create bipartite graph
-    G.create( nrORs(), nrIRs(), edges.begin(), edges.end() );
+    G.construct( nrORs(), nrIRs(), edges.begin(), edges.end() );
 
     // Check counting numbers
     Check_Counting_Numbers();
     
     // Create messages and beliefs
-    _Qa.clear();
-    _Qa.reserve( nrORs() );
+    Qa.clear();
+    Qa.reserve( nrORs() );
     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
-        _Qa.push_back( OR(alpha) );
+        Qa.push_back( OR(alpha) );
 
-    _Qb.clear();
-    _Qb.reserve( nrIRs() );
+    Qb.clear();
+    Qb.reserve( nrIRs() );
     for( size_t beta = 0; beta < nrIRs(); beta++ ) 
-        _Qb.push_back( Factor( IR(beta), 1.0 ) );
+        Qb.push_back( Factor( IR(beta), 1.0 ) );
 
     // DIFF with JTree::GenerateJT:  no messages
     
@@ -344,27 +339,7 @@ void TreeEP::ConstructRG( const DEdgeVec &tree ) {
             //subTree.resize( subTreeSize );  // FIXME
 //          cout << "subtree " << I << " has size " << subTreeSize << endl;
 
-/*
-            char fn[30];
-            sprintf( fn, "/tmp/subtree_%d.dot", I );
-            std::ofstream dots(fn);
-            dots << "graph G {" << endl;
-            dots << "graph[size=\"9,9\"];" << endl;
-            dots << "node[shape=circle,width=0.4,fixedsize=true];" << endl;
-            for( size_t i = 0; i < nrVars(); i++ )
-                dots << "\tx" << var(i).label() << ((factor(I).vars() >> var(i)) ? "[color=blue];" : ";") << endl;
-            dots << "node[shape=box,style=filled,color=lightgrey,width=0.3,height=0.3,fixedsize=true];" << endl;
-            for( size_t J = 0; J < nrFactors(); J++ )
-                dots << "\tp" << J << ";" << endl;
-            for( size_t iI = 0; iI < FactorGraph::nr_edges(); iI++ )
-                dots << "\tx" << var(FactorGraph::edge(iI).first).label() << " -- p" << FactorGraph::edge(iI).second << ";" << endl;
-            for( size_t a = 0; a < tree.size(); a++ )
-                dots << "\tx" << var(tree[a].n1).label() << " -- x" << var(tree[a].n2).label() << " [color=red];" << endl;
-            dots << "}" << endl;
-            dots.close();
-*/
-            
-            TreeEPSubTree QI( subTree, _RTree, _Qa, _Qb, &factor(I) );
+            TreeEPSubTree QI( subTree, RTree, Qa, Qb, &factor(I) );
             _Q[I] = QI;
         }
     // Previous root of first off-tree factor should be the root of the last off-tree factor
@@ -376,27 +351,23 @@ void TreeEP::ConstructRG( const DEdgeVec &tree ) {
             //subTree.resize( subTreeSize ); // FIXME
 //          cout << "subtree " << I << " has size " << subTreeSize << endl;
 
-            TreeEPSubTree QI( subTree, _RTree, _Qa, _Qb, &factor(I) );
+            TreeEPSubTree QI( subTree, RTree, Qa, Qb, &factor(I) );
             _Q[I] = QI;
             break;
         }
 
-    if( Verbose() >= 3 ) {
+    if( props.verbose >= 3 ) {
         cout << "Resulting regiongraph: " << *this << endl;
     }
 }
 
 
 string TreeEP::identify() const { 
-    stringstream result (stringstream::out);
-    result << Name << GetProperties();
-    return result.str();
+    return string(Name) + printProperties();
 }
 
 
 void TreeEP::init() {
-    assert( checkProperties() );
-
     runHUGIN();
 
     // Init factor approximations
@@ -407,12 +378,12 @@ void TreeEP::init() {
 
 
 double TreeEP::run() {
-    if( Verbose() >= 1 )
+    if( props.verbose >= 1 )
         cout << "Starting " << identify() << "...";
-    if( Verbose() >= 3)
+    if( props.verbose >= 3)
         cout << endl;
 
-    clock_t tic = toc();
+    double tic = toc();
     Diffs diffs(nrVars(), 1.0);
 
     vector<Factor> old_beliefs;
@@ -420,16 +391,14 @@ double TreeEP::run() {
     for( size_t i = 0; i < nrVars(); i++ )
         old_beliefs.push_back(belief(var(i)));
 
-    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( _iters=0; _iters < props.maxiter && diffs.maxDiff() > props.tol; _iters++ ) {
         for( size_t I = 0; I < nrFactors(); I++ )
             if( offtree(I) ) {  
-                _Q[I].InvertAndMultiply( _Qa, _Qb );
-                _Q[I].HUGIN_with_I( _Qa, _Qb );
-                _Q[I].InvertAndMultiply( _Qa, _Qb );
+                _Q[I].InvertAndMultiply( Qa, Qb );
+                _Q[I].HUGIN_with_I( Qa, Qb );
+                _Q[I].InvertAndMultiply( Qa, Qb );
             }
 
         // calculate new beliefs and compare with old ones
@@ -439,45 +408,46 @@ double TreeEP::run() {
             old_beliefs[i] = nb;
         }
 
-        if( Verbose() >= 3 )
-            cout << "TreeEP::run:  maxdiff " << diffs.max() << " after " << iter+1 << " passes" << endl;
+        if( props.verbose >= 3 )
+            cout << Name << "::run:  maxdiff " << diffs.maxDiff() << " after " << _iters+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 << "TreeEP::run:  WARNING: not converged within " << MaxIter() << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.max() << endl;
+            cout << Name << "::run:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << diffs.maxDiff() << endl;
         } else {
-            if( Verbose() >= 3 )
-                cout << "TreeEP::run:  ";
-            cout << "converged in " << iter << " passes (" << toc() - tic << " clocks)." << endl;
+            if( props.verbose >= 3 )
+                cout << Name << "::run:  ";
+            cout << "converged in " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
         }
     }
 
-    return diffs.max();
+    return diffs.maxDiff();
 }
 
 
-Complex TreeEP::logZ() const {
+Real TreeEP::logZ() const {
     double sum = 0.0;
 
     // entropy of the tree
     for( size_t beta = 0; beta < nrIRs(); beta++ )
-        sum -= real(_Qb[beta].entropy());
+        sum -= Qb[beta].entropy();
     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
-        sum += real(_Qa[alpha].entropy());
+        sum += Qa[alpha].entropy();
 
     // energy of the on-tree factors
     for( size_t alpha = 0; alpha < nrORs(); alpha++ )
-        sum += (OR(alpha).log0() * _Qa[alpha]).totalSum();
+        sum += (OR(alpha).log0() * Qa[alpha]).totalSum();
 
     // energy of the off-tree factors
     for( size_t I = 0; I < nrFactors(); I++ )
         if( offtree(I) )
-            sum += (_Q.find(I))->second.logZ( _Qa, _Qb );
+            sum += (_Q.find(I))->second.logZ( Qa, Qb );
     
     return sum;
 }