ceph/src/boost/libs/graph_parallel/include/boost/graph/distributed/betweenness_centrality.hpp

   1 // Copyright 2004 The Trustees of Indiana University.
   2
   3 // Distributed under the Boost Software License, Version 1.0.
   4 // (See accompanying file LICENSE_1_0.txt or copy at
   5 // http://www.boost.org/LICENSE_1_0.txt)
   6
   7 //  Authors: Douglas Gregor
   8 //           Andrew Lumsdaine
   9 #ifndef BOOST_GRAPH_PARALLEL_BRANDES_BETWEENNESS_CENTRALITY_HPP
  10 #define BOOST_GRAPH_PARALLEL_BRANDES_BETWEENNESS_CENTRALITY_HPP
  11
  12 #ifndef BOOST_GRAPH_USE_MPI
  13 #error "Parallel BGL files should not be included unless <boost/graph/use_mpi.hpp> has been included"
  14 #endif
  15
  16 // #define COMPUTE_PATH_COUNTS_INLINE
  17
  18 #include <boost/graph/betweenness_centrality.hpp>
  19 #include <boost/graph/overloading.hpp>
  20 #include <boost/graph/distributed/concepts.hpp>
  21 #include <boost/graph/graph_traits.hpp>
  22 #include <boost/config.hpp>
  23 #include <boost/assert.hpp>
  24
  25 // For additive_reducer
  26 #include <boost/graph/distributed/distributed_graph_utility.hpp>
  27 #include <boost/type_traits/is_convertible.hpp>
  28 #include <boost/type_traits/is_same.hpp>
  29 #include <boost/property_map/property_map.hpp>
  30 #include <boost/graph/named_function_params.hpp>
  31
  32 #include <boost/property_map/parallel/distributed_property_map.hpp>
  33 #include <boost/graph/distributed/detail/dijkstra_shortest_paths.hpp>
  34 #include <boost/tuple/tuple.hpp>
  35
  36 // NGE - Needed for minstd_rand at L807, should pass vertex list
  37 //       or generator instead
  38 #include <boost/random/linear_congruential.hpp>
  39
  40 #include <algorithm>
  41 #include <stack>
  42 #include <vector>
  43
  44 // Appending reducer
  45 template <typename T>
  46 struct append_reducer {
  47   BOOST_STATIC_CONSTANT(bool, non_default_resolver = true);
  48
  49   template<typename K>
  50   T operator()(const K&) const { return T(); }
  51
  52   template<typename K>
  53   T operator()(const K&, const T& x, const T& y) const
  54   {
  55     T z(x.begin(), x.end());
  56     for (typename T::const_iterator iter = y.begin(); iter != y.end(); ++iter)
  57       if (std::find(z.begin(), z.end(), *iter) == z.end())
  58         z.push_back(*iter);
  59
  60     return z;
  61   }
  62 };
  63
  64 namespace boost {
  65
  66   namespace serialization {
  67
  68     // TODO(nge): Write generalized serialization for tuples
  69     template<typename Archive, typename T1, typename T2, typename T3,
  70              typename T4>
  71     void serialize(Archive & ar,
  72                    boost::tuple<T1,T2,T3, T4>& t,
  73                    const unsigned int)
  74     {
  75       ar & boost::tuples::get<0>(t);
  76       ar & boost::tuples::get<1>(t);
  77       ar & boost::tuples::get<2>(t);
  78       ar & boost::tuples::get<3>(t);
  79     }
  80
  81   } // serialization
  82
  83   template <typename OwnerMap, typename Tuple>
  84   class get_owner_of_first_tuple_element {
  85
  86   public:
  87     typedef typename property_traits<OwnerMap>::value_type owner_type;
  88
  89     get_owner_of_first_tuple_element(OwnerMap owner) : owner(owner) { }
  90
  91     owner_type get_owner(Tuple t) { return get(owner, boost::tuples::get<0>(t)); }
  92
  93   private:
  94     OwnerMap owner;
  95   };
  96
  97   template <typename OwnerMap, typename Tuple>
  98   typename get_owner_of_first_tuple_element<OwnerMap, Tuple>::owner_type
  99   get(get_owner_of_first_tuple_element<OwnerMap, Tuple> o, Tuple t)
 100   { return o.get_owner(t); }
 101
 102   template <typename OwnerMap>
 103   class get_owner_of_first_pair_element {
 104
 105   public:
 106     typedef typename property_traits<OwnerMap>::value_type owner_type;
 107
 108     get_owner_of_first_pair_element(OwnerMap owner) : owner(owner) { }
 109
 110     template <typename Vertex, typename T>
 111     owner_type get_owner(std::pair<Vertex, T> p) { return get(owner, p.first); }
 112
 113   private:
 114     OwnerMap owner;
 115   };
 116
 117   template <typename OwnerMap, typename Vertex, typename T>
 118   typename get_owner_of_first_pair_element<OwnerMap>::owner_type
 119   get(get_owner_of_first_pair_element<OwnerMap> o, std::pair<Vertex, T> p)
 120   { return o.get_owner(p); }
 121
 122   namespace graph { namespace parallel { namespace detail {
 123
 124   template<typename DistanceMap, typename IncomingMap>
 125   class betweenness_centrality_msg_value
 126   {
 127     typedef typename property_traits<DistanceMap>::value_type distance_type;
 128     typedef typename property_traits<IncomingMap>::value_type incoming_type;
 129     typedef typename incoming_type::value_type incoming_value_type;
 130
 131   public:
 132     typedef std::pair<distance_type, incoming_value_type> type;
 133
 134     static type create(distance_type dist, incoming_value_type source)
 135     { return std::make_pair(dist, source); }
 136   };
 137
 138
 139   /************************************************************************/
 140   /* Delta-stepping Betweenness Centrality                                */
 141   /************************************************************************/
 142
 143   template<typename Graph, typename DistanceMap, typename IncomingMap,
 144            typename EdgeWeightMap, typename PathCountMap
 145 #ifdef COMPUTE_PATH_COUNTS_INLINE
 146            , typename IsSettledMap, typename VertexIndexMap
 147 #endif
 148            >
 149   class betweenness_centrality_delta_stepping_impl {
 150     // Could inherit from delta_stepping_impl to get run() method
 151     // but for the time being it's just reproduced here
 152
 153     typedef typename graph_traits<Graph>::vertex_descriptor Vertex;
 154     typedef typename graph_traits<Graph>::degree_size_type Degree;
 155     typedef typename property_traits<EdgeWeightMap>::value_type Dist;
 156     typedef typename property_traits<IncomingMap>::value_type IncomingType;
 157     typedef typename boost::graph::parallel::process_group_type<Graph>::type
 158       ProcessGroup;
 159
 160     typedef std::list<Vertex> Bucket;
 161     typedef typename Bucket::iterator BucketIterator;
 162     typedef typename std::vector<Bucket*>::size_type BucketIndex;
 163
 164     typedef betweenness_centrality_msg_value<DistanceMap, IncomingMap>
 165       MessageValue;
 166
 167     enum {
 168       // Relax a remote vertex. The message contains a pair<Vertex,
 169       // MessageValue>, the first part of which is the vertex whose
 170       // tentative distance is being relaxed and the second part
 171       // contains either the new distance (if there is no predecessor
 172       // map) or a pair with the distance and predecessor.
 173       msg_relax
 174     };
 175
 176   public:
 177
 178     // Must supply delta, ctor that guesses delta removed
 179     betweenness_centrality_delta_stepping_impl(const Graph& g,
 180                                                DistanceMap distance,
 181                                                IncomingMap incoming,
 182                                                EdgeWeightMap weight,
 183                                                PathCountMap path_count,
 184 #ifdef COMPUTE_PATH_COUNTS_INLINE
 185                                                IsSettledMap is_settled,
 186                                                VertexIndexMap vertex_index,
 187 #endif
 188                                                Dist delta);
 189
 190     void run(Vertex s);
 191
 192   private:
 193     // Relax the edge (u, v), creating a new best path of distance x.
 194     void relax(Vertex u, Vertex v, Dist x);
 195
 196     // Synchronize all of the processes, by receiving all messages that
 197     // have not yet been received.
 198     void synchronize()
 199     {
 200       using boost::parallel::synchronize;
 201       synchronize(pg);
 202     }
 203
 204     // Setup triggers for msg_relax messages
 205     void setup_triggers()
 206     {
 207       using boost::parallel::simple_trigger;
 208       simple_trigger(pg, msg_relax, this,
 209                      &betweenness_centrality_delta_stepping_impl::handle_msg_relax);
 210     }
 211
 212     void handle_msg_relax(int /*source*/, int /*tag*/,
 213                           const std::pair<Vertex, typename MessageValue::type>& data,
 214                           boost::parallel::trigger_receive_context)
 215     { relax(data.second.second, data.first, data.second.first); }
 216
 217     const Graph& g;
 218     IncomingMap incoming;
 219     DistanceMap distance;
 220     EdgeWeightMap weight;
 221     PathCountMap path_count;
 222 #ifdef COMPUTE_PATH_COUNTS_INLINE
 223     IsSettledMap is_settled;
 224     VertexIndexMap vertex_index;
 225 #endif
 226     Dist delta;
 227     ProcessGroup pg;
 228     typename property_map<Graph, vertex_owner_t>::const_type owner;
 229     typename property_map<Graph, vertex_local_t>::const_type local;
 230
 231     // A "property map" that contains the position of each vertex in
 232     // whatever bucket it resides in.
 233     std::vector<BucketIterator> position_in_bucket;
 234
 235     // Bucket data structure. The ith bucket contains all local vertices
 236     // with (tentative) distance in the range [i*delta,
 237     // (i+1)*delta).
 238     std::vector<Bucket*> buckets;
 239
 240     // This "dummy" list is used only so that we can initialize the
 241     // position_in_bucket property map with non-singular iterators. This
 242     // won't matter for most implementations of the C++ Standard
 243     // Library, but it avoids undefined behavior and allows us to run
 244     // with library "debug modes".
 245     std::list<Vertex> dummy_list;
 246
 247     // A "property map" that states which vertices have been deleted
 248     // from the bucket in this iteration.
 249     std::vector<bool> vertex_was_deleted;
 250   };
 251
 252   template<typename Graph, typename DistanceMap, typename IncomingMap,
 253            typename EdgeWeightMap, typename PathCountMap
 254 #ifdef COMPUTE_PATH_COUNTS_INLINE
 255            , typename IsSettledMap, typename VertexIndexMap
 256 #endif
 257            >
 258   betweenness_centrality_delta_stepping_impl<
 259     Graph, DistanceMap, IncomingMap, EdgeWeightMap, PathCountMap
 260 #ifdef COMPUTE_PATH_COUNTS_INLINE
 261            , IsSettledMap, VertexIndexMap
 262 #endif
 263     >::
 264   betweenness_centrality_delta_stepping_impl(const Graph& g,
 265                                              DistanceMap distance,
 266                                              IncomingMap incoming,
 267                                              EdgeWeightMap weight,
 268                                              PathCountMap path_count,
 269 #ifdef COMPUTE_PATH_COUNTS_INLINE
 270                                              IsSettledMap is_settled,
 271                                              VertexIndexMap vertex_index,
 272 #endif
 273                                              Dist delta)
 274     : g(g),
 275       incoming(incoming),
 276       distance(distance),
 277       weight(weight),
 278       path_count(path_count),
 279 #ifdef COMPUTE_PATH_COUNTS_INLINE
 280       is_settled(is_settled),
 281       vertex_index(vertex_index),
 282 #endif
 283       delta(delta),
 284       pg(boost::graph::parallel::process_group_adl(g), boost::parallel::attach_distributed_object()),
 285       owner(get(vertex_owner, g)),
 286       local(get(vertex_local, g))
 287
 288   { setup_triggers(); }
 289
 290   template<typename Graph, typename DistanceMap, typename IncomingMap,
 291            typename EdgeWeightMap, typename PathCountMap
 292 #ifdef COMPUTE_PATH_COUNTS_INLINE
 293            , typename IsSettledMap, typename VertexIndexMap
 294 #endif
 295            >
 296   void
 297   betweenness_centrality_delta_stepping_impl<
 298     Graph, DistanceMap, IncomingMap, EdgeWeightMap, PathCountMap
 299 #ifdef COMPUTE_PATH_COUNTS_INLINE
 300            , IsSettledMap, VertexIndexMap
 301 #endif
 302     >::
 303   run(Vertex s)
 304   {
 305     typedef typename boost::graph::parallel::process_group_type<Graph>::type
 306       process_group_type;
 307     typename process_group_type::process_id_type id = process_id(pg);
 308
 309     Dist inf = (std::numeric_limits<Dist>::max)();
 310
 311     // None of the vertices are stored in the bucket.
 312     position_in_bucket.clear();
 313     position_in_bucket.resize(num_vertices(g), dummy_list.end());
 314
 315     // None of the vertices have been deleted
 316     vertex_was_deleted.clear();
 317     vertex_was_deleted.resize(num_vertices(g), false);
 318
 319     // No path from s to any other vertex, yet
 320     BGL_FORALL_VERTICES_T(v, g, Graph)
 321       put(distance, v, inf);
 322
 323     // The distance to the starting node is zero
 324     if (get(owner, s) == id)
 325       // Put "s" into its bucket (bucket 0)
 326       relax(s, s, 0);
 327     else
 328       // Note that we know the distance to s is zero
 329       cache(distance, s, 0);
 330
 331 #ifdef COMPUTE_PATH_COUNTS_INLINE
 332     // Synchronize here to deliver initial relaxation since we don't
 333     // synchronize at the beginning of the inner loop any more
 334     synchronize();
 335
 336     // Incoming edge count map is an implementation detail and should
 337     // be freed as soon as possible so build it here
 338     typedef typename graph_traits<Graph>::edges_size_type edges_size_type;
 339
 340     std::vector<edges_size_type> incoming_edge_countS(num_vertices(g));
 341     iterator_property_map<typename std::vector<edges_size_type>::iterator, VertexIndexMap>
 342       incoming_edge_count(incoming_edge_countS.begin(), vertex_index);
 343 #endif
 344
 345     BucketIndex max_bucket = (std::numeric_limits<BucketIndex>::max)();
 346     BucketIndex current_bucket = 0;
 347     do {
 348 #ifdef COMPUTE_PATH_COUNTS_INLINE
 349       // We need to clear the outgoing map after every bucket so just build it here
 350       std::vector<IncomingType> outgoingS(num_vertices(g));
 351       IncomingMap outgoing(outgoingS.begin(), vertex_index);
 352
 353       outgoing.set_reduce(append_reducer<IncomingType>());
 354 #else
 355       // Synchronize with all of the other processes.
 356       synchronize();
 357 #endif
 358
 359       // Find the next bucket that has something in it.
 360       while (current_bucket < buckets.size()
 361              && (!buckets[current_bucket] || buckets[current_bucket]->empty()))
 362         ++current_bucket;
 363       if (current_bucket >= buckets.size())
 364         current_bucket = max_bucket;
 365
 366       // Find the smallest bucket (over all processes) that has vertices
 367       // that need to be processed.
 368       using boost::parallel::all_reduce;
 369       using boost::parallel::minimum;
 370       current_bucket = all_reduce(pg, current_bucket, minimum<BucketIndex>());
 371
 372       if (current_bucket == max_bucket)
 373         // There are no non-empty buckets in any process; exit.
 374         break;
 375
 376       // Contains the set of vertices that have been deleted in the
 377       // relaxation of "light" edges. Note that we keep track of which
 378       // vertices were deleted with the property map
 379       // "vertex_was_deleted".
 380       std::vector<Vertex> deleted_vertices;
 381
 382       // Repeatedly relax light edges
 383       bool nonempty_bucket;
 384       do {
 385         // Someone has work to do in this bucket.
 386
 387         if (current_bucket < buckets.size() && buckets[current_bucket]) {
 388           Bucket& bucket = *buckets[current_bucket];
 389           // For each element in the bucket
 390           while (!bucket.empty()) {
 391             Vertex u = bucket.front();
 392
 393             // Remove u from the front of the bucket
 394             bucket.pop_front();
 395
 396             // Insert u into the set of deleted vertices, if it hasn't
 397             // been done already.
 398             if (!vertex_was_deleted[get(local, u)]) {
 399               vertex_was_deleted[get(local, u)] = true;
 400               deleted_vertices.push_back(u);
 401             }
 402
 403             // Relax each light edge.
 404             Dist u_dist = get(distance, u);
 405             BGL_FORALL_OUTEDGES_T(u, e, g, Graph)
 406               if (get(weight, e) <= delta) // light edge
 407                 relax(u, target(e, g), u_dist + get(weight, e));
 408           }
 409         }
 410
 411         // Synchronize with all of the other processes.
 412         synchronize();
 413
 414         // Is the bucket empty now?
 415         nonempty_bucket = (current_bucket < buckets.size()
 416                            && buckets[current_bucket]
 417                            && !buckets[current_bucket]->empty());
 418       } while (all_reduce(pg, nonempty_bucket, std::logical_or<bool>()));
 419
 420       // Relax heavy edges for each of the vertices that we previously
 421       // deleted.
 422       for (typename std::vector<Vertex>::iterator iter = deleted_vertices.begin();
 423            iter != deleted_vertices.end(); ++iter) {
 424         // Relax each heavy edge.
 425         Vertex u = *iter;
 426         Dist u_dist = get(distance, u);
 427         BGL_FORALL_OUTEDGES_T(u, e, g, Graph)
 428           if (get(weight, e) > delta) // heavy edge
 429             relax(u, target(e, g), u_dist + get(weight, e));
 430
 431 #ifdef COMPUTE_PATH_COUNTS_INLINE
 432         // Set outgoing paths
 433         IncomingType in = get(incoming, u);
 434         for (typename IncomingType::iterator pred = in.begin(); pred != in.end(); ++pred)
 435           if (get(owner, *pred) == id) {
 436             IncomingType x = get(outgoing, *pred);
 437             if (std::find(x.begin(), x.end(), u) == x.end())
 438               x.push_back(u);
 439             put(outgoing, *pred, x);
 440           } else {
 441             IncomingType in;
 442             in.push_back(u);
 443             put(outgoing, *pred, in);
 444           }
 445
 446         // Set incoming edge counts
 447         put(incoming_edge_count, u, in.size());
 448 #endif
 449       }
 450
 451 #ifdef COMPUTE_PATH_COUNTS_INLINE
 452       synchronize();  // Deliver heavy edge relaxations and outgoing paths
 453
 454       // Build Queue
 455       typedef typename property_traits<PathCountMap>::value_type PathCountType;
 456       typedef std::pair<Vertex, PathCountType> queue_value_type;
 457       typedef typename property_map<Graph, vertex_owner_t>::const_type OwnerMap;
 458       typedef typename get_owner_of_first_pair_element<OwnerMap> IndirectOwnerMap;
 459
 460       typedef boost::queue<queue_value_type> local_queue_type;
 461       typedef boost::graph::distributed::distributed_queue<process_group_type,
 462                                                            IndirectOwnerMap,
 463                                                            local_queue_type> dist_queue_type;
 464
 465       IndirectOwnerMap indirect_owner(owner);
 466       dist_queue_type Q(pg, indirect_owner);
 467
 468       // Find sources to initialize queue
 469       BGL_FORALL_VERTICES_T(v, g, Graph) {
 470         if (get(is_settled, v) && !(get(outgoing, v).empty())) {
 471           put(incoming_edge_count, v, 1);
 472           Q.push(std::make_pair(v, 0)); // Push this vertex with no additional path count
 473         }
 474       }
 475
 476       // Set path counts for vertices in this bucket
 477       while (!Q.empty()) {
 478         queue_value_type t = Q.top(); Q.pop();
 479         Vertex v = t.first;
 480         PathCountType p = t.second;
 481
 482         put(path_count, v, get(path_count, v) + p);
 483         put(incoming_edge_count, v, get(incoming_edge_count, v) - 1);
 484
 485         if (get(incoming_edge_count, v) == 0) {
 486           IncomingType out = get(outgoing, v);
 487           for (typename IncomingType::iterator iter = out.begin(); iter != out.end(); ++iter)
 488             Q.push(std::make_pair(*iter, get(path_count, v)));
 489         }
 490       }
 491
 492       // Mark the vertices in this bucket settled
 493       for (typename std::vector<Vertex>::iterator iter = deleted_vertices.begin();
 494            iter != deleted_vertices.end(); ++iter)
 495         put(is_settled, *iter, true);
 496
 497       // No need to clear path count map as it is never read/written remotely
 498       // No need to clear outgoing map as it is re-alloced every bucket
 499 #endif
 500
 501       // Go to the next bucket: the current bucket must already be empty.
 502       ++current_bucket;
 503     } while (true);
 504
 505     // Delete all of the buckets.
 506     for (typename std::vector<Bucket*>::iterator iter = buckets.begin();
 507          iter != buckets.end(); ++iter) {
 508       if (*iter) {
 509         delete *iter;
 510         *iter = 0;
 511       }
 512     }
 513   }
 514
 515   template<typename Graph, typename DistanceMap, typename IncomingMap,
 516            typename EdgeWeightMap, typename PathCountMap
 517 #ifdef COMPUTE_PATH_COUNTS_INLINE
 518            , typename IsSettledMap, typename VertexIndexMap
 519 #endif
 520            >
 521   void
 522   betweenness_centrality_delta_stepping_impl<
 523     Graph, DistanceMap, IncomingMap, EdgeWeightMap, PathCountMap
 524 #ifdef COMPUTE_PATH_COUNTS_INLINE
 525            , IsSettledMap, VertexIndexMap
 526 #endif
 527     >::
 528   relax(Vertex u, Vertex v, Dist x)
 529   {
 530
 531     if (x <= get(distance, v)) {
 532
 533       // We're relaxing the edge to vertex v.
 534       if (get(owner, v) == process_id(pg)) {
 535         if (x < get(distance, v)) {
 536           // Compute the new bucket index for v
 537           BucketIndex new_index = static_cast<BucketIndex>(x / delta);
 538
 539           // Make sure there is enough room in the buckets data structure.
 540           if (new_index >= buckets.size()) buckets.resize(new_index + 1, 0);
 541
 542           // Make sure that we have allocated the bucket itself.
 543           if (!buckets[new_index]) buckets[new_index] = new Bucket;
 544
 545           if (get(distance, v) != (std::numeric_limits<Dist>::max)()
 546               && !vertex_was_deleted[get(local, v)]) {
 547             // We're moving v from an old bucket into a new one. Compute
 548             // the old index, then splice it in.
 549             BucketIndex old_index
 550               = static_cast<BucketIndex>(get(distance, v) / delta);
 551             buckets[new_index]->splice(buckets[new_index]->end(),
 552                                        *buckets[old_index],
 553                                        position_in_bucket[get(local, v)]);
 554           } else {
 555             // We're inserting v into a bucket for the first time. Put it
 556             // at the end.
 557             buckets[new_index]->push_back(v);
 558           }
 559
 560           // v is now at the last position in the new bucket
 561           position_in_bucket[get(local, v)] = buckets[new_index]->end();
 562           --position_in_bucket[get(local, v)];
 563
 564           // Update tentative distance information and incoming, path_count
 565           if (u != v) put(incoming, v, IncomingType(1, u));
 566           put(distance, v, x);
 567         }        // u != v covers initial source relaxation and self-loops
 568         else if (x == get(distance, v) && u != v) {
 569
 570           // Add incoming edge if it's not already in the list
 571           IncomingType in = get(incoming, v);
 572           if (std::find(in.begin(), in.end(), u) == in.end()) {
 573             in.push_back(u);
 574             put(incoming, v, in);
 575           }
 576         }
 577       } else {
 578         // The vertex is remote: send a request to the vertex's owner
 579         send(pg, get(owner, v), msg_relax,
 580              std::make_pair(v, MessageValue::create(x, u)));
 581
 582         // Cache tentative distance information
 583         cache(distance, v, x);
 584       }
 585     }
 586   }
 587
 588   /************************************************************************/
 589   /* Shortest Paths function object for betweenness centrality            */
 590   /************************************************************************/
 591
 592   template<typename WeightMap>
 593   struct brandes_shortest_paths {
 594     typedef typename property_traits<WeightMap>::value_type weight_type;
 595
 596     brandes_shortest_paths()
 597       : weight(1), delta(0)  { }
 598     brandes_shortest_paths(weight_type delta)
 599       : weight(1), delta(delta)  { }
 600     brandes_shortest_paths(WeightMap w)
 601       : weight(w), delta(0)  { }
 602     brandes_shortest_paths(WeightMap w, weight_type delta)
 603       : weight(w), delta(delta)  { }
 604
 605     template<typename Graph, typename IncomingMap, typename DistanceMap,
 606              typename PathCountMap
 607 #ifdef COMPUTE_PATH_COUNTS_INLINE
 608              , typename IsSettledMap, typename VertexIndexMap
 609 #endif
 610
 611              >
 612     void
 613     operator()(Graph& g,
 614                typename graph_traits<Graph>::vertex_descriptor s,
 615                IncomingMap incoming,
 616                DistanceMap distance,
 617                PathCountMap path_count
 618 #ifdef COMPUTE_PATH_COUNTS_INLINE
 619                , IsSettledMap is_settled,
 620                VertexIndexMap vertex_index
 621 #endif
 622                )
 623     {
 624       // The "distance" map needs to act like one, retrieving the default
 625       // value of infinity.
 626       set_property_map_role(vertex_distance, distance);
 627
 628       // Only calculate delta the first time operator() is called
 629       // This presumes g is the same every time, but so does the fact
 630       // that we're reusing the weight map
 631       if (delta == 0)
 632         set_delta(g);
 633
 634       // TODO (NGE): Restructure the code so we don't have to construct
 635       //             impl every time?
 636       betweenness_centrality_delta_stepping_impl<
 637           Graph, DistanceMap, IncomingMap, WeightMap, PathCountMap
 638 #ifdef COMPUTE_PATH_COUNTS_INLINE
 639           , IsSettledMap, VertexIndexMap
 640 #endif
 641             >
 642         impl(g, distance, incoming, weight, path_count,
 643 #ifdef COMPUTE_PATH_COUNTS_INLINE
 644              is_settled, vertex_index,
 645 #endif
 646              delta);
 647
 648       impl.run(s);
 649     }
 650
 651   private:
 652
 653     template <typename Graph>
 654     void
 655     set_delta(const Graph& g)
 656     {
 657       using boost::parallel::all_reduce;
 658       using boost::parallel::maximum;
 659       using std::max;
 660
 661       typedef typename graph_traits<Graph>::degree_size_type Degree;
 662       typedef weight_type Dist;
 663
 664       // Compute the maximum edge weight and degree
 665       Dist max_edge_weight = 0;
 666       Degree max_degree = 0;
 667       BGL_FORALL_VERTICES_T(u, g, Graph) {
 668         max_degree = max BOOST_PREVENT_MACRO_SUBSTITUTION (max_degree, out_degree(u, g));
 669         BGL_FORALL_OUTEDGES_T(u, e, g, Graph)
 670           max_edge_weight = max BOOST_PREVENT_MACRO_SUBSTITUTION (max_edge_weight, get(weight, e));
 671       }
 672
 673       max_edge_weight = all_reduce(process_group(g), max_edge_weight, maximum<Dist>());
 674       max_degree = all_reduce(process_group(g), max_degree, maximum<Degree>());
 675
 676       // Take a guess at delta, based on what works well for random
 677       // graphs.
 678       delta = max_edge_weight / max_degree;
 679       if (delta == 0)
 680         delta = 1;
 681     }
 682
 683     WeightMap     weight;
 684     weight_type   delta;
 685   };
 686
 687   // Perform a single SSSP from the specified vertex and update the centrality map(s)
 688   template<typename Graph, typename CentralityMap, typename EdgeCentralityMap,
 689            typename IncomingMap, typename DistanceMap, typename DependencyMap,
 690            typename PathCountMap,
 691 #ifdef COMPUTE_PATH_COUNTS_INLINE
 692            typename IsSettledMap,
 693 #endif
 694            typename VertexIndexMap, typename ShortestPaths>
 695   void
 696   do_brandes_sssp(const Graph& g,
 697                   CentralityMap centrality,
 698                   EdgeCentralityMap edge_centrality_map,
 699                   IncomingMap incoming,
 700                   DistanceMap distance,
 701                   DependencyMap dependency,
 702                   PathCountMap path_count,
 703 #ifdef COMPUTE_PATH_COUNTS_INLINE
 704                   IsSettledMap is_settled,
 705 #endif
 706                   VertexIndexMap vertex_index,
 707                   ShortestPaths shortest_paths,
 708                   typename graph_traits<Graph>::vertex_descriptor s)
 709   {
 710     using boost::detail::graph::update_centrality;
 711     using boost::graph::parallel::process_group;
 712
 713     typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
 714     typedef typename graph_traits<Graph>::edges_size_type edges_size_type;
 715
 716     typedef typename property_traits<IncomingMap>::value_type incoming_type;
 717     typedef typename property_traits<DistanceMap>::value_type distance_type;
 718     typedef typename property_traits<DependencyMap>::value_type dependency_type;
 719     typedef typename property_traits<PathCountMap>::value_type path_count_type;
 720
 721     typedef typename incoming_type::iterator incoming_iterator;
 722
 723     typedef typename property_map<Graph, vertex_owner_t>::const_type OwnerMap;
 724     OwnerMap owner = get(vertex_owner, g);
 725
 726     typedef typename boost::graph::parallel::process_group_type<Graph>::type
 727       process_group_type;
 728     process_group_type pg = process_group(g);
 729     typename process_group_type::process_id_type id = process_id(pg);
 730
 731     // TODO: Is it faster not to clear some of these maps?
 732     // Initialize for this iteration
 733     distance.clear();
 734     incoming.clear();
 735     path_count.clear();
 736     dependency.clear();
 737     BGL_FORALL_VERTICES_T(v, g, Graph) {
 738       put(path_count, v, 0);
 739       put(dependency, v, 0);
 740     }
 741
 742     if (get(owner, s) == id) {
 743       put(incoming, s, incoming_type());
 744 #ifdef COMPUTE_PATH_COUNTS_INLINE
 745       put(path_count, s, 1);
 746       put(is_settled, s, true);
 747 #endif
 748     }
 749
 750     // Execute the shortest paths algorithm. This will be either
 751     // a weighted or unweighted customized breadth-first search,
 752     shortest_paths(g, s, incoming, distance, path_count
 753 #ifdef COMPUTE_PATH_COUNTS_INLINE
 754                    , is_settled, vertex_index
 755 #endif
 756                    );
 757
 758 #ifndef COMPUTE_PATH_COUNTS_INLINE
 759
 760     //
 761     // TODO: Optimize case where source has no out-edges
 762     //
 763
 764     // Count of incoming edges to tell when all incoming edges have been relaxed in
 765     // the induced shortest paths DAG
 766     std::vector<edges_size_type> incoming_edge_countS(num_vertices(g));
 767     iterator_property_map<typename std::vector<edges_size_type>::iterator, VertexIndexMap>
 768       incoming_edge_count(incoming_edge_countS.begin(), vertex_index);
 769
 770     BGL_FORALL_VERTICES_T(v, g, Graph) {
 771       put(incoming_edge_count, v, get(incoming, v).size());
 772     }
 773
 774     if (get(owner, s) == id) {
 775       put(incoming_edge_count, s, 1);
 776       put(incoming, s, incoming_type());
 777     }
 778
 779     std::vector<incoming_type> outgoingS(num_vertices(g));
 780     iterator_property_map<typename std::vector<incoming_type>::iterator, VertexIndexMap>
 781       outgoing(outgoingS.begin(), vertex_index);
 782
 783     outgoing.set_reduce(append_reducer<incoming_type>());
 784
 785     // Mark forward adjacencies in DAG of shortest paths
 786
 787     // TODO: It's possible to do this using edge flags but it's not currently done this way
 788     //       because during traversal of the DAG we would have to examine all out edges
 789     //       which would lead to more memory accesses and a larger cache footprint.
 790     //
 791     //       In the bidirectional graph case edge flags would be an excellent way of marking
 792     //       edges in the DAG of shortest paths
 793     BGL_FORALL_VERTICES_T(v, g, Graph) {
 794       incoming_type i = get(incoming, v);
 795       for (typename incoming_type::iterator iter = i.begin(); iter != i.end(); ++iter) {
 796         if (get(owner, *iter) == id) {
 797           incoming_type x = get(outgoing, *iter);
 798           if (std::find(x.begin(), x.end(), v) == x.end())
 799             x.push_back(v);
 800           put(outgoing, *iter, x);
 801         } else {
 802           incoming_type in;
 803           in.push_back(v);
 804           put(outgoing, *iter, in);
 805         }
 806       }
 807     }
 808
 809     synchronize(pg);
 810
 811     // Traverse DAG induced by forward edges in dependency order and compute path counts
 812     {
 813       typedef std::pair<vertex_descriptor, path_count_type> queue_value_type;
 814       typedef get_owner_of_first_pair_element<OwnerMap> IndirectOwnerMap;
 815
 816       typedef boost::queue<queue_value_type> local_queue_type;
 817       typedef boost::graph::distributed::distributed_queue<process_group_type,
 818                                                            IndirectOwnerMap,
 819                                                            local_queue_type> dist_queue_type;
 820
 821       IndirectOwnerMap indirect_owner(owner);
 822       dist_queue_type Q(pg, indirect_owner);
 823
 824       if (get(owner, s) == id)
 825         Q.push(std::make_pair(s, 1));
 826
 827       while (!Q.empty()) {
 828         queue_value_type t = Q.top(); Q.pop();
 829         vertex_descriptor v = t.first;
 830         path_count_type p = t.second;
 831
 832         put(path_count, v, get(path_count, v) + p);
 833         put(incoming_edge_count, v, get(incoming_edge_count, v) - 1);
 834
 835         if (get(incoming_edge_count, v) == 0) {
 836           incoming_type out = get(outgoing, v);
 837           for (typename incoming_type::iterator iter = out.begin(); iter != out.end(); ++iter)
 838             Q.push(std::make_pair(*iter, get(path_count, v)));
 839         }
 840       }
 841     }
 842
 843 #endif // COMPUTE_PATH_COUNTS_INLINE
 844
 845     //
 846     // Compute dependencies
 847     //
 848
 849
 850     // Build the distributed_queue
 851     // Value type consists of 1) target of update 2) source of update
 852     // 3) dependency of source 4) path count of source
 853     typedef boost::tuple<vertex_descriptor, vertex_descriptor, dependency_type, path_count_type>
 854       queue_value_type;
 855     typedef get_owner_of_first_tuple_element<OwnerMap, queue_value_type> IndirectOwnerMap;
 856
 857     typedef boost::queue<queue_value_type> local_queue_type;
 858     typedef boost::graph::distributed::distributed_queue<process_group_type,
 859                                                          IndirectOwnerMap,
 860                                                          local_queue_type> dist_queue_type;
 861
 862     IndirectOwnerMap indirect_owner(owner);
 863     dist_queue_type Q(pg, indirect_owner);
 864
 865     // Calculate number of vertices each vertex depends on, when a vertex has been pushed
 866     // that number of times then we will update it
 867     // AND Request path counts of sources of incoming edges
 868     std::vector<dependency_type> dependency_countS(num_vertices(g), 0);
 869     iterator_property_map<typename std::vector<dependency_type>::iterator, VertexIndexMap>
 870       dependency_count(dependency_countS.begin(), vertex_index);
 871
 872     dependency_count.set_reduce(boost::graph::distributed::additive_reducer<dependency_type>());
 873
 874     path_count.set_max_ghost_cells(0);
 875
 876     BGL_FORALL_VERTICES_T(v, g, Graph) {
 877       if (get(distance, v) < (std::numeric_limits<distance_type>::max)()) {
 878         incoming_type el = get(incoming, v);
 879         for (incoming_iterator vw = el.begin(); vw != el.end(); ++vw) {
 880           if (get(owner, *vw) == id)
 881             put(dependency_count, *vw, get(dependency_count, *vw) + 1);
 882           else {
 883             put(dependency_count, *vw, 1);
 884
 885             // Request path counts
 886             get(path_count, *vw);
 887           }
 888
 889           // request() doesn't work here, perhaps because we don't have a copy of this
 890           // ghost cell already?
 891         }
 892       }
 893     }
 894
 895     synchronize(pg);
 896
 897     // Push vertices with non-zero distance/path count and zero dependency count
 898     BGL_FORALL_VERTICES_T(v, g, Graph) {
 899       if (get(distance, v) < (std::numeric_limits<distance_type>::max)()
 900           && get(dependency_count, v) == 0)
 901         Q.push(boost::make_tuple(v, v, get(dependency, v), get(path_count, v)));
 902     }
 903
 904     dependency.set_max_ghost_cells(0);
 905     while(!Q.empty()) {
 906
 907       queue_value_type x = Q.top(); Q.pop();
 908       vertex_descriptor w = boost::tuples::get<0>(x);
 909       vertex_descriptor source = boost::tuples::get<1>(x);
 910       dependency_type dep = boost::tuples::get<2>(x);
 911       path_count_type pc = boost::tuples::get<3>(x);
 912
 913       cache(dependency, source, dep);
 914       cache(path_count, source, pc);
 915
 916       if (get(dependency_count, w) != 0)
 917         put(dependency_count, w, get(dependency_count, w) - 1);
 918
 919       if (get(dependency_count, w) == 0) {
 920
 921         // Update dependency and centrality of sources of incoming edges
 922         incoming_type el = get(incoming, w);
 923         for (incoming_iterator vw = el.begin(); vw != el.end(); ++vw) {
 924           vertex_descriptor v = *vw;
 925
 926           BOOST_ASSERT(get(path_count, w) != 0);
 927
 928           dependency_type factor = dependency_type(get(path_count, v))
 929             / dependency_type(get(path_count, w));
 930           factor *= (dependency_type(1) + get(dependency, w));
 931
 932           if (get(owner, v) == id)
 933             put(dependency, v, get(dependency, v) + factor);
 934           else
 935             put(dependency, v, factor);
 936
 937           update_centrality(edge_centrality_map, v, factor);
 938         }
 939
 940         if (w != s)
 941           update_centrality(centrality, w, get(dependency, w));
 942
 943         // Push sources of edges in incoming edge list
 944         for (incoming_iterator vw = el.begin(); vw != el.end(); ++vw)
 945           Q.push(boost::make_tuple(*vw, w, get(dependency, w), get(path_count, w)));
 946       }
 947     }
 948   }
 949
 950   template<typename Graph, typename CentralityMap, typename EdgeCentralityMap,
 951            typename IncomingMap, typename DistanceMap, typename DependencyMap,
 952            typename PathCountMap, typename VertexIndexMap, typename ShortestPaths,
 953            typename Buffer>
 954   void
 955   brandes_betweenness_centrality_impl(const Graph& g,
 956                                       CentralityMap centrality,
 957                                       EdgeCentralityMap edge_centrality_map,
 958                                       IncomingMap incoming,
 959                                       DistanceMap distance,
 960                                       DependencyMap dependency,
 961                                       PathCountMap path_count,
 962                                       VertexIndexMap vertex_index,
 963                                       ShortestPaths shortest_paths,
 964                                       Buffer sources)
 965   {
 966     using boost::detail::graph::init_centrality_map;
 967     using boost::detail::graph::divide_centrality_by_two;
 968     using boost::graph::parallel::process_group;
 969
 970     typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
 971
 972     typedef typename property_traits<DistanceMap>::value_type distance_type;
 973     typedef typename property_traits<DependencyMap>::value_type dependency_type;
 974
 975     // Initialize centrality
 976     init_centrality_map(vertices(g), centrality);
 977     init_centrality_map(edges(g), edge_centrality_map);
 978
 979     // Set the reduction operation on the dependency map to be addition
 980     dependency.set_reduce(boost::graph::distributed::additive_reducer<dependency_type>());
 981     distance.set_reduce(boost::graph::distributed::choose_min_reducer<distance_type>());
 982
 983     // Don't allow remote procs to write incoming or path_count maps
 984     // updating them is handled inside the betweenness_centrality_queue
 985     incoming.set_consistency_model(0);
 986     path_count.set_consistency_model(0);
 987
 988     typedef typename boost::graph::parallel::process_group_type<Graph>::type
 989       process_group_type;
 990     process_group_type pg = process_group(g);
 991
 992 #ifdef COMPUTE_PATH_COUNTS_INLINE
 993     // Build is_settled maps
 994     std::vector<bool> is_settledS(num_vertices(g));
 995     typedef iterator_property_map<std::vector<bool>::iterator, VertexIndexMap>
 996       IsSettledMap;
 997
 998     IsSettledMap is_settled(is_settledS.begin(), vertex_index);
 999 #endif
1000
1001     if (!sources.empty()) {
1002       // DO SSSPs
1003       while (!sources.empty()) {
1004         do_brandes_sssp(g, centrality, edge_centrality_map, incoming, distance,
1005                         dependency, path_count,
1006 #ifdef COMPUTE_PATH_COUNTS_INLINE
1007                         is_settled,
1008 #endif
1009                         vertex_index, shortest_paths, sources.top());
1010         sources.pop();
1011       }
1012     } else { // Exact Betweenness Centrality
1013       typedef typename graph_traits<Graph>::vertices_size_type vertices_size_type;
1014       vertices_size_type n = num_vertices(g);
1015       n = boost::parallel::all_reduce(pg, n, std::plus<vertices_size_type>());
1016
1017       for (vertices_size_type i = 0; i < n; ++i) {
1018         vertex_descriptor v = vertex(i, g);
1019
1020         do_brandes_sssp(g, centrality, edge_centrality_map, incoming, distance,
1021                         dependency, path_count,
1022 #ifdef COMPUTE_PATH_COUNTS_INLINE
1023                         is_settled,
1024 #endif
1025                         vertex_index, shortest_paths, v);
1026       }
1027     }
1028
1029     typedef typename graph_traits<Graph>::directed_category directed_category;
1030     const bool is_undirected =
1031       is_convertible<directed_category*, undirected_tag*>::value;
1032     if (is_undirected) {
1033       divide_centrality_by_two(vertices(g), centrality);
1034       divide_centrality_by_two(edges(g), edge_centrality_map);
1035     }
1036   }
1037
1038   template<typename Graph, typename CentralityMap, typename EdgeCentralityMap,
1039            typename IncomingMap, typename DistanceMap, typename DependencyMap,
1040            typename PathCountMap, typename VertexIndexMap, typename ShortestPaths,
1041            typename Stack>
1042   void
1043   do_sequential_brandes_sssp(const Graph& g,
1044                              CentralityMap centrality,
1045                              EdgeCentralityMap edge_centrality_map,
1046                              IncomingMap incoming,
1047                              DistanceMap distance,
1048                              DependencyMap dependency,
1049                              PathCountMap path_count,
1050                              VertexIndexMap vertex_index,
1051                              ShortestPaths shortest_paths,
1052                              Stack& ordered_vertices,
1053                              typename graph_traits<Graph>::vertex_descriptor v)
1054   {
1055     using boost::detail::graph::update_centrality;
1056
1057     typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
1058
1059     // Initialize for this iteration
1060     BGL_FORALL_VERTICES_T(w, g, Graph) {
1061       // put(path_count, w, 0);
1062       incoming[w].clear();
1063       put(dependency, w, 0);
1064     }
1065
1066     put(path_count, v, 1);
1067     incoming[v].clear();
1068
1069     // Execute the shortest paths algorithm. This will be either
1070     // Dijkstra's algorithm or a customized breadth-first search,
1071     // depending on whether the graph is weighted or unweighted.
1072     shortest_paths(g, v, ordered_vertices, incoming, distance,
1073                    path_count, vertex_index);
1074
1075     while (!ordered_vertices.empty()) {
1076       vertex_descriptor w = ordered_vertices.top();
1077       ordered_vertices.pop();
1078
1079       typedef typename property_traits<IncomingMap>::value_type
1080             incoming_type;
1081       typedef typename incoming_type::iterator incoming_iterator;
1082       typedef typename property_traits<DependencyMap>::value_type
1083         dependency_type;
1084
1085       for (incoming_iterator vw = incoming[w].begin();
1086            vw != incoming[w].end(); ++vw) {
1087         vertex_descriptor v = source(*vw, g);
1088         dependency_type factor = dependency_type(get(path_count, v))
1089           / dependency_type(get(path_count, w));
1090         factor *= (dependency_type(1) + get(dependency, w));
1091         put(dependency, v, get(dependency, v) + factor);
1092         update_centrality(edge_centrality_map, *vw, factor);
1093       }
1094
1095       if (w != v) {
1096         update_centrality(centrality, w, get(dependency, w));
1097       }
1098     }
1099   }
1100
1101   // Betweenness Centrality variant that duplicates graph across processors
1102   // and parallizes SSSPs
1103   // This function expects a non-distributed graph and property-maps
1104   template<typename ProcessGroup, typename Graph,
1105            typename CentralityMap, typename EdgeCentralityMap,
1106            typename IncomingMap, typename DistanceMap,
1107            typename DependencyMap, typename PathCountMap,
1108            typename VertexIndexMap, typename ShortestPaths,
1109            typename Buffer>
1110   void
1111   non_distributed_brandes_betweenness_centrality_impl(const ProcessGroup& pg,
1112                                                       const Graph& g,
1113                                                       CentralityMap centrality,
1114                                                       EdgeCentralityMap edge_centrality_map,
1115                                                       IncomingMap incoming, // P
1116                                                       DistanceMap distance,         // d
1117                                                       DependencyMap dependency,     // delta
1118                                                       PathCountMap path_count,      // sigma
1119                                                       VertexIndexMap vertex_index,
1120                                                       ShortestPaths shortest_paths,
1121                                                       Buffer sources)
1122   {
1123     using boost::detail::graph::init_centrality_map;
1124     using boost::detail::graph::divide_centrality_by_two;
1125     using boost::graph::parallel::process_group;
1126
1127     typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
1128
1129     typedef ProcessGroup process_group_type;
1130
1131     typename process_group_type::process_id_type id = process_id(pg);
1132     typename process_group_type::process_size_type p = num_processes(pg);
1133
1134     // Initialize centrality
1135     init_centrality_map(vertices(g), centrality);
1136     init_centrality_map(edges(g), edge_centrality_map);
1137
1138     std::stack<vertex_descriptor> ordered_vertices;
1139
1140     if (!sources.empty()) {
1141       std::vector<vertex_descriptor> local_sources;
1142
1143       for (int i = 0; i < id; ++i) if (!sources.empty()) sources.pop();
1144       while (!sources.empty()) {
1145         local_sources.push_back(sources.top());
1146
1147         for (int i = 0; i < p; ++i) if (!sources.empty()) sources.pop();
1148       }
1149
1150       // DO SSSPs
1151       for(size_t i = 0; i < local_sources.size(); ++i)
1152         do_sequential_brandes_sssp(g, centrality, edge_centrality_map, incoming,
1153                                    distance, dependency, path_count, vertex_index,
1154                                    shortest_paths, ordered_vertices, local_sources[i]);
1155
1156     } else { // Exact Betweenness Centrality
1157       typedef typename graph_traits<Graph>::vertices_size_type vertices_size_type;
1158       vertices_size_type n = num_vertices(g);
1159
1160       for (vertices_size_type i = id; i < n; i += p) {
1161         vertex_descriptor v = vertex(i, g);
1162
1163         do_sequential_brandes_sssp(g, centrality, edge_centrality_map, incoming,
1164                                    distance, dependency, path_count, vertex_index,
1165                                    shortest_paths, ordered_vertices, v);
1166       }
1167     }
1168
1169     typedef typename graph_traits<Graph>::directed_category directed_category;
1170     const bool is_undirected =
1171       is_convertible<directed_category*, undirected_tag*>::value;
1172     if (is_undirected) {
1173       divide_centrality_by_two(vertices(g), centrality);
1174       divide_centrality_by_two(edges(g), edge_centrality_map);
1175     }
1176
1177     // Merge the centrality maps by summing the values at each vertex)
1178     // TODO(nge): this copy-out, reduce, copy-in is lame
1179     typedef typename property_traits<CentralityMap>::value_type centrality_type;
1180     typedef typename property_traits<EdgeCentralityMap>::value_type edge_centrality_type;
1181
1182     std::vector<centrality_type> centrality_v(num_vertices(g));
1183     std::vector<edge_centrality_type> edge_centrality_v;
1184     edge_centrality_v.reserve(num_edges(g));
1185
1186     BGL_FORALL_VERTICES_T(v, g, Graph) {
1187       centrality_v[get(vertex_index, v)] = get(centrality, v);
1188     }
1189
1190     // Skip when EdgeCentralityMap is a dummy_property_map
1191     if (!is_same<EdgeCentralityMap, dummy_property_map>::value) {
1192       BGL_FORALL_EDGES_T(e, g, Graph) {
1193         edge_centrality_v.push_back(get(edge_centrality_map, e));
1194       }
1195       // NGE: If we trust that the order of elements in the vector isn't changed in the
1196       //      all_reduce below then this method avoids the need for an edge index map
1197     }
1198
1199     using boost::parallel::all_reduce;
1200
1201     all_reduce(pg, &centrality_v[0], &centrality_v[centrality_v.size()],
1202                &centrality_v[0], std::plus<centrality_type>());
1203
1204     if (edge_centrality_v.size())
1205       all_reduce(pg, &edge_centrality_v[0], &edge_centrality_v[edge_centrality_v.size()],
1206                  &edge_centrality_v[0], std::plus<edge_centrality_type>());
1207
1208     BGL_FORALL_VERTICES_T(v, g, Graph) {
1209       put(centrality, v, centrality_v[get(vertex_index, v)]);
1210     }
1211
1212     // Skip when EdgeCentralityMap is a dummy_property_map
1213     if (!is_same<EdgeCentralityMap, dummy_property_map>::value) {
1214       int i = 0;
1215       BGL_FORALL_EDGES_T(e, g, Graph) {
1216         put(edge_centrality_map, e, edge_centrality_v[i]);
1217         ++i;
1218       }
1219     }
1220   }
1221
1222 } } } // end namespace graph::parallel::detail
1223
1224 template<typename Graph, typename CentralityMap, typename EdgeCentralityMap,
1225          typename IncomingMap, typename DistanceMap, typename DependencyMap,
1226          typename PathCountMap, typename VertexIndexMap, typename Buffer>
1227 void
1228 brandes_betweenness_centrality(const Graph& g,
1229                                CentralityMap centrality,
1230                                EdgeCentralityMap edge_centrality_map,
1231                                IncomingMap incoming,
1232                                DistanceMap distance,
1233                                DependencyMap dependency,
1234                                PathCountMap path_count,
1235                                VertexIndexMap vertex_index,
1236                                Buffer sources,
1237                                typename property_traits<DistanceMap>::value_type delta
1238                                BOOST_GRAPH_ENABLE_IF_MODELS_PARM(Graph,distributed_graph_tag))
1239 {
1240   typedef typename property_traits<DistanceMap>::value_type distance_type;
1241   typedef static_property_map<distance_type> WeightMap;
1242
1243   graph::parallel::detail::brandes_shortest_paths<WeightMap>
1244     shortest_paths(delta);
1245
1246   graph::parallel::detail::brandes_betweenness_centrality_impl(g, centrality,
1247                                                                edge_centrality_map,
1248                                                                incoming, distance,
1249                                                                dependency, path_count,
1250                                                                vertex_index,
1251                                                                shortest_paths,
1252                                                                sources);
1253 }
1254
1255 template<typename Graph, typename CentralityMap, typename EdgeCentralityMap,
1256          typename IncomingMap, typename DistanceMap, typename DependencyMap,
1257          typename PathCountMap, typename VertexIndexMap, typename WeightMap,
1258          typename Buffer>
1259 void
1260 brandes_betweenness_centrality(const Graph& g,
1261                                CentralityMap centrality,
1262                                EdgeCentralityMap edge_centrality_map,
1263                                IncomingMap incoming,
1264                                DistanceMap distance,
1265                                DependencyMap dependency,
1266                                PathCountMap path_count,
1267                                VertexIndexMap vertex_index,
1268                                Buffer sources,
1269                                typename property_traits<WeightMap>::value_type delta,
1270                                WeightMap weight_map
1271                                BOOST_GRAPH_ENABLE_IF_MODELS_PARM(Graph,distributed_graph_tag))
1272 {
1273   graph::parallel::detail::brandes_shortest_paths<WeightMap> shortest_paths(weight_map, delta);
1274
1275   graph::parallel::detail::brandes_betweenness_centrality_impl(g, centrality,
1276                                                                edge_centrality_map,
1277                                                                incoming, distance,
1278                                                                dependency, path_count,
1279                                                                vertex_index,
1280                                                                shortest_paths,
1281                                                                sources);
1282 }
1283
1284 namespace graph { namespace parallel { namespace detail {
1285   template<typename Graph, typename CentralityMap, typename EdgeCentralityMap,
1286            typename WeightMap, typename VertexIndexMap, typename Buffer>
1287   void
1288   brandes_betweenness_centrality_dispatch2(const Graph& g,
1289                                            CentralityMap centrality,
1290                                            EdgeCentralityMap edge_centrality_map,
1291                                            WeightMap weight_map,
1292                                            VertexIndexMap vertex_index,
1293                                            Buffer sources,
1294                                            typename property_traits<WeightMap>::value_type delta)
1295   {
1296     typedef typename graph_traits<Graph>::degree_size_type degree_size_type;
1297     typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
1298     typedef typename mpl::if_c<(is_same<CentralityMap,
1299                                         dummy_property_map>::value),
1300                                          EdgeCentralityMap,
1301                                CentralityMap>::type a_centrality_map;
1302     typedef typename property_traits<a_centrality_map>::value_type
1303       centrality_type;
1304
1305     typename graph_traits<Graph>::vertices_size_type V = num_vertices(g);
1306
1307     std::vector<std::vector<vertex_descriptor> > incoming(V);
1308     std::vector<centrality_type> distance(V);
1309     std::vector<centrality_type> dependency(V);
1310     std::vector<degree_size_type> path_count(V);
1311
1312     brandes_betweenness_centrality(
1313       g, centrality, edge_centrality_map,
1314       make_iterator_property_map(incoming.begin(), vertex_index),
1315       make_iterator_property_map(distance.begin(), vertex_index),
1316       make_iterator_property_map(dependency.begin(), vertex_index),
1317       make_iterator_property_map(path_count.begin(), vertex_index),
1318       vertex_index, unwrap_ref(sources), delta,
1319       weight_map);
1320   }
1321
1322   // TODO: Should the type of the distance and dependency map depend on the
1323   //       value type of the centrality map?
1324   template<typename Graph, typename CentralityMap, typename EdgeCentralityMap,
1325            typename VertexIndexMap, typename Buffer>
1326   void
1327   brandes_betweenness_centrality_dispatch2(const Graph& g,
1328                                            CentralityMap centrality,
1329                                            EdgeCentralityMap edge_centrality_map,
1330                                            VertexIndexMap vertex_index,
1331                                            Buffer sources,
1332                                            typename graph_traits<Graph>::edges_size_type delta)
1333   {
1334     typedef typename graph_traits<Graph>::degree_size_type degree_size_type;
1335     typedef typename graph_traits<Graph>::edges_size_type edges_size_type;
1336     typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
1337
1338     typename graph_traits<Graph>::vertices_size_type V = num_vertices(g);
1339
1340     std::vector<std::vector<vertex_descriptor> > incoming(V);
1341     std::vector<edges_size_type> distance(V);
1342     std::vector<edges_size_type> dependency(V);
1343     std::vector<degree_size_type> path_count(V);
1344
1345     brandes_betweenness_centrality(
1346       g, centrality, edge_centrality_map,
1347       make_iterator_property_map(incoming.begin(), vertex_index),
1348       make_iterator_property_map(distance.begin(), vertex_index),
1349       make_iterator_property_map(dependency.begin(), vertex_index),
1350       make_iterator_property_map(path_count.begin(), vertex_index),
1351       vertex_index, unwrap_ref(sources), delta);
1352   }
1353
1354   template<typename WeightMap>
1355   struct brandes_betweenness_centrality_dispatch1
1356   {
1357     template<typename Graph, typename CentralityMap, typename EdgeCentralityMap,
1358              typename VertexIndexMap, typename Buffer>
1359     static void
1360     run(const Graph& g, CentralityMap centrality, EdgeCentralityMap edge_centrality_map,
1361         VertexIndexMap vertex_index, Buffer sources,
1362         typename property_traits<WeightMap>::value_type delta, WeightMap weight_map)
1363     {
1364       boost::graph::parallel::detail::brandes_betweenness_centrality_dispatch2(
1365        g, centrality, edge_centrality_map, weight_map, vertex_index, sources, delta);
1366     }
1367   };
1368
1369   template<>
1370   struct brandes_betweenness_centrality_dispatch1<boost::param_not_found>
1371   {
1372     template<typename Graph, typename CentralityMap, typename EdgeCentralityMap,
1373              typename VertexIndexMap, typename Buffer>
1374     static void
1375     run(const Graph& g, CentralityMap centrality, EdgeCentralityMap edge_centrality_map,
1376         VertexIndexMap vertex_index, Buffer sources,
1377         typename graph_traits<Graph>::edges_size_type delta,
1378         boost::param_not_found)
1379     {
1380       boost::graph::parallel::detail::brandes_betweenness_centrality_dispatch2(
1381        g, centrality, edge_centrality_map, vertex_index, sources, delta);
1382     }
1383   };
1384
1385 } } } // end namespace graph::parallel::detail
1386
1387 template<typename Graph, typename Param, typename Tag, typename Rest>
1388 void
1389 brandes_betweenness_centrality(const Graph& g,
1390                                const bgl_named_params<Param,Tag,Rest>& params
1391                                BOOST_GRAPH_ENABLE_IF_MODELS_PARM(Graph,distributed_graph_tag))
1392 {
1393   typedef bgl_named_params<Param,Tag,Rest> named_params;
1394
1395   typedef queue<typename graph_traits<Graph>::vertex_descriptor> queue_t;
1396   queue_t q;
1397
1398   typedef typename get_param_type<edge_weight_t, named_params>::type ew_param;
1399   typedef typename detail::choose_impl_result<mpl::true_, Graph, ew_param, edge_weight_t>::type ew;
1400   graph::parallel::detail::brandes_betweenness_centrality_dispatch1<ew>::run(
1401     g,
1402     choose_param(get_param(params, vertex_centrality),
1403                  dummy_property_map()),
1404     choose_param(get_param(params, edge_centrality),
1405                  dummy_property_map()),
1406     choose_const_pmap(get_param(params, vertex_index), g, vertex_index),
1407     choose_param(get_param(params, buffer_param_t()), boost::ref(q)),
1408     choose_param(get_param(params, lookahead_t()), 0),
1409     choose_const_pmap(get_param(params, edge_weight), g, edge_weight));
1410 }
1411
1412 template<typename Graph, typename CentralityMap>
1413 void
1414 brandes_betweenness_centrality(const Graph& g, CentralityMap centrality
1415                                BOOST_GRAPH_ENABLE_IF_MODELS_PARM(Graph,distributed_graph_tag))
1416 {
1417   typedef queue<typename graph_traits<Graph>::vertex_descriptor> queue_t;
1418   queue_t q;
1419
1420   boost::graph::parallel::detail::brandes_betweenness_centrality_dispatch2(
1421     g, centrality, dummy_property_map(), get(vertex_index, g), boost::ref(q), 0);
1422 }
1423
1424 template<typename Graph, typename CentralityMap, typename EdgeCentralityMap>
1425 void
1426 brandes_betweenness_centrality(const Graph& g, CentralityMap centrality,
1427                                EdgeCentralityMap edge_centrality_map
1428                                BOOST_GRAPH_ENABLE_IF_MODELS_PARM(Graph,distributed_graph_tag))
1429 {
1430   typedef queue<int> queue_t;
1431   queue_t q;
1432
1433   boost::graph::parallel::detail::brandes_betweenness_centrality_dispatch2(
1434     g, centrality, edge_centrality_map, get(vertex_index, g), boost::ref(q), 0);
1435 }
1436
1437 template<typename ProcessGroup, typename Graph, typename CentralityMap,
1438          typename EdgeCentralityMap, typename IncomingMap, typename DistanceMap,
1439          typename DependencyMap, typename PathCountMap, typename VertexIndexMap,
1440          typename Buffer>
1441 void
1442 non_distributed_brandes_betweenness_centrality(const ProcessGroup& pg,
1443                                                const Graph& g,
1444                                                CentralityMap centrality,
1445                                                EdgeCentralityMap edge_centrality_map,
1446                                                IncomingMap incoming,
1447                                                DistanceMap distance,
1448                                                DependencyMap dependency,
1449                                                PathCountMap path_count,
1450                                                VertexIndexMap vertex_index,
1451                                                Buffer sources)
1452 {
1453   detail::graph::brandes_unweighted_shortest_paths shortest_paths;
1454
1455   graph::parallel::detail::non_distributed_brandes_betweenness_centrality_impl(pg, g, centrality,
1456                                                                                edge_centrality_map,
1457                                                                                incoming, distance,
1458                                                                                dependency, path_count,
1459                                                                                vertex_index,
1460                                                                                shortest_paths,
1461                                                                                sources);
1462 }
1463
1464 template<typename ProcessGroup, typename Graph, typename CentralityMap,
1465          typename EdgeCentralityMap, typename IncomingMap, typename DistanceMap,
1466          typename DependencyMap, typename PathCountMap, typename VertexIndexMap,
1467          typename WeightMap, typename Buffer>
1468 void
1469 non_distributed_brandes_betweenness_centrality(const ProcessGroup& pg,
1470                                                const Graph& g,
1471                                                CentralityMap centrality,
1472                                                EdgeCentralityMap edge_centrality_map,
1473                                                IncomingMap incoming,
1474                                                DistanceMap distance,
1475                                                DependencyMap dependency,
1476                                                PathCountMap path_count,
1477                                                VertexIndexMap vertex_index,
1478                                                WeightMap weight_map,
1479                                                Buffer sources)
1480 {
1481   detail::graph::brandes_dijkstra_shortest_paths<WeightMap> shortest_paths(weight_map);
1482
1483   graph::parallel::detail::non_distributed_brandes_betweenness_centrality_impl(pg, g, centrality,
1484                                                                                edge_centrality_map,
1485                                                                                incoming, distance,
1486                                                                                dependency, path_count,
1487                                                                                vertex_index,
1488                                                                                shortest_paths,
1489                                                                                sources);
1490 }
1491
1492 namespace detail { namespace graph {
1493   template<typename ProcessGroup, typename Graph, typename CentralityMap,
1494            typename EdgeCentralityMap, typename WeightMap, typename VertexIndexMap,
1495            typename Buffer>
1496   void
1497   non_distributed_brandes_betweenness_centrality_dispatch2(const ProcessGroup& pg,
1498                                                            const Graph& g,
1499                                                            CentralityMap centrality,
1500                                                            EdgeCentralityMap edge_centrality_map,
1501                                                            WeightMap weight_map,
1502                                                            VertexIndexMap vertex_index,
1503                                                            Buffer sources)
1504   {
1505     typedef typename graph_traits<Graph>::degree_size_type degree_size_type;
1506     typedef typename graph_traits<Graph>::edge_descriptor edge_descriptor;
1507     typedef typename mpl::if_c<(is_same<CentralityMap,
1508                                         dummy_property_map>::value),
1509                                          EdgeCentralityMap,
1510                                CentralityMap>::type a_centrality_map;
1511     typedef typename property_traits<a_centrality_map>::value_type
1512       centrality_type;
1513
1514     typename graph_traits<Graph>::vertices_size_type V = num_vertices(g);
1515
1516     std::vector<std::vector<edge_descriptor> > incoming(V);
1517     std::vector<centrality_type> distance(V);
1518     std::vector<centrality_type> dependency(V);
1519     std::vector<degree_size_type> path_count(V);
1520
1521     non_distributed_brandes_betweenness_centrality(
1522       pg, g, centrality, edge_centrality_map,
1523       make_iterator_property_map(incoming.begin(), vertex_index),
1524       make_iterator_property_map(distance.begin(), vertex_index),
1525       make_iterator_property_map(dependency.begin(), vertex_index),
1526       make_iterator_property_map(path_count.begin(), vertex_index),
1527       vertex_index, weight_map, unwrap_ref(sources));
1528   }
1529
1530
1531   template<typename ProcessGroup, typename Graph, typename CentralityMap,
1532            typename EdgeCentralityMap, typename VertexIndexMap, typename Buffer>
1533   void
1534   non_distributed_brandes_betweenness_centrality_dispatch2(const ProcessGroup& pg,
1535                                                            const Graph& g,
1536                                                            CentralityMap centrality,
1537                                                            EdgeCentralityMap edge_centrality_map,
1538                                                            VertexIndexMap vertex_index,
1539                                                            Buffer sources)
1540   {
1541     typedef typename graph_traits<Graph>::degree_size_type degree_size_type;
1542     typedef typename graph_traits<Graph>::edge_descriptor edge_descriptor;
1543     typedef typename mpl::if_c<(is_same<CentralityMap,
1544                                         dummy_property_map>::value),
1545                                          EdgeCentralityMap,
1546                                CentralityMap>::type a_centrality_map;
1547     typedef typename property_traits<a_centrality_map>::value_type
1548       centrality_type;
1549
1550     typename graph_traits<Graph>::vertices_size_type V = num_vertices(g);
1551
1552     std::vector<std::vector<edge_descriptor> > incoming(V);
1553     std::vector<centrality_type> distance(V);
1554     std::vector<centrality_type> dependency(V);
1555     std::vector<degree_size_type> path_count(V);
1556
1557     non_distributed_brandes_betweenness_centrality(
1558       pg, g, centrality, edge_centrality_map,
1559       make_iterator_property_map(incoming.begin(), vertex_index),
1560       make_iterator_property_map(distance.begin(), vertex_index),
1561       make_iterator_property_map(dependency.begin(), vertex_index),
1562       make_iterator_property_map(path_count.begin(), vertex_index),
1563       vertex_index, unwrap_ref(sources));
1564   }
1565
1566   template<typename WeightMap>
1567   struct non_distributed_brandes_betweenness_centrality_dispatch1
1568   {
1569     template<typename ProcessGroup, typename Graph, typename CentralityMap,
1570              typename EdgeCentralityMap, typename VertexIndexMap, typename Buffer>
1571     static void
1572     run(const ProcessGroup& pg, const Graph& g, CentralityMap centrality,
1573         EdgeCentralityMap edge_centrality_map, VertexIndexMap vertex_index,
1574         Buffer sources, WeightMap weight_map)
1575     {
1576       non_distributed_brandes_betweenness_centrality_dispatch2(pg, g, centrality, edge_centrality_map,
1577                                                                weight_map, vertex_index, sources);
1578     }
1579   };
1580
1581   template<>
1582   struct non_distributed_brandes_betweenness_centrality_dispatch1<param_not_found>
1583   {
1584     template<typename ProcessGroup, typename Graph, typename CentralityMap,
1585              typename EdgeCentralityMap, typename VertexIndexMap, typename Buffer>
1586     static void
1587     run(const ProcessGroup& pg, const Graph& g, CentralityMap centrality,
1588         EdgeCentralityMap edge_centrality_map, VertexIndexMap vertex_index,
1589         Buffer sources, param_not_found)
1590     {
1591       non_distributed_brandes_betweenness_centrality_dispatch2(pg, g, centrality, edge_centrality_map,
1592                                                                vertex_index, sources);
1593     }
1594   };
1595
1596 } } // end namespace detail::graph
1597
1598 template<typename ProcessGroup, typename Graph, typename Param, typename Tag, typename Rest>
1599 void
1600 non_distributed_brandes_betweenness_centrality(const ProcessGroup& pg, const Graph& g,
1601                                                const bgl_named_params<Param,Tag,Rest>& params)
1602 {
1603   typedef bgl_named_params<Param,Tag,Rest> named_params;
1604
1605   typedef queue<int> queue_t;
1606   queue_t q;
1607
1608   typedef typename get_param_type<edge_weight_t, named_params>::type ew_param;
1609   typedef typename detail::choose_impl_result<mpl::true_, Graph, ew_param, edge_weight_t>::type ew;
1610   detail::graph::non_distributed_brandes_betweenness_centrality_dispatch1<ew>::run(
1611     pg, g,
1612     choose_param(get_param(params, vertex_centrality),
1613                  dummy_property_map()),
1614     choose_param(get_param(params, edge_centrality),
1615                  dummy_property_map()),
1616     choose_const_pmap(get_param(params, vertex_index), g, vertex_index),
1617     choose_param(get_param(params, buffer_param_t()),  boost::ref(q)),
1618     choose_const_pmap(get_param(params, edge_weight), g, edge_weight));
1619 }
1620
1621 template<typename ProcessGroup, typename Graph, typename CentralityMap>
1622 void
1623 non_distributed_brandes_betweenness_centrality(const ProcessGroup& pg, const Graph& g,
1624                                                CentralityMap centrality)
1625 {
1626   typedef queue<int> queue_t;
1627   queue_t q;
1628
1629   detail::graph::non_distributed_brandes_betweenness_centrality_dispatch2(
1630     pg, g, centrality, dummy_property_map(), get(vertex_index, g), boost::ref(q));
1631 }
1632
1633 template<typename ProcessGroup, typename Graph, typename CentralityMap,
1634          typename Buffer>
1635 void
1636 non_distributed_brandes_betweenness_centrality(const ProcessGroup& pg, const Graph& g,
1637                                                CentralityMap centrality, Buffer sources)
1638 {
1639   detail::graph::non_distributed_brandes_betweenness_centrality_dispatch2(
1640     pg, g, centrality, dummy_property_map(), get(vertex_index, g), sources);
1641 }
1642
1643 template<typename ProcessGroup, typename Graph, typename CentralityMap,
1644          typename EdgeCentralityMap, typename Buffer>
1645 void
1646 non_distributed_brandes_betweenness_centrality(const ProcessGroup& pg, const Graph& g,
1647                                                CentralityMap centrality,
1648                                                EdgeCentralityMap edge_centrality_map,
1649                                                Buffer sources)
1650 {
1651   detail::graph::non_distributed_brandes_betweenness_centrality_dispatch2(
1652     pg, g, centrality, edge_centrality_map, get(vertex_index, g), sources);
1653 }
1654
1655 // Compute the central point dominance of a graph.
1656 // TODO: Make sure central point dominance works in parallel case
1657 template<typename Graph, typename CentralityMap>
1658 typename property_traits<CentralityMap>::value_type
1659 central_point_dominance(const Graph& g, CentralityMap centrality
1660                         BOOST_GRAPH_ENABLE_IF_MODELS_PARM(Graph,distributed_graph_tag))
1661 {
1662   using std::max;
1663
1664   typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;
1665   typedef typename property_traits<CentralityMap>::value_type centrality_type;
1666   typedef typename graph_traits<Graph>::vertices_size_type vertices_size_type;
1667
1668   typedef typename boost::graph::parallel::process_group_type<Graph>::type
1669     process_group_type;
1670   process_group_type pg = boost::graph::parallel::process_group(g);
1671
1672   vertices_size_type n = num_vertices(g);
1673
1674   using boost::parallel::all_reduce;
1675   n = all_reduce(pg, n, std::plus<vertices_size_type>());
1676
1677   // Find max centrality
1678   centrality_type max_centrality(0);
1679   vertex_iterator v, v_end;
1680   for (boost::tie(v, v_end) = vertices(g); v != v_end; ++v) {
1681     max_centrality = (max)(max_centrality, get(centrality, *v));
1682   }
1683
1684   // All reduce to get global max centrality
1685   max_centrality = all_reduce(pg, max_centrality, boost::parallel::maximum<centrality_type>());
1686
1687   // Compute central point dominance
1688   centrality_type sum(0);
1689   for (boost::tie(v, v_end) = vertices(g); v != v_end; ++v) {
1690     sum += (max_centrality - get(centrality, *v));
1691   }
1692
1693   sum = all_reduce(pg, sum, std::plus<centrality_type>());
1694
1695   return sum/(n-1);
1696 }
1697
1698 } // end namespace boost
1699
1700 #endif // BOOST_GRAPH_PARALLEL_BRANDES_BETWEENNESS_CENTRALITY_HPP