Modify shortest_distances() and include shortest_path()

This implements a shortest_path() function which returns the shortest
path between two vertices.

This also modifies shortest_distances, such that a predecessor map is
optionally returned.

src/graph/topology/graph_distance.cc

@@ -37,10 +37,9 @@ 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)
+    void tree_edge(typename graph_traits<Graph>::edge_descriptor e,
+                   Graph& g)
     {
         _pred[target(e,g)] = source(e,g);
     }
@@ -49,7 +48,7 @@ public:
     void discover_vertex(typename graph_traits<Graph>::vertex_descriptor v,
                          Graph& g)
     {
-        if (_pred[v] == v)
+        if (size_t(_pred[v]) == v)
             return;
         size_t dist = _dist_map[_pred[v]] + 1;
         if (dist > _max_dist)
@@ -93,9 +92,10 @@ private:
 
 struct do_bfs_search
 {
-    template <class Graph, class VertexIndexMap, class DistMap>
+    template <class Graph, class VertexIndexMap, class DistMap, class PredMap>
     void operator()(const Graph& g, size_t source, VertexIndexMap vertex_index,
-                    DistMap dist_map, long double max_dist) const
+                    DistMap dist_map, PredMap pred_map, long double max_dist)
+        const
     {
         typedef typename property_traits<DistMap>::value_type dist_t;
         dist_t max_d = (max_dist > 0) ?
@@ -112,17 +112,13 @@ struct do_bfs_search
         }
         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>
+                                 visitor(bfs_max_visitor<DistMap, PredMap>
                                          (dist_map, pred_map, max_d)).
                                  vertex_index_map(vertex_index).
                                  color_map(color_map));
@@ -133,10 +129,11 @@ struct do_bfs_search
 
 struct do_djk_search
 {
-    template <class Graph, class VertexIndexMap, class DistMap, class WeightMap>
+    template <class Graph, class VertexIndexMap, class DistMap, class PredMap,
+              class WeightMap>
     void operator()(const Graph& g, size_t source, VertexIndexMap vertex_index,
-                    DistMap dist_map, WeightMap weight, long double max_dist)
-        const
+                    DistMap dist_map, PredMap pred_map, WeightMap weight,
+                    long double max_dist) const
     {
         typedef typename property_traits<DistMap>::value_type dist_t;
         dist_t max_d = (max_dist > 0) ?
@@ -148,6 +145,7 @@ struct do_djk_search
                                     weight_map(weight).
                                     distance_map(dist_map).
                                     vertex_index_map(vertex_index).
+                                    predecessor_map(pred_map).
                                     visitor(djk_max_visitor<DistMap>
                                             (dist_map, max_d)));
         }
@@ -156,21 +154,29 @@ struct do_djk_search
 };
 
 void get_dists(GraphInterface& gi, size_t source, boost::any dist_map,
-               boost::any weight, long double max_dist)
+               boost::any weight, boost::any pred_map, long double max_dist)
 {
+    typedef property_map_type
+        ::apply<int64_t, GraphInterface::vertex_index_map_t>::type pred_map_t;
+
+    pred_map_t pmap = any_cast<pred_map_t>(pred_map);
+
     if (weight.empty())
     {
         run_action<>()
             (gi, bind<void>(do_bfs_search(), _1, source, gi.GetVertexIndex(),
-                            _2, max_dist),
-             writable_vertex_scalar_properties())
+                            _2, pmap.get_unchecked(num_vertices(gi.GetGraph())),
+                            max_dist),
+             writable_vertex_scalar_properties(),
+             mpl::vector<pred_map_t>())
             (dist_map);
     }
     else
     {
         run_action<>()
             (gi, bind<void>(do_djk_search(), _1, source, gi.GetVertexIndex(),
-                            _2, _3, max_dist),
+                            _2, pmap.get_unchecked(num_vertices(gi.GetGraph())),
+                            _3, max_dist),
              writable_vertex_scalar_properties(),
              edge_scalar_properties())
             (dist_map, weight);

src/graph_tool/topology/__init__.py

@@ -51,7 +51,8 @@ import random, sys, numpy, weakref
 __all__ = ["isomorphism", "subgraph_isomorphism", "mark_subgraph",
            "min_spanning_tree", "dominator_tree", "topological_sort",
            "transitive_closure", "label_components",
-           "label_biconnected_components", "shortest_distance", "is_planar"]
+           "label_biconnected_components", "shortest_distance",
+           "shortest_path", "is_planar"]
 
 def isomorphism(g1, g2, isomap=False):
     r"""Check whether two graphs are isomorphic.
@@ -165,7 +166,7 @@ def mark_subgraph(g, sub, vmap, emap, vmask=None, emask=None):
 
     This returns a vertex and an edge property map, with value type 'bool',
     indicating whether or not a vertex/edge in `g` corresponds to the subgraph
-    `sub`. 
+    `sub`.
     """
     if vmask == None:
         vmask = g.new_vertex_property("bool")
