Add two-phase option to shortest augmenting path maxflow algorithm

networkx/algorithms/flow/shortest_augmenting_path.py

@@ -18,35 +18,29 @@
            'shortest_augmenting_path_flow']
 
 
-def shortest_augmenting_path_impl(G, s, t, capacity):
+def shortest_augmenting_path_impl(G, s, t, capacity, two_phase):
     """Implementation of the shortest augmenting path algorithm.
     """
     R = _build_residual_network(G, s, t, capacity)
 
-    def reverse_bfs():
-        """Perform a reverse breadth-first search from t in the residual
-        network.
-        """
-        heights = {t: 0}
-        q = deque([(t, 0)])
-        while q:
-            u, height = q.popleft()
-            height += 1
-            for v, u, attr in R.in_edges_iter(u, data=True):
-                if attr['flow'] < attr['capacity'] and v not in heights:
-                    heights[v] = height
-                    q.append((v, height))
-        return heights
-
     # Initialize heights of the nodes.
-    heights = reverse_bfs()
+    heights = {t: 0}
+    q = deque([(t, 0)])
+    while q:
+        u, height = q.popleft()
+        height += 1
+        for v, u, attr in R.in_edges_iter(u, data=True):
+            if attr['flow'] < attr['capacity'] and v not in heights:
+                heights[v] = height
+                q.append((v, height))
 
     if s not in heights:
         # t is not reachable from s in the residual network. The maximum flow
         # must be zero.
         return R
 
     n = len(G)
+    m = R.size() / 2
 
     # Initialize heights and 'current edge' data structures of the nodes.
     for u in R:
@@ -90,9 +84,13 @@ def relabel(u):
                 height = min(height, R.node[v]['height'])
         return height + 1
 
+    # Phase 1: Look for shortest augmenting paths using depth-first search.
+
     path = [s]
     u = s
-    while True:
+    d = n if not two_phase else int(min(m ** 0.5, 2 * n ** (2. / 3)))
+    done = False
+    while not done:
         height = R.node[u]['height']
         curr_edge = R.node[u]['curr_edge']
         # Depth-first search for the next node on the path to t.
@@ -114,10 +112,15 @@ def relabel(u):
                     # can now be terminated.
                     return R
                 height = relabel(u)
-                if u == s and height >= n:
-                    # t is disconnected from s in the residual network. No more
-                    # augmenting paths exist.
-                    return R
+                if u == s and height >= d:
+                    if not two_phase:
+                        # t is disconnected from s in the residual network. No
+                        # more augmenting paths exist.
+                        return R
+                    else:
+                        # t is at least d steps away from s. End of phase 1.
+                        done = True
+                        break
                 counts[height] += 1
                 R.node[u]['height'] = height
                 if u != s:
@@ -133,8 +136,36 @@ def relabel(u):
             path = [s]
             u = s
 
+    # Phase 2: Look for shortest augmenting paths using breadth-first search.
+    while True:
+        pred = {s: None}
+        q = deque([s])
+        done = False
+        while not done and q:
+            u = q.popleft()
+            for v, attr in R[u].items():
+                if v not in pred and attr['flow'] < attr['capacity']:
+                    pred[v] = u
+                    if v == t:
+                        done = True
+                        break
+                    q.append(v)
+        if not done:
+            # No augmenting path found.
+            return R
+        # Trace a s-t path using the predecessor array.
+        path = [t]
+        u = t
+        while True:
+            u = pred[u]
+            path.append(u)
+            if u == s:
+                break
+        path.reverse()
+        augment(path)
+
 
-def shortest_augmenting_path(G, s, t, capacity='capacity'):
+def shortest_augmenting_path(G, s, t, capacity='capacity', two_phase=False):
     """Find a maximum single-commodity flow using the shortest augmenting path
     algorithm.
 
@@ -161,6 +192,11 @@ def shortest_augmenting_path(G, s, t, capacity='capacity'):
         attribute is not present, the edge is considered to have
         infinite capacity. Default value: 'capacity'.
 
+    two_phase : bool
+        If True, a two-phase variant is used. The two-phase variant improves
+        the running time on unit-capacity networks from `O(nm)` to
+        `O(\min(n^{2/3}, m^{1/2}) m)`. Default value: False.
+
     Returns
     -------
     flow_value : integer, float
@@ -198,11 +234,12 @@ def shortest_augmenting_path(G, s, t, capacity='capacity'):
     >>> flow, F['a']['c']
     (3.0, 2.0)
     """
