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