// Copyright 2004 The Trustees of Indiana University.

// Distributed under the Boost Software License, Version 1.0.

// (See accompanying file LICENSE_1_0.txt or copy at

// http://www.boost.org/LICENSE_1_0.txt)

//  Authors: Douglas Gregor

//           Andrew Lumsdaine

#ifndef BOOST_GRAPH_BRANDES_BETWEENNESS_CENTRALITY_HPP

#define BOOST_GRAPH_BRANDES_BETWEENNESS_CENTRALITY_HPP

#include <stack>

#include <vector>

#include <boost/graph/dijkstra_shortest_paths.hpp>

#include <boost/graph/breadth_first_search.hpp>

#include <boost/graph/relax.hpp>

#include <boost/graph/graph_traits.hpp>

#include <boost/tuple/tuple.hpp>

#include <boost/type_traits/is_convertible.hpp>

#include <boost/type_traits/is_same.hpp>

#include <boost/mpl/if.hpp>

#include <boost/property_map.hpp>

#include <boost/graph/named_function_params.hpp>

#include <algorithm>

namespace boost {

namespace detail { namespace graph {

  /**

   * Customized visitor passed to Dijkstra's algorithm by Brandes'

   * betweenness centrality algorithm. This visitor is responsible for

   * keeping track of the order in which vertices are discovered, the

   * predecessors on the shortest path(s) to a vertex, and the number

   * of shortest paths.

   */

  template<typename Graph, typename WeightMap, typename IncomingMap,

           typename DistanceMap, typename PathCountMap>

  struct brandes_dijkstra_visitor : public bfs_visitor<>

  {

    typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;

    typedef typename graph_traits<Graph>::edge_descriptor edge_descriptor;

    brandes_dijkstra_visitor(std::stack<vertex_descriptor>& ordered_vertices,

                             WeightMap weight,

                             IncomingMap incoming,

                             DistanceMap distance,

                             PathCountMap path_count)

      : ordered_vertices(ordered_vertices), weight(weight), 

        incoming(incoming), distance(distance),

        path_count(path_count)

    { }

    /**

     * Whenever an edge e = (v, w) is relaxed, the incoming edge list

     * for w is set to {(v, w)} and the shortest path count of w is set to

     * the number of paths that reach {v}.

     */

    void edge_relaxed(edge_descriptor e, const Graph& g) 

    { 

      vertex_descriptor v = source(e, g), w = target(e, g);

      incoming[w].clear();

      incoming[w].push_back(e);

      put(path_count, w, get(path_count, v));

    }

    /**

     * If an edge e = (v, w) was not relaxed, it may still be the case

     * that we've found more equally-short paths, so include {(v, w)} in the

     * incoming edges of w and add all of the shortest paths to v to the

     * shortest path count of w.

     */

    void edge_not_relaxed(edge_descriptor e, const Graph& g) 

    {

      typedef typename property_traits<WeightMap>::value_type weight_type;

      typedef typename property_traits<DistanceMap>::value_type distance_type;

      vertex_descriptor v = source(e, g), w = target(e, g);

      distance_type d_v = get(distance, v), d_w = get(distance, w);

      weight_type w_e = get(weight, e);

      closed_plus<distance_type> combine;

      if (d_w == combine(d_v, w_e)) {

        put(path_count, w, get(path_count, w) + get(path_count, v));

        incoming[w].push_back(e);

      }

    }

    /// Keep track of vertices as they are reached

    void examine_vertex(vertex_descriptor w, const Graph&) 

    { 

      ordered_vertices.push(w);

    }

  private:

    std::stack<vertex_descriptor>& ordered_vertices;

    WeightMap weight;

    IncomingMap incoming;

    DistanceMap distance;

    PathCountMap path_count;

  };

  /**

   * Function object that calls Dijkstra's shortest paths algorithm

   * using the Dijkstra visitor for the Brandes betweenness centrality

   * algorithm.

   */

  template<typename WeightMap>

  struct brandes_dijkstra_shortest_paths

  {

    brandes_dijkstra_shortest_paths(WeightMap weight_map) 

      : weight_map(weight_map) { }

    template<typename Graph, typename IncomingMap, typename DistanceMap, 

             typename PathCountMap, typename VertexIndexMap>

    void 

    operator()(Graph& g, 

               typename graph_traits<Graph>::vertex_descriptor s,

               std::stack<typename graph_traits<Graph>::vertex_descriptor>& ov,

               IncomingMap incoming,

               DistanceMap distance,

               PathCountMap path_count,

               VertexIndexMap vertex_index)

