Merge pull request #13 from Fkawala/master

update: topology.shortest_distance and topology.shortest_path to sear…

src/graph/topology/graph_distance.cc

@@ -22,7 +22,7 @@
 
 #include <boost/graph/breadth_first_search.hpp>
 #include <boost/graph/dijkstra_shortest_paths.hpp>
-
+#include <boost/python/stl_iterator.hpp>
 #include <boost/python.hpp>
 
 using namespace std;
@@ -75,6 +75,57 @@ private:
     size_t _dist;
 };
 
+template <class DistMap, class PredMap>
+class bfs_max_multiple_targets_visitor:
+    public boost::bfs_visitor<null_visitor>
+{
+public:
+    bfs_max_multiple_targets_visitor(DistMap dist_map, PredMap pred, size_t max_dist, std::unordered_set<std::size_t> target)
+        : _dist_map(dist_map), _pred(pred), _max_dist(max_dist), _target(target),
+          _dist(0) {}
+
+    template <class Graph>
+    void tree_edge(typename graph_traits<Graph>::edge_descriptor e,
+                   Graph& g)
+    {
+        _pred[target(e,g)] = source(e,g);
+    }
+
+    template <class Graph>
+    void examine_vertex(typename graph_traits<Graph>::vertex_descriptor v,
+                        Graph&)
+    {
+        typedef typename property_traits<DistMap>::value_type val_t;
+        if ( _dist_map[v] > val_t(_max_dist))
+            throw stop_search();
+    }
+
+    template <class Graph>
+    void discover_vertex(typename graph_traits<Graph>::vertex_descriptor v,
+                         Graph&)
+    {
+        if (size_t(_pred[v]) == v)
+            return;
+        _dist_map[v] = _dist_map[_pred[v]] + 1;
+
+        auto search = _target.find(v);
+        if (search != _target.end())
+        {
+            _target.erase(*search);
+            if (_target.empty())
+                throw stop_search();
+        };
+    }
+
+private:
+    DistMap _dist_map;
+    PredMap _pred;
+    size_t _max_dist;
+    std::unordered_set<std::size_t> _target;
+    size_t _dist;
+};
+
+
 template <class DistMap>
 class djk_max_visitor:
     public boost::dijkstra_visitor<null_visitor>
@@ -104,14 +155,52 @@ private:
 };
 
 
+template <class DistMap>
+class djk_max_multiple_targets_visitor:
+    public boost::dijkstra_visitor<null_visitor>
+{
+public:
+    djk_max_multiple_targets_visitor(DistMap dist_map,
+                                     typename property_traits<DistMap>::value_type max_dist, 
+                                     std::unordered_set<std::size_t> target)
+        : _dist_map(dist_map), _max_dist(max_dist), _target(target) {}
+
+
+    template <class Graph>
+    void examine_vertex(typename graph_traits<Graph>::vertex_descriptor u,
+                        Graph&)
+    {
+        if (_dist_map[u] > _max_dist)
+            throw stop_search();
+        
+        auto search = _target.find(u);
+        if (search != _target.end())
+        {
+            _target.erase(*search);
+            if (_target.empty())
+                throw stop_search();
+        };
+    }
+
+
+private:
+    DistMap _dist_map;
+    typename property_traits<DistMap>::value_type _max_dist;
+    std::unordered_set<std::size_t> _target;
+};
+
+
 struct do_bfs_search
 {
     template <class Graph, class VertexIndexMap, class DistMap, class PredMap>
-    void operator()(const Graph& g, size_t source, size_t target,
+    void operator()(const Graph& g, size_t source, boost::python::list target_list,
                     VertexIndexMap vertex_index, DistMap dist_map,
                     PredMap pred_map, long double max_dist) const
     {
         typedef typename property_traits<DistMap>::value_type dist_t;
+
+        boost::python::stl_input_iterator<size_t> begin(target_list), end;
+        std::unordered_set<std::size_t> tgt(begin, end);        
         dist_t max_d = (max_dist > 0) ?
             max_dist : numeric_limits<dist_t>::max();
 
@@ -123,14 +212,28 @@ struct do_bfs_search
 
         pred_map[vertex(source, g)] = vertex(source, g);
         unchecked_vector_property_map<boost::default_color_type, VertexIndexMap>
-            color_map(vertex_index, num_vertices(g));
+        color_map(vertex_index, num_vertices(g));
+
         try
-        {
-            breadth_first_search(g, vertex(source, g),
-                                 visitor(bfs_max_visitor<DistMap, PredMap>
-                                         (dist_map, pred_map, max_d, target)).
-                                 vertex_index_map(vertex_index).
-                                 color_map(color_map));
+        {   
+            if (tgt.size() <= 1)
+            {
+                size_t target = tgt.empty() ? graph_traits<GraphInterface::multigraph_t>::null_vertex() : *tgt.begin();
+                breadth_first_search(g, vertex(source, g),
+                                     visitor(bfs_max_visitor<DistMap, PredMap>
+                                             (dist_map, pred_map, max_d, target)).
+                                     vertex_index_map(vertex_index).
+                                     color_map(color_map));
+            }
+            else
+            {
+                breadth_first_search(g, vertex(source, g),
+                                     visitor(bfs_max_multiple_targets_visitor<DistMap, PredMap>
+                                             (dist_map, pred_map, max_d, tgt)).
+                                     vertex_index_map(vertex_index).
+                                     color_map(color_map));
+            }
+
         }
         catch (stop_search&) {}
     }
@@ -140,13 +243,15 @@ struct do_djk_search
 {
     template <class Graph, class VertexIndexMap, class DistMap, class PredMap,
               class WeightMap>
-    void operator()(const Graph& g, size_t source, size_t target,
+    void operator()(const Graph& g, size_t source, boost::python::list target_list,
                     VertexIndexMap vertex_index, 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) ?
-            max_dist : numeric_limits<dist_t>::max();
+        max_dist : numeric_limits<dist_t>::max();
+        boost::python::stl_input_iterator<size_t> begin(target_list), end;
+        std::unordered_set<std::size_t> tgt(begin, end);
 
         int i, N = num_vertices(g);
         #pragma omp parallel for default(shared) private(i) schedule(runtime) if (N > 100)
@@ -155,33 +260,46 @@ struct do_djk_search
         dist_map[source] = 0;
 
         try
-        {
-            dijkstra_shortest_paths(g, vertex(source, g),
-                                    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, target)));
+        {   
+            if (tgt.size() <= 1)
+            {
+                size_t target = tgt.empty() ? graph_traits<GraphInterface::multigraph_t>::null_vertex() : *tgt.begin();
+                dijkstra_shortest_paths(g, vertex(source, g),
+                                        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, target)));
+            }
+            else
+            {
+                dijkstra_shortest_paths(g, vertex(source, g),
+                                        weight_map(weight).
+                                        distance_map(dist_map).
+                                        vertex_index_map(vertex_index).
+                                        predecessor_map(pred_map).
+                                        visitor(djk_max_multiple_targets_visitor<DistMap>
+                                                (dist_map, max_d, tgt)));
+            }
+
         }
         catch (stop_search&) {}
     }
 };
 
