Small changes
[libdai.git] / include / dai / weightedgraph.h
1 /* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
2 Radboud University Nijmegen, The Netherlands
3
4 This file is part of libDAI.
5
6 libDAI is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 libDAI is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with libDAI; if not, write to the Free Software
18 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
19 */
20
21
22 #ifndef __defined_libdai_weightedgraph_h
23 #define __defined_libdai_weightedgraph_h
24
25
26 #include <vector>
27 #include <map>
28 #include <iostream>
29 #include <set>
30 #include <cassert>
31 #include <limits>
32
33 #include <boost/graph/adjacency_list.hpp>
34 #include <boost/graph/prim_minimum_spanning_tree.hpp>
35
36
37
38 namespace dai {
39
40
41 /// Directed edge
42 class DEdge {
43 public:
44 size_t n1, n2;
45
46 DEdge() {}
47 DEdge( size_t m1, size_t m2 ) : n1(m1), n2(m2) {}
48 bool operator==( const DEdge &x ) const {
49 return ((n1 == x.n1) && (n2 == x.n2));
50 }
51 bool operator!=( const DEdge &x ) const {
52 return !(*this == x);
53 }
54 bool operator<( const DEdge &x ) const {
55 return( (n1 < x.n1) || ((n1 == x.n1) && (n2 < x.n2)) );
56 }
57 friend std::ostream & operator << (std::ostream & os, const DEdge & e) {
58 os << "(" << e.n1 << "," << e.n2 << ")";
59 return os;
60 }
61 };
62
63
64 /// Undirected edge
65 class UEdge {
66 public:
67 size_t n1, n2;
68
69 UEdge() {}
70 UEdge( size_t m1, size_t m2 ) : n1(m1), n2(m2) {}
71 UEdge( const DEdge & e ) : n1(e.n1), n2(e.n2) {}
72 bool operator==( const UEdge &x ) {
73 return ((n1 == x.n1) && (n2 == x.n2)) || ((n1 == x.n2) && (n2 == x.n1));
74 }
75 bool operator<( const UEdge &x ) const {
76 size_t s = n1, l = n2;
77 if( s > l )
78 std::swap( s, l );
79 size_t xs = x.n1, xl = x.n2;
80 if( xs > xl )
81 std::swap( xs, xl );
82 return( (s < xs) || ((s == xs) && (l < xl)) );
83 }
84 friend std::ostream & operator << (std::ostream & os, const UEdge & e) {
85 if( e.n1 < e.n2 )
86 os << "{" << e.n1 << "," << e.n2 << "}";
87 else
88 os << "{" << e.n2 << "," << e.n1 << "}";
89 return os;
90 }
91 };
92
93
94 typedef std::vector<UEdge> UEdgeVec;
95 typedef std::vector<DEdge> DEdgeVec;
96 template<class T> class WeightedGraph : public std::map<UEdge, T> {};
97 typedef std::set<UEdge> Graph;
98
99
100 /// Use Prim's algorithm to construct a minimal spanning tree from the weighted graph Graph
101 /// The weights should be positive. Use implementation in Boost Graph Library.
102 template<typename T> DEdgeVec MinSpanningTreePrims( const WeightedGraph<T> & Graph ) {
103 DEdgeVec result;
104 if( Graph.size() > 0 ) {
105 using namespace boost;
106 using namespace std;
107 typedef adjacency_list< vecS, vecS, undirectedS, property<vertex_distance_t, int>, property<edge_weight_t, double> > boostGraph;
108 typedef pair<size_t, size_t> E;
109
110 set<size_t> nodes;
111 vector<E> edges;
112 vector<double> weights;
113 edges.reserve( Graph.size() );
114 weights.reserve( Graph.size() );
115 for( typename WeightedGraph<T>::const_iterator e = Graph.begin(); e != Graph.end(); e++ ) {
116 weights.push_back( e->second );
117 edges.push_back( E( e->first.n1, e->first.n2 ) );
118 nodes.insert( e->first.n1 );
119 nodes.insert( e->first.n2 );
120 }
121
122 boostGraph g( edges.begin(), edges.end(), weights.begin(), nodes.size() );
123 vector< graph_traits< boostGraph >::vertex_descriptor > p( num_vertices(g) );
124 prim_minimum_spanning_tree( g, &(p[0]) );
125
126 // Store tree edges in result
127 result.reserve( nodes.size() - 1 );
128 size_t root = 0;
129 for( size_t i = 0; i != p.size(); i++ )
130 if( p[i] != i )
131 result.push_back( DEdge( p[i], i ) );
132 else
133 root = i;
134
135 // We have to store the minimum spanning tree in the right
136 // order, such that for all (i1, j1), (i2, j2) in result,
137 // if j1 == i2 then (i1, j1) comes before (i2, j2) in result.
138 // We do this by reordering the contents of result, effectively
139 // growing the tree starting at the root. At each step,
140 // result[0..N-1] are the edges already added to the tree,
141 // whereas the other elements of result still have to be added.
142 // The elements of nodes are the vertices that still have to
143 // be added to the tree.
144
145 // Start with the root
146 nodes.erase( root );
147 size_t N = 0;
148
149 // Iteratively add edges and nodes to the growing tree
150 while( N != result.size() ) {
151 for( size_t e = N; e != result.size(); e++ ) {
152 bool e1_in_tree = !nodes.count( result[e].n1 );
153 if( e1_in_tree ) {
154 nodes.erase( result[e].n2 );
155 swap( result[N], result[e] );
156 N++;
157 break;
158 }
159 }
160 }
161 }
162
163 return result;
164 }
165
166
167 /// Use Prim's algorithm to construct a maximal spanning tree from the weighted graph Graph
168 /// The weights should be positive. Use implementation in Boost Graph Library.
169 template<typename T> DEdgeVec MaxSpanningTreePrims( const WeightedGraph<T> & Graph ) {
170 T maxweight = Graph.begin()->second;
171 for( typename WeightedGraph<T>::const_iterator it = Graph.begin(); it != Graph.end(); it++ )
172 if( it->second > maxweight )
173 maxweight = it->second;
174 // make a copy of the graph
175 WeightedGraph<T> gr( Graph );
176 // invoke MinSpanningTreePrims with negative weights
177 // (which have to be shifted to satisfy positivity criterion)
178 for( typename WeightedGraph<T>::iterator it = gr.begin(); it != gr.end(); it++ )
179 it->second = maxweight - it->second;
180 return MinSpanningTreePrims( gr );
181 }
182
183
184 /// Calculate rooted tree from a tree T and a root
185 DEdgeVec GrowRootedTree( const Graph & T, size_t Root );
186
187
188 UEdgeVec RandomDRegularGraph( size_t N, size_t d );
189
190
191 /// Use Dijkstra's algorithm to solve the shortest path problem for the weighted graph Graph and source vertex index s
192 template<typename T>
193 std::map<size_t,size_t> DijkstraShortestPaths( const WeightedGraph<T> & Graph, size_t s ) {
194 /*
195 * from wikipedia
196 *
197 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.
198
199 function Dijkstra(G, w, s)
200 for each vertex v in V[G] // Initializations
201 d[v] := infinity // Unknown distance function from s to v
202 previous[v] := undefined
203 d[s] := 0 // Distance from s to s
204 S := empty set // Set of all visited vertices
205 Q := V[G] // Set of all unvisited vertices
206 while Q is not an empty set // The algorithm itself
207 u := Extract_Min(Q) // Remove best vertex from priority queue
208 S := S union {u} // Mark it 'visited'
209 for each edge (u,v) outgoing from u
210 if d[u] + w(u,v) < d[v] // Relax (u,v)
211 d[v] := d[u] + w(u,v)
212 previous[v] := u
213
214 To find the shortest path from s to t:
215
216 u := t
217 while defined previous[u]
218 insert u to the beginning of S
219 u := previous[u]
220
221 */
222
223 // Calculate set of nodes in G
224 std::set<size_t> nodes;
225 for( typename WeightedGraph<T>::const_iterator e = Graph.begin(); e != Graph.end(); e++ ) {
226 nodes.insert( e->n1.n1 );
227 nodes.insert( e->n1.n2 );
228 }
229 if( !nodes.count( s ) )
230 return std::map<size_t, size_t>();
231
232 std::map<size_t, double> d;
233 std::map<size_t, size_t> previous;
234 for( std::set<size_t>::const_iterator n = nodes.begin(); n != nodes.end(); n++ )
235 d[*n] = std::numeric_limits<double>::infinity();
236
237 d[s] = 0;
238 std::set<size_t> S;
239 std::set<size_t> Q = nodes;
240
241 while( Q.size() ) {
242 double least = d[*Q.begin()];
243 std::set<size_t>::iterator u_least = Q.begin();
244 for( std::set<size_t>::iterator _u = Q.begin(); _u != Q.end(); _u++ )
245 if( d[*_u] < least ) {
246 u_least = _u;
247 least = d[*_u];
248 }
249 size_t u = *u_least;
250 Q.erase( u_least );
251
252 S.insert( u );
253 for( typename WeightedGraph<T>::const_iterator e = Graph.begin(); e != Graph.end(); e++ ) {
254 size_t v;
255 if( e->n1.n1 == u )
256 v = e->n1.n2;
257 else if( e->n1.n2 == u )
258 v = e->n1.n1;
259 else
260 continue;
261 if( d[u] + e->n2 < d[v] ) {
262 d[v] = d[u] + e->n2;
263 previous[v] = u;
264 }
265 }
266 }
267
268 return previous;
269 }
270
271
272 } // end of namespace dai
273
274
275 #endif