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"""
2Capacity scaling minimum cost flow algorithm.
3"""
5__all__ = ["capacity_scaling"]
7from itertools import chain
8from math import log
10import networkx as nx
12from ...utils import BinaryHeap, arbitrary_element, not_implemented_for
15def _detect_unboundedness(R):
16 """Detect infinite-capacity negative cycles."""
17 G = nx.DiGraph()
18 G.add_nodes_from(R)
20 # Value simulating infinity.
21 inf = R.graph["inf"]
22 # True infinity.
23 f_inf = float("inf")
24 for u in R:
25 for v, e in R[u].items():
26 # Compute the minimum weight of infinite-capacity (u, v) edges.
27 w = f_inf
28 for k, e in e.items():
29 if e["capacity"] == inf:
30 w = min(w, e["weight"])
31 if w != f_inf:
32 G.add_edge(u, v, weight=w)
34 if nx.negative_edge_cycle(G):
35 raise nx.NetworkXUnbounded(
36 "Negative cost cycle of infinite capacity found. "
37 "Min cost flow may be unbounded below."
38 )
41@not_implemented_for("undirected")
42def _build_residual_network(G, demand, capacity, weight):
43 """Build a residual network and initialize a zero flow."""
44 if sum(G.nodes[u].get(demand, 0) for u in G) != 0:
45 raise nx.NetworkXUnfeasible("Sum of the demands should be 0.")
47 R = nx.MultiDiGraph()
48 R.add_nodes_from(
49 (u, {"excess": -G.nodes[u].get(demand, 0), "potential": 0}) for u in G
50 )
52 inf = float("inf")
53 # Detect selfloops with infinite capacities and negative weights.
54 for u, v, e in nx.selfloop_edges(G, data=True):
55 if e.get(weight, 0) < 0 and e.get(capacity, inf) == inf:
56 raise nx.NetworkXUnbounded(
57 "Negative cost cycle of infinite capacity found. "
58 "Min cost flow may be unbounded below."
59 )
61 # Extract edges with positive capacities. Self loops excluded.
62 if G.is_multigraph():
63 edge_list = [
64 (u, v, k, e)
65 for u, v, k, e in G.edges(data=True, keys=True)
66 if u != v and e.get(capacity, inf) > 0
67 ]
68 else:
69 edge_list = [
70 (u, v, 0, e)
71 for u, v, e in G.edges(data=True)
72 if u != v and e.get(capacity, inf) > 0
73 ]
74 # Simulate infinity with the larger of the sum of absolute node imbalances
75 # the sum of finite edge capacities or any positive value if both sums are
76 # zero. This allows the infinite-capacity edges to be distinguished for
77 # unboundedness detection and directly participate in residual capacity
78 # calculation.
79 inf = (
80 max(
81 sum(abs(R.nodes[u]["excess"]) for u in R),
82 2
83 * sum(
84 e[capacity]
85 for u, v, k, e in edge_list
86 if capacity in e and e[capacity] != inf
87 ),
88 )
89 or 1
90 )
91 for u, v, k, e in edge_list:
92 r = min(e.get(capacity, inf), inf)
93 w = e.get(weight, 0)
94 # Add both (u, v) and (v, u) into the residual network marked with the
95 # original key. (key[1] == True) indicates the (u, v) is in the
96 # original network.
97 R.add_edge(u, v, key=(k, True), capacity=r, weight=w, flow=0)
98 R.add_edge(v, u, key=(k, False), capacity=0, weight=-w, flow=0)
100 # Record the value simulating infinity.
101 R.graph["inf"] = inf
103 _detect_unboundedness(R)
105 return R
108def _build_flow_dict(G, R, capacity, weight):
109 """Build a flow dictionary from a residual network."""
110 inf = float("inf")
111 flow_dict = {}
112 if G.is_multigraph():
113 for u in G:
114 flow_dict[u] = {}
115 for v, es in G[u].items():
116 flow_dict[u][v] = {
117 # Always saturate negative selfloops.
118 k: (
119 0
120 if (
121 u != v or e.get(capacity, inf) <= 0 or e.get(weight, 0) >= 0
122 )
123 else e[capacity]
124 )
125 for k, e in es.items()
126 }
127 for v, es in R[u].items():
128 if v in flow_dict[u]:
129 flow_dict[u][v].update(
130 (k[0], e["flow"]) for k, e in es.items() if e["flow"] > 0
131 )
132 else:
133 for u in G:
134 flow_dict[u] = {
135 # Always saturate negative selfloops.
136 v: (
137 0
138 if (u != v or e.get(capacity, inf) <= 0 or e.get(weight, 0) >= 0)
139 else e[capacity]
140 )
141 for v, e in G[u].items()
142 }
143 flow_dict[u].update(
144 (v, e["flow"])
145 for v, es in R[u].items()
146 for e in es.values()
147 if e["flow"] > 0
148 )
149 return flow_dict
152@nx._dispatchable(
153 node_attrs="demand", edge_attrs={"capacity": float("inf"), "weight": 0}
154)
155def capacity_scaling(
156 G, demand="demand", capacity="capacity", weight="weight", heap=BinaryHeap
157):
158 r"""Find a minimum cost flow satisfying all demands in digraph G.
160 This is a capacity scaling successive shortest augmenting path algorithm.
162 G is a digraph with edge costs and capacities and in which nodes
163 have demand, i.e., they want to send or receive some amount of
164 flow. A negative demand means that the node wants to send flow, a
165 positive demand means that the node want to receive flow. A flow on
166 the digraph G satisfies all demand if the net flow into each node
167 is equal to the demand of that node.
169 Parameters
170 ----------
171 G : NetworkX graph
172 DiGraph or MultiDiGraph on which a minimum cost flow satisfying all
173 demands is to be found.
175 demand : string
176 Nodes of the graph G are expected to have an attribute demand
177 that indicates how much flow a node wants to send (negative
178 demand) or receive (positive demand). Note that the sum of the
179 demands should be 0 otherwise the problem in not feasible. If
180 this attribute is not present, a node is considered to have 0
181 demand. Default value: 'demand'.
183 capacity : string
184 Edges of the graph G are expected to have an attribute capacity
185 that indicates how much flow the edge can support. If this
186 attribute is not present, the edge is considered to have
187 infinite capacity. Default value: 'capacity'.
189 weight : string
190 Edges of the graph G are expected to have an attribute weight
191 that indicates the cost incurred by sending one unit of flow on
192 that edge. If not present, the weight is considered to be 0.
193 Default value: 'weight'.
195 heap : class
196 Type of heap to be used in the algorithm. It should be a subclass of
197 :class:`MinHeap` or implement a compatible interface.
199 If a stock heap implementation is to be used, :class:`BinaryHeap` is
200 recommended over :class:`PairingHeap` for Python implementations without
201 optimized attribute accesses (e.g., CPython) despite a slower
202 asymptotic running time. For Python implementations with optimized
203 attribute accesses (e.g., PyPy), :class:`PairingHeap` provides better
204 performance. Default value: :class:`BinaryHeap`.
206 Returns
207 -------
208 flowCost : integer
209 Cost of a minimum cost flow satisfying all demands.
211 flowDict : dictionary
212 If G is a digraph, a dict-of-dicts keyed by nodes such that
213 flowDict[u][v] is the flow on edge (u, v).
214 If G is a MultiDiGraph, a dict-of-dicts-of-dicts keyed by nodes
215 so that flowDict[u][v][key] is the flow on edge (u, v, key).
217 Raises
218 ------
219 NetworkXError
220 This exception is raised if the input graph is not directed,
221 not connected.
223 NetworkXUnfeasible
224 This exception is raised in the following situations:
226 * The sum of the demands is not zero. Then, there is no
227 flow satisfying all demands.
228 * There is no flow satisfying all demand.
230 NetworkXUnbounded
231 This exception is raised if the digraph G has a cycle of
232 negative cost and infinite capacity. Then, the cost of a flow
233 satisfying all demands is unbounded below.
235 Notes
236 -----
237 This algorithm does not work if edge weights are floating-point numbers.
239 See also
240 --------
241 :meth:`network_simplex`
243 Examples
244 --------
245 A simple example of a min cost flow problem.
247 >>> G = nx.DiGraph()
248 >>> G.add_node("a", demand=-5)
249 >>> G.add_node("d", demand=5)
250 >>> G.add_edge("a", "b", weight=3, capacity=4)
251 >>> G.add_edge("a", "c", weight=6, capacity=10)
252 >>> G.add_edge("b", "d", weight=1, capacity=9)
253 >>> G.add_edge("c", "d", weight=2, capacity=5)
254 >>> flowCost, flowDict = nx.capacity_scaling(G)
255 >>> flowCost
256 24
257 >>> flowDict
258 {'a': {'b': 4, 'c': 1}, 'd': {}, 'b': {'d': 4}, 'c': {'d': 1}}
260 It is possible to change the name of the attributes used for the
261 algorithm.
263 >>> G = nx.DiGraph()
264 >>> G.add_node("p", spam=-4)
265 >>> G.add_node("q", spam=2)
266 >>> G.add_node("a", spam=-2)
267 >>> G.add_node("d", spam=-1)
268 >>> G.add_node("t", spam=2)
269 >>> G.add_node("w", spam=3)
270 >>> G.add_edge("p", "q", cost=7, vacancies=5)
271 >>> G.add_edge("p", "a", cost=1, vacancies=4)
272 >>> G.add_edge("q", "d", cost=2, vacancies=3)
273 >>> G.add_edge("t", "q", cost=1, vacancies=2)
274 >>> G.add_edge("a", "t", cost=2, vacancies=4)
275 >>> G.add_edge("d", "w", cost=3, vacancies=4)
276 >>> G.add_edge("t", "w", cost=4, vacancies=1)
277 >>> flowCost, flowDict = nx.capacity_scaling(
278 ... G, demand="spam", capacity="vacancies", weight="cost"
279 ... )
280 >>> flowCost
281 37
282 >>> flowDict
283 {'p': {'q': 2, 'a': 2}, 'q': {'d': 1}, 'a': {'t': 4}, 'd': {'w': 2}, 't': {'q': 1, 'w': 1}, 'w': {}}
284 """
285 R = _build_residual_network(G, demand, capacity, weight)
287 inf = float("inf")
288 # Account cost of negative selfloops.
289 flow_cost = sum(
290 0
291 if e.get(capacity, inf) <= 0 or e.get(weight, 0) >= 0
292 else e[capacity] * e[weight]
293 for u, v, e in nx.selfloop_edges(G, data=True)
294 )
296 # Determine the maximum edge capacity.
297 wmax = max(chain([-inf], (e["capacity"] for u, v, e in R.edges(data=True))))
298 if wmax == -inf:
299 # Residual network has no edges.
300 return flow_cost, _build_flow_dict(G, R, capacity, weight)
302 R_nodes = R.nodes
303 R_succ = R.succ
305 delta = 2 ** int(log(wmax, 2))
306 while delta >= 1:
307 # Saturate Δ-residual edges with negative reduced costs to achieve
308 # Δ-optimality.
309 for u in R:
310 p_u = R_nodes[u]["potential"]
311 for v, es in R_succ[u].items():
312 for k, e in es.items():
313 flow = e["capacity"] - e["flow"]
314 if e["weight"] - p_u + R_nodes[v]["potential"] < 0:
315 flow = e["capacity"] - e["flow"]
316 if flow >= delta:
317 e["flow"] += flow
318 R_succ[v][u][(k[0], not k[1])]["flow"] -= flow
319 R_nodes[u]["excess"] -= flow
320 R_nodes[v]["excess"] += flow
321 # Determine the Δ-active nodes.
322 S = set()
323 T = set()
324 S_add = S.add
325 S_remove = S.remove
326 T_add = T.add
327 T_remove = T.remove
328 for u in R:
329 excess = R_nodes[u]["excess"]
330 if excess >= delta:
331 S_add(u)
332 elif excess <= -delta:
333 T_add(u)
334 # Repeatedly augment flow from S to T along shortest paths until
335 # Δ-feasibility is achieved.
336 while S and T:
337 s = arbitrary_element(S)
338 t = None
339 # Search for a shortest path in terms of reduce costs from s to
340 # any t in T in the Δ-residual network.
341 d = {}
342 pred = {s: None}
343 h = heap()
344 h_insert = h.insert
345 h_get = h.get
346 h_insert(s, 0)
347 while h:
348 u, d_u = h.pop()
349 d[u] = d_u
350 if u in T:
351 # Path found.
352 t = u
353 break
354 p_u = R_nodes[u]["potential"]
355 for v, es in R_succ[u].items():
356 if v in d:
357 continue
358 wmin = inf
359 # Find the minimum-weighted (u, v) Δ-residual edge.
360 for k, e in es.items():
361 if e["capacity"] - e["flow"] >= delta:
362 w = e["weight"]
363 if w < wmin:
364 wmin = w
365 kmin = k
366 emin = e
367 if wmin == inf:
368 continue
369 # Update the distance label of v.
370 d_v = d_u + wmin - p_u + R_nodes[v]["potential"]
371 if h_insert(v, d_v):
372 pred[v] = (u, kmin, emin)
373 if t is not None:
374 # Augment Δ units of flow from s to t.
375 while u != s:
376 v = u
377 u, k, e = pred[v]
378 e["flow"] += delta
379 R_succ[v][u][(k[0], not k[1])]["flow"] -= delta
380 # Account node excess and deficit.
381 R_nodes[s]["excess"] -= delta
382 R_nodes[t]["excess"] += delta
383 if R_nodes[s]["excess"] < delta:
384 S_remove(s)
385 if R_nodes[t]["excess"] > -delta:
386 T_remove(t)
387 # Update node potentials.
388 d_t = d[t]
389 for u, d_u in d.items():
390 R_nodes[u]["potential"] -= d_u - d_t
391 else:
392 # Path not found.
393 S_remove(s)
394 delta //= 2
396 if any(R.nodes[u]["excess"] != 0 for u in R):
397 raise nx.NetworkXUnfeasible("No flow satisfying all demands.")
399 # Calculate the flow cost.
400 for u in R:
401 for v, es in R_succ[u].items():
402 for e in es.values():
403 flow = e["flow"]
404 if flow > 0:
405 flow_cost += flow * e["weight"]
407 return flow_cost, _build_flow_dict(G, R, capacity, weight)