Include minimum distance algorithms

This includes the shortest_distance() function in the topology module.

src/graph/topology/Makefile.am

@@ -21,7 +21,9 @@ libgraph_tool_topology_la_SOURCES = \
     graph_dominator_tree.cc \
     graph_topological_sort.cc \
     graph_transitive_closure.cc \
-    graph_components.cc
+    graph_components.cc \
+    graph_distance.cc \
+    graph_all_distances.cc
 
 libgraph_tool_topology_la_include_HEADERS = \
     graph_components.hh

src/graph/topology/graph_all_distances.cc

+// Copyright (C) 2008  Tiago de Paula Peixoto <tiago@forked.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.hh"
+#include "graph_filtering.hh"
+#include "graph_properties.hh"
+#include "graph_selectors.hh"
+
+#include <boost/python.hpp>
+
+#include <boost/graph/johnson_all_pairs_shortest.hpp>
+#include <boost/graph/floyd_warshall_shortest.hpp>
+
+using namespace std;
+using namespace boost;
+using namespace graph_tool;
+
+struct do_all_pairs_search
+{
+    template <class Graph, class VertexIndexMap, class DistMap, class WeightMap>
+    void operator()(const Graph& g, VertexIndexMap vertex_index,
+                    DistMap dist_map, WeightMap weight, bool dense) const
+    {
+        typedef typename property_traits<DistMap>::value_type::value_type
+            dist_t;
+
+        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;
+            dist_map[v].clear();
+            dist_map[v].resize(num_vertices(g), 0);
+        }
+
+        if (dense)
+        {
+            floyd_warshall_all_pairs_shortest_paths
+                (g, dist_map,
+                 weight_map(ConvertedPropertyMap<WeightMap,dist_t>(weight)).
+                 vertex_index_map(vertex_index));
+        }
+        else
+        {
+            johnson_all_pairs_shortest_paths
+                (g, dist_map,
+                 weight_map(ConvertedPropertyMap<WeightMap,dist_t>(weight)).
+                 vertex_index_map(vertex_index));
+        }
+    }
+};
+
+void get_all_dists(GraphInterface& gi, boost::any dist_map, boost::any weight,
+                   bool dense)
+{
+    typedef ConstantPropertyMap<size_t,GraphInterface::edge_t> cweight_map_t;
+
+    if (weight.empty())
+        weight = boost::any(cweight_map_t(1));
+
+    run_action<>()
+        (gi,
+         bind<void>(do_all_pairs_search(), _1, gi.GetVertexIndex(),
+                    _2, _3, dense),
+         vertex_scalar_vector_properties(),
+         mpl::push_back<edge_scalar_properties,cweight_map_t>::type())
+        (dist_map, weight);
+}
+
+void export_all_dists()
+{
+    python::def("get_all_dists", &get_all_dists);
+};
+

src/graph/topology/graph_distance.cc

