Merged prob.h, factorgraph.h, factograph.cpp from SVN head (broken!)

diff --git a/ChangeLog b/ChangeLog
index 310044a..e003b28 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -84,6 +84,11 @@ libDAI-0.2.2 (2008-??-??)
       moved isConnected() to BipartiteGraph
       removed updatedFactor(size_t)
       removed _normtype and NormType()
       moved isConnected() to BipartiteGraph
       removed updatedFactor(size_t)
       removed _normtype and NormType()

+      WriteToDotFile(const char *filename) -> display( std::ostream& os )
+      undoProb(size_t) -> restoreFactor(size_t)
+      saveProb(size_t) -> backupFactor(size_t)
+      undoProbs(const VarSet &) -> restoreFactors(const VarSet &)
+      saveProbs(const VarSet &) -> backupFactors(const VarSet &)

   - RegionGraph::
       nr_ORs() -> nrORs()
       nr_IRs() -> nrIRs()
   - RegionGraph::
       nr_ORs() -> nrORs()
       nr_IRs() -> nrIRs()

diff --git a/STATUS b/STATUS
index 97c59a1..db2f534 100644 (file)
--- a/STATUS
+++ b/STATUS
@@ -53,19 +53,19 @@ alldai.cpp
 properties.h
 properties.cpp
 factor.h
 properties.h
 properties.cpp
 factor.h

+prob.h
+factorgraph.h
+factorgraph.cpp

 
 FILES IN SVN HEAD THAT ARE STILL RELEVANT:
 ChangeLog
 README
 TODO
 
 FILES IN SVN HEAD THAT ARE STILL RELEVANT:
 ChangeLog
 README
 TODO

-prob.h
-factorgraph.h
-factorgraph.cpp

 regiongraph.h
 regiongraph.cpp
 regiongraph.h
 regiongraph.cpp

-

 daialg.h
 daialg.cpp
 daialg.h
 daialg.cpp

+

 bp.h
 bp.cpp
 hak.h
 bp.h
 bp.cpp
 hak.h

