Add shortest_simple_paths function.

This function is based on Yen's algorithm for finding k shortest paths.
It generates all simple paths between source and target starting from
the shortest one. See #762 and #793 for prior work on this and discussion.

This implementation is based on Andrey's work on #762. It generalizes his
approach there to the weighted case, as in Greg's implementation posted
on #793.

Regarding the discussion on #762 about how to do the filtering of nodes
and edges efficiently, this PR adds private modifications of the functions
`bidirectional_shortest_path` and `bidirectional_dijkstra` (that accept
containers with nodes and edges to exclude during the SP search) for its
exclusive use.

This still needs work, especially on documentation and tests. I'm pushing
it so we can have the discussion with the code available.

One thing to discuss is whether or not to include an one liner function
to get the k shortest paths. Could be something like this:

```python
from itertools import islice
def k_shortest_paths(G, source, target, k, weight=None):
        return list(islice(nx.shortest_simple_paths(G, source, target, weight=weight), k))
```

I think that, given that `all_simple_paths` is way faster in computing
all paths, the principal use of `shortest_simple_paths` will be to get
the k best/shortest paths. So it will be a more user friendly interface
if we add it. However, we could also add this function as example in the
`shortest_simple_paths` docstring. I'm ok either way. Thoughts?

1 files changed, 473 insertions, 2 deletions

diff --git a/networkx/algorithms/simple_paths.py b/networkx/algorithms/simple_paths.py
index f72d4d25..9f02bfb1 100644
--- a/networkx/algorithms/simple_paths.py
+++ b/networkx/algorithms/simple_paths.py
@@ -3,10 +3,21 @@
 #    Sergio Nery Simoes <sergionery@gmail.com>
 #    All rights reserved.
 #    BSD license.
+from heapq import heappush, heappop
+from itertools import count
+
 import networkx as nx
