## iof-tools / networkxMiCe / networkx-master / networkx / algorithms / flow / networksimplex.py @ 5cef0f13

History | View | Annotate | Download (20.7 KB)

1 | 5cef0f13 | tiamilani | ```
# -*- coding: utf-8 -*-
``` |
---|---|---|---|

2 | ```
"""
``` |
||

3 | ```
Minimum cost flow algorithms on directed connected graphs.
``` |
||

4 | ```
"""
``` |
||

5 | |||

6 | ```
__author__ = """Loïc Séguin-C. <loicseguin@gmail.com>"""
``` |
||

7 | ```
# Copyright (C) 2010 Loïc Séguin-C. <loicseguin@gmail.com>
``` |
||

8 | ```
# All rights reserved.
``` |
||

9 | ```
# BSD license.
``` |
||

10 | |||

11 | ```
__all__ = ['network_simplex']
``` |
||

12 | |||

13 | from itertools import chain, islice, repeat |
||

14 | from math import ceil, sqrt |
||

15 | import networkx as nx |
||

16 | from networkx.utils import not_implemented_for |
||

17 | |||

18 | ```
try:
``` |
||

19 | from itertools import izip as zip |
||

20 | except ImportError: |
||

21 | ```
pass
``` |
||

22 | ```
try:
``` |
||

23 | ```
range = xrange
``` |
||

24 | except NameError: |
||

25 | ```
pass
``` |
||

26 | |||

27 | |||

28 | @not_implemented_for('undirected') |
||

29 | def network_simplex(G, demand='demand', capacity='capacity', weight='weight'): |
||

30 | ```
r"""Find a minimum cost flow satisfying all demands in digraph G.
``` |
||

31 | ```
``` |
||

32 | ```
This is a primal network simplex algorithm that uses the leaving
``` |
||

33 | ```
arc rule to prevent cycling.
``` |
||

34 | ```
``` |
||

35 | ```
G is a digraph with edge costs and capacities and in which nodes
``` |
||

36 | ```
have demand, i.e., they want to send or receive some amount of
``` |
||

37 | ```
flow. A negative demand means that the node wants to send flow, a
``` |
||

38 | ```
positive demand means that the node want to receive flow. A flow on
``` |
||

39 | ```
the digraph G satisfies all demand if the net flow into each node
``` |
||

40 | ```
is equal to the demand of that node.
``` |
||

41 | ```
``` |
||

42 | ```
Parameters
``` |
||

43 | ```
----------
``` |
||

44 | ```
G : NetworkX graph
``` |
||

45 | ```
DiGraph on which a minimum cost flow satisfying all demands is
``` |
||

46 | ```
to be found.
``` |
||

47 | ```
``` |
||

48 | ```
demand : string
``` |
||

49 | ```
Nodes of the graph G are expected to have an attribute demand
``` |
||

50 | ```
that indicates how much flow a node wants to send (negative
``` |
||

51 | ```
demand) or receive (positive demand). Note that the sum of the
``` |
||

52 | ```
demands should be 0 otherwise the problem in not feasible. If
``` |
||

53 | ```
this attribute is not present, a node is considered to have 0
``` |
||

54 | ```
demand. Default value: 'demand'.
``` |
||

55 | ```
``` |
||

56 | ```
capacity : string
``` |
||

57 | ```
Edges of the graph G are expected to have an attribute capacity
``` |
||

58 | ```
that indicates how much flow the edge can support. If this
``` |
||

59 | ```
attribute is not present, the edge is considered to have
``` |
||

60 | ```
infinite capacity. Default value: 'capacity'.
``` |
||

61 | ```
``` |
||

62 | ```
weight : string
``` |
||

63 | ```
Edges of the graph G are expected to have an attribute weight
``` |
||

64 | ```
that indicates the cost incurred by sending one unit of flow on
``` |
||

65 | ```
that edge. If not present, the weight is considered to be 0.
``` |
||

66 | ```
Default value: 'weight'.
``` |
||

67 | ```
``` |
||

68 | ```
Returns
``` |
||

69 | ```
-------
``` |
||

70 | ```
flowCost : integer, float
``` |
||

