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