+// Copyright (C) 2008  Tiago de Paula Peixoto <tiago@forked.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.hh"
+#include "graph_filtering.hh"
+#include "graph_properties.hh"
+#include "graph_selectors.hh"
+
+#include <boost/graph/breadth_first_search.hpp>
+#include <boost/graph/dijkstra_shortest_paths.hpp>
+
+#include <boost/python.hpp>
+
+using namespace std;
+using namespace boost;
+using namespace graph_tool;
+
+struct stop_search {};
+
+template <class DistMap, class PredMap>
+class bfs_max_visitor:
+    public boost::bfs_visitor<null_visitor>
+{
+public:
+    bfs_max_visitor(DistMap dist_map, PredMap pred, size_t max_dist)
+        : _dist_map(dist_map), _pred(pred), _max_dist(max_dist), _dist(0) {}
+
+
+    template <class Graph>
+    void examine_edge(typename graph_traits<Graph>::edge_descriptor e,
+                      Graph& g)
+    {
+        _pred[target(e,g)] = source(e,g);
+    }
+
+    template <class Graph>
+    void discover_vertex(typename graph_traits<Graph>::vertex_descriptor v,
+                         Graph& g)
+    {
+        if (_pred[v] == v)
+            return;
+        size_t dist = _dist_map[_pred[v]] + 1;
+        if (dist > _max_dist)
+            throw stop_search();
+        _dist_map[v] = dist;
+    }
+
+private:
+    DistMap _dist_map;
+    PredMap _pred;
+    size_t _max_dist;
+    size_t _dist;
+};
+
+template <class DistMap, class PredMap>
+class djk_max_visitor:
+    public boost::dijkstra_visitor<null_visitor>
+{
+public:
+    djk_max_visitor(DistMap dist_map, PredMap pred_map,
+                    typename property_traits<DistMap>::value_type max_dist)
+        : _dist_map(dist_map), _max_dist(max_dist) {}
+
+    template <class Graph>
+    void examine_vertex(typename graph_traits<Graph>::vertex_descriptor u,
+                        Graph& g)
+    {
+        if (_dist_map[_pred_map[u]] >= _max_dist)
+            throw stop_search();
+    }
+
+private:
+    DistMap _dist_map;
+    PredMap _pred_map;
+    typename property_traits<DistMap>::value_type _max_dist;
+};
+
+
+struct do_bfs_search
+{
+    template <class Graph, class VertexIndexMap, class DistMap>
+    void operator()(const Graph& g, size_t source, VertexIndexMap vertex_index,
+                    DistMap dist_map, long double max_dist) const
+    {
+        typedef typename property_traits<DistMap>::value_type dist_t;
+        dist_t max_d = (max_dist > 0) ?
+            max_dist : numeric_limits<dist_t>::max();
+
+        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;
+            dist_map[v] = numeric_limits<dist_t>::max();
+        }
+        dist_map[vertex(source,g)] = 0;
+
+        typedef unchecked_vector_property_map
+            <typename graph_traits<Graph>::vertex_descriptor, VertexIndexMap>
+            pred_map_t;
+        pred_map_t pred_map(vertex_index, num_vertices(g));
+        pred_map[vertex(source, g)] = vertex(source, g);
+        unchecked_vector_property_map<boost::default_color_type, VertexIndexMap>
+            color_map(vertex_index, num_vertices(g));
+        try
+        {
+            breadth_first_search(g, vertex(source, g),
+                                 visitor(bfs_max_visitor<DistMap, pred_map_t>
+                                         (dist_map, pred_map, max_d)).
+                                 vertex_index_map(vertex_index).
+                                 color_map(color_map));
+        }
+        catch (stop_search&) {}
+    }
+};
+
+struct do_djk_search
+{
+    template <class Graph, class VertexIndexMap, class DistMap, class WeightMap>
+    void operator()(const Graph& g, size_t source, VertexIndexMap vertex_index,
+                    DistMap dist_map, WeightMap weight, long double max_dist)
+        const
+    {
+        typedef typename property_traits<DistMap>::value_type dist_t;
+        dist_t max_d = (max_dist > 0) ?
+            max_dist : numeric_limits<dist_t>::max();
+
+        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;
+            dist_map[v] = numeric_limits<dist_t>::max();
+        }
+        dist_map[vertex(source,g)] = 0;
+        typedef unchecked_vector_property_map
+            <typename graph_traits<Graph>::vertex_descriptor, VertexIndexMap>
+            pred_map_t;
+        pred_map_t pred_map(vertex_index, num_vertices(g));
+        unchecked_vector_property_map<boost::default_color_type, VertexIndexMap>
+            color_map(vertex_index, num_vertices(g));
+        try
+        {
+            dijkstra_shortest_paths(g, vertex(source, g),
+                                    visitor(djk_max_visitor<DistMap, pred_map_t>
+                                            (dist_map, pred_map, max_d)).
+                                    weight_map(weight).
+                                    distance_map(dist_map).
+                                    predecessor_map(pred_map).
+                                    vertex_index_map(vertex_index).
+                                    color_map(color_map));
+        }
+        catch (stop_search&) {}
+    }
+};
+
+void get_dists(GraphInterface& gi, size_t source, boost::any dist_map,
+               boost::any weight, long double max_dist)
+{
+    if (weight.empty())
+    {
+        run_action<>()
+            (gi, bind<void>(do_bfs_search(), _1, source, gi.GetVertexIndex(),
+                            _2, max_dist),
+             writable_vertex_scalar_properties())
+            (dist_map);
+    }
+    else
+    {
+        run_action<>()
+            (gi, bind<void>(do_djk_search(), _1, source, gi.GetVertexIndex(),
+                            _2, _3, max_dist),
+             writable_vertex_scalar_properties(),
+             edge_scalar_properties())
+            (dist_map, weight);
+
+    }
+}
+
+void export_dists()
+{
+    python::def("get_dists", &get_dists);
+};
+