71 | ```
Cost of a minimum cost flow satisfying all demands.
``` |
||

72 | ```
``` |
||

73 | ```
flowDict : dictionary
``` |
||

74 | ```
Dictionary of dictionaries keyed by nodes such that
``` |
||

75 | ```
flowDict[u][v] is the flow edge (u, v).
``` |
||

76 | ```
``` |
||

77 | ```
Raises
``` |
||

78 | ```
------
``` |
||

79 | ```
NetworkXError
``` |
||

80 | ```
This exception is raised if the input graph is not directed,
``` |
||

81 | ```
not connected or is a multigraph.
``` |
||

82 | ```
``` |
||

83 | ```
NetworkXUnfeasible
``` |
||

84 | ```
This exception is raised in the following situations:
``` |
||

85 | ```
``` |
||

86 | ```
* The sum of the demands is not zero. Then, there is no
``` |
||

87 | ```
flow satisfying all demands.
``` |
||

88 | ```
* There is no flow satisfying all demand.
``` |
||

89 | ```
``` |
||

90 | ```
NetworkXUnbounded
``` |
||

91 | ```
This exception is raised if the digraph G has a cycle of
``` |
||

92 | ```
negative cost and infinite capacity. Then, the cost of a flow
``` |
||

93 | ```
satisfying all demands is unbounded below.
``` |
||

94 | ```
``` |
||

95 | ```
Notes
``` |
||

96 | ```
-----
``` |
||

97 | ```
This algorithm is not guaranteed to work if edge weights or demands
``` |
||

98 | ```
are floating point numbers (overflows and roundoff errors can
``` |
||

99 | ```
cause problems). As a workaround you can use integer numbers by
``` |
||

100 | ```
multiplying the relevant edge attributes by a convenient
``` |
||

101 | ```
constant factor (eg 100).
``` |
||

102 | ```
``` |
||

103 | ```
See also
``` |
||

104 | ```
--------
``` |
||

105 | ```
cost_of_flow, max_flow_min_cost, min_cost_flow, min_cost_flow_cost
``` |
||

106 | ```
``` |
||

107 | ```
Examples
``` |
||

108 | ```
--------
``` |
||

109 | ```
A simple example of a min cost flow problem.
``` |
||

110 | ```
``` |
||

111 | ```
>>> import networkx as nx
``` |
||

112 | ```
>>> G = nx.DiGraph()
``` |
||

113 | ```
>>> G.add_node('a', demand=-5)
``` |
||

114 | ```
>>> G.add_node('d', demand=5)
``` |
||

115 | ```
>>> G.add_edge('a', 'b', weight=3, capacity=4)
``` |
||

116 | ```
>>> G.add_edge('a', 'c', weight=6, capacity=10)
``` |
||

117 | ```
>>> G.add_edge('b', 'd', weight=1, capacity=9)
``` |
||

118 | ```
>>> G.add_edge('c', 'd', weight=2, capacity=5)
``` |
||

119 | ```
>>> flowCost, flowDict = nx.network_simplex(G)
``` |
||

120 | ```
>>> flowCost
``` |
||

121 | ```
24
``` |
||

122 | ```
>>> flowDict # doctest: +SKIP
``` |
||

123 | ```
{'a': {'c': 1, 'b': 4}, 'c': {'d': 1}, 'b': {'d': 4}, 'd': {}}
``` |
||

124 | ```
``` |
||

125 | ```
The mincost flow algorithm can also be used to solve shortest path
``` |
||

126 | ```
problems. To find the shortest path between two nodes u and v,
``` |
||

127 | ```
give all edges an infinite capacity, give node u a demand of -1 and
``` |
||

128 | ```
node v a demand a 1. Then run the network simplex. The value of a
``` |
||

129 | ```
min cost flow will be the distance between u and v and edges
``` |
||

130 | ```
carrying positive flow will indicate the path.
``` |
||

131 | ```
``` |
||

132 | ```
>>> G=nx.DiGraph()
``` |
||

133 | ```
>>> G.add_weighted_edges_from([('s', 'u' ,10), ('s' ,'x' ,5),
``` |
||

134 | ```
... ('u', 'v' ,1), ('u' ,'x' ,2),
``` |
||

