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"""
2Minimum cost flow algorithms on directed connected graphs.
3"""
5__all__ = ["network_simplex"]
7from itertools import chain, islice, repeat
8from math import ceil, sqrt
10import networkx as nx
11from networkx.utils import not_implemented_for
14class _DataEssentialsAndFunctions:
15 def __init__(
16 self, G, multigraph, demand="demand", capacity="capacity", weight="weight"
17 ):
18 # Number all nodes and edges and hereafter reference them using ONLY their numbers
19 self.node_list = list(G) # nodes
20 self.node_indices = {u: i for i, u in enumerate(self.node_list)} # node indices
21 self.node_demands = [
22 G.nodes[u].get(demand, 0) for u in self.node_list
23 ] # node demands
25 self.edge_sources = [] # edge sources
26 self.edge_targets = [] # edge targets
27 if multigraph:
28 self.edge_keys = [] # edge keys
29 self.edge_indices = {} # edge indices
30 self.edge_capacities = [] # edge capacities
31 self.edge_weights = [] # edge weights
33 if not multigraph:
34 edges = G.edges(data=True)
35 else:
36 edges = G.edges(data=True, keys=True)
38 inf = float("inf")
39 edges = (e for e in edges if e[0] != e[1] and e[-1].get(capacity, inf) != 0)
40 for i, e in enumerate(edges):
41 self.edge_sources.append(self.node_indices[e[0]])
42 self.edge_targets.append(self.node_indices[e[1]])
43 if multigraph:
44 self.edge_keys.append(e[2])
45 self.edge_indices[e[:-1]] = i
46 self.edge_capacities.append(e[-1].get(capacity, inf))
47 self.edge_weights.append(e[-1].get(weight, 0))
49 # spanning tree specific data to be initialized
51 self.edge_count = None # number of edges
52 self.edge_flow = None # edge flows
53 self.node_potentials = None # node potentials
54 self.parent = None # parent nodes
55 self.parent_edge = None # edges to parents
56 self.subtree_size = None # subtree sizes
57 self.next_node_dft = None # next nodes in depth-first thread
58 self.prev_node_dft = None # previous nodes in depth-first thread
59 self.last_descendent_dft = None # last descendants in depth-first thread
60 self._spanning_tree_initialized = (
61 False # False until initialize_spanning_tree() is called
62 )
64 def initialize_spanning_tree(self, n, faux_inf):
65 self.edge_count = len(self.edge_indices) # number of edges
66 self.edge_flow = list(
67 chain(repeat(0, self.edge_count), (abs(d) for d in self.node_demands))
68 ) # edge flows
69 self.node_potentials = [
70 faux_inf if d <= 0 else -faux_inf for d in self.node_demands
71 ] # node potentials
72 self.parent = list(chain(repeat(-1, n), [None])) # parent nodes
73 self.parent_edge = list(
74 range(self.edge_count, self.edge_count + n)
75 ) # edges to parents
76 self.subtree_size = list(chain(repeat(1, n), [n + 1])) # subtree sizes
77 self.next_node_dft = list(
78 chain(range(1, n), [-1, 0])
79 ) # next nodes in depth-first thread
80 self.prev_node_dft = list(range(-1, n)) # previous nodes in depth-first thread
81 self.last_descendent_dft = list(
82 chain(range(n), [n - 1])
83 ) # last descendants in depth-first thread
84 self._spanning_tree_initialized = True # True only if all the assignments pass
86 def find_apex(self, p, q):
87 """
88 Find the lowest common ancestor of nodes p and q in the spanning tree.
89 """
90 size_p = self.subtree_size[p]
91 size_q = self.subtree_size[q]
92 while True:
93 while size_p < size_q:
94 p = self.parent[p]
95 size_p = self.subtree_size[p]
96 while size_p > size_q:
97 q = self.parent[q]
98 size_q = self.subtree_size[q]
99 if size_p == size_q:
100 if p != q:
101 p = self.parent[p]
102 size_p = self.subtree_size[p]
103 q = self.parent[q]
104 size_q = self.subtree_size[q]
105 else:
106 return p
108 def trace_path(self, p, w):
109 """
110 Returns the nodes and edges on the path from node p to its ancestor w.
111 """
112 Wn = [p]
113 We = []
114 while p != w:
115 We.append(self.parent_edge[p])
116 p = self.parent[p]
117 Wn.append(p)
118 return Wn, We
120 def find_cycle(self, i, p, q):
121 """
122 Returns the nodes and edges on the cycle containing edge i == (p, q)
123 when the latter is added to the spanning tree.
125 The cycle is oriented in the direction from p to q.
126 """
127 w = self.find_apex(p, q)
128 Wn, We = self.trace_path(p, w)
129 Wn.reverse()
130 We.reverse()
131 if We != [i]:
132 We.append(i)
133 WnR, WeR = self.trace_path(q, w)
134 del WnR[-1]
135 Wn += WnR
136 We += WeR
137 return Wn, We
139 def augment_flow(self, Wn, We, f):
140 """
141 Augment f units of flow along a cycle represented by Wn and We.
142 """
143 for i, p in zip(We, Wn):
144 if self.edge_sources[i] == p:
145 self.edge_flow[i] += f
146 else:
147 self.edge_flow[i] -= f
149 def trace_subtree(self, p):
150 """
151 Yield the nodes in the subtree rooted at a node p.
152 """
153 yield p
154 l = self.last_descendent_dft[p]
155 while p != l:
156 p = self.next_node_dft[p]
157 yield p
159 def remove_edge(self, s, t):
160 """
161 Remove an edge (s, t) where parent[t] == s from the spanning tree.
162 """
163 size_t = self.subtree_size[t]
164 prev_t = self.prev_node_dft[t]
165 last_t = self.last_descendent_dft[t]
166 next_last_t = self.next_node_dft[last_t]
167 # Remove (s, t).
168 self.parent[t] = None
169 self.parent_edge[t] = None
170 # Remove the subtree rooted at t from the depth-first thread.
171 self.next_node_dft[prev_t] = next_last_t
172 self.prev_node_dft[next_last_t] = prev_t
173 self.next_node_dft[last_t] = t
174 self.prev_node_dft[t] = last_t
175 # Update the subtree sizes and last descendants of the (old) ancestors
176 # of t.
177 while s is not None:
178 self.subtree_size[s] -= size_t
179 if self.last_descendent_dft[s] == last_t:
180 self.last_descendent_dft[s] = prev_t
181 s = self.parent[s]
183 def make_root(self, q):
184 """
185 Make a node q the root of its containing subtree.
186 """
187 ancestors = []
188 while q is not None:
189 ancestors.append(q)
190 q = self.parent[q]
191 ancestors.reverse()
192 for p, q in zip(ancestors, islice(ancestors, 1, None)):
193 size_p = self.subtree_size[p]
194 last_p = self.last_descendent_dft[p]
195 prev_q = self.prev_node_dft[q]
196 last_q = self.last_descendent_dft[q]
197 next_last_q = self.next_node_dft[last_q]
198 # Make p a child of q.
199 self.parent[p] = q
200 self.parent[q] = None
201 self.parent_edge[p] = self.parent_edge[q]
202 self.parent_edge[q] = None
203 self.subtree_size[p] = size_p - self.subtree_size[q]
204 self.subtree_size[q] = size_p
205 # Remove the subtree rooted at q from the depth-first thread.
206 self.next_node_dft[prev_q] = next_last_q
207 self.prev_node_dft[next_last_q] = prev_q
208 self.next_node_dft[last_q] = q
209 self.prev_node_dft[q] = last_q
210 if last_p == last_q:
211 self.last_descendent_dft[p] = prev_q
212 last_p = prev_q
213 # Add the remaining parts of the subtree rooted at p as a subtree
214 # of q in the depth-first thread.
215 self.prev_node_dft[p] = last_q
216 self.next_node_dft[last_q] = p
217 self.next_node_dft[last_p] = q
218 self.prev_node_dft[q] = last_p
219 self.last_descendent_dft[q] = last_p
221 def add_edge(self, i, p, q):
222 """
223 Add an edge (p, q) to the spanning tree where q is the root of a subtree.
224 """
225 last_p = self.last_descendent_dft[p]
226 next_last_p = self.next_node_dft[last_p]
227 size_q = self.subtree_size[q]
228 last_q = self.last_descendent_dft[q]
229 # Make q a child of p.
230 self.parent[q] = p
231 self.parent_edge[q] = i
232 # Insert the subtree rooted at q into the depth-first thread.
233 self.next_node_dft[last_p] = q
234 self.prev_node_dft[q] = last_p
235 self.prev_node_dft[next_last_p] = last_q
236 self.next_node_dft[last_q] = next_last_p
237 # Update the subtree sizes and last descendants of the (new) ancestors
238 # of q.
239 while p is not None:
240 self.subtree_size[p] += size_q
241 if self.last_descendent_dft[p] == last_p:
242 self.last_descendent_dft[p] = last_q
243 p = self.parent[p]
245 def update_potentials(self, i, p, q):
246 """
247 Update the potentials of the nodes in the subtree rooted at a node
248 q connected to its parent p by an edge i.
249 """
250 if q == self.edge_targets[i]:
251 d = self.node_potentials[p] - self.edge_weights[i] - self.node_potentials[q]
252 else:
253 d = self.node_potentials[p] + self.edge_weights[i] - self.node_potentials[q]
254 for q in self.trace_subtree(q):
255 self.node_potentials[q] += d
257 def reduced_cost(self, i):
258 """Returns the reduced cost of an edge i."""
259 c = (
260 self.edge_weights[i]
261 - self.node_potentials[self.edge_sources[i]]
262 + self.node_potentials[self.edge_targets[i]]
263 )
264 return c if self.edge_flow[i] == 0 else -c
266 def find_entering_edges(self):
267 """Yield entering edges until none can be found."""
268 if self.edge_count == 0:
269 return
271 # Entering edges are found by combining Dantzig's rule and Bland's
272 # rule. The edges are cyclically grouped into blocks of size B. Within
273 # each block, Dantzig's rule is applied to find an entering edge. The
274 # blocks to search is determined following Bland's rule.
275 B = int(ceil(sqrt(self.edge_count))) # pivot block size
276 M = (self.edge_count + B - 1) // B # number of blocks needed to cover all edges
277 m = 0 # number of consecutive blocks without eligible
278 # entering edges
279 f = 0 # first edge in block
280 while m < M:
281 # Determine the next block of edges.
282 l = f + B
283 if l <= self.edge_count:
284 edges = range(f, l)
285 else:
286 l -= self.edge_count
287 edges = chain(range(f, self.edge_count), range(l))
288 f = l
289 # Find the first edge with the lowest reduced cost.
290 i = min(edges, key=self.reduced_cost)
291 c = self.reduced_cost(i)
292 if c >= 0:
293 # No entering edge found in the current block.
294 m += 1
295 else:
296 # Entering edge found.
297 if self.edge_flow[i] == 0:
298 p = self.edge_sources[i]
299 q = self.edge_targets[i]
300 else:
301 p = self.edge_targets[i]
302 q = self.edge_sources[i]
303 yield i, p, q
304 m = 0
305 # All edges have nonnegative reduced costs. The current flow is
306 # optimal.
308 def residual_capacity(self, i, p):
309 """Returns the residual capacity of an edge i in the direction away
310 from its endpoint p.
311 """
312 return (
313 self.edge_capacities[i] - self.edge_flow[i]
314 if self.edge_sources[i] == p
315 else self.edge_flow[i]
316 )
318 def find_leaving_edge(self, Wn, We):
319 """Returns the leaving edge in a cycle represented by Wn and We."""
320 j, s = min(
321 zip(reversed(We), reversed(Wn)),
322 key=lambda i_p: self.residual_capacity(*i_p),
323 )
324 t = self.edge_targets[j] if self.edge_sources[j] == s else self.edge_sources[j]
325 return j, s, t
328@not_implemented_for("undirected")
329@nx._dispatchable(
330 node_attrs="demand", edge_attrs={"capacity": float("inf"), "weight": 0}
331)
332def network_simplex(G, demand="demand", capacity="capacity", weight="weight"):
333 r"""Find a minimum cost flow satisfying all demands in digraph G.
335 This is a primal network simplex algorithm that uses the leaving
336 arc rule to prevent cycling.
338 G is a digraph with edge costs and capacities and in which nodes
339 have demand, i.e., they want to send or receive some amount of
340 flow. A negative demand means that the node wants to send flow, a
341 positive demand means that the node want to receive flow. A flow on
342 the digraph G satisfies all demand if the net flow into each node
343 is equal to the demand of that node.
345 Parameters
346 ----------
347 G : NetworkX graph
348 DiGraph on which a minimum cost flow satisfying all demands is
349 to be found.
351 demand : string
352 Nodes of the graph G are expected to have an attribute demand
353 that indicates how much flow a node wants to send (negative
354 demand) or receive (positive demand). Note that the sum of the
355 demands should be 0 otherwise the problem in not feasible. If
356 this attribute is not present, a node is considered to have 0
357 demand. Default value: 'demand'.
359 capacity : string
360 Edges of the graph G are expected to have an attribute capacity
361 that indicates how much flow the edge can support. If this
362 attribute is not present, the edge is considered to have
363 infinite capacity. Default value: 'capacity'.
365 weight : string
366 Edges of the graph G are expected to have an attribute weight
367 that indicates the cost incurred by sending one unit of flow on
368 that edge. If not present, the weight is considered to be 0.
369 Default value: 'weight'.
371 Returns
372 -------
373 flowCost : integer, float
374 Cost of a minimum cost flow satisfying all demands.
376 flowDict : dictionary
377 Dictionary of dictionaries keyed by nodes such that
378 flowDict[u][v] is the flow edge (u, v).
380 Raises
381 ------
382 NetworkXError
383 This exception is raised if the input graph is not directed or
384 not connected.
386 NetworkXUnfeasible
387 This exception is raised in the following situations:
389 * The sum of the demands is not zero. Then, there is no
390 flow satisfying all demands.
391 * There is no flow satisfying all demand.
393 NetworkXUnbounded
394 This exception is raised if the digraph G has a cycle of
395 negative cost and infinite capacity. Then, the cost of a flow
396 satisfying all demands is unbounded below.
398 Notes
399 -----
400 This algorithm is not guaranteed to work if edge weights or demands
401 are floating point numbers (overflows and roundoff errors can
402 cause problems). As a workaround you can use integer numbers by
403 multiplying the relevant edge attributes by a convenient
404 constant factor (eg 100).
406 See also
407 --------
408 cost_of_flow, max_flow_min_cost, min_cost_flow, min_cost_flow_cost
410 Examples
411 --------
412 A simple example of a min cost flow problem.
414 >>> G = nx.DiGraph()
415 >>> G.add_node("a", demand=-5)
416 >>> G.add_node("d", demand=5)
417 >>> G.add_edge("a", "b", weight=3, capacity=4)
418 >>> G.add_edge("a", "c", weight=6, capacity=10)
419 >>> G.add_edge("b", "d", weight=1, capacity=9)
420 >>> G.add_edge("c", "d", weight=2, capacity=5)
421 >>> flowCost, flowDict = nx.network_simplex(G)
422 >>> flowCost
423 24
424 >>> flowDict
425 {'a': {'b': 4, 'c': 1}, 'd': {}, 'b': {'d': 4}, 'c': {'d': 1}}
427 The mincost flow algorithm can also be used to solve shortest path
428 problems. To find the shortest path between two nodes u and v,
429 give all edges an infinite capacity, give node u a demand of -1 and
430 node v a demand a 1. Then run the network simplex. The value of a
431 min cost flow will be the distance between u and v and edges
432 carrying positive flow will indicate the path.
434 >>> G = nx.DiGraph()
435 >>> G.add_weighted_edges_from(
436 ... [
437 ... ("s", "u", 10),
438 ... ("s", "x", 5),
439 ... ("u", "v", 1),
440 ... ("u", "x", 2),
441 ... ("v", "y", 1),
442 ... ("x", "u", 3),
443 ... ("x", "v", 5),
444 ... ("x", "y", 2),
445 ... ("y", "s", 7),
446 ... ("y", "v", 6),
447 ... ]
448 ... )
449 >>> G.add_node("s", demand=-1)
450 >>> G.add_node("v", demand=1)
451 >>> flowCost, flowDict = nx.network_simplex(G)
452 >>> flowCost == nx.shortest_path_length(G, "s", "v", weight="weight")
453 True
454 >>> sorted([(u, v) for u in flowDict for v in flowDict[u] if flowDict[u][v] > 0])
455 [('s', 'x'), ('u', 'v'), ('x', 'u')]
456 >>> nx.shortest_path(G, "s", "v", weight="weight")
457 ['s', 'x', 'u', 'v']
459 It is possible to change the name of the attributes used for the
460 algorithm.
462 >>> G = nx.DiGraph()
463 >>> G.add_node("p", spam=-4)
464 >>> G.add_node("q", spam=2)
465 >>> G.add_node("a", spam=-2)
466 >>> G.add_node("d", spam=-1)
467 >>> G.add_node("t", spam=2)
468 >>> G.add_node("w", spam=3)
469 >>> G.add_edge("p", "q", cost=7, vacancies=5)
470 >>> G.add_edge("p", "a", cost=1, vacancies=4)
471 >>> G.add_edge("q", "d", cost=2, vacancies=3)
472 >>> G.add_edge("t", "q", cost=1, vacancies=2)
473 >>> G.add_edge("a", "t", cost=2, vacancies=4)
474 >>> G.add_edge("d", "w", cost=3, vacancies=4)
475 >>> G.add_edge("t", "w", cost=4, vacancies=1)
476 >>> flowCost, flowDict = nx.network_simplex(
477 ... G, demand="spam", capacity="vacancies", weight="cost"
478 ... )
479 >>> flowCost
480 37
481 >>> flowDict
482 {'p': {'q': 2, 'a': 2}, 'q': {'d': 1}, 'a': {'t': 4}, 'd': {'w': 2}, 't': {'q': 1, 'w': 1}, 'w': {}}
484 References
485 ----------
486 .. [1] Z. Kiraly, P. Kovacs.
487 Efficient implementation of minimum-cost flow algorithms.
488 Acta Universitatis Sapientiae, Informatica 4(1):67--118. 2012.
489 .. [2] R. Barr, F. Glover, D. Klingman.
490 Enhancement of spanning tree labeling procedures for network
491 optimization.
492 INFOR 17(1):16--34. 1979.
493 """
494 ###########################################################################
495 # Problem essentials extraction and sanity check
496 ###########################################################################
498 if len(G) == 0:
499 raise nx.NetworkXError("graph has no nodes")
501 multigraph = G.is_multigraph()
503 # extracting data essential to problem
504 DEAF = _DataEssentialsAndFunctions(
505 G, multigraph, demand=demand, capacity=capacity, weight=weight
506 )
508 ###########################################################################
509 # Quick Error Detection
510 ###########################################################################
512 inf = float("inf")
513 for u, d in zip(DEAF.node_list, DEAF.node_demands):
514 if abs(d) == inf:
515 raise nx.NetworkXError(f"node {u!r} has infinite demand")
516 for e, w in zip(DEAF.edge_indices, DEAF.edge_weights):
517 if abs(w) == inf:
518 raise nx.NetworkXError(f"edge {e!r} has infinite weight")
519 if not multigraph:
520 edges = nx.selfloop_edges(G, data=True)
521 else:
522 edges = nx.selfloop_edges(G, data=True, keys=True)
523 for e in edges:
524 if abs(e[-1].get(weight, 0)) == inf:
525 raise nx.NetworkXError(f"edge {e[:-1]!r} has infinite weight")
527 ###########################################################################
528 # Quick Infeasibility Detection
529 ###########################################################################
531 if sum(DEAF.node_demands) != 0:
532 raise nx.NetworkXUnfeasible("total node demand is not zero")
533 for e, c in zip(DEAF.edge_indices, DEAF.edge_capacities):
534 if c < 0:
535 raise nx.NetworkXUnfeasible(f"edge {e!r} has negative capacity")
536 if not multigraph:
537 edges = nx.selfloop_edges(G, data=True)
538 else:
539 edges = nx.selfloop_edges(G, data=True, keys=True)
540 for e in edges:
541 if e[-1].get(capacity, inf) < 0:
542 raise nx.NetworkXUnfeasible(f"edge {e[:-1]!r} has negative capacity")
544 ###########################################################################
545 # Initialization
546 ###########################################################################
548 # Add a dummy node -1 and connect all existing nodes to it with infinite-
549 # capacity dummy edges. Node -1 will serve as the root of the
550 # spanning tree of the network simplex method. The new edges will used to
551 # trivially satisfy the node demands and create an initial strongly
552 # feasible spanning tree.
553 for i, d in enumerate(DEAF.node_demands):
554 # Must be greater-than here. Zero-demand nodes must have
555 # edges pointing towards the root to ensure strong feasibility.
556 if d > 0:
557 DEAF.edge_sources.append(-1)
558 DEAF.edge_targets.append(i)
559 else:
560 DEAF.edge_sources.append(i)
561 DEAF.edge_targets.append(-1)
562 faux_inf = (
563 3
564 * max(
565 sum(c for c in DEAF.edge_capacities if c < inf),
566 sum(abs(w) for w in DEAF.edge_weights),
567 sum(abs(d) for d in DEAF.node_demands),
568 )
569 or 1
570 )
572 n = len(DEAF.node_list) # number of nodes
573 DEAF.edge_weights.extend(repeat(faux_inf, n))
574 DEAF.edge_capacities.extend(repeat(faux_inf, n))
576 # Construct the initial spanning tree.
577 DEAF.initialize_spanning_tree(n, faux_inf)
579 ###########################################################################
580 # Pivot loop
581 ###########################################################################
583 for i, p, q in DEAF.find_entering_edges():
584 Wn, We = DEAF.find_cycle(i, p, q)
585 j, s, t = DEAF.find_leaving_edge(Wn, We)
586 DEAF.augment_flow(Wn, We, DEAF.residual_capacity(j, s))
587 # Do nothing more if the entering edge is the same as the leaving edge.
588 if i != j:
589 if DEAF.parent[t] != s:
590 # Ensure that s is the parent of t.
591 s, t = t, s
592 if We.index(i) > We.index(j):
593 # Ensure that q is in the subtree rooted at t.
594 p, q = q, p
595 DEAF.remove_edge(s, t)
596 DEAF.make_root(q)
597 DEAF.add_edge(i, p, q)
598 DEAF.update_potentials(i, p, q)
600 ###########################################################################
601 # Infeasibility and unboundedness detection
602 ###########################################################################
604 if any(DEAF.edge_flow[i] != 0 for i in range(-n, 0)):
605 raise nx.NetworkXUnfeasible("no flow satisfies all node demands")
607 if any(DEAF.edge_flow[i] * 2 >= faux_inf for i in range(DEAF.edge_count)) or any(
608 e[-1].get(capacity, inf) == inf and e[-1].get(weight, 0) < 0
609 for e in nx.selfloop_edges(G, data=True)
610 ):
611 raise nx.NetworkXUnbounded("negative cycle with infinite capacity found")
613 ###########################################################################
614 # Flow cost calculation and flow dict construction
615 ###########################################################################
617 del DEAF.edge_flow[DEAF.edge_count :]
618 flow_cost = sum(w * x for w, x in zip(DEAF.edge_weights, DEAF.edge_flow))
619 flow_dict = {n: {} for n in DEAF.node_list}
621 def add_entry(e):
622 """Add a flow dict entry."""
623 d = flow_dict[e[0]]
624 for k in e[1:-2]:
625 try:
626 d = d[k]
627 except KeyError:
628 t = {}
629 d[k] = t
630 d = t
631 d[e[-2]] = e[-1]
633 DEAF.edge_sources = (
634 DEAF.node_list[s] for s in DEAF.edge_sources
635 ) # Use original nodes.
636 DEAF.edge_targets = (
637 DEAF.node_list[t] for t in DEAF.edge_targets
638 ) # Use original nodes.
639 if not multigraph:
640 for e in zip(DEAF.edge_sources, DEAF.edge_targets, DEAF.edge_flow):
641 add_entry(e)
642 edges = G.edges(data=True)
643 else:
644 for e in zip(
645 DEAF.edge_sources, DEAF.edge_targets, DEAF.edge_keys, DEAF.edge_flow
646 ):
647 add_entry(e)
648 edges = G.edges(data=True, keys=True)
649 for e in edges:
650 if e[0] != e[1]:
651 if e[-1].get(capacity, inf) == 0:
652 add_entry(e[:-1] + (0,))
653 else:
654 w = e[-1].get(weight, 0)
655 if w >= 0:
656 add_entry(e[:-1] + (0,))
657 else:
658 c = e[-1][capacity]
659 flow_cost += w * c
660 add_entry(e[:-1] + (c,))
662 return flow_cost, flow_dict