Add Boykov Kolmogorov algorithm for maximum flow problems.

This algorithm has worse case complexity `O(n^2m|C|)` where |C| is the
cost of the minimum cut [1], but it is faster than other algorithms with
theoretically better worse case complexities for some problems, such as
vision and image processing, and maybe others.

For the case of connectivity problems in networks with higly skewed
degree distributions (eg social networks) this algorithm doesn't run
faster than Edmonds-Karp but after reading about it I thought that it
was really cool and decided to implement it.

This is still WIP as the marking heuristic described in [2] is
not yet implemented. According to the authors it speeds up the running
time considerably.

[1] Boykov, Y., & Kolmogorov, V. (2004). An experimental comparison of
min-cut/max-flow algorithms for energy minimization in vision.
http://www.csd.uwo.ca/~yuri/Papers/pami04.pdf
[2] Vladimir Kolmogorov. Graph-based Algorithms for Multi-camera
Reconstruction Problem. PhD thesis, Cornell University, CS Department,
2003. pp. 109-114.
https://pub.ist.ac.at/~vnk/papers/thesis.pdf

10 files changed, 341 insertions, 3 deletions

diff --git a/doc/source/reference/algorithms.flow.rst b/doc/source/reference/algorithms.flow.rst
index 8a813c0b..c8390f0a 100644
--- a/doc/source/reference/algorithms.flow.rst
+++ b/doc/source/reference/algorithms.flow.rst
@@ -48,6 +48,14 @@ Dinitz
    dinitz
 
 
+Boykov-Kolmogorov
+-----------------
+.. autosummary::
+   :toctree: generated/
+
+   boykov_kolmogorov
+
+
 Utils
 -----
 .. autosummary::
diff --git a/networkx/algorithms/connectivity/connectivity.py b/networkx/algorithms/connectivity/connectivity.py
index 6d18b462..b10afd63 100644
--- a/networkx/algorithms/connectivity/connectivity.py
+++ b/networkx/algorithms/connectivity/connectivity.py
@@ -10,6 +10,7 @@ from operator import itemgetter
 import networkx as nx
 # Define the default maximum flow function to use in all flow based
 # connectivity algorithms.
+from networkx.algorithms.flow import boykov_kolmogorov
 from networkx.algorithms.flow import dinitz
 from networkx.algorithms.flow import edmonds_karp
 from networkx.algorithms.flow import shortest_augmenting_path