135 | ```
... ('v', 'y' ,1), ('x' ,'u' ,3),
``` |
||

136 | ```
... ('x', 'v' ,5), ('x' ,'y' ,2),
``` |
||

137 | ```
... ('y', 's' ,7), ('y' ,'v' ,6)])
``` |
||

138 | ```
>>> G.add_node('s', demand = -1)
``` |
||

139 | ```
>>> G.add_node('v', demand = 1)
``` |
||

140 | ```
>>> flowCost, flowDict = nx.network_simplex(G)
``` |
||

141 | ```
>>> flowCost == nx.shortest_path_length(G, 's', 'v', weight='weight')
``` |
||

142 | ```
True
``` |
||

143 | ```
>>> sorted([(u, v) for u in flowDict for v in flowDict[u] if flowDict[u][v] > 0])
``` |
||

144 | ```
[('s', 'x'), ('u', 'v'), ('x', 'u')]
``` |
||

145 | ```
>>> nx.shortest_path(G, 's', 'v', weight = 'weight')
``` |
||

146 | ```
['s', 'x', 'u', 'v']
``` |
||

147 | ```
``` |
||

148 | ```
It is possible to change the name of the attributes used for the
``` |
||

149 | ```
algorithm.
``` |
||

150 | ```
``` |
||

151 | ```
>>> G = nx.DiGraph()
``` |
||

152 | ```
>>> G.add_node('p', spam=-4)
``` |
||

153 | ```
>>> G.add_node('q', spam=2)
``` |
||

154 | ```
>>> G.add_node('a', spam=-2)
``` |
||

155 | ```
>>> G.add_node('d', spam=-1)
``` |
||

156 | ```
>>> G.add_node('t', spam=2)
``` |
||

157 | ```
>>> G.add_node('w', spam=3)
``` |
||

158 | ```
>>> G.add_edge('p', 'q', cost=7, vacancies=5)
``` |
||

159 | ```
>>> G.add_edge('p', 'a', cost=1, vacancies=4)
``` |
||

160 | ```
>>> G.add_edge('q', 'd', cost=2, vacancies=3)
``` |
||

161 | ```
>>> G.add_edge('t', 'q', cost=1, vacancies=2)
``` |
||

162 | ```
>>> G.add_edge('a', 't', cost=2, vacancies=4)
``` |
||

163 | ```
>>> G.add_edge('d', 'w', cost=3, vacancies=4)
``` |
||

164 | ```
>>> G.add_edge('t', 'w', cost=4, vacancies=1)
``` |
||

165 | ```
>>> flowCost, flowDict = nx.network_simplex(G, demand='spam',
``` |
||

166 | ```
... capacity='vacancies',
``` |
||

167 | ```
... weight='cost')
``` |
||

168 | ```
>>> flowCost
``` |
||

169 | ```
37
``` |
||

170 | ```
>>> flowDict # doctest: +SKIP
``` |
||

171 | ```
{'a': {'t': 4}, 'd': {'w': 2}, 'q': {'d': 1}, 'p': {'q': 2, 'a': 2}, 't': {'q': 1, 'w': 1}, 'w': {}}
``` |
||

172 | ```
``` |
||

173 | ```
References
``` |
||

174 | ```
----------
``` |
||

175 | ```
.. [1] Z. Kiraly, P. Kovacs.
``` |
||

176 | ```
Efficient implementation of minimum-cost flow algorithms.
``` |
||

177 | ```
Acta Universitatis Sapientiae, Informatica 4(1):67--118. 2012.
``` |
||

178 | ```
.. [2] R. Barr, F. Glover, D. Klingman.
``` |
||

179 | ```
Enhancement of spanning tree labeling procedures for network
``` |
||

180 | ```
optimization.
``` |
||

181 | ```
INFOR 17(1):16--34. 1979.
``` |
||

182 | ```
"""
``` |
||

183 | ```
###########################################################################
``` |
||

184 | ```
# Problem essentials extraction and sanity check
``` |
||

185 | ```
###########################################################################
``` |
||

186 | |||

187 | if len(G) == 0: |
||

188 | raise nx.NetworkXError('graph has no nodes') |
||

