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"""
2Boykov-Kolmogorov algorithm for maximum flow problems.
3"""
5from collections import deque
6from operator import itemgetter
8import networkx as nx
9from networkx.algorithms.flow.utils import build_residual_network
11__all__ = ["boykov_kolmogorov"]
14@nx._dispatchable(
15 edge_attrs={"capacity": float("inf")}, returns_graph=True, preserve_edge_attrs=True
16)
17def boykov_kolmogorov(
18 G, s, t, capacity="capacity", residual=None, value_only=False, cutoff=None
19):
20 r"""Find a maximum single-commodity flow using Boykov-Kolmogorov algorithm.
22 This function returns the residual network resulting after computing
23 the maximum flow. See below for details about the conventions
24 NetworkX uses for defining residual networks.
26 This algorithm has worse case complexity $O(n^2 m |C|)$ for $n$ nodes, $m$
27 edges, and $|C|$ the cost of the minimum cut [1]_. This implementation
28 uses the marking heuristic defined in [2]_ which improves its running
29 time in many practical problems.
31 Parameters
32 ----------
33 G : NetworkX graph
34 Edges of the graph are expected to have an attribute called
35 'capacity'. If this attribute is not present, the edge is
36 considered to have infinite capacity.
38 s : node
39 Source node for the flow.
41 t : node
42 Sink node for the flow.
44 capacity : string or function (default= 'capacity')
45 If this is a string, then edge capacity will be accessed via the
46 edge attribute with this key (that is, the capacity of the edge
47 joining `u` to `v` will be ``G.edges[u, v][capacity]``). If no
48 such edge attribute exists, the capacity of the edge is assumed to
49 be infinite.
51 If this is a function, the capacity of an edge is the value
52 returned by the function. The function must accept exactly three
53 positional arguments: the two endpoints of an edge and the
54 dictionary of edge attributes for that edge. The function must
55 return a number or None to indicate a hidden edge.
57 residual : NetworkX graph
58 Residual network on which the algorithm is to be executed. If None, a
59 new residual network is created. Default value: None.
61 value_only : bool
62 If True compute only the value of the maximum flow. This parameter
63 will be ignored by this algorithm because it is not applicable.
65 cutoff : integer, float
66 If specified, the algorithm will terminate when the flow value reaches
67 or exceeds the cutoff. In this case, it may be unable to immediately
68 determine a minimum cut. Default value: None.
70 Returns
71 -------
72 R : NetworkX DiGraph
73 Residual network after computing the maximum flow.
75 Raises
76 ------
77 NetworkXError
78 The algorithm does not support MultiGraph and MultiDiGraph. If
79 the input graph is an instance of one of these two classes, a
80 NetworkXError is raised.
82 NetworkXUnbounded
83 If the graph has a path of infinite capacity, the value of a
84 feasible flow on the graph is unbounded above and the function
85 raises a NetworkXUnbounded.
87 See also
88 --------
89 :meth:`maximum_flow`
90 :meth:`minimum_cut`
91 :meth:`preflow_push`
92 :meth:`shortest_augmenting_path`
94 Notes
95 -----
96 The residual network :samp:`R` from an input graph :samp:`G` has the
97 same nodes as :samp:`G`. :samp:`R` is a DiGraph that contains a pair
98 of edges :samp:`(u, v)` and :samp:`(v, u)` iff :samp:`(u, v)` is not a
99 self-loop, and at least one of :samp:`(u, v)` and :samp:`(v, u)` exists
100 in :samp:`G`.
102 For each edge :samp:`(u, v)` in :samp:`R`, :samp:`R[u][v]['capacity']`
103 is equal to the capacity of :samp:`(u, v)` in :samp:`G` if it exists
104 in :samp:`G` or zero otherwise. If the capacity is infinite,
105 :samp:`R[u][v]['capacity']` will have a high arbitrary finite value
106 that does not affect the solution of the problem. This value is stored in
107 :samp:`R.graph['inf']`. For each edge :samp:`(u, v)` in :samp:`R`,
108 :samp:`R[u][v]['flow']` represents the flow function of :samp:`(u, v)` and
109 satisfies :samp:`R[u][v]['flow'] == -R[v][u]['flow']`.
111 The flow value, defined as the total flow into :samp:`t`, the sink, is
112 stored in :samp:`R.graph['flow_value']`. If :samp:`cutoff` is not
113 specified, reachability to :samp:`t` using only edges :samp:`(u, v)` such
114 that :samp:`R[u][v]['flow'] < R[u][v]['capacity']` induces a minimum
115 :samp:`s`-:samp:`t` cut.
117 Examples
118 --------
119 >>> from networkx.algorithms.flow import boykov_kolmogorov
121 The functions that implement flow algorithms and output a residual
122 network, such as this one, are not imported to the base NetworkX
123 namespace, so you have to explicitly import them from the flow package.
125 >>> G = nx.DiGraph()
126 >>> G.add_edge("x", "a", capacity=3.0)
127 >>> G.add_edge("x", "b", capacity=1.0)
128 >>> G.add_edge("a", "c", capacity=3.0)
129 >>> G.add_edge("b", "c", capacity=5.0)
130 >>> G.add_edge("b", "d", capacity=4.0)
131 >>> G.add_edge("d", "e", capacity=2.0)
132 >>> G.add_edge("c", "y", capacity=2.0)
133 >>> G.add_edge("e", "y", capacity=3.0)
134 >>> R = boykov_kolmogorov(G, "x", "y")
135 >>> flow_value = nx.maximum_flow_value(G, "x", "y")
136 >>> flow_value
137 3.0
138 >>> flow_value == R.graph["flow_value"]
139 True
141 A nice feature of the Boykov-Kolmogorov algorithm is that a partition
142 of the nodes that defines a minimum cut can be easily computed based
143 on the search trees used during the algorithm. These trees are stored
144 in the graph attribute `trees` of the residual network.
146 >>> source_tree, target_tree = R.graph["trees"]
147 >>> partition = (set(source_tree), set(G) - set(source_tree))
149 Or equivalently:
151 >>> partition = (set(G) - set(target_tree), set(target_tree))
153 References
154 ----------
155 .. [1] Boykov, Y., & Kolmogorov, V. (2004). An experimental comparison
156 of min-cut/max-flow algorithms for energy minimization in vision.
157 Pattern Analysis and Machine Intelligence, IEEE Transactions on,
158 26(9), 1124-1137.
159 https://doi.org/10.1109/TPAMI.2004.60
161 .. [2] Vladimir Kolmogorov. Graph-based Algorithms for Multi-camera
162 Reconstruction Problem. PhD thesis, Cornell University, CS Department,
163 2003. pp. 109-114.
164 https://web.archive.org/web/20170809091249/https://pub.ist.ac.at/~vnk/papers/thesis.pdf
166 """
167 R = boykov_kolmogorov_impl(G, s, t, capacity, residual, cutoff)
168 R.graph["algorithm"] = "boykov_kolmogorov"
169 nx._clear_cache(R)
170 return R
173def boykov_kolmogorov_impl(G, s, t, capacity, residual, cutoff):
174 if s not in G:
175 raise nx.NetworkXError(f"node {str(s)} not in graph")
176 if t not in G:
177 raise nx.NetworkXError(f"node {str(t)} not in graph")
178 if s == t:
179 raise nx.NetworkXError("source and sink are the same node")
181 if residual is None:
182 R = build_residual_network(G, capacity)
183 else:
184 R = residual
186 # Initialize/reset the residual network.
187 # This is way too slow
188 # nx.set_edge_attributes(R, 0, 'flow')
189 for u in R:
190 for e in R[u].values():
191 e["flow"] = 0
193 # Use an arbitrary high value as infinite. It is computed
194 # when building the residual network.
195 INF = R.graph["inf"]
197 if cutoff is None:
198 cutoff = INF
200 R_succ = R.succ
201 R_pred = R.pred
203 def grow():
204 """Bidirectional breadth-first search for the growth stage.
206 Returns a connecting edge, that is and edge that connects
207 a node from the source search tree with a node from the
208 target search tree.
209 The first node in the connecting edge is always from the
210 source tree and the last node from the target tree.
211 """
212 while active:
213 u = active[0]
214 if u in source_tree:
215 this_tree = source_tree
216 other_tree = target_tree
217 neighbors = R_succ
218 else:
219 this_tree = target_tree
220 other_tree = source_tree
221 neighbors = R_pred
222 for v, attr in neighbors[u].items():
223 if attr["capacity"] - attr["flow"] > 0:
224 if v not in this_tree:
225 if v in other_tree:
226 return (u, v) if this_tree is source_tree else (v, u)
227 this_tree[v] = u
228 dist[v] = dist[u] + 1
229 timestamp[v] = timestamp[u]
230 active.append(v)
231 elif v in this_tree and _is_closer(u, v):
232 this_tree[v] = u
233 dist[v] = dist[u] + 1
234 timestamp[v] = timestamp[u]
235 _ = active.popleft()
236 return None, None
238 def augment(u, v):
239 """Augmentation stage.
241 Reconstruct path and determine its residual capacity.
242 We start from a connecting edge, which links a node
243 from the source tree to a node from the target tree.
244 The connecting edge is the output of the grow function
245 and the input of this function.
246 """
247 attr = R_succ[u][v]
248 flow = min(INF, attr["capacity"] - attr["flow"])
249 path = [u]
250 # Trace a path from u to s in source_tree.
251 w = u
252 while w != s:
253 n = w
254 w = source_tree[n]
255 attr = R_pred[n][w]
256 flow = min(flow, attr["capacity"] - attr["flow"])
257 path.append(w)
258 path.reverse()
259 # Trace a path from v to t in target_tree.
260 path.append(v)
261 w = v
262 while w != t:
263 n = w
264 w = target_tree[n]
265 attr = R_succ[n][w]
266 flow = min(flow, attr["capacity"] - attr["flow"])
267 path.append(w)
268 # Augment flow along the path and check for saturated edges.
269 it = iter(path)
270 u = next(it)
271 these_orphans = []
272 for v in it:
273 R_succ[u][v]["flow"] += flow
274 R_succ[v][u]["flow"] -= flow
275 if R_succ[u][v]["flow"] == R_succ[u][v]["capacity"]:
276 if v in source_tree:
277 source_tree[v] = None
278 these_orphans.append(v)
279 if u in target_tree:
280 target_tree[u] = None
281 these_orphans.append(u)
282 u = v
283 orphans.extend(sorted(these_orphans, key=dist.get))
284 return flow
286 def adopt():
287 """Adoption stage.
289 Reconstruct search trees by adopting or discarding orphans.
290 During augmentation stage some edges got saturated and thus
291 the source and target search trees broke down to forests, with
292 orphans as roots of some of its trees. We have to reconstruct
293 the search trees rooted to source and target before we can grow
294 them again.
295 """
296 while orphans:
297 u = orphans.popleft()
298 if u in source_tree:
299 tree = source_tree
300 neighbors = R_pred
301 else:
302 tree = target_tree
303 neighbors = R_succ
304 nbrs = ((n, attr, dist[n]) for n, attr in neighbors[u].items() if n in tree)
305 for v, attr, d in sorted(nbrs, key=itemgetter(2)):
306 if attr["capacity"] - attr["flow"] > 0:
307 if _has_valid_root(v, tree):
308 tree[u] = v
309 dist[u] = dist[v] + 1
310 timestamp[u] = time
311 break
312 else:
313 nbrs = (
314 (n, attr, dist[n]) for n, attr in neighbors[u].items() if n in tree
315 )
316 for v, attr, d in sorted(nbrs, key=itemgetter(2)):
317 if attr["capacity"] - attr["flow"] > 0:
318 if v not in active:
319 active.append(v)
320 if tree[v] == u:
321 tree[v] = None
322 orphans.appendleft(v)
323 if u in active:
324 active.remove(u)
325 del tree[u]
327 def _has_valid_root(n, tree):
328 path = []
329 v = n
330 while v is not None:
331 path.append(v)
332 if v in (s, t):
333 base_dist = 0
334 break
335 elif timestamp[v] == time:
336 base_dist = dist[v]
337 break
338 v = tree[v]
339 else:
340 return False
341 length = len(path)
342 for i, u in enumerate(path, 1):
343 dist[u] = base_dist + length - i
344 timestamp[u] = time
345 return True
347 def _is_closer(u, v):
348 return timestamp[v] <= timestamp[u] and dist[v] > dist[u] + 1
350 source_tree = {s: None}
351 target_tree = {t: None}
352 active = deque([s, t])
353 orphans = deque()
354 flow_value = 0
355 # data structures for the marking heuristic
356 time = 1
357 timestamp = {s: time, t: time}
358 dist = {s: 0, t: 0}
359 while flow_value < cutoff:
360 # Growth stage
361 u, v = grow()
362 if u is None:
363 break
364 time += 1
365 # Augmentation stage
366 flow_value += augment(u, v)
367 # Adoption stage
368 adopt()
370 if flow_value * 2 > INF:
371 raise nx.NetworkXUnbounded("Infinite capacity path, flow unbounded above.")
373 # Add source and target tree in a graph attribute.
374 # A partition that defines a minimum cut can be directly
375 # computed from the search trees as explained in the docstrings.
376 R.graph["trees"] = (source_tree, target_tree)
377 # Add the standard flow_value graph attribute.
378 R.graph["flow_value"] = flow_value
379 return R