graph/test/betweenness_centrality_test.cpp

// Copyright 2004 The Trustees of Indiana University.

// Use, modification and distribution is subject to 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
#include <boost/graph/betweenness_centrality.hpp>

#include <boost/graph/adjacency_list.hpp>
#include <vector>
#include <stack>
#include <queue>
#include <boost/property_map/property_map.hpp>
#include <boost/core/lightweight_test.hpp>
#include <boost/random/uniform_01.hpp>
#include <boost/random/linear_congruential.hpp>
#include <boost/lexical_cast.hpp>

using namespace boost;

const double error_tolerance = 0.001;

typedef property< edge_weight_t, double, property< edge_index_t, std::size_t > >
    EdgeProperties;

struct weighted_edge
{
    int source, target;
    double weight;
};

template < typename Graph >
void run_weighted_test(Graph*, int V, weighted_edge edge_init[], int E,
    double correct_centrality[])
{
    Graph g(V);
    typedef typename graph_traits< Graph >::vertex_descriptor Vertex;
    typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator;
    typedef typename graph_traits< Graph >::edge_descriptor Edge;

    std::vector< Vertex > vertices(V);
    {
        vertex_iterator v, v_end;
        int index = 0;
        for (boost::tie(v, v_end) = boost::vertices(g); v != v_end;
             ++v, ++index)
        {
            put(vertex_index, g, *v, index);
            vertices[index] = *v;
        }
    }

    std::vector< Edge > edges(E);
    for (int e = 0; e < E; ++e)
    {
        edges[e] = add_edge(
            vertices[edge_init[e].source], vertices[edge_init[e].target], g)
                       .first;
        put(edge_weight, g, edges[e], 1.0);
    }

    std::vector< double > centrality(V);
    brandes_betweenness_centrality(g,
        centrality_map(make_iterator_property_map(
                           centrality.begin(), get(vertex_index, g), double()))
            .vertex_index_map(get(vertex_index, g))
            .weight_map(get(edge_weight, g)));

    for (int v = 0; v < V; ++v)
    {
        BOOST_TEST(centrality[v] == correct_centrality[v]);
    }
}

struct unweighted_edge
{
    int source, target;
};

template < typename Graph >
void run_unweighted_test(Graph*, int V, unweighted_edge edge_init[], int E,
    double correct_centrality[], double* correct_edge_centrality = 0)
{
    Graph g(V);
    typedef typename graph_traits< Graph >::vertex_descriptor Vertex;
    typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator;
    typedef typename graph_traits< Graph >::edge_descriptor Edge;

    std::vector< Vertex > vertices(V);
    {
        vertex_iterator v, v_end;
        int index = 0;
        for (boost::tie(v, v_end) = boost::vertices(g); v != v_end;
             ++v, ++index)
        {
            put(vertex_index, g, *v, index);
            vertices[index] = *v;
        }
    }

    std::vector< Edge > edges(E);
    for (int e = 0; e < E; ++e)
    {
        edges[e] = add_edge(
            vertices[edge_init[e].source], vertices[edge_init[e].target], g)
                       .first;
        put(edge_weight, g, edges[e], 1.0);
        put(edge_index, g, edges[e], e);
    }

    std::vector< double > centrality(V);
    std::vector< double > edge_centrality1(E);

    brandes_betweenness_centrality(g,
        centrality_map(make_iterator_property_map(
                           centrality.begin(), get(vertex_index, g), double()))
            .edge_centrality_map(make_iterator_property_map(
                edge_centrality1.begin(), get(edge_index, g), double()))
            .vertex_index_map(get(vertex_index, g)));

    std::vector< double > centrality2(V);
    std::vector< double > edge_centrality2(E);
    brandes_betweenness_centrality(g,
        vertex_index_map(get(vertex_index, g))
            .weight_map(get(edge_weight, g))
            .centrality_map(make_iterator_property_map(
                centrality2.begin(), get(vertex_index, g), double()))
            .edge_centrality_map(make_iterator_property_map(
                edge_centrality2.begin(), get(edge_index, g), double())));

    std::vector< double > edge_centrality3(E);
    brandes_betweenness_centrality(g,
        edge_centrality_map(make_iterator_property_map(
            edge_centrality3.begin(), get(edge_index, g), double())));

    for (int v = 0; v < V; ++v)
    {
        BOOST_TEST(centrality[v] == centrality2[v]);

        double relative_error = correct_centrality[v] == 0.0
            ? centrality[v]
            : (centrality[v] - correct_centrality[v]) / correct_centrality[v];
        if (relative_error < 0)
            relative_error = -relative_error;
        BOOST_TEST(relative_error < error_tolerance);
    }

    for (int e = 0; e < E; ++e)
    {
        BOOST_TEST(edge_centrality1[e] == edge_centrality2[e]);
        BOOST_TEST(edge_centrality1[e] == edge_centrality3[e]);

        if (correct_edge_centrality)
        {
            double relative_error = correct_edge_centrality[e] == 0.0
                ? edge_centrality1[e]
                : (edge_centrality1[e] - correct_edge_centrality[e])
                    / correct_edge_centrality[e];
            if (relative_error < 0)
                relative_error = -relative_error;
            BOOST_TEST(relative_error < error_tolerance);

            if (relative_error >= error_tolerance)
            {
                std::cerr << "Edge " << e << " has edge centrality "
                          << edge_centrality1[e] << ", should be "
                          << correct_edge_centrality[e] << std::endl;
            }
        }
    }
}