189 | |||

190 | ```
# Number all nodes and edges and hereafter reference them using ONLY their
``` |
||

191 | ```
# numbers
``` |
||

192 | |||

193 | N = list(G) # nodes |
||

194 | I = {u: i for i, u in enumerate(N)} # node indices |
||

195 | D = [G.nodes[u].get(demand, 0) for u in N] # node demands |
||

196 | |||

197 | inf = float('inf') |
||

198 | for p, b in zip(N, D): |
||

199 | if abs(b) == inf: |
||

200 | raise nx.NetworkXError('node %r has infinite demand' % (p,)) |
||

201 | |||

202 | multigraph = G.is_multigraph() |
||

203 | ```
S = [] # edge sources
``` |
||

204 | ```
T = [] # edge targets
``` |
||

205 | ```
if multigraph:
``` |
||

206 | ```
K = [] # edge keys
``` |
||

207 | ```
E = {} # edge indices
``` |
||

208 | ```
U = [] # edge capacities
``` |
||

209 | ```
C = [] # edge weights
``` |
||

210 | |||

211 | if not multigraph: |
||

212 | ```
edges = G.edges(data=True)
``` |
||

213 | ```
else:
``` |
||

214 | edges = G.edges(data=True, keys=True) |
||

215 | edges = (e for e in edges |
||

216 | if e[0] != e[1] and e[-1].get(capacity, inf) != 0) |
||

217 | for i, e in enumerate(edges): |
||

218 | ```
S.append(I[e[0]])
``` |
||

219 | ```
T.append(I[e[1]])
``` |
||

220 | ```
if multigraph:
``` |
||

221 | ```
K.append(e[2])
``` |
||

222 | ```
E[e[:-1]] = i
``` |
||

223 | ```
U.append(e[-1].get(capacity, inf))
``` |
||

224 | C.append(e[-1].get(weight, 0)) |
||

225 | |||

226 | for e, c in zip(E, C): |
||

227 | if abs(c) == inf: |
||

228 | raise nx.NetworkXError('edge %r has infinite weight' % (e,)) |
||

229 | if not multigraph: |
||

230 | ```
edges = nx.selfloop_edges(G, data=True)
``` |
||

231 | ```
else:
``` |
||

232 | edges = nx.selfloop_edges(G, data=True, keys=True) |
||

233 | for e in edges: |
||

234 | if abs(e[-1].get(weight, 0)) == inf: |
||

235 | raise nx.NetworkXError('edge %r has infinite weight' % (e[:-1],)) |
||

236 | |||

237 | ```
###########################################################################
``` |
||

238 | ```
# Quick infeasibility detection
``` |
||

239 | ```
###########################################################################
``` |
||

240 | |||

241 | if sum(D) != 0: |
||

242 | raise nx.NetworkXUnfeasible('total node demand is not zero') |
||

243 | for e, u in zip(E, U): |
||

244 | if u < 0: |
||

245 | raise nx.NetworkXUnfeasible('edge %r has negative capacity' % (e,)) |
||

246 | if not multigraph: |
||

247 | ```
edges = nx.selfloop_edges(G, data=True)
``` |
||

248 | ```
else:
``` |
||

249 | edges = nx.selfloop_edges(G, data=True, keys=True) |
||

250 | for e in edges: |
||

251 | if e[-1].get(capacity, inf) < 0: |
||

252 | ```
raise nx.NetworkXUnfeasible(
``` |
||

253 | 'edge %r has negative capacity' % (e[:-1],)) |
||

254 | |||

255 | ```
###########################################################################
``` |
||

256 | ```
# Initialization
``` |
||

257 | ```
###########################################################################
``` |
||

258 | |||

259 | ```
# Add a dummy node -1 and connect all existing nodes to it with infinite-
``` |
||

260 | ```
# capacity dummy edges. Node -1 will serve as the root of the
``` |
||

261 | ```
# spanning tree of the network simplex method. The new edges will used to
``` |
||

262 | ```
# trivially satisfy the node demands and create an initial strongly
``` |
||

263 | ```
# feasible spanning tree.
``` |
||

264 | n = len(N) # number of nodes |
||

265 | for p, d in enumerate(D): |
||

