面白かった!! 上手に変数変換することで「2 変数劣モジュラ関数の和の最小化」になるタイプの問題だった。
問題概要
考えたこと
一目見て、2 変数劣モジュラ関数の最小化 (燃やす埋める) っぽいと感じた。値が 500 以下というのも怪しい。
次のような変数 を考えたくなった。ここで、 値をとる整数変数を 個の 0-1 変数で表現するテクを用いている。
そして、次のような制約条件が成り立つ。
しかし、上はともかく、下は劣モジュラ関数にはならない。そこで、変数 の定義をビット反転して、以下のように変更することにした。
これで、制約条件は次のように記述できることとなった。
総合すると、各制約条件は次のように書ける (これは実装中のメモ)。
コード
#include <bits/stdc++.h>
using namespace std;
template<class COST> struct TwoVariableSubmodularOpt {
TwoVariableSubmodularOpt() : N(2), S(0), T(0), OFFSET(0) {}
TwoVariableSubmodularOpt(int n, COST inf = 0)
: N(n), S(n), T(n + 1), OFFSET(0), INF(inf), list(n + 2) {}
void init(int n, COST inf = 0) {
N = n, S = n, T = n + 1;
OFFSET = 0, INF = inf;
list.assign(N + 2, Edge());
pos.clear();
}
void add_single_cost(int xi, COST false_cost, COST true_cost) {
assert(0 <= xi && xi < N);
if (false_cost >= true_cost) {
OFFSET += true_cost;
add_edge(S, xi, false_cost - true_cost);
} else {
OFFSET += false_cost;
add_edge(xi, T, true_cost - false_cost);
}
}
void add_psp_constraint(int xi, int xj) {
assert(0 <= xi && xi < N);
assert(0 <= xj && xj < N);
add_edge(xi, xj, INF);
}
void add_psp_penalty(int xi, int xj, COST C) {
assert(0 <= xi && xi < N);
assert(0 <= xj && xj < N);
assert(C >= 0);
add_edge(xi, xj, C);
}
void add_both_true_profit(int xi, int xj, COST P) {
assert(0 <= xi && xi < N);
assert(0 <= xj && xj < N);
assert(P >= 0);
OFFSET -= P;
add_edge(S, xi, P);
add_edge(xi, xj, P);
}
void add_both_false_profit(int xi, int xj, COST P) {
assert(0 <= xi && xi < N);
assert(0 <= xj && xj < N);
assert(P >= 0);
OFFSET -= P;
add_edge(xj, T, P);
add_edge(xi, xj, P);
}
void add_submodular_function(int xi, int xj, COST A, COST B, COST C, COST D) {
assert(0 <= xi && xi < N);
assert(0 <= xj && xj < N);
assert(B + C >= A + D);
OFFSET += A;
add_single_cost(xi, 0, D - B);
add_single_cost(xj, 0, B - A);
add_psp_penalty(xi, xj, B + C - A - D);
}
void add_all_true_profit(const vector<int> &xs, COST P) {
assert(P >= 0);
int y = (int)list.size();
list.resize(y + 1);
OFFSET -= P;
add_edge(S, y, P);
for (auto xi : xs) {
assert(xi >= 0 && xi < N);
add_edge(y, xi, INF);
}
}
void add_all_false_profit(const vector<int> &xs, COST P) {
assert(P >= 0);
int y = (int)list.size();
list.resize(y + 1);
OFFSET -= P;
add_edge(y, T, P);
for (auto xi : xs) {
assert(xi >= 0 && xi < N);
add_edge(xi, y, INF);
}
}
COST solve() {
return dinic() + OFFSET;
}
vector<bool> reconstruct() {
vector<bool> res(N, false), seen(list.size(), false);
queue<int> que;
seen[S] = true;
que.push(S);
while (!que.empty()) {
int v = que.front();
que.pop();
for (const auto &e : list[v]) {
if (e.cap && !seen[e.to]) {
if (e.to < N) res[e.to] = true;
seen[e.to] = true;
que.push(e.to);
}
}
}
return res;
}
friend ostream& operator << (ostream& s, const TwoVariableSubmodularOpt &G) {
const auto &edges = G.get_edges();
for (const auto &e : edges) s << e << endl;
return s;
}
private:
struct Edge {
int rev, from, to;
COST cap, icap, flow;
Edge(int r, int f, int t, COST c)
: rev(r), from(f), to(t), cap(c), icap(c), flow(0) {}
void reset() { cap = icap, flow = 0; }
friend ostream& operator << (ostream& s, const Edge& E) {
return s << E.from << "->" << E.to << '(' << E.flow << '/' << E.icap << ')';
}
};
int N, S, T;
COST OFFSET, INF;
vector<vector<Edge>> list;
vector<pair<int,int>> pos;
Edge &get_rev_edge(const Edge &e) {
if (e.from != e.to) return list[e.to][e.rev];
else return list[e.to][e.rev + 1];
}
Edge &get_edge(int i) {
return list[pos[i].first][pos[i].second];
}
const Edge &get_edge(int i) const {
return list[pos[i].first][pos[i].second];
}
vector<Edge> get_edges() const {
vector<Edge> edges;
for (int i = 0; i < (int)pos.size(); ++i) {
edges.push_back(get_edge(i));
}
return edges;
}
void add_edge(int from, int to, COST cap) {
if (!cap) return;
pos.emplace_back(from, (int)list[from].size());
list[from].push_back(Edge((int)list[to].size(), from, to, cap));
list[to].push_back(Edge((int)list[from].size() - 1, to, from, 0));
}
COST dinic(COST limit_flow) {
COST current_flow = 0;
vector<int> level((int)list.size(), -1), iter((int)list.size(), 0);
auto bfs = [&]() -> void {
level.assign((int)list.size(), -1);
level[S] = 0;
queue<int> que;
que.push(S);
while (!que.empty()) {
int v = que.front();
que.pop();
for (const Edge &e : list[v]) {
if (level[e.to] < 0 && e.cap > 0) {
level[e.to] = level[v] + 1;
if (e.to == T) return;
que.push(e.to);
}
}
}
};
auto dfs = [&](auto self, int v, COST up_flow) {
if (v == T) return up_flow;
COST res_flow = 0;
for (int &i = iter[v]; i < (int)list[v].size(); ++i) {
Edge &e = list[v][i], &re = get_rev_edge(e);
if (level[v] >= level[e.to] || e.cap == 0) continue;
COST flow = self(self, e.to, min(up_flow - res_flow, e.cap));
if (flow <= 0) continue;
res_flow += flow;
e.cap -= flow, e.flow += flow;
re.cap += flow, re.flow -= flow;
if (res_flow == up_flow) break;
}
return res_flow;
};
while (current_flow < limit_flow) {
bfs();
if (level[T] < 0) break;
iter.assign((int)iter.size(), 0);
while (current_flow < limit_flow) {
COST flow = dfs(dfs, S, limit_flow - current_flow);
if (!flow) break;
current_flow += flow;
}
}
return current_flow;
};
COST dinic() {
return dinic(numeric_limits<COST>::max());
}
};
int main() {
const int VAL = 500;
int N, M;
cin >> N >> M;
vector<int> X(N), Y(N);
for (int i = 0; i < N; ++i) cin >> X[i];
for (int i = 0; i < N; ++i) cin >> Y[i];
vector A(M, vector(N, 0));
for (int i = 0; i < M; ++i) for (int j = 0; j < N; ++j) cin >> A[i][j];
const long long INF = 1LL<<50;
TwoVariableSubmodularOpt<long long> tvs(N * VAL * 2 + M, INF);
for (int i = 0; i < N; ++i) {
for (int j = 0; j < VAL; ++j) {
tvs.add_single_cost(i * VAL + j, 0, 1);
tvs.add_single_cost(N * VAL + i * VAL + j, 1, 0);
if (j + 1 < VAL) {
tvs.add_psp_constraint(i * VAL + j + 1, i * VAL + j);
tvs.add_psp_constraint(N * VAL + i * VAL + j, N * VAL + i * VAL + j + 1);
}
}
tvs.add_single_cost(i * VAL + X[i] - 1, INF, 0);
tvs.add_single_cost(N * VAL + i * VAL + Y[i] - 1, 0, INF);
}
for (int j = 0; j < M; ++j) {
for (int i = 0; i < N; ++i) {
tvs.add_psp_constraint(N * VAL * 2 + j, i * VAL + A[j][i] - 1);
tvs.add_psp_constraint(N * VAL + i * VAL + A[j][i] - 1, N * VAL * 2 + j);
}
}
cout << tvs.solve() << endl;
}