@@ -203,6 +204,8 @@ def local_node_connectivity(G, s, t, flow_func=None, auxiliary=None,
         kwargs['cutoff'] = cutoff
     elif flow_func is dinitz:
         kwargs['cutoff'] = cutoff
+    elif flow_func is boykov_kolmogorov:
+        kwargs['cutoff'] = cutoff
 
     return nx.maximum_flow_value(H, '%sB' % mapping[s], '%sA' % mapping[t], **kwargs)
 
@@ -638,6 +641,8 @@ def local_edge_connectivity(G, u, v, flow_func=None, auxiliary=None,
         kwargs['cutoff'] = cutoff
     elif flow_func is dinitz:
         kwargs['cutoff'] = cutoff
+    elif flow_func is boykov_kolmogorov:
+        kwargs['cutoff'] = cutoff
 
     return nx.maximum_flow_value(H, u, v, **kwargs)
 
diff --git a/networkx/algorithms/connectivity/tests/test_connectivity.py b/networkx/algorithms/connectivity/tests/test_connectivity.py
index 33ea72dd..88f73b51 100644
--- a/networkx/algorithms/connectivity/tests/test_connectivity.py
+++ b/networkx/algorithms/connectivity/tests/test_connectivity.py
@@ -7,6 +7,7 @@ from networkx.algorithms.connectivity import local_edge_connectivity
 from networkx.algorithms.connectivity import local_node_connectivity
 
 flow_funcs = [
+    flow.boykov_kolmogorov,
     flow.dinitz,
     flow.edmonds_karp,
     flow.preflow_push,
diff --git a/networkx/algorithms/connectivity/tests/test_cuts.py b/networkx/algorithms/connectivity/tests/test_cuts.py
index 2d12c41b..b75e3e80 100644
--- a/networkx/algorithms/connectivity/tests/test_cuts.py
+++ b/networkx/algorithms/connectivity/tests/test_cuts.py
@@ -7,6 +7,7 @@ from networkx.algorithms.connectivity import minimum_st_node_cut
 from networkx.utils import arbitrary_element
 
 flow_funcs = [
+    flow.boykov_kolmogorov,
     flow.dinitz,
     flow.edmonds_karp,
     flow.preflow_push,
diff --git a/networkx/algorithms/connectivity/tests/test_kcutsets.py b/networkx/algorithms/connectivity/tests/test_kcutsets.py
index 66ae3204..0ecd5e63 100644
--- a/networkx/algorithms/connectivity/tests/test_kcutsets.py
+++ b/networkx/algorithms/connectivity/tests/test_kcutsets.py
@@ -9,6 +9,7 @@ from networkx.algorithms.connectivity.kcutsets import _is_separating_set
 
 
 flow_funcs = [
+    flow.boykov_kolmogorov,
     flow.dinitz,
     flow.edmonds_karp,
     flow.preflow_push,
diff --git a/networkx/algorithms/flow/__init__.py b/networkx/algorithms/flow/__init__.py
index ffb54cb1..561440bd 100644
--- a/networkx/algorithms/flow/__init__.py
+++ b/networkx/algorithms/flow/__init__.py
@@ -1,5 +1,6 @@
 from .maxflow import *
 from .mincost import *
+from .boykovkolmogorov import *
 from .dinitz_alg import *
 from .edmondskarp import *
 from .preflowpush import *
diff --git a/networkx/algorithms/flow/boykovkolmogorov.py b/networkx/algorithms/flow/boykovkolmogorov.py
new file mode 100644
index 00000000..736360cf
--- /dev/null
+++ b/networkx/algorithms/flow/boykovkolmogorov.py
@@ -0,0 +1,312 @@
+# boykovkolmogorov.py - Boykov Kolmogorov algorithm for maximum flow problems.
+#
+# Copyright 2016 NetworkX developers.
+#
+# This file is part of NetworkX.
+#
+# NetworkX is distributed under a BSD license; see LICENSE.txt for more
+# information.
+#
+# Author: Jordi Torrents <jordi.t21@gmail.com>
+"""
+Boykov-Kolmogorov algorithm for maximum flow problems.
+"""
+from collections import deque
+
+import networkx as nx
+from networkx.algorithms.flow.utils import build_residual_network
+
+__all__ = ['boykov_kolmogorov']
+
+def boykov_kolmogorov(G, s, t, capacity='capacity', residual=None, value_only=False, cutoff=None):
+    r"""Find a maximum single-commodity flow using Boykov-Kolmogorov algorithm.
+
+    This function returns the residual network resulting after computing
+    the maximum flow. See below for details about the conventions
+    NetworkX uses for defining residual networks.
+
+    This algorithm has worse case complexity `O(n^2 m |C|)` for `n` nodes, `m`
+    edges, and `|C|` the cost of the minimum cut [1]_.
+
+    Parameters
+    ----------
+    G : NetworkX graph
+        Edges of the graph are expected to have an attribute called
+        'capacity'. If this attribute is not present, the edge is
+        considered to have infinite capacity.
+
+    s : node
+        Source node for the flow.
+
+    t : node
+        Sink node for the flow.
+
+    capacity : string
+        Edges of the graph G are expected to have an attribute capacity
+        that indicates how much flow the edge can support. If this
+        attribute is not present, the edge is considered to have
+        infinite capacity. Default value: 'capacity'.
+
+    residual : NetworkX graph
+        Residual network on which the algorithm is to be executed. If None, a
+        new residual network is created. Default value: None.
+
+    value_only : bool
+        If True compute only the value of the maximum flow. This parameter
+        will be ignored by this algorithm because it is not applicable.
+
+    cutoff : integer, float
+        If specified, the algorithm will terminate when the flow value reaches
+        or exceeds the cutoff. In this case, it may be unable to immediately
+        determine a minimum cut. Default value: None.
+
+    Returns
+    -------
+    R : NetworkX DiGraph
+        Residual network after computing the maximum flow.
+
+    Raises
+    ------
+    NetworkXError
+        The algorithm does not support MultiGraph and MultiDiGraph. If
+        the input graph is an instance of one of these two classes, a
+        NetworkXError is raised.
+
+    NetworkXUnbounded
+        If the graph has a path of infinite capacity, the value of a
+        feasible flow on the graph is unbounded above and the function
+        raises a NetworkXUnbounded.
+
+    See also
+    --------
+    :meth:`maximum_flow`
+    :meth:`minimum_cut`
+    :meth:`preflow_push`
+    :meth:`shortest_augmenting_path`
+
+    Notes
+    -----
+    The residual network :samp:`R` from an input graph :samp:`G` has the
+    same nodes as :samp:`G`. :samp:`R` is a DiGraph that contains a pair
+    of edges :samp:`(u, v)` and :samp:`(v, u)` iff :samp:`(u, v)` is not a
+    self-loop, and at least one of :samp:`(u, v)` and :samp:`(v, u)` exists
+    in :samp:`G`.
+
+    For each edge :samp:`(u, v)` in :samp:`R`, :samp:`R[u][v]['capacity']`
+    is equal to the capacity of :samp:`(u, v)` in :samp:`G` if it exists
+    in :samp:`G` or zero otherwise. If the capacity is infinite,
+    :samp:`R[u][v]['capacity']` will have a high arbitrary finite value
+    that does not affect the solution of the problem. This value is stored in
+    :samp:`R.graph['inf']`. For each edge :samp:`(u, v)` in :samp:`R`,
+    :samp:`R[u][v]['flow']` represents the flow function of :samp:`(u, v)` and
+    satisfies :samp:`R[u][v]['flow'] == -R[v][u]['flow']`.
+
+    The flow value, defined as the total flow into :samp:`t`, the sink, is
+    stored in :samp:`R.graph['flow_value']`. If :samp:`cutoff` is not
+    specified, reachability to :samp:`t` using only edges :samp:`(u, v)` such
+    that :samp:`R[u][v]['flow'] < R[u][v]['capacity']` induces a minimum
+    :samp:`s`-:samp:`t` cut.
+
+    Examples
+    --------
+    >>> import networkx as nx
+    >>> from networkx.algorithms.flow import boykov_kolmogorov
+
+    The functions that implement flow algorithms and output a residual
+    network, such as this one, are not imported to the base NetworkX
+    namespace, so you have to explicitly import them from the flow package.
+
+    >>> G = nx.DiGraph()
+    >>> G.add_edge('x','a', capacity=3.0)
+    >>> G.add_edge('x','b', capacity=1.0)
+    >>> G.add_edge('a','c', capacity=3.0)
+    >>> G.add_edge('b','c', capacity=5.0)
+    >>> G.add_edge('b','d', capacity=4.0)
+    >>> G.add_edge('d','e', capacity=2.0)
+    >>> G.add_edge('c','y', capacity=2.0)
+    >>> G.add_edge('e','y', capacity=3.0)
+    >>> R = boykov_kolmogorov(G, 'x', 'y')
+    >>> flow_value = nx.maximum_flow_value(G, 'x', 'y')
+    >>> flow_value
+    3.0
+    >>> flow_value == R.graph['flow_value']
+    True
+
+    References
+    ----------
+    .. [1] Boykov, Y., & Kolmogorov, V. (2004). An experimental comparison
+           of min-cut/max-flow algorithms for energy minimization in vision.
+           Pattern Analysis and Machine Intelligence, IEEE Transactions on,
+           26(9), 1124-1137.
+
+    """
+    R = boykov_kolmogorov_impl(G, s, t, capacity, residual, cutoff)
+    R.graph['algorithm'] = 'boykov_kolmogorov'
+    return R
+
+
+def boykov_kolmogorov_impl(G, s, t, capacity, residual, cutoff):
+    if s not in G:
+        raise nx.NetworkXError('node %s not in graph' % str(s))
+    if t not in G:
+        raise nx.NetworkXError('node %s not in graph' % str(t))
+    if s == t:
+        raise nx.NetworkXError('source and sink are the same node')
+
+    if residual is None:
+        R = build_residual_network(G, capacity)
+    else:
+        R = residual
+
+    # Initialize/reset the residual network.
+    nx.set_edge_attributes(R, 'flow', 0)
+
+    # Use an arbitrary high value as infinite. It is computed
+    # when building the residual network.
+    INF = R.graph['inf']
+
+    if cutoff is None:
+        cutoff = INF
+
+    R_succ = R.succ
+    R_pred = R.pred
+
+
+    def grow():
+        """Bidirectional breadth-first search for the growth stage.
+
+           Returns a connecting edge, that is and edge that connects
+           a node from the source search tree with a node from the
+           target search tree.
+           The first node in the connecting edge is always from the
+           source tree and the last node from the target tree.
+        """
+        while active:
+            u = active[0]
+            if u in source_tree:
+                for v, attr in R_succ[u].items():
+                    if v not in source_tree and attr['capacity'] - attr['flow'] > 0:
+                        if v in target_tree:
+                            return u, v
+                        source_tree[v] = u
+                        active.append(v)
+                _ = active.popleft()
+            elif u in target_tree:
+                for v, attr in R_pred[u].items():
+                    if v not in target_tree and attr['capacity'] - attr['flow'] > 0:
+                        if v in source_tree:
+                            return v, u
+                        target_tree[v] = u
+                        active.append(v)
+                _ = active.popleft()
+        return None, None
+
+
+    def augment(u, v):
+        """Augmentation stage.
+        """
+        # Reconstruct path and determine its residual capacity.
+        # We start from a connecting edge, which links a node
+        # from the source tree to a node from the target tree.
+        # The connecting edge is the output of the grow function.
+        attr = R_succ[u][v]
+        flow = min(INF, attr['capacity'] - attr['flow'])
+        path = [u]
+        # Trace a path from u to s in source_tree.
+        w = u
+        while w != s:
+            n = w
+            w = source_tree[n]
+            attr = R_pred[n][w]
+            flow = min(flow, attr['capacity'] - attr['flow'])
+            path.append(w)
+        path.reverse()
+        # Trace a path from v to t in target_tree.
+        path.append(v)
+        w = v
+        while w != t:
+            n = w
+            w = target_tree[n]
+            attr = R_succ[n][w]
+            flow = min(flow, attr['capacity'] - attr['flow'])
+            path.append(w)
+        # Augment flow along the path and check for saturated edges.
+        it = iter(path)
+        u = next(it)
+        for v in it:
+            R_succ[u][v]['flow'] += flow
+            R_succ[v][u]['flow'] -= flow
+            if R_succ[u][v]['flow'] == R_succ[u][v]['capacity']:
+                if v in source_tree:
+                    source_tree[v] = None
+                    orphans.append(v)
+                if u in target_tree:
+                    target_tree[u] = None
+                    orphans.append(u)
+            u = v
+        return flow
+
+
+    def adopt():
+        """Reconstruct the trees by adopting or discarding orphans.
+        """
+        while orphans:
+            u = orphans.popleft()
+            tree, neighbors = _get_tree_and_neighbors(u)
+            nbrs = ((n, attr) for n, attr in neighbors[u].items() if n in tree)
+            for v, attr in nbrs:
+                if attr['capacity'] - attr['flow'] > 0:
+                    if _has_valid_root(v, tree):
+                        tree[u] = v
+                        break
+            else:
+                nbrs = ((n, attr) for n, attr in neighbors[u].items() if n in tree)
+                for v, attr in nbrs:
+                    if attr['capacity'] - attr['flow'] > 0:
+                        if v not in active:
+                            active.append(v)
+                    if tree[v] == u:
+                        tree[v] = None
+                        orphans.appendleft(v)
+                if u in active:
+                    active.remove(u)
+                del tree[u]
+
+
+    def _has_valid_root(n, tree):
+        v = n
+        while v is not None:
+            if v == s or v == t:
+                return True
+            v = tree[v]
+        return False
+
+
+    def _get_tree_and_neighbors(n):
+        if n in source_tree:
+            return source_tree, R_pred
+        elif n in target_tree:
+            return target_tree, R_succ
+
+
+    source_tree = {s: None}
+    target_tree = {t: None}
+    active = deque([s, t])
+    orphans = deque()
+    flow_value = 0
+    while flow_value < cutoff:
+        # Growth stage
+        u, v = grow()
+        if u is None:
+            break
+        # Augmentation stage
+        flow_value += augment(u, v)
+        # Adoption stage
+        adopt()
+
+
+    if flow_value * 2 > INF:
+        raise nx.NetworkXUnbounded('Infinite capacity path, flow unbounded above.')
+
+    R.graph['flow_value'] = flow_value
+    return R
diff --git a/networkx/algorithms/flow/maxflow.py b/networkx/algorithms/flow/maxflow.py
index 9cccf011..9d31b517 100644
--- a/networkx/algorithms/flow/maxflow.py
+++ b/networkx/algorithms/flow/maxflow.py
@@ -4,6 +4,7 @@ Maximum flow (and minimum cut) algorithms on capacitated graphs.
 """
 import networkx as nx
 
+from .boykovkolmogorov import boykov_kolmogorov
 from .dinitz_alg import dinitz
 from .edmondskarp import edmonds_karp
 from .preflowpush import preflow_push
@@ -12,7 +13,13 @@ from .utils import build_flow_dict
 # Define the default flow function for computing maximum flow.
 default_flow_func = preflow_push
 # Functions that don't support cutoff for minimum cut computations.
-flow_funcs = [dinitz, edmonds_karp, preflow_push, shortest_augmenting_path]
+flow_funcs = [
+    boykov_kolmogorov,
+    dinitz,
+    edmonds_karp,
+    preflow_push,
+    shortest_augmenting_path,
+]
 
 __all__ = ['maximum_flow',
            'maximum_flow_value',
diff --git a/networkx/algorithms/flow/tests/test_maxflow.py b/networkx/algorithms/flow/tests/test_maxflow.py
index 30c0603a..fc1bbbcd 100644
--- a/networkx/algorithms/flow/tests/test_maxflow.py
+++ b/networkx/algorithms/flow/tests/test_maxflow.py
@@ -5,12 +5,13 @@ from nose.tools import *
 
 import networkx as nx
 from networkx.algorithms.flow import build_flow_dict, build_residual_network
+from networkx.algorithms.flow import boykov_kolmogorov
 from networkx.algorithms.flow import edmonds_karp
 from networkx.algorithms.flow import preflow_push
 from networkx.algorithms.flow import shortest_augmenting_path
 from networkx.algorithms.flow import dinitz
 
-flow_funcs = [dinitz, edmonds_karp, preflow_push, shortest_augmenting_path]
+flow_funcs = [boykov_kolmogorov, dinitz, edmonds_karp, preflow_push, shortest_augmenting_path]
 max_min_funcs = [nx.maximum_flow, nx.minimum_cut]
 flow_value_funcs = [nx.maximum_flow_value, nx.minimum_cut_value]
 interface_funcs = sum([max_min_funcs, flow_value_funcs], [])
diff --git a/networkx/algorithms/flow/tests/test_maxflow_large_graph.py b/networkx/algorithms/flow/tests/test_maxflow_large_graph.py
index 89e4f595..f9db6778 100644
--- a/networkx/algorithms/flow/tests/test_maxflow_large_graph.py
+++ b/networkx/algorithms/flow/tests/test_maxflow_large_graph.py
@@ -11,6 +11,7 @@ from nose.tools import *
 
 import networkx as nx
 from networkx.algorithms.flow import build_flow_dict, build_residual_network
+from networkx.algorithms.flow import boykov_kolmogorov
 from networkx.algorithms.flow import edmonds_karp
 from networkx.algorithms.flow import preflow_push
 from networkx.algorithms.flow import shortest_augmenting_path
@@ -18,7 +19,7 @@ from networkx.algorithms.flow import dinitz
 
 # Dinitz algorithm is too slow in big and dense networks to test it here.
 # It is alredy tested in test_maxflow.py
-flow_funcs = [edmonds_karp, preflow_push, shortest_augmenting_path]
+flow_funcs = [boykov_kolmogorov, edmonds_karp, preflow_push, shortest_augmenting_path]
 
 msg = "Assertion failed in function: {0}"