    {

      typedef brandes_dijkstra_visitor<Graph, WeightMap, IncomingMap, 

                                       DistanceMap, PathCountMap> visitor_type;

      visitor_type visitor(ov, weight_map, incoming, distance, path_count);

      dijkstra_shortest_paths(g, s, 

                              boost::weight_map(weight_map)

                              .vertex_index_map(vertex_index)

                              .distance_map(distance)

                              .visitor(visitor));

    }

  private:

    WeightMap weight_map;

  };

  /**

   * Function object that invokes breadth-first search for the

   * unweighted form of the Brandes betweenness centrality algorithm.

   */

  struct brandes_unweighted_shortest_paths

  {

    /**

     * Customized visitor passed to breadth-first search, which

     * records predecessor and the number of shortest paths to each

     * vertex.

     */

    template<typename Graph, typename IncomingMap, typename DistanceMap, 

             typename PathCountMap>

    struct visitor_type : public bfs_visitor<>

    {

      typedef typename graph_traits<Graph>::edge_descriptor edge_descriptor;

      typedef typename graph_traits<Graph>::vertex_descriptor 

        vertex_descriptor;

      visitor_type(IncomingMap incoming, DistanceMap distance, 

                   PathCountMap path_count, 

                   std::stack<vertex_descriptor>& ordered_vertices)

        : incoming(incoming), distance(distance), 

          path_count(path_count), ordered_vertices(ordered_vertices) { }

      /// Keep track of vertices as they are reached

      void examine_vertex(vertex_descriptor v, Graph&)

      {

        ordered_vertices.push(v);

      }

      /**

       * Whenever an edge e = (v, w) is labelled a tree edge, the

       * incoming edge list for w is set to {(v, w)} and the shortest

       * path count of w is set to the number of paths that reach {v}.

       */

      void tree_edge(edge_descriptor e, Graph& g)

      {

        vertex_descriptor v = source(e, g);

        vertex_descriptor w = target(e, g);

        put(distance, w, get(distance, v) + 1);

        put(path_count, w, get(path_count, v));

        incoming[w].push_back(e);

      }

      /**

       * If an edge e = (v, w) is not a tree edge, it may still be the

       * case that we've found more equally-short paths, so include (v, w)

       * in the incoming edge list of w and add all of the shortest

       * paths to v to the shortest path count of w.

       */

      void non_tree_edge(edge_descriptor e, Graph& g)

      {

        vertex_descriptor v = source(e, g);

        vertex_descriptor w = target(e, g);

        if (get(distance, w) == get(distance, v) + 1) {

          put(path_count, w, get(path_count, w) + get(path_count, v));

          incoming[w].push_back(e);

        }

      }

    private:

      IncomingMap incoming;

      DistanceMap distance;

      PathCountMap path_count;

      std::stack<vertex_descriptor>& ordered_vertices;

    };

    template<typename Graph, typename IncomingMap, typename DistanceMap, 

             typename PathCountMap, typename VertexIndexMap>

    void 

    operator()(Graph& g, 

               typename graph_traits<Graph>::vertex_descriptor s,

               std::stack<typename graph_traits<Graph>::vertex_descriptor>& ov,

               IncomingMap incoming,

               DistanceMap distance,

               PathCountMap path_count,

               VertexIndexMap vertex_index)

    {

      typedef typename graph_traits<Graph>::vertex_descriptor

        vertex_descriptor;

      visitor_type<Graph, IncomingMap, DistanceMap, PathCountMap>

        visitor(incoming, distance, path_count, ov);

      std::vector<default_color_type> 

        colors(num_vertices(g), color_traits<default_color_type>::white());

      boost::queue<vertex_descriptor> Q;

      breadth_first_visit(g, s, Q, visitor, 

                          make_iterator_property_map(colors.begin(), 

                                                     vertex_index));

    }

  };

  // When the edge centrality map is a dummy property map, no

  // initialization is needed.

  template<typename Iter>

  inline void 

  init_centrality_map(std::pair<Iter, Iter>, dummy_property_map) { }

  // When we have a real edge centrality map, initialize all of the

  // centralities to zero.

  template<typename Iter, typename Centrality>

  void 

  init_centrality_map(std::pair<Iter, Iter> keys, Centrality centrality_map)

