# min cost max flow blog

a guest
Jul 18th, 2022
1,053
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. #include <iostream>
2. #include <vector>
3. #include <queue>
4. #include <tuple>
5.
6. using namespace std;
7.
8. typedef long long ll;
9.
10. template<typename T>
11. using inverse_priority_queue = priority_queue<T, vector<T>, greater<T>>;
12.
13. const ll INF = 1e18 + 1000;
14.
15. // A reference implementation of Dinitz's algorithm for minimum-cost flows.
16. // Note: this implementation has been created with educational
17. // value in mind, not efficiency. Other implementations with
18. // better constant factors exist; if you're looking for a black
19. // box to use in a contest, you're better off looking elsewhere.
20.
21. // This struct represents the state of an edge, together with its
22. // reverse edge that may appear in some G_f, during the execution
23. // of the algorithm.
24. struct Edge {
25. int from, to; // The edge is from->to in the original graph.
26. ll capac, flow;
27.
28. // The cost of the edge. In this implementation, we add potentials
29. // to the edges directly (it is easy to see that if we add potentials,
30. // the reverse edge will still have cost equal to -(forward edge cost).
31. ll _cost;
32.
33. // Gets the endpoint other than u of the edge.
34. int oth (int u) {
35. return u ^ from ^ to; // xor trick
36. }
37.
38. // Gets the capacity of the edge u->oth(u) in G_f.
39. ll capacity (int u) {
40. if (u == from) {
41. return capac - flow;
42. } else {
43. return flow;
44. }
45. }
46.
47. // Gets the cost of the edge u->oth(u) in G_f
48. ll cost (int u) {
49. if (u == from) {
50. return _cost;
51. } else {
52. return -_cost;
53. }
54. }
55.
56. // Returns whether the edge u->oth(u) exists in u.
57. bool exists (int u) {
58. return capacity(u) != 0;
59. }
60.
61. // Adds df flow to the edge in G if u == from,
62. // subtracts it otherwise.
63. void add_flow (int u, ll df) {
64. if (u == from) {
65. flow += df;
66. } else {
67. flow -= df;
68. }
69. }
70. };
71.
72. // At any given time, represents one particular G_f. We modify this class on the
73. // fly when adding a flow.
74. class MaxFlow {
75. // number of vertices
76. int n;
77. int source, sink;
78.
79. // the current potential of the sink; the only potential which we actually need
80. ll sink_potential;
81.
82. // adj[u] is the set of edges incident to u
84.
85. // we also keep a list of all edges, for easy retrieval of flow on edge
86. vector<Edge*> edges;
87.
88. void recalc_edge_costs (const vector<ll> &dist) {
89. for (auto e : edges) {
90. if (dist[e->from] >= INF || dist[e->to] >= INF) {
91. // one of these vertices is no longer reachable. we can ignore these edges entirely
92. continue;
93. }
94.
95. e->_cost += dist[e->from];
96. e->_cost -= dist[e->to];
97. }
98.
99. sink_potential += dist[sink];
100. }
101.
102. // returns true if sink is reachable from source
103. bool recalc_dist_bellman_ford () {
104. vector<ll> dist (n, INF);
105.
106. dist[source] = 0;
107. for (int k = 0; k < n; k++) {
108. for (auto e : edges) {
109. if (e->from == e->to) {
110. // see the constructor. we pretend those don't exist.
111. continue;
112. }
113.
114. // this is calculated at the beginning, so we can use "from" and "to", no
115. // flow has been added yet.
116. dist[e->to] = min(dist[e->to], dist[e->from] + e->_cost);
117. }
118. }
119.
120. if (dist[sink] >= INF) {
121. return false;
122. }
123.
124. recalc_edge_costs(dist);
125. return true;
126. }
127.
128. bool recalc_dist_dijkstra () {
129. vector<ll> dist (n, INF);
130. vector<int> processed (n, 0);
131.
132. inverse_priority_queue<pair<ll, int>> Q;
133. dist[source] = 0;
134. Q.push({dist[source], source});
135.
136. while (!Q.empty()) {
137. int u;
138. ll d;
139. tie(d, u) = Q.top();
140. Q.pop();
141.
142. if (processed[u]) {
143. continue;
144. }
145. processed[u] = true;
146.
147. for (auto e : adj[u]) {
148. if (!e->exists(u)) {
149. continue;
150. }
151.
152. if (dist[u] + e->cost(u) < dist[e->oth(u)]) {
153. dist[e->oth(u)] = dist[u] + e->cost(u);
154. Q.push({dist[e->oth(u)], e->oth(u)});
155. }
156. }
157. }
158.
159. if (dist[sink] >= INF) {
160. return false;
161. }
162.
163. recalc_edge_costs(dist);
164. return true;
165. }
166.
167. // lvl[u] is the layer (or distance from s) at the current DAG.
168. // updated via dinitz_bfs
169. vector<int> lvl;
170.
171. void dinitz_bfs () {
172. // we use a value greater than n to denote "unexplored"
173. // or what is called "undefined" in the blog. no vertex actually
174. // has a distance greater than n.
175. lvl = vector<int> (n, n + 10);
176.
177. queue<int> Q;
178. Q.push(source);
179. lvl[source] = 0;
180.
181. while (!Q.empty()) {
182. int u = Q.front();
183. Q.pop();
184.
185. for (auto e : adj[u]) {
186. if (!e->exists(u)) {
187. continue;
188. }
189.
190. // this is supposed to be run in the shortest-paths DAG.
191. // if we use shortest-path potentials, then the shortest-paths DAG is
192. // exactly the edges that have cost 0.
193. if (e->cost(u) != 0) {
194. continue;
195. }
196.
197. int v = e->oth(u);
198. if (lvl[v] > n) {
199. lvl[v] = lvl[u] + 1;
200. Q.push(v);
201. }
202.
203. // in the actual implementation, we don't explicitly add edges to the
204. // DAG. we simply say that an edge u->v is in the DAG if lvl[v] = lvl[u] + 1
205. }
206. }
207. }
208.
209. bool in_current_dag (int u, int v) {
210. return lvl[v] == lvl[u] + 1;
211. }
212.
213. // as in the blog, tries to push F flow from u
214. // returns how much flow could actually be pushed
215. // the third parameter is used to maintain what is
216. // called v.blocked in the blog; note that it is a
217. // reference.
218. ll dinitz_dfs (int u, ll F, vector<int> &blocked) {
219. if (u == sink) {
220. return F;
221. }
222.
223. ll flow_pushed = 0;
224. bool all_blocked = true;
225. for (auto e : adj[u]) {
226. int v = e->oth(u);
227. if (!in_current_dag(u, v)) {
228. continue;
229. }
230.
231. if (e->exists(u) && e->cost(u) == 0 && !blocked[v]) {
232. // note that here we use e->capacity(u) instead of
233. // the e.capacity - e.flow in the blog, as the Edge
234. // struct automatically modifies the capacities when
236. ll dF = dinitz_dfs(v, min(F, e->capacity(u)), blocked);
237.
238. flow_pushed += dF;
239. F -= dF;
241. }
242.
243. if (e->exists(u) && e->cost(u) == 0 && !blocked[v]) {
244. all_blocked = false;
245. }
246. }
247.
248. if (all_blocked) {
249. blocked[u] = true;
250. }
251.
252. return flow_pushed;
253. }
254.
255. ll dinitz_dfs () {
256. vector<int> blocked (n, false);
257. return dinitz_dfs(source, INF, blocked);
258. }
259.
260. public:
261. MaxFlow (int _source, int _sink, const vector<tuple<int, int, ll, ll>> &_edges)
262. : source(_source), sink(_sink), sink_potential(0) {
263. n = max(1 + source, 1 + sink);
264. for (auto e : _edges) {
265. n = max(n, max(1 + get<0>(e), 1 + get<1>(e)));
266. }
267.
268. adj = vector<vector<Edge*>> (n, vector<Edge*> (0));
269. for (auto e : _edges) {
270. auto ee = new Edge();
271. ee->from = get<0>(e);
272. ee->to = get<1>(e);
273. ee->flow = 0;
274. ee->capac = get<2>(e);
275. ee->_cost = get<3>(e);
276.
277. edges.push_back(ee);
278. if (ee->from != ee->to) {
279. // don't add self-loops to the adjacency list; they don't do anything
280. // useful to the graph and may mess up the algorithm.
281. // if they are positive, then they are useless
282. // if they are negative, then the assumption is violated
285. }
286. }
287. }
288.
289. // returns a pair <flow, cost>
290. pair<ll, ll> calc_max_flow () {
291. bool any_path = recalc_dist_bellman_ford();
292. if (!any_path) {
293. return {0LL, 0LL};
294. }
295.
296. ll tot_flow = 0, tot_cost = 0;
297. while (true) {
298. any_path = recalc_dist_dijkstra();
299. if (!any_path) {
300. break;
301. }
302.
303. while (true) {
304. dinitz_bfs();
305. if (lvl[sink] > n) {
306. // t is not reachable from s anymore.
307. break;
308. }
309.
310. ll cur_flow = dinitz_dfs();
311. tot_flow += cur_flow;
312. tot_cost += cur_flow * sink_potential;
313. }
314. }
315.
316. return {tot_flow, tot_cost};
317. }
318.
319. ll flow_on_edge (int idx) {
320. return edges[idx]->flow;
321. }
322. };
323.
324. // convenience class for building the graph without awkward make_tuples
325. class GraphBuilder {
326. int source, sink;
327. vector<tuple<int, int, ll, ll>> edges;
328.
329. public:
330. GraphBuilder (int _source, int _sink) : source(_source), sink(_sink), edges(0) {}
331.
332. void add_edge (int u, int v, ll capacity, ll cost) {
333. edges.push_back(make_tuple(u, v, capacity, cost));
334. }
335.
336. MaxFlow build () {
337. return MaxFlow(source, sink, edges);
338. }
339. };
340.
341. int main () {
342. ios::sync_with_stdio(false);
343. cin.tie(0);
344.
345. int n, m;
346. cin >> n >> m;
347.
348. GraphBuilder gb (1, n);
349. for (int i = 0; i < m; i++) {
350. int u, v;
351. ll capac, cost;
352. cin >> u >> v >> capac >> cost;
353.