template < typename Graph > void run_wheel_test(Graph*, int V)
{
    typedef typename graph_traits< Graph >::vertex_descriptor Vertex;
    typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator;
    typedef typename graph_traits< Graph >::edge_descriptor Edge;

    Graph g(V);
    Vertex center = *boost::vertices(g).first;

    std::vector< Vertex > vertices(V);
    {
        vertex_iterator v, v_end;
        int index = 0;
        for (boost::tie(v, v_end) = boost::vertices(g); v != v_end;
             ++v, ++index)
        {
            put(vertex_index, g, *v, index);
            vertices[index] = *v;
            if (*v != center)
            {
                Edge e = add_edge(*v, center, g).first;
                put(edge_weight, g, e, 1.0);
            }
        }
    }

    std::vector< double > centrality(V);
    brandes_betweenness_centrality(g,
        make_iterator_property_map(
            centrality.begin(), get(vertex_index, g), double()));

    std::vector< double > centrality2(V);
    brandes_betweenness_centrality(g,
        centrality_map(make_iterator_property_map(
                           centrality2.begin(), get(vertex_index, g), double()))
            .vertex_index_map(get(vertex_index, g))
            .weight_map(get(edge_weight, g)));

    relative_betweenness_centrality(g,
        make_iterator_property_map(
            centrality.begin(), get(vertex_index, g), double()));

    relative_betweenness_centrality(g,
        make_iterator_property_map(
            centrality2.begin(), get(vertex_index, g), double()));

    for (int v = 0; v < V; ++v)
    {
        BOOST_TEST(centrality[v] == centrality2[v]);
        BOOST_TEST(
            (v == 0 && centrality[v] == 1) || (v != 0 && centrality[v] == 0));
    }

    double dominance = central_point_dominance(g,
        make_iterator_property_map(
            centrality2.begin(), get(vertex_index, g), double()));
    BOOST_TEST(dominance == 1.0);
}

template < typename MutableGraph >
void randomly_add_edges(MutableGraph& g, double edge_probability)
{
    typedef typename graph_traits< MutableGraph >::directed_category
        directed_category;

    minstd_rand gen;
    uniform_01< minstd_rand, double > rand_gen(gen);

    typedef typename graph_traits< MutableGraph >::vertex_descriptor vertex;
    typename graph_traits< MutableGraph >::vertex_iterator vi, vi_end;
    for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
    {
        vertex v = *vi;
        typename graph_traits< MutableGraph >::vertex_iterator wi
            = is_same< directed_category, undirected_tag >::value
            ? vi
            : vertices(g).first;
        while (wi != vi_end)
        {
            vertex w = *wi++;
            if (v != w)
            {
                if (rand_gen() < edge_probability)
                    add_edge(v, w, g);
            }
        }
    }
}

