capacityscaling.py 14 KB

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