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(
292 edge_attrs={"capacity": float("inf")}, returns_graph=True, preserve_edge_attrs=True
293)
294def preflow_push(
295 G, s, t, capacity="capacity", residual=None, global_relabel_freq=1, value_only=False
296):
297 r"""Find a maximum single-commodity flow using the highest-label
298 preflow-push algorithm.
300 This function returns the residual network resulting after computing
301 the maximum flow. See below for details about the conventions
302 NetworkX uses for defining residual networks.
304 This algorithm has a running time of $O(n^2 \sqrt{m})$ for $n$ nodes and
305 $m$ edges.
308 Parameters
309 ----------
310 G : NetworkX graph
311 Edges of the graph are expected to have an attribute called
312 'capacity'. If this attribute is not present, the edge is
313 considered to have infinite capacity.
315 s : node
316 Source node for the flow.
318 t : node
319 Sink node for the flow.
321 capacity : string or function (default= 'capacity')
322 If this is a string, then edge capacity will be accessed via the
323 edge attribute with this key (that is, the capacity of the edge
324 joining `u` to `v` will be ``G.edges[u, v][capacity]``). If no
325 such edge attribute exists, the capacity of the edge is assumed to
326 be infinite.
328 If this is a function, the capacity of an edge is the value
329 returned by the function. The function must accept exactly three
330 positional arguments: the two endpoints of an edge and the
331 dictionary of edge attributes for that edge. The function must
332 return a number or None to indicate a hidden edge.
334 residual : NetworkX graph
335 Residual network on which the algorithm is to be executed. If None, a
336 new residual network is created. Default value: None.
338 global_relabel_freq : integer, float
339 Relative frequency of applying the global relabeling heuristic to speed
340 up the algorithm. If it is None, the heuristic is disabled. Default
341 value: 1.
343 value_only : bool
344 If False, compute a maximum flow; otherwise, compute a maximum preflow
345 which is enough for computing the maximum flow value. Default value:
346 False.
348 Returns
349 -------
350 R : NetworkX DiGraph
351 Residual network after computing the maximum flow.
353 Raises
354 ------
355 NetworkXError
356 The algorithm does not support MultiGraph and MultiDiGraph. If
357 the input graph is an instance of one of these two classes, a
358 NetworkXError is raised.
360 NetworkXUnbounded
361 If the graph has a path of infinite capacity, the value of a
362 feasible flow on the graph is unbounded above and the function
363 raises a NetworkXUnbounded.
365 See also
366 --------
367 :meth:`maximum_flow`
368 :meth:`minimum_cut`
369 :meth:`edmonds_karp`
370 :meth:`shortest_augmenting_path`
372 Notes
373 -----
374 The residual network :samp:`R` from an input graph :samp:`G` has the
375 same nodes as :samp:`G`. :samp:`R` is a DiGraph that contains a pair
376 of edges :samp:`(u, v)` and :samp:`(v, u)` iff :samp:`(u, v)` is not a
377 self-loop, and at least one of :samp:`(u, v)` and :samp:`(v, u)` exists
378 in :samp:`G`. For each node :samp:`u` in :samp:`R`,
379 :samp:`R.nodes[u]['excess']` represents the difference between flow into
380 :samp:`u` and flow out of :samp:`u`.
382 For each edge :samp:`(u, v)` in :samp:`R`, :samp:`R[u][v]['capacity']`
383 is equal to the capacity of :samp:`(u, v)` in :samp:`G` if it exists
384 in :samp:`G` or zero otherwise. If the capacity is infinite,
385 :samp:`R[u][v]['capacity']` will have a high arbitrary finite value
386 that does not affect the solution of the problem. This value is stored in
387 :samp:`R.graph['inf']`. For each edge :samp:`(u, v)` in :samp:`R`,
388 :samp:`R[u][v]['flow']` represents the flow function of :samp:`(u, v)` and
389 satisfies :samp:`R[u][v]['flow'] == -R[v][u]['flow']`.
391 The flow value, defined as the total flow into :samp:`t`, the sink, is
392 stored in :samp:`R.graph['flow_value']`. Reachability to :samp:`t` using
393 only edges :samp:`(u, v)` such that
394 :samp:`R[u][v]['flow'] < R[u][v]['capacity']` induces a minimum
395 :samp:`s`-:samp:`t` cut.
397 Examples
398 --------
399 >>> from networkx.algorithms.flow import preflow_push
401 The functions that implement flow algorithms and output a residual
402 network, such as this one, are not imported to the base NetworkX
403 namespace, so you have to explicitly import them from the flow package.
405 >>> G = nx.DiGraph()
406 >>> G.add_edge("x", "a", capacity=3.0)
407 >>> G.add_edge("x", "b", capacity=1.0)
408 >>> G.add_edge("a", "c", capacity=3.0)
409 >>> G.add_edge("b", "c", capacity=5.0)
410 >>> G.add_edge("b", "d", capacity=4.0)
411 >>> G.add_edge("d", "e", capacity=2.0)
412 >>> G.add_edge("c", "y", capacity=2.0)
413 >>> G.add_edge("e", "y", capacity=3.0)
414 >>> R = preflow_push(G, "x", "y")
415 >>> flow_value = nx.maximum_flow_value(G, "x", "y")
416 >>> flow_value == R.graph["flow_value"]
417 True
418 >>> # preflow_push also stores the maximum flow value
419 >>> # in the excess attribute of the sink node t
420 >>> flow_value == R.nodes["y"]["excess"]
421 True
422 >>> # For some problems, you might only want to compute a
423 >>> # maximum preflow.
424 >>> R = preflow_push(G, "x", "y", value_only=True)
425 >>> flow_value == R.graph["flow_value"]
426 True
427 >>> flow_value == R.nodes["y"]["excess"]
428 True
430 """
431 R = preflow_push_impl(G, s, t, capacity, residual, global_relabel_freq, value_only)
432 R.graph["algorithm"] = "preflow_push"
433 nx._clear_cache(R)
434 return R