Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2Highest-label preflow-push algorithm for maximum flow problems.
3"""
5from collections import deque
6from itertools import islice
8import networkx as nx
10from ...utils import arbitrary_element
11from .utils import (
12 CurrentEdge,
13 GlobalRelabelThreshold,
14 Level,
15 build_residual_network,
16 detect_unboundedness,
17)
19__all__ = ["preflow_push"]
22def preflow_push_impl(G, s, t, capacity, residual, global_relabel_freq, value_only):
23 """Implementation of the highest-label preflow-push algorithm."""
24 if s not in G:
25 raise nx.NetworkXError(f"node {str(s)} not in graph")
26 if t not in G:
27 raise nx.NetworkXError(f"node {str(t)} not in graph")
28 if s == t:
29 raise nx.NetworkXError("source and sink are the same node")
31 if global_relabel_freq is None:
32 global_relabel_freq = 0
33 if global_relabel_freq < 0:
34 raise nx.NetworkXError("global_relabel_freq must be nonnegative.")
36 if residual is None:
37 R = build_residual_network(G, capacity)
38 else:
39 R = residual
41 detect_unboundedness(R, s, t)
43 R_nodes = R.nodes
44 R_pred = R.pred
45 R_succ = R.succ
47 # Initialize/reset the residual network.
48 for u in R:
49 R_nodes[u]["excess"] = 0
50 for e in R_succ[u].values():
51 e["flow"] = 0
53 def reverse_bfs(src):
54 """Perform a reverse breadth-first search from src in the residual
55 network.
56 """
57 heights = {src: 0}
58 q = deque([(src, 0)])
59 while q:
60 u, height = q.popleft()
61 height += 1
62 for v, attr in R_pred[u].items():
63 if v not in heights and attr["flow"] < attr["capacity"]:
64 heights[v] = height
65 q.append((v, height))
66 return heights
68 # Initialize heights of the nodes.
69 heights = reverse_bfs(t)
71 if s not in heights:
72 # t is not reachable from s in the residual network. The maximum flow
73 # must be zero.
74 R.graph["flow_value"] = 0
75 return R
77 n = len(R)
78 # max_height represents the height of the highest level below level n with
79 # at least one active node.
80 max_height = max(heights[u] for u in heights if u != s)
81 heights[s] = n
83 grt = GlobalRelabelThreshold(n, R.size(), global_relabel_freq)
85 # Initialize heights and 'current edge' data structures of the nodes.
86 for u in R:
87 R_nodes[u]["height"] = heights[u] if u in heights else n + 1
88 R_nodes[u]["curr_edge"] = CurrentEdge(R_succ[u])
90 def push(u, v, flow):
91 """Push flow units of flow from u to v."""
92 R_succ[u][v]["flow"] += flow
93 R_succ[v][u]["flow"] -= flow
94 R_nodes[u]["excess"] -= flow
95 R_nodes[v]["excess"] += flow
97 # The maximum flow must be nonzero now. Initialize the preflow by
98 # saturating all edges emanating from s.
99 for u, attr in R_succ[s].items():
100 flow = attr["capacity"]
101 if flow > 0:
102 push(s, u, flow)
104 # Partition nodes into levels.
105 levels = [Level() for i in range(2 * n)]
106 for u in R:
107 if u != s and u != t:
108 level = levels[R_nodes[u]["height"]]
109 if R_nodes[u]["excess"] > 0:
110 level.active.add(u)
111 else:
112 level.inactive.add(u)
114 def activate(v):
115 """Move a node from the inactive set to the active set of its level."""
116 if v != s and v != t:
117 level = levels[R_nodes[v]["height"]]
118 if v in level.inactive:
119 level.inactive.remove(v)
120 level.active.add(v)
122 def relabel(u):
123 """Relabel a node to create an admissible edge."""
124 grt.add_work(len(R_succ[u]))
125 return (
126 min(
127 R_nodes[v]["height"]
128 for v, attr in R_succ[u].items()
129 if attr["flow"] < attr["capacity"]
130 )
131 + 1
132 )
134 def discharge(u, is_phase1):
135 """Discharge a node until it becomes inactive or, during phase 1 (see
136 below), its height reaches at least n. The node is known to have the
137 largest height among active nodes.
138 """
139 height = R_nodes[u]["height"]
140 curr_edge = R_nodes[u]["curr_edge"]
141 # next_height represents the next height to examine after discharging
142 # the current node. During phase 1, it is capped to below n.
143 next_height = height
144 levels[height].active.remove(u)
145 while True:
146 v, attr = curr_edge.get()
147 if height == R_nodes[v]["height"] + 1 and attr["flow"] < attr["capacity"]:
148 flow = min(R_nodes[u]["excess"], attr["capacity"] - attr["flow"])
149 push(u, v, flow)
150 activate(v)
151 if R_nodes[u]["excess"] == 0:
152 # The node has become inactive.
153 levels[height].inactive.add(u)
154 break
155 try:
156 curr_edge.move_to_next()
157 except StopIteration:
158 # We have run off the end of the adjacency list, and there can
159 # be no more admissible edges. Relabel the node to create one.
160 height = relabel(u)
161 if is_phase1 and height >= n - 1:
162 # Although the node is still active, with a height at least
163 # n - 1, it is now known to be on the s side of the minimum
164 # s-t cut. Stop processing it until phase 2.
165 levels[height].active.add(u)
166 break
167 # The first relabel operation after global relabeling may not
168 # increase the height of the node since the 'current edge' data
169 # structure is not rewound. Use height instead of (height - 1)
170 # in case other active nodes at the same level are missed.
171 next_height = height
172 R_nodes[u]["height"] = height
173 return next_height
175 def gap_heuristic(height):
176 """Apply the gap heuristic."""
177 # Move all nodes at levels (height + 1) to max_height to level n + 1.
178 for level in islice(levels, height + 1, max_height + 1):
179 for u in level.active:
180 R_nodes[u]["height"] = n + 1
181 for u in level.inactive:
182 R_nodes[u]["height"] = n + 1
183 levels[n + 1].active.update(level.active)
184 level.active.clear()
185 levels[n + 1].inactive.update(level.inactive)
186 level.inactive.clear()
188 def global_relabel(from_sink):
189 """Apply the global relabeling heuristic."""
190 src = t if from_sink else s
191 heights = reverse_bfs(src)
192 if not from_sink:
193 # s must be reachable from t. Remove t explicitly.
194 del heights[t]
195 max_height = max(heights.values())
196 if from_sink:
197 # Also mark nodes from which t is unreachable for relabeling. This
198 # serves the same purpose as the gap heuristic.
199 for u in R:
200 if u not in heights and R_nodes[u]["height"] < n:
201 heights[u] = n + 1
202 else:
203 # Shift the computed heights because the height of s is n.
204 for u in heights:
205 heights[u] += n
206 max_height += n
207 del heights[src]
208 for u, new_height in heights.items():
209 old_height = R_nodes[u]["height"]
210 if new_height != old_height:
211 if u in levels[old_height].active:
212 levels[old_height].active.remove(u)
213 levels[new_height].active.add(u)
214 else:
215 levels[old_height].inactive.remove(u)
216 levels[new_height].inactive.add(u)
217 R_nodes[u]["height"] = new_height
218 return max_height
220 # Phase 1: Find the maximum preflow by pushing as much flow as possible to
221 # t.
223 height = max_height
224 while height > 0:
225 # Discharge active nodes in the current level.
226 while True:
227 level = levels[height]
228 if not level.active:
229 # All active nodes in the current level have been discharged.
230 # Move to the next lower level.
231 height -= 1
232 break
233 # Record the old height and level for the gap heuristic.
234 old_height = height
235 old_level = level
236 u = arbitrary_element(level.active)
237 height = discharge(u, True)
238 if grt.is_reached():
239 # Global relabeling heuristic: Recompute the exact heights of
240 # all nodes.
241 height = global_relabel(True)
242 max_height = height
243 grt.clear_work()
244 elif not old_level.active and not old_level.inactive:
245 # Gap heuristic: If the level at old_height is empty (a 'gap'),
246 # a minimum cut has been identified. All nodes with heights
247 # above old_height can have their heights set to n + 1 and not
248 # be further processed before a maximum preflow is found.
249 gap_heuristic(old_height)
250 height = old_height - 1
251 max_height = height
252 else:
253 # Update the height of the highest level with at least one
254 # active node.
255 max_height = max(max_height, height)
257 # A maximum preflow has been found. The excess at t is the maximum flow
258 # value.
259 if value_only:
260 R.graph["flow_value"] = R_nodes[t]["excess"]
261 return R
263 # Phase 2: Convert the maximum preflow into a maximum flow by returning the
264 # excess to s.
266 # Relabel all nodes so that they have accurate heights.
267 height = global_relabel(False)
268 grt.clear_work()
270 # Continue to discharge the active nodes.
271 while height > n:
272 # Discharge active nodes in the current level.
273 while True:
274 level = levels[height]
275 if not level.active:
276 # All active nodes in the current level have been discharged.
277 # Move to the next lower level.
278 height -= 1
279 break
280 u = arbitrary_element(level.active)
281 height = discharge(u, False)
282 if grt.is_reached():
283 # Global relabeling heuristic.
284 height = global_relabel(False)
285 grt.clear_work()
287 R.graph["flow_value"] = R_nodes[t]["excess"]
288 return R
291@nx._dispatchable(edge_attrs={"capacity": float("inf")}, returns_graph=True)
292def preflow_push(
293 G, s, t, capacity="capacity", residual=None, global_relabel_freq=1, value_only=False
294):
295 r"""Find a maximum single-commodity flow using the highest-label
296 preflow-push algorithm.
298 This function returns the residual network resulting after computing
299 the maximum flow. See below for details about the conventions
300 NetworkX uses for defining residual networks.
302 This algorithm has a running time of $O(n^2 \sqrt{m})$ for $n$ nodes and
303 $m$ edges.
306 Parameters
307 ----------
308 G : NetworkX graph
309 Edges of the graph are expected to have an attribute called
310 'capacity'. If this attribute is not present, the edge is
311 considered to have infinite capacity.
313 s : node
314 Source node for the flow.
316 t : node
317 Sink node for the flow.
319 capacity : string
320 Edges of the graph G are expected to have an attribute capacity
321 that indicates how much flow the edge can support. If this
322 attribute is not present, the edge is considered to have
323 infinite capacity. Default value: 'capacity'.
325 residual : NetworkX graph
326 Residual network on which the algorithm is to be executed. If None, a
327 new residual network is created. Default value: None.
329 global_relabel_freq : integer, float
330 Relative frequency of applying the global relabeling heuristic to speed
331 up the algorithm. If it is None, the heuristic is disabled. Default
332 value: 1.
334 value_only : bool
335 If False, compute a maximum flow; otherwise, compute a maximum preflow
336 which is enough for computing the maximum flow value. Default value:
337 False.
339 Returns
340 -------
341 R : NetworkX DiGraph
342 Residual network after computing the maximum flow.
344 Raises
345 ------
346 NetworkXError
347 The algorithm does not support MultiGraph and MultiDiGraph. If
348 the input graph is an instance of one of these two classes, a
349 NetworkXError is raised.
351 NetworkXUnbounded
352 If the graph has a path of infinite capacity, the value of a
353 feasible flow on the graph is unbounded above and the function
354 raises a NetworkXUnbounded.
356 See also
357 --------
358 :meth:`maximum_flow`
359 :meth:`minimum_cut`
360 :meth:`edmonds_karp`
361 :meth:`shortest_augmenting_path`
363 Notes
364 -----
365 The residual network :samp:`R` from an input graph :samp:`G` has the
366 same nodes as :samp:`G`. :samp:`R` is a DiGraph that contains a pair
367 of edges :samp:`(u, v)` and :samp:`(v, u)` iff :samp:`(u, v)` is not a
368 self-loop, and at least one of :samp:`(u, v)` and :samp:`(v, u)` exists
369 in :samp:`G`. For each node :samp:`u` in :samp:`R`,
370 :samp:`R.nodes[u]['excess']` represents the difference between flow into
371 :samp:`u` and flow out of :samp:`u`.
373 For each edge :samp:`(u, v)` in :samp:`R`, :samp:`R[u][v]['capacity']`
374 is equal to the capacity of :samp:`(u, v)` in :samp:`G` if it exists
375 in :samp:`G` or zero otherwise. If the capacity is infinite,
376 :samp:`R[u][v]['capacity']` will have a high arbitrary finite value
377 that does not affect the solution of the problem. This value is stored in
378 :samp:`R.graph['inf']`. For each edge :samp:`(u, v)` in :samp:`R`,
379 :samp:`R[u][v]['flow']` represents the flow function of :samp:`(u, v)` and
380 satisfies :samp:`R[u][v]['flow'] == -R[v][u]['flow']`.
382 The flow value, defined as the total flow into :samp:`t`, the sink, is
383 stored in :samp:`R.graph['flow_value']`. Reachability to :samp:`t` using
384 only edges :samp:`(u, v)` such that
385 :samp:`R[u][v]['flow'] < R[u][v]['capacity']` induces a minimum
386 :samp:`s`-:samp:`t` cut.
388 Examples
389 --------
390 >>> from networkx.algorithms.flow import preflow_push
392 The functions that implement flow algorithms and output a residual
393 network, such as this one, are not imported to the base NetworkX
394 namespace, so you have to explicitly import them from the flow package.
396 >>> G = nx.DiGraph()
397 >>> G.add_edge("x", "a", capacity=3.0)
398 >>> G.add_edge("x", "b", capacity=1.0)
399 >>> G.add_edge("a", "c", capacity=3.0)
400 >>> G.add_edge("b", "c", capacity=5.0)
401 >>> G.add_edge("b", "d", capacity=4.0)
402 >>> G.add_edge("d", "e", capacity=2.0)
403 >>> G.add_edge("c", "y", capacity=2.0)
404 >>> G.add_edge("e", "y", capacity=3.0)
405 >>> R = preflow_push(G, "x", "y")
406 >>> flow_value = nx.maximum_flow_value(G, "x", "y")
407 >>> flow_value == R.graph["flow_value"]
408 True
409 >>> # preflow_push also stores the maximum flow value
410 >>> # in the excess attribute of the sink node t
411 >>> flow_value == R.nodes["y"]["excess"]
412 True
413 >>> # For some problems, you might only want to compute a
414 >>> # maximum preflow.
415 >>> R = preflow_push(G, "x", "y", value_only=True)
416 >>> flow_value == R.graph["flow_value"]
417 True
418 >>> flow_value == R.nodes["y"]["excess"]
419 True
421 """
422 R = preflow_push_impl(G, s, t, capacity, residual, global_relabel_freq, value_only)
423 R.graph["algorithm"] = "preflow_push"
424 nx._clear_cache(R)
425 return R