From f93a3b5a3afb239711fe4542a473b64a2aaf68f8 Mon Sep 17 00:00:00 2001 From: Kohei Morita Date: Mon, 14 Dec 2020 01:01:10 +0900 Subject: [PATCH 1/2] optimize mcf --- atcoder/internal_csr | 1 + atcoder/internal_csr.hpp | 33 +++++++ atcoder/internal_scc.hpp | 20 +---- atcoder/mincostflow.hpp | 186 ++++++++++++++++++++++++--------------- 4 files changed, 153 insertions(+), 87 deletions(-) create mode 100644 atcoder/internal_csr create mode 100644 atcoder/internal_csr.hpp diff --git a/atcoder/internal_csr b/atcoder/internal_csr new file mode 100644 index 0000000..d6a98f5 --- /dev/null +++ b/atcoder/internal_csr @@ -0,0 +1 @@ +#include "atcoder/internal_csr.hpp" diff --git a/atcoder/internal_csr.hpp b/atcoder/internal_csr.hpp new file mode 100644 index 0000000..c102c0b --- /dev/null +++ b/atcoder/internal_csr.hpp @@ -0,0 +1,33 @@ +#ifndef ATCODER_INTERNAL_CSR_HPP +#define ATCODER_INTERNAL_CSR_HPP 1 + +#include +#include +#include + +namespace atcoder { +namespace internal { + +template struct csr { + std::vector start; + std::vector elist; + csr(int n, const std::vector>& edges) + : start(n + 1), elist(edges.size()) { + for (auto e : edges) { + start[e.first + 1]++; + } + for (int i = 1; i <= n; i++) { + start[i] += start[i - 1]; + } + auto counter = start; + for (auto e : edges) { + elist[counter[e.first]++] = e.second; + } + } +}; + +} // namespace internal + +} // namespace atcoder + +#endif // ATCODER_INTERNAL_CSR_HPP diff --git a/atcoder/internal_scc.hpp b/atcoder/internal_scc.hpp index 4fc811e..ab5ed4e 100644 --- a/atcoder/internal_scc.hpp +++ b/atcoder/internal_scc.hpp @@ -5,27 +5,11 @@ #include #include +#include "atcoder/internal_csr" + namespace atcoder { namespace internal { -template struct csr { - std::vector start; - std::vector elist; - csr(int n, const std::vector>& edges) - : start(n + 1), elist(edges.size()) { - for (auto e : edges) { - start[e.first + 1]++; - } - for (int i = 1; i <= n; i++) { - start[i] += start[i - 1]; - } - auto counter = start; - for (auto e : edges) { - elist[counter[e.first]++] = e.second; - } - } -}; - // Reference: // R. Tarjan, // Depth-First Search and Linear Graph Algorithms diff --git a/atcoder/mincostflow.hpp b/atcoder/mincostflow.hpp index d532aac..3808853 100644 --- a/atcoder/mincostflow.hpp +++ b/atcoder/mincostflow.hpp @@ -7,25 +7,23 @@ #include #include +#include "atcoder/internal_csr" +#include "atcoder/internal_queue" + namespace atcoder { template struct mcf_graph { public: mcf_graph() {} - mcf_graph(int n) : _n(n), g(n) {} + mcf_graph(int n) : _n(n) {} int add_edge(int from, int to, Cap cap, Cost cost) { assert(0 <= from && from < _n); assert(0 <= to && to < _n); assert(0 <= cap); assert(0 <= cost); - int m = int(pos.size()); - pos.push_back({from, int(g[from].size())}); - int from_id = int(g[from].size()); - int to_id = int(g[to].size()); - if (from == to) to_id++; - g[from].push_back(_edge{to, to_id, cap, cost}); - g[to].push_back(_edge{from, from_id, 0, -cost}); + int m = int(_edges.size()); + _edges.push_back({from, to, cap, 0, cost}); return m; } @@ -36,28 +34,17 @@ template struct mcf_graph { }; edge get_edge(int i) { - int m = int(pos.size()); + int m = int(_edges.size()); assert(0 <= i && i < m); - auto _e = g[pos[i].first][pos[i].second]; - auto _re = g[_e.to][_e.rev]; - return edge{ - pos[i].first, _e.to, _e.cap + _re.cap, _re.cap, _e.cost, - }; - } - std::vector edges() { - int m = int(pos.size()); - std::vector result(m); - for (int i = 0; i < m; i++) { - result[i] = get_edge(i); - } - return result; + return _edges[i]; } + std::vector edges() { return _edges; } std::pair flow(int s, int t) { return flow(s, t, std::numeric_limits::max()); } std::pair flow(int s, int t, Cap flow_limit) { - return slope(s, t, flow_limit).back(); + return slope(s, t, flow_limit).back(); } std::vector> slope(int s, int t) { return slope(s, t, std::numeric_limits::max()); @@ -66,49 +53,122 @@ template struct mcf_graph { assert(0 <= s && s < _n); assert(0 <= t && t < _n); assert(s != t); + + int m = int(_edges.size()); + std::vector edge_idx(m); + + auto g = [&]() { + std::vector degree(_n), redge_idx(m); + std::vector> elist; + elist.reserve(2 * m); + for (int i = 0; i < m; i++) { + auto e = _edges[i]; + edge_idx[i] = degree[e.from]++; + redge_idx[i] = degree[e.to]++; + elist.push_back({e.from, {e.to, -1, e.cap - e.flow, e.cost}}); + elist.push_back({e.to, {e.from, -1, e.flow, -e.cost}}); + } + auto _g = internal::csr<_edge>(_n, elist); + for (int i = 0; i < m; i++) { + auto e = _edges[i]; + edge_idx[i] += _g.start[e.from]; + redge_idx[i] += _g.start[e.to]; + _g.elist[edge_idx[i]].rev = redge_idx[i]; + _g.elist[redge_idx[i]].rev = edge_idx[i]; + } + return _g; + }(); + + auto result = slope(g, s, t, flow_limit); + + for (int i = 0; i < m; i++) { + auto e = g.elist[edge_idx[i]]; + _edges[i].flow = _edges[i].cap - e.cap; + } + + return result; + } + + private: + int _n; + std::vector _edges; + + // inside edge + struct _edge { + int to, rev; + Cap cap; + Cost cost; + }; + + std::vector> slope(internal::csr<_edge>& g, + int s, + int t, + Cap flow_limit) { // variants (C = maxcost): // -(n-1)C <= dual[s] <= dual[i] <= dual[t] = 0 // reduced cost (= e.cost + dual[e.from] - dual[e.to]) >= 0 for all edge - std::vector dual(_n, 0), dist(_n); - std::vector pv(_n), pe(_n); + + // dual_dist[i] = (dual[i], dist[i]) + std::vector> dual_dist(_n); + std::vector prev_e(_n); std::vector vis(_n); struct Q { Cost key; int to; bool operator<(Q r) const { return key > r.key; } }; + std::vector que_min; std::vector que; auto dual_ref = [&]() { - std::fill(dist.begin(), dist.end(), - std::numeric_limits::max()); + for (int i = 0; i < _n; i++) { + dual_dist[i].second = std::numeric_limits::max(); + } std::fill(vis.begin(), vis.end(), false); + que_min.clear(); que.clear(); - dist[s] = 0; - que.push_back(Q{0, s}); - std::push_heap(que.begin(), que.end()); - while (!que.empty()) { - int v = que.front().to; - std::pop_heap(que.begin(), que.end()); - que.pop_back(); + // que[0..heap_r) was heapified + size_t heap_r = 0; + + dual_dist[s].second = 0; + que_min.push_back(s); + while (!que_min.empty() || !que.empty()) { + int v; + if (!que_min.empty()) { + v = que_min.back(); + que_min.pop_back(); + } else { + while (heap_r < que.size()) { + heap_r++; + std::push_heap(que.begin(), que.begin() + heap_r); + } + v = que.front().to; + std::pop_heap(que.begin(), que.end()); + que.pop_back(); + heap_r--; + } if (vis[v]) continue; vis[v] = true; if (v == t) break; // dist[v] = shortest(s, v) + dual[s] - dual[v] // dist[v] >= 0 (all reduced cost are positive) // dist[v] <= (n-1)C - for (int i = 0; i < int(g[v].size()); i++) { - auto e = g[v][i]; - if (vis[e.to] || !e.cap) continue; + Cost dual_v = dual_dist[v].first, dist_v = dual_dist[v].second; + for (int i = g.start[v]; i < g.start[v + 1]; i++) { + auto e = g.elist[i]; + if (!e.cap) continue; // |-dual[e.to] + dual[v]| <= (n-1)C // cost <= C - -(n-1)C + 0 = nC - Cost cost = e.cost - dual[e.to] + dual[v]; - if (dist[e.to] - dist[v] > cost) { - dist[e.to] = dist[v] + cost; - pv[e.to] = v; - pe[e.to] = i; - que.push_back(Q{dist[e.to], e.to}); - std::push_heap(que.begin(), que.end()); + Cost cost = e.cost - dual_dist[e.to].first + dual_v; + if (dual_dist[e.to].second - dist_v > cost) { + Cost dist_to = dist_v + cost; + dual_dist[e.to].second = dist_to; + prev_e[e.to] = e.rev; + if (dist_to == dist_v) { + que_min.push_back(e.to); + } else { + que.push_back(Q{dist_to, e.to}); + } } } } @@ -119,29 +179,29 @@ template struct mcf_graph { for (int v = 0; v < _n; v++) { if (!vis[v]) continue; // dual[v] = dual[v] - dist[t] + dist[v] - // = dual[v] - (shortest(s, t) + dual[s] - dual[t]) + (shortest(s, v) + dual[s] - dual[v]) - // = - shortest(s, t) + dual[t] + shortest(s, v) - // = shortest(s, v) - shortest(s, t) >= 0 - (n-1)C - dual[v] -= dist[t] - dist[v]; + // = dual[v] - (shortest(s, t) + dual[s] - dual[t]) + + // (shortest(s, v) + dual[s] - dual[v]) = - shortest(s, + // t) + dual[t] + shortest(s, v) = shortest(s, v) - + // shortest(s, t) >= 0 - (n-1)C + dual_dist[v].first -= dual_dist[t].second - dual_dist[v].second; } return true; }; Cap flow = 0; Cost cost = 0, prev_cost_per_flow = -1; - std::vector> result; - result.push_back({flow, cost}); + std::vector> result = {{Cap(0), Cost(0)}}; while (flow < flow_limit) { if (!dual_ref()) break; Cap c = flow_limit - flow; - for (int v = t; v != s; v = pv[v]) { - c = std::min(c, g[pv[v]][pe[v]].cap); + for (int v = t; v != s; v = g.elist[prev_e[v]].to) { + c = std::min(c, g.elist[g.elist[prev_e[v]].rev].cap); } - for (int v = t; v != s; v = pv[v]) { - auto& e = g[pv[v]][pe[v]]; - e.cap -= c; - g[v][e.rev].cap += c; + for (int v = t; v != s; v = g.elist[prev_e[v]].to) { + auto& e = g.elist[prev_e[v]]; + e.cap += c; + g.elist[e.rev].cap -= c; } - Cost d = -dual[s]; + Cost d = -dual_dist[s].first; flow += c; cost += c * d; if (prev_cost_per_flow == d) { @@ -152,18 +212,6 @@ template struct mcf_graph { } return result; } - - private: - int _n; - - struct _edge { - int to, rev; - Cap cap; - Cost cost; - }; - - std::vector> pos; - std::vector> g; }; } // namespace atcoder From 6c48f276d20a2290237e94a073f5ab005f927e9e Mon Sep 17 00:00:00 2001 From: Kohei Morita Date: Mon, 14 Dec 2020 01:12:31 +0900 Subject: [PATCH 2/2] remove space --- atcoder/mincostflow.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/atcoder/mincostflow.hpp b/atcoder/mincostflow.hpp index 3808853..4525a13 100644 --- a/atcoder/mincostflow.hpp +++ b/atcoder/mincostflow.hpp @@ -44,7 +44,7 @@ template struct mcf_graph { return flow(s, t, std::numeric_limits::max()); } std::pair flow(int s, int t, Cap flow_limit) { - return slope(s, t, flow_limit).back(); + return slope(s, t, flow_limit).back(); } std::vector> slope(int s, int t) { return slope(s, t, std::numeric_limits::max());