Advertisement
Guest User

min cost max flow blog

a guest
Jul 18th, 2022
865
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 8.96 KB | None | 0 0
  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
  83. vector<vector<Edge*>> adj;
  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
  235. // we add flow.
  236. ll dF = dinitz_dfs(v, min(F, e->capacity(u)), blocked);
  237.  
  238. flow_pushed += dF;
  239. F -= dF;
  240. e->add_flow(u, 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
  283. adj[ee->from].push_back(ee);
  284. adj[ee->to].push_back(ee);
  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.  
  354. gb.add_edge(u, v, capac, cost);
  355. }
  356.  
  357. auto mf = gb.build();
  358. auto ans = mf.calc_max_flow();
  359. cout << ans.first << " " << ans.second << '\n';
  360. }
  361.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement