### comparison env/lib/python3.9/site-packages/networkx/algorithms/flow/preflowpush.py @ 0:4f3585e2f14bdraftdefaulttip

author shellac Mon, 22 Mar 2021 18:12:50 +0000
comparison
equal inserted replaced
-1:000000000000 0:4f3585e2f14b
1 """
2 Highest-label preflow-push algorithm for maximum flow problems.
3 """
4
5 from collections import deque
6 from itertools import islice
7 import networkx as nx
8 from ...utils import arbitrary_element
9 from .utils import build_residual_network
10 from .utils import CurrentEdge
11 from .utils import detect_unboundedness
12 from .utils import GlobalRelabelThreshold
13 from .utils import Level
14
15 __all__ = ["preflow_push"]
16
17
18 def preflow_push_impl(G, s, t, capacity, residual, global_relabel_freq, value_only):
19 """Implementation of the highest-label preflow-push algorithm.
20 """
21 if s not in G:
22 raise nx.NetworkXError(f"node {str(s)} not in graph")
23 if t not in G:
24 raise nx.NetworkXError(f"node {str(t)} not in graph")
25 if s == t:
26 raise nx.NetworkXError("source and sink are the same node")
27
28 if global_relabel_freq is None:
29 global_relabel_freq = 0
30 if global_relabel_freq < 0:
31 raise nx.NetworkXError("global_relabel_freq must be nonnegative.")
32
33 if residual is None:
34 R = build_residual_network(G, capacity)
35 else:
36 R = residual
37
38 detect_unboundedness(R, s, t)
39
40 R_nodes = R.nodes
41 R_pred = R.pred
42 R_succ = R.succ
43
44 # Initialize/reset the residual network.
45 for u in R:
46 R_nodes[u]["excess"] = 0
47 for e in R_succ[u].values():
48 e["flow"] = 0
49
50 def reverse_bfs(src):
51 """Perform a reverse breadth-first search from src in the residual
52 network.
53 """
54 heights = {src: 0}
55 q = deque([(src, 0)])
56 while q:
57 u, height = q.popleft()
58 height += 1
59 for v, attr in R_pred[u].items():
60 if v not in heights and attr["flow"] < attr["capacity"]:
61 heights[v] = height
62 q.append((v, height))
63 return heights
64
65 # Initialize heights of the nodes.
66 heights = reverse_bfs(t)
67
68 if s not in heights:
69 # t is not reachable from s in the residual network. The maximum flow
70 # must be zero.
71 R.graph["flow_value"] = 0
72 return R
73
74 n = len(R)
75 # max_height represents the height of the highest level below level n with
76 # at least one active node.
77 max_height = max(heights[u] for u in heights if u != s)
78 heights[s] = n
79
80 grt = GlobalRelabelThreshold(n, R.size(), global_relabel_freq)
81
82 # Initialize heights and 'current edge' data structures of the nodes.
83 for u in R:
84 R_nodes[u]["height"] = heights[u] if u in heights else n + 1
85 R_nodes[u]["curr_edge"] = CurrentEdge(R_succ[u])
86
87 def push(u, v, flow):
88 """Push flow units of flow from u to v.
89 """
90 R_succ[u][v]["flow"] += flow
91 R_succ[v][u]["flow"] -= flow
92 R_nodes[u]["excess"] -= flow
93 R_nodes[v]["excess"] += flow
94
95 # The maximum flow must be nonzero now. Initialize the preflow by
96 # saturating all edges emanating from s.
97 for u, attr in R_succ[s].items():
98 flow = attr["capacity"]
99 if flow > 0:
100 push(s, u, flow)
101
102 # Partition nodes into levels.
103 levels = [Level() for i in range(2 * n)]
104 for u in R:
105 if u != s and u != t:
106 level = levels[R_nodes[u]["height"]]
107 if R_nodes[u]["excess"] > 0:
109 else:
111
112 def activate(v):
113 """Move a node from the inactive set to the active set of its level.
114 """
115 if v != s and v != t:
116 level = levels[R_nodes[v]["height"]]
117 if v in level.inactive:
118 level.inactive.remove(v)
120
121 def relabel(u):
122 """Relabel a node to create an admissible edge.
123 """
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 )
133
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.
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.
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
174
175 def gap_heuristic(height):
176 """Apply the gap heuristic.
177 """
178 # Move all nodes at levels (height + 1) to max_height to level n + 1.
179 for level in islice(levels, height + 1, max_height + 1):
180 for u in level.active:
181 R_nodes[u]["height"] = n + 1
182 for u in level.inactive:
183 R_nodes[u]["height"] = n + 1
184 levels[n + 1].active.update(level.active)
185 level.active.clear()
186 levels[n + 1].inactive.update(level.inactive)
187 level.inactive.clear()
188
189 def global_relabel(from_sink):
190 """Apply the global relabeling heuristic.
191 """
192 src = t if from_sink else s
193 heights = reverse_bfs(src)
194 if not from_sink:
195 # s must be reachable from t. Remove t explicitly.
196 del heights[t]
197 max_height = max(heights.values())
198 if from_sink:
199 # Also mark nodes from which t is unreachable for relabeling. This
200 # serves the same purpose as the gap heuristic.
201 for u in R:
202 if u not in heights and R_nodes[u]["height"] < n:
203 heights[u] = n + 1
204 else:
205 # Shift the computed heights because the height of s is n.
206 for u in heights:
207 heights[u] += n
208 max_height += n
209 del heights[src]
210 for u, new_height in heights.items():
211 old_height = R_nodes[u]["height"]
212 if new_height != old_height:
213 if u in levels[old_height].active:
214 levels[old_height].active.remove(u)
216 else:
217 levels[old_height].inactive.remove(u)
219 R_nodes[u]["height"] = new_height
220 return max_height
221
222 # Phase 1: Find the maximum preflow by pushing as much flow as possible to
223 # t.
224
225 height = max_height
226 while height > 0:
227 # Discharge active nodes in the current level.
228 while True:
229 level = levels[height]
230 if not level.active:
231 # All active nodes in the current level have been discharged.
232 # Move to the next lower level.
233 height -= 1
234 break
235 # Record the old height and level for the gap heuristic.
236 old_height = height
237 old_level = level
238 u = arbitrary_element(level.active)
239 height = discharge(u, True)
240 if grt.is_reached():
241 # Global relabeling heuristic: Recompute the exact heights of
242 # all nodes.
243 height = global_relabel(True)
244 max_height = height
245 grt.clear_work()
246 elif not old_level.active and not old_level.inactive:
247 # Gap heuristic: If the level at old_height is empty (a 'gap'),
248 # a minimum cut has been identified. All nodes with heights
249 # above old_height can have their heights set to n + 1 and not
250 # be further processed before a maximum preflow is found.
251 gap_heuristic(old_height)
252 height = old_height - 1
253 max_height = height
254 else:
255 # Update the height of the highest level with at least one
256 # active node.
257 max_height = max(max_height, height)
258
259 # A maximum preflow has been found. The excess at t is the maximum flow
260 # value.
261 if value_only:
262 R.graph["flow_value"] = R_nodes[t]["excess"]
263 return R
264
265 # Phase 2: Convert the maximum preflow into a maximum flow by returning the
266 # excess to s.
267
268 # Relabel all nodes so that they have accurate heights.
269 height = global_relabel(False)
270 grt.clear_work()
271
272 # Continue to discharge the active nodes.
273 while height > n:
274 # Discharge active nodes in the current level.
275 while True:
276 level = levels[height]
277 if not level.active:
278 # All active nodes in the current level have been discharged.
279 # Move to the next lower level.
280 height -= 1
281 break
282 u = arbitrary_element(level.active)
283 height = discharge(u, False)
284 if grt.is_reached():
285 # Global relabeling heuristic.
286 height = global_relabel(False)
287 grt.clear_work()
288
289 R.graph["flow_value"] = R_nodes[t]["excess"]
290 return R
291
292
293 def preflow_push(
294 G, s, t, capacity="capacity", residual=None, global_relabel_freq=1, value_only=False
295 ):
296 r"""Find a maximum single-commodity flow using the highest-label
297 preflow-push algorithm.
298
299 This function returns the residual network resulting after computing
300 the maximum flow. See below for details about the conventions
301 NetworkX uses for defining residual networks.
302
303 This algorithm has a running time of $O(n^2 \sqrt{m})$ for $n$ nodes and
304 $m$ edges.
305
306
307 Parameters
308 ----------
309 G : NetworkX graph
310 Edges of the graph are expected to have an attribute called
311 'capacity'. If this attribute is not present, the edge is
312 considered to have infinite capacity.
313
314 s : node
315 Source node for the flow.
316
317 t : node
318 Sink node for the flow.
319
320 capacity : string
321 Edges of the graph G are expected to have an attribute capacity
322 that indicates how much flow the edge can support. If this
323 attribute is not present, the edge is considered to have
324 infinite capacity. Default value: 'capacity'.
325
326 residual : NetworkX graph
327 Residual network on which the algorithm is to be executed. If None, a
328 new residual network is created. Default value: None.
329
330 global_relabel_freq : integer, float
331 Relative frequency of applying the global relabeling heuristic to speed
332 up the algorithm. If it is None, the heuristic is disabled. Default
333 value: 1.
334
335 value_only : bool
336 If False, compute a maximum flow; otherwise, compute a maximum preflow
337 which is enough for computing the maximum flow value. Default value:
338 False.
339
340 Returns
341 -------
342 R : NetworkX DiGraph
343 Residual network after computing the maximum flow.
344
345 Raises
346 ------
347 NetworkXError
348 The algorithm does not support MultiGraph and MultiDiGraph. If
349 the input graph is an instance of one of these two classes, a
350 NetworkXError is raised.
351
352 NetworkXUnbounded
353 If the graph has a path of infinite capacity, the value of a
354 feasible flow on the graph is unbounded above and the function
355 raises a NetworkXUnbounded.
356
358 --------
359 :meth:maximum_flow
360 :meth:minimum_cut
361 :meth:edmonds_karp
362 :meth:shortest_augmenting_path
363
364 Notes
365 -----
366 The residual network :samp:R from an input graph :samp:G has the
367 same nodes as :samp:G. :samp:R is a DiGraph that contains a pair
368 of edges :samp:(u, v) and :samp:(v, u) iff :samp:(u, v) is not a
369 self-loop, and at least one of :samp:(u, v) and :samp:(v, u) exists
370 in :samp:G. For each node :samp:u in :samp:R,
371 :samp:R.nodes[u]['excess'] represents the difference between flow into
372 :samp:u and flow out of :samp:u.
373
374 For each edge :samp:(u, v) in :samp:R, :samp:R[u][v]['capacity']
375 is equal to the capacity of :samp:(u, v) in :samp:G if it exists
376 in :samp:G or zero otherwise. If the capacity is infinite,
377 :samp:R[u][v]['capacity'] will have a high arbitrary finite value
378 that does not affect the solution of the problem. This value is stored in
379 :samp:R.graph['inf']. For each edge :samp:(u, v) in :samp:R,
380 :samp:R[u][v]['flow'] represents the flow function of :samp:(u, v) and
381 satisfies :samp:R[u][v]['flow'] == -R[v][u]['flow'].
382
383 The flow value, defined as the total flow into :samp:t, the sink, is
384 stored in :samp:R.graph['flow_value']. Reachability to :samp:t using
385 only edges :samp:(u, v) such that
386 :samp:R[u][v]['flow'] < R[u][v]['capacity'] induces a minimum
387 :samp:s-:samp:t cut.
388
389 Examples
390 --------
391 >>> from networkx.algorithms.flow import preflow_push
392
393 The functions that implement flow algorithms and output a residual
394 network, such as this one, are not imported to the base NetworkX
395 namespace, so you have to explicitly import them from the flow package.
396
397 >>> G = nx.DiGraph()
406 >>> R = preflow_push(G, "x", "y")
407 >>> flow_value = nx.maximum_flow_value(G, "x", "y")
408 >>> flow_value == R.graph["flow_value"]
409 True
410 >>> # preflow_push also stores the maximum flow value
411 >>> # in the excess attribute of the sink node t
412 >>> flow_value == R.nodes["y"]["excess"]
413 True
414 >>> # For some problems, you might only want to compute a
415 >>> # maximum preflow.
416 >>> R = preflow_push(G, "x", "y", value_only=True)
417 >>> flow_value == R.graph["flow_value"]
418 True
419 >>> flow_value == R.nodes["y"]["excess"]
420 True
421
422 """
423 R = preflow_push_impl(G, s, t, capacity, residual, global_relabel_freq, value_only)
424 R.graph["algorithm"] = "preflow_push"
425 return R