266 | if d > 0: # Must be greater-than here. Zero-demand nodes must have |
||

267 | ```
# edges pointing towards the root to ensure strong
``` |
||

268 | ```
# feasibility.
``` |
||

269 | ```
S.append(-1)
``` |
||

270 | T.append(p) |
||

271 | ```
else:
``` |
||

272 | S.append(p) |
||

273 | ```
T.append(-1)
``` |
||

274 | faux_inf = 3 * max(chain([sum(u for u in U if u < inf), |
||

275 | sum(abs(c) for c in C)], |
||

276 | (abs(d) for d in D))) or 1 |
||

277 | C.extend(repeat(faux_inf, n)) |
||

278 | U.extend(repeat(faux_inf, n)) |
||

279 | |||

280 | ```
# Construct the initial spanning tree.
``` |
||

281 | e = len(E) # number of edges |
||

282 | x = list(chain(repeat(0, e), (abs(d) for d in D))) # edge flows |
||

283 | pi = [faux_inf if d <= 0 else -faux_inf for d in D] # node potentials |
||

284 | parent = list(chain(repeat(-1, n), [None])) # parent nodes |
||

285 | edge = list(range(e, e + n)) # edges to parents |
||

286 | size = list(chain(repeat(1, n), [n + 1])) # subtree sizes |
||

287 | next = list(chain(range(1, n), [-1, 0])) # next nodes in depth-first thread |
||

288 | prev = list(range(-1, n)) # previous nodes in depth-first thread |
||

289 | last = list(chain(range(n), [n - 1])) # last descendants in depth-first thread |
||

290 | |||

291 | ```
###########################################################################
``` |
||

292 | ```
# Pivot loop
``` |
||

293 | ```
###########################################################################
``` |
||

294 | |||

295 | def reduced_cost(i): |
||

296 | ```
"""Returns the reduced cost of an edge i.
``` |
||

297 | ```
"""
``` |
||

298 | c = C[i] - pi[S[i]] + pi[T[i]] |
||

299 | return c if x[i] == 0 else -c |
||

300 | |||

301 | def find_entering_edges(): |
||

302 | ```
"""Yield entering edges until none can be found.
``` |
||

303 | ```
"""
``` |
||

304 | if e == 0: |
||

305 | ```
return
``` |
||

306 | |||

307 | ```
# Entering edges are found by combining Dantzig's rule and Bland's
``` |
||

308 | ```
# rule. The edges are cyclically grouped into blocks of size B. Within
``` |
||

309 | ```
# each block, Dantzig's rule is applied to find an entering edge. The
``` |
||

310 | ```
# blocks to search is determined following Bland's rule.
``` |
||

311 | B = int(ceil(sqrt(e))) # pivot block size |
||

312 | M = (e + B - 1) // B # number of blocks needed to cover all edges |
||

313 | m = 0 # number of consecutive blocks without eligible |
||

314 | ```
# entering edges
``` |
||

315 | f = 0 # first edge in block |
||

316 | ```
while m < M:
``` |
||

317 | ```
# Determine the next block of edges.
``` |
||

318 | l = f + B |
||

319 | ```
if l <= e:
``` |
||

320 | ```
edges = range(f, l)
``` |
||

321 | ```
else:
``` |
||

322 | l -= e |
||

323 | edges = chain(range(f, e), range(l)) |
||

324 | f = l |
||

325 | ```
# Find the first edge with the lowest reduced cost.
``` |
||

326 | ```
i = min(edges, key=reduced_cost)
``` |
||

327 | c = reduced_cost(i) |
||

328 | if c >= 0: |
||

329 | ```
# No entering edge found in the current block.
``` |
||

330 | ```
m += 1
``` |
||

331 | ```
else:
``` |
||

332 | ```
# Entering edge found.
``` |
||

333 | if x[i] == 0: |
||

334 | p = S[i] |
||

335 | q = T[i] |
||

336 | ```
else:
``` |
||

337 | p = T[i] |
||

338 | q = S[i] |
||

339 | ```
yield i, p, q
``` |
||

340 | ```
m = 0
``` |
||

341 | ```
# All edges have nonnegative reduced costs. The current flow is
``` |
||