src/graph/topology/graph_topology.cc

@@ -24,20 +24,16 @@ using namespace graph_tool;
 
 bool check_isomorphism(GraphInterface& gi1, GraphInterface& gi2,
                        boost::any iso_map);
-
 void get_kruskal_spanning_tree(GraphInterface& gi, boost::any weight_map,
                                boost::any tree_map);
-
 void get_prim_spanning_tree(GraphInterface& gi, size_t root,
                             boost::any weight_map, boost::any tree_map);
-
 void topological_sort(GraphInterface& gi, vector<int32_t>& sort);
-
 void dominator_tree(GraphInterface& gi, size_t entry, boost::any pred_map);
-
 void transitive_closure(GraphInterface& gi, GraphInterface& tcgi);
-
 void export_components();
+void export_dists();
+void export_all_dists();
 
 BOOST_PYTHON_MODULE(libgraph_tool_topology)
 {
@@ -48,4 +44,6 @@ BOOST_PYTHON_MODULE(libgraph_tool_topology)
     def("dominator_tree", &dominator_tree);
     def("transitive_closure", &transitive_closure);
     export_components();
+    export_dists();
+    export_all_dists();
 }

src/graph_tool/topology/__init__.py

@@ -26,6 +26,7 @@ Summary
 .. autosummary::
    :nosignatures:
 
+   shortest_distance
    isomorphism
    min_spanning_tree
    dominator_tree
@@ -42,11 +43,11 @@ from .. dl_import import dl_import
 dl_import("import libgraph_tool_topology")
 
 from .. core import _prop, Vector_int32_t, _check_prop_writable, \
