Replaced sub_nb class in mr.h by boost::dynamic_bitset
[libdai.git] / src / jtree.cpp
index b82b0e5..0f28d5b 100644 (file)
@@ -32,53 +32,64 @@ using namespace std;
 const char *JTree::Name = "JTREE";
 
 
-bool JTree::checkProperties() {
-    if (!HasProperty("verbose") )
-        return false;
-    if( !HasProperty("updates") )
-        return false;
+void JTree::setProperties( const PropertySet &opts ) {
+    assert( opts.hasKey("verbose") );
+    assert( opts.hasKey("updates") );
     
-    ConvertPropertyTo<size_t>("verbose");
-    ConvertPropertyTo<UpdateType>("updates");
+    props.verbose = opts.getStringAs<size_t>("verbose");
+    props.updates = opts.getStringAs<Properties::UpdateType>("updates");
+}
+
+
+PropertySet JTree::getProperties() const {
+    PropertySet opts;
+    opts.Set( "verbose", props.verbose );
+    opts.Set( "updates", props.updates );
+    return opts;
+}
+
 
-    return true;
+string JTree::printProperties() const {
+    stringstream s( stringstream::out );
+    s << "[";
+    s << "verbose=" << props.verbose << ",";
+    s << "updates=" << props.updates << "]";
+    return s.str();
 }
 
 
-JTree::JTree( const FactorGraph &fg, const Properties &opts, bool automatic) : DAIAlgRG(fg, opts), _RTree(), _Qa(), _Qb(), _mes(), _logZ() {
-    assert( checkProperties() );
+JTree::JTree( const FactorGraph &fg, const PropertySet &opts, bool automatic ) : DAIAlgRG(fg), _RTree(), _Qa(), _Qb(), _mes(), _logZ(), props() {
+    setProperties( opts );
 
-    if( automatic ) {
-        ClusterGraph _cg;
+    if( !isConnected() ) 
+       DAI_THROW(FACTORGRAPH_NOT_CONNECTED); 
 
-        // Copy factors
+    if( automatic ) {
+        // Create ClusterGraph which contains factors as clusters
+        vector<VarSet> cl;
+        cl.reserve( fg.nrFactors() );
         for( size_t I = 0; I < nrFactors(); I++ )
-            _cg.insert( factor(I).vars() );
-        if( Verbose() >= 3 )
+            cl.push_back( factor(I).vars() );
+        ClusterGraph _cg( cl );
+
+        if( props.verbose >= 3 )
             cout << "Initial clusters: " << _cg << endl;
 
         // Retain only maximal clusters
         _cg.eraseNonMaximal();
-        if( Verbose() >= 3 )
+        if( props.verbose >= 3 )
             cout << "Maximal clusters: " << _cg << endl;
 
         vector<VarSet> ElimVec = _cg.VarElim_MinFill().eraseNonMaximal().toVector();
-        if( Verbose() >= 3 ) {
-            cout << "VarElim_MinFill result: {" << endl;
-            for( size_t i = 0; i < ElimVec.size(); i++ ) {
-                if( i != 0 )
-                    cout << ", ";
-                cout << ElimVec[i];
-            }
-            cout << "}" << endl;
-        }
+        if( props.verbose >= 3 )
+            cout << "VarElim_MinFill result: " << ElimVec << endl;
 
         GenerateJT( ElimVec );
     }
 }
 
 
-void JTree::GenerateJT( const vector<VarSet> &Cliques ) {
+void JTree::GenerateJT( const std::vector<VarSet> &Cliques ) {
     // Construct a weighted graph (each edge is weighted with the cardinality 
     // of the intersection of the nodes, where the nodes are the elements of
     // Cliques).
@@ -86,75 +97,78 @@ void JTree::GenerateJT( const vector<VarSet> &Cliques ) {
     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
 
     // Create outer regions
-    ORs().reserve( Cliques.size() );
+    ORs.reserve( Cliques.size() );
     for( size_t i = 0; i < Cliques.size(); i++ )
-        ORs().push_back( FRegion( Factor(Cliques[i], 1.0), 1.0 ) );
+        ORs.push_back( FRegion( Factor(Cliques[i], 1.0), 1.0 ) );
 
     // For each factor, find an outer region that subsumes that factor.
     // Then, multiply the outer region with that factor.
     for( size_t I = 0; I < nrFactors(); I++ ) {
         size_t alpha;
-        for( alpha = 0; alpha < nr_ORs(); alpha++ )
+        for( alpha = 0; alpha < nrORs(); alpha++ )
             if( OR(alpha).vars() >> factor(I).vars() ) {
-//              OR(alpha) *= factor(I);
-                _fac2OR[I] = alpha;
+                fac2OR.push_back( alpha );
                 break;
             }
-        assert( alpha != nr_ORs() );
+        assert( alpha != nrORs() );
     }
     RecomputeORs();
 
     // Create inner regions and edges
-    IRs().reserve( _RTree.size() );
-    Redges().reserve( 2 * _RTree.size() );
+    IRs.reserve( _RTree.size() );
+    vector<Edge> edges;
+    edges.reserve( 2 * _RTree.size() );
     for( size_t i = 0; i < _RTree.size(); i++ ) {
-        Redges().push_back( R_edge_t( _RTree[i].n1, IRs().size() ) );
-        Redges().push_back( R_edge_t( _RTree[i].n2, IRs().size() ) );
+        edges.push_back( Edge( _RTree[i].n1, nrIRs() ) );
+        edges.push_back( Edge( _RTree[i].n2, nrIRs() ) );
         // 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 ) );
     }
 
-    // Regenerate BipartiteGraph internals
-    Regenerate();
+    // create bipartite graph
+    G.construct( nrORs(), nrIRs(), edges.begin(), edges.end() );
 
     // Create messages and beliefs
     _Qa.clear();
-    _Qa.reserve( nr_ORs() );
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+    _Qa.reserve( nrORs() );
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ )
         _Qa.push_back( OR(alpha) );
 
     _Qb.clear();
-    _Qb.reserve( nr_IRs() );
-    for( size_t beta = 0; beta < nr_IRs(); beta++ ) 
+    _Qb.reserve( nrIRs() );
+    for( size_t beta = 0; beta < nrIRs(); beta++ ) 
         _Qb.push_back( Factor( IR(beta), 1.0 ) );
 
     _mes.clear();
-    _mes.reserve( nr_Redges() );
-    for( size_t e = 0; e < nr_Redges(); e++ )
-        _mes.push_back( Factor( IR(Redge(e).second), 1.0 ) );
+    _mes.reserve( nrORs() );
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
+        _mes.push_back( vector<Factor>() );
+        _mes[alpha].reserve( nbOR(alpha).size() );
+        foreach( const Neighbor &beta, nbOR(alpha) )
+            _mes[alpha].push_back( Factor( IR(beta), 1.0 ) );
+    }
 
     // Check counting numbers
     Check_Counting_Numbers();
 
-    if( Verbose() >= 3 ) {
+    if( props.verbose >= 3 ) {
         cout << "Resulting regiongraph: " << *this << endl;
     }
 }
 
 
 string JTree::identify() const {
-    stringstream result (stringstream::out);
-    result << Name << GetProperties();
-    return result.str();
+    return string(Name) + printProperties();
 }
 
 
@@ -178,9 +192,9 @@ Factor JTree::belief( const VarSet &ns ) const {
 
 vector<Factor> JTree::beliefs() const {
     vector<Factor> result;
-    for( size_t beta = 0; beta < nr_IRs(); beta++ )
+    for( size_t beta = 0; beta < nrIRs(); beta++ )
         result.push_back( _Qb[beta] );
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ )
         result.push_back( _Qa[alpha] );
     return result;
 }
@@ -193,10 +207,10 @@ Factor JTree::belief( const Var &n ) const {
 
 // Needs no init
 void JTree::runHUGIN() {
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ )
         _Qa[alpha] = OR(alpha);
 
-    for( size_t beta = 0; beta < nr_IRs(); beta++ )
+    for( size_t beta = 0; beta < nrIRs(); beta++ )
         _Qb[beta].fill( 1.0 );
 
     // CollectEvidence
@@ -204,15 +218,15 @@ void JTree::runHUGIN() {
     for( size_t i = _RTree.size(); (i--) != 0; ) {
 //      Make outer region _RTree[i].n1 consistent with outer region _RTree[i].n2
 //      IR(i) = seperator OR(_RTree[i].n1) && OR(_RTree[i].n2)
-        Factor new_Qb = _Qa[_RTree[i].n2].part_sum( IR( i ) );
-        _logZ += log(new_Qb.normalize( Prob::NORMPROB ));
+        Factor new_Qb = _Qa[_RTree[i].n2].partSum( IR( i ) );
+        _logZ += log(new_Qb.normalize());
         _Qa[_RTree[i].n1] *= new_Qb.divided_by( _Qb[i] ); 
         _Qb[i] = new_Qb;
     }
     if( _RTree.empty() )
-        _logZ += log(_Qa[0].normalize( Prob::NORMPROB ) );
+        _logZ += log(_Qa[0].normalize() );
     else
-        _logZ += log(_Qa[_RTree[0].n1].normalize( Prob::NORMPROB ));
+        _logZ += log(_Qa[_RTree[0].n1].normalize());
 
     // DistributeEvidence
     for( size_t i = 0; i < _RTree.size(); i++ ) {
@@ -224,8 +238,8 @@ void JTree::runHUGIN() {
     }
 
     // Normalize
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
-        _Qa[alpha].normalize( Prob::NORMPROB );
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ )
+        _Qa[alpha].normalize();
 }
 
 
@@ -233,69 +247,71 @@ void JTree::runHUGIN() {
 void JTree::runShaferShenoy() {
     // First pass
     _logZ = 0.0;
-    for( size_t e = _RTree.size(); (e--) != 0; ) {
+    for( size_t e = nrIRs(); (e--) != 0; ) {
         // send a message from _RTree[e].n2 to _RTree[e].n1
         // or, actually, from the seperator IR(e) to _RTree[e].n1
 
-        size_t i = _RTree[e].n2;
-        size_t j = _RTree[e].n1;
+        size_t i = nbIR(e)[1].node; // = _RTree[e].n2
+        size_t j = nbIR(e)[0].node; // = _RTree[e].n1
+        size_t _e = nbIR(e)[0].dual;
         
         Factor piet = OR(i);
-        for( R_nb_cit k = nbOR(i).begin(); k != nbOR(i).end(); k++ )
-            if( *k != e )
-                piet *= message( i, *k );
-        message( j, e ) = piet.part_sum( IR(e) );
-        _logZ += log( message(j,e).normalize( Prob::NORMPROB ) );
+        foreach( const Neighbor &k, nbOR(i) )
+            if( k != e ) 
+                piet *= message( i, k.iter );
+        message( j, _e ) = piet.partSum( IR(e) );
+        _logZ += log( message(j,_e).normalize() );
     }
 
     // Second pass
-    for( size_t e = 0; e < _RTree.size(); e++ ) {
-        size_t i = _RTree[e].n1;
-        size_t j = _RTree[e].n2;
+    for( size_t e = 0; e < nrIRs(); e++ ) {
+        size_t i = nbIR(e)[0].node; // = _RTree[e].n1
+        size_t j = nbIR(e)[1].node; // = _RTree[e].n2
+        size_t _e = nbIR(e)[1].dual;
         
         Factor piet = OR(i);
-        for( R_nb_cit k = nbOR(i).begin(); k != nbOR(i).end(); k++ )
-            if( *k != e )
-                piet *= message( i, *k );
-        message( j, e ) = piet.marginal( IR(e) );
+        foreach( const Neighbor &k, nbOR(i) )
+            if( k != e )
+                piet *= message( i, k.iter );
+        message( j, _e ) = piet.marginal( IR(e) );
     }
 
     // Calculate beliefs
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
         Factor piet = OR(alpha);
-        for( R_nb_cit k = nbOR(alpha).begin(); k != nbOR(alpha).end(); k++ )
-            piet *= message( alpha, *k );
-        if( _RTree.empty() ) {
-            _logZ += log( piet.normalize( Prob::NORMPROB ) );
+        foreach( const Neighbor &k, nbOR(alpha) )
+            piet *= message( alpha, k.iter );
+        if( nrIRs() == 0 ) {
+            _logZ += log( piet.normalize() );
             _Qa[alpha] = piet;
-        } else if( alpha == _RTree[0].n1 ) {
-            _logZ += log( piet.normalize( Prob::NORMPROB ) );
+        } else if( alpha == nbIR(0)[0].node /*_RTree[0].n1*/ ) {
+            _logZ += log( piet.normalize() );
             _Qa[alpha] = piet;
         } else
-            _Qa[alpha] = piet.normalized( Prob::NORMPROB );
+            _Qa[alpha] = piet.normalized();
     }
 
     // Only for logZ (and for belief)...
-    for( size_t beta = 0; beta < nr_IRs(); beta++ ) 
-        _Qb[beta] = _Qa[nbIR(beta)[0]].marginal( IR(beta) );
+    for( size_t beta = 0; beta < nrIRs(); beta++ ) 
+        _Qb[beta] = _Qa[nbIR(beta)[0].node].marginal( IR(beta) );
 }
 
 
 double JTree::run() {
-    if( Updates() == UpdateType::HUGIN )
+    if( props.updates == Properties::UpdateType::HUGIN )
         runHUGIN();
-    else if( Updates() == UpdateType::SHSH )
+    else if( props.updates == Properties::UpdateType::SHSH )
         runShaferShenoy();
     return 0.0;
 }
 
 
-Complex JTree::logZ() const {
-    Complex sum = 0.0;
-    for( size_t beta = 0; beta < nr_IRs(); beta++ )
-        sum += Complex(IR(beta).c()) * _Qb[beta].entropy();
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
-        sum += Complex(OR(alpha).c()) * _Qa[alpha].entropy();
+Real JTree::logZ() const {
+    Real sum = 0.0;
+    for( size_t beta = 0; beta < nrIRs(); beta++ )
+        sum += IR(beta).c() * _Qb[beta].entropy();
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
+        sum += OR(alpha).c() * _Qa[alpha].entropy();
         sum += (OR(alpha).log0() * _Qa[alpha]).totalSum();
     }
     return sum;
@@ -306,37 +322,29 @@ Complex JTree::logZ() const {
 size_t JTree::findEfficientTree( const VarSet& ns, DEdgeVec &Tree, size_t PreviousRoot ) const {
     // find new root clique (the one with maximal statespace overlap with ns)
     size_t maxval = 0, maxalpha = 0;
-    for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
-        size_t val = (ns & OR(alpha).vars()).stateSpace();
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
+        size_t val = (ns & OR(alpha).vars()).states();
         if( val > maxval ) {
             maxval = val;
             maxalpha = alpha;
         }
     }
 
-//  for( size_t e = 0; e < _RTree.size(); e++ )
-//      cout << OR(_RTree[e].n1).vars() << "->" << OR(_RTree[e].n2).vars() << ",  ";
-//  cout << endl;
     // grow new tree
     Graph oldTree;
     for( DEdgeVec::const_iterator e = _RTree.begin(); e != _RTree.end(); e++ )
         oldTree.insert( UEdge(e->n1, e->n2) );
     DEdgeVec newTree = GrowRootedTree( oldTree, maxalpha );
-//  cout << ns << ": ";
-//  for( size_t e = 0; e < newTree.size(); e++ )
-//      cout << OR(newTree[e].n1).vars() << "->" << OR(newTree[e].n2).vars() << ",  ";
-//  cout << endl;
     
     // identify subtree that contains variables of ns which are not in the new root
     VarSet nsrem = ns / OR(maxalpha).vars();
-//  cout << "nsrem:" << nsrem << endl;
     set<DEdge> subTree;
     // for each variable in ns that is not in the root clique
     for( VarSet::const_iterator n = nsrem.begin(); n != nsrem.end(); n++ ) {
         // find first occurence of *n in the tree, which is closest to the root
         size_t e = 0;
         for( ; e != newTree.size(); e++ ) {
-            if( OR(newTree[e].n2).vars() && *n )
+            if( OR(newTree[e].n2).vars().contains( *n ) )
                 break;
         }
         assert( e != newTree.size() );
@@ -368,10 +376,6 @@ size_t JTree::findEfficientTree( const VarSet& ns, DEdgeVec &Tree, size_t Previo
                 pos = newTree[e-1].n1;
             }
     }
-//  cout << "subTree: " << endl;
-//  for( set<DEdge>::const_iterator sTi = subTree.begin(); sTi != subTree.end(); sTi++ )
-//      cout << OR(sTi->n1).vars() << "->" << OR(sTi->n2).vars() << ",  ";
-//  cout << endl;
 
     // Resulting Tree is a reordered copy of newTree
     // First add edges in subTree to Tree
@@ -379,9 +383,7 @@ size_t JTree::findEfficientTree( const VarSet& ns, DEdgeVec &Tree, size_t Previo
     for( DEdgeVec::const_iterator e = newTree.begin(); e != newTree.end(); e++ )
         if( subTree.count( *e ) ) {
             Tree.push_back( *e );
-//          cout << OR(e->n1).vars() << "->" << OR(e->n2).vars() << ",  ";
         }
-//  cout << endl;
     // Then add edges pointing away from nsrem
     // FIXME
 /*  for( DEdgeVec::const_iterator e = newTree.begin(); e != newTree.end(); e++ )
@@ -390,7 +392,6 @@ size_t JTree::findEfficientTree( const VarSet& ns, DEdgeVec &Tree, size_t Previo
                 if( e->n1 == sTi->n1 || e->n1 == sTi->n2 ||
                     e->n2 == sTi->n1 || e->n2 == sTi->n2 ) {
                     Tree.push_back( *e );
-//                  cout << OR(e->n1).vars() << "->" << OR(e->n2).vars() << ",  ";
                 }
             }*/
     // FIXME
@@ -404,10 +405,8 @@ size_t JTree::findEfficientTree( const VarSet& ns, DEdgeVec &Tree, size_t Previo
                 }
             if( found ) {
                 Tree.push_back( *e );
-                cout << OR(e->n1).vars() << "->" << OR(e->n2).vars() << ",  ";
             }
-        }
-    cout << endl;*/
+        }*/
     size_t subTreeSize = Tree.size();
     // Then add remaining edges
     for( DEdgeVec::const_iterator e = newTree.begin(); e != newTree.end(); e++ )
@@ -443,8 +442,6 @@ Factor JTree::calcMarginal( const VarSet& ns ) {
             VarSet nsrem = ns / OR(T.front().n1).vars();
             Factor Pns (ns, 0.0);
             
-            multind mi( nsrem );
-
             // Save _Qa and _Qb on the subtree
             map<size_t,Factor> _Qa_old;
             map<size_t,Factor> _Qb_old;
@@ -453,10 +450,10 @@ Factor JTree::calcMarginal( const VarSet& ns ) {
                 size_t alpha1 = T[i].n1;
                 size_t alpha2 = T[i].n2;
                 size_t beta;
-                for( beta = 0; beta < nr_IRs(); beta++ )
+                for( beta = 0; beta < nrIRs(); beta++ )
                     if( UEdge( _RTree[beta].n1, _RTree[beta].n2 ) == UEdge( alpha1, alpha2 ) )
                         break;
-                assert( beta != nr_IRs() );
+                assert( beta != nrIRs() );
                 b[i] = beta;
 
                 if( !_Qa_old.count( alpha1 ) )
@@ -468,33 +465,30 @@ Factor JTree::calcMarginal( const VarSet& ns ) {
             }
                 
             // 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++ ) {
                 // CollectEvidence
                 double logZ = 0.0;
                 for( size_t i = Tsize; (i--) != 0; ) {
-            //      Make outer region T[i].n1 consistent with outer region T[i].n2
-            //      IR(i) = seperator OR(T[i].n1) && OR(T[i].n2)
+                // Make outer region T[i].n1 consistent with outer region T[i].n2
+                // IR(i) = seperator OR(T[i].n1) && OR(T[i].n2)
 
-                    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[T[i].n2].vars() >> *n ) {
                             Factor piet( *n, 0.0 );
-                            piet[vi[k]] = 1.0;
+                            piet[s(*n)] = 1.0;
                             _Qa[T[i].n2] *= piet; 
                         }
 
-                    Factor new_Qb = _Qa[T[i].n2].part_sum( IR( b[i] ) );
-                    logZ += log(new_Qb.normalize( Prob::NORMPROB ));
+                    Factor new_Qb = _Qa[T[i].n2].partSum( IR( b[i] ) );
+                    logZ += log(new_Qb.normalize());
                     _Qa[T[i].n1] *= new_Qb.divided_by( _Qb[b[i]] ); 
                     _Qb[b[i]] = new_Qb;
                 }
-                logZ += log(_Qa[T[0].n1].normalize( Prob::NORMPROB ));
+                logZ += log(_Qa[T[0].n1].normalize());
 
                 Factor piet( nsrem, 0.0 );
-                piet[j] = exp(logZ);
-                Pns += piet * _Qa[T[0].n1].part_sum( ns / nsrem );      // OPTIMIZE ME
+                piet[s] = exp(logZ);
+                Pns += piet * _Qa[T[0].n1].partSum( ns / nsrem );      // OPTIMIZE ME
 
                 // Restore clamped beliefs
                 for( map<size_t,Factor>::const_iterator alpha = _Qa_old.begin(); alpha != _Qa_old.end(); alpha++ )
@@ -503,10 +497,39 @@ Factor JTree::calcMarginal( const VarSet& ns ) {
                     _Qb[beta->first] = beta->second;
             }
 
-            return( Pns.normalized(Prob::NORMPROB) );
+            return( Pns.normalized() );
         }
     }
 }
 
 
+// first return value is treewidth
+// second return value is number of states in largest clique
+pair<size_t,size_t> treewidth( const FactorGraph & fg ) {
+    ClusterGraph _cg;
+
+    // Copy factors
+    for( size_t I = 0; I < fg.nrFactors(); I++ )
+        _cg.insert( fg.factor(I).vars() );
+
+    // Retain only maximal clusters
+    _cg.eraseNonMaximal();
+
+    // Obtain elimination sequence
+    vector<VarSet> ElimVec = _cg.VarElim_MinFill().eraseNonMaximal().toVector();
+
+    // Calculate treewidth
+    size_t treewidth = 0;
+    size_t nrstates = 0;
+    for( size_t i = 0; i < ElimVec.size(); i++ ) {
+        if( ElimVec[i].size() > treewidth )
+            treewidth = ElimVec[i].size();
+        if( ElimVec[i].states() > nrstates )
+            nrstates = ElimVec[i].states();
+    }
+
+    return pair<size_t,size_t>(treewidth, nrstates);
+}
+
+
 } // end of namespace dai