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
5 This file is part of libDAI.
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.
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.
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
23 #ifndef __defined_libdai_weightedgraph_h
24 #define __defined_libdai_weightedgraph_h
34 #include <boost/graph/adjacency_list.hpp>
35 #include <boost/graph/prim_minimum_spanning_tree.hpp>
48 DEdge( size_t m1
, size_t m2
) : n1(m1
), n2(m2
) {}
49 bool operator==( const DEdge
&x
) const {
50 return ((n1
== x
.n1
) && (n2
== x
.n2
));
52 bool operator!=( const DEdge
&x
) const {
55 bool operator<( const DEdge
&x
) const {
56 return( (n1
< x
.n1
) || ((n1
== x
.n1
) && (n2
< x
.n2
)) );
58 friend std::ostream
& operator << (std::ostream
& os
, const DEdge
& e
) {
59 os
<< "(" << e
.n1
<< "," << e
.n2
<< ")";
71 UEdge( size_t m1
, size_t m2
) : n1(m1
), n2(m2
) {}
72 UEdge( const DEdge
& e
) : n1(e
.n1
), n2(e
.n2
) {}
73 bool operator==( const UEdge
&x
) {
74 return ((n1
== x
.n1
) && (n2
== x
.n2
)) || ((n1
== x
.n2
) && (n2
== x
.n1
));
76 bool operator<( const UEdge
&x
) const {
77 size_t s
= n1
, l
= n2
;
80 size_t xs
= x
.n1
, xl
= x
.n2
;
83 return( (s
< xs
) || ((s
== xs
) && (l
< xl
)) );
85 friend std::ostream
& operator << (std::ostream
& os
, const UEdge
& e
) {
87 os
<< "{" << e
.n1
<< "," << e
.n2
<< "}";
89 os
<< "{" << e
.n2
<< "," << e
.n1
<< "}";
95 typedef std::vector
<UEdge
> UEdgeVec
;
96 typedef std::vector
<DEdge
> DEdgeVec
;
97 template<class T
> class WeightedGraph
: public std::map
<UEdge
, T
> {};
98 typedef std::set
<UEdge
> Graph
;
101 /// Use Prim's algorithm to construct a minimal spanning tree from the weighted graph Graph
102 /// The weights should be positive. Use implementation in Boost Graph Library.
103 template<typename T
> DEdgeVec
MinSpanningTreePrims( const WeightedGraph
<T
> & Graph
) {
105 if( Graph
.size() > 0 ) {
106 using namespace boost
;
108 typedef adjacency_list
< vecS
, vecS
, undirectedS
, property
<vertex_distance_t
, int>, property
<edge_weight_t
, double> > boostGraph
;
109 typedef pair
<size_t, size_t> E
;
113 vector
<double> weights
;
114 edges
.reserve( Graph
.size() );
115 weights
.reserve( Graph
.size() );
116 for( typename WeightedGraph
<T
>::const_iterator e
= Graph
.begin(); e
!= Graph
.end(); e
++ ) {
117 weights
.push_back( e
->second
);
118 edges
.push_back( E( e
->first
.n1
, e
->first
.n2
) );
119 nodes
.insert( e
->first
.n1
);
120 nodes
.insert( e
->first
.n2
);
123 boostGraph
g( edges
.begin(), edges
.end(), weights
.begin(), nodes
.size() );
124 vector
< graph_traits
< boostGraph
>::vertex_descriptor
> p( num_vertices(g
) );
125 prim_minimum_spanning_tree( g
, &(p
[0]) );
127 // Store tree edges in result
128 result
.reserve( nodes
.size() - 1 );
130 for( size_t i
= 0; i
!= p
.size(); i
++ )
132 result
.push_back( DEdge( p
[i
], i
) );
136 // We have to store the minimum spanning tree in the right
137 // order, such that for all (i1, j1), (i2, j2) in result,
138 // if j1 == i2 then (i1, j1) comes before (i2, j2) in result.
139 // We do this by reordering the contents of result, effectively
140 // growing the tree starting at the root. At each step,
141 // result[0..N-1] are the edges already added to the tree,
142 // whereas the other elements of result still have to be added.
143 // The elements of nodes are the vertices that still have to
144 // be added to the tree.
146 // Start with the root
150 // Iteratively add edges and nodes to the growing tree
151 while( N
!= result
.size() ) {
152 for( size_t e
= N
; e
!= result
.size(); e
++ ) {
153 bool e1_in_tree
= !nodes
.count( result
[e
].n1
);
155 nodes
.erase( result
[e
].n2
);
156 swap( result
[N
], result
[e
] );
168 /// Use Prim's algorithm to construct a maximal spanning tree from the weighted graph Graph
169 /// The weights should be positive. Use implementation in Boost Graph Library.
170 template<typename T
> DEdgeVec
MaxSpanningTreePrims( const WeightedGraph
<T
> & Graph
) {
171 T maxweight
= Graph
.begin()->second
;
172 for( typename WeightedGraph
<T
>::const_iterator it
= Graph
.begin(); it
!= Graph
.end(); it
++ )
173 if( it
->second
> maxweight
)
174 maxweight
= it
->second
;
175 // make a copy of the graph
176 WeightedGraph
<T
> gr( Graph
);
177 // invoke MinSpanningTreePrims with negative weights
178 // (which have to be shifted to satisfy positivity criterion)
179 for( typename WeightedGraph
<T
>::iterator it
= gr
.begin(); it
!= gr
.end(); it
++ )
180 it
->second
= maxweight
- it
->second
;
181 return MinSpanningTreePrims( gr
);
185 /// Calculate rooted tree from a tree T and a root
186 DEdgeVec
GrowRootedTree( const Graph
& T
, size_t Root
);
189 UEdgeVec
RandomDRegularGraph( size_t N
, size_t d
);
192 /// Use Dijkstra's algorithm to solve the shortest path problem for the weighted graph Graph and source vertex index s
194 std::map
<size_t,size_t> DijkstraShortestPaths( const WeightedGraph
<T
> & Graph
, size_t s
) {
198 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.
200 function Dijkstra(G, w, s)
201 for each vertex v in V[G] // Initializations
202 d[v] := infinity // Unknown distance function from s to v
203 previous[v] := undefined
204 d[s] := 0 // Distance from s to s
205 S := empty set // Set of all visited vertices
206 Q := V[G] // Set of all unvisited vertices
207 while Q is not an empty set // The algorithm itself
208 u := Extract_Min(Q) // Remove best vertex from priority queue
209 S := S union {u} // Mark it 'visited'
210 for each edge (u,v) outgoing from u
211 if d[u] + w(u,v) < d[v] // Relax (u,v)
212 d[v] := d[u] + w(u,v)
215 To find the shortest path from s to t:
218 while defined previous[u]
219 insert u to the beginning of S
224 // Calculate set of nodes in G
225 std::set
<size_t> nodes
;
226 for( typename WeightedGraph
<T
>::const_iterator e
= Graph
.begin(); e
!= Graph
.end(); e
++ ) {
227 nodes
.insert( e
->n1
.n1
);
228 nodes
.insert( e
->n1
.n2
);
230 if( !nodes
.count( s
) )
231 return std::map
<size_t, size_t>();
233 std::map
<size_t, double> d
;
234 std::map
<size_t, size_t> previous
;
235 for( std::set
<size_t>::const_iterator n
= nodes
.begin(); n
!= nodes
.end(); n
++ )
236 d
[*n
] = std::numeric_limits
<double>::infinity();
240 std::set
<size_t> Q
= nodes
;
243 double least
= d
[*Q
.begin()];
244 std::set
<size_t>::iterator u_least
= Q
.begin();
245 for( std::set
<size_t>::iterator _u
= Q
.begin(); _u
!= Q
.end(); _u
++ )
246 if( d
[*_u
] < least
) {
254 for( typename WeightedGraph
<T
>::const_iterator e
= Graph
.begin(); e
!= Graph
.end(); e
++ ) {
258 else if( e
->n1
.n2
== u
)
262 if( d
[u
] + e
->n2
< d
[v
] ) {
273 } // end of namespace dai