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 #ifndef __defined_libdai_weightedgraph_h
  24 #define __defined_libdai_weightedgraph_h
  25
  26
  27 #include <vector>
  28 #include <map>
  29 #include <iostream>
  30 #include <set>
  31 #include <cassert>
  32 #include <limits>
  33
  34 #include <boost/graph/adjacency_list.hpp>
  35 #include <boost/graph/prim_minimum_spanning_tree.hpp>
  36
  37
  38
  39 namespace dai {
  40
  41
  42 /// Directed edge
  43 class DEdge {
  44     public:
  45         size_t  n1, n2;
  46
  47         DEdge() {}
  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));
  51         }
  52         bool operator!=( const DEdge &x ) const {
  53             return !(*this == x);
  54         }
  55         bool operator<( const DEdge &x ) const {
  56             return( (n1 < x.n1) || ((n1 == x.n1) && (n2 < x.n2)) );
  57         }
  58         friend std::ostream & operator << (std::ostream & os, const DEdge & e) {
  59             os << "(" << e.n1 << "," << e.n2 << ")";
  60             return os;
  61         }
  62 };
  63
  64
  65 /// Undirected edge
  66 class UEdge {
  67     public:
  68         size_t  n1, n2;
  69
  70         UEdge() {}
  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));
  75         }
  76         bool operator<( const UEdge &x ) const {
  77             size_t s = n1, l = n2;
  78             if( s > l )
  79                 std::swap( s, l );
  80             size_t xs = x.n1, xl = x.n2;
  81             if( xs > xl )
  82                 std::swap( xs, xl );
  83             return( (s < xs) || ((s == xs) && (l < xl)) );
  84         }
  85         friend std::ostream & operator << (std::ostream & os, const UEdge & e) {
  86             if( e.n1 < e.n2 )
  87                 os << "{" << e.n1 << "," << e.n2 << "}";
  88             else
  89                 os << "{" << e.n2 << "," << e.n1 << "}";
  90             return os;
  91         }
  92 };
  93
  94
  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;
  99
 100
 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 ) {
 104     DEdgeVec result;
 105     if( Graph.size() > 0 ) {
 106         using namespace boost;
 107         using namespace std;
 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;
 110
 111         set<size_t> nodes;
 112         vector<E> edges;
 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 );
 121         }
 122
 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]) );
 126
 127         // Store tree edges in result
 128         result.reserve( nodes.size() - 1 );
 129         size_t root = 0;
 130         for( size_t i = 0; i != p.size(); i++ )
 131             if( p[i] != i )
 132                 result.push_back( DEdge( p[i], i ) );
 133             else
 134                 root = i;
 135
 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.
 145
 146         // Start with the root
 147         nodes.erase( root );
 148         size_t N = 0;
 149
 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 );
 154                 if( e1_in_tree ) {
 155                     nodes.erase( result[e].n2 );
 156                     swap( result[N], result[e] );
 157                     N++;
 158                     break;
 159                 }
 160             }
 161         }
 162     }
 163
 164     return result;
 165 }
 166
 167
 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 );
 182 }
 183
 184
 185 /// Calculate rooted tree from a tree T and a root
 186 DEdgeVec GrowRootedTree( const Graph & T, size_t Root );
 187
 188
 189 UEdgeVec RandomDRegularGraph( size_t N, size_t d );
 190
 191
 192 /// Use Dijkstra's algorithm to solve the shortest path problem for the weighted graph Graph and source vertex index s
 193 template<typename T>
 194 std::map<size_t,size_t> DijkstraShortestPaths( const WeightedGraph<T> & Graph, size_t s ) {
 195 /*
 196  * from wikipedia
 197  *
 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.
 199
 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)
 213                         previous[v] := u
 214
 215 To find the shortest path from s to t:
 216
 217 u := t
 218 while defined previous[u]
 219       insert u to the beginning of S
 220       u := previous[u]
 221
 222 */
 223
 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 );
 229     }
 230     if( !nodes.count( s ) )
 231         return std::map<size_t, size_t>();
 232
 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();
 237
 238     d[s] = 0;
 239     std::set<size_t> S;
 240     std::set<size_t> Q = nodes;
 241
 242     while( Q.size() ) {
 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 ) {
 247                 u_least = _u;
 248                 least = d[*_u];
 249             }
 250         size_t u = *u_least;
 251         Q.erase( u_least );
 252
 253         S.insert( u );
 254         for( typename WeightedGraph<T>::const_iterator e = Graph.begin(); e != Graph.end(); e++ ) {
 255             size_t v;
 256             if( e->n1.n1 == u )
 257                 v = e->n1.n2;
 258             else if( e->n1.n2 == u )
 259                 v = e->n1.n1;
 260             else
 261                 continue;
 262             if( d[u] + e->n2 < d[v] ) {
 263                 d[v] = d[u] + e->n2;
 264                 previous[v] = u;
 265             }
 266         }
 267     }
 268
 269     return previous;
 270 }
 271
 272
 273 } // end of namespace dai
 274
 275
 276 #endif