342 | ```
# optimal.
``` |
||

343 | |||

344 | def find_apex(p, q): |
||

345 | ```
"""Find the lowest common ancestor of nodes p and q in the spanning
``` |
||

346 | ```
tree.
``` |
||

347 | ```
"""
``` |
||

348 | size_p = size[p] |
||

349 | size_q = size[q] |
||

350 | while True: |
||

351 | ```
while size_p < size_q:
``` |
||

352 | p = parent[p] |
||

353 | size_p = size[p] |
||

354 | ```
while size_p > size_q:
``` |
||

355 | q = parent[q] |
||

356 | size_q = size[q] |
||

357 | ```
if size_p == size_q:
``` |
||

358 | ```
if p != q:
``` |
||

359 | p = parent[p] |
||

360 | size_p = size[p] |
||

361 | q = parent[q] |
||

362 | size_q = size[q] |
||

363 | ```
else:
``` |
||

364 | ```
return p
``` |
||

365 | |||

366 | def trace_path(p, w): |
||

367 | ```
"""Returns the nodes and edges on the path from node p to its ancestor
``` |
||

368 | ```
w.
``` |
||

369 | ```
"""
``` |
||

370 | Wn = [p] |
||

371 | We = [] |
||

372 | ```
while p != w:
``` |
||

373 | We.append(edge[p]) |
||

374 | p = parent[p] |
||

375 | Wn.append(p) |
||

376 | ```
return Wn, We
``` |
||

377 | |||

378 | def find_cycle(i, p, q): |
||

379 | ```
"""Returns the nodes and edges on the cycle containing edge i == (p, q)
``` |
||

380 | ```
when the latter is added to the spanning tree.
``` |
||

381 | ```
``` |
||

382 | ```
The cycle is oriented in the direction from p to q.
``` |
||

383 | ```
"""
``` |
||

384 | w = find_apex(p, q) |
||

385 | Wn, We = trace_path(p, w) |
||

386 | Wn.reverse() |
||

387 | We.reverse() |
||

388 | ```
if We != [i]:
``` |
||

389 | We.append(i) |
||

390 | WnR, WeR = trace_path(q, w) |
||

391 | del WnR[-1] |
||

392 | Wn += WnR |
||

393 | We += WeR |
||

394 | ```
return Wn, We
``` |
||

395 | |||

396 | def residual_capacity(i, p): |
||

397 | ```
"""Returns the residual capacity of an edge i in the direction away
``` |
||

398 | ```
from its endpoint p.
``` |
||

399 | ```
"""
``` |
||

400 | return U[i] - x[i] if S[i] == p else x[i] |
||

401 | |||

402 | def find_leaving_edge(Wn, We): |
||

403 | ```
"""Returns the leaving edge in a cycle represented by Wn and We.
``` |
||

404 | ```
"""
``` |
||