diff --git a/include/dai/factorgraph.h b/include/dai/factorgraph.h
index 567d9bb..0659ded 100644 (file)
--- a/include/dai/factorgraph.h
+++ b/include/dai/factorgraph.h
@@ -32,10 +32,6 @@
 namespace dai {
 
 
 namespace dai {
 
 

-bool hasShortLoops( const std::vector<Factor> &P );
-void RemoveShortLoops( std::vector<Factor> &P );
-
-

 class FactorGraph {
     public:
         BipartiteGraph         G;
 class FactorGraph {
     public:
         BipartiteGraph         G;


@@ -45,14 +41,14 @@ class FactorGraph {
         typedef BipartiteGraph::Neighbors Neighbors;
         typedef BipartiteGraph::Edge      Edge;
 
         typedef BipartiteGraph::Neighbors Neighbors;
         typedef BipartiteGraph::Edge      Edge;
 

-    protected:
-        std::map<size_t,Prob>  _undoProbs;
+    private:
+        std::map<size_t,Factor>  _backupFactors;

 
     public:
         /// Default constructor
 
     public:
         /// Default constructor

-        FactorGraph() : G(), vars(), factors(), _undoProbs() {};
+        FactorGraph() : G(), vars(), factors(), _backupFactors() {}

         /// Copy constructor
         /// Copy constructor

-        FactorGraph(const FactorGraph & x) : G(x.G), vars(x.vars), factors(x.factors), _undoProbs(x._undoProbs) {};
+        FactorGraph(const FactorGraph & x) : G(x.G), vars(x.vars), factors(x.factors), _backupFactors(x._backupFactors) {}

         /// Construct FactorGraph from vector of Factors
         FactorGraph(const std::vector<Factor> &P);
         // Construct a FactorGraph from given factor and variable iterators
         /// Construct FactorGraph from vector of Factors
         FactorGraph(const std::vector<Factor> &P);
         // Construct a FactorGraph from given factor and variable iterators


@@ -65,19 +61,33 @@ class FactorGraph {
                 G          = x.G;
                 vars       = x.vars;
                 factors    = x.factors;
                 G          = x.G;
                 vars       = x.vars;
                 factors    = x.factors;

-                _undoProbs = x._undoProbs;
+                _backupFactors = x._backupFactors;

             }
             return *this;
         }
         virtual ~FactorGraph() {}
 
             }
             return *this;
         }
         virtual ~FactorGraph() {}
 

+        /// Create (virtual default constructor)
+        virtual FactorGraph* create() const {
+            return new FactorGraph(*this);
+        }
+
+        /// Clone (virtual copy constructor)
+        virtual FactorGraph* clone() const {
+            return new FactorGraph();
+        }
+

         // aliases
         Var & var(size_t i) { return vars[i]; }
         // aliases
         Var & var(size_t i) { return vars[i]; }

+        /// Get const reference to i'th variable

         const Var & var(size_t i) const { return vars[i]; }
         Factor & factor(size_t I) { return factors[I]; }
         const Var & var(size_t i) const { return vars[i]; }
         Factor & factor(size_t I) { return factors[I]; }

+        /// Get const reference to I'th factor

         const Factor & factor(size_t I) const { return factors[I]; }
 
         const Factor & factor(size_t I) const { return factors[I]; }
 

+        /// Get number of variables

         size_t nrVars() const { return vars.size(); }
         size_t nrVars() const { return vars.size(); }

+        /// Get number of factors

         size_t nrFactors() const { return factors.size(); }
         size_t nrEdges() const { return G.nrEdges(); }
 
         size_t nrFactors() const { return factors.size(); }
         size_t nrEdges() const { return G.nrEdges(); }
 


@@ -98,11 +108,22 @@ class FactorGraph {
         /// Provides full access to neighbor of factor
         Neighbor & nbF( size_t I, size_t _i ) { return G.nb2(I)[_i]; }
 
         /// Provides full access to neighbor of factor
         Neighbor & nbF( size_t I, size_t _i ) { return G.nb2(I)[_i]; }
 

-        size_t findVar(const Var & n) const {
+        /// Get index of variable n
+        size_t findVar( const Var & n ) const {

             size_t i = find( vars.begin(), vars.end(), n ) - vars.begin();
             assert( i != nrVars() );
             return i;
         }
             size_t i = find( vars.begin(), vars.end(), n ) - vars.begin();
             assert( i != nrVars() );
             return i;
         }

+
+        /// Get set of indexes for set of variables
+        std::set<size_t> findVars( VarSet &ns ) const {
+            std::set<size_t> indexes;
+            for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++ )
+                indexes.insert( findVar( *n ) );
+            return indexes;
+        }
+
+        /// Get index of first factor involving ns

         size_t findFactor(const VarSet &ns) const {
             size_t I;
             for( I = 0; I < nrFactors(); I++ )
         size_t findFactor(const VarSet &ns) const {
             size_t I;
             for( I = 0; I < nrFactors(); I++ )


@@ -112,29 +133,77 @@ class FactorGraph {
             return I;
         }
 
             return I;
         }
 

-        friend std::ostream& operator << (std::ostream& os, const FactorGraph& fg);
-        friend std::istream& operator >> (std::istream& is, FactorGraph& fg);
+        /// Return all variables that occur in a factor involving the i'th variable, itself included
+        VarSet Delta( unsigned i ) const;

 
 

+        /// Return all variables that occur in a factor involving some variable in ns, ns itself included
+        VarSet Delta( const VarSet &ns ) const;
+
+        /// Return all variables that occur in a factor involving the i'th variable, n itself excluded

         VarSet delta( unsigned i ) const;
         VarSet delta( unsigned i ) const;

-        VarSet Delta( unsigned i ) const;
-        virtual void makeCavity( unsigned i );

 
 

-        long ReadFromFile(const char *filename);
-        long WriteToFile(const char *filename) const;
-        long WriteToDotFile(const char *filename) const;
+        /// Return all variables that occur in a factor involving some variable in ns, ns itself excluded
+        VarSet delta( const VarSet & ns ) const {
+            return Delta( ns ) / ns;
+        }

 
 

-        virtual void clamp( const Var & n, size_t i );
-        
-        bool hasNegatives() const;
+        /// Set the content of the I'th factor and make a backup of its old content if backup == true
+        virtual void setFactor( size_t I, const Factor &newFactor, bool backup = false ) {
+            assert( newFactor.vars() == factors[I].vars() ); 
+            if( backup )
+                backupFactor( I );
+            factors[I] = newFactor; 
+        }
+
+        /// Set the contents of all factors as specified by facs and make a backup of the old contents if backup == true
+        virtual void setFactors( const std::map<size_t, Factor> & facs, bool backup = false ) {
+            for( std::map<size_t, Factor>::const_iterator fac = facs.begin(); fac != facs.end(); fac++ ) {
+                if( backup )
+                    backupFactor( fac->first );
+                setFactor( fac->first, fac->second );
+            }
+        }
+
+        /// Clamp variable n to value i (i.e. multiply with a Kronecker delta \f$\delta_{x_n, i}\f$);
+        /// If backup == true, make a backup of all factors that are changed
+        virtual void clamp( const Var & n, size_t i, bool backup = false );
+
+        /// Set all factors interacting with the i'th variable 1
+        virtual void makeCavity( unsigned i, bool backup = false );
+
+        /// Backup the factors specified by indices in facs
+        virtual void backupFactors( const std::set<size_t> & facs );
+
+        /// Restore all factors to the backup copies
+        virtual void restoreFactors();
+
+        bool isConnected() const { return G.isConnected(); }
+        bool isTree() const { return G.isTree(); }
+
+        friend std::ostream& operator << (std::ostream& os, const FactorGraph& fg);
+        friend std::istream& operator >> (std::istream& is, FactorGraph& fg);
+
+        void ReadFromFile(const char *filename);
+        void WriteToFile(const char *filename) const;
+        void display( std::ostream& os ) const;

         
         std::vector<VarSet> Cliques() const;
 
         
         std::vector<VarSet> Cliques() const;
 

-        virtual void undoProbs( const VarSet &ns );
-        void saveProbs( const VarSet &ns );
-        virtual void undoProb( size_t I );
-        void saveProb( size_t I );
+        // Clamp variable v_i to value state (i.e. multiply with a Kronecker delta \f$\delta_{x_{v_i},x}\f$);
+        // This version changes the factor graph structure and thus returns a newly constructed FactorGraph
+        // and keeps the current one constant, contrary to clamp()
+        FactorGraph clamped( const Var & v_i, size_t x ) const;
+
+        FactorGraph maximalFactors() const;
+
+        bool isPairwise() const;
+        bool isBinary() const;

 
     private:
 
     private:

+        void restoreFactor( size_t I );
+        void backupFactor( size_t I );
+        void restoreFactors( const VarSet &ns );
+        void backupFactors( const VarSet &ns );

         /// Part of constructors (creates edges, neighbors and adjacency matrix)
         void createGraph( size_t nrEdges );
 };
         /// Part of constructors (creates edges, neighbors and adjacency matrix)
         void createGraph( size_t nrEdges );
 };


@@ -142,7 +211,7 @@ class FactorGraph {
 
 // assumes that the set of variables in [var_begin,var_end) is the union of the variables in the factors in [fact_begin, fact_end)
 template<typename FactorInputIterator, typename VarInputIterator>
 
 // assumes that the set of variables in [var_begin,var_end) is the union of the variables in the factors in [fact_begin, fact_end)
 template<typename FactorInputIterator, typename VarInputIterator>

-FactorGraph::FactorGraph(FactorInputIterator fact_begin, FactorInputIterator fact_end, VarInputIterator var_begin, VarInputIterator var_end, size_t nr_fact_hint, size_t nr_var_hint ) : G(), _undoProbs() {
+FactorGraph::FactorGraph(FactorInputIterator fact_begin, FactorInputIterator fact_end, VarInputIterator var_begin, VarInputIterator var_end, size_t nr_fact_hint, size_t nr_var_hint ) : G(), _backupFactors() {

     // add factors
     size_t nrEdges = 0;
     factors.reserve( nr_fact_hint );
     // add factors
     size_t nrEdges = 0;
     factors.reserve( nr_fact_hint );


@@ -150,7 +219,7 @@ FactorGraph::FactorGraph(FactorInputIterator fact_begin, FactorInputIterator fac
         factors.push_back( *p2 );
         nrEdges += p2->vars().size();
     }
         factors.push_back( *p2 );
         nrEdges += p2->vars().size();
     }

- 
+

     // add variables
     vars.reserve( nr_var_hint );
     for( VarInputIterator p1 = var_begin; p1 != var_end; ++p1 )
     // add variables
     vars.reserve( nr_var_hint );
     for( VarInputIterator p1 = var_begin; p1 != var_end; ++p1 )

diff --git a/include/dai/prob.h b/include/dai/prob.h
index 6037339..b0466d2 100644 (file)
--- a/include/dai/prob.h
+++ b/include/dai/prob.h
@@ -42,6 +42,11 @@ template<typename T> class      TProb;
 typedef TProb<Real>             Prob;
 
 
 typedef TProb<Real>             Prob;
 
 

+// predefine friends
+template<typename T> TProb<T> min( const TProb<T> &a, const TProb<T> &b );
+template<typename T> TProb<T> max( const TProb<T> &a, const TProb<T> &b );
+
+

 /// TProb<T> implements a probability vector of type T.
 /// T should be castable from and to double.
 template <typename T> class TProb {
 /// TProb<T> implements a probability vector of type T.
 /// T should be castable from and to double.
 template <typename T> class TProb {


@@ -63,7 +68,7 @@ template <typename T> class TProb {
         typedef enum { DISTL1, DISTLINF, DISTTV } DistType;
         
         /// Default constructor
         typedef enum { DISTL1, DISTLINF, DISTTV } DistType;
         
         /// Default constructor

-        TProb() {}
+        TProb() : _p() {}

         
         /// Construct uniform distribution of given length
         explicit TProb( size_t n ) : _p(std::vector<T>(n, 1.0 / n)) {}
         
         /// Construct uniform distribution of given length
         explicit TProb( size_t n ) : _p(std::vector<T>(n, 1.0 / n)) {}


@@ -81,7 +86,13 @@ template <typename T> class TProb {
         std::vector<T> & p() { return _p; }
         
         /// Provide read access to ith element of _p
         std::vector<T> & p() { return _p; }
         
         /// Provide read access to ith element of _p

-        T operator[]( size_t i ) const { return _p[i]; }
+        T operator[]( size_t i ) const { 
+#ifdef DAI_DEBUG
+            return _p.at(i);
+#else
+            return _p[i];
+#endif
+        }

         
         /// Provide full access to ith element of _p
         T& operator[]( size_t i ) { return _p[i]; }
         
         /// Provide full access to ith element of _p
         T& operator[]( size_t i ) { return _p[i]; }


@@ -112,6 +123,14 @@ template <typename T> class TProb {
             return *this;
         }
 
             return *this;
         }
 

+        /// Make entries epsilon if they are smaller than epsilon
+        TProb<T>& makePositive (Real epsilon) {
+            for( size_t i = 0; i < size(); i++ )
+                if( (0 < (Real)_p[i]) && ((Real)_p[i] < epsilon) )
+                    _p[i] = epsilon;
+            return *this;
+        }
+

         /// Multiplication with T x
         TProb<T>& operator*= (T x) {
             std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::multiplies<T>(), x) );
         /// Multiplication with T x
         TProb<T>& operator*= (T x) {
             std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::multiplies<T>(), x) );


@@ -140,7 +159,7 @@ template <typename T> class TProb {
             quot /= x;
             return quot;
         }
             quot /= x;
             return quot;
         }

-        
+

         /// addition of x
         TProb<T>& operator+= (T x) {
             std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::plus<T>(), x ) );
         /// addition of x
         TProb<T>& operator+= (T x) {
             std::transform( _p.begin(), _p.end(), _p.begin(), std::bind2nd( std::plus<T>(), x ) );


@@ -167,6 +186,17 @@ template <typename T> class TProb {
             return diff;
         }
 
             return diff;
         }
 

+        /// Pointwise comparison
+        bool operator<= (const TProb<T> & q) const {
+#ifdef DAI_DEBUG
+            assert( size() == q.size() );
+#endif
+            for( size_t i = 0; i < size(); i++ )
+                if( !(_p[i] <= q[i]) )
+                    return false;
+            return true;
+        }
+

         /// Pointwise multiplication with q
         TProb<T>& operator*= (const TProb<T> & q) {
 #ifdef DAI_DEBUG
         /// Pointwise multiplication with q
         TProb<T>& operator*= (const TProb<T> & q) {
 #ifdef DAI_DEBUG


@@ -288,6 +318,30 @@ template <typename T> class TProb {
             return power;
         }
 
             return power;
         }
 

+        /// Pointwise signum
+        TProb<T> sgn() const {
+            TProb<T> x;
+            x._p.reserve( size() );
+            for( size_t i = 0; i < size(); i++ ) {
+                T s = 0;
+                if( _p[i] > 0 )
+                    s = 1;
+                else if( _p[i] < 0 )
+                    s = -1;
+                x._p.push_back( s );
+            }
+            return x;
+        }
+
+        /// Pointwise absolute value
+        TProb<T> abs() const {
+            TProb<T> x;
+            x._p.reserve( size() );
+            for( size_t i = 0; i < size(); i++ )
+                x._p.push_back( _p[i] < 0 ? (-p[i]) : p[i] );
+            return x;
+        }
+

         /// Pointwise exp
         const TProb<T>& takeExp() {
             std::transform( _p.begin(), _p.end(), _p.begin(),  std::ptr_fun<T, T>(std::exp) );
         /// Pointwise exp
         const TProb<T>& takeExp() {
             std::transform( _p.begin(), _p.end(), _p.begin(),  std::ptr_fun<T, T>(std::exp) );


@@ -396,8 +450,14 @@ template <typename T> class TProb {
             return Z;
         }
 
             return Z;
         }
 

+        /// Returns minimum value
+        T minVal() const {
+            T Z = *std::min_element( _p.begin(), _p.end() );
+            return Z;
+        }
+

         /// Normalize, using the specified norm
         /// Normalize, using the specified norm

-        T normalize( NormType norm ) {
+        T normalize( NormType norm = NORMPROB ) {

             T Z = 0.0;
             if( norm == NORMPROB )
                 Z = totalSum();
             T Z = 0.0;
             if( norm == NORMPROB )
                 Z = totalSum();


@@ -411,7 +471,7 @@ template <typename T> class TProb {
         }
 
         /// Return normalized copy of *this, using the specified norm
         }
 
         /// Return normalized copy of *this, using the specified norm

-        TProb<T> normalized( NormType norm ) const {
+        TProb<T> normalized( NormType norm = NORMPROB ) const {

             TProb<T> result(*this);
             result.normalize( norm );
             return result;
             TProb<T> result(*this);
             result.normalize( norm );
             return result;


@@ -440,14 +500,45 @@ template <typename T> class TProb {
             return S;
         }
 
             return S;
         }
 

+        /// Returns TProb<T> containing the pointwise minimum of a and b (which should have equal size)
+        friend TProb<T> min <> ( const TProb<T> &a, const TProb<T> &b );
+
+        /// Returns TProb<T> containing the pointwise maximum of a and b (which should have equal size)
+        friend TProb<T> max <> ( const TProb<T> &a, const TProb<T> &b );
+

         friend std::ostream& operator<< (std::ostream& os, const TProb<T>& P) {
         friend std::ostream& operator<< (std::ostream& os, const TProb<T>& P) {

+            os << "[";

             std::copy( P._p.begin(), P._p.end(), std::ostream_iterator<T>(os, " ") );
             std::copy( P._p.begin(), P._p.end(), std::ostream_iterator<T>(os, " ") );

-            os << std::endl;
+            os << "]";

             return os;
         }
 };
 
 
             return os;
         }
 };
 
 

+template<typename T> TProb<T> min( const TProb<T> &a, const TProb<T> &b ) {
+    assert( a.size() == b.size() );
+    TProb<T> result( a.size() );
+    for( size_t i = 0; i < a.size(); i++ )
+        if( a[i] < b[i] )
+            result[i] = a[i];
+        else
+            result[i] = b[i];
+    return result;
+}
+
+
+template<typename T> TProb<T> max( const TProb<T> &a, const TProb<T> &b ) {
+    assert( a.size() == b.size() );
+    TProb<T> result( a.size() );
+    for( size_t i = 0; i < a.size(); i++ )
+        if( a[i] > b[i] )
+            result[i] = a[i];
+        else
+            result[i] = b[i];
+    return result;
+}
+
+

 } // end of namespace dai
 
 
 } // end of namespace dai
 
 

diff --git a/include/dai/regiongraph.h b/include/dai/regiongraph.h
index 16916c0..6268d76 100644 (file)
--- a/include/dai/regiongraph.h
+++ b/include/dai/regiongraph.h
@@ -184,10 +184,10 @@ class RegionGraph : public FactorGraph {
         void makeCavity( size_t i ) { FactorGraph::makeCavity( i ); RecomputeORs( var(i) ); }
 
         /// We have to overload FactorGraph::undoProbs because the corresponding outer regions have to be recomputed
         void makeCavity( size_t i ) { FactorGraph::makeCavity( i ); RecomputeORs( var(i) ); }
 
         /// We have to overload FactorGraph::undoProbs because the corresponding outer regions have to be recomputed

-        void undoProbs( const VarSet &ns ) { FactorGraph::undoProbs( ns ); RecomputeORs( ns ); }
+        void undoProbs( const VarSet &ns ) { FactorGraph::restoreFactors( ns ); RecomputeORs( ns ); }

 
         /// We have to overload FactorGraph::undoProb because the corresponding outer regions have to be recomputed
 
         /// We have to overload FactorGraph::undoProb because the corresponding outer regions have to be recomputed

-        void undoProb( size_t I ) { FactorGraph::undoProb( I ); RecomputeOR( I ); }
+        void undoProb( size_t I ) { FactorGraph::restoreFactor( I ); RecomputeOR( I ); }

 
         /// Send RegionGraph to output stream
         friend std::ostream & operator << ( std::ostream & os, const RegionGraph & rg );
 
         /// Send RegionGraph to output stream
         friend std::ostream & operator << ( std::ostream & os, const RegionGraph & rg );

diff --git a/src/factorgraph.cpp b/src/factorgraph.cpp
index 4df27c8..1a25f8d 100644 (file)
--- a/src/factorgraph.cpp
+++ b/src/factorgraph.cpp
@@ -53,7 +53,7 @@ FactorGraph::FactorGraph( const std::vector<Factor> &P ) : G(), _undoProbs() {
     vars.reserve( _vars.size() );
     for( set<Var>::const_iterator p1 = _vars.begin(); p1 != _vars.end(); p1++ )
         vars.push_back( *p1 );
     vars.reserve( _vars.size() );
     for( set<Var>::const_iterator p1 = _vars.begin(); p1 != _vars.end(); p1++ )
         vars.push_back( *p1 );

-    
+

     // create graph structure
     createGraph( nrEdges );
 }
     // create graph structure
     createGraph( nrEdges );
 }


@@ -115,22 +115,23 @@ istream& operator >> (istream& is, FactorGraph& fg) {
 
     try {
         vector<Factor> factors;
 
     try {
         vector<Factor> factors;

-        size_t nr_f;
+        size_t nr_Factors;

         string line;
         
         while( (is.peek()) == '#' )
             getline(is,line);
         string line;
         
         while( (is.peek()) == '#' )
             getline(is,line);

-        is >> nr_f;
+        is >> nr_Factors;

         if( is.fail() )
             DAI_THROW(INVALID_FACTORGRAPH_FILE);
         if( verbose >= 2 )
         if( is.fail() )
             DAI_THROW(INVALID_FACTORGRAPH_FILE);
         if( verbose >= 2 )

-            cout << "Reading " << nr_f << " factors..." << endl;
+            cout << "Reading " << nr_Factors << " factors..." << endl;

 
         getline (is,line);
         if( is.fail() )
             DAI_THROW(INVALID_FACTORGRAPH_FILE);
 
 
         getline (is,line);
         if( is.fail() )
             DAI_THROW(INVALID_FACTORGRAPH_FILE);
 

-        for( size_t I = 0; I < nr_f; I++ ) {
+        map<long,size_t> vardims;
+        for( size_t I = 0; I < nr_Factors; I++ ) {

             if( verbose >= 3 )
                 cout << "Reading factor " << I << "..." << endl;
             size_t nr_members;
             if( verbose >= 3 )
                 cout << "Reading factor " << I << "..." << endl;
             size_t nr_members;


@@ -148,11 +149,8 @@ istream& operator >> (istream& is, FactorGraph& fg) {
                 is >> mi_label;
                 labels.push_back(mi_label);
             }
                 is >> mi_label;
                 labels.push_back(mi_label);
             }

-            if( verbose >= 3 ) {
-                cout << "  labels: ";
-                copy (labels.begin(), labels.end(), ostream_iterator<int>(cout, " "));
-                cout << endl;
-            }
+            if( verbose >= 3 )
+                cout << "  labels: " << labels << endl;

 
             vector<size_t> dims;
             for( size_t mi = 0; mi < nr_members; mi++ ) {
 
             vector<size_t> dims;
             for( size_t mi = 0; mi < nr_members; mi++ ) {


@@ -162,17 +160,22 @@ istream& operator >> (istream& is, FactorGraph& fg) {
                 is >> mi_dim;
                 dims.push_back(mi_dim);
             }
                 is >> mi_dim;
                 dims.push_back(mi_dim);
             }

-            if( verbose >= 3 ) {
-                cout << "  dimensions: ";
-                copy (dims.begin(), dims.end(), ostream_iterator<int>(cout, " "));
-                cout << endl;
-            }
+            if( verbose >= 3 )
+                cout << "  dimensions: " << dims << endl;

 
             // add the Factor
             VarSet I_vars;
 
             // add the Factor
             VarSet I_vars;

-            for( size_t mi = 0; mi < nr_members; mi++ )
+            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_THROW(INVALID_FACTORGRAPH_FILE);
+                } else
+                    vardims[labels[mi]] = dims[mi];

                 I_vars |= Var(labels[mi], dims[mi]);
                 I_vars |= Var(labels[mi], dims[mi]);

-            factors.push_back(Factor(I_vars,0.0));
+            }
+            factors.push_back( Factor( I_vars, 0.0 ) );

             
             // calculate permutation sigma (internally, members are sorted)
             vector<size_t> sigma(nr_members,0);
             
             // calculate permutation sigma (internally, members are sorted)
             vector<size_t> sigma(nr_members,0);


@@ -182,11 +185,8 @@ istream& operator >> (istream& is, FactorGraph& fg) {
                 vector<long>::iterator j_loc = find(labels.begin(),labels.end(),search_for);
                 sigma[mi] = j_loc - labels.begin();
             }
                 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;
-            }
+            if( verbose >= 3 )
+                cout << "  sigma: " << sigma << endl;

 
             // calculate multindices
             Permute permindex( dims, sigma );
 
             // calculate multindices
             Permute permindex( dims, sigma );


@@ -214,10 +214,8 @@ istream& operator >> (istream& is, FactorGraph& fg) {
             }
         }
 
             }
         }
 

-        if( verbose >= 3 ) {
-            cout << "factors:" << endl;
-            copy(factors.begin(), factors.end(), ostream_iterator<Factor>(cout,"\n"));
-        }
+        if( verbose >= 3 )
+            cout << "factors:" << factors << endl;

 
         fg = FactorGraph(factors);
     } catch (char *e) {
 
         fg = FactorGraph(factors);
     } catch (char *e) {


@@ -229,141 +227,72 @@ istream& operator >> (istream& is, FactorGraph& fg) {
 
 
 VarSet FactorGraph::delta( unsigned i ) const {
 
 
 VarSet FactorGraph::delta( unsigned i ) const {

+    return( Delta(i) / var(i) );
+}
+
+
+VarSet FactorGraph::Delta( unsigned i ) const {

     // calculate Markov Blanket
     // calculate Markov Blanket

-    VarSet del;
+    VarSet Del;

     foreach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
         foreach( const Neighbor &j, nbF(I) ) // for all neighboring variables j of I
     foreach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
         foreach( const Neighbor &j, nbF(I) ) // for all neighboring variables j of I

-            if( j != i )
-                del |= var(j);
+            Del |= var(j);

 
 

-    return del;
+    return Del;

 }
 
 
 }
 
 

-VarSet FactorGraph::Delta( unsigned i ) const {
-    return( delta(i) | var(i) );
+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( unsigned i ) {
+void FactorGraph::makeCavity( unsigned i, bool backup ) {

     // fills all Factors that include var(i) with ones
     // fills all Factors that include var(i) with ones

+    map<size_t,Factor> newFacs;

     foreach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
     foreach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i

-        factor(I).fill( 1.0 );
+        newFacs[I] = Factor(factor(I).vars(), 1.0);
+    setFactors( newFacs, backup );

 }
 
 
 }
 
 

-bool FactorGraph::hasNegatives() const {
-    bool result = false;
-    for( size_t I = 0; I < nrFactors() && !result; I++ )
-        if( factor(I).hasNegatives() )
-            result = true;
-    return result;
-}
- 
-
-long FactorGraph::ReadFromFile(const char *filename) {
+void FactorGraph::ReadFromFile( const char *filename ) {

     ifstream infile;
     ifstream infile;

-    infile.open (filename);
-    if (infile.is_open()) {
+    infile.open( filename );
+    if( infile.is_open() ) {

         infile >> *this;
         infile.close();
         infile >> *this;
         infile.close();

-        return 0;
-    } else {
-        cout << "ERROR OPENING FILE" << endl;
-        return 1;
-    }
+    } else
+        DAI_THROW(CANNOT_READ_FILE);

 }
 
 
 }
 
 

-long FactorGraph::WriteToFile(const char *filename) const {
+void FactorGraph::WriteToFile( const char *filename ) const {

     ofstream outfile;
     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 << *this;

         outfile.close();
         outfile.close();

-        return 0;
-    } else {
-        cout << "ERROR OPENING FILE" << endl;
-        return 1;
-    }
-}
-
-
-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 i = 0; i < nrVars(); i++ )
-                foreach( const Neighbor &I, nbV(i) )  // for all neighboring factors I of i
-                    outfile << "\tx" << var(i).label() << " -- p" << I << ";" << endl;
-            outfile << "}" << endl;
-        } catch (char *e) {
-            cout << e << endl;
-            return 1;
-        }
-        outfile.close();
-        return 0;
-    } else {
-        cout << "ERROR OPENING FILE" << endl;
-        return 1;
-    }
-}
-
-
-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;
-            }
-        if( found )
-            break;
-    }
-    return found;
+    } else
+        DAI_THROW(CANNOT_WRITE_FILE);

 }
 
 
 }
 
 

