Implement vertex/edge_percolation()

src/graph/topology/Makefile.am

@@ -27,6 +27,7 @@ libgraph_tool_topology_la_SOURCES = \
     graph_maximal_planar.cc \
     graph_maximal_vertex_set.cc \
     graph_minimum_spanning_tree.cc \
+    graph_percolation.cc \
     graph_planar.cc \
     graph_random_matching.cc \
     graph_random_spanning_tree.cc \
@@ -45,5 +46,6 @@ libgraph_tool_topology_la_SOURCES = \
 libgraph_tool_topology_la_include_HEADERS = \
     graph_components.hh \
     graph_kcore.hh \
+    graph_percolation.hh \
     graph_similarity.hh \
     graph_vertex_similarity.hh

src/graph/topology/graph_percolation.cc

+// graph-tool -- a general graph modification and manipulation thingy
+//
+// Copyright (C) 2006-2017 Tiago de Paula Peixoto <tiago@skewed.de>
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License
+// as published by the Free Software Foundation; either version 3
+// of the License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+#include "graph_tool.hh"
+#include "numpy_bind.hh"
+
+#include "graph_percolation.hh"
+
+using namespace std;
+using namespace boost;
+using namespace graph_tool;
+
+void percolate_edge(GraphInterface& gi, boost::any tree, boost::any size,
+                    python::object edges, python::object max_size)
+{
+    typedef property_map_type::apply<int64_t,
+                                     GraphInterface::vertex_index_map_t>::type
+        tree_t;
+    tree_t tree_map;
+    try
+    {
+        tree_map = any_cast<tree_t>(tree);
+    }
+    catch (bad_any_cast&)
+    {
+        throw GraphException("tree map must be a vertex property map of value type int64_t");
+    }
+
+    tree_t size_map;
+    try
+    {
+        size_map = any_cast<tree_t>(size);
+    }
+    catch (bad_any_cast&)
+    {
+        throw GraphException("size map must be a vertex property map of value type int64_t");
+    }
+
+    multi_array_ref<uint64_t, 2> es = get_array<uint64_t, 2>(edges);
+    multi_array_ref<uint64_t, 1> ms = get_array<uint64_t, 1>(max_size);
+
+    run_action<graph_tool::detail::never_directed>()
+        (gi, [&](auto& g){ edge_percolate(g, tree_map, size_map, ms, es); })();
+}
+
+
+void percolate_vertex(GraphInterface& gi, boost::any tree, boost::any size,
+                      python::object vertices, python::object max_size)
+{
+    typedef property_map_type::apply<int64_t,
+                                     GraphInterface::vertex_index_map_t>::type
+        tree_t;
+    tree_t tree_map;
+    try
+    {
+        tree_map = any_cast<tree_t>(tree);
+    }
+    catch (bad_any_cast&)
+    {
+        throw GraphException("tree map must be a vertex property map of value type int64_t");
+    }
+
+    tree_t size_map;
+    try
+    {
+        size_map = any_cast<tree_t>(size);
+    }
+    catch (bad_any_cast&)
+    {
+        throw GraphException("size map must be a vertex property map of value type int64_t");
+    }
+
+    multi_array_ref<uint64_t, 1> vs = get_array<uint64_t, 1>(vertices);
+    multi_array_ref<uint64_t, 1> ms = get_array<uint64_t, 1>(max_size);
+
+    run_action<graph_tool::detail::never_directed>()
+        (gi, [&](auto& g){ vertex_percolate(g, tree_map, size_map, ms, vs); })();
+}
+
+#include <boost/python.hpp>
+
+void export_percolation()
+{
+    using namespace boost::python;
+
+    def("percolate_edge", percolate_edge);
+    def("percolate_vertex", percolate_vertex);
+};

src/graph/topology/graph_percolation.hh

