Implement random_spanning_tree() and label_out_component()

src/graph/topology/Makefile.am

@@ -26,6 +26,7 @@ libgraph_tool_topology_la_SOURCES = \
     graph_minimum_spanning_tree.cc \
     graph_planar.cc \
     graph_random_matching.cc \
+    graph_random_spanning_tree.cc \
     graph_reciprocity.cc \
     graph_similarity.cc \
     graph_subgraph_isomorphism.cc \

src/graph/topology/graph_components.cc

@@ -50,9 +50,16 @@ do_label_biconnected_components(GraphInterface& gi, boost::any comp,
     return wrap_vector_owned(hist);
 }
 
+void do_label_out_component(GraphInterface& gi, size_t root, boost::any prop)
+{
+    run_action<>()(gi, bind<void>(label_out_component(), _1, _2, root),
+                   writable_vertex_scalar_properties())(prop);
+}
+
 void export_components()
 {
     python::def("label_components", &do_label_components);
     python::def("label_biconnected_components",
                 &do_label_biconnected_components);
+    python::def("label_out_component", &do_label_out_component);
 };

src/graph/topology/graph_components.hh

@@ -63,7 +63,7 @@ public:
         boost::put(_base_map, k, v);
 
         vector<size_t>& h = _hist;
-        size_t bin = v;        
+        size_t bin = v;
         if (bin > _max)
             return;
         if (bin >= h.size())
@@ -160,6 +160,33 @@ struct label_biconnected_components
 };
 
 
+struct label_out_component
+{
+    template <class CompMap>
+    class marker_visitor: public bfs_visitor<>
+    {
+    public:
+        marker_visitor() { }
+        marker_visitor(CompMap comp) : _comp(comp) { }
+
+        template <class Vertex, class Graph>
+        void discover_vertex(Vertex u, const Graph& g)
+        {
+            _comp[u] = true;
+        }
+    private:
+        CompMap _comp;
+    };
+
+    template <class Graph, class CompMap>
+    void operator()(const Graph& g, CompMap comp_map, size_t root) const
+    {
+        marker_visitor<CompMap> marker(comp_map);
+        breadth_first_search(g, vertex(root, g), visitor(marker));
+    }
+};
+
+
 } // graph_tool namespace
 
 #endif // GRAPH_COMPONENTS_HH

src/graph/topology/graph_random_spanning_tree.cc

+// graph-tool -- a general graph modification and manipulation thingy
+//
+// Copyright (C) 2007-2012 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_filtering.hh"
+#include "graph.hh"
+#include "graph_properties.hh"
+
+#include <boost/graph/random_spanning_tree.hpp>
+
+#if (GCC_VERSION >= 40400)
+#   include <tr1/random>
+#else
+#   include <boost/tr1/random.hpp>
+#endif
+
+using namespace std;
+using namespace boost;
+using namespace graph_tool;
+
+typedef tr1::mt19937 rng_t;
+
+struct get_random_span_tree
+{
+    template <class Graph, class IndexMap, class WeightMap, class TreeMap,
+              class RNG>
+    void operator()(const Graph& g, size_t root, IndexMap vertex_index,
+                    WeightMap weights, TreeMap tree_map, RNG& rng) const
+    {
+        typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
+        typedef typename graph_traits<Graph>::edge_descriptor edge_t;
+
+        unchecked_vector_property_map<vertex_t,IndexMap>
+            pred_map(vertex_index, num_vertices(g));
+        random_spanning_tree(g, rng, predecessor_map(pred_map).
+                             root_vertex(vertex(root, g)).
+                             weight_map(weights).
+                             vertex_index_map(vertex_index));
+
+        // convert the predecessor map to a tree map, and avoid trouble with
+        // parallel edges
+        int i, N = num_vertices(g);
+        #pragma omp parallel for default(shared) private(i) schedule(dynamic)
+        for (i = 0; i < N; ++i)
+        {
+            typename graph_traits<Graph>::vertex_descriptor v = vertex(i, g);
+            if (v == graph_traits<Graph>::null_vertex())
+                continue;
+
+            vector<edge_t> edges;
+            vector<typename property_traits<WeightMap>::value_type> ws;
+            typename graph_traits<Graph>::out_edge_iterator e, e_end;
+            for (tie(e,e_end) = out_edges(v, g); e != e_end; ++e)
+            {
+                if (target(*e,g) == pred_map[v])
+                {
+                    edges.push_back(*e);
+                    ws.push_back(weights[*e]);
+                }
+            }
+            if (!edges.empty())
+            {
+                edge_t e = *(edges.begin() +
+                             size_t(min_element(ws.begin(),
+                                                ws.end())-ws.begin()));
+                tree_map[e] = 1;
+            }
+        }
+    }
+};
+
+typedef property_map_types::apply<mpl::vector<uint8_t>,
+                                  GraphInterface::edge_index_map_t,
+                                  mpl::bool_<false> >::type
+    tree_properties;
+
+void get_random_spanning_tree(GraphInterface& gi, size_t root,
+                              boost::any weight_map, boost::any tree_map,
+                              size_t seed)
+{
+    rng_t rng(static_cast<rng_t::result_type>(seed));
+
+    typedef ConstantPropertyMap<size_t,GraphInterface::edge_t> cweight_t;
+
+    if (weight_map.empty())
+        weight_map = cweight_t(1);
+
+    typedef mpl::push_back<writable_edge_scalar_properties, cweight_t>::type
+        weight_maps;
+
+    run_action<>()
+        (gi, bind<void>(get_random_span_tree(), _1, root, gi.GetVertexIndex(),
+                        _2, _3, ref(rng)),
+         weight_maps(), tree_properties())(weight_map, tree_map);
+}

