ioftools / networkxMiCe / networkxmaster / networkx / algorithms / flow / boykovkolmogorov.py @ 5cef0f13
History  View  Annotate  Download (13.3 KB)
1 
# boykovkolmogorov.py  Boykov Kolmogorov algorithm for maximum flow problems.


2 
#

3 
# Copyright 20162019 NetworkX developers.

4 
#

5 
# This file is part of NetworkX.

6 
#

7 
# NetworkX is distributed under a BSD license; see LICENSE.txt for more

8 
# information.

9 
#

10 
# Author: Jordi Torrents <jordi.t21@gmail.com>

11 
"""

12 
BoykovKolmogorov algorithm for maximum flow problems.

13 
"""

14 
from collections import deque 
15 
from operator import itemgetter 
16  
17 
import networkx as nx 
18 
from networkx.algorithms.flow.utils import build_residual_network 
19  
20 
__all__ = ['boykov_kolmogorov']

21  
22  
23 
def boykov_kolmogorov(G, s, t, capacity='capacity', residual=None, 
24 
value_only=False, cutoff=None): 
25 
r"""Find a maximum singlecommodity flow using BoykovKolmogorov algorithm.

26 

27 
This function returns the residual network resulting after computing

28 
the maximum flow. See below for details about the conventions

29 
NetworkX uses for defining residual networks.

30 

31 
This algorithm has worse case complexity $O(n^2 m C)$ for $n$ nodes, $m$

32 
edges, and $C$ the cost of the minimum cut [1]_. This implementation

33 
uses the marking heuristic defined in [2]_ which improves its running

34 
time in many practical problems.

35 

36 
Parameters

37 


38 
G : NetworkX graph

39 
Edges of the graph are expected to have an attribute called

40 
'capacity'. If this attribute is not present, the edge is

41 
considered to have infinite capacity.

42 

43 
s : node

44 
Source node for the flow.

45 

46 
t : node

47 
Sink node for the flow.

48 

49 
capacity : string

50 
Edges of the graph G are expected to have an attribute capacity

51 
that indicates how much flow the edge can support. If this

52 
attribute is not present, the edge is considered to have

53 
infinite capacity. Default value: 'capacity'.

54 

55 
residual : NetworkX graph

56 
Residual network on which the algorithm is to be executed. If None, a

57 
new residual network is created. Default value: None.

58 

59 
value_only : bool

60 
If True compute only the value of the maximum flow. This parameter

61 
will be ignored by this algorithm because it is not applicable.

62 

63 
cutoff : integer, float

64 
If specified, the algorithm will terminate when the flow value reaches

65 
or exceeds the cutoff. In this case, it may be unable to immediately

66 
determine a minimum cut. Default value: None.

67 

68 
Returns

69 


70 
R : NetworkX DiGraph

71 
Residual network after computing the maximum flow.

72 

73 
Raises

74 


75 
NetworkXError

76 
The algorithm does not support MultiGraph and MultiDiGraph. If

77 
the input graph is an instance of one of these two classes, a

78 
NetworkXError is raised.

79 

80 
NetworkXUnbounded

81 
If the graph has a path of infinite capacity, the value of a

82 
feasible flow on the graph is unbounded above and the function

83 
raises a NetworkXUnbounded.

84 

85 
See also

86 


87 
:meth:`maximum_flow`

88 
:meth:`minimum_cut`

89 
:meth:`preflow_push`

90 
:meth:`shortest_augmenting_path`

91 

92 
Notes

93 


94 
The residual network :samp:`R` from an input graph :samp:`G` has the

95 
same nodes as :samp:`G`. :samp:`R` is a DiGraph that contains a pair

96 
of edges :samp:`(u, v)` and :samp:`(v, u)` iff :samp:`(u, v)` is not a

97 
selfloop, and at least one of :samp:`(u, v)` and :samp:`(v, u)` exists

98 
in :samp:`G`.

99 

100 
For each edge :samp:`(u, v)` in :samp:`R`, :samp:`R[u][v]['capacity']`

101 
is equal to the capacity of :samp:`(u, v)` in :samp:`G` if it exists

102 
in :samp:`G` or zero otherwise. If the capacity is infinite,

103 
:samp:`R[u][v]['capacity']` will have a high arbitrary finite value

104 
that does not affect the solution of the problem. This value is stored in

105 
:samp:`R.graph['inf']`. For each edge :samp:`(u, v)` in :samp:`R`,

106 
:samp:`R[u][v]['flow']` represents the flow function of :samp:`(u, v)` and

107 
satisfies :samp:`R[u][v]['flow'] == R[v][u]['flow']`.

108 

109 
The flow value, defined as the total flow into :samp:`t`, the sink, is

110 
stored in :samp:`R.graph['flow_value']`. If :samp:`cutoff` is not

111 
specified, reachability to :samp:`t` using only edges :samp:`(u, v)` such

112 
that :samp:`R[u][v]['flow'] < R[u][v]['capacity']` induces a minimum

113 
:samp:`s`:samp:`t` cut.

114 

115 
Examples

116 


117 
>>> import networkx as nx

118 
>>> from networkx.algorithms.flow import boykov_kolmogorov

119 

120 
The functions that implement flow algorithms and output a residual

121 
network, such as this one, are not imported to the base NetworkX

122 
namespace, so you have to explicitly import them from the flow package.

123 

124 
>>> G = nx.DiGraph()

125 
>>> G.add_edge('x','a', capacity=3.0)

126 
>>> G.add_edge('x','b', capacity=1.0)

127 
>>> G.add_edge('a','c', capacity=3.0)

128 
>>> G.add_edge('b','c', capacity=5.0)

129 
>>> G.add_edge('b','d', capacity=4.0)

130 
>>> G.add_edge('d','e', capacity=2.0)

131 
>>> G.add_edge('c','y', capacity=2.0)

132 
>>> G.add_edge('e','y', capacity=3.0)

133 
>>> R = boykov_kolmogorov(G, 'x', 'y')

134 
>>> flow_value = nx.maximum_flow_value(G, 'x', 'y')

135 
>>> flow_value

136 
3.0

137 
>>> flow_value == R.graph['flow_value']

138 
True

139 

140 
A nice feature of the BoykovKolmogorov algorithm is that a partition

141 
of the nodes that defines a minimum cut can be easily computed based

142 
on the search trees used during the algorithm. These trees are stored

143 
in the graph attribute `trees` of the residual network.

144 

145 
>>> source_tree, target_tree = R.graph['trees']

146 
>>> partition = (set(source_tree), set(G)  set(source_tree))

147 

148 
Or equivalently:

149 

150 
>>> partition = (set(G)  set(target_tree), set(target_tree))

151 

152 
References

153 


154 
.. [1] Boykov, Y., & Kolmogorov, V. (2004). An experimental comparison

155 
of mincut/maxflow algorithms for energy minimization in vision.

156 
Pattern Analysis and Machine Intelligence, IEEE Transactions on,

157 
26(9), 11241137.

158 
http://www.csd.uwo.ca/~yuri/Papers/pami04.pdf

159 

160 
.. [2] Vladimir Kolmogorov. Graphbased Algorithms for Multicamera

161 
Reconstruction Problem. PhD thesis, Cornell University, CS Department,

162 
2003. pp. 109114.

163 
https://pub.ist.ac.at/~vnk/papers/thesis.pdf

164 

165 
"""

