New git HEAD version
[libdai.git] / src / factorgraph.cpp
index fcbf626..3cff7eb 100644 (file)
@@ -1,25 +1,13 @@
-/*  Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]
-    Radboud University Nijmegen, The Netherlands
-    
-    This file is part of libDAI.
-
-    libDAI is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    libDAI is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with libDAI; if not, write to the Free Software
-    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
-*/
+/*  This file is part of libDAI - http://www.libdai.org/
+ *
+ *  Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
+ *
+ *  Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
+ */
 
 
 #include <iostream>
+#include <iomanip>
 #include <iterator>
 #include <map>
 #include <set>
 #include <string>
 #include <algorithm>
 #include <functional>
-#include <tr1/unordered_map>
 #include <dai/factorgraph.h>
+#include <dai/util.h>
+#include <dai/exceptions.h>
+#include <boost/lexical_cast.hpp>
 
 
 namespace dai {
@@ -37,66 +27,50 @@ namespace dai {
 using namespace std;
 
 
-FactorGraph::FactorGraph( const vector<Factor> &P ) : BipartiteGraph<Var,Factor>(), _undoProbs(), _normtype(Prob::NORMPROB) {
+FactorGraph::FactorGraph( const std::vector<Factor> &P ) : _G(), _backup() {
     // add factors, obtain variables
-    set<Var> _vars;
-    V2s().reserve( P.size() );
+    set<Var> varset;
+    _factors.reserve( P.size() );
     size_t nrEdges = 0;
     for( vector<Factor>::const_iterator p2 = P.begin(); p2 != P.end(); p2++ ) {
-        V2s().push_back( *p2 );
-       copy( p2->vars().begin(), p2->vars().end(), inserter( _vars, _vars.begin() ) );
-       nrEdges += p2->vars().size();
+        _factors.push_back( *p2 );
+        copy( p2->vars().begin(), p2->vars().end(), inserter( varset, varset.begin() ) );
+        nrEdges += p2->vars().size();
     }
 
-    // add _vars
-    V1s().reserve( _vars.size() );
-    for( set<Var>::const_iterator p1 = _vars.begin(); p1 != _vars.end(); p1++ )
-        V1s().push_back( *p1 );
-    
+    // add vars
+    _vars.reserve( varset.size() );
+    for( set<Var>::const_iterator p1 = varset.begin(); p1 != varset.end(); p1++ )
+        _vars.push_back( *p1 );
+
     // create graph structure
-    createGraph( nrEdges );
+    constructGraph( nrEdges );
 }
 
 
-/// Part of constructors (creates edges, neighbours and adjacency matrix)
-void FactorGraph::createGraph( size_t nrEdges ) {
+void FactorGraph::constructGraph( size_t nrEdges ) {
     // create a mapping for indices
-    std::tr1::unordered_map<size_t, size_t> hashmap;
-    
-    for( size_t i = 0; i < vars().size(); i++ )        
-       hashmap[vars()[i].label()] = i;
-    
-    // create edges
-    edges().reserve( nrEdges );
+    hash_map<size_t, size_t> hashmap;
+
+    for( size_t i = 0; i < vars().size(); i++ )
+        hashmap[var(i).label()] = i;
+
+    // create edge list
+    vector<Edge> edges;
+    edges.reserve( nrEdges );
     for( size_t i2 = 0; i2 < nrFactors(); i2++ ) {
         const VarSet& ns = factor(i2).vars();
         for( VarSet::const_iterator q = ns.begin(); q != ns.end(); q++ )
-            edges().push_back(_edge_t(hashmap[q->label()], i2));               
+            edges.push_back( Edge(hashmap[q->label()], i2) );
     }
 
-    // calc neighbours and adjacency matrix
-    Regenerate();
+    // create bipartite graph
+    _G.construct( nrVars(), nrFactors(), edges.begin(), edges.end() );
 }
 
 
-/*FactorGraph& FactorGraph::addFactor( const Factor &I ) {
-    // add Factor
-    _V2.push_back( I );
-
-    // add new vars in Factor
-    for( VarSet::const_iterator i = I.vars().begin(); i != I.vars().end(); i++ ) {
-        size_t i_ind = find(vars().begin(), vars().end(), *i) - vars().begin();
-        if( i_ind == vars().size() )
-            _V1.push_back( *i );
-        _E12.push_back( _edge_t( i_ind, nrFactors() - 1 ) );
-    }
-
-    Regenerate();
-    return(*this);
-}*/
-
-
-ostream& operator << (ostream& os, const FactorGraph& fg) {
+/// Writes a FactorGraph to an output stream
+std::ostream& operator<< ( std::ostream &os, const FactorGraph &fg ) {
     os << fg.nrFactors() << endl;
 
     for( size_t I = 0; I < fg.nrFactors(); I++ ) {
@@ -109,521 +83,461 @@ ostream& operator << (ostream& os, const FactorGraph& fg) {
             os << i->states() << " ";
         os << endl;
         size_t nr_nonzeros = 0;
-        for( size_t k = 0; k < fg.factor(I).stateSpace(); k++ )
-            if( fg.factor(I)[k] != 0.0 )
+        for( size_t k = 0; k < fg.factor(I).nrStates(); k++ )
+            if( fg.factor(I)[k] != (Real)0 )
                 nr_nonzeros++;
         os << nr_nonzeros << endl;
-        for( size_t k = 0; k < fg.factor(I).stateSpace(); k++ )
-            if( fg.factor(I)[k] != 0.0 ) {
-                char buf[20];
-                sprintf(buf,"%18.14g", fg.factor(I)[k]);
-                os << k << " " << buf << endl;
-            }
+        for( size_t k = 0; k < fg.factor(I).nrStates(); k++ )
+            if( fg.factor(I)[k] != (Real)0 )
+                os << k << " " << setw(os.precision()+4) << fg.factor(I)[k] << endl;
     }
 
     return(os);
 }
 
 
-istream& operator >> (istream& is, FactorGraph& fg) {
+/// Reads a FactorGraph from an input stream
+std::istream& operator>> ( std::istream& is, FactorGraph &fg ) {
     long verbose = 0;
 
-    try {
-        vector<Factor> factors;
-        size_t nr_f;
-        string line;
-        
+    vector<Factor> facs;
+    size_t nr_Factors;
+    string line;
+
+    while( (is.peek()) == '#' )
+        getline(is,line);
+    is >> nr_Factors;
+    if( is.fail() )
+        DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Cannot read number of factors");
+    if( verbose >= 1 )
+        cerr << "Reading " << nr_Factors << " factors..." << endl;
+
+    getline (is,line);
+    if( is.fail() || line.size() > 0 )
+        DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Expecting empty line");
+
+    map<long,size_t> vardims;
+    for( size_t I = 0; I < nr_Factors; I++ ) {
+        if( verbose >= 2 )
+            cerr << "Reading factor " << I << "..." << endl;
+        size_t nr_members;
         while( (is.peek()) == '#' )
             getline(is,line);
-        is >> nr_f;
-        if( is.fail() )
-            throw "ReadFromFile: unable to read number of Factors";
+        is >> nr_members;
         if( verbose >= 2 )
-            cout << "Reading " << nr_f << " factors..." << endl;
+            cerr << "  nr_members: " << nr_members << endl;
 
-        getline (is,line);
-        if( is.fail() )
-            throw "ReadFromFile: empty line expected";
-
-        for( size_t I = 0; I < nr_f; I++ ) {
-            if( verbose >= 3 )
-                cout << "Reading factor " << I << "..." << endl;
-            size_t nr_members;
+        vector<long> labels;
+        for( size_t mi = 0; mi < nr_members; mi++ ) {
+            long mi_label;
             while( (is.peek()) == '#' )
                 getline(is,line);
-            is >> nr_members;
-            if( verbose >= 3 )
-                cout << "  nr_members: " << nr_members << endl;
-
-            vector<long> labels;
-            for( size_t mi = 0; mi < nr_members; mi++ ) {
-                long mi_label;
-                while( (is.peek()) == '#' )
-                    getline(is,line);
-                is >> mi_label;
-                labels.push_back(mi_label);
-            }
-            if( verbose >= 3 ) {
-                cout << "  labels: ";
-                copy (labels.begin(), labels.end(), ostream_iterator<int>(cout, " "));
-                cout << endl;
-            }
+            is >> mi_label;
+            labels.push_back(mi_label);
+        }
+        if( verbose >= 2 )
+            cerr << "  labels: " << labels << endl;
 
-            vector<size_t> dims;
-            for( size_t mi = 0; mi < nr_members; mi++ ) {
-                size_t mi_dim;
-                while( (is.peek()) == '#' )
-                    getline(is,line);
-                is >> mi_dim;
-                dims.push_back(mi_dim);
-            }
-            if( verbose >= 3 ) {
-                cout << "  dimensions: ";
-                copy (dims.begin(), dims.end(), ostream_iterator<int>(cout, " "));
-                cout << endl;
-            }
+        vector<size_t> dims;
+        for( size_t mi = 0; mi < nr_members; mi++ ) {
+            size_t mi_dim;
+            while( (is.peek()) == '#' )
+                getline(is,line);
+            is >> mi_dim;
+            dims.push_back(mi_dim);
+        }
+        if( verbose >= 2 )
+            cerr << "  dimensions: " << dims << endl;
+
+        // add the Factor
+        vector<Var> Ivars;
+        Ivars.reserve( nr_members );
+        for( size_t mi = 0; mi < nr_members; mi++ ) {
+            map<long,size_t>::iterator vdi = vardims.find( labels[mi] );
+            if( vdi != vardims.end() ) {
+                // check whether dimensions are consistent
+                if( vdi->second != dims[mi] )
+                    DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Variable with label " + boost::lexical_cast<string>(labels[mi]) + " has inconsistent dimensions.");
+            } else
+                vardims[labels[mi]] = dims[mi];
+            Ivars.push_back( Var(labels[mi], dims[mi]) );
+        }
+        facs.push_back( Factor( VarSet( Ivars.begin(), Ivars.end(), Ivars.size() ), (Real)0 ) );
+        if( verbose >= 2 )
+            cerr << "  vardims: " << vardims << endl;
 
-            // add the Factor
-            VarSet I_vars;
-            for( size_t mi = 0; mi < nr_members; mi++ )
-                I_vars.insert( Var(labels[mi], dims[mi]) );
-            factors.push_back(Factor(I_vars,0.0));
-            
-            // calculate permutation sigma (internally, members are sorted)
-            vector<long> sigma(nr_members,0);
-            VarSet::iterator j = I_vars.begin();
-            for( size_t mi = 0; mi < nr_members; mi++,j++ ) {
-                long search_for = j->label();
-                vector<long>::iterator j_loc = find(labels.begin(),labels.end(),search_for);
-                sigma[mi] = j_loc - labels.begin();
-            }
-            if( verbose >= 3 ) {
-                cout << "  sigma: ";
-                copy( sigma.begin(), sigma.end(), ostream_iterator<int>(cout," "));
-                cout << endl;
-            }
+        // calculate permutation object
+        Permute permindex( Ivars );
 
-            // calculate multindices
-            vector<size_t> sdims(nr_members,0);
-            for( size_t k = 0; k < nr_members; k++ ) {
-                sdims[k] = dims[sigma[k]];
-            }
-            multind mi(dims);
-            multind smi(sdims);
-            if( verbose >= 3 ) {
-                cout << "  mi.max(): " << mi.max() << endl;
-                cout << "       ";
-                for( size_t k=0; k < nr_members; k++ ) 
-                    cout << labels[k] << " ";
-                cout << "   ";
-                for( size_t k=0; k < nr_members; k++ ) 
-                    cout << labels[sigma[k]] << " ";
-                cout << endl;
-            }
-            
-            // read values
-            size_t nr_nonzeros;
+        // read values
+        size_t nr_nonzeros;
+        while( (is.peek()) == '#' )
+            getline(is,line);
+        is >> nr_nonzeros;
+        if( verbose >= 2 )
+            cerr << "  nonzeroes: " << nr_nonzeros << endl;
+        for( size_t k = 0; k < nr_nonzeros; k++ ) {
+            size_t li;
+            Real val;
             while( (is.peek()) == '#' )
                 getline(is,line);
-            is >> nr_nonzeros;
-            if( verbose >= 3 ) 
-                cout << "  nonzeroes: " << nr_nonzeros << endl;
-            for( size_t k = 0; k < nr_nonzeros; k++ ) {
-                size_t li;
-                double val;
-                while( (is.peek()) == '#' )
-                    getline(is,line);
-                is >> li;
-                while( (is.peek()) == '#' )
-                    getline(is,line);
-                is >> val;
-
-                vector<size_t> vi = mi.vi(li);
-                vector<size_t> svi(vi.size(),0);
-                for( size_t k = 0; k < vi.size(); k++ )
-                    svi[k] = vi[sigma[k]];
-                size_t sli = smi.li(svi);
-                if( verbose >= 3 ) {
-                    cout << "    " << li << ": ";
-                    copy( vi.begin(), vi.end(), ostream_iterator<size_t>(cout," "));
-                    cout << "-> ";
-                    copy( svi.begin(), svi.end(), ostream_iterator<size_t>(cout," "));
-                    cout << ": " << sli << endl;
-                }
-                factors.back()[sli] = val;
-            }
-        }
+            is >> li;
+            while( (is.peek()) == '#' )
+                getline(is,line);
+            is >> val;
 
-        if( verbose >= 3 ) {
-            cout << "factors:" << endl;
-            copy(factors.begin(), factors.end(), ostream_iterator<Factor>(cout,"\n"));
+            // store value, but permute indices first according to internal representation
+            facs.back().set( permindex.convertLinearIndex( li ), val );
         }
-
-        fg = FactorGraph(factors);
-    } catch (char *e) {
-        cout << e << endl;
     }
 
+    if( verbose >= 3 )
+        cerr << "factors:" << facs << endl;
+
+    fg = FactorGraph(facs);
+
     return is;
 }
 
 
-VarSet FactorGraph::delta(const Var & n) const {
+VarSet FactorGraph::Delta( size_t i ) const {
     // calculate Markov Blanket
-    size_t i = findVar( n );
-
-    VarSet del;
-    for( _nb_cit I = nb1(i).begin(); I != nb1(i).end(); ++I )
-        for( _nb_cit j = nb2(*I).begin(); j != nb2(*I).end(); ++j )
-            if( *j != i )
-                del |= var(*j);
+    VarSet Del;
+    bforeach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
+        bforeach( const Neighbor &j, nbF(I) ) // for all neighboring variables j of I
+            Del |= var(j);
 
-    return del;
+    return Del;
 }
 
 
-VarSet FactorGraph::Delta(const Var & n) const {
-    return( delta(n) | n );
-}
-
-
-void FactorGraph::makeFactorCavity(size_t I) {
-    // fill Factor I with ones
-    factor(I).fill(1.0);
+VarSet FactorGraph::Delta( const VarSet &ns ) const {
+    VarSet result;
+    for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++ )
+        result |= Delta( findVar(*n) );
+    return result;
 }
 
 
-void FactorGraph::makeCavity(const Var & n) {
-    // fills all Factors that include Var n with ones
-    size_t i = findVar( n );
+SmallSet<size_t> FactorGraph::Deltai( size_t i ) const {
+    // calculate Markov Blanket
+    SmallSet<size_t> Del;
+    bforeach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
+        bforeach( const Neighbor &j, nbF(I) ) // for all neighboring variables j of I
+            Del |= j;
 
-    for( _nb_cit I = nb1(i).begin(); I != nb1(i).end(); ++I )
-        factor(*I).fill(1.0);
+    return Del;
 }
 
 
-bool FactorGraph::hasNegatives() const {
-    bool result = false;
-    for( size_t I = 0; I < nrFactors() && !result; I++ )
-        if( factor(I).hasNegatives() )
-           result = true;
-    return result;
-}
-
-/*FactorGraph & FactorGraph::DeleteFactor(size_t I) {
-    // Go through all edges
-    for( vector<_edge_t>::iterator edge = _E12.begin(); edge != _E12.end(); edge++ )
-        if( edge->second >= I ) {
-            if( edge->second == I )
-                edge->second = -1UL;
-            else 
-                (edge->second)--;
-        }
-    // Remove all edges containing I
-    for( vector<_edge_t>::iterator edge = _E12.begin(); edge != _E12.end(); edge++ )
-        if( edge->second == -1UL )
-            edge = _E12.erase( edge );
-//  vector<_edge_t>::iterator new_end = _E12.remove_if( _E12.begin(), _E12.end(), compose1( bind2nd(equal_to<size_t>(), -1), select2nd<_edge_t>() ) );
-//  _E12.erase( new_end, _E12.end() );
-
-    // Erase the factor
-    _V2.erase( _V2.begin() + I );
-    
-    Regenerate();
-
-    return *this;
+void FactorGraph::makeCavity( size_t i, bool backup ) {
+    // fills all Factors that include var(i) with ones
+    map<size_t,Factor> newFacs;
+    bforeach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
+        newFacs[I] = Factor( factor(I).vars(), (Real)1 );
+    setFactors( newFacs, backup );
 }
 
 
-FactorGraph & FactorGraph::DeleteVar(size_t i) {
-    // Go through all edges
-    for( vector<_edge_t>::iterator edge = _E12.begin(); edge != _E12.end(); edge++ )
-        if( edge->first >= i ) {
-            if( edge->first == i )
-                edge->first = -1UL;
-            else 
-                (edge->first)--;
-        }
-    // Remove all edges containing i
-    for( vector<_edge_t>::iterator edge = _E12.begin(); edge != _E12.end(); edge++ )
-        if( edge->first == -1UL )
-            edge = _E12.erase( edge );
-                
-//  vector<_edge_t>::iterator new_end = _E12.remove_if( _E12.begin(), _E12.end(), compose1( bind2nd(equal_to<size_t>(), -1), select1st<_edge_t>() ) );
-//  _E12.erase( new_end, _E12.end() );
-
-    // Erase the variable
-    _V1.erase( _V1.begin() + i );
-    
-    Regenerate();
-
-    return *this;
-}*/
+void FactorGraph::makeRegionCavity( std::vector<size_t> facInds, bool backup ) {
+    map<size_t,Factor> newFacs;
+    for( size_t I = 0; I < facInds.size(); I++ )
+       newFacs[facInds[I]] = Factor(factor(facInds[I]).vars(), (Real)1);
+    setFactors( newFacs, backup );
+}
 
 
-long FactorGraph::ReadFromFile(const char *filename) {
+void FactorGraph::ReadFromFile( const char *filename ) {
     ifstream infile;
-    infile.open (filename);
-    if (infile.is_open()) {
+    infile.open( filename );
+    if( infile.is_open() ) {
         infile >> *this;
         infile.close();
-        return 0;
-    } else {
-        cout << "ERROR OPENING FILE" << endl;
-        return 1;
-    }
+    } else
+        DAI_THROWE(CANNOT_READ_FILE,"Cannot read from file " + std::string(filename));
 }
 
 
-long FactorGraph::WriteToFile(const char *filename) const {
+void FactorGraph::WriteToFile( const char *filename, size_t precision ) const {
     ofstream outfile;
-    outfile.open (filename);
-    if (outfile.is_open()) {
-        try {
-            outfile << *this;
-        } catch (char *e) {
-            cout << e << endl;
-            return 1;
-        }
+    outfile.open( filename );
+    if( outfile.is_open() ) {
+        outfile.precision( precision );
+        outfile << *this;
         outfile.close();
-        return 0;
-    } else {
-        cout << "ERROR OPENING FILE" << endl;
-        return 1;
-    }
+    } else
+        DAI_THROWE(CANNOT_WRITE_FILE,"Cannot write to file " + std::string(filename));
 }
 
 
-long FactorGraph::WriteToDotFile(const char *filename) const {
-    ofstream outfile;
-    outfile.open (filename);
-    if (outfile.is_open()) {
-        try {
-            outfile << "graph G {" << endl;
-            outfile << "graph[size=\"9,9\"];" << endl;
-            outfile << "node[shape=circle,width=0.4,fixedsize=true];" << endl;
-            for( size_t i = 0; i < nrVars(); i++ )
-                outfile << "\tx" << var(i).label() << ";" << endl;
-            outfile << "node[shape=box,style=filled,color=lightgrey,width=0.3,height=0.3,fixedsize=true];" << endl;
-            for( size_t I = 0; I < nrFactors(); I++ )
-                outfile << "\tp" << I << ";" << endl;
-            for( size_t iI = 0; iI < nr_edges(); iI++ )
-                outfile << "\tx" << var(edge(iI).first).label() << " -- p" << edge(iI).second << ";" << endl;
-            outfile << "}" << endl;
-        } catch (char *e) {
-            cout << e << endl;
-            return 1;
-        }
-        outfile.close();
-        return 0;
-    } else {
-        cout << "ERROR OPENING FILE" << endl;
-        return 1;
-    }
+void FactorGraph::printDot( std::ostream &os ) const {
+    os << "graph FactorGraph {" << endl;
+    os << "node[shape=circle,width=0.4,fixedsize=true];" << endl;
+    for( size_t i = 0; i < nrVars(); i++ )
+        os << "\tv" << var(i).label() << ";" << endl;
+    os << "node[shape=box,width=0.3,height=0.3,fixedsize=true];" << endl;
+    for( size_t I = 0; I < nrFactors(); I++ )
+        os << "\tf" << I << ";" << endl;
+    for( size_t i = 0; i < nrVars(); i++ )
+        bforeach( const Neighbor &I, nbV(i) )  // for all neighboring factors I of i
+            os << "\tv" << var(i).label() << " -- f" << I << ";" << endl;
+    os << "}" << endl;
+}
+
+
+GraphAL FactorGraph::MarkovGraph() const {
+    GraphAL G( nrVars() );
+    for( size_t i = 0; i < nrVars(); i++ )
+        bforeach( const Neighbor &I, nbV(i) )
+            bforeach( const Neighbor &j, nbF(I) )
+                if( i < j )
+                    G.addEdge( i, j, true );
+    return G;
 }
 
 
-bool hasShortLoops( const vector<Factor> &P ) {
-    bool found = false;
-    vector<Factor>::const_iterator I, J;
-    for( I = P.begin(); I != P.end(); I++ ) {
-        J = I;
-        J++;
-        for( ; J != P.end(); J++ )
-            if( (I->vars() & J->vars()).size() >= 2 ) {
-                found = true;
-                break;
+bool FactorGraph::isMaximal( size_t I ) const {
+    const VarSet& I_vars = factor(I).vars();
+    size_t I_size = I_vars.size();
+
+    if( I_size == 0 ) {
+        for( size_t J = 0; J < nrFactors(); J++ ) 
+            if( J != I )
+                if( factor(J).vars().size() > 0 )
+                    return false;
+        return true;
+    } else {
+        bforeach( const Neighbor& i, nbF(I) ) {
+            bforeach( const Neighbor& J, nbV(i) ) {
+                if( J != I )
+                    if( (factor(J).vars() >> I_vars) && (factor(J).vars().size() != I_size) )
+                        return false;
             }
-        if( found )
-            break;
+        }
+        return true;
     }
-    return found;
 }
 
 
-void RemoveShortLoops(vector<Factor> &P) {
-    bool found = true;
-    while( found ) {
-        found = false;
-        vector<Factor>::iterator I, J;
-        for( I = P.begin(); I != P.end(); I++ ) {
-            J = I;
-            J++;
-            for( ; J != P.end(); J++ )
-                if( (I->vars() & J->vars()).size() >= 2 ) {
-                    found = true;
-                    break;
-                }
-            if( found )
-                break;
-        }
-        if( found ) {
-            cout << "Merging factors " << I->vars() << " and " << J->vars() << endl;
-            *I *= *J;
-            P.erase(J);
+size_t FactorGraph::maximalFactor( size_t I ) const {
+    const VarSet& I_vars = factor(I).vars();
+    size_t I_size = I_vars.size();
+
+    if( I_size == 0 ) {
+        for( size_t J = 0; J < nrFactors(); J++ )
+            if( J != I )
+                if( factor(J).vars().size() > 0 )
+                    return maximalFactor( J );
+        return I;
+    } else {
+        bforeach( const Neighbor& i, nbF(I) ) {
+            bforeach( const Neighbor& J, nbV(i) ) {
+                if( J != I )
+                    if( (factor(J).vars() >> I_vars) && (factor(J).vars().size() != I_size) )
+                        return maximalFactor( J );
+            }
         }
+        return I;
     }
 }
 
 
-Factor FactorGraph::ExactMarginal(const VarSet & x) const {
-    Factor P;
+vector<VarSet> FactorGraph::maximalFactorDomains() const {
+    vector<VarSet> result;
+
     for( size_t I = 0; I < nrFactors(); I++ )
-        P *= factor(I);
-    return P.marginal(x);
+        if( isMaximal( I ) )
+            result.push_back( factor(I).vars() );
+
+    if( result.size() == 0 )
+        result.push_back( VarSet() );
+    return result;
 }
 
 
-Real FactorGraph::ExactlogZ() const {
-    Factor P;
+Real FactorGraph::logScore( const std::vector<size_t>& statevec ) const {
+    // Construct a State object that represents statevec
+    // This decouples the representation of the joint state in statevec from the factor graph
+    map<Var, size_t> statemap;
+    for( size_t i = 0; i < statevec.size(); i++ )
+        statemap[var(i)] = statevec[i];
+    State S(statemap);
+
+    // Evaluate the log probability of the joint configuration in statevec
+    // by summing the log factor entries of the factors that correspond to this joint configuration
+    Real lS = 0.0;
     for( size_t I = 0; I < nrFactors(); I++ )
-        P *= factor(I);
-    return std::log(P.totalSum());
+        lS += dai::log( factor(I)[BigInt_size_t(S(factor(I).vars()))] );
+    return lS;
 }
 
 
-vector<VarSet> FactorGraph::Cliques() const {
-    vector<VarSet> result;
-    
-    for( size_t I = 0; I < nrFactors(); I++ ) {
-        bool maximal = true;
-        for( size_t J = 0; (J < nrFactors()) && maximal; J++ )
-            if( (factor(J).vars() >> factor(I).vars()) && !(factor(J).vars() == factor(I).vars()) )
-                maximal = false;
-        
-        if( maximal )
-            result.push_back( factor(I).vars() );
-    }
+void FactorGraph::clamp( size_t i, size_t x, bool backup ) {
+    DAI_ASSERT( x <= var(i).states() );
+    Factor mask( var(i), (Real)0 );
+    mask.set( x, (Real)1 );
 
-    return result;
-}
+    map<size_t, Factor> newFacs;
+    bforeach( const Neighbor &I, nbV(i) )
+        newFacs[I] = factor(I) * mask;
+    setFactors( newFacs, backup );
 
+    return;
+}
 
-void FactorGraph::clamp( const Var & n, size_t i ) {
-    assert( i <= n.states() );
-
-/*  if( do_surgery ) {
-        size_t ni = find( vars().begin(), vars().end(), n) - vars().begin();
-
-        if( ni != nrVars() ) {
-            for( _nb_cit I = nb1(ni).begin(); I != nb1(ni).end(); I++ ) {
-                if( factor(*I).size() == 1 )
-                    // Remove this single-variable factor
-    //              I = (_V2.erase(I))--;
-                    _E12.erase( _E12.begin() + VV2E(ni, *I) );
-                else {
-                    // Replace it by the slice
-                    Index ind_I_min_n( factor(*I), factor(*I) / n );
-                    Index ind_n( factor(*I), n );
-                    Factor slice_I( factor(*I) / n );
-                    for( size_t ind_I = 0; ind_I < factor(*I).stateSpace(); ++ind_I, ++ind_I_min_n, ++ind_n )
-                        if( ind_n == i )
-                            slice_I[ind_I_min_n] = factor(*I)[ind_I];
-                    factor(*I) = slice_I;
-
-                    // Remove the edge between n and I
-                    _E12.erase( _E12.begin() + VV2E(ni, *I) );
-                }
-            }
 
-            Regenerate();
-            
-            // remove all unconnected factors
-            for( size_t I = 0; I < nrFactors(); I++ )
-                if( nb2(I).size() == 0 )
-                    DeleteFactor(I--);
+void FactorGraph::clampVar( size_t i, const vector<size_t> &is, bool backup ) {
+    Var n = var(i);
+    Factor mask_n( n, (Real)0 );
 
-            DeleteVar( ni );
+    bforeach( size_t i, is ) {
+        DAI_ASSERT( i <= n.states() );
+        mask_n.set( i, (Real)1 );
+    }
 
-            // FIXME
-        }
-    } */
+    map<size_t, Factor> newFacs;
+    bforeach( const Neighbor &I, nbV(i) )
+        newFacs[I] = factor(I) * mask_n;
+    setFactors( newFacs, backup );
+}
 
-    // The cheap solution (at least in terms of coding time) is to multiply every factor
-    // that contains the variable with a delta function
 
-    Factor delta_n_i(n,0.0);
-    delta_n_i[i] = 1.0;
+void FactorGraph::clampFactor( size_t I, const vector<size_t> &is, bool backup ) {
+    size_t st = factor(I).nrStates();
+    Factor newF( factor(I).vars(), (Real)0 );
 
-    // For all factors that contain n
-    for( size_t I = 0; I < nrFactors(); I++ ) 
-        if( factor(I).vars() && n )
-            // Multiply it with a delta function
-            factor(I) *= delta_n_i;
+    bforeach( size_t i, is ) {
+        DAI_ASSERT( i <= st );
+        newF.set( i, factor(I)[i] );
+    }
 
-    return;
+    setFactor( I, newF, backup );
 }
 
 
-void FactorGraph::saveProb( size_t I ) {
-    map<size_t,Prob>::iterator it = _undoProbs.find( I );
-    if( it != _undoProbs.end() )
-        cout << "FactorGraph::saveProb:  WARNING: _undoProbs[I] already defined!" << endl;
-    _undoProbs[I] = factor(I).p();
+void FactorGraph::backupFactor( size_t I ) {
+    map<size_t,Factor>::iterator it = _backup.find( I );
+    if( it != _backup.end() )
+        DAI_THROW(MULTIPLE_UNDO);
+    _backup[I] = factor(I);
 }
 
 
-void FactorGraph::undoProb( size_t I ) {
-    map<size_t,Prob>::iterator it = _undoProbs.find( I );
-    if( it != _undoProbs.end() ) {
-        factor(I).p() = (*it).second;
-        _undoProbs.erase(it);
-    }
+void FactorGraph::restoreFactor( size_t I ) {
+    map<size_t,Factor>::iterator it = _backup.find( I );
+    if( it != _backup.end() ) {
+        setFactor(I, it->second);
+        _backup.erase(it);
+    } else
+        DAI_THROW(OBJECT_NOT_FOUND);
 }
 
 
-void FactorGraph::saveProbs( const VarSet &ns ) {
-    if( !_undoProbs.empty() )
-        cout << "FactorGraph::saveProbs:  WARNING: _undoProbs not empy!" << endl;
+void FactorGraph::backupFactors( const VarSet &ns ) {
     for( size_t I = 0; I < nrFactors(); I++ )
-        if( factor(I).vars() && ns )
-            _undoProbs[I] = factor(I).p();
+        if( factor(I).vars().intersects( ns ) )
+            backupFactor( I );
 }
 
 
-void FactorGraph::undoProbs( const VarSet &ns ) {
-    for( map<size_t,Prob>::iterator uI = _undoProbs.begin(); uI != _undoProbs.end(); ) {
-        if( factor((*uI).first).vars() && ns ) {
-//          cout << "undoing " << factor((*uI).first).vars() << endl;
-//          cout << "from " << factor((*uI).first).p() << " to " << (*uI).second << endl;
-            factor((*uI).first).p() = (*uI).second;
-            _undoProbs.erase(uI++);
+void FactorGraph::restoreFactors( const VarSet &ns ) {
+    map<size_t,Factor> facs;
+    for( map<size_t,Factor>::iterator uI = _backup.begin(); uI != _backup.end(); ) {
+        if( factor(uI->first).vars().intersects( ns ) ) {
+            facs.insert( *uI );
+            _backup.erase(uI++);
         } else
             uI++;
     }
+    setFactors( facs );
+}
+
+
+void FactorGraph::restoreFactors() {
+    setFactors( _backup );
+    _backup.clear();
+}
+
+
+void FactorGraph::backupFactors( const std::set<size_t> & facs ) {
+    for( std::set<size_t>::const_iterator fac = facs.begin(); fac != facs.end(); fac++ )
+        backupFactor( *fac );
 }
 
 
-bool FactorGraph::isConnected() const {
-    if( nrVars() == 0 )
-        return false;
-    else {
-        Var n = var( 0 );
+bool FactorGraph::isPairwise() const {
+    bool pairwise = true;
+    for( size_t I = 0; I < nrFactors() && pairwise; I++ )
+        if( factor(I).vars().size() > 2 )
+            pairwise = false;
+    return pairwise;
+}
 
-        VarSet component = n;
 
-        VarSet remaining;
-        for( size_t i = 1; i < nrVars(); i++ )
-            remaining |= var(i);
+bool FactorGraph::isBinary() const {
+    bool binary = true;
+    for( size_t i = 0; i < nrVars() && binary; i++ )
+        if( var(i).states() > 2 )
+            binary = false;
+    return binary;
+}
 
-        bool found_new_vars = true;
-        while( found_new_vars ) {
-            VarSet new_vars;
-            for( VarSet::const_iterator m = remaining.begin(); m != remaining.end(); m++ )
-                if( delta(*m) && component )
-                    new_vars |= *m;
 
-            if( new_vars.empty() )
-                found_new_vars = false;
-            else 
-                found_new_vars = true;
+FactorGraph FactorGraph::clamped( size_t i, size_t state ) const {
+    Var v = var( i );
+    Real zeroth_order = (Real)1;
+    vector<Factor> clamped_facs;
+    clamped_facs.push_back( createFactorDelta( v, state ) );
+    for( size_t I = 0; I < nrFactors(); I++ ) {
+        VarSet v_I = factor(I).vars();
+        Factor new_factor;
+        if( v_I.intersects( v ) )
+            new_factor = factor(I).slice( v, state );
+        else
+            new_factor = factor(I);
+
+        if( new_factor.vars().size() != 0 ) {
+            size_t J = 0;
+            // if it can be merged with a previous one, do that
+            for( J = 0; J < clamped_facs.size(); J++ )
+                if( clamped_facs[J].vars() == new_factor.vars() ) {
+                    clamped_facs[J] *= new_factor;
+                    break;
+                }
+            // otherwise, push it back
+            if( J == clamped_facs.size() || clamped_facs.size() == 0 )
+                clamped_facs.push_back( new_factor );
+        } else
+            zeroth_order *= new_factor[0];
+    }
+    *(clamped_facs.begin()) *= zeroth_order;
+    return FactorGraph( clamped_facs );
+}
 
-            component |= new_vars;
-            remaining /= new_vars;
-        };
-        return remaining.empty();
+
+FactorGraph FactorGraph::maximalFactors() const {
+    vector<size_t> maxfac( nrFactors() );
+    map<size_t,size_t> newindex;
+    size_t nrmax = 0;
+    for( size_t I = 0; I < nrFactors(); I++ ) {
+        maxfac[I] = I;
+        VarSet maxfacvars = factor(maxfac[I]).vars();
+        for( size_t J = 0; J < nrFactors(); J++ ) {
+            VarSet Jvars = factor(J).vars();
+            if( Jvars >> maxfacvars && (Jvars != maxfacvars) ) {
+                maxfac[I] = J;
+                maxfacvars = factor(maxfac[I]).vars();
+            }
+        }
+        if( maxfac[I] == I )
+            newindex[I] = nrmax++;
     }
+
+    vector<Factor> facs( nrmax );
+    for( size_t I = 0; I < nrFactors(); I++ )
+        facs[newindex[maxfac[I]]] *= factor(I);
+
+    return FactorGraph( facs.begin(), facs.end(), vars().begin(), vars().end(), facs.size(), nrVars() );
 }