  {

    typedef typename property_traits<Centrality>::value_type 

      centrality_type;

    while (keys.first != keys.second) {

      put(centrality_map, *keys.first, centrality_type(0));

      ++keys.first;

    }

  }

  // When the edge centrality map is a dummy property map, no update

  // is performed.

  template<typename Key, typename T>

  inline void 

  update_centrality(dummy_property_map, const Key&, const T&) { }

  // When we have a real edge centrality map, add the value to the map

  template<typename CentralityMap, typename Key, typename T>

  inline void 

  update_centrality(CentralityMap centrality_map, Key k, const T& x)

  { put(centrality_map, k, get(centrality_map, k) + x); }

  template<typename Iter>

  inline void 

  divide_centrality_by_two(std::pair<Iter, Iter>, dummy_property_map) {}

  template<typename Iter, typename CentralityMap>

  inline void

  divide_centrality_by_two(std::pair<Iter, Iter> keys, 

                           CentralityMap centrality_map)

  {

    typename property_traits<CentralityMap>::value_type two(2);

    while (keys.first != keys.second) {

      put(centrality_map, *keys.first, get(centrality_map, *keys.first) / two);

      ++keys.first;

    }

  }

  template<typename Graph, typename CentralityMap, typename EdgeCentralityMap,

           typename IncomingMap, typename DistanceMap, 

           typename DependencyMap, typename PathCountMap,

           typename VertexIndexMap, typename ShortestPaths>

  void 

  brandes_betweenness_centrality_impl(const Graph& g, 

                                      CentralityMap centrality,     // C_B

                                      EdgeCentralityMap edge_centrality_map,

                                      IncomingMap incoming, // P

                                      DistanceMap distance,         // d

                                      DependencyMap dependency,     // delta

                                      PathCountMap path_count,      // sigma

                                      VertexIndexMap vertex_index,

                                      ShortestPaths shortest_paths)

  {

    typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;

    typedef typename graph_traits<Graph>::edge_iterator edge_iterator;

    typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;

    // Initialize centrality

    init_centrality_map(vertices(g), centrality);

    init_centrality_map(edges(g), edge_centrality_map);

    std::stack<vertex_descriptor> ordered_vertices;

    vertex_iterator s, s_end;

    for (tie(s, s_end) = vertices(g); s != s_end; ++s) {

      // Initialize for this iteration

      vertex_iterator w, w_end;

      for (tie(w, w_end) = vertices(g); w != w_end; ++w) {

        incoming[*w].clear();

        put(path_count, *w, 0);

        put(dependency, *w, 0);

      }

      put(path_count, *s, 1);

      // Execute the shortest paths algorithm. This will be either

      // Dijkstra's algorithm or a customized breadth-first search,

      // depending on whether the graph is weighted or unweighted.

      shortest_paths(g, *s, ordered_vertices, incoming, distance,

                     path_count, vertex_index);

      while (!ordered_vertices.empty()) {

        vertex_descriptor w = ordered_vertices.top();

        ordered_vertices.pop();

        typedef typename property_traits<IncomingMap>::value_type

          incoming_type;

        typedef typename incoming_type::iterator incoming_iterator;

        typedef typename property_traits<DependencyMap>::value_type 

          dependency_type;

        for (incoming_iterator vw = incoming[w].begin();

             vw != incoming[w].end(); ++vw) {

          vertex_descriptor v = source(*vw, g);

          dependency_type factor = dependency_type(get(path_count, v))

            / dependency_type(get(path_count, w));

          factor *= (dependency_type(1) + get(dependency, w));

          put(dependency, v, get(dependency, v) + factor);

          update_centrality(edge_centrality_map, *vw, factor);

        }

        if (w != *s) {

          update_centrality(centrality, w, get(dependency, w));

        }

      }

    }

    typedef typename graph_traits<Graph>::directed_category directed_category;

    const bool is_undirected = 

      is_convertible<directed_category*, undirected_tag*>::value;

    if (is_undirected) {

      divide_centrality_by_two(vertices(g), centrality);

      divide_centrality_by_two(edges(g), edge_centrality_map);

    }

  }

} } // end namespace detail::graph

template<typename Graph, typename CentralityMap, typename EdgeCentralityMap,

         typename IncomingMap, typename DistanceMap, 

         typename DependencyMap, typename PathCountMap, 