166 
R = boykov_kolmogorov_impl(G, s, t, capacity, residual, cutoff) 
167 
R.graph['algorithm'] = 'boykov_kolmogorov' 
168 
return R

169  
170  
171 
def boykov_kolmogorov_impl(G, s, t, capacity, residual, cutoff): 
172 
if s not in G: 
173 
raise nx.NetworkXError('node %s not in graph' % str(s)) 
174 
if t not in G: 
175 
raise nx.NetworkXError('node %s not in graph' % str(t)) 
176 
if s == t:

177 
raise nx.NetworkXError('source and sink are the same node') 
178  
179 
if residual is None: 
180 
R = build_residual_network(G, capacity) 
181 
else:

182 
R = residual 
183  
184 
# Initialize/reset the residual network.

185 
# This is way too slow

186 
#nx.set_edge_attributes(R, 0, 'flow')

187 
for u in R: 
188 
for e in R[u].values(): 
189 
e['flow'] = 0 
190  
191 
# Use an arbitrary high value as infinite. It is computed

192 
# when building the residual network.

193 
INF = R.graph['inf']

194  
195 
if cutoff is None: 
196 
cutoff = INF 
197  
198 
R_succ = R.succ 
199 
R_pred = R.pred 
200  
201 
def grow(): 
202 
"""Bidirectional breadthfirst search for the growth stage.

203 

204 
Returns a connecting edge, that is and edge that connects

205 
a node from the source search tree with a node from the

206 
target search tree.

207 
The first node in the connecting edge is always from the

208 
source tree and the last node from the target tree.

209 
"""

210 
while active:

211 
u = active[0]

212 
if u in source_tree: 
213 
this_tree = source_tree 
214 
other_tree = target_tree 
215 
neighbors = R_succ 
216 
else:

217 
this_tree = target_tree 
218 
other_tree = source_tree 
219 
neighbors = R_pred 
220 
for v, attr in neighbors[u].items(): 
221 
if attr['capacity']  attr['flow'] > 0: 
222 
if v not in this_tree: 
223 
if v in other_tree: 
224 
return (u, v) if this_tree is source_tree else (v, u) 
225 
this_tree[v] = u 
226 
dist[v] = dist[u] + 1

227 
timestamp[v] = timestamp[u] 
228 
active.append(v) 
229 
elif v in this_tree and _is_closer(u, v): 
230 
this_tree[v] = u 
231 
dist[v] = dist[u] + 1

232 
timestamp[v] = timestamp[u] 
233 
_ = active.popleft() 
234 
return None, None 
235  
236 
def augment(u, v): 
237 
"""Augmentation stage.

238 

239 
Reconstruct path and determine its residual capacity.

240 
We start from a connecting edge, which links a node

241 
from the source tree to a node from the target tree.

242 
The connecting edge is the output of the grow function

243 
and the input of this function.

244 
"""

