Mercurial > repos > shellac > sam_consensus_v3
comparison env/lib/python3.9/site-packages/networkx/algorithms/flow/networksimplex.py @ 0:4f3585e2f14b draft default tip
"planemo upload commit 60cee0fc7c0cda8592644e1aad72851dec82c959"
author | shellac |
---|---|
date | Mon, 22 Mar 2021 18:12:50 +0000 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:4f3585e2f14b |
---|---|
1 """ | |
2 Minimum cost flow algorithms on directed connected graphs. | |
3 """ | |
4 | |
5 __all__ = ["network_simplex"] | |
6 | |
7 from itertools import chain, islice, repeat | |
8 from math import ceil, sqrt | |
9 import networkx as nx | |
10 from networkx.utils import not_implemented_for | |
11 | |
12 | |
13 @not_implemented_for("undirected") | |
14 def network_simplex(G, demand="demand", capacity="capacity", weight="weight"): | |
15 r"""Find a minimum cost flow satisfying all demands in digraph G. | |
16 | |
17 This is a primal network simplex algorithm that uses the leaving | |
18 arc rule to prevent cycling. | |
19 | |
20 G is a digraph with edge costs and capacities and in which nodes | |
21 have demand, i.e., they want to send or receive some amount of | |
22 flow. A negative demand means that the node wants to send flow, a | |
23 positive demand means that the node want to receive flow. A flow on | |
24 the digraph G satisfies all demand if the net flow into each node | |
25 is equal to the demand of that node. | |
26 | |
27 Parameters | |
28 ---------- | |
29 G : NetworkX graph | |
30 DiGraph on which a minimum cost flow satisfying all demands is | |
31 to be found. | |
32 | |
33 demand : string | |
34 Nodes of the graph G are expected to have an attribute demand | |
35 that indicates how much flow a node wants to send (negative | |
36 demand) or receive (positive demand). Note that the sum of the | |
37 demands should be 0 otherwise the problem in not feasible. If | |
38 this attribute is not present, a node is considered to have 0 | |
39 demand. Default value: 'demand'. | |
40 | |
41 capacity : string | |
42 Edges of the graph G are expected to have an attribute capacity | |
43 that indicates how much flow the edge can support. If this | |
44 attribute is not present, the edge is considered to have | |
45 infinite capacity. Default value: 'capacity'. | |
46 | |
47 weight : string | |
48 Edges of the graph G are expected to have an attribute weight | |
49 that indicates the cost incurred by sending one unit of flow on | |
50 that edge. If not present, the weight is considered to be 0. | |
51 Default value: 'weight'. | |
52 | |
53 Returns | |
54 ------- | |
55 flowCost : integer, float | |
56 Cost of a minimum cost flow satisfying all demands. | |
57 | |
58 flowDict : dictionary | |
59 Dictionary of dictionaries keyed by nodes such that | |
60 flowDict[u][v] is the flow edge (u, v). | |
61 | |
62 Raises | |
63 ------ | |
64 NetworkXError | |
65 This exception is raised if the input graph is not directed, | |
66 not connected or is a multigraph. | |
67 | |
68 NetworkXUnfeasible | |
69 This exception is raised in the following situations: | |
70 | |
71 * The sum of the demands is not zero. Then, there is no | |
72 flow satisfying all demands. | |
73 * There is no flow satisfying all demand. | |
74 | |
75 NetworkXUnbounded | |
76 This exception is raised if the digraph G has a cycle of | |
77 negative cost and infinite capacity. Then, the cost of a flow | |
78 satisfying all demands is unbounded below. | |
79 | |
80 Notes | |
81 ----- | |
82 This algorithm is not guaranteed to work if edge weights or demands | |
83 are floating point numbers (overflows and roundoff errors can | |
84 cause problems). As a workaround you can use integer numbers by | |
85 multiplying the relevant edge attributes by a convenient | |
86 constant factor (eg 100). | |
87 | |
88 See also | |
89 -------- | |
90 cost_of_flow, max_flow_min_cost, min_cost_flow, min_cost_flow_cost | |
91 | |
92 Examples | |
93 -------- | |
94 A simple example of a min cost flow problem. | |
95 | |
96 >>> G = nx.DiGraph() | |
97 >>> G.add_node("a", demand=-5) | |
98 >>> G.add_node("d", demand=5) | |
99 >>> G.add_edge("a", "b", weight=3, capacity=4) | |
100 >>> G.add_edge("a", "c", weight=6, capacity=10) | |
101 >>> G.add_edge("b", "d", weight=1, capacity=9) | |
102 >>> G.add_edge("c", "d", weight=2, capacity=5) | |
103 >>> flowCost, flowDict = nx.network_simplex(G) | |
104 >>> flowCost | |
105 24 | |
106 >>> flowDict # doctest: +SKIP | |
107 {'a': {'c': 1, 'b': 4}, 'c': {'d': 1}, 'b': {'d': 4}, 'd': {}} | |
108 | |
109 The mincost flow algorithm can also be used to solve shortest path | |
110 problems. To find the shortest path between two nodes u and v, | |
111 give all edges an infinite capacity, give node u a demand of -1 and | |
112 node v a demand a 1. Then run the network simplex. The value of a | |
113 min cost flow will be the distance between u and v and edges | |
114 carrying positive flow will indicate the path. | |
115 | |
116 >>> G = nx.DiGraph() | |
117 >>> G.add_weighted_edges_from( | |
118 ... [ | |
119 ... ("s", "u", 10), | |
120 ... ("s", "x", 5), | |
121 ... ("u", "v", 1), | |
122 ... ("u", "x", 2), | |
123 ... ("v", "y", 1), | |
124 ... ("x", "u", 3), | |
125 ... ("x", "v", 5), | |
126 ... ("x", "y", 2), | |
127 ... ("y", "s", 7), | |
128 ... ("y", "v", 6), | |
129 ... ] | |
130 ... ) | |
131 >>> G.add_node("s", demand=-1) | |
132 >>> G.add_node("v", demand=1) | |
133 >>> flowCost, flowDict = nx.network_simplex(G) | |
134 >>> flowCost == nx.shortest_path_length(G, "s", "v", weight="weight") | |
135 True | |
136 >>> sorted([(u, v) for u in flowDict for v in flowDict[u] if flowDict[u][v] > 0]) | |
137 [('s', 'x'), ('u', 'v'), ('x', 'u')] | |
138 >>> nx.shortest_path(G, "s", "v", weight="weight") | |
139 ['s', 'x', 'u', 'v'] | |
140 | |
141 It is possible to change the name of the attributes used for the | |
142 algorithm. | |
143 | |
144 >>> G = nx.DiGraph() | |
145 >>> G.add_node("p", spam=-4) | |
146 >>> G.add_node("q", spam=2) | |
147 >>> G.add_node("a", spam=-2) | |
148 >>> G.add_node("d", spam=-1) | |
149 >>> G.add_node("t", spam=2) | |
150 >>> G.add_node("w", spam=3) | |
151 >>> G.add_edge("p", "q", cost=7, vacancies=5) | |
152 >>> G.add_edge("p", "a", cost=1, vacancies=4) | |
153 >>> G.add_edge("q", "d", cost=2, vacancies=3) | |
154 >>> G.add_edge("t", "q", cost=1, vacancies=2) | |
155 >>> G.add_edge("a", "t", cost=2, vacancies=4) | |
156 >>> G.add_edge("d", "w", cost=3, vacancies=4) | |
157 >>> G.add_edge("t", "w", cost=4, vacancies=1) | |
158 >>> flowCost, flowDict = nx.network_simplex( | |
159 ... G, demand="spam", capacity="vacancies", weight="cost" | |
160 ... ) | |
161 >>> flowCost | |
162 37 | |
163 >>> flowDict # doctest: +SKIP | |
164 {'a': {'t': 4}, 'd': {'w': 2}, 'q': {'d': 1}, 'p': {'q': 2, 'a': 2}, 't': {'q': 1, 'w': 1}, 'w': {}} | |
165 | |
166 References | |
167 ---------- | |
168 .. [1] Z. Kiraly, P. Kovacs. | |
169 Efficient implementation of minimum-cost flow algorithms. | |
170 Acta Universitatis Sapientiae, Informatica 4(1):67--118. 2012. | |
171 .. [2] R. Barr, F. Glover, D. Klingman. | |
172 Enhancement of spanning tree labeling procedures for network | |
173 optimization. | |
174 INFOR 17(1):16--34. 1979. | |
175 """ | |
176 ########################################################################### | |
177 # Problem essentials extraction and sanity check | |
178 ########################################################################### | |
179 | |
180 if len(G) == 0: | |
181 raise nx.NetworkXError("graph has no nodes") | |
182 | |
183 # Number all nodes and edges and hereafter reference them using ONLY their | |
184 # numbers | |
185 | |
186 N = list(G) # nodes | |
187 I = {u: i for i, u in enumerate(N)} # node indices | |
188 D = [G.nodes[u].get(demand, 0) for u in N] # node demands | |
189 | |
190 inf = float("inf") | |
191 for p, b in zip(N, D): | |
192 if abs(b) == inf: | |
193 raise nx.NetworkXError(f"node {p!r} has infinite demand") | |
194 | |
195 multigraph = G.is_multigraph() | |
196 S = [] # edge sources | |
197 T = [] # edge targets | |
198 if multigraph: | |
199 K = [] # edge keys | |
200 E = {} # edge indices | |
201 U = [] # edge capacities | |
202 C = [] # edge weights | |
203 | |
204 if not multigraph: | |
205 edges = G.edges(data=True) | |
206 else: | |
207 edges = G.edges(data=True, keys=True) | |
208 edges = (e for e in edges if e[0] != e[1] and e[-1].get(capacity, inf) != 0) | |
209 for i, e in enumerate(edges): | |
210 S.append(I[e[0]]) | |
211 T.append(I[e[1]]) | |
212 if multigraph: | |
213 K.append(e[2]) | |
214 E[e[:-1]] = i | |
215 U.append(e[-1].get(capacity, inf)) | |
216 C.append(e[-1].get(weight, 0)) | |
217 | |
218 for e, c in zip(E, C): | |
219 if abs(c) == inf: | |
220 raise nx.NetworkXError(f"edge {e!r} has infinite weight") | |
221 if not multigraph: | |
222 edges = nx.selfloop_edges(G, data=True) | |
223 else: | |
224 edges = nx.selfloop_edges(G, data=True, keys=True) | |
225 for e in edges: | |
226 if abs(e[-1].get(weight, 0)) == inf: | |
227 raise nx.NetworkXError(f"edge {e[:-1]!r} has infinite weight") | |
228 | |
229 ########################################################################### | |
230 # Quick infeasibility detection | |
231 ########################################################################### | |
232 | |
233 if sum(D) != 0: | |
234 raise nx.NetworkXUnfeasible("total node demand is not zero") | |
235 for e, u in zip(E, U): | |
236 if u < 0: | |
237 raise nx.NetworkXUnfeasible(f"edge {e!r} has negative capacity") | |
238 if not multigraph: | |
239 edges = nx.selfloop_edges(G, data=True) | |
240 else: | |
241 edges = nx.selfloop_edges(G, data=True, keys=True) | |
242 for e in edges: | |
243 if e[-1].get(capacity, inf) < 0: | |
244 raise nx.NetworkXUnfeasible(f"edge {e[:-1]!r} has negative capacity") | |
245 | |
246 ########################################################################### | |
247 # Initialization | |
248 ########################################################################### | |
249 | |
250 # Add a dummy node -1 and connect all existing nodes to it with infinite- | |
251 # capacity dummy edges. Node -1 will serve as the root of the | |
252 # spanning tree of the network simplex method. The new edges will used to | |
253 # trivially satisfy the node demands and create an initial strongly | |
254 # feasible spanning tree. | |
255 n = len(N) # number of nodes | |
256 for p, d in enumerate(D): | |
257 # Must be greater-than here. Zero-demand nodes must have | |
258 # edges pointing towards the root to ensure strong | |
259 # feasibility. | |
260 if d > 0: | |
261 S.append(-1) | |
262 T.append(p) | |
263 else: | |
264 S.append(p) | |
265 T.append(-1) | |
266 faux_inf = ( | |
267 3 | |
268 * max( | |
269 chain( | |
270 [sum(u for u in U if u < inf), sum(abs(c) for c in C)], | |
271 (abs(d) for d in D), | |
272 ) | |
273 ) | |
274 or 1 | |
275 ) | |
276 C.extend(repeat(faux_inf, n)) | |
277 U.extend(repeat(faux_inf, n)) | |
278 | |
279 # Construct the initial spanning tree. | |
280 e = len(E) # number of edges | |
281 x = list(chain(repeat(0, e), (abs(d) for d in D))) # edge flows | |
282 pi = [faux_inf if d <= 0 else -faux_inf for d in D] # node potentials | |
283 parent = list(chain(repeat(-1, n), [None])) # parent nodes | |
284 edge = list(range(e, e + n)) # edges to parents | |
285 size = list(chain(repeat(1, n), [n + 1])) # subtree sizes | |
286 next = list(chain(range(1, n), [-1, 0])) # next nodes in depth-first thread | |
287 prev = list(range(-1, n)) # previous nodes in depth-first thread | |
288 last = list(chain(range(n), [n - 1])) # last descendants in depth-first thread | |
289 | |
290 ########################################################################### | |
291 # Pivot loop | |
292 ########################################################################### | |
293 | |
294 def reduced_cost(i): | |
295 """Returns the reduced cost of an edge i. | |
296 """ | |
297 c = C[i] - pi[S[i]] + pi[T[i]] | |
298 return c if x[i] == 0 else -c | |
299 | |
300 def find_entering_edges(): | |
301 """Yield entering edges until none can be found. | |
302 """ | |
303 if e == 0: | |
304 return | |
305 | |
306 # Entering edges are found by combining Dantzig's rule and Bland's | |
307 # rule. The edges are cyclically grouped into blocks of size B. Within | |
308 # each block, Dantzig's rule is applied to find an entering edge. The | |
309 # blocks to search is determined following Bland's rule. | |
310 B = int(ceil(sqrt(e))) # pivot block size | |
311 M = (e + B - 1) // B # number of blocks needed to cover all edges | |
312 m = 0 # number of consecutive blocks without eligible | |
313 # entering edges | |
314 f = 0 # first edge in block | |
315 while m < M: | |
316 # Determine the next block of edges. | |
317 l = f + B | |
318 if l <= e: | |
319 edges = range(f, l) | |
320 else: | |
321 l -= e | |
322 edges = chain(range(f, e), range(l)) | |
323 f = l | |
324 # Find the first edge with the lowest reduced cost. | |
325 i = min(edges, key=reduced_cost) | |
326 c = reduced_cost(i) | |
327 if c >= 0: | |
328 # No entering edge found in the current block. | |
329 m += 1 | |
330 else: | |
331 # Entering edge found. | |
332 if x[i] == 0: | |
333 p = S[i] | |
334 q = T[i] | |
335 else: | |
336 p = T[i] | |
337 q = S[i] | |
338 yield i, p, q | |
339 m = 0 | |
340 # All edges have nonnegative reduced costs. The current flow is | |
341 # optimal. | |
342 | |
343 def find_apex(p, q): | |
344 """Find the lowest common ancestor of nodes p and q in the spanning | |
345 tree. | |
346 """ | |
347 size_p = size[p] | |
348 size_q = size[q] | |
349 while True: | |
350 while size_p < size_q: | |
351 p = parent[p] | |
352 size_p = size[p] | |
353 while size_p > size_q: | |
354 q = parent[q] | |
355 size_q = size[q] | |
356 if size_p == size_q: | |
357 if p != q: | |
358 p = parent[p] | |
359 size_p = size[p] | |
360 q = parent[q] | |
361 size_q = size[q] | |
362 else: | |
363 return p | |
364 | |
365 def trace_path(p, w): | |
366 """Returns the nodes and edges on the path from node p to its ancestor | |
367 w. | |
368 """ | |
369 Wn = [p] | |
370 We = [] | |
371 while p != w: | |
372 We.append(edge[p]) | |
373 p = parent[p] | |
374 Wn.append(p) | |
375 return Wn, We | |
376 | |
377 def find_cycle(i, p, q): | |
378 """Returns the nodes and edges on the cycle containing edge i == (p, q) | |
379 when the latter is added to the spanning tree. | |
380 | |
381 The cycle is oriented in the direction from p to q. | |
382 """ | |
383 w = find_apex(p, q) | |
384 Wn, We = trace_path(p, w) | |
385 Wn.reverse() | |
386 We.reverse() | |
387 if We != [i]: | |
388 We.append(i) | |
389 WnR, WeR = trace_path(q, w) | |
390 del WnR[-1] | |
391 Wn += WnR | |
392 We += WeR | |
393 return Wn, We | |
394 | |
395 def residual_capacity(i, p): | |
396 """Returns the residual capacity of an edge i in the direction away | |
397 from its endpoint p. | |
398 """ | |
399 return U[i] - x[i] if S[i] == p else x[i] | |
400 | |
401 def find_leaving_edge(Wn, We): | |
402 """Returns the leaving edge in a cycle represented by Wn and We. | |
403 """ | |
404 j, s = min( | |
405 zip(reversed(We), reversed(Wn)), key=lambda i_p: residual_capacity(*i_p) | |
406 ) | |
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 # Do nothing more if the entering edge is the same as the leaving edge. | |
529 if i != j: | |
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 any( | |
549 e[-1].get(capacity, inf) == inf and e[-1].get(weight, 0) < 0 | |
550 for e in nx.selfloop_edges(G, data=True) | |
551 ): | |
552 raise nx.NetworkXUnbounded("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 |