src/graph/topology/graph_topology.cc

@@ -39,6 +39,9 @@ void subgraph_isomorphism(GraphInterface& gi1, GraphInterface& gi2,
                           size_t n_max, size_t seed);
 double reciprocity(GraphInterface& gi);
 bool is_bipartite(GraphInterface& gi, boost::any part_map);
+void get_random_spanning_tree(GraphInterface& gi, size_t root,
+                              boost::any weight_map, boost::any tree_map,
+                              size_t seed);
 
 void export_components();
 void export_similarity();
@@ -61,6 +64,7 @@ BOOST_PYTHON_MODULE(libgraph_tool_topology)
     def("is_planar", &is_planar);
     def("reciprocity", &reciprocity);
     def("is_bipartite", &is_bipartite);
+    def("random_spanning_tree", &get_random_spanning_tree);
     export_components();
     export_similarity();
     export_dists();

src/graph_tool/topology/__init__.py

@@ -38,12 +38,14 @@ Summary
    max_cardinality_matching
    max_independent_vertex_set
    min_spanning_tree
+   random_spanning_tree
    dominator_tree
    topological_sort
    transitive_closure
    label_components
    label_biconnected_components
    label_largest_component
+   label_out_component
    is_bipartite
    is_planar
    edge_reciprocity
@@ -64,10 +66,12 @@ from .. flow import libgraph_tool_flow
 import random, sys, numpy
 __all__ = ["isomorphism", "subgraph_isomorphism", "mark_subgraph",
            "max_cardinality_matching", "max_independent_vertex_set",
-           "min_spanning_tree", "dominator_tree", "topological_sort",
-           "transitive_closure", "label_components", "label_largest_component",
-           "label_biconnected_components", "shortest_distance", "shortest_path",
-           "pseudo_diameter", "is_bipartite", "is_planar", "similarity", "edge_reciprocity"]
+           "min_spanning_tree", "random_spanning_tree", "dominator_tree",
+           "topological_sort", "transitive_closure", "label_components",
+           "label_largest_component", "label_biconnected_components",
+           "label_out_component", "shortest_distance", "shortest_path",
+           "pseudo_diameter", "is_bipartite", "is_planar", "similarity",
+           "edge_reciprocity"]
 
 
 def similarity(g1, g2, label1=None, label2=None, norm=True):
@@ -284,13 +288,13 @@ def min_spanning_tree(g, weights=None, root=None, tree_map=None):
     ----------
     g : :class:`~graph_tool.Graph`
         Graph to be used.
-    weights : :class:`~graph_tool.PropertyMap` (optional, default: None)
+    weights : :class:`~graph_tool.PropertyMap` (optional, default: `None`)
         The edge weights. If provided, the minimum spanning tree will minimize
         the edge weights.
-    root : :class:`~graph_tool.Vertex` (optional, default: None)
+    root : :class:`~graph_tool.Vertex` (optional, default: `None`)
         Root of the minimum spanning tree. If this is provided, Prim's algorithm
         is used. Otherwise, Kruskal's algorithm is used.
-    tree_map : :class:`~graph_tool.PropertyMap` (optional, default: None)
+    tree_map : :class:`~graph_tool.PropertyMap` (optional, default: `None`)
         If provided, the edge tree map will be written in this property map.
 
     Returns
@@ -361,6 +365,86 @@ def min_spanning_tree(g, weights=None, root=None, tree_map=None):
     return tree_map
 
 
