#define __defined_libdai_bipgraph_h
+#include <ostream>
#include <vector>
-#include <algorithm>
+#include <cassert>
#include <dai/util.h>
/// 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() {}
/// (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 );
/// (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 );
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;
}
}
+ /// 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;
+ }
};
#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
std::vector<Factor> factors;
typedef BipartiteGraph::Neighbor Neighbor;
typedef BipartiteGraph::Neighbors Neighbors;
+ typedef BipartiteGraph::Edge Edge;
protected:
std::map<size_t,Prob> _undoProbs;
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 ) {}
#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 {
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 );
}
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
double tic = toc();
Diffs diffs(nrVars(), 1.0);
- typedef pair<size_t,size_t> Edge;
vector<Edge> update_seq;
vector<Factor> old_beliefs;
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;
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++ ) {
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;
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 );
}
}
// Construct maximal spanning tree using Prim's algorithm
- _RTree = MaxSpanningTreePrim( JuncGraph );
+ _RTree = MaxSpanningTreePrims( JuncGraph );
// Construct corresponding region graph
// 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++ ) {
}
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 )
// 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.
// 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 );
}
}
// 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;
}
// find maximal spanning tree
- ConstructRG( MaxSpanningTreePrim( wg ) );
+ ConstructRG( MaxSpanningTreePrims( wg ) );
} else {
assert( 0 == 1 );
}
}
// Construct maximal spanning tree using Prim's algorithm
- _RTree = MaxSpanningTreePrim( JuncGraph );
+ _RTree = MaxSpanningTreePrims( JuncGraph );
// Construct corresponding region graph
// 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++ ) {
projects / libdai.git / commitdiff