-     _check_prop_scalar, Graph
+     _check_prop_scalar,  _check_prop_vector, Graph
 import random, sys, numpy
 __all__ = ["isomorphism", "min_spanning_tree", "dominator_tree",
            "topological_sort", "transitive_closure", "label_components",
-           "label_biconnected_components"]
+           "label_biconnected_components", "shortest_distance"]
 
 def isomorphism(g1, g2, isomap=False):
     """Check whether two graphs are isomorphisms. If `isomap` is True, a vertex
@@ -413,3 +414,149 @@ def label_biconnected_components(g, eprop=None, vprop=None):
     finally:
         g.pop_filter(directed=True)
     return eprop, vprop, nc
+
+def shortest_distance(g, source=None, weights=None, max_dist=None,
+                      directed=None, dense=False, dist_map=None):
+    """
+    Calculate the distance of all vertices from a given source, or the all pairs
+    shortest paths, if the source is not specified.
+
+    Parameters
+    ----------
+    g : :class:`~graph_tool.Graph`
+        Graph to be used.
+    source : :class:`~graph_tool.Vertex` (optional, default: None)
+        Vertex source of the search. If unspecified, the all pairs shortest
+        distances are computed.
+    weights : :class:`~graph_tool.PropertyMap` (optional, default: None)
+        The edge weights. If provided, the minimum spanning tree will minimize
+        the edge weights.
+    max_dist : scalar value (optional, default: None)
+        If specified, this limits the maximum distance of the vertices
+        are searched. This parameter has no effect if source == None.
+    directed : bool (optional, default:None)
+        Treat graph as directed or not, independently of its actual
+        directionality.
+    dense : bool (optional, default: False)
+        If true, and source == None, the Floyd-Warshall algorithm is used,
+        otherwise the Johnson algorithm is used. If source != None, this option
+        has no effect.
+    dist_map : :class:`~graph_tool.PropertyMap` (optional, default: None)
+        Vertex property to store the distances. If none is supplied, one
+        is created.
+
+    Returns
+    -------
+    dist_map : :class:`~graph_tool.PropertyMap`
+        Vertex property map with the distances from source. If source is 'None',
+        it will have a vector value type, with the distances to every vertex.
+
+    Notes
+    -----
+
+    If a source is given, the distances are calculated with a breadth-first
+    search (BFS) or Dijkstra's algorithm [dijkstra]_, if weights are given. If
+    source is not given, the distances are calculated with Johnson's algorithm
+    [johnson-apsp]_. If dense=True, the Floyd-Warshall algorithm
+    [floyd-warshall-apsp]_ is used instead.
+
+    If source is specified, the algorithm runs in :math:`O(V + E)` time, or
+    :math:`O(V \log V)` if weights are given. If source is not specified, it
+    runs in :math:`O(VE\log V)` time, or :math:`O(V^3)` if dense == True.
+
+    Examples
+    --------
+    >>> from numpy.random import seed, poisson
+    >>> seed(42)
+    >>> g = gt.random_graph(100, lambda: (poisson(3), poisson(3)))
+    >>> dist = gt.shortest_distance(g, source=g.vertex(0))
+    >>> print dist.get_array()
+    [         0          2 2147483647          4 2147483647          6
+              4          4          3          4          4          5
+              5          2          4          5          5          5
+              4          4          6          6 2147483647          5
+              4          4          4          6          4          4
+              5          5          3          3          4          4
+              2          3          3          4 2147483647 2147483647
+              4          4          3          3          1          5
+              5          4          5          2          4          4
+              4          1          3          2          3          4
+              3          5          5          1          3 2147483647
+              5          5          5          4          3 2147483647
+              3          2          3          3          3          5
+              4          4          4          3          5 2147483647
+              5          6          4          5          3          5
+              5          4          5          4          1          6
+              4          3          3          4]
+    >>> dist = gt.shortest_distance(g)
+    >>> print array(dist[g.vertex(0)])
+    [         0          2 2147483647          4 2147483647          6
+              4          4          3          4          4          5
+              5          2          4          5          5          5
+              4          4          6          6 2147483647          5
+              4          4          4          6          4          4
+              5          5          3          3          4          4
+              2          3          3          4 2147483647 2147483647
+              4          4          3          3          1          5
+              5          4          5          2          4          4
+              4          1          3          2          3          4
+              3          5          5          1          3 2147483647
+              5          5          5          4          3 2147483647
+              3          2          3          3          3          5
+              4          4          4          3          5 2147483647
+              5          6          4          5          3          5
+              5          4          5          4          1          6
+              4          3          3          4]
+
+
+    References
+    ----------
+    .. [bfs] Edward Moore, "The shortest path through a maze", International
+       Symposium on the Theory of Switching (1959), Harvard University
+       Press;http://www.boost.org/libs/graph/doc/breadth_first_search.html
+    .. [dijkstra] E. Dijkstra, "A note on two problems in connexion with
+       graphs." Numerische Mathematik, 1:269-271, 1959.
+       http://www.boost.org/libs/graph/doc/dijkstra_shortest_paths.html
+    .. [johnson-apsp] http://www.boost.org/libs/graph/doc/johnson_all_pairs_shortest.html
+    .. [floyd-warshall-apsp] http://www.boost.org/libs/graph/doc/floyd_warshall_shortest.html
+    """
+
+    if weights == None:
+        dist_type = 'int32_t'
+    else:
+        dist_type = weights.value_type()
+
+    if dist_map == None:
+        if source != None:
+            dist_map = g.new_vertex_property(dist_type)
+        else:
+            dist_map = g.new_vertex_property("vector<%s>" % dist_type)
+
+    _check_prop_writable(dist_map, name="dist_map")
+    if source != None:
+        _check_prop_scalar(dist_map, name="dist_map")
+    else:
+        _check_prop_vector(dist_map, name="dist_map")
+
+    if max_dist == None:
+        max_dist = 0
+
+    if directed != None:
+        g.stash_filter(directed=True)
+        g.set_directed(directed)
+
+    try:
+        if source != None:
+            libgraph_tool_topology.get_dists(g._Graph__graph, int(source),
+                                             _prop("v", g, dist_map),
+                                             _prop("e", g, weights),
+                                             float(max_dist))
+        else:
+            libgraph_tool_topology.get_all_dists(g._Graph__graph,
+                                                 _prop("v", g, dist_map),
+                                                 _prop("e", g, weights), dense)
+
+    finally:
+        if directed != None:
+            g.pop_filter(directed=True)
+    return dist_map