         typename VertexIndexMap>

void 

brandes_betweenness_centrality(const Graph& g, 

                               CentralityMap centrality,     // C_B

                               EdgeCentralityMap edge_centrality_map,

                               IncomingMap incoming, // P

                               DistanceMap distance,         // d

                               DependencyMap dependency,     // delta

                               PathCountMap path_count,      // sigma

                               VertexIndexMap vertex_index)

{

  detail::graph::brandes_unweighted_shortest_paths shortest_paths;

  detail::graph::brandes_betweenness_centrality_impl(g, centrality, 

                                                     edge_centrality_map,

                                                     incoming, distance,

                                                     dependency, path_count,

                                                     vertex_index, 

                                                     shortest_paths);

}

template<typename Graph, typename CentralityMap, typename EdgeCentralityMap, 

         typename IncomingMap, typename DistanceMap, 

         typename DependencyMap, typename PathCountMap, 

         typename VertexIndexMap, typename WeightMap>    

void 

brandes_betweenness_centrality(const Graph& g, 

                               CentralityMap centrality,     // C_B

                               EdgeCentralityMap edge_centrality_map,

                               IncomingMap incoming, // P

                               DistanceMap distance,         // d

                               DependencyMap dependency,     // delta

                               PathCountMap path_count,      // sigma

                               VertexIndexMap vertex_index,

                               WeightMap weight_map)

{

  detail::graph::brandes_dijkstra_shortest_paths<WeightMap>

    shortest_paths(weight_map);

  detail::graph::brandes_betweenness_centrality_impl(g, centrality, 

                                                     edge_centrality_map,

                                                     incoming, distance,

                                                     dependency, path_count,

                                                     vertex_index, 

                                                     shortest_paths);

}

namespace detail { namespace graph {

  template<typename Graph, typename CentralityMap, typename EdgeCentralityMap,

           typename WeightMap, typename VertexIndexMap>

  void 

  brandes_betweenness_centrality_dispatch2(const Graph& g,

                                           CentralityMap centrality,

                                           EdgeCentralityMap edge_centrality_map,

                                           WeightMap weight_map,

                                           VertexIndexMap vertex_index)

  {

    typedef typename graph_traits<Graph>::degree_size_type degree_size_type;

    typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;

    typedef typename graph_traits<Graph>::edge_descriptor edge_descriptor;

    typedef typename mpl::if_c<(is_same<CentralityMap, 

                                        dummy_property_map>::value),

                                         EdgeCentralityMap, 

                               CentralityMap>::type a_centrality_map;

    typedef typename property_traits<a_centrality_map>::value_type 

      centrality_type;

    typename graph_traits<Graph>::vertices_size_type V = num_vertices(g);

    std::vector<std::vector<edge_descriptor> > incoming(V);

    std::vector<centrality_type> distance(V);

    std::vector<centrality_type> dependency(V);

    std::vector<degree_size_type> path_count(V);

    brandes_betweenness_centrality(

      g, centrality, edge_centrality_map,

      make_iterator_property_map(incoming.begin(), vertex_index),

      make_iterator_property_map(distance.begin(), vertex_index),

      make_iterator_property_map(dependency.begin(), vertex_index),

      make_iterator_property_map(path_count.begin(), vertex_index),

      vertex_index,

      weight_map);

  }

  template<typename Graph, typename CentralityMap, typename EdgeCentralityMap,

           typename VertexIndexMap>

  void 

  brandes_betweenness_centrality_dispatch2(const Graph& g,

                                           CentralityMap centrality,

                                           EdgeCentralityMap edge_centrality_map,

                                           VertexIndexMap vertex_index)

  {

    typedef typename graph_traits<Graph>::degree_size_type degree_size_type;

    typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;

    typedef typename graph_traits<Graph>::edge_descriptor edge_descriptor;

    typedef typename mpl::if_c<(is_same<CentralityMap, 

                                        dummy_property_map>::value),

                                         EdgeCentralityMap, 

                               CentralityMap>::type a_centrality_map;

    typedef typename property_traits<a_centrality_map>::value_type 

      centrality_type;

    typename graph_traits<Graph>::vertices_size_type V = num_vertices(g);

    std::vector<std::vector<edge_descriptor> > incoming(V);

    std::vector<centrality_type> distance(V);

    std::vector<centrality_type> dependency(V);

    std::vector<degree_size_type> path_count(V);

    brandes_betweenness_centrality(

      g, centrality, edge_centrality_map,

      make_iterator_property_map(incoming.begin(), vertex_index),

      make_iterator_property_map(distance.begin(), vertex_index),

      make_iterator_property_map(dependency.begin(), vertex_index),

      make_iterator_property_map(path_count.begin(), vertex_index),

      vertex_index);

  }

