Improved ClusterGraph implementation and MaxSpanningTreePrims implementation.
authorJoris Mooij <jorism@osun.tuebingen.mpg.de>
Thu, 11 Sep 2008 07:36:23 +0000 (09:36 +0200)
committerJoris Mooij <jorism@osun.tuebingen.mpg.de>
Thu, 11 Sep 2008 07:36:23 +0000 (09:36 +0200)
- Added several methods to BipartiteGraph
- Improved ClusterGraph implementation (taken from SVN head)
- Improved MaxSpanningTree algorithm (using boost::graph library, taken from SVN head)

12 files changed:
include/dai/bipgraph.h
include/dai/clustergraph.h
include/dai/factorgraph.h
include/dai/varset.h
include/dai/weightedgraph.h
src/bp.cpp
src/clustergraph.cpp
src/factorgraph.cpp
src/jtree.cpp
src/lc.cpp
src/regiongraph.cpp
src/treeep.cpp

index ef8e4cb..6c4f3f4 100644 (file)
@@ -23,8 +23,9 @@
 #define __defined_libdai_bipgraph_h
 
 
+#include <ostream>
 #include <vector>
-#include <algorithm>
+#include <cassert>
 #include <dai/util.h>
 
 
@@ -62,12 +63,21 @@ class BipartiteGraph {
         /// Neighbors is a vector of Neigbor entries
         typedef std::vector<Neighbor> Neighbors;
 
+        /// Edge is used as index of an edge
+        typedef std::pair<size_t,size_t> Edge;
+
     private:
         /// _nb1 contains for each node of the first kind a list of its neighbors
         std::vector<Neighbors> _nb1;
         /// _nb2 contains for each node of the second kind a list of its neighbors
         std::vector<Neighbors> _nb2;
 
+        /// Used internally by isTree()
+        struct levelType {
+            std::vector<size_t> ind1;       // indices of vertices of type 1
+            std::vector<size_t> ind2;       // indices of vertices of type 2
+        };
+
     public:
         /// Default constructor
         BipartiteGraph() : _nb1(), _nb2() {}
@@ -88,7 +98,7 @@ class BipartiteGraph {
         /// (more precisely, a std::pair<unsigned, unsigned> where the first integer corresponds
         /// to the node of the first type, and the second integer corresponds to the node of the
         /// second type). nr1 is the number of nodes of the first type, nr2 the number of nodes
-        /// of the second type.
+        /// of the second type. The values between the iterators begin and end should be of type Edge.
         template<typename EdgeInputIterator>
         void create( size_t nr1, size_t nr2, EdgeInputIterator begin, EdgeInputIterator end );
 
@@ -96,7 +106,7 @@ class BipartiteGraph {
         /// (more precisely, a std::pair<unsigned, unsigned> where the first integer corresponds
         /// to the node of the first type, and the second integer corresponds to the node of the
         /// second type). nr1 is the number of nodes of the first type, nr2 the number of nodes
-        /// of the second type.
+        /// of the second type. The values between the iterators begin and end should be of type Edge.
         template<typename EdgeInputIterator>
         BipartiteGraph( size_t nr1, size_t nr2, EdgeInputIterator begin, EdgeInputIterator end ) : _nb1( nr1 ), _nb2( nr2 ) {
             create( nr1, nr2, begin, end );
@@ -134,9 +144,150 @@ class BipartiteGraph {
                 sum += nb1(i1).size();
             return sum;
         }
+        
+        /// Add node of first type without neighbors.
+        void add1() {
+            _nb1.push_back( Neighbors() );
+        }
+        
+        /// Add node of first type without neighbors.
+        void add2() {
+            _nb2.push_back( Neighbors() );
+        }
+
+        /// Add node of first type with neighbors specified by a range.
+        /// For improved efficiency, the size of the range may be specified by sizeHint.
+        /// *NodeInputIterator should be a size_t.
+        template <typename NodeInputIterator>
+        void add1( NodeInputIterator begin, NodeInputIterator end, size_t sizeHint = 0 ) {
+            Neighbors nbs1new;
+            nbs1new.reserve( sizeHint );
+            size_t iter = 0;
+            for( NodeInputIterator it = begin; it != end; ++it ) {
+                assert( *it < nr2() );
+
+                Neighbor nb1new;
+                nb1new.node = *it;
+                nb1new.iter = iter;
+                nb1new.dual = nb2(*it).size();
+
+                Neighbor nb2new;
+                nb2new.node = nr1();
+                nb2new.iter = nb2(*it).size();
+                nb2new.dual = iter++;
+
+                nbs1new.push_back( nb1new );
+                nb2( *it ).push_back( nb2new );
+            }
+            _nb1.push_back( nbs1new );
+        }
+
+        /// Add node of second type with neighbors specified by a range.
+        /// For improved efficiency, the size of the range may be specified by sizeHint.
+        /// *NodeInputIterator should be a size_t.
+        template <typename NodeInputIterator>
+        void add2( NodeInputIterator begin, NodeInputIterator end, size_t sizeHint = 0 ) {
+            Neighbors nbs2new;
+            nbs2new.reserve( sizeHint );
+            size_t iter = 0;
+            for( NodeInputIterator it = begin; it != end; ++it ) {
+                assert( *it < nr1() );
+
+                Neighbor nb2new;
+                nb2new.node = *it;
+                nb2new.iter = iter;
+                nb2new.dual = nb1(*it).size();
+
+                Neighbor nb1new;
+                nb1new.node = nr2();
+                nb1new.iter = nb1(*it).size();
+                nb1new.dual = iter++;
+
+                nbs2new.push_back( nb2new );
+                nb1( *it ).push_back( nb1new );
+            }
+            _nb2.push_back( nbs2new );
+        }
+
+        /// Remove node of first type
+        void erase1( size_t n1 ) {
+            assert( n1 < nr1() );
+            // Erase neighbor entry of node n1
+            _nb1.erase( _nb1.begin() + n1 );
+            // Adjust neighbor entries of nodes of second type
+            for( size_t n2 = 0; n2 < nr2(); n2++ )
+                for( size_t iter = 0; iter < nb2(n2).size(); ) {
+                    if( nb2(n2, iter).node == n1 ) {
+                        // delete this entry, because it points to the deleted node
+                        nb2(n2).erase( nb2(n2).begin() + iter );
+                        // adjust all subsequent entries:
+                        // update their iter and the corresponding dual of the neighboring node of the first kind
+                        for( size_t newiter = iter; newiter < nb2(n2).size(); newiter++ ) {
+                            nb2( n2, newiter ).iter = newiter;
+                            nb1( nb2(n2, newiter).node, nb2(n2, newiter).dual ).dual = newiter;
+                        }
+                    } else if( nb2(n2, iter).node > n1 ) {
+                        nb2(n2, iter).node--;
+                        iter++;
+                    } else
+                        iter++;
+                }
+        }
+
+        /// Remove node of second type
+        void erase2( size_t n2 ) {
+            assert( n2 < nr2() );
+            // Erase neighbor entry of node n2
+            _nb2.erase( _nb2.begin() + n2 );
+            // Adjust neighbor entries of nodes of first type
+            for( size_t n1 = 0; n1 < nr1(); n1++ )
+                for( size_t iter = 0; iter < nb1(n1).size(); ) {
+                    if( nb1(n1, iter).node == n2 ) {
+                        // delete this entry, because it points to the deleted node
+                        nb1(n1).erase( nb1(n1).begin() + iter );
+                        // adjust all subsequent entries:
+                        // update their iter and the corresponding dual of the neighboring node of the first kind
+                        for( size_t newiter = iter; newiter < nb1(n1).size(); newiter++ ) {
+                            nb1( n1, newiter ).iter = newiter;
+                            nb2( nb1(n1, newiter).node, nb1(n1, newiter).dual ).dual = newiter;
+                        }
+                    } else if( nb1(n1, iter).node > n2 ) {
+                        nb1(n1, iter).node--;
+                        iter++;
+                    } else
+                        iter++;
+                }
+        }
+
+        /// Calculate neighbors of neighbors of node n1 of first type.
+        /// If include == true, include n1 itself, otherwise exclude n1.
+        std::vector<size_t> delta1( size_t n1, bool include = false ) const {
+            std::vector<size_t> result;
+            foreach( const Neighbor &n2, nb1(n1) )
+                foreach( const Neighbor &m1, nb2(n2) )
+                    if( include || (m1 != n1) )
+                        result.push_back( m1 );
+            // remove duplicates
+            std::vector<size_t>::iterator it = unique( result.begin(), result.end() );
+            result.erase( it, result.end() );
+            return result;
+        }
+
+        /// Calculate neighbors of neighbors of node n2 of second type
+        /// If include == true, include n2 itself, otherwise exclude n2.
+        std::vector<size_t> delta2( size_t n2, bool include = false ) const {
+            std::vector<size_t> result;
+            foreach( const Neighbor &n1, nb2(n2) )
+                foreach( const Neighbor &m2, nb1(n1) )
+                    if( include || (m2 != n2) )
+                        result.push_back( m2 );
+            // remove duplicates
+            std::vector<size_t>::iterator it = unique( result.begin(), result.end() );
+            result.erase( it, result.end() );
+            return result;
+        }
 
         /// Returns true if the graph is connected
-        /// FIXME: this should be optimized
         bool isConnected() const {
             if( nr1() == 0 ) {
                 return true;
@@ -187,6 +338,91 @@ class BipartiteGraph {
             }
         }
 
+        /// Returns true if the graph is a tree, i.e., if it is singly connected and connected.
+        /// This is equivalent to whether for each pair of vertices in the graph, there exists
+        /// a unique path in the graph that starts at the first and ends at the second vertex.
+        bool isTree() const {
+            using namespace std;
+            vector<levelType> levels;
+
+            bool foundCycle = false;
+            size_t nr_1 = 0;
+            size_t nr_2 = 0;
+
+            if( nr1() == 0 || nr2() == 0 )
+                return true;
+            else {
+                levelType newLevel;
+                do {
+                    newLevel.ind1.clear();
+                    newLevel.ind2.clear();
+                    if( levels.size() == 0 ) {
+                        size_t n1 = 0;
+                        // add n1 to ind1
+                        newLevel.ind1 = vector<size_t>( 1, n1 );
+                        // add all neighbors of n1 to ind2
+                        newLevel.ind2.reserve( nb1(n1).size() );
+                        foreach( const Neighbor &n2, nb1(n1) )
+                            newLevel.ind2.push_back( n2 );
+                    } else {
+                        const levelType &prevLevel = levels.back();
+                        // build newLevel.ind1
+                        foreach( size_t n2, prevLevel.ind2 ) { // for all n2 in the previous level
+                            foreach( const Neighbor &n1, nb2(n2) ) { // for all neighbors n1 of n2
+                                if( find( prevLevel.ind1.begin(), prevLevel.ind1.end(), n1 ) == prevLevel.ind1.end() ) { // n1 not in previous level
+                                    if( find( newLevel.ind1.begin(), newLevel.ind1.end(), n1 ) != newLevel.ind1.end() )
+                                        foundCycle = true; // n1 already in new level: we found a cycle
+                                    else
+                                        newLevel.ind1.push_back( n1 ); // add n1 to new level
+                                }
+                                if( foundCycle )
+                                    break;
+                            }
+                            if( foundCycle )
+                                break;
+                        }
+                        // build newLevel.ind2
+                        foreach( size_t n1, newLevel.ind1 ) { // for all n1 in this level
+                            foreach( const Neighbor &n2, nb1(n1) ) { // for all neighbors n2 of n1
+                                if( find( prevLevel.ind2.begin(), prevLevel.ind2.end(), n2 ) == prevLevel.ind2.end() ) { // n2 not in previous level
+                                    if( find( newLevel.ind2.begin(), newLevel.ind2.end(), n2 ) != newLevel.ind2.end() )
+                                        foundCycle = true; // n2 already in new level: we found a cycle
+                                    else
+                                        newLevel.ind2.push_back( n2 ); // add n2 to new level
+                                }
+                                if( foundCycle )
+                                    break;
+                            }
+                            if( foundCycle )
+                                break;
+                        }
+                    }
+                    levels.push_back( newLevel );
+                    nr_1 += newLevel.ind1.size();
+                    nr_2 += newLevel.ind2.size();
+                } while( ((newLevel.ind1.size() != 0) || (newLevel.ind2.size() != 0)) && !foundCycle );
+                if( nr_1 == nr1() && nr_2 == nr2() && !foundCycle )
+                    return true;
+                else
+                    return false;
+            }
+        }
+
+        /// Stream to output stream os in graphviz .dot syntax
+        void display( std::ostream& os ) const {
+            using namespace std;
+            os << "graph G {" << endl;
+            os << "node[shape=circle,width=0.4,fixedsize=true];" << endl;
+            for( size_t n1 = 0; n1 < nr1(); n1++ )
+                os << "\tx" << n1 << ";" << endl;
+            os << "node[shape=box,width=0.3,height=0.3,fixedsize=true];" << endl;
+            for( size_t n2 = 0; n2 < nr2(); n2++ )
+                os << "\ty" << n2 << ";" << endl;
+            for( size_t n1 = 0; n1 < nr1(); n1++ )
+                foreach( const Neighbor &n2, nb1(n1) )
+                    os << "\tx" << n1 << " -- y" << n2 << ";" << endl;
+            os << "}" << endl;
+        }
 };
 
 
index 1de3e0d..9d2772e 100644 (file)
 #include <set>
 #include <vector>
 #include <dai/varset.h>
+#include <dai/bipgraph.h>
 
 
 namespace dai {
 
 
-/// A ClusterGraph is a hypergraph with VarSets as nodes.
-/// It is implemented as a set<VarSet> in which the adjacency
-/// relationship is computed realtime. It may be better to
-/// implement it as a bipartitegraph, though. The additional
-/// functionality compared to a simple set<VarSet> is
-/// finding maximal clusters, finding cliques, etc...
-class ClusterGraph : public std::set<VarSet> {
-    public:
-        /// Default constructor
-        ClusterGraph() : std::set<VarSet>() {}
-
-        /// Construct from vector<VarSet>
-        ClusterGraph( const std::vector<VarSet> & cls ) {
-            insert( cls.begin(), cls.end() );
-        }
-        
-        /// Copy constructor
-        ClusterGraph(const ClusterGraph &x) : std::set<VarSet>(x) {}
-
-        /// Assignment operator
-        ClusterGraph& operator=( const ClusterGraph &x ) {
-            if( this != &x ) {
-                std::set<VarSet>::operator=( x );
+    /// A ClusterGraph is a hypergraph with VarSets as nodes.
+    /// It is implemented as bipartite graph with variable nodes
+    /// and cluster (VarSet) nodes. The additional
+    /// functionality compared to a simple set<VarSet> is
+    /// finding maximal clusters, finding cliques, etc...
+    class ClusterGraph {
+        public:
+            BipartiteGraph       G;
+            std::vector<Var>     vars;
+            std::vector<VarSet>  clusters;
+
+            typedef BipartiteGraph::Neighbor Neighbor;
+            typedef BipartiteGraph::Edge     Edge;
+
+        public:
+            /// Default constructor
+            ClusterGraph() : G(), vars(), clusters() {}
+
+            /// Construct from vector<VarSet>
+            ClusterGraph( const std::vector<VarSet> & cls );
+            
+            /// Copy constructor
+            ClusterGraph( const ClusterGraph &x ) : G(x.G), vars(x.vars), clusters(x.clusters) {}
+
+            /// Assignment operator
+            ClusterGraph& operator=( const ClusterGraph &x ) {
+                if( this != &x ) {
+                    G = x.G;
+                    vars = x.vars;
+                    clusters = x.clusters;
+                }
+                return *this;
             }
-            return *this;
-        }
 
-        /// Returns true if ns is a maximal member of *this under inclusion (VarSet::operator<<)
-        bool isMaximal( const VarSet &ns ) const {
-            if( count( ns ) ) {
-                // ns is a member
+            /// Returns true if cluster I is not contained in a larger cluster
+            bool isMaximal( size_t I ) const {
+#ifdef DAI_DEBUG
+                assert( I < G.nr2() );
+#endif
+                const VarSet & clI = clusters[I];
                 bool maximal = true;
-                for( const_iterator x = begin(); x != end() && maximal; x++ )
-                    if( (ns << *x) && (ns != *x) )
-                        maximal = false;
+                // The following may not be optimal, since it may repeatedly test the same cluster *J
+                foreach( const Neighbor &i, G.nb2(I) ) {
+                    foreach( const Neighbor &J, G.nb1(i) )
+                        if( (J != I) && (clI << clusters[J]) ) {
+                            maximal = false;
+                            break;
+                        }
+                    if( !maximal )
+                        break;
+                }
                 return maximal;
-            } else
-                return false;
-        }
-
-        /// Erase all VarSets that are not maximal
-        ClusterGraph& eraseNonMaximal() {
-            for( iterator x = begin(); x != end(); )
-                if( !isMaximal(*x) )
-                    erase(x++);
-                else
-                    x++;
-            return *this;
-        }
-
-        /// Return union of all members
-        VarSet vars() const {
-            VarSet result;
-            for( const_iterator x = begin(); x != end(); x++ )
-                result |= *x;
-            return result;
-        }
-
-        /// Returns true if n1 and n2 are adjacent, i.e.\ by
-        /// definition, are both contained in some cluster in *this
-        bool adj( const Var& n1, const Var& n2 ) {
-            bool result = false;
-            for( iterator x = begin(); (x != end()) && (!result); x++ )
-                if( (*x && n1) && (*x && n2) )
-                    result = true;
-            return result;
-        }
-        
-        /// Returns union of clusters that contain n, minus n
-        VarSet Delta( const Var& n ) const {
-            VarSet result;
-            for( const_iterator x = begin(); x != end(); x++ )
-                if( (*x && n) )
-                    result |= *x;
-            return result;
-        }
-
-        /// Returns union of clusters that contain n, minus n
-        VarSet delta( const Var& n ) const {
-            return Delta( n ) / n;
-        }
-        
-        /// Erases all members that contain n
-        ClusterGraph& eraseSubsuming( const Var& n ) {
-            for( iterator x = begin(); x != end(); )
-                if( (*x && n) )
-                    erase(x++);
-                else
-                    x++;
-            return *this;
-        }
-        
-        /// Send to output stream
-        friend std::ostream & operator << ( std::ostream & os, const ClusterGraph & cl ) {
-            os << "{";
-            ClusterGraph::const_iterator x = cl.begin();
-            if( x != cl.end() )
-                os << *(x++);
-            for( ; x != cl.end(); x++ )
-                os << ", " << *x;
-            os << "}";
-            return os;
-        }
-
-        /// Convert to vector<VarSet>
-        std::vector<VarSet> toVector() const {
-            std::vector<VarSet> result;
-            result.reserve( size() );
-            for( const_iterator x = begin(); x != end(); x++ )
-                result.push_back( *x );
-            return result;
-        }
-
-        /// Calculate cost of eliminating variable n
-        /// using as a measure "number of added edges in the adjacency graph"
-        /// where the adjacency graph has the variables as its nodes and
-        /// connects nodes i1 and i2 iff i1 and i2 occur in some common factor I
-        size_t eliminationCost( const Var& n ) {
-            VarSet d_n = delta( n );
-            size_t cost = 0;
-
-            // for each unordered pair {i1,i2} adjacent to n
-            for( VarSet::const_iterator i1 = d_n.begin(); i1 != d_n.end(); i1++ ) {
-                VarSet d_i1 = delta( *i1 );
-                for( VarSet::const_iterator i2 = i1; (++i2) != d_n.end(); ) {
-                    // if i1 and i2 are not adjacent, eliminating n would make them adjacent
-                    if( !adj(*i1, *i2) )
-                        cost++;
+            }
+
+            /// Erase all VarSets that are not maximal
+            ClusterGraph& eraseNonMaximal() {
+                for( size_t I = 0; I < G.nr2(); ) {
+                    if( !isMaximal(I) ) {
+                        clusters.erase( clusters.begin() + I );
+                        G.erase2(I);
+                    } else
+                        I++;
+                }
+                return *this;
+            }
+
+            /// Return number of clusters
+            size_t size() const {
+                return G.nr2();
+            }
+
+            /// Return index of variable n
+            size_t findVar( const Var &n ) const {
+                return find( vars.begin(), vars.end(), n ) - vars.begin();
+            }
+
+            /// Returns true if vars with indices i1 and i2 are adjacent, i.e., both contained in the same cluster
+            bool adj( size_t i1, size_t i2 ) {
+                bool result = false;
+                foreach( const Neighbor &I, G.nb1(i1) )
+                    if( find( G.nb2(I).begin(), G.nb2(I).end(), i2 ) != G.nb2(I).end() ) {
+                        result = true;
+                        break;
+                    }
+                return result;
+            }
+            
+            /// Returns union of clusters that contain the variable with index i
+            VarSet Delta( size_t i ) const {
+                VarSet result;
+                foreach( const Neighbor &I, G.nb1(i) )
+                    result |= clusters[I];
+                return result;
+            }
+
+            /// Inserts a cluster (if it does not already exist)
+            void insert( const VarSet &cl ) {
+                if( find( clusters.begin(), clusters.end(), cl ) == clusters.end() ) {
+                    clusters.push_back( cl );
+                    // add variables (if necessary) and calculate neighborhood of new cluster
+                    std::vector<size_t> nbs;
+                    for( VarSet::const_iterator n = cl.begin(); n != cl.end(); n++ ) {
+                        size_t iter = find( vars.begin(), vars.end(), *n ) - vars.begin();
+                        nbs.push_back( iter );
+                        if( iter == vars.size() ) {
+                            G.add1();
+                            vars.push_back( *n );
+                        }
+                    }
+                    G.add2( nbs.begin(), nbs.end(), nbs.size() );
+                }
+            }
+
+            /// Returns union of clusters that contain variable with index i, minus this variable
+            VarSet delta( size_t i ) const {
+                return Delta( i ) / vars[i];
+            }
+            
+            /// Erases all clusters that contain n where n is the variable with index i
+            ClusterGraph& eraseSubsuming( size_t i ) {
+                while( G.nb1(i).size() ) {
+                    clusters.erase( clusters.begin() + G.nb1(i)[0] );
+                    G.erase2( G.nb1(i)[0] );
                 }
+                return *this;
+            }
+            
+            /// Provide read access to clusters
+            const std::vector<VarSet> & toVector() const { return clusters; }
+
+            /// Calculate cost of eliminating the variable with index i,
+            /// using as a measure "number of added edges in the adjacency graph"
+            /// where the adjacency graph has the variables as its nodes and
+            /// connects nodes i1 and i2 iff i1 and i2 occur in some common cluster
+            size_t eliminationCost( size_t i ) {
+                std::vector<size_t> id_n = G.delta1( i );
+
+                size_t cost = 0;
+
+                // for each unordered pair {i1,i2} adjacent to n
+                for( size_t _i1 = 0; _i1 < id_n.size(); _i1++ )
+                    for( size_t _i2 = _i1 + 1; _i2 < id_n.size(); _i2++ ) {
+                        // if i1 and i2 are not adjacent, eliminating n would make them adjacent
+                        if( !adj(id_n[_i1], id_n[_i2]) )
+                            cost++;
+                    }
+
+                return cost;
+            }
+
+            /// Perform Variable Elimination without Probs, i.e. only keeping track of
+            /// the interactions that are created along the way.
+            /// Input:  a set of outer clusters and an elimination sequence
+            /// Output: a set of elimination "cliques"
+            ClusterGraph VarElim( const std::vector<Var> &ElimSeq ) const;
+
+            /// Perform Variable Eliminiation using the MinFill heuristic
+            ClusterGraph VarElim_MinFill() const;
+
+            /// Send to output stream
+            friend std::ostream & operator << ( std::ostream & os, const ClusterGraph & cl ) {
+                os << cl.toVector();
+                return os;
             }
-            return cost;
-        }
-
-        /// Perform Variable Elimination without Probs, i.e. only keeping track of
-        /// the interactions that are created along the way.
-        /// Input:  a set of outer clusters and an elimination sequence
-        /// Output: a set of elimination "cliques"
-        ClusterGraph VarElim( const std::vector<Var> &ElimSeq ) const;
-
-        /// As Taylan does it
-        ClusterGraph VarElim_MinFill() const;
-};
+    };
 
 
 } // end of namespace dai
index adb6722..77f5f8a 100644 (file)
@@ -44,6 +44,7 @@ class FactorGraph {
         std::vector<Factor>    factors;
         typedef BipartiteGraph::Neighbor  Neighbor;
         typedef BipartiteGraph::Neighbors Neighbors;
+        typedef BipartiteGraph::Edge      Edge;
 
     protected:
         std::map<size_t,Prob>  _undoProbs;
index 9450637..606c3f8 100644 (file)
@@ -78,6 +78,13 @@ class VarSet : private std::set<Var> {
             calcStateSpace();
         }
 
+        /// Construct from a vector<Var>
+        VarSet( const std::vector<Var> &ns ) {
+            for( std::vector<Var>::const_iterator n = ns.begin(); n != ns.end(); n++ )
+                insert( *n );
+            calcStateSpace();
+        }
+
         /// Copy constructor
         VarSet( const VarSet &x ) : std::set<Var>( x ), _statespace( x._statespace ) {}
 
index 990d409..9423059 100644 (file)
 #include <map>
 #include <iostream>
 #include <set>
+#include <cassert>
+#include <limits>
+
+#include <boost/graph/adjacency_list.hpp>
+#include <boost/graph/prim_minimum_spanning_tree.hpp>
+
 
 
 namespace dai {
@@ -92,71 +98,87 @@ template<class T> class WeightedGraph : public std::map<UEdge, T> {};
 typedef std::set<UEdge>     Graph;
 
 
-/// Use Prim's algorithm to construct a maximal spanning tree from the weighted graph Graph
-template<typename T> DEdgeVec MaxSpanningTreePrim( const WeightedGraph<T> & Graph ) {
-    const long verbose = 0;
-
+/// Use Prim's algorithm to construct a minimal spanning tree from the weighted graph Graph
+/// The weights should be positive. Use implementation in Boost Graph Library.
+template<typename T> DEdgeVec MinSpanningTreePrims( const WeightedGraph<T> & Graph ) {
     DEdgeVec result;
-    if( Graph.size() == 0 )
-        return result;
-    else {
-        // Make a copy
-        WeightedGraph<T> Gr = Graph;
-
-        // Nodes in the tree
-        std::set<size_t> treeV;
-
-        // Start with one node
-        treeV.insert( Gr.begin()->first.n1 );
-        
-        // Perform Prim's algorithm
-        while( Gr.size() ) {
-            typename WeightedGraph<T>::iterator largest = Gr.end();
-            
-            for( typename WeightedGraph<T>::iterator e = Gr.begin(); e != Gr.end(); ) {
-                if( verbose >= 1 )
-                    std::cout << "considering edge " << e->first << "...";
-                bool e1_in_treeV = treeV.count( e->first.n1 );
-                bool e2_in_treeV = treeV.count( e->first.n2 );
-                if( e1_in_treeV && e2_in_treeV ) {
-                    if( verbose >= 1 )
-                        std::cout << "red";
-                    Gr.erase( e++ );    // Nice trick! 
-                } else if( e1_in_treeV || e2_in_treeV ) {
-                    if( verbose >= 1 )
-                        std::cout << e->second;
-                    if( (largest == Gr.end()) || (e->second > largest->second) ) {
-                        largest = e;    // largest edge connected to the tree (until now)
-                        if( verbose >= 1 )
-                            std::cout << " and largest!";
-                    } 
-                    e++;
-                } else {
-                    if( verbose >= 1 )
-                        std::cout << "out of reach";
-                    e++;
-                }
-                if( verbose >= 1 )
-                    std::cout << std::endl;
-            }
+    if( Graph.size() > 0 ) {
+        using namespace boost;
+        using namespace std;
+        typedef adjacency_list< vecS, vecS, undirectedS, property<vertex_distance_t, int>, property<edge_weight_t, double> > boostGraph;
+        typedef pair<size_t, size_t> E;
+
+        set<size_t> nodes;
+        vector<E> edges;
+        vector<double> weights;
+        edges.reserve( Graph.size() );
+        weights.reserve( Graph.size() );
+        for( typename WeightedGraph<T>::const_iterator e = Graph.begin(); e != Graph.end(); e++ ) {
+            weights.push_back( e->second );
+            edges.push_back( E( e->first.n1, e->first.n2 ) );
+            nodes.insert( e->first.n1 );
+            nodes.insert( e->first.n2 );
+        }
+
+        boostGraph g( edges.begin(), edges.end(), weights.begin(), nodes.size() );
+        vector< graph_traits< boostGraph >::vertex_descriptor > p( num_vertices(g) );
+        prim_minimum_spanning_tree( g, &(p[0]) );
 
-            if( largest != Gr.end() ) {
-                if( verbose >= 1 )
-                    std::cout << "largest = " << largest->first << std::endl;
-                // Add directed edge, pointing away from the root
-                if( treeV.count( largest->first.n1 ) ) {
-                    result.push_back( DEdge( largest->first.n1, largest->first.n2 ) );
-                    treeV.insert( largest->first.n2 );
-                } else {
-                    result.push_back( DEdge( largest->first.n2, largest->first.n1 ) );
-                    treeV.insert( largest->first.n1 );
+        // Store tree edges in result
+        result.reserve( nodes.size() - 1 );
+        size_t root = 0;
+        for( size_t i = 0; i != p.size(); i++ )
+            if( p[i] != i )
+                result.push_back( DEdge( p[i], i ) );
+            else
+                root = i;
+
+        // We have to store the minimum spanning tree in the right
+        // order, such that for all (i1, j1), (i2, j2) in result,
+        // if j1 == i2 then (i1, j1) comes before (i2, j2) in result.
+        // We do this by reordering the contents of result, effectively
+        // growing the tree starting at the root. At each step, 
+        // result[0..N-1] are the edges already added to the tree,
+        // whereas the other elements of result still have to be added.
+        // The elements of nodes are the vertices that still have to
+        // be added to the tree.
+
+        // Start with the root
+        nodes.erase( root );
+        size_t N = 0;
+
+        // Iteratively add edges and nodes to the growing tree
+        while( N != result.size() ) {
+            for( size_t e = N; e != result.size(); e++ ) {
+                bool e1_in_tree = !nodes.count( result[e].n1 );
+                if( e1_in_tree ) {
+                    nodes.erase( result[e].n2 );
+                    swap( result[N], result[e] );
+                    N++;
+                    break;
                 }
-                Gr.erase( largest );
             }
         }
-
-        return result;
     }
+
+    return result;
+}
+
+
+/// Use Prim's algorithm to construct a maximal spanning tree from the weighted graph Graph
+/// The weights should be positive. Use implementation in Boost Graph Library.
+template<typename T> DEdgeVec MaxSpanningTreePrims( const WeightedGraph<T> & Graph ) {
+    T maxweight = Graph.begin()->second;
+    for( typename WeightedGraph<T>::const_iterator it = Graph.begin(); it != Graph.end(); it++ )
+        if( it->second > maxweight )
+            maxweight = it->second;
+    // make a copy of the graph
+    WeightedGraph<T> gr( Graph );
+    // invoke MinSpanningTreePrims with negative weights
+    // (which have to be shifted to satisfy positivity criterion)
+    for( typename WeightedGraph<T>::iterator it = gr.begin(); it != gr.end(); it++ )
+        it->second = maxweight - it->second;
+    return MinSpanningTreePrims( gr );
 }
 
 
@@ -167,6 +189,87 @@ DEdgeVec GrowRootedTree( const Graph & T, size_t Root );
 UEdgeVec RandomDRegularGraph( size_t N, size_t d );
 
 
+/// Use Dijkstra's algorithm to solve the shortest path problem for the weighted graph Graph and source vertex index s
+template<typename T>
+std::map<size_t,size_t> DijkstraShortestPaths( const WeightedGraph<T> & Graph, size_t s ) {
+/*
+ * from wikipedia 
+ *
+In the following algorithm, u := Extract_Min(Q) searches for the vertex u in the vertex set Q that has the least d[u] value. That vertex is removed from the set Q and returned to the user.
+
+  function Dijkstra(G, w, s)
+     for each vertex v in V[G]                        // Initializations
+           d[v] := infinity                                 // Unknown distance function from s to v
+           previous[v] := undefined
+     d[s] := 0                                        // Distance from s to s
+     S := empty set                                   // Set of all visited vertices
+     Q := V[G]                                        // Set of all unvisited vertices
+     while Q is not an empty set                      // The algorithm itself
+           u := Extract_Min(Q)                        // Remove best vertex from priority queue
+           S := S union {u}                           // Mark it 'visited'
+           for each edge (u,v) outgoing from u
+                  if d[u] + w(u,v) < d[v]             // Relax (u,v)
+                        d[v] := d[u] + w(u,v) 
+                        previous[v] := u
+
+To find the shortest path from s to t:
+
+u := t
+while defined previous[u]
+      insert u to the beginning of S
+      u := previous[u]
+       
+*/
+
+    // Calculate set of nodes in G
+    std::set<size_t> nodes;
+    for( typename WeightedGraph<T>::const_iterator e = Graph.begin(); e != Graph.end(); e++ ) {
+        nodes.insert( e->n1.n1 );
+        nodes.insert( e->n1.n2 );
+    }
+    if( !nodes.count( s ) )
+        return std::map<size_t, size_t>();
+
+    std::map<size_t, double> d;
+    std::map<size_t, size_t> previous;
+    for( std::set<size_t>::const_iterator n = nodes.begin(); n != nodes.end(); n++ )
+        d[*n] = std::numeric_limits<double>::infinity();
+
+    d[s] = 0;
+    std::set<size_t> S;
+    std::set<size_t> Q = nodes;
+
+    while( Q.size() ) {
+        double least = d[*Q.begin()];
+        std::set<size_t>::iterator u_least = Q.begin();
+        for( std::set<size_t>::iterator _u = Q.begin(); _u != Q.end(); _u++ )
+            if( d[*_u] < least ) {
+                u_least = _u;
+                least = d[*_u];
+            }
+        size_t u = *u_least;
+        Q.erase( u_least );
+            
+        S.insert( u );
+        for( typename WeightedGraph<T>::const_iterator e = Graph.begin(); e != Graph.end(); e++ ) {
+            size_t v;
+            if( e->n1.n1 == u )
+                v = e->n1.n2;
+            else if( e->n1.n2 == u )
+                v = e->n1.n1;
+            else
+                continue;
+            if( d[u] + e->n2 < d[v] ) {
+                d[v] = d[u] + e->n2;
+                previous[v] = u;
+            }
+        }
+    }
+
+    return previous;
+}
+
+
 } // end of namespace dai
 
 
index 1045ded..79d751f 100644 (file)
@@ -191,7 +191,6 @@ double BP::run() {
     double tic = toc();
     Diffs diffs(nrVars(), 1.0);
     
-    typedef pair<size_t,size_t> Edge;
     vector<Edge> update_seq;
 
     vector<Factor> old_beliefs;
index 80e107e..4b0695e 100644 (file)
@@ -32,87 +32,84 @@ namespace dai {
 using namespace std;
 
 
-ClusterGraph ClusterGraph::VarElim( const std::vector<Var> & ElimSeq ) const {
-    const long verbose = 0;
+ClusterGraph::ClusterGraph( const std::vector<VarSet> & cls ) : G(), vars(), clusters() {
+    // construct vars, clusters and edge list
+    vector<Edge> edges;
+    foreach( const VarSet &cl, cls ) {
+        if( find( clusters.begin(), clusters.end(), cl ) == clusters.end() ) {
+            // add cluster
+            size_t n2 = clusters.size();
+            clusters.push_back( cl );
+            for( VarSet::const_iterator n = cl.begin(); n != cl.end(); n++ ) {
+                size_t n1 = find( vars.begin(), vars.end(), *n ) - vars.begin();
+                if( n1 == vars.size() )
+                    // add variable
+                    vars.push_back( *n );
+                edges.push_back( Edge( n1, n2 ) );
+            }
+        } // disregard duplicate clusters
+    }
+
+    // Create bipartite graph
+    G.create( vars.size(), clusters.size(), edges.begin(), edges.end() );
+}
+
 
+ClusterGraph ClusterGraph::VarElim_MinFill() const {
     // Make a copy
-    ClusterGraph _Cl(*this);
+    ClusterGraph cl(*this);
+    cl.eraseNonMaximal();
 
     ClusterGraph result;
-    _Cl.eraseNonMaximal();
+
+    // Construct set of variable indices
+    set<size_t> varindices;
+    for( size_t i = 0; i < vars.size(); ++i )
+        varindices.insert( i );
     
     // Do variable elimination
-    for( vector<Var>::const_iterator n = ElimSeq.begin(); n != ElimSeq.end(); n++ ) {
-        assert( _Cl.vars() && *n );
-
-        if( verbose >= 1 )
-            cout << "Cost of eliminating " << *n << ": " << _Cl.eliminationCost( *n ) << " new edges" << endl;
-        
-        result.insert( _Cl.Delta(*n) );
-
-        if( verbose >= 1 )
-            cout << "_Cl = " << _Cl << endl;
-
-        if( verbose >= 1 )
-            cout << "After inserting " << _Cl.delta(*n) << ", _Cl = ";
-        _Cl.insert( _Cl.delta(*n) );
-        if( verbose >= 1 )
-            cout << _Cl << endl;
-
-        if( verbose >= 1 )
-            cout << "After erasing clusters that contain " << *n <<  ", _Cl = ";
-        _Cl.eraseSubsuming( *n );
-        if( verbose >= 1 )
-            cout << _Cl << endl;
-
-        if( verbose >= 1 )
-            cout << "After erasing nonmaximal clusters, _Cl = ";
-        _Cl.eraseNonMaximal();
-        if( verbose >= 1 )
-            cout << _Cl << endl;
+    while( !varindices.empty() ) {
+        set<size_t>::const_iterator lowest = varindices.end();
+        size_t lowest_cost = -1UL;
+        for( set<size_t>::const_iterator i = varindices.begin(); i != varindices.end(); i++ ) {
+            size_t cost = cl.eliminationCost( *i );
+            if( lowest == varindices.end() || lowest_cost > cost ) {
+                lowest = i;
+                lowest_cost = cost;
+            }
+        }
+        size_t i = *lowest;
+
+        result.insert( cl.Delta( i ) );
+
+        cl.insert( cl.delta( i ) );
+        cl.eraseSubsuming( i );
+        cl.eraseNonMaximal();
+        varindices.erase( i );
     }
 
     return result;
 }
 
 
-ClusterGraph ClusterGraph::VarElim_MinFill() const {
-    const long verbose = 0;
 
+ClusterGraph ClusterGraph::VarElim( const std::vector<Var> & ElimSeq ) const {
     // Make a copy
-    ClusterGraph _Cl(*this);
-    VarSet _vars( vars() );
+    ClusterGraph cl(*this);
+    cl.eraseNonMaximal();
 
     ClusterGraph result;
-    _Cl.eraseNonMaximal();
     
     // Do variable elimination
-    while( !_vars.empty() ) {
-        if( verbose >= 1 )
-            cout << "Var  Eliminiation cost" << endl;
-        VarSet::const_iterator lowest = _vars.end();
-        size_t lowest_cost = -1UL;
-        for( VarSet::const_iterator n = _vars.begin(); n != _vars.end(); n++ ) {
-            size_t cost = _Cl.eliminationCost( *n );
-            if( verbose >= 1 )
-                cout << *n << "  " << cost << endl;
-            if( lowest == _vars.end() || lowest_cost > cost ) {
-                lowest = n;
-                lowest_cost = cost;
-            }
-        }
-        Var n = *lowest;
-
-        if( verbose >= 1 )
-            cout << "Lowest: " << n << " (" << lowest_cost << ")" << endl;
-
-        result.insert( _Cl.Delta(n) );
+    for( vector<Var>::const_iterator n = ElimSeq.begin(); n != ElimSeq.end(); n++ ) {
+        size_t i = cl.findVar( *n );
+        assert( i != cl.vars.size() );
 
-        _Cl.insert( _Cl.delta(n) );
-        _Cl.eraseSubsuming( n );
-        _Cl.eraseNonMaximal();
-        _vars /= n;
+        result.insert( cl.Delta(i) );
 
+        cl.insert( cl.delta(i) );
+        cl.eraseSubsuming( i );
+        cl.eraseNonMaximal();
     }
 
     return result;
index f2227bf..5d650e8 100644 (file)
@@ -68,7 +68,6 @@ void FactorGraph::createGraph( size_t nrEdges ) {
         hashmap[var(i).label()] = i;
     
     // create edge list
-    typedef pair<unsigned,unsigned> Edge;
     vector<Edge> edges;
     edges.reserve( nrEdges );
     for( size_t i2 = 0; i2 < nrFactors(); i2++ ) {
index 247635f..b846d35 100644 (file)
@@ -49,11 +49,13 @@ JTree::JTree( const FactorGraph &fg, const Properties &opts, bool automatic ) :
     assert( checkProperties() );
 
     if( automatic ) {
-        ClusterGraph _cg;
-
-        // Copy factors
+        // Copy VarSets of factors
+        vector<VarSet> cl;
+        cl.reserve( fg.nrFactors() );
         for( size_t I = 0; I < nrFactors(); I++ )
-            _cg.insert( factor(I).vars() );
+            cl.push_back( factor(I).vars() );
+        ClusterGraph _cg( cl );
+
         if( Verbose() >= 3 )
             cout << "Initial clusters: " << _cg << endl;
 
@@ -63,15 +65,8 @@ JTree::JTree( const FactorGraph &fg, const Properties &opts, bool automatic ) :
             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( Verbose() >= 3 )
+            cout << "VarElim_MinFill result: " << ElimVec << endl;
 
         GenerateJT( ElimVec );
     }
@@ -90,7 +85,7 @@ void JTree::GenerateJT( const std::vector<VarSet> &Cliques ) {
         }
     
     // Construct maximal spanning tree using Prim's algorithm
-    _RTree = MaxSpanningTreePrim( JuncGraph );
+    _RTree = MaxSpanningTreePrims( JuncGraph );
 
     // Construct corresponding region graph
 
@@ -115,7 +110,6 @@ void JTree::GenerateJT( const std::vector<VarSet> &Cliques ) {
 
     // Create inner regions and edges
     IRs.reserve( _RTree.size() );
-    typedef pair<size_t,size_t> Edge;
     vector<Edge> edges;
     edges.reserve( 2 * _RTree.size() );
     for( size_t i = 0; i < _RTree.size(); i++ ) {
index d0b45a2..57c84d5 100644 (file)
@@ -286,7 +286,6 @@ double LC::run() {
     }
 
     size_t nredges = nrEdges();
-    typedef pair<size_t,size_t> Edge;
     vector<Edge> update_seq;
     update_seq.reserve( nredges );
     for( size_t i = 0; i < nrVars(); ++i )
index 8449745..72b3eb9 100644 (file)
@@ -71,8 +71,8 @@ RegionGraph::RegionGraph( const FactorGraph &fg, const std::vector<VarSet> &cl )
     
     // Create outer regions, giving them counting number 1.0
     ORs.reserve( cg.size() );
-    for( ClusterGraph::const_iterator alpha = cg.begin(); alpha != cg.end(); alpha++ )
-        ORs.push_back( FRegion(Factor(*alpha, 1.0), 1.0) );
+    foreach( const VarSet &ns, cg.clusters )
+        ORs.push_back( FRegion(Factor(ns, 1.0), 1.0) );
 
     // For each factor, find an outer regions that subsumes that factor.
     // Then, multiply the outer region with that factor.
@@ -91,9 +91,9 @@ RegionGraph::RegionGraph( const FactorGraph &fg, const std::vector<VarSet> &cl )
     
     // Create inner regions - first pass
     set<VarSet> betas;
-    for( ClusterGraph::const_iterator alpha = cg.begin(); alpha != cg.end(); alpha++ )
-        for( ClusterGraph::const_iterator alpha2 = alpha; (++alpha2) != cg.end(); ) {
-            VarSet intersect = (*alpha) & (*alpha2);
+    for( size_t alpha = 0; alpha < cg.clusters.size(); alpha++ )
+        for( size_t alpha2 = alpha; (++alpha2) != cg.clusters.size(); ) {
+            VarSet intersect = cg.clusters[alpha] & cg.clusters[alpha2];
             if( intersect.size() > 0 )
                 betas.insert( intersect );
         }
index f076aa8..d5c4614 100644 (file)
@@ -217,10 +217,10 @@ TreeEP::TreeEP( const FactorGraph &fg, const Properties &opts ) : JTree(fg, opts
             }
 
             // find maximal spanning tree
-            ConstructRG( MaxSpanningTreePrim( wg ) );
+            ConstructRG( MaxSpanningTreePrims( wg ) );
 
 //            cout << "Constructing maximum spanning tree..." << endl;
-//            DEdgeVec MST = MaxSpanningTreePrim( wg );
+//            DEdgeVec MST = MaxSpanningTreePrims( wg );
 //            cout << "Maximum spanning tree:" << endl;
 //            for( DEdgeVec::const_iterator e = MST.begin(); e != MST.end(); e++ )
 //                cout << *e << endl; 
@@ -245,7 +245,7 @@ TreeEP::TreeEP( const FactorGraph &fg, const Properties &opts ) : JTree(fg, opts
             }
 
             // find maximal spanning tree
-            ConstructRG( MaxSpanningTreePrim( wg ) );
+            ConstructRG( MaxSpanningTreePrims( wg ) );
         } else {
             assert( 0 == 1 );
         }
@@ -269,7 +269,7 @@ void TreeEP::ConstructRG( const DEdgeVec &tree ) {
         }
     
     // Construct maximal spanning tree using Prim's algorithm
-    _RTree = MaxSpanningTreePrim( JuncGraph );
+    _RTree = MaxSpanningTreePrims( JuncGraph );
 
     // Construct corresponding region graph
 
@@ -297,7 +297,6 @@ void TreeEP::ConstructRG( const DEdgeVec &tree ) {
 
     // Create inner regions and edges
     IRs.reserve( _RTree.size() );
-    typedef pair<size_t,size_t> Edge;
     vector<Edge> edges;
     edges.reserve( 2 * _RTree.size() );
     for( size_t i = 0; i < _RTree.size(); i++ ) {