-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);
-        }
-    }
+void FactorGraph::display( ostream &os ) const {
+    os << "graph G {" << 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++ )
+        foreach( const Neighbor &I, nbV(i) )  // for all neighboring factors I of i
+            os << "\tv" << var(i).label() << " -- f" << I << ";" << endl;
+    os << "}" << endl;

 }
 
 
 }
 
 


@@ -373,7 +302,7 @@ vector<VarSet> FactorGraph::Cliques() const {
     for( size_t I = 0; I < nrFactors(); I++ ) {
         bool maximal = true;
         for( size_t J = 0; (J < nrFactors()) && maximal; J++ )
     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()) )
+            if( (factor(J).vars() >> factor(I).vars()) && (factor(J).vars() != factor(I).vars()) )

                 maximal = false;
         
         if( maximal )
                 maximal = false;
         
         if( maximal )


@@ -384,7 +313,7 @@ vector<VarSet> FactorGraph::Cliques() const {
 }
 
 
 }
 
 

-void FactorGraph::clamp( const Var & n, size_t i ) {
+void FactorGraph::clamp( const Var & n, size_t i, bool backup ) {

     assert( i <= n.states() );
 
     // Multiply each factor that contains the variable with a delta function
     assert( i <= n.states() );
 
     // Multiply each factor that contains the variable with a delta function


@@ -392,52 +321,137 @@ void FactorGraph::clamp( const Var & n, size_t i ) {
     Factor delta_n_i(n,0.0);
     delta_n_i[i] = 1.0;
 
     Factor delta_n_i(n,0.0);
     delta_n_i[i] = 1.0;
 

+    map<size_t, Factor> newFacs;

     // For all factors that contain n
     for( size_t I = 0; I < nrFactors(); I++ ) 
         if( factor(I).vars().contains( n ) )
             // Multiply it with a delta function
     // For all factors that contain n
     for( size_t I = 0; I < nrFactors(); I++ ) 
         if( factor(I).vars().contains( n ) )
             // Multiply it with a delta function

-            factor(I) *= delta_n_i;
+            newFacs[I] = factor(I) * delta_n_i;
+    setFactors( newFacs, backup );

 
     return;
 }
 
 
 
     return;
 }
 
 

-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 = _backupFactors.find( I );
+    if( it != _backupFactors.end() )
+        DAI_THROW( MULTIPLE_UNDO );
+    _backupFactors[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 = _backupFactors.find( I );
+    if( it != _backupFactors.end() ) {
+        setFactor(I, it->second);
+        _backupFactors.erase(it);

     }
 }
 
 
     }
 }
 
 