  template<typename WeightMap>

  struct brandes_betweenness_centrality_dispatch1

  {

    template<typename Graph, typename CentralityMap, 

             typename EdgeCentralityMap, typename VertexIndexMap>

    static void 

    run(const Graph& g, CentralityMap centrality, 

        EdgeCentralityMap edge_centrality_map, VertexIndexMap vertex_index,

        WeightMap weight_map)

    {

      brandes_betweenness_centrality_dispatch2(g, centrality, edge_centrality_map,

                                               weight_map, vertex_index);

    }

  };

  template<>

  struct brandes_betweenness_centrality_dispatch1<error_property_not_found>

  {

    template<typename Graph, typename CentralityMap, 

             typename EdgeCentralityMap, typename VertexIndexMap>

    static void 

    run(const Graph& g, CentralityMap centrality, 

        EdgeCentralityMap edge_centrality_map, VertexIndexMap vertex_index,

        error_property_not_found)

    {

      brandes_betweenness_centrality_dispatch2(g, centrality, edge_centrality_map,

                                               vertex_index);

    }

  };

} } // end namespace detail::graph

template<typename Graph, typename Param, typename Tag, typename Rest>

void 

brandes_betweenness_centrality(const Graph& g, 

                               const bgl_named_params<Param,Tag,Rest>& params)

{

  typedef bgl_named_params<Param,Tag,Rest> named_params;

  typedef typename property_value<named_params, edge_weight_t>::type ew;

  detail::graph::brandes_betweenness_centrality_dispatch1<ew>::run(

    g, 

    choose_param(get_param(params, vertex_centrality), 

                 dummy_property_map()),

    choose_param(get_param(params, edge_centrality), 

                 dummy_property_map()),

    choose_const_pmap(get_param(params, vertex_index), g, vertex_index),

    get_param(params, edge_weight));

}

template<typename Graph, typename CentralityMap>

void 

brandes_betweenness_centrality(const Graph& g, CentralityMap centrality)

{

  detail::graph::brandes_betweenness_centrality_dispatch2(

    g, centrality, dummy_property_map(), get(vertex_index, g));

}

template<typename Graph, typename CentralityMap, typename EdgeCentralityMap>

void 

brandes_betweenness_centrality(const Graph& g, CentralityMap centrality,

                               EdgeCentralityMap edge_centrality_map)

{

  detail::graph::brandes_betweenness_centrality_dispatch2(

    g, centrality, edge_centrality_map, get(vertex_index, g));

}

/**

 * Converts "absolute" betweenness centrality (as computed by the

 * brandes_betweenness_centrality algorithm) in the centrality map

 * into "relative" centrality. The result is placed back into the

 * given centrality map.

 */

template<typename Graph, typename CentralityMap>

void 

relative_betweenness_centrality(const Graph& g, CentralityMap centrality)

{

  typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;

  typedef typename property_traits<CentralityMap>::value_type centrality_type;

  typename graph_traits<Graph>::vertices_size_type n = num_vertices(g);

  centrality_type factor = centrality_type(2)/centrality_type(n*n - 3*n + 2);

  vertex_iterator v, v_end;

  for (tie(v, v_end) = vertices(g); v != v_end; ++v) {

    put(centrality, *v, factor * get(centrality, *v));

  }

}

// Compute the central point dominance of a graph.

template<typename Graph, typename CentralityMap>

typename property_traits<CentralityMap>::value_type

central_point_dominance(const Graph& g, CentralityMap centrality)

{

  using std::max;

  typedef typename graph_traits<Graph>::vertex_iterator vertex_iterator;

  typedef typename property_traits<CentralityMap>::value_type centrality_type;

  typename graph_traits<Graph>::vertices_size_type n = num_vertices(g);

  // Find max centrality

  centrality_type max_centrality(0);

  vertex_iterator v, v_end;

  for (tie(v, v_end) = vertices(g); v != v_end; ++v) {

    max_centrality = (max)(max_centrality, get(centrality, *v));

  }

  // Compute central point dominance

  centrality_type sum(0);

  for (tie(v, v_end) = vertices(g); v != v_end; ++v) {

    sum += (max_centrality - get(centrality, *v));

  }

  return sum/(n-1);

}

} // end namespace boost

#endif // BOOST_GRAPH_BRANDES_BETWEENNESS_CENTRALITY_HPP