New git HEAD version
[libdai.git] / include / dai / weightedgraph.h
1 /* This file is part of libDAI - http://www.libdai.org/
2 *
3 * Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
4 *
5 * Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
6 */
7
8
9 /** \file
10 * \brief Defines some utility functions for (weighted) undirected graphs, trees and rooted trees.
11 * \todo Improve general support for graphs and trees (in particular, a good tree implementation is needed).
12 */
13
14
15 #ifndef __defined_libdai_weightedgraph_h
16 #define __defined_libdai_weightedgraph_h
17
18
19 #include <vector>
20 #include <map>
21 #include <iostream>
22 #include <set>
23 #include <limits>
24 #include <climits> // Work-around for bug in boost graph library
25
26 #include <boost/graph/adjacency_list.hpp>
27 #include <boost/graph/prim_minimum_spanning_tree.hpp>
28 #include <boost/graph/kruskal_min_spanning_tree.hpp>
29
30 #include <dai/util.h>
31 #include <dai/exceptions.h>
32 #include <dai/graph.h>
33
34
35 namespace dai {
36
37
38 /// Represents a directed edge
39 class DEdge {
40 public:
41 /// First node index (source of edge)
42 size_t first;
43
44 /// Second node index (target of edge)
45 size_t second;
46
47 /// Default constructor
48 DEdge() : first(0), second(0) {}
49
50 /// Constructs a directed edge pointing from \a m1 to \a m2
51 DEdge( size_t m1, size_t m2 ) : first(m1), second(m2) {}
52
53 /// Tests for equality
54 bool operator==( const DEdge &x ) const { return ((first == x.first) && (second == x.second)); }
55
56 /// Smaller-than operator (performs lexicographical comparison)
57 bool operator<( const DEdge &x ) const {
58 return( (first < x.first) || ((first == x.first) && (second < x.second)) );
59 }
60
61 /// Writes a directed edge to an output stream
62 friend std::ostream & operator << (std::ostream & os, const DEdge & e) {
63 os << "(" << e.first << "->" << e.second << ")";
64 return os;
65 }
66
67 /// Formats a directed edge as a string
68 std::string toString() const {
69 std::stringstream ss;
70 ss << *this;
71 return ss.str();
72 }
73 };
74
75
76 /// Represents an undirected edge
77 class UEdge {
78 public:
79 /// First node index
80 size_t first;
81
82 /// Second node index
83 size_t second;
84
85 /// Default constructor
86 UEdge() : first(0), second(0) {}
87
88 /// Constructs an undirected edge between \a m1 and \a m2
89 UEdge( size_t m1, size_t m2 ) : first(m1), second(m2) {}
90
91 /// Construct from DEdge
92 UEdge( const DEdge &e ) : first(e.first), second(e.second) {}
93
94 /// Tests for inequality (disregarding the ordering of the nodes)
95 bool operator==( const UEdge &x ) {
96 return ((first == x.first) && (second == x.second)) || ((first == x.second) && (second == x.first));
97 }
98
99 /// Smaller-than operator
100 bool operator<( const UEdge &x ) const {
101 size_t s = std::min( first, second );
102 size_t l = std::max( first, second );
103 size_t xs = std::min( x.first, x.second );
104 size_t xl = std::max( x.first, x.second );
105 return( (s < xs) || ((s == xs) && (l < xl)) );
106 }
107
108 /// Writes an undirected edge to an output stream
109 friend std::ostream & operator << (std::ostream & os, const UEdge & e) {
110 if( e.first < e.second )
111 os << "{" << e.first << "--" << e.second << "}";
112 else
113 os << "{" << e.second << "--" << e.first << "}";
114 return os;
115 }
116
117 /// Formats an undirected edge as a string
118 std::string toString() const {
119 std::stringstream ss;
120 ss << *this;
121 return ss.str();
122 }
123 };
124
125
126 /// Represents an undirected graph, implemented as a std::set of undirected edges
127 class GraphEL : public std::set<UEdge> {
128 public:
129 /// Default constructor
130 GraphEL() {}
131
132 /// Construct from range of objects that can be cast to UEdge
133 template <class InputIterator>
134 GraphEL( InputIterator begin, InputIterator end ) {
135 insert( begin, end );
136 }
137
138 /// Construct from GraphAL
139 GraphEL( const GraphAL& G ) {
140 for( size_t n1 = 0; n1 < G.nrNodes(); n1++ )
141 bforeach( const Neighbor n2, G.nb(n1) )
142 if( n1 < n2 )
143 insert( UEdge( n1, n2 ) );
144 }
145 };
146
147
148 /// Represents an undirected weighted graph, with weights of type \a T, implemented as a std::map mapping undirected edges to weights
149 template<class T> class WeightedGraph : public std::map<UEdge, T> {};
150
151
152 /// Represents a rooted tree, implemented as a vector of directed edges
153 /** By convention, the edges are stored such that they point away from
154 * the root and such that edges nearer to the root come before edges
155 * farther away from the root.
156 */
157 class RootedTree : public std::vector<DEdge> {
158 public:
159 /// Default constructor
160 RootedTree() {}
161
162 /// Constructs a rooted tree from a tree and a root
163 /** \pre T has no cycles and contains node \a Root
164 */
165 RootedTree( const GraphEL &T, size_t Root );
166 };
167
168
169 /// Constructs a minimum spanning tree from the (non-negatively) weighted graph \a G.
170 /** \param G Weighted graph that should have non-negative weights.
171 * \param usePrim If true, use Prim's algorithm (complexity O(E log(V))), otherwise, use Kruskal's algorithm (complexity O(E log(E))).
172 * \note Uses implementation from Boost Graph Library.
173 * \note The vertices of \a G must be in the range [0,N) where N is the number of vertices of \a G.
174 */
175 template<typename T> RootedTree MinSpanningTree( const WeightedGraph<T> &G, bool usePrim ) {
176 RootedTree result;
177 if( G.size() > 0 ) {
178 using namespace boost;
179 using namespace std;
180 typedef adjacency_list< vecS, vecS, undirectedS, no_property, property<edge_weight_t, double> > boostGraph;
181
182 set<size_t> nodes;
183 vector<UEdge> edges;
184 vector<double> weights;
185 edges.reserve( G.size() );
186 weights.reserve( G.size() );
187 for( typename WeightedGraph<T>::const_iterator e = G.begin(); e != G.end(); e++ ) {
188 weights.push_back( e->second );
189 edges.push_back( e->first );
190 nodes.insert( e->first.first );
191 nodes.insert( e->first.second );
192 }
193
194 size_t N = nodes.size();
195 for( set<size_t>::const_iterator it = nodes.begin(); it != nodes.end(); it++ )
196 if( *it >= N )
197 DAI_THROWE(RUNTIME_ERROR,"Vertices must be in range [0..N) where N is the number of vertices.");
198
199 boostGraph g( edges.begin(), edges.end(), weights.begin(), nodes.size() );
200 size_t root = *(nodes.begin());
201 GraphEL tree;
202 if( usePrim ) {
203 // Prim's algorithm
204 vector< graph_traits< boostGraph >::vertex_descriptor > p(N);
205 prim_minimum_spanning_tree( g, &(p[0]) );
206
207 // Store tree edges in result
208 for( size_t i = 0; i != p.size(); i++ ) {
209 if( p[i] != i )
210 tree.insert( UEdge( p[i], i ) );
211 }
212 } else {
213 // Kruskal's algorithm
214 vector< graph_traits< boostGraph >::edge_descriptor > t;
215 t.reserve( N - 1 );
216 kruskal_minimum_spanning_tree( g, std::back_inserter(t) );
217
218 // Store tree edges in result
219 for( size_t i = 0; i != t.size(); i++ ) {
220 size_t v1 = source( t[i], g );
221 size_t v2 = target( t[i], g );
222 if( v1 != v2 )
223 tree.insert( UEdge( v1, v2 ) );
224 }
225 }
226
227 // Direct edges in order to obtain a rooted tree
228 result = RootedTree( tree, root );
229 }
230 return result;
231 }
232
233
234 /// Constructs a minimum spanning tree from the (non-negatively) weighted graph \a G.
235 /** \param G Weighted graph that should have non-negative weights.
236 * \param usePrim If true, use Prim's algorithm (complexity O(E log(V))), otherwise, use Kruskal's algorithm (complexity O(E log(E))).
237 * \note Uses implementation from Boost Graph Library.
238 * \note The vertices of \a G must be in the range [0,N) where N is the number of vertices of \a G.
239 */
240 template<typename T> RootedTree MaxSpanningTree( const WeightedGraph<T> &G, bool usePrim ) {
241 if( G.size() == 0 )
242 return RootedTree();
243 else {
244 T maxweight = G.begin()->second;
245 for( typename WeightedGraph<T>::const_iterator it = G.begin(); it != G.end(); it++ )
246 if( it->second > maxweight )
247 maxweight = it->second;
248 // make a copy of the graph
249 WeightedGraph<T> gr( G );
250 // invoke MinSpanningTree with negative weights
251 // (which have to be shifted to satisfy positivity criterion)
252 for( typename WeightedGraph<T>::iterator it = gr.begin(); it != gr.end(); it++ )
253 it->second = maxweight - it->second;
254 return MinSpanningTree( gr, usePrim );
255 }
256 }
257
258
259 } // end of namespace dai
260
261
262 #endif