@@ -537,7 +538,8 @@ def label_biconnected_components(g, eprop=None, vprop=None):
     return eprop, vprop, nc
 
 def shortest_distance(g, source=None, weights=None, max_dist=None,
-                      directed=None, dense=False, dist_map=None):
+                      directed=None, dense=False, dist_map=None,
+                      pred_map=False):
     """
     Calculate the distance of all vertices from a given source, or the all pairs
     shortest paths, if the source is not specified.
@@ -547,7 +549,7 @@ def shortest_distance(g, source=None, weights=None, max_dist=None,
     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
+        Source vertex 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
@@ -565,6 +567,9 @@ def shortest_distance(g, source=None, weights=None, max_dist=None,
     dist_map : :class:`~graph_tool.PropertyMap` (optional, default: None)
         Vertex property to store the distances. If none is supplied, one
         is created.
+    pred_map : bool (optional, default: False)
+        If true, a vertex property map with the predecessors is returned.
+        Ignored if source=None.
 
     Returns
     -------
@@ -668,9 +673,11 @@ def shortest_distance(g, source=None, weights=None, max_dist=None,
 
     try:
         if source != None:
+            pmap = g.copy_property(g.vertex_index, value_type="int64_t")
             libgraph_tool_topology.get_dists(g._Graph__graph, int(source),
                                              _prop("v", g, dist_map),
                                              _prop("e", g, weights),
+                                             _prop("v", g, pmap),
                                              float(max_dist))
         else:
             libgraph_tool_topology.get_all_dists(g._Graph__graph,
@@ -680,7 +687,104 @@ def shortest_distance(g, source=None, weights=None, max_dist=None,
     finally:
         if directed != None:
             g.pop_filter(directed=True)
-    return dist_map
+    if source != None and pred_map:
+        return dist_map, pmap
+    else:
+        return dist_map
+
+def shortest_path(g, source, target, weights=None, pred_map=None):
+    """
+    Return the shortest path from `source` to `target`.
+
+    Parameters
+    ----------
+    g : :class:`~graph_tool.Graph`
+        Graph to be used.
+    source : :class:`~graph_tool.Vertex`
+        Source vertex of the search.
+    source : :class:`~graph_tool.Vertex`
+        Target vertex of the search.
+    weights : :class:`~graph_tool.PropertyMap` (optional, default: None)
+        The edge weights. If provided, the minimum spanning tree will minimize
+        the edge weights.
+    pred_map :  :class:`~graph_tool.PropertyMap` (optional, default: None)
+        Vertex property map with the predecessors in the search tree. If this is
+        provided, the shortest paths are not computed, and are obtained directly
+        from this map.
+
+    Returns
+    -------
+    vertex_list : list of :class:`~graph_tool.Vertex`
+        List of vertices from `source` to `target` in the shortest path.
+    edge_list : list of :class:`~graph_tool.Edge`
+        List of edges from `source` to `target` in the shortest path.
+
+    Notes
+    -----
+
+    The paths are computed with a breadth-first search (BFS) or Dijkstra's
+    algorithm [dijkstra]_, if weights are given.
+
+    The algorithm runs in :math:`O(V + E)` time, or :math:`O(V \log V)` if
+    weights are given.
+
+    Examples
+    --------
+    >>> from numpy.random import seed, poisson
+    >>> seed(42)
+    >>> g = gt.random_graph(300, lambda: (poisson(3), poisson(3)))
+    >>> vlist, elist = gt.shortest_path(g, g.vertex(10), g.vertex(11))
+    >>> print [str(v) for v in vlist]
+    ['10', '226', '288', '94', '221', '169', '203', '38', '281', '11']
+    >>> print [str(e) for e in elist]
+    ['(10,226)', '(226,288)', '(288,94)', '(94,221)', '(221,169)', '(169,203)', '(203,38)', '(38,281)', '(281,11)']
+
+    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
+    """
+
+    if pred_map == None:
+        dists, pred_map = shortest_distance(g, source, weights=weights,
+                                            pred_map=True)
+
+    if pred_map[target] == int(target): # no path to source
+        return [], []
+
+    vlist = [target]
+    elist = []
+
+    if weights != None:
+        max_w = weights.a.max() + 1
+    else:
+        max_w = None
+
+    v = target
+    while v != source:
+        p = g.vertex(pred_map[v])
+        min_w = max_w
+        pe = None
+        s = None
+        for e in v.in_edges() if g.is_directed() else v.out_edges():
+            s = e.source() if g.is_directed() else e.target()
+            if s == p:
+                if weights != None:
+                    if weights[e] < min_w:
+                        min_w = weights[e]
+                        pe = e
+                else:
+                    pe = e
+                    break
+        elist.insert(0, pe)
+        vlist.insert(0, p)
+        v = p
+    return vlist, elist
+
 
 def is_planar(g, embedding=False, kuratowski=False):
     """