-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().intersects( ns ) )
     for( size_t I = 0; I < nrFactors(); I++ )
         if( factor(I).vars().intersects( ns ) )

-            _undoProbs[I] = factor(I).p();
+            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().intersects( 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 = _backupFactors.begin(); uI != _backupFactors.end(); ) {
+        if( factor(uI->first).vars().intersects( ns ) ) {
+            facs.insert( *uI );
+            _backupFactors.erase(uI++);

         } else
             uI++;
     }
         } else
             uI++;
     }

+    setFactors( facs );
+}
+
+
+void FactorGraph::restoreFactors() {
+    setFactors( _backupFactors );
+    _backupFactors.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::isPairwise() const {
+    bool pairwise = true;
+    for( size_t I = 0; I < nrFactors() && pairwise; I++ )
+        if( factor(I).vars().size() > 2 )
+            pairwise = false;
+    return pairwise;
+}
+
+
+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;
+}
+
+
+FactorGraph FactorGraph::clamped( const Var & v_i, size_t state ) const {
+    Real zeroth_order = 1.0;
+    vector<Factor> clamped_facs;
+    for( size_t I = 0; I < nrFactors(); I++ ) {
+        VarSet v_I = factor(I).vars();
+        Factor new_factor;
+        if( v_I.intersects( v_i ) )
+            new_factor = factor(I).slice( v_i, 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 );
+}
+
+
+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() );

 }
 
 
 }