+
 __author__ = """\n""".join(['Sérgio Nery Simões <sergionery@gmail.com>',
-                            'Aric Hagberg <aric.hagberg@gmail.com>'])
-__all__ = ['all_simple_paths']
+                            'Aric Hagberg <aric.hagberg@gmail.com>',
+                            'Andrey Paramonov',
+                            'Jordi Torrents <jordi.t21@gmail.com>'])
+
+__all__ = [
+    'all_simple_paths',
+    'shortest_simple_paths',
+]
+
 
 def all_simple_paths(G, source, target, cutoff=None):
     """Generate all simple paths in the graph G from source to target.
@@ -75,6 +86,7 @@ def all_simple_paths(G, source, target, cutoff=None):
     else:
         return _all_simple_paths_graph(G, source, target, cutoff=cutoff)
 
+
 def _all_simple_paths_graph(G, source, target, cutoff=None):
     if cutoff < 1:
         return
@@ -122,3 +134,462 @@ def _all_simple_paths_multigraph(G, source, target, cutoff=None):
                 yield visited + [target]
             stack.pop()
             visited.pop()
+
+
+def shortest_simple_paths(G, source, target, weight=None):
+    """Generate all simple paths in the graph G from source to target,
+    starting from shortest ones.
+
+    A simple path is a path with no repeated nodes.
+
+    If a weighted shortest path search is to be used, no negative weights
+    are allawed.
+
+    Parameters
+    ----------
+    G : NetworkX graph
+
+    source : node
+       Starting node for path
+
+    target : node
+       Ending node for path
+
+    weight : string
+        Name of the edge attribute to be used as a weight. If None all
+        edges are considered to have unit weight. Default value None.
+
+    Returns
+    -------
+    path_generator: generator
+       A generator that produces lists of simple paths, in order from
+       shortest to longest.
+
+    Notes
+    -----
+    This procedure is based on algorithm by Jin Y. Yen [1]_.  Finding
+    the first K paths requires O(KN^3) operations.
+
+    References
+    ----------
+    .. [1] Jin Y. Yen, "Finding the K Shortest Loopless Paths in a
+       Network", Management Science, Vol. 17, No. 11, Theory Series
+       (Jul., 1971), pp. 712-716.
+
+    See Also
+    --------
+    all_shortest_paths, shortest_path
+
+    """
+    if weight is None:
+        length_func = len
+        shortest_path_func = _bidirectional_shortest_path
+    else:
+        def length_func(path):
+            return sum(G.edge[u][v][weight] for (u, v) in zip(path, path[1:]))
+        shortest_path_func = _bidirectional_dijkstra
+
+    listA = list()
+    listB = PathBuffer()
+    prev_path = None
+    while True:
+        if not prev_path:
+            length, path = shortest_path_func(G, source, target)
+            listB.push(length, path)
+        else:
+            ignore_nodes = set()
+            ignore_edges = set()
+            for i in range(1, len(prev_path)):
+                root = prev_path[:i]
+                for path in listA:
+                    if path[:i] == root:
+                        ignore_edges.add((path[i-1], path[i]))
+                ignore_nodes.add(root[-1])
+                try:
+                    length, spur = shortest_path_func(G, root[-1], target,
+                                                      ignore_nodes=ignore_nodes,
+                                                      ignore_edges=ignore_edges)
+                    assert root[-1] == spur[0]
+                    assert not set(root[:-1]).intersection(set(spur))
+                    path = root[:-1] + spur
+                    listB.push(length_func(path), path)
+                except nx.NetworkXNoPath:
+                    pass
+
+        if listB:
+            path = listB.pop()
+            assert path not in listA
+            yield path
+            listA.append(path)
+            prev_path = path
+        else:
+            break
+
+
+class PathBuffer:
+    def __init__(self):
+        self.paths = set()
+        self.sortedpaths = list()
+        self.counter = count()
+
+    def __len__(self):
+        return len(self.sortedpaths)
+
+    def push(self, cost, path):
+        hashable_path = tuple(path)
+        if hashable_path not in self.paths:
+            heappush(self.sortedpaths, (cost, next(self.counter), path))
+            self.paths.add(hashable_path)
+
+    def pop(self):
+        (cost, num, path) = heappop(self.sortedpaths)
+        hashable_path = tuple(path)
+        self.paths.remove(hashable_path)
+        return path
+
+
+def _bidirectional_shortest_path(G, source, target,
+                                 ignore_nodes=None,
+                                 ignore_edges=None):
+    """Return the shortest path between source and target ignoring nodes and 
+       edges in the containers ignore_nodes and ignore_edges.
+
+    This is a custom modification of the standard bidirectional shortest path implementation
+    at networkx.algorithms.unweighted
+
+    Parameters
+    ----------
+    G : NetworkX graph
+
+    source : node
+       starting node for path
+
+    target : node
+       ending node for path
+
+    ignore_nodes : container of nodes
+       nodes to ignore, optional
+
+    ignore_edges : container of edges
+       edges to ignore, optional
+
+    Returns
+    -------
+    path: list
+       List of nodes in a path from source to target.
+
+    Raises
+    ------
+    NetworkXNoPath
+       If no path exists between source and target.
+
+    See Also
+    --------
+    shortest_path
+
+    """
+    # call helper to do the real work
+    results=_bidirectional_pred_succ(G,source,target,ignore_nodes,ignore_edges)
+    pred,succ,w=results
+
+    # build path from pred+w+succ
+    path=[]
+    # from w to target
+    while w is not None:
+        path.append(w)
+        w=succ[w]
+    # from source to w
+    w=pred[path[0]]
+    while w is not None:
+        path.insert(0,w)
+        w=pred[w]
+
+    return len(path), path
+
+
+def _bidirectional_pred_succ(G, source, target, ignore_nodes=None, ignore_edges=None):
+    """Bidirectional shortest path helper.
+       Returns (pred,succ,w) where
+       pred is a dictionary of predecessors from w to the source, and
+       succ is a dictionary of successors from w to the target.
+    """
+    # does BFS from both source and target and meets in the middle
+    if target == source:
+        return ({target:None},{source:None},source)
+
+    # handle either directed or undirected
+    if G.is_directed():
+        Gpred=G.predecessors_iter
+        Gsucc=G.successors_iter
+    else:
+        Gpred=G.neighbors_iter
+        Gsucc=G.neighbors_iter
+
+    # support optional nodes filter
+    if ignore_nodes:
+        def filter_iter(nodes_iter):
+            def iterate(v):
+                for w in nodes_iter(v):
+                    if w not in ignore_nodes:
+                        yield w
+            return iterate
+
+        Gpred=filter_iter(Gpred)
+        Gsucc=filter_iter(Gsucc)
+
+    # support optional edges filter
+    if ignore_edges:
+        if G.is_directed():
+            def filter_pred_iter(pred_iter):
+                def iterate(v):
+                    for w in pred_iter(v):
+                        if (w, v) not in ignore_edges:
+                            yield w
+                return iterate
+
+            def filter_succ_iter(succ_iter):
+                def iterate(v):
+                    for w in succ_iter(v):
+                        if (v, w) not in ignore_edges:
+                            yield w
+                return iterate
+
+            Gpred=filter_pred_iter(Gpred)
+            Gsucc=filter_succ_iter(Gsucc)
+
+        else:
+            def filter_iter(nodes_iter):
+                def iterate(v):
+                    for w in nodes_iter(v):
+                        if (v, w) not in ignore_edges \
+                                and (w, v) not in ignore_edges:
+                            yield w
+                return iterate
+
+            Gpred=filter_iter(Gpred)
+            Gsucc=filter_iter(Gsucc)
+
+    # predecesssor and successors in search
+    pred={source:None}
+    succ={target:None}
+
+    # initialize fringes, start with forward
+    forward_fringe=[source]
+    reverse_fringe=[target]
+
+    while forward_fringe and reverse_fringe:
+        if len(forward_fringe) <= len(reverse_fringe):
+            this_level=forward_fringe
+            forward_fringe=[]
+            for v in this_level:
+                for w in Gsucc(v):
+                    if w not in pred:
+                        forward_fringe.append(w)
+                        pred[w]=v
+                    if w in succ:
+                        # found path
+                        return pred,succ,w
+        else:
+            this_level=reverse_fringe
+            reverse_fringe=[]
+            for v in this_level:
+                for w in Gpred(v):
+                    if w not in succ:
+                        succ[w]=v
+                        reverse_fringe.append(w)
+                    if w in pred:
+                        # found path
+                        return pred,succ,w
+
+    raise nx.NetworkXNoPath("No path between %s and %s." % (source, target))
+
+
+def _bidirectional_dijkstra(G, source, target, weight='weight',
+                            ignore_nodes=None, ignore_edges=None):
+    """Dijkstra's algorithm for shortest paths using bidirectional search.
+
+    Parameters
+    ----------
+    G : NetworkX graph
+
+    source : node
+       Starting node.
+
+    target : node
+       Ending node.
+
+    weight: string, optional (default='weight')
+       Edge data key corresponding to the edge weight
+
+    Returns
+    -------
+    length : number
+        Shortest path length.
+
+    Returns a tuple of two dictionaries keyed by node.
+    The first dictionary stores distance from the source.
+    The second stores the path from the source to that node.
+
+    Raises
+    ------
+    NetworkXNoPath
+        If no path exists between source and target.
+
+    Examples
+    --------
+    >>> G=nx.path_graph(5)
+    >>> length,path=nx.bidirectional_dijkstra(G,0,4)
+    >>> print(length)
+    4
+    >>> print(path)
+    [0, 1, 2, 3, 4]
+
+    Notes
+    -----
+    Edge weight attributes must be numerical.
+    Distances are calculated as sums of weighted edges traversed.
+
+    In practice  bidirectional Dijkstra is much more than twice as fast as
+    ordinary Dijkstra.
+
+    Ordinary Dijkstra expands nodes in a sphere-like manner from the
+    source. The radius of this sphere will eventually be the length
+    of the shortest path. Bidirectional Dijkstra will expand nodes
+    from both the source and the target, making two spheres of half
+    this radius. Volume of the first sphere is pi*r*r while the
+    others are 2*pi*r/2*r/2, making up half the volume.
+
+    This algorithm is not guaranteed to work if edge weights
+    are negative or are floating point numbers
+    (overflows and roundoff errors can cause problems).
+
+    See Also
+    --------
+    shortest_path
+    shortest_path_length
+    """
+    if source == target:
+        return (0, [source])
+
+    # handle either directed or undirected
+    if G.is_directed():
+        Gpred=G.predecessors_iter
+        Gsucc=G.successors_iter
+    else:
+        Gpred=G.neighbors_iter
+        Gsucc=G.neighbors_iter
+
+    # support optional nodes filter
+    if ignore_nodes:
+        def filter_iter(nodes_iter):
+            def iterate(v):
+                for w in nodes_iter(v):
+                    if w not in ignore_nodes:
+                        yield w
+            return iterate
+
+        Gpred=filter_iter(Gpred)
+        Gsucc=filter_iter(Gsucc)
+
+    # support optional edges filter
+    if ignore_edges:
+        if G.is_directed():
+            def filter_pred_iter(pred_iter):
+                def iterate(v):
+                    for w in pred_iter(v):
+                        if (w, v) not in ignore_edges:
+                            yield w
+                return iterate
+
+            def filter_succ_iter(succ_iter):
+                def iterate(v):
+                    for w in succ_iter(v):
+                        if (v, w) not in ignore_edges:
+                            yield w
+                return iterate
+
+            Gpred=filter_pred_iter(Gpred)
+            Gsucc=filter_succ_iter(Gsucc)
+
+        else:
+            def filter_iter(nodes_iter):
+                def iterate(v):
+                    for w in nodes_iter(v):
+                        if (v, w) not in ignore_edges \
+                                and (w, v) not in ignore_edges:
+                            yield w
+                return iterate
+
+            Gpred=filter_iter(Gpred)
+            Gsucc=filter_iter(Gsucc)
+
+    push = heappush
+    pop = heappop
+    # Init:   Forward             Backward
+    dists  = [{},                {}]  # dictionary of final distances
+    paths  = [{source: [source]}, {target: [target]}]  # dictionary of paths
+    fringe = [[],                []]  # heap of (distance, node) tuples for
+                                      # extracting next node to expand
+    seen   = [{source: 0},        {target: 0}]  # dictionary of distances to
+                                                # nodes seen
+    c = count()
+    # initialize fringe heap
+    push(fringe[0], (0, next(c), source))
+    push(fringe[1], (0, next(c), target))
+    # neighs for extracting correct neighbor information
+    neighs = [Gsucc, Gpred]
+    # variables to hold shortest discovered path
+    #finaldist = 1e30000
+    finalpath = []
+    dir = 1
+    while fringe[0] and fringe[1]:
+        # choose direction
+        # dir == 0 is forward direction and dir == 1 is back
+        dir = 1 - dir
+        # extract closest to expand
+        (dist, _, v) = pop(fringe[dir])
+        if v in dists[dir]:
+            # Shortest path to v has already been found
+            continue
+        # update distance
+        dists[dir][v] = dist  # equal to seen[dir][v]
+        if v in dists[1 - dir]:
+            # if we have scanned v in both directions we are done
+            # we have now discovered the shortest path
+            return (finaldist, finalpath)
+
+        for w in neighs[dir](v):
+            if(dir == 0):  # forward
+                if G.is_multigraph():
+                    minweight = min((dd.get(weight, 1)
+                                     for k, dd in G[v][w].items()))
+                else:
+                    minweight = G[v][w].get(weight, 1)
+                vwLength = dists[dir][v] + minweight  # G[v][w].get(weight,1)
+            else:  # back, must remember to change v,w->w,v
+                if G.is_multigraph():
+                    minweight = min((dd.get(weight, 1)
+                                     for k, dd in G[w][v].items()))
+                else:
+                    minweight = G[w][v].get(weight, 1)
+                vwLength = dists[dir][v] + minweight  # G[w][v].get(weight,1)
+
+            if w in dists[dir]:
+                if vwLength < dists[dir][w]:
+                    raise ValueError(
+                        "Contradictory paths found: negative weights?")
+            elif w not in seen[dir] or vwLength < seen[dir][w]:
+                # relaxing
+                seen[dir][w] = vwLength
+                push(fringe[dir], (vwLength, next(c), w))
+                paths[dir][w] = paths[dir][v] + [w]
+                if w in seen[0] and w in seen[1]:
+                    # see if this path is better than than the already
+                    # discovered shortest path
+                    totaldist = seen[0][w] + seen[1][w]
+                    if finalpath == [] or finaldist > totaldist:
+                        finaldist = totaldist
+                        revpath = paths[1][w][:]
+                        revpath.reverse()
+                        finalpath = paths[0][w] + revpath[1:]
+    raise nx.NetworkXNoPath("No path between %s and %s." % (source, target))