-void get_dists(GraphInterface& gi, size_t source, int tgt, boost::any dist_map,
-               boost::any weight, boost::any pred_map, long double max_dist)
+void get_dists(GraphInterface& gi, size_t source, boost::python::list tgt, 
+               boost::any dist_map, 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);
 
-    size_t target = tgt < 0 ? graph_traits<GraphInterface::multigraph_t>::null_vertex() : tgt;
-
     if (weight.empty())
     {
         run_action<>()
-            (gi, std::bind(do_bfs_search(), placeholders::_1, source, target, gi.GetVertexIndex(),
+            (gi, std::bind(do_bfs_search(), placeholders::_1, source, tgt, gi.GetVertexIndex(),
                            placeholders::_2, pmap.get_unchecked(num_vertices(gi.GetGraph())),
                            max_dist),
              writable_vertex_scalar_properties())
@@ -190,7 +308,7 @@ void get_dists(GraphInterface& gi, size_t source, int tgt, boost::any dist_map,
     else
     {
         run_action<>()
-            (gi, std::bind(do_djk_search(), placeholders::_1, source, target, gi.GetVertexIndex(),
+            (gi, std::bind(do_djk_search(), placeholders::_1, source, tgt, gi.GetVertexIndex(),
                            placeholders::_2, pmap.get_unchecked(num_vertices(gi.GetGraph())),
                            placeholders::_3, max_dist),
              writable_vertex_scalar_properties(),

src/graph_tool/topology/__init__.py

@@ -69,7 +69,7 @@ from .. import _prop, Vector_int32_t, _check_prop_writable, \
      _check_prop_scalar, _check_prop_vector, Graph, PropertyMap, GraphView,\
      libcore, _get_rng, _degree, perfect_prop_hash
 from .. stats import label_self_loops
-import random, sys, numpy
+import random, sys, numpy, collections
 
 __all__ = ["isomorphism", "subgraph_isomorphism", "mark_subgraph",
            "max_cardinality_matching", "max_independent_vertex_set",
@@ -1110,8 +1110,8 @@ def shortest_distance(g, source=None, target=None, weights=None, max_dist=None,
     source : :class:`~graph_tool.Vertex` (optional, default: None)
         Source vertex of the search. If unspecified, the all pairs shortest
         distances are computed.
-    target : :class:`~graph_tool.Vertex` (optional, default: None)
-        Target vertex of the search. If unspecified, the distance to all
+    target : :class:`~graph_tool.Vertex` or iterable of such objects (optional, default: None)
+        Target vertex/vertices of the search. If unspecified, the distance to all
         vertices from the source will be computed.
     weights : :class:`~graph_tool.PropertyMap` (optional, default: None)
         The edge weights. If provided, the minimum spanning tree will minimize
@@ -1216,6 +1216,12 @@ def shortest_distance(g, source=None, target=None, weights=None, max_dist=None,
     .. [floyd-warshall-apsp] http://www.boost.org/libs/graph/doc/floyd_warshall_shortest.html
     """
 
+    if not isinstance(target, list):
+        if isinstance(target, collections.Iterable):
+            target = list(target)
+        else:
+            target = [target]
+
     if weights is None:
         dist_type = 'int32_t'
     else:
@@ -1241,14 +1247,11 @@ def shortest_distance(g, source=None, target=None, weights=None, max_dist=None,
     else:
         u = g
 
-    if target is None:
-        target = -1
-
     if source is not None:
         pmap = g.copy_property(u.vertex_index, value_type="int64_t")
         libgraph_tool_topology.get_dists(g._Graph__graph,
                                          int(source),
-                                         int(target),
+                                         list(target),
                                          _prop("v", g, dist_map),
                                          _prop("e", g, weights),
                                          _prop("v", g, pmap),
@@ -1258,9 +1261,12 @@ def shortest_distance(g, source=None, target=None, weights=None, max_dist=None,
                                              _prop("v", g, dist_map),
                                              _prop("e", g, weights), dense)
 
-
-    if source is not None and target != -1:
-        dist_map = dist_map[target]
+    if source is not None:
+        if len(target) > 1:
+            dist_map = numpy.array([dist_map[target] for target in target])
+        else:
+            # Standard behaviour when there is one single target
+            dist_map = dist_map[target[0]]
 
     if source is not None and pred_map:
         return dist_map, pmap
@@ -1278,8 +1284,8 @@ def shortest_path(g, source, target, weights=None, pred_map=None):
         Graph to be used.
     source : :class:`~graph_tool.Vertex`
         Source vertex of the search.
-    target : :class:`~graph_tool.Vertex`
-        Target vertex of the search.
+    target : :class:`~graph_tool.Vertex` or iterable of such objects
+        Target vertex or vertices of the search.
     weights : :class:`~graph_tool.PropertyMap` (optional, default: None)
         The edge weights.
     pred_map :  :class:`~graph_tool.PropertyMap` (optional, default: None)
@@ -1330,43 +1336,67 @@ def shortest_path(g, source, target, weights=None, pred_map=None):
        graphs." Numerische Mathematik, 1:269-271, 1959.
     .. [dijkstra-boost] http://www.boost.org/libs/graph/doc/dijkstra_shortest_paths.html
     """
+    
+    if not isinstance(target, collections.Iterable):
+        target = [target]
+    
+    vlists, elists = {}, {}
+    targets = set(target)
 
     if pred_map is None:
-        pred_map = shortest_distance(g, source, target,
-                                     weights=weights,
-                                     pred_map=True)[1]
-
-    if pred_map[target] == int(target):  # no path to target
-        return [], []
+        _, pred_map = shortest_distance(g, source, target,
+                                        weights=weights,
+                                        pred_map=True)
 
-    vlist = [target]
-    elist = []
+    if any(pred_map[tgt] != int(tgt) for tgt in target):
+        # there is a path to at least one of the targets 
 
-    if weights is not 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 is not 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
+        if weights is not None:
+            max_w = weights.a.max() + 1
+        else:
+            max_w = None
+
+        for tgt in target:
+
+            if tgt in vlists:
+                # the current target has already been found
+                continue
+
+            v = tgt
+            vlist, elist = [tgt], []
+            targets_indexes = {}
+
+            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 is not None:
+                            if weights[e] < min_w:
+                                min_w = weights[e]
+                                pe = e
+                        else:
+                            pe = e
+                            break
+                elist.insert(0, pe)
+                vlist.insert(0, p)
+
+                if v in targets:
+                    targets_indexes[v] = len(elist) - 1
+                
+                v = p
+
+            # the paths are updated for every encountered target
+            path_len = len(elist)
+            vlists.update({tgt: vlist[:path_len-idx+1] for tgt, idx in targets_indexes.iteritems()})
+            elists.update({tgt: elist[:path_len-idx] for tgt, idx in targets_indexes.iteritems()})
+
+    if len(target) == 1:
+        return elists[target[0]], vlists[target[0]]
+    return [elists[tgt] for tgt in target], [vlists[tgt] for tgt in target]
 
 
 def pseudo_diameter(g, source=None, weights=None):
@@ -1670,6 +1700,7 @@ def make_maximal_planar(g, unfilter=False):
     g = GraphView(g, directed=False)
     libgraph_tool_topology.maximal_planar(g._Graph__graph)
 
+
 def is_DAG(g):
     """
     Return `True` if the graph is a directed acyclic graph (DAG).