graph/distributed/delta_stepping_shortest_paths.hpp

0001 // Copyright (C) 2007 The Trustees of Indiana University.
0002
0003 // Use, modification and distribution is subject to the Boost Software
0004 // License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0005 // http://www.boost.org/LICENSE_1_0.txt)
0006
0007 //  Authors: Douglas Gregor
0008 //           Andrew Lumsdaine
0009
0010 /**************************************************************************
0011  * This source file implements the Delta-stepping algorithm:              *
0012  *                                                                        *
0013  *   Ulrich Meyer and Peter Sanders. Parallel Shortest Path for Arbitrary *
0014  *   Graphs. In Proceedings from the 6th International Euro-Par           *
0015  *   Conference on Parallel Processing, pages 461--470, 2000.             *
0016  *                                                                        *
0017  *   Ulrich Meyer, Peter Sanders: [Delta]-stepping: A Parallelizable      *
0018  *   Shortest Path Algorithm. J. Algorithms 49(1): 114-152, 2003.         *
0019  *                                                                        *
0020  * There are several potential optimizations we could still perform for   *
0021  * this algorithm implementation:                                         *
0022  *                                                                        *
0023  *   - Implement "shortcuts", which bound the number of reinsertions      *
0024  *     in a single bucket (to one). The computation of shortcuts looks    *
0025  *     expensive in a distributed-memory setting, but it could be         *
0026  *     ammortized over many queries.                                      *
0027  *                                                                        *
0028  *   - The size of the "buckets" data structure can be bounded to         *
0029  *     max_edge_weight/delta buckets, if we treat the buckets as a        *
0030  *     circular array.                                                    *
0031  *                                                                        *
0032  *   - If we partition light/heavy edges ahead of time, we could improve  *
0033  *     relaxation performance by only iterating over the right portion    *
0034  *     of the out-edge list each time.                                    *
0035  *                                                                        *
0036  *   - Implement a more sophisticated algorithm for guessing "delta",     *
0037  *     based on the shortcut-finding algorithm.                           *
0038  **************************************************************************/
0039 #ifndef BOOST_GRAPH_DELTA_STEPPING_SHORTEST_PATHS_HPP
0040 #define BOOST_GRAPH_DELTA_STEPPING_SHORTEST_PATHS_HPP
0041
0042 #ifndef BOOST_GRAPH_USE_MPI
0043 #error "Parallel BGL files should not be included unless <boost/graph/use_mpi.hpp> has been included"
0044 #endif
0045
0046 #include <boost/config.hpp>
0047 #include <boost/graph/graph_traits.hpp>
0048 #include <boost/property_map/property_map.hpp>
0049 #include <boost/property_map/parallel/parallel_property_maps.hpp>
0050 #include <boost/graph/iteration_macros.hpp>
0051 #include <limits>
0052 #include <list>
0053 #include <vector>
0054 #include <boost/graph/parallel/container_traits.hpp>
0055 #include <boost/graph/parallel/properties.hpp>
0056 #include <boost/graph/distributed/detail/dijkstra_shortest_paths.hpp>
0057 #include <utility> // for std::pair
0058 #include <functional> // for std::logical_or
0059 #include <boost/graph/parallel/algorithm.hpp> // for all_reduce
0060 #include <cassert>
0061 #include <algorithm> // for std::min, std::max
0062 #include <boost/graph/parallel/simple_trigger.hpp>
0063
0064 #ifdef PBGL_DELTA_STEPPING_DEBUG
0065 #  include <iostream> // for std::cerr
0066 #endif
0067
0068 namespace boost { namespace graph { namespace distributed {
0069
0070 template<typename Graph, typename PredecessorMap, typename DistanceMap,
0071          typename EdgeWeightMap>
0072 class delta_stepping_impl {
0073   typedef typename graph_traits<Graph>::vertex_descriptor Vertex;
0074   typedef typename graph_traits<Graph>::degree_size_type Degree;
0075   typedef typename property_traits<EdgeWeightMap>::value_type Dist;
0076   typedef typename boost::graph::parallel::process_group_type<Graph>::type
0077     ProcessGroup;
0078
0079   typedef std::list<Vertex> Bucket;
0080   typedef typename Bucket::iterator BucketIterator;
0081   typedef typename std::vector<Bucket*>::size_type BucketIndex;
0082
0083   typedef detail::dijkstra_msg_value<DistanceMap, PredecessorMap> MessageValue;
0084
0085   enum {
0086     // Relax a remote vertex. The message contains a pair<Vertex,
0087     // MessageValue>, the first part of which is the vertex whose
0088     // tentative distance is being relaxed and the second part
0089     // contains either the new distance (if there is no predecessor
0090     // map) or a pair with the distance and predecessor.
0091     msg_relax
0092   };
0093
0094 public:
0095   delta_stepping_impl(const Graph& g,
0096                       PredecessorMap predecessor,
0097                       DistanceMap distance,
0098                       EdgeWeightMap weight,
0099                       Dist delta);
0100
0101   delta_stepping_impl(const Graph& g,
0102                       PredecessorMap predecessor,
0103                       DistanceMap distance,
0104                       EdgeWeightMap weight);
0105
0106   void run(Vertex s);
0107
0108 private:
0109   // Relax the edge (u, v), creating a new best path of distance x.
0110   void relax(Vertex u, Vertex v, Dist x);
0111
0112   // Synchronize all of the processes, by receiving all messages that
0113   // have not yet been received.
0114   void synchronize();
0115
0116   // Handle a relax message that contains only the target and
0117   // distance. This kind of message will be received when the
0118   // predecessor map is a dummy_property_map.
0119   void handle_relax_message(Vertex v, Dist x) { relax(v, v, x); }
0120
0121   // Handle a relax message that contains the source (predecessor),
0122   // target, and distance. This kind of message will be received when
0123   // the predecessor map is not a dummy_property_map.
0124   void handle_relax_message(Vertex v, const std::pair<Dist, Vertex>& p)
0125   { relax(p.second, v, p.first); }
0126
0127   // Setup triggers for msg_relax messages
0128   void setup_triggers();
0129
0130   void handle_msg_relax(int /*source*/, int /*tag*/,
0131                         const std::pair<Vertex, typename MessageValue::type>& data,
0132                         trigger_receive_context)
0133   { handle_relax_message(data.first, data.second); }
0134
0135   const Graph& g;
0136   PredecessorMap predecessor;
0137   DistanceMap distance;
0138   EdgeWeightMap weight;
0139   Dist delta;
0140   ProcessGroup pg;
0141   typename property_map<Graph, vertex_owner_t>::const_type owner;
0142   typename property_map<Graph, vertex_local_t>::const_type local;
0143
0144   // A "property map" that contains the position of each vertex in
0145   // whatever bucket it resides in.
0146   std::vector<BucketIterator> position_in_bucket;
0147
0148   // Bucket data structure. The ith bucket contains all local vertices
0149   // with (tentative) distance in the range [i*delta,
0150   // (i+1)*delta).
0151   std::vector<Bucket*> buckets;
0152
0153   // This "dummy" list is used only so that we can initialize the
0154   // position_in_bucket property map with non-singular iterators. This
0155   // won't matter for most implementations of the C++ Standard
0156   // Library, but it avoids undefined behavior and allows us to run
0157   // with library "debug modes".
0158   std::list<Vertex> dummy_list;
0159
0160   // A "property map" that states which vertices have been deleted
0161   // from the bucket in this iteration.
0162   std::vector<bool> vertex_was_deleted;
0163 };
0164
0165 template<typename Graph, typename PredecessorMap, typename DistanceMap,
0166          typename EdgeWeightMap>
0167 delta_stepping_impl<Graph, PredecessorMap, DistanceMap, EdgeWeightMap>::
0168 delta_stepping_impl(const Graph& g,
0169                     PredecessorMap predecessor,
0170                     DistanceMap distance,
0171                     EdgeWeightMap weight,
0172                     Dist delta)
0173     : g(g),
0174       predecessor(predecessor),
0175       distance(distance),
0176       weight(weight),
0177       delta(delta),
0178       pg(boost::graph::parallel::process_group_adl(g), attach_distributed_object()),
0179       owner(get(vertex_owner, g)),
0180       local(get(vertex_local, g))
0181 {
0182   setup_triggers();
0183 }
0184
0185 template<typename Graph, typename PredecessorMap, typename DistanceMap,
0186          typename EdgeWeightMap>
0187 delta_stepping_impl<Graph, PredecessorMap, DistanceMap, EdgeWeightMap>::
0188 delta_stepping_impl(const Graph& g,
0189                     PredecessorMap predecessor,
0190                     DistanceMap distance,
0191                     EdgeWeightMap weight)
0192     : g(g),
0193       predecessor(predecessor),
0194       distance(distance),
0195       weight(weight),
0196       pg(boost::graph::parallel::process_group_adl(g), attach_distributed_object()),
0197       owner(get(vertex_owner, g)),
0198       local(get(vertex_local, g))
0199 {
0200   using boost::parallel::all_reduce;
0201   using boost::parallel::maximum;
0202   using std::max;
0203
0204   // Compute the maximum edge weight and degree
0205   Dist max_edge_weight = 0;
0206   Degree max_degree = 0;
0207   BGL_FORALL_VERTICES_T(u, g, Graph) {
0208     max_degree = max BOOST_PREVENT_MACRO_SUBSTITUTION (max_degree, out_degree(u, g));
0209     BGL_FORALL_OUTEDGES_T(u, e, g, Graph)
0210       max_edge_weight = max BOOST_PREVENT_MACRO_SUBSTITUTION (max_edge_weight, get(weight, e));
0211   }
0212
0213   max_edge_weight = all_reduce(pg, max_edge_weight, maximum<Dist>());
0214   max_degree = all_reduce(pg, max_degree, maximum<Degree>());
0215
0216   // Take a guess at delta, based on what works well for random
0217   // graphs.
0218   delta = max_edge_weight / max_degree;
0219   if (delta == 0)
0220     delta = 1;
0221
0222   setup_triggers();
0223 }
0224
0225 template<typename Graph, typename PredecessorMap, typename DistanceMap,
0226          typename EdgeWeightMap>
0227 void
0228 delta_stepping_impl<Graph, PredecessorMap, DistanceMap, EdgeWeightMap>::
0229 run(Vertex s)
0230 {
0231   Dist inf = (std::numeric_limits<Dist>::max)();
0232
0233   // None of the vertices are stored in the bucket.
0234   position_in_bucket.clear();
0235   position_in_bucket.resize(num_vertices(g), dummy_list.end());
0236
0237   // None of the vertices have been deleted
0238   vertex_was_deleted.clear();
0239   vertex_was_deleted.resize(num_vertices(g), false);
0240
0241   // No path from s to any other vertex, yet
0242   BGL_FORALL_VERTICES_T(v, g, Graph)
0243     put(distance, v, inf);
0244
0245   // The distance to the starting node is zero
0246   if (get(owner, s) == process_id(pg))
0247     // Put "s" into its bucket (bucket 0)
0248     relax(s, s, 0);
0249   else
0250     // Note that we know the distance to s is zero
0251     cache(distance, s, 0);
0252
0253   BucketIndex max_bucket = (std::numeric_limits<BucketIndex>::max)();
0254   BucketIndex current_bucket = 0;
0255   do {
0256     // Synchronize with all of the other processes.
0257     synchronize();
0258
0259     // Find the next bucket that has something in it.
0260     while (current_bucket < buckets.size()
0261            && (!buckets[current_bucket] || buckets[current_bucket]->empty()))
0262       ++current_bucket;
0263     if (current_bucket >= buckets.size())
0264       current_bucket = max_bucket;
0265
0266 #ifdef PBGL_DELTA_STEPPING_DEBUG
0267     std::cerr << "#" << process_id(pg) << ": lowest bucket is #"
0268               << current_bucket << std::endl;
0269 #endif
0270     // Find the smallest bucket (over all processes) that has vertices
0271     // that need to be processed.
0272     using boost::parallel::all_reduce;
0273     using boost::parallel::minimum;
0274     current_bucket = all_reduce(pg, current_bucket, minimum<BucketIndex>());
0275
0276     if (current_bucket == max_bucket)
0277       // There are no non-empty buckets in any process; exit.
0278       break;
0279
0280 #ifdef PBGL_DELTA_STEPPING_DEBUG
0281     if (process_id(pg) == 0)
0282       std::cerr << "Processing bucket #" << current_bucket << std::endl;
0283 #endif
0284
0285     // Contains the set of vertices that have been deleted in the
0286     // relaxation of "light" edges. Note that we keep track of which
0287     // vertices were deleted with the property map
0288     // "vertex_was_deleted".
0289     std::vector<Vertex> deleted_vertices;
0290
0291     // Repeatedly relax light edges
0292     bool nonempty_bucket;
0293     do {
0294       // Someone has work to do in this bucket.
0295
0296       if (current_bucket < buckets.size() && buckets[current_bucket]) {
0297         Bucket& bucket = *buckets[current_bucket];
0298         // For each element in the bucket
0299         while (!bucket.empty()) {
0300           Vertex u = bucket.front();
0301
0302 #ifdef PBGL_DELTA_STEPPING_DEBUG
0303           std::cerr << "#" << process_id(pg) << ": processing vertex "
0304                     << get(vertex_global, g, u).second << "@"
0305                     << get(vertex_global, g, u).first
0306                     << std::endl;
0307 #endif
0308
0309           // Remove u from the front of the bucket
0310           bucket.pop_front();
0311
0312           // Insert u into the set of deleted vertices, if it hasn't
0313           // been done already.
0314           if (!vertex_was_deleted[get(local, u)]) {
0315             vertex_was_deleted[get(local, u)] = true;
0316             deleted_vertices.push_back(u);
0317           }
0318
0319           // Relax each light edge.
0320           Dist u_dist = get(distance, u);
0321           BGL_FORALL_OUTEDGES_T(u, e, g, Graph)
0322             if (get(weight, e) <= delta) // light edge
0323               relax(u, target(e, g), u_dist + get(weight, e));
0324         }
0325       }
0326
0327       // Synchronize with all of the other processes.
0328       synchronize();
0329
0330       // Is the bucket empty now?
0331       nonempty_bucket = (current_bucket < buckets.size()
0332                          && buckets[current_bucket]
0333                          && !buckets[current_bucket]->empty());
0334      } while (all_reduce(pg, nonempty_bucket, std::logical_or<bool>()));
0335
0336     // Relax heavy edges for each of the vertices that we previously
0337     // deleted.
0338     for (typename std::vector<Vertex>::iterator iter = deleted_vertices.begin();
0339          iter != deleted_vertices.end(); ++iter) {
0340       // Relax each heavy edge.
0341       Vertex u = *iter;
0342       Dist u_dist = get(distance, u);
0343       BGL_FORALL_OUTEDGES_T(u, e, g, Graph)
0344         if (get(weight, e) > delta) // heavy edge
0345           relax(u, target(e, g), u_dist + get(weight, e));
0346     }
0347
0348     // Go to the next bucket: the current bucket must already be empty.
0349     ++current_bucket;
0350   } while (true);
0351
0352   // Delete all of the buckets.
0353   for (typename std::vector<Bucket*>::iterator iter = buckets.begin();
0354        iter != buckets.end(); ++iter) {
0355     if (*iter) {
0356       delete *iter;
0357       *iter = 0;
0358     }
0359   }
0360 }
0361
0362 template<typename Graph, typename PredecessorMap, typename DistanceMap,
0363          typename EdgeWeightMap>
0364 void
0365 delta_stepping_impl<Graph, PredecessorMap, DistanceMap, EdgeWeightMap>::
0366 relax(Vertex u, Vertex v, Dist x)
0367 {
0368 #ifdef PBGL_DELTA_STEPPING_DEBUG
0369   std::cerr << "#" << process_id(pg) << ": relax("
0370             << get(vertex_global, g, u).second << "@"
0371             << get(vertex_global, g, u).first << ", "
0372             << get(vertex_global, g, v).second << "@"
0373             << get(vertex_global, g, v).first << ", "
0374             << x << ")" << std::endl;
0375 #endif
0376
0377   if (x < get(distance, v)) {
0378     // We're relaxing the edge to vertex v.
0379     if (get(owner, v) == process_id(pg)) {
0380       // Compute the new bucket index for v
0381       BucketIndex new_index = static_cast<BucketIndex>(x / delta);
0382
0383       // Make sure there is enough room in the buckets data structure.
0384       if (new_index >= buckets.size()) buckets.resize(new_index + 1, 0);
0385
0386       // Make sure that we have allocated the bucket itself.
0387       if (!buckets[new_index]) buckets[new_index] = new Bucket;
0388
0389       if (get(distance, v) != (std::numeric_limits<Dist>::max)()
0390           && !vertex_was_deleted[get(local, v)]) {
0391         // We're moving v from an old bucket into a new one. Compute
0392         // the old index, then splice it in.
0393         BucketIndex old_index
0394           = static_cast<BucketIndex>(get(distance, v) / delta);
0395         buckets[new_index]->splice(buckets[new_index]->end(),
0396                                    *buckets[old_index],
0397                                    position_in_bucket[get(local, v)]);
0398       } else {
0399         // We're inserting v into a bucket for the first time. Put it
0400         // at the end.
0401         buckets[new_index]->push_back(v);
0402       }
0403
0404       // v is now at the last position in the new bucket
0405       position_in_bucket[get(local, v)] = buckets[new_index]->end();
0406       --position_in_bucket[get(local, v)];
0407
0408       // Update predecessor and tentative distance information
0409       put(predecessor, v, u);
0410       put(distance, v, x);
0411     } else {
0412 #ifdef PBGL_DELTA_STEPPING_DEBUG
0413       std::cerr << "#" << process_id(pg) << ": sending relax("
0414                 << get(vertex_global, g, u).second << "@"
0415                 << get(vertex_global, g, u).first << ", "
0416                 << get(vertex_global, g, v).second << "@"
0417                 << get(vertex_global, g, v).first << ", "
0418             << x << ") to " << get(owner, v) << std::endl;
0419
0420 #endif
0421       // The vertex is remote: send a request to the vertex's owner
0422       send(pg, get(owner, v), msg_relax,
0423            std::make_pair(v, MessageValue::create(x, u)));
0424
0425       // Cache tentative distance information
0426       cache(distance, v, x);
0427     }
0428   }
0429 }
0430
0431 template<typename Graph, typename PredecessorMap, typename DistanceMap,
0432          typename EdgeWeightMap>
0433 void
0434 delta_stepping_impl<Graph, PredecessorMap, DistanceMap, EdgeWeightMap>::
0435 synchronize()
0436 {
0437   using boost::parallel::synchronize;
0438
0439   // Synchronize at the process group level.
0440   synchronize(pg);
0441
0442   // Receive any relaxation request messages.
0443 //   typedef typename ProcessGroup::process_id_type process_id_type;
0444 //   while (optional<std::pair<process_id_type, int> > stp = probe(pg)) {
0445 //     // Receive the relaxation message
0446 //     assert(stp->second == msg_relax);
0447 //     std::pair<Vertex, typename MessageValue::type> data;
0448 //     receive(pg, stp->first, stp->second, data);
0449
0450 //     // Turn it into a "relax" call
0451 //     handle_relax_message(data.first, data.second);
0452 //   }
0453 }
0454
0455 template<typename Graph, typename PredecessorMap, typename DistanceMap,
0456          typename EdgeWeightMap>
0457 void
0458 delta_stepping_impl<Graph, PredecessorMap, DistanceMap, EdgeWeightMap>::
0459 setup_triggers()
0460 {
0461   using boost::graph::parallel::simple_trigger;
0462
0463   simple_trigger(pg, msg_relax, this,
0464                  &delta_stepping_impl::handle_msg_relax);
0465 }
0466
0467 template<typename Graph, typename PredecessorMap, typename DistanceMap,
0468          typename EdgeWeightMap>
0469 void
0470 delta_stepping_shortest_paths
0471   (const Graph& g,
0472    typename graph_traits<Graph>::vertex_descriptor s,
0473    PredecessorMap predecessor, DistanceMap distance, EdgeWeightMap weight,
0474    typename property_traits<EdgeWeightMap>::value_type delta)
0475 {
0476   // The "distance" map needs to act like one, retrieving the default
0477   // value of infinity.
0478   set_property_map_role(vertex_distance, distance);
0479
0480   // Construct the implementation object, which will perform all of
0481   // the actual work.
0482   delta_stepping_impl<Graph, PredecessorMap, DistanceMap, EdgeWeightMap>
0483     impl(g, predecessor, distance, weight, delta);
0484
0485   // Run the delta-stepping algorithm. The results will show up in
0486   // "predecessor" and "weight".
0487   impl.run(s);
0488 }
0489
0490 template<typename Graph, typename PredecessorMap, typename DistanceMap,
0491          typename EdgeWeightMap>
0492 void
0493 delta_stepping_shortest_paths
0494   (const Graph& g,
0495    typename graph_traits<Graph>::vertex_descriptor s,
0496    PredecessorMap predecessor, DistanceMap distance, EdgeWeightMap weight)
0497 {
0498   // The "distance" map needs to act like one, retrieving the default
0499   // value of infinity.
0500   set_property_map_role(vertex_distance, distance);
0501
0502   // Construct the implementation object, which will perform all of
0503   // the actual work.
0504   delta_stepping_impl<Graph, PredecessorMap, DistanceMap, EdgeWeightMap>
0505     impl(g, predecessor, distance, weight);
0506
0507   // Run the delta-stepping algorithm. The results will show up in
0508   // "predecessor" and "weight".
0509   impl.run(s);
0510 }
0511
0512 } } } // end namespace boost::graph::distributed
0513
0514 #endif // BOOST_GRAPH_DELTA_STEPPING_SHORTEST_PATHS_HPP