-    R = shortest_augmenting_path_impl(G, s, t, capacity)
+    R = shortest_augmenting_path_impl(G, s, t, capacity, two_phase)
     return (R.node[t]['excess'], _build_flow_dict(G, R))
 
 
-def shortest_augmenting_path_value(G, s, t, capacity='capacity'):
+def shortest_augmenting_path_value(G, s, t, capacity='capacity',
+                                   two_phase=False):
     """Find a maximum single-commodity flow using the shortest augmenting path
     algorithm.
 
@@ -229,6 +266,11 @@ def shortest_augmenting_path_value(G, s, t, capacity='capacity'):
         attribute is not present, the edge is considered to have
         infinite capacity. Default value: 'capacity'.
 
+    two_phase : bool
+        If True, a two-phase variant is used. The two-phase variant improves
+        the running time on unit-capacity networks from `O(nm)` to
+        `O(\min(n^{2/3}, m^{1/2}) m)`. Default value: False.
+
     Returns
     -------
     flow_value : integer, float
@@ -262,11 +304,12 @@ def shortest_augmenting_path_value(G, s, t, capacity='capacity'):
     >>> flow
     3.0
     """
-    R = shortest_augmenting_path_impl(G, s, t, capacity)
+    R = shortest_augmenting_path_impl(G, s, t, capacity, two_phase)
     return R.node[t]['excess']
 
 
-def shortest_augmenting_path_flow(G, s, t, capacity='capacity'):
+def shortest_augmenting_path_flow(G, s, t, capacity='capacity',
+                                  two_phase=False):
     """Find a maximum single-commodity flow using the shortest augmenting path
     algorithm.
 
@@ -293,6 +336,11 @@ def shortest_augmenting_path_flow(G, s, t, capacity='capacity'):
         attribute is not present, the edge is considered to have
         infinite capacity. Default value: 'capacity'.
 
+    two_phase : bool
+        If True, a two-phase variant is used. The two-phase variant improves
+        the running time on unit-capacity networks from `O(nm)` to
+        `O(\min(n^{2/3}, m^{1/2}) m)`. Default value: False.
+
     Returns
     -------
     flow_dict : dictionary
@@ -327,5 +375,5 @@ def shortest_augmenting_path_flow(G, s, t, capacity='capacity'):
     >>> F['a']['c']
     2.0
     """
-    R = shortest_augmenting_path_impl(G, s, t, capacity)
+    R = shortest_augmenting_path_impl(G, s, t, capacity, two_phase)
     return _build_flow_dict(G, R)

networkx/algorithms/flow/tests/test_maxflow.py

@@ -42,7 +42,12 @@ def compare_flows(G, s, t, solnFlows, solnValue, capacity = 'capacity'):
     flowValue, flowDict = nx.preflow_push(G, s, t, capacity)
     assert_equal(flowValue, solnValue)
     validate_flows(G, s, t, flowDict, solnValue, capacity)
-    flowValue, flowDict = nx.shortest_augmenting_path(G, s, t, capacity)
+    flowValue, flowDict = nx.shortest_augmenting_path(G, s, t, capacity,
+                                                      two_phase=False)
+    assert_equal(flowValue, solnValue)
+    validate_flows(G, s, t, flowDict, solnValue, capacity)
+    flowValue, flowDict = nx.shortest_augmenting_path(G, s, t, capacity,
+                                                      two_phase=True)
     assert_equal(flowValue, solnValue)
     validate_flows(G, s, t, flowDict, solnValue, capacity)
     assert_equal(nx.min_cut(G, s, t, capacity), solnValue)

networkx/algorithms/flow/tests/test_maxflow_large_graph.py

@@ -72,15 +72,21 @@ def test_complete_graph(self):
             G[u][v]['capacity'] = 5
         assert_equal(nx.ford_fulkerson(G, 1, 2)[0], 5 * (N - 1))
         assert_equal(nx.preflow_push(G, 1, 2)[0], 5 * (N - 1))