245 
attr = R_succ[u][v] 
246 
flow = min(INF, attr['capacity']  attr['flow']) 
247 
path = [u] 
248 
# Trace a path from u to s in source_tree.

249 
w = u 
250 
while w != s:

251 
n = w 
252 
w = source_tree[n] 
253 
attr = R_pred[n][w] 
254 
flow = min(flow, attr['capacity']  attr['flow']) 
255 
path.append(w) 
256 
path.reverse() 
257 
# Trace a path from v to t in target_tree.

258 
path.append(v) 
259 
w = v 
260 
while w != t:

261 
n = w 
262 
w = target_tree[n] 
263 
attr = R_succ[n][w] 
264 
flow = min(flow, attr['capacity']  attr['flow']) 
265 
path.append(w) 
266 
# Augment flow along the path and check for saturated edges.

267 
it = iter(path)

268 
u = next(it)

269 
these_orphans = [] 
270 
for v in it: 
271 
R_succ[u][v]['flow'] += flow

272 
R_succ[v][u]['flow'] = flow

273 
if R_succ[u][v]['flow'] == R_succ[u][v]['capacity']: 
274 
if v in source_tree: 
275 
source_tree[v] = None

276 
these_orphans.append(v) 
277 
if u in target_tree: 
278 
target_tree[u] = None

279 
these_orphans.append(u) 
280 
u = v 
281 
orphans.extend(sorted(these_orphans, key=dist.get))

282 
return flow

283  
284 
def adopt(): 
285 
"""Adoption stage.

286 

287 
Reconstruct search trees by adopting or discarding orphans.

288 
During augmentation stage some edges got saturated and thus

289 
the source and target search trees broke down to forests, with

290 
orphans as roots of some of its trees. We have to reconstruct

291 
the search trees rooted to source and target before we can grow

292 
them again.

293 
"""

294 
while orphans:

295 
u = orphans.popleft() 
296 
if u in source_tree: 
297 
tree = source_tree 
298 
neighbors = R_pred 
299 
else:

300 
tree = target_tree 
301 
neighbors = R_succ 
302 
nbrs = ((n, attr, dist[n]) for n, attr in neighbors[u].items() 
303 
if n in tree) 
304 
for v, attr, d in sorted(nbrs, key=itemgetter(2)): 
305 
if attr['capacity']  attr['flow'] > 0: 
306 
if _has_valid_root(v, tree):

307 
tree[u] = v 
308 
dist[u] = dist[v] + 1

309 
timestamp[u] = time 
310 
break

311 
else:

312 
nbrs = ((n, attr, dist[n]) for n, attr in neighbors[u].items() 
313 
if n in tree) 
314 
for v, attr, d in sorted(nbrs, key=itemgetter(2)): 
315 
if attr['capacity']  attr['flow'] > 0: 
316 
if v not in active: 
317 
active.append(v) 
318 
if tree[v] == u:

319 
tree[v] = None

320 
orphans.appendleft(v) 
321 
if u in active: 
322 
active.remove(u) 
323 
del tree[u]

324  
325 
def _has_valid_root(n, tree): 
326 
path = [] 
327 
v = n 
328 
while v is not None: 
329 
path.append(v) 
330 
if v == s or v == t: 
331 
base_dist = 0

332 
break

333 
elif timestamp[v] == time:

334 
base_dist = dist[v] 
335 
break

336 
v = tree[v] 
337 
else:

338 
return False 
339 
length = len(path)

340 
for i, u in enumerate(path, 1): 
341 
dist[u] = base_dist + length  i 
342 
timestamp[u] = time 
343 
return True 
344  
345 
def _is_closer(u, v): 
346 
return timestamp[v] <= timestamp[u] and dist[v] > dist[u] + 1 
347  
348 
source_tree = {s: None}

349 
target_tree = {t: None}

350 
active = deque([s, t]) 
351 
orphans = deque() 
352 
flow_value = 0

353 
# data structures for the marking heuristic

354 
time = 1

355 
timestamp = {s: time, t: time} 
356 
dist = {s: 0, t: 0} 
357 
while flow_value < cutoff:

358 
# Growth stage

359 
u, v = grow() 
360 
if u is None: 
361 
break

362 
time += 1

363 
# Augmentation stage

364 
flow_value += augment(u, v) 
365 
# Adoption stage

366 
adopt() 
367  
368 
if flow_value * 2 > INF: 
369 
raise nx.NetworkXUnbounded('Infinite capacity path, flow unbounded above.') 
370  
371 
# Add source and target tree in a graph attribute.

372 
# A partition that defines a minimum cut can be directly

373 
# computed from the search trees as explained in the docstrings.

374 
R.graph['trees'] = (source_tree, target_tree)

375 
# Add the standard flow_value graph attribute.

376 
R.graph['flow_value'] = flow_value

377 
return R
