Merge branch 'pletscher'
[libdai.git] / include / dai / weightedgraph.h
1 /* This file is part of libDAI - http://www.libdai.org/
2 *
3 * libDAI is licensed under the terms of the GNU General Public License version
4 * 2, or (at your option) any later version. libDAI is distributed without any
5 * warranty. See the file COPYING for more details.
6 *
7 * Copyright (C) 2006-2009 Joris Mooij [joris dot mooij at libdai dot org]
8 * Copyright (C) 2006-2007 Radboud University Nijmegen, The Netherlands
9 */
10
11
12 /** \file
13 * \brief Defines some utility functions for weighted graphs
14 * \todo Improve documentation
15 * \todo Improve general support for graphs and trees.
16 */
17
18
19 #ifndef __defined_libdai_weightedgraph_h
20 #define __defined_libdai_weightedgraph_h
21
22
23 #include <vector>
24 #include <map>
25 #include <iostream>
26 #include <set>
27 #include <limits>
28 #include <climits> // Work-around for bug in boost graph library
29
30 #include <boost/graph/adjacency_list.hpp>
31 #include <boost/graph/prim_minimum_spanning_tree.hpp>
32
33
34 namespace dai {
35
36
37 /// Represents a directed edge pointing from n1 to n2
38 class DEdge {
39 public:
40 size_t n1; ///< First node index
41 size_t n2; ///< Second node index
42
43 /// Default constructor
44 DEdge() {}
45
46 /// Constructor
47 DEdge( size_t m1, size_t m2 ) : n1(m1), n2(m2) {}
48
49 /// Tests for equality
50 bool operator==( const DEdge &x ) const { return ((n1 == x.n1) && (n2 == x.n2)); }
51
52 /// Tests for inequality
53 bool operator!=( const DEdge &x ) const { return !(*this == x); }
54
55 /// Smaller-than operator (performs lexicographical comparison)
56 bool operator<( const DEdge &x ) const {
57 return( (n1 < x.n1) || ((n1 == x.n1) && (n2 < x.n2)) );
58 }
59
60 /// Writes a DEdge to an output stream
61 friend std::ostream & operator << (std::ostream & os, const DEdge & e) {
62 os << "(" << e.n1 << "," << e.n2 << ")";
63 return os;
64 }
65 };
66
67
68 /// Undirected edge between nodes n1 and n2
69 class UEdge {
70 public:
71 size_t n1; ///< First node index
72 size_t n2; ///< Second node index
73
74 /// Default constructor
75 UEdge() {}
76
77 /// Constructor
78 UEdge( size_t m1, size_t m2 ) : n1(m1), n2(m2) {}
79
80 /// Construct from DEdge
81 UEdge( const DEdge & e ) : n1(e.n1), n2(e.n2) {}
82
83 /// Tests for inequality (disregarding the ordering of n1 and n2)
84 bool operator==( const UEdge &x ) {
85 return ((n1 == x.n1) && (n2 == x.n2)) || ((n1 == x.n2) && (n2 == x.n1));
86 }
87
88 /// Smaller-than operator
89 bool operator<( const UEdge &x ) const {
90 size_t s = n1, l = n2;
91 if( s > l )
92 std::swap( s, l );
93 size_t xs = x.n1, xl = x.n2;
94 if( xs > xl )
95 std::swap( xs, xl );
96 return( (s < xs) || ((s == xs) && (l < xl)) );
97 }
98
99 /// Writes a UEdge to an output stream
100 friend std::ostream & operator << (std::ostream & os, const UEdge & e) {
101 if( e.n1 < e.n2 )
102 os << "{" << e.n1 << "," << e.n2 << "}";
103 else
104 os << "{" << e.n2 << "," << e.n1 << "}";
105 return os;
106 }
107 };
108
109
110 /// Vector of UEdge
111 typedef std::vector<UEdge> UEdgeVec;
112
113 /// Vector of DEdge
114 typedef std::vector<DEdge> DEdgeVec;
115
116 /// Represents an undirected weighted graph, with weights of type T
117 template<class T> class WeightedGraph : public std::map<UEdge, T> {};
118
119 /// Represents an undirected graph
120 typedef std::set<UEdge> Graph;
121
122
123 /// Uses Prim's algorithm to construct a minimal spanning tree from the (positively) weighted graph G.
124 /** Uses implementation in Boost Graph Library.
125 */
126 template<typename T> DEdgeVec MinSpanningTreePrims( const WeightedGraph<T> &G ) {
127 DEdgeVec result;
128 if( G.size() > 0 ) {
129 using namespace boost;
130 using namespace std;
131 typedef adjacency_list< vecS, vecS, undirectedS, property<vertex_distance_t, int>, property<edge_weight_t, double> > boostGraph;
132 typedef pair<size_t, size_t> E;
133
134 set<size_t> nodes;
135 vector<E> edges;
136 vector<double> weights;
137 edges.reserve( G.size() );
138 weights.reserve( G.size() );
139 for( typename WeightedGraph<T>::const_iterator e = G.begin(); e != G.end(); e++ ) {
140 weights.push_back( e->second );
141 edges.push_back( E( e->first.n1, e->first.n2 ) );
142 nodes.insert( e->first.n1 );
143 nodes.insert( e->first.n2 );
144 }
145
146 boostGraph g( edges.begin(), edges.end(), weights.begin(), nodes.size() );
147 vector< graph_traits< boostGraph >::vertex_descriptor > p( num_vertices(g) );
148 prim_minimum_spanning_tree( g, &(p[0]) );
149
150 // Store tree edges in result
151 result.reserve( nodes.size() - 1 );
152 size_t root = 0;
153 for( size_t i = 0; i != p.size(); i++ )
154 if( p[i] != i )
155 result.push_back( DEdge( p[i], i ) );
156 else
157 root = i;
158
159 // We have to store the minimum spanning tree in the right
160 // order, such that for all (i1, j1), (i2, j2) in result,
161 // if j1 == i2 then (i1, j1) comes before (i2, j2) in result.
162 // We do this by reordering the contents of result, effectively
163 // growing the tree starting at the root. At each step,
164 // result[0..N-1] are the edges already added to the tree,
165 // whereas the other elements of result still have to be added.
166 // The elements of nodes are the vertices that still have to
167 // be added to the tree.
168
169 // Start with the root
170 nodes.erase( root );
171 size_t N = 0;
172
173 // Iteratively add edges and nodes to the growing tree
174 while( N != result.size() ) {
175 for( size_t e = N; e != result.size(); e++ ) {
176 bool e1_in_tree = !nodes.count( result[e].n1 );
177 if( e1_in_tree ) {
178 nodes.erase( result[e].n2 );
179 swap( result[N], result[e] );
180 N++;
181 break;
182 }
183 }
184 }
185 }
186
187 return result;
188 }
189
190
191 /// Use Prim's algorithm to construct a minimal spanning tree from the (positively) weighted graph G.
192 /** Uses implementation in Boost Graph Library.
193 */
194 template<typename T> DEdgeVec MaxSpanningTreePrims( const WeightedGraph<T> & Graph ) {
195 if( Graph.size() == 0 )
196 return DEdgeVec();
197 else {
198 T maxweight = Graph.begin()->second;
199 for( typename WeightedGraph<T>::const_iterator it = Graph.begin(); it != Graph.end(); it++ )
200 if( it->second > maxweight )
201 maxweight = it->second;
202 // make a copy of the graph
203 WeightedGraph<T> gr( Graph );
204 // invoke MinSpanningTreePrims with negative weights
205 // (which have to be shifted to satisfy positivity criterion)
206 for( typename WeightedGraph<T>::iterator it = gr.begin(); it != gr.end(); it++ )
207 it->second = maxweight - it->second;
208 return MinSpanningTreePrims( gr );
209 }
210 }
211
212
213 /// Constructs a rooted tree from a tree and a root
214 DEdgeVec GrowRootedTree( const Graph & T, size_t Root );
215
216
217 /// Constructs a random undirected graph of N nodes, where each node has connectivity d
218 UEdgeVec RandomDRegularGraph( size_t N, size_t d );
219
220
221 } // end of namespace dai
222
223
224 #endif