-        assert_equal(nx.shortest_augmenting_path(G, 1, 2)[0], 5 * (N - 1))
+        assert_equal(nx.shortest_augmenting_path(G, 1, 2, two_phase=False)[0],
+                     5 * (N - 1))
+        assert_equal(nx.shortest_augmenting_path(G, 1, 2, two_phase=True)[0],
+                     5 * (N - 1))
 
     def test_pyramid(self):
         N = 10
 #        N = 100 # this gives a graph with 5051 nodes
         G = gen_pyramid(N)
         assert_almost_equal(nx.ford_fulkerson(G, (0, 0), 't')[0], 1.)
         assert_almost_equal(nx.preflow_push(G, (0, 0), 't')[0], 1.)
-        assert_almost_equal(nx.shortest_augmenting_path(G, (0, 0), 't')[0], 1.)
+        assert_almost_equal(nx.shortest_augmenting_path(
+            G, (0, 0), 't', two_phase=False)[0], 1.)
+        assert_almost_equal(nx.shortest_augmenting_path(
+            G, (0, 0), 't', two_phase=True)[0], 1.)
 
     def test_gl1(self):
         G = read_graph('gl1')
@@ -89,27 +95,45 @@ def test_gl1(self):
         validate_flows(G, s, t, 156545, *nx.ford_fulkerson(G, s, t))
         validate_flows(G, s, t, 156545, nx.preflow_push_value(G, s, t),
                        nx.preflow_push_flow(G, s, t))
+        validate_flows(
+            G, s, t, 156545,
+            nx.shortest_augmenting_path_value(G, s, t, two_phase=False),
+            nx.shortest_augmenting_path_flow(G, s, t, two_phase=False))
         validate_flows(G, s, t, 156545,
-                       nx.shortest_augmenting_path_value(G, s, t),
-                       nx.shortest_augmenting_path_flow(G, s, t))
+            nx.shortest_augmenting_path_value(G, s, t, two_phase=True),
+            nx.shortest_augmenting_path_flow(G, s, t, two_phase=True))
 
     def test_gw1(self):
         G = read_graph('gw1')
         s = 1
         t = len(G)
         validate_flows(G, s, t, 1202018, *nx.ford_fulkerson(G, s, t))
-        validate_flows(G, s, t, 1202018,
-                       nx.shortest_augmenting_path_value(G, s, t),
-                       nx.shortest_augmenting_path_flow(G, s, t))
+        validate_flows(G, s, t, 1202018, nx.preflow_push_value(G, s, t),
+                       nx.preflow_push_flow(G, s, t))
+        validate_flows(
+            G, s, t, 1202018,
+            nx.shortest_augmenting_path_value(G, s, t, two_phase=False),
+            nx.shortest_augmenting_path_flow(G, s, t, two_phase=False))
+        validate_flows(
+            G, s, t, 1202018,
+            nx.shortest_augmenting_path_value(G, s, t, two_phase=True),
+            nx.shortest_augmenting_path_flow(G, s, t, two_phase=True))
 
     def test_wlm3(self):
         G = read_graph('wlm3')
         s = 1
         t = len(G)
         validate_flows(G, s, t, 11875108, *nx.ford_fulkerson(G, s, t))
-        validate_flows(G, s, t, 11875108,
-                       nx.shortest_augmenting_path_value(G, s, t),
-                       nx.shortest_augmenting_path_flow(G, s, t))
+        validate_flows(G, s, t, 11875108, nx.preflow_push_value(G, s, t),
+                       nx.preflow_push_flow(G, s, t))
+        validate_flows(
+            G, s, t, 11875108,
+            nx.shortest_augmenting_path_value(G, s, t, two_phase=False),
+            nx.shortest_augmenting_path_flow(G, s, t, two_phase=False))
+        validate_flows(
+            G, s, t, 11875108,
+            nx.shortest_augmenting_path_value(G, s, t, two_phase=True),
+            nx.shortest_augmenting_path_flow(G, s, t, two_phase=True))
 
     def test_preflow_push_global_relabel(self):
         G = read_graph('gw1')