Skip to content

optimize mcf #82

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Dec 14, 2020
Merged

optimize mcf #82

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
f93a3b5
optimize mcf
yosupo06 Dec 13, 2020
6c48f27
remove space
yosupo06 Dec 13, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions atcoder/internal_csr
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
#include "atcoder/internal_csr.hpp"
33 changes: 33 additions & 0 deletions atcoder/internal_csr.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,33 @@
#ifndef ATCODER_INTERNAL_CSR_HPP
#define ATCODER_INTERNAL_CSR_HPP 1

#include <algorithm>
#include <utility>
#include <vector>

namespace atcoder {
namespace internal {

template <class E> struct csr {
std::vector<int> start;
std::vector<E> elist;
csr(int n, const std::vector<std::pair<int, E>>& 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
20 changes: 2 additions & 18 deletions atcoder/internal_scc.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,27 +5,11 @@
#include <utility>
#include <vector>

#include "atcoder/internal_csr"

namespace atcoder {
namespace internal {

template <class E> struct csr {
std::vector<int> start;
std::vector<E> elist;
csr(int n, const std::vector<std::pair<int, E>>& 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
Expand Down
184 changes: 116 additions & 68 deletions atcoder/mincostflow.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,25 +7,23 @@
#include <queue>
#include <vector>

#include "atcoder/internal_csr"
#include "atcoder/internal_queue"

namespace atcoder {

template <class Cap, class Cost> 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;
}

Expand All @@ -36,22 +34,11 @@ template <class Cap, class Cost> 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<edge> edges() {
int m = int(pos.size());
std::vector<edge> result(m);
for (int i = 0; i < m; i++) {
result[i] = get_edge(i);
}
return result;
return _edges[i];
}
std::vector<edge> edges() { return _edges; }

std::pair<Cap, Cost> flow(int s, int t) {
return flow(s, t, std::numeric_limits<Cap>::max());
Expand All @@ -66,49 +53,122 @@ template <class Cap, class Cost> struct mcf_graph {
assert(0 <= s && s < _n);
assert(0 <= t && t < _n);
assert(s != t);

int m = int(_edges.size());
std::vector<int> edge_idx(m);

auto g = [&]() {
std::vector<int> degree(_n), redge_idx(m);
std::vector<std::pair<int, _edge>> 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<edge> _edges;

// inside edge
struct _edge {
int to, rev;
Cap cap;
Cost cost;
};

std::vector<std::pair<Cap, Cost>> 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<Cost> dual(_n, 0), dist(_n);
std::vector<int> pv(_n), pe(_n);

// dual_dist[i] = (dual[i], dist[i])
std::vector<std::pair<Cost, Cost>> dual_dist(_n);
std::vector<int> prev_e(_n);
std::vector<bool> vis(_n);
struct Q {
Cost key;
int to;
bool operator<(Q r) const { return key > r.key; }
};
std::vector<int> que_min;
std::vector<Q> que;
auto dual_ref = [&]() {
std::fill(dist.begin(), dist.end(),
std::numeric_limits<Cost>::max());
for (int i = 0; i < _n; i++) {
dual_dist[i].second = std::numeric_limits<Cost>::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});
}
}
}
}
Expand All @@ -119,29 +179,29 @@ template <class Cap, class Cost> 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<std::pair<Cap, Cost>> result;
result.push_back({flow, cost});
std::vector<std::pair<Cap, Cost>> 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) {
Expand All @@ -152,18 +212,6 @@ template <class Cap, class Cost> struct mcf_graph {
}
return result;
}

private:
int _n;

struct _edge {
int to, rev;
Cap cap;
Cost cost;
};

std::vector<std::pair<int, int>> pos;
std::vector<std::vector<_edge>> g;
};

} // namespace atcoder
Expand Down