template < typename Graph, typename VertexIndexMap, typename CentralityMap >
void simple_unweighted_betweenness_centrality(
    const Graph& g, VertexIndexMap index, CentralityMap centrality)
{
    typedef typename boost::graph_traits< Graph >::vertex_descriptor vertex;
    typedef
        typename boost::graph_traits< Graph >::vertex_iterator vertex_iterator;
    typedef typename boost::graph_traits< Graph >::adjacency_iterator
        adjacency_iterator;
    typedef typename boost::graph_traits< Graph >::vertices_size_type
        vertices_size_type;
    typedef typename boost::property_traits< CentralityMap >::value_type
        centrality_type;

    vertex_iterator vi, vi_end;
    for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi)
        put(centrality, *vi, 0);

    vertex_iterator si, si_end;
    for (boost::tie(si, si_end) = vertices(g); si != si_end; ++si)
    {
        vertex s = *si;

        // S <-- empty stack
        std::stack< vertex > S;

        // P[w] <-- empty list, w \in V
        typedef std::vector< vertex > Predecessors;
        std::vector< Predecessors > predecessors(num_vertices(g));

        // sigma[t] <-- 0, t \in V
        std::vector< vertices_size_type > sigma(num_vertices(g), 0);

        // sigma[s] <-- 1
        sigma[get(index, s)] = 1;

        // d[t] <-- -1, t \in V
        std::vector< int > d(num_vertices(g), -1);

        // d[s] <-- 0
        d[get(index, s)] = 0;

        // Q <-- empty queue
        std::queue< vertex > Q;

        // enqueue s --> Q
        Q.push(s);

        while (!Q.empty())
        {
            // dequeue v <-- Q
            vertex v = Q.front();
            Q.pop();

            // push v --> S
            S.push(v);

            adjacency_iterator wi, wi_end;
            for (boost::tie(wi, wi_end) = adjacent_vertices(v, g); wi != wi_end;
                 ++wi)
            {
                vertex w = *wi;

                // w found for the first time?
                if (d[get(index, w)] < 0)
                {
                    // enqueue w --> Q
                    Q.push(w);

                    // d[w] <-- d[v] + 1
                    d[get(index, w)] = d[get(index, v)] + 1;
                }

                // shortest path to w via v?
                if (d[get(index, w)] == d[get(index, v)] + 1)
                {
                    // sigma[w] = sigma[w] + sigma[v]
                    sigma[get(index, w)] += sigma[get(index, v)];

                    // append v --> P[w]
                    predecessors[get(index, w)].push_back(v);
                }
            }
        }

        // delta[v] <-- 0, v \in V
        std::vector< centrality_type > delta(num_vertices(g), 0);

        // S returns vertices in order of non-increasing distance from s
        while (!S.empty())
        {
            // pop w <-- S
            vertex w = S.top();
            S.pop();

            const Predecessors& w_preds = predecessors[get(index, w)];
            for (typename Predecessors::const_iterator vi = w_preds.begin();
                 vi != w_preds.end(); ++vi)
            {
                vertex v = *vi;
                // delta[v] <-- delta[v] + (sigma[v]/sigma[w])*(1 + delta[w])
                delta[get(index, v)] += ((centrality_type)sigma[get(index, v)]
                                            / sigma[get(index, w)])
                    * (1 + delta[get(index, w)]);
            }

            if (w != s)
            {
                // C_B[w] <-- C_B[w] + delta[w]
                centrality[w] += delta[get(index, w)];
            }
        }
    }

    typedef typename graph_traits< Graph >::directed_category directed_category;
    const bool is_undirected
        = is_same< directed_category, undirected_tag >::value;
    if (is_undirected)
    {
        vertex_iterator v, v_end;
        for (boost::tie(v, v_end) = vertices(g); v != v_end; ++v)
        {
            put(centrality, *v, get(centrality, *v) / centrality_type(2));
        }
    }
}