405 | j, s = min(zip(reversed(We), reversed(Wn)), |
||

406 | ```
key=lambda i_p: residual_capacity(*i_p))
``` |
||

407 | t = T[j] if S[j] == s else S[j] |
||

408 | ```
return j, s, t
``` |
||

409 | |||

410 | def augment_flow(Wn, We, f): |
||

411 | ```
"""Augment f units of flow along a cycle represented by Wn and We.
``` |
||

412 | ```
"""
``` |
||

413 | for i, p in zip(We, Wn): |
||

414 | ```
if S[i] == p:
``` |
||

415 | x[i] += f |
||

416 | ```
else:
``` |
||

417 | x[i] -= f |
||

418 | |||

419 | def trace_subtree(p): |
||

420 | ```
"""Yield the nodes in the subtree rooted at a node p.
``` |
||

421 | ```
"""
``` |
||

422 | ```
yield p
``` |
||

423 | l = last[p] |
||

424 | ```
while p != l:
``` |
||

425 | ```
p = next[p]
``` |
||

426 | ```
yield p
``` |
||

427 | |||

428 | def remove_edge(s, t): |
||

429 | ```
"""Remove an edge (s, t) where parent[t] == s from the spanning tree.
``` |
||

430 | ```
"""
``` |
||

431 | size_t = size[t] |
||

432 | prev_t = prev[t] |
||

433 | last_t = last[t] |
||

434 | ```
next_last_t = next[last_t]
``` |
||

435 | ```
# Remove (s, t).
``` |
||

436 | ```
parent[t] = None
``` |
||

437 | ```
edge[t] = None
``` |
||

438 | ```
# Remove the subtree rooted at t from the depth-first thread.
``` |
||

439 | ```
next[prev_t] = next_last_t
``` |
||

440 | prev[next_last_t] = prev_t |
||

441 | ```
next[last_t] = t
``` |
||

442 | prev[t] = last_t |
||

443 | ```
# Update the subtree sizes and last descendants of the (old) acenstors
``` |
||

444 | ```
# of t.
``` |
||

445 | while s is not None: |
||

446 | size[s] -= size_t |
||

447 | ```
if last[s] == last_t:
``` |
||

448 | last[s] = prev_t |
||

449 | s = parent[s] |
||

450 | |||

451 | def make_root(q): |
||

452 | ```
"""Make a node q the root of its containing subtree.
``` |
||

453 | ```
"""
``` |
||

454 | ancestors = [] |
||

455 | while q is not None: |
||

456 | ancestors.append(q) |
||

457 | q = parent[q] |
||

458 | ancestors.reverse() |
||

459 | for p, q in zip(ancestors, islice(ancestors, 1, None)): |
||

460 | size_p = size[p] |
||

461 | last_p = last[p] |
||

462 | prev_q = prev[q] |
||

463 | last_q = last[q] |
||

464 | ```
next_last_q = next[last_q]
``` |
||

465 | ```
# Make p a child of q.
``` |
||

466 | parent[p] = q |
||

467 | ```
parent[q] = None
``` |
||

468 | edge[p] = edge[q] |
||

469 | ```
edge[q] = None
``` |
||

470 | size[p] = size_p - size[q] |
||

471 | size[q] = size_p |
||

472 | ```
# Remove the subtree rooted at q from the depth-first thread.
``` |
||

473 | ```
next[prev_q] = next_last_q
``` |
||

474 | prev[next_last_q] = prev_q |
||

475 | ```
next[last_q] = q
``` |
||

476 | prev[q] = last_q |
||

477 | ```
if last_p == last_q:
``` |
||

478 | last[p] = prev_q |
||

479 | last_p = prev_q |
||

480 | ```
# Add the remaining parts of the subtree rooted at p as a subtree
``` |
||

481 | ```
# of q in the depth-first thread.
``` |
||

482 | prev[p] = last_q |
||

483 | ```
next[last_q] = p
``` |
||

484 | ```
next[last_p] = q
``` |
||

485 | prev[q] = last_p |
||

486 | last[q] = last_p |
||

487 | |||

488 | def add_edge(i, p, q): |
||

489 | ```
"""Add an edge (p, q) to the spanning tree where q is the root of a
``` |
||

490 | ```
subtree.
``` |
||

491 | ```
"""
``` |
||

492 | last_p = last[p] |
||

493 | ```
next_last_p = next[last_p]
``` |
||

494 | size_q = size[q] |
||

495 | last_q = last[q] |
||

496 | ```
# Make q a child of p.
``` |
||

497 | parent[q] = p |
||

498 | edge[q] = i |
||

499 | ```
# Insert the subtree rooted at q into the depth-first thread.
``` |
||

500 | ```
next[last_p] = q
``` |
||

501 | prev[q] = last_p |
||

502 | prev[next_last_p] = last_q |
||

503 | ```
next[last_q] = next_last_p
``` |
||

504 | ```
# Update the subtree sizes and last descendants of the (new) ancestors
``` |
||

505 | ```
# of q.
``` |
||

506 | while p is not None: |
||

507 | size[p] += size_q |
||

508 | ```
if last[p] == last_p:
``` |
||

509 | last[p] = last_q |
||

510 | p = parent[p] |
||

511 | |||

512 | def update_potentials(i, p, q): |
||

513 | ```
"""Update the potentials of the nodes in the subtree rooted at a node
``` |
||

514 | ```
q connected to its parent p by an edge i.
``` |
||

515 | ```
"""
``` |
||

516 | ```
if q == T[i]:
``` |
||

517 | d = pi[p] - C[i] - pi[q] |
||

518 | ```
else:
``` |
||

519 | d = pi[p] + C[i] - pi[q] |
||

520 | for q in trace_subtree(q): |
||

521 | pi[q] += d |
||

522 | |||

523 | ```
# Pivot loop
``` |
||

524 | for i, p, q in find_entering_edges(): |
||

525 | Wn, We = find_cycle(i, p, q) |
||

526 | j, s, t = find_leaving_edge(Wn, We) |
||

527 | augment_flow(Wn, We, residual_capacity(j, s)) |
||

528 | if i != j: # Do nothing more if the entering edge is the same as the |
||

529 | ```
# the leaving edge.
``` |
||

530 | ```
if parent[t] != s:
``` |
||

531 | ```
# Ensure that s is the parent of t.
``` |
||

532 | s, t = t, s |
||

533 | ```
if We.index(i) > We.index(j):
``` |
||

534 | ```
# Ensure that q is in the subtree rooted at t.
``` |
||

535 | p, q = q, p |
||

536 | remove_edge(s, t) |
||

537 | make_root(q) |
||

538 | add_edge(i, p, q) |
||

539 | update_potentials(i, p, q) |
||

540 | |||

541 | ```
###########################################################################
``` |
||

542 | ```
# Infeasibility and unboundedness detection
``` |
||

543 | ```
###########################################################################
``` |
||

544 | |||

545 | if any(x[i] != 0 for i in range(-n, 0)): |
||

546 | raise nx.NetworkXUnfeasible('no flow satisfies all node demands') |
||

547 | |||

548 | if (any(x[i] * 2 >= faux_inf for i in range(e)) or |
||

549 | any(e[-1].get(capacity, inf) == inf and e[-1].get(weight, 0) < 0 |
||

550 | for e in nx.selfloop_edges(G, data=True))): |
||

551 | ```
raise nx.NetworkXUnbounded(
``` |
||

552 | ```
'negative cycle with infinite capacity found')
``` |
||

553 | |||

554 | ```
###########################################################################
``` |
||

555 | ```
# Flow cost calculation and flow dict construction
``` |
||

556 | ```
###########################################################################
``` |
||

557 | |||

558 | ```
del x[e:]
``` |
||

559 | flow_cost = sum(c * x for c, x in zip(C, x)) |
||

560 | flow_dict = {n: {} for n in N} |
||

561 | |||

562 | def add_entry(e): |
||

563 | ```
"""Add a flow dict entry.
``` |
||

564 | ```
"""
``` |
||

565 | ```
d = flow_dict[e[0]]
``` |
||

566 | for k in e[1:-2]: |
||

567 | ```
try:
``` |
||

568 | d = d[k] |
||

569 | except KeyError: |
||

570 | t = {} |
||

571 | d[k] = t |
||

572 | d = t |
||

573 | d[e[-2]] = e[-1] |
||

574 | |||

575 | S = (N[s] for s in S) # Use original nodes. |
||

576 | T = (N[t] for t in T) # Use original nodes. |
||

577 | if not multigraph: |
||

578 | for e in zip(S, T, x): |
||

579 | add_entry(e) |
||

580 | ```
edges = G.edges(data=True)
``` |
||

581 | ```
else:
``` |
||

582 | for e in zip(S, T, K, x): |
||

583 | add_entry(e) |
||

584 | edges = G.edges(data=True, keys=True) |
||

585 | for e in edges: |
||

586 | if e[0] != e[1]: |
||

587 | if e[-1].get(capacity, inf) == 0: |
||

588 | add_entry(e[:-1] + (0,)) |
||

589 | ```
else:
``` |
||

590 | c = e[-1].get(weight, 0) |
||

591 | if c >= 0: |
||

592 | add_entry(e[:-1] + (0,)) |
||

593 | ```
else:
``` |
||

594 | ```
u = e[-1][capacity]
``` |
||

595 | flow_cost += c * u |
||

596 | ```
add_entry(e[:-1] + (u,))
``` |
||

597 | |||

598 | ` return flow_cost, flow_dict` |