+def random_spanning_tree(g, weights=None, root=None, tree_map=None):
+    """
+    Return a random spanning tree of a given graph, which can be directed or
+    undirected.
+
+    Parameters
+    ----------
+    g : :class:`~graph_tool.Graph`
+        Graph to be used.
+    weights : :class:`~graph_tool.PropertyMap` (optional, default: `None`)
+        The edge weights. If provided, the probability of a particular spanning
+        tree being selected is the product of its edge weights.
+    root : :class:`~graph_tool.Vertex` (optional, default: `None`)
+        Root of the spanning tree. If not provided, it will be selected randomly.
+    tree_map : :class:`~graph_tool.PropertyMap` (optional, default: `None`)
+        If provided, the edge tree map will be written in this property map.
+
+    Returns
+    -------
+    tree_map : :class:`~graph_tool.PropertyMap`
+        Edge property map with mark the tree edges: 1 for tree edge, 0
+        otherwise.
+
+    Notes
+    -----
+    The typical running time for random graphs is :math:`O(N\log N)`.
+
+    Examples
+    --------
+    >>> from numpy.random import seed, random
+    >>> seed(42)
+    >>> g, pos = gt.triangulation(random((400, 2)) * 10, type="delaunay")
+    >>> weight = g.new_edge_property("double")
+    >>> for e in g.edges():
+    ...    weight[e] = linalg.norm(pos[e.target()].a - pos[e.source()].a)
+    >>> tree = gt.random_spanning_tree(g, weights=weight)
+    >>> gt.graph_draw(g, pos=pos, output="rtriang_orig.pdf")
+    <...>
+    >>> g.set_edge_filter(tree)
+    >>> gt.graph_draw(g, pos=pos, output="triang_min_span_tree.pdf")
+    <...>
+
+
+    .. image:: rtriang_orig.*
+        :width: 400px
+    .. image:: triang_random_span_tree.*
+        :width: 400px
+
+    *Left:* Original graph, *Right:* A random spanning tree.
+
+    References
+    ----------
+
+    .. [wilson-generating-1996] David Bruce Wilson, "Generating random spanning
+       trees more quickly than the cover time", Proceedings of the twenty-eighth
+       annual ACM symposium on Theory of computing, Pages 296-303, ACM New York,
+       1996, :doi:`10.1145/237814.237880`
+    .. [boost-rst] http://www.boost.org/libs/graph/doc/random_spanning_tree.html
+    """
+    if tree_map is None:
+        tree_map = g.new_edge_property("bool")
+    if tree_map.value_type() != "bool":
+        raise ValueError("edge property 'tree_map' must be of value type bool.")
+
+    if root is None:
+        root = g.vertex(numpy.random.randint(0, g.num_vertices()),
+                        use_index=False)
+
+    # we need to restrict ourselves to the in-component of root
+    l = label_out_component(GraphView(g, reversed=True), root)
+    g = GraphView(g, vfilt=l)
+
+    seed = numpy.random.randint(0, sys.maxsize)
+    libgraph_tool_topology.\
+        random_spanning_tree(g._Graph__graph, int(root),
+                             _prop("e", g, weights),
+                             _prop("e", g, tree_map), seed)
+    return tree_map
+
+
 def dominator_tree(g, root, dom_map=None):
     """Return a vertex property map the dominator vertices for each vertex.
 
@@ -600,6 +684,54 @@ def label_largest_component(g, directed=None):
     return label
 
 
+def label_out_component(g, root):
+    """
+    Label the out-component (or simply the component for undirected graphs) of a
+    root vertex.
+
+    Parameters
+    ----------
+    g : :class:`~graph_tool.Graph`
+        Graph to be used.
+    root : :class:`~graph_tool.Vertex`
+        The root vertex.
+
+    Returns
+    -------
+    comp : :class:`~graph_tool.PropertyMap`
+         Boolean vertex property map which labels the out-component.
+
+    Notes
+    -----
+    The algorithm runs in :math:`O(V + E)` time.
+
+    Examples
+    --------
+    >>> from numpy.random import seed, poisson
+    >>> seed(43)
+    >>> g = gt.random_graph(100, lambda: poisson(1), directed=False)
+    >>> l = gt.label_out_component(g, g.vertex(0))
+    >>> print(l.a)
+    [1 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 0 0 1 0 0 1 1 1 1 0 0 0 1 0 0 0 0 0 0 0 1
+     1 1 0 0 0 0 1 0 1 1 0 0 0 1 1 0 0 1 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0
+     0 0 0 1 1 0 1 1 0 0 0 0 0 1 1 0 1 0 1 0 1 0 0 0 0 0]
+
+    The in-component can be obtained by reversing the graph.
+
+    >>> l = gt.label_out_component(GraphView(g, reversed=True), g.vertex(0))
+    >>> print(l.a)
+    [1 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 0 0 1 0 0 1 1 1 1 0 0 0 1 0 0 0 0 0 0 0 1
+     1 1 0 0 0 0 1 0 1 1 0 0 0 1 1 0 0 1 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0
+     0 0 0 1 1 0 1 1 0 0 0 0 0 1 1 0 1 0 1 0 1 0 0 0 0 0]
+    """
+
+    label = g.new_vertex_property("bool")
+    libgraph_tool_topology.\
+             label_out_component(g._Graph__graph, int(root),
+                                 _prop("v", g, label))
+    return label
+
+
 def label_biconnected_components(g, eprop=None, vprop=None):
     """
     Label the edges of biconnected components, and the vertices which are