template < typename Graph > void random_unweighted_test(Graph*, int n)
{
    Graph g(n);

    {
        typename graph_traits< Graph >::vertex_iterator v, v_end;
        int index = 0;
        for (boost::tie(v, v_end) = boost::vertices(g); v != v_end;
             ++v, ++index)
        {
            put(vertex_index, g, *v, index);
        }
    }

    randomly_add_edges(g, 0.20);

    std::cout << "Random graph with " << n << " vertices and " << num_edges(g)
              << " edges.\n";

    std::cout << "  Direct translation of Brandes' algorithm...";
    std::vector< double > centrality(n);
    simple_unweighted_betweenness_centrality(g, get(vertex_index, g),
        make_iterator_property_map(
            centrality.begin(), get(vertex_index, g), double()));
    std::cout << "DONE.\n";

    std::cout << "  Real version, unweighted...";
    std::vector< double > centrality2(n);
    brandes_betweenness_centrality(g,
        make_iterator_property_map(
            centrality2.begin(), get(vertex_index, g), double()));
    std::cout << "DONE.\n";

    if (!std::equal(centrality.begin(), centrality.end(), centrality2.begin()))
    {
        for (std::size_t v = 0; v < centrality.size(); ++v)
        {
            double relative_error = centrality[v] == 0.0
                ? centrality2[v]
                : (centrality2[v] - centrality[v]) / centrality[v];
            if (relative_error < 0)
                relative_error = -relative_error;
            BOOST_TEST(relative_error < error_tolerance);
        }
    }

    std::cout << "  Real version, weighted...";
    std::vector< double > centrality3(n);

    for (typename graph_traits< Graph >::edge_iterator ei = edges(g).first;
         ei != edges(g).second; ++ei)
        put(edge_weight, g, *ei, 1);

    brandes_betweenness_centrality(g,
        weight_map(get(edge_weight, g))
            .centrality_map(make_iterator_property_map(
                centrality3.begin(), get(vertex_index, g), double())));
    std::cout << "DONE.\n";

    if (!std::equal(centrality.begin(), centrality.end(), centrality3.begin()))
    {
        for (std::size_t v = 0; v < centrality.size(); ++v)
        {
            double relative_error = centrality[v] == 0.0
                ? centrality3[v]
                : (centrality3[v] - centrality[v]) / centrality[v];
            if (relative_error < 0)
                relative_error = -relative_error;
            BOOST_TEST(relative_error < error_tolerance);
        }
    }
}

int main(int argc, char* argv[])
{
    int random_test_num_vertices = 300;
    if (argc >= 2)
        random_test_num_vertices = boost::lexical_cast< int >(argv[1]);
    typedef adjacency_list< listS, listS, undirectedS,
        property< vertex_index_t, int >, EdgeProperties >
        Graph;
    typedef adjacency_list< listS, listS, directedS,
        property< vertex_index_t, int >, EdgeProperties >
        Digraph;

    struct unweighted_edge ud_edge_init1[5]
        = { { 0, 1 }, { 0, 3 }, { 1, 2 }, { 3, 2 }, { 2, 4 } };
    double ud_centrality1[5] = { 0.5, 1.0, 3.5, 1.0, 0.0 };
    run_unweighted_test((Graph*)0, 5, ud_edge_init1, 5, ud_centrality1);

    // Example borrowed from the JUNG test suite
    struct unweighted_edge ud_edge_init2[10] = {
        { 0, 1 },
        { 0, 6 },
        { 1, 2 },
        { 1, 3 },
        { 2, 4 },
        { 3, 4 },
        { 4, 5 },
        { 5, 8 },
        { 7, 8 },
        { 6, 7 },
    };
    double ud_centrality2[9]
        = { 0.2142 * 28, 0.2797 * 28, 0.0892 * 28, 0.0892 * 28, 0.2797 * 28,
              0.2142 * 28, 0.1666 * 28, 0.1428 * 28, 0.1666 * 28 };
    double ud_edge_centrality2[10] = { 10.66666, 9.33333, 6.5, 6.5, 6.5, 6.5,
        10.66666, 9.33333, 8.0, 8.0 };

    run_unweighted_test(
        (Graph*)0, 9, ud_edge_init2, 10, ud_centrality2, ud_edge_centrality2);

    weighted_edge dw_edge_init1[6] = { { 0, 1, 1.0 }, { 0, 3, 1.0 },
        { 1, 2, 0.5 }, { 3, 1, 1.0 }, { 3, 4, 1.0 }, { 4, 2, 0.5 } };
    double dw_centrality1[5] = { 0.0, 1.5, 0.0, 1.0, 0.5 };
    run_weighted_test((Digraph*)0, 5, dw_edge_init1, 6, dw_centrality1);

    run_wheel_test((Graph*)0, 15);

    random_unweighted_test((Graph*)0, random_test_num_vertices);

    return boost::report_errors();
}