+// graph-tool -- a general graph modification and manipulation thingy
+//
+// Copyright (C) 2006-2017 Tiago de Paula Peixoto <tiago@skewed.de>
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License
+// as published by the Free Software Foundation; either version 3
+// of the License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef GRAPH_PERCOLATION_HH
+#define GRAPH_PERCOLATION_HH
+
+namespace graph_tool
+{
+using namespace std;
+using namespace boost;
+
+template <class Graph, class TreeMap>
+typename graph_traits<Graph>::vertex_descriptor
+find_root(size_t vi, TreeMap tree, Graph& g,
+          vector<size_t>& temp)
+{
+    auto parent = vertex(vi, g);
+    temp.clear();
+    while(size_t(tree[vertex(parent, g)]) != parent)
+    {
+        temp.push_back(parent);
+        parent = tree[vertex(parent, g)];
+    }
+    // path compression
+    for (auto v: temp)
+        tree[vertex(v, g)] = parent;
+    return vertex(parent, g);
+}
+
+
+template <class Graph, class TreeMap, class SizeMap>
+void join_cluster(const pair<size_t, size_t>& e, TreeMap tree, SizeMap size,
+                  size_t& max_size, Graph& g, vector<size_t>& temp)
+{
+    auto rs = find_root(e.first, tree, g, temp);
+    auto rt = find_root(e.second, tree, g, temp);
+    if (rt != rs)
+    {
+        if (size[rs] < size[rt])
+            swap(rs, rt);
+        tree[rt] = rs;
+        size[rs] += size[rt];
+        max_size = max(size_t(size[rs]), max_size);
+    }
+}
+
+template <class Graph, class TreeMap,
+          class SizeMap, class MaxSize, class Edges>
+void edge_percolate(Graph& g, TreeMap tree, SizeMap size, MaxSize& max_size,
+                    Edges& edges)
+{
+    vector<size_t> temp;
+    size_t s = 0;
+    for (size_t i = 0; i < edges.size(); ++i)
+    {
+        join_cluster(make_pair(size_t(edges[i][0]), size_t(edges[i][1])),
+                     tree, size, s, g, temp);
+        max_size[i] = s;
+    }
+
+    boost::multi_array_ref<typename Edges::element, 1>
+        vertices(edges.data(),
+                 boost::extents[edges.num_elements()]);
+
+    for (auto& v : vertices)
+    {
+        auto root = find_root(v, tree, g, temp);
+        size[v] = size[root];
+    }
+};
+
+template <class Graph,  class TreeMap, class SizeMap, class MaxSize,
+          class Vertices>
+void vertex_percolate(Graph& g, TreeMap tree, SizeMap size, MaxSize& max_size,
+                      Vertices& vertices)
+{
+    vector<size_t> temp;
+    size_t s = 0;
+    for (size_t i = 0; i < vertices.size(); ++i)
+    {
+        auto v = vertex(vertices[i], g);
+        if (v == graph_traits<Graph>::null_vertex())
+        {
+            max_size[i] = s;
+            continue;
+        }
+
+        for (auto a : adjacent_vertices_range(v, g))
+        {
+            join_cluster(make_pair(size_t(v), size_t(a)), tree, size, s, g,
+                         temp);
+        }
+        max_size[i] = s;
+    }
+
+    //flatten tree
+    for (auto vi : vertices)
+    {
+        auto v = vertex(vi, g);
+        if (v == graph_traits<Graph>::null_vertex())
+            continue;
+        auto root = find_root(vi, tree, g, temp);
+        size[v] = size[root];
+    }
+}
+
+} // graph_tool namespace
+
+#endif // GRAPH_PERCOLATION_HH

src/graph/topology/graph_topology.cc

@@ -55,6 +55,7 @@ vector<int32_t> get_tsp(GraphInterface& gi, size_t src, boost::any weight_map);
 
 void export_components();
 void export_kcore();
+void export_percolation();
 void export_similarity();
 void export_dists();
 void export_all_dists();
@@ -83,6 +84,7 @@ BOOST_PYTHON_MODULE(libgraph_tool_topology)
     def("get_tsp", &get_tsp);
     export_components();
     export_kcore();
+    export_percolation();
     export_similarity();
     export_dists();
     export_all_dists();

src/graph_tool/topology/__init__.py

@@ -53,6 +53,8 @@ Summary
    label_biconnected_components
    label_largest_component
    label_out_component
+   vertex_percolation
+   edge_percolation
    kcore_decomposition
    is_bipartite
    is_DAG
@@ -82,9 +84,10 @@ __all__ = ["isomorphism", "subgraph_isomorphism", "mark_subgraph",
            "topological_sort", "transitive_closure", "tsp_tour",
            "sequential_vertex_coloring", "label_components",
            "label_largest_component", "label_biconnected_components",
-           "label_out_component", "kcore_decomposition", "shortest_distance",
-           "shortest_path", "all_shortest_paths", "all_predecessors",
-           "all_paths", "all_circuits", "pseudo_diameter", "is_bipartite", "is_DAG",
+           "label_out_component", "vertex_percolation", "edge_percolation",
+           "kcore_decomposition", "shortest_distance", "shortest_path",
+           "all_shortest_paths", "all_predecessors", "all_paths",
+           "all_circuits", "pseudo_diameter", "is_bipartite", "is_DAG",
            "is_planar", "make_maximal_planar", "similarity", "vertex_similarity",
            "edge_reciprocity"]
 
@@ -1255,6 +1258,169 @@ def label_biconnected_components(g, eprop=None, vprop=None):
                                           _prop("v", g, vprop))
     return eprop, vprop, hist
 
+def vertex_percolation(g, vertices):
+    """Compute the size of the largest component as vertices are (virtually)
+    removed from the graph.
+
+    Parameters
+    ----------
+    g : :class:`~graph_tool.Graph`
+        Graph to be used.
+    vertices : :class:`numpy.ndarray` or iterable of ints
+        List of vertices in reversed order of removal.
+
+    Returns
+    -------
+    size : :class:`numpy.ndarray`
+        Size of the largest component prior to removal of each vertex.
+    comp : :class:`~graph_tool.PropertyMap`
+        Vertex property map with component labels.
+
+    Notes
+    -----
+
+    The algorithm runs in :math:`O(V + E)` time.
+
+    Examples
+    --------
+    .. testcode::
+       :hide:
+
+       import numpy.random
+       numpy.random.seed(42)
+       gt.seed_rng(42)
+
+    >>> g = gt.random_graph(10000, lambda: geometric(1./4) + 1, directed=False)
+    >>> vertices = sorted([v for v in g.vertices()], key=lambda v: v.out_degree())
+    >>> sizes, comp = gt.vertex_percolation(g, vertices)
+    >>> numpy.random.shuffle(vertices)
+    >>> sizes2, comp = gt.vertex_percolation(g, vertices)
+    >>> figure()
+    <...>
+    >>> plot(sizes, label="Targeted")
+    [...]
+    >>> plot(sizes2, label="Random")
+    [...]
+    >>> xlabel("Vertices remaining")
+    <...>
+    >>> ylabel("Size of largest component")
+    <...>
+    >>> legend(loc="lower right")
+    <...>
+    >>> savefig("vertex-percolation.svg")
+
+    .. figure:: vertex-percolation.*
+        :align: center
+
+        Targeted and random vertex percolation of a random graph with an
+        exponential degree distribution.
+
+    References
+    ----------
+    .. [newman-ziff] M. E. J. Newman, R. M. Ziff, "A fast Monte Carlo algorithm
+       for site or bond percolation", Phys. Rev. E 64, 016706 (2001)
+       :doi:`10.1103/PhysRevE.64.016706`, :arxiv:`cond-mat/0101295`
+
+    """
+    vertices = numpy.asarray(vertices, dtype="uint64")
+
+    tree = g.vertex_index.copy("int64_t")
+    size = g.new_vertex_property("int64_t", 1)
+    max_size = numpy.zeros(len(vertices), dtype="uint64")
+
+    u = GraphView(g, directed=False)
+
+    libgraph_tool_topology.\
+        percolate_vertex(u._Graph__graph,
+                         _prop("v", u, tree),
+                         _prop("v", u, size),
+                         vertices, max_size)
+
+    return max_size, tree
+
+def edge_percolation(g, edges):
+    """Compute the size of the largest component as edges are (virtually)
+    removed from the graph.
+
+    Parameters
+    ----------
+    g : :class:`~graph_tool.Graph`
+        Graph to be used.
+    edges : :class:`numpy.ndarray` or iterable of pairs of ints
+        List of edges in reversed order of removal. If the type is
+        :class:`numpy.ndarray`, it should have a shape ``(E, 2)``, where ``E``
+        is the number of edges, such that ``edges[i,0]`` and ``edges[i,1]`` are
+        the both endpoints of edge ``i``.
+
+    Returns
+    -------
+    size : :class:`numpy.ndarray`
+        Size of the largest component prior to removal of each edge.
+    comp : :class:`~graph_tool.PropertyMap`
+        Vertex property map with component labels.
+
+    Notes
+    -----
+
+    The algorithm runs in :math:`O(E)` time.
+
+    Examples
+    --------
+    .. testcode::
+       :hide:
+
+       import numpy.random
+       numpy.random.seed(42)
+       gt.seed_rng(42)
+
+    >>> g = gt.random_graph(10000, lambda: geometric(1./4) + 1, directed=False)
+    >>> edges = sorted([(e.source(), e.target()) for e in g.edges()],
+    ...                key=lambda e: e[0].out_degree() * e[1].out_degree())
+    >>> sizes, comp = gt.edge_percolation(g, edges)
+    >>> numpy.random.shuffle(edges)
+    >>> sizes2, comp = gt.edge_percolation(g, edges)
+    >>> figure()
+    <...>
+    >>> plot(sizes, label="Targeted")
+    [...]
+    >>> plot(sizes2, label="Random")
+    [...]
+    >>> xlabel("Edges remaining")
+    <...>
+    >>> ylabel("Size of largest component")
+    <...>
+    >>> legend(loc="lower right")
+    <...>
+    >>> savefig("edge-percolation.svg")
+
+    .. figure:: edge-percolation.*
+        :align: center
+
+        Targeted and random edge percolation of a random graph with an
+        exponential degree distribution.
+
+    References
+    ----------
+    .. [newman-ziff] M. E. J. Newman, R. M. Ziff, "A fast Monte Carlo algorithm
+       for site or bond percolation", Phys. Rev. E 64, 016706 (2001)
+       :doi:`10.1103/PhysRevE.64.016706`, :arxiv:`cond-mat/0101295`
+
+    """
+    edges = numpy.asarray(edges, dtype="uint64")
+
+    tree = g.vertex_index.copy("int64_t")
+    size = g.new_vertex_property("int64_t", 1)
+    max_size = numpy.zeros(len(edges), dtype="uint64")
+
+    u = GraphView(g, directed=False)
+
+    libgraph_tool_topology.\
+        percolate_edge(u._Graph__graph,
+                       _prop("v", u, tree),
+                       _prop("v", u, size),
+                       edges, max_size)
+    return max_size, tree
+
 def kcore_decomposition(g, deg="out", vprop=None):
     """
     Perform a k-core decomposition of the given graph.