Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/networkx/algorithms/flow/preflowpush.py: 5%
173 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-10-20 07:00 +0000
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._dispatch(
292 graphs={"G": 0, "residual?": 4},
293 edge_attrs={"capacity": float("inf")},
294 preserve_edge_attrs={"residual": {"capacity": float("inf")}},
295 preserve_graph_attrs={"residual"},
296)
297def preflow_push(
298 G, s, t, capacity="capacity", residual=None, global_relabel_freq=1, value_only=False
299):
300 r"""Find a maximum single-commodity flow using the highest-label
301 preflow-push algorithm.
303 This function returns the residual network resulting after computing
304 the maximum flow. See below for details about the conventions
305 NetworkX uses for defining residual networks.
307 This algorithm has a running time of $O(n^2 \sqrt{m})$ for $n$ nodes and
308 $m$ edges.
311 Parameters
312 ----------
313 G : NetworkX graph
314 Edges of the graph are expected to have an attribute called
315 'capacity'. If this attribute is not present, the edge is
316 considered to have infinite capacity.
318 s : node
319 Source node for the flow.
321 t : node
322 Sink node for the flow.
324 capacity : string
325 Edges of the graph G are expected to have an attribute capacity
326 that indicates how much flow the edge can support. If this
327 attribute is not present, the edge is considered to have
328 infinite capacity. Default value: 'capacity'.
330 residual : NetworkX graph
331 Residual network on which the algorithm is to be executed. If None, a
332 new residual network is created. Default value: None.
334 global_relabel_freq : integer, float
335 Relative frequency of applying the global relabeling heuristic to speed
336 up the algorithm. If it is None, the heuristic is disabled. Default
337 value: 1.
339 value_only : bool
340 If False, compute a maximum flow; otherwise, compute a maximum preflow
341 which is enough for computing the maximum flow value. Default value:
342 False.
344 Returns
345 -------
346 R : NetworkX DiGraph
347 Residual network after computing the maximum flow.
349 Raises
350 ------
351 NetworkXError
352 The algorithm does not support MultiGraph and MultiDiGraph. If
353 the input graph is an instance of one of these two classes, a
354 NetworkXError is raised.
356 NetworkXUnbounded
357 If the graph has a path of infinite capacity, the value of a
358 feasible flow on the graph is unbounded above and the function
359 raises a NetworkXUnbounded.
361 See also
362 --------
363 :meth:`maximum_flow`
364 :meth:`minimum_cut`
365 :meth:`edmonds_karp`
366 :meth:`shortest_augmenting_path`
368 Notes
369 -----
370 The residual network :samp:`R` from an input graph :samp:`G` has the
371 same nodes as :samp:`G`. :samp:`R` is a DiGraph that contains a pair
372 of edges :samp:`(u, v)` and :samp:`(v, u)` iff :samp:`(u, v)` is not a
373 self-loop, and at least one of :samp:`(u, v)` and :samp:`(v, u)` exists
374 in :samp:`G`. For each node :samp:`u` in :samp:`R`,
375 :samp:`R.nodes[u]['excess']` represents the difference between flow into
376 :samp:`u` and flow out of :samp:`u`.
378 For each edge :samp:`(u, v)` in :samp:`R`, :samp:`R[u][v]['capacity']`
379 is equal to the capacity of :samp:`(u, v)` in :samp:`G` if it exists
380 in :samp:`G` or zero otherwise. If the capacity is infinite,
381 :samp:`R[u][v]['capacity']` will have a high arbitrary finite value
382 that does not affect the solution of the problem. This value is stored in
383 :samp:`R.graph['inf']`. For each edge :samp:`(u, v)` in :samp:`R`,
384 :samp:`R[u][v]['flow']` represents the flow function of :samp:`(u, v)` and
385 satisfies :samp:`R[u][v]['flow'] == -R[v][u]['flow']`.
387 The flow value, defined as the total flow into :samp:`t`, the sink, is
388 stored in :samp:`R.graph['flow_value']`. Reachability to :samp:`t` using
389 only edges :samp:`(u, v)` such that
390 :samp:`R[u][v]['flow'] < R[u][v]['capacity']` induces a minimum
391 :samp:`s`-:samp:`t` cut.
393 Examples
394 --------
395 >>> from networkx.algorithms.flow import preflow_push
397 The functions that implement flow algorithms and output a residual
398 network, such as this one, are not imported to the base NetworkX
399 namespace, so you have to explicitly import them from the flow package.
401 >>> G = nx.DiGraph()
402 >>> G.add_edge("x", "a", capacity=3.0)
403 >>> G.add_edge("x", "b", capacity=1.0)
404 >>> G.add_edge("a", "c", capacity=3.0)
405 >>> G.add_edge("b", "c", capacity=5.0)
406 >>> G.add_edge("b", "d", capacity=4.0)
407 >>> G.add_edge("d", "e", capacity=2.0)
408 >>> G.add_edge("c", "y", capacity=2.0)
409 >>> G.add_edge("e", "y", capacity=3.0)
410 >>> R = preflow_push(G, "x", "y")
411 >>> flow_value = nx.maximum_flow_value(G, "x", "y")
412 >>> flow_value == R.graph["flow_value"]
413 True
414 >>> # preflow_push also stores the maximum flow value
415 >>> # in the excess attribute of the sink node t
416 >>> flow_value == R.nodes["y"]["excess"]
417 True
418 >>> # For some problems, you might only want to compute a
419 >>> # maximum preflow.
420 >>> R = preflow_push(G, "x", "y", value_only=True)
421 >>> flow_value == R.graph["flow_value"]
422 True
423 >>> flow_value == R.nodes["y"]["excess"]
424 True
426 """
427 R = preflow_push_impl(G, s, t, capacity, residual, global_relabel_freq, value_only)
428 R.graph["algorithm"] = "preflow_push"
429 return R