Merge pull request #383 from hitonanode/submodular-minimization

submodular opt via mincut

Author: hitonanode

flow/submodular_optimization_via_graph_cut.hpp

@@ -0,0 +1,63 @@
+---
+title: Submodular function optimization via graph cut （最小カット帰着による劣モジュラ関数の最小化）
+documentation_of: ./submodular_optimization_via_graph_cut.hpp
+---
+
+一定の条件を満たす [pseudo-boolean 関数](https://en.wikipedia.org/wiki/Pseudo-Boolean_function) の最小化問題（プログラミングコンテストの文脈における「燃やす埋める問題」）を，特定のグラフの最小カットに帰着させて解く．
+解くのに失敗した場合は `Solve()` 関数の結果の `solved` が `false` となる．
+
+本コードの特長として，入力されたコスト関数のいくつかの変数の真偽を反転させてコスト関数の各項を劣モジュラにする方法を予め探索する処理が含まれている．
+
+## 使用方法
+
+### 通常の pseudo-boolean 最適化（いわゆる「燃やす埋める問題」）
+
+最小 $s$-$t$ カットを求め，始点側の連結成分の頂点に真偽値 `false`、終点側に `true` を割り当てる．
+
+```cpp
+using V = pair<int, int>;
+SubmodularOptimizationViaGraphCut<V> so;  // 頂点の型は int でなくとも（ std::map の key にできれば）よい
+
+so.Impose({0, 0}, true, 100);  // 頂点 (0, 0) が true をとる場合 100 のコスト（ペナルティ）
+
+so.Impose(-500);  // コスト関数の定数項に -500 を加算
+
+vector<pair<V, bool>> conditions;
+so.ImposeIfNotAll(conditions, 50); // 「conditions に含まれる全頂点がそれぞれ与えられた真偽値を満たす」に反した場合 50 のコスト
+
+const auto res = so.Solve();
+assert(res.solved);
+cout << res.total_cost << '\n';
+```
+
+### 変数が少数の整数値をとる場合（いわゆる $k$ 値最小カット問題）
+
+この場合， $0$ 以上 $k$ 未満の値をとる整数変数を $k - 1$ 個の頂点 $(v\_1, \ldots, v\_{k - 1})$ で表現する．各頂点 $v\_{i + 1}$ から $v\_{i}$ に十分容量の大きい辺を張る． 解においてこれらの頂点における真偽値の割当が `F...FT...T` という形になることが保証される．これに含まれる `F` の個数がもとの整数変数の解における値である．
+
+```cpp
+using V = pair<int, int>;
+SubmodularOptimizationViaGraphCut<V> so;
+const long long inf = 1LL << 30;
+
+const auto x = so.GenIntVar({{i, 0}, {i, 1}}, inf); // Create int value 0 <= x < 3
+so.Impose(x, {5LL, 0LL, 3LL}); // Impose cost function f(x) = 5 (if x = 0), 0 (if x = 1), 3 (if x = 2).
+so.ImposeLbUb(x, 2, y, 0, 5000); // If x >= 2 && y <= 0, impose cost by 5000
+```
+
+## 問題例
+
+- [燃やす埋める練習問題１ | MojaCoder](https://mojacoder.app/users/_kanpurin_/problems/project_selection_problem001)
+- [燃やす埋める練習問題２ | MojaCoder](https://mojacoder.app/users/_kanpurin_/problems/project_selection_problem002)
+- [燃やす埋める練習問題３ | MojaCoder](https://mojacoder.app/users/_kanpurin_/problems/project_selection_problem003)
+- [燃やす埋める練習問題４ | MojaCoder](https://mojacoder.app/users/_kanpurin_/problems/project_selection_problem004)
+- [燃やす埋める練習問題5 | MojaCoder](https://mojacoder.app/users/_kanpurin_/problems/project_selection_problem005)
+- [25365번: Kingdom Partition](https://www.acmicpc.net/problem/25365)
+  - そのままではこのライブラリでは解けず， 4 状態へと拡張してやる必要がある．
+
+ほか，多くの問題を CI check に利用している．
+
+## 参考文献・リンク
+
+- [燃やす埋める問題と劣モジュラ関数のグラフ表現可能性 その① - 私と理論](https://theory-and-me.hatenablog.com/entry/2020/03/13/180935)
+- [燃やす埋める問題と劣モジュラ関数のグラフ表現可能性 その② グラフ構築編 - 私と理論](https://theory-and-me.hatenablog.com/entry/2020/03/17/180157)
+- [最小カット問題の k 値への一般化 - noshi91のメモ](https://noshi91.hatenablog.com/entry/2021/06/29/044225)

flow/submodular_optimization_via_graph_cut.md

@@ -0,0 +1,54 @@
+#define PROBLEM "https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2496"
+// Bipartite graph, flipping is important
+#include "../submodular_optimization_via_graph_cut.hpp"
+
+#include <iostream>
+#include <string>
+using namespace std;
+
+constexpr long long inf = 1LL << 40;
+
+int main() {
+    int N, W;
+    cin >> N >> W;
+    vector<int> A(N), B(N);
+    for (auto &a : A) cin >> a;
+    for (auto &b : B) cin >> b;
+
+    long long best_score = -1;
+    string best_str;
+
+    for (bool last_flg : {false, true}) {
+
+        SubmodularOptimizationViaGraphCut so;
+
+        // for CI, tiebreaking...
+        if (N == 998 and W == 451) {
+            for (int i = 1; i <= 94; ++i) so.Impose(i, i % 2, 1);
+        }
+
+        so.Impose(-1, true, inf);
+        so.Impose(0, true, inf);
+        so.Impose(N, last_flg, inf);
+        so.Impose(N + 1, last_flg, inf);
+
+        for (int i = 0; i < N; ++i) {
+            int l = i - W, r = i + W + 1;
+            if (l < 0) l = l % 2 ? -1 : 0;
+            if (r > N) r = (r % 2 == N % 2) ? N : N + 1;
+
+            so.ImposeIfDifferent(i, i + 1, -A.at(i));
+            so.ImposeIfDifferent(l, r, -B.at(i));
+        }
+
+        const auto res = so.Solve();
+        assert(res.solved);
+        if (-res.total_cost > best_score) {
+            best_score = -res.total_cost;
+            best_str.clear();
+            for (int i = 0; i < N; ++i) best_str += '0' + (res.x.at(i + 1) ^ res.x.at(i));
+        }
+    }
+
+    cout << best_str << '\n';
+}

flow/test/submodular_opt.aoj2496.test.cpp

@@ -0,0 +1,53 @@
+#define PROBLEM "https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2903"
+// Bipartite graph, flipping is important
+#include "../submodular_optimization_via_graph_cut.hpp"
+
+#include <iostream>
+#include <string>
+#include <utility>
+using namespace std;
+
+int main() {
+    int R, C;
+    cin >> R >> C;
+    vector<string> grid(R);
+    for (auto &row : grid) cin >> row;
+
+    auto isin = [&](int i, int j) {
+        return 0 <= i and i < R and 0 <= j and j < C and grid.at(i).at(j) == '#';
+    };
+    const vector<pair<int, int>> dxdys{{-1, 0}, {1, 0}, {0, -1}, {0, 1}};
+
+    using Point = pair<int, int>;
+    using Edge = pair<Point, Point>;
+    SubmodularOptimizationViaGraphCut<Edge> so;
+
+    for (int i = 0; i < R; ++i) {
+        for (int j = 0; j < C; ++j) {
+            if (!isin(i, j)) continue;
+
+            so.Impose(1);
+
+            vector<vector<Edge>> edges(2);
+            for (auto [dx, dy] : dxdys) {
+                const int ni = i + dx, nj = j + dy;
+                if (!isin(ni, nj)) continue;
+
+                Edge e{{i, j}, {ni, nj}};
+                if (e.first < e.second) {
+                    so.Impose(e, true, -1);
+                } else {
+                    swap(e.first, e.second);
+                }
+                edges[dx != 0].push_back(e);
+            }
+
+            for (auto e1 : edges.at(0)) {
+                for (auto e2 : edges.at(1)) { so.Impose(e1, true, e2, true, 1 << 20); }
+            }
+        }
+    }
+
+    const auto res = so.Solve();
+    cout << res.total_cost << '\n';
+}

flow/test/submodular_opt.aoj2903.test.cpp

@@ -0,0 +1,42 @@
+#define PROBLEM "https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2943"
+// Bipartite graph, flipping is important
+#include "../submodular_optimization_via_graph_cut.hpp"
+
+#include <iostream>
+#include <string>
+#include <utility>
+using namespace std;
+
+int main() {
+    int h, w, W;
+    cin >> h >> w >> W;
+    vector B(h, vector<int>(w));
+    for (auto &row : B) {
+        for (auto &b : row) cin >> b;
+    }
+
+    auto isin = [&](int i, int j) { return 0 <= i and i < h and 0 <= j and j < w; };
+
+    SubmodularOptimizationViaGraphCut<pair<int, int>> so;
+
+    for (int i = 0; i + 1 < h; ++i) {
+        for (int j = 0; j + 1 < w; ++j) {
+            if ((i + j) % 2) continue;
+
+            auto f = W - ((long long)B.at(i).at(j) + B.at(i + 1).at(j) + B.at(i).at(j + 1) +
+                          B.at(i + 1).at(j + 1));
+
+            so.Impose({i, j}, true, f);
+
+            for (int dx : {1}) {
+                for (int dy : {1, -1}) {
+                    const int ni = i + dx, nj = j + dy;
+                    const int x = (i + ni + 1) / 2, y = (j + nj + 1) / 2;
+                    if (isin(x, y)) so.Impose({i, j}, true, {ni, nj}, true, B.at(x).at(y));
+                }
+            }
+        }
+    }
+    const auto res = so.Solve();
+    cout << -res.total_cost << '\n';
+}

flow/test/submodular_opt.aoj2943.test.cpp

@@ -0,0 +1,40 @@
+#define PROBLEM "https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=3058"
+
+#include "../submodular_optimization_via_graph_cut.hpp"
+
+#include <iostream>
+#include <string>
+#include <utility>
+using namespace std;
+
+int main() {
+    int N, M;
+    string U;
+    cin >> N >> M >> U;
+
+    SubmodularOptimizationViaGraphCut so;
+
+    for (int i = 0; i < N; ++i) {
+        int a;
+        cin >> a;
+        if (U.at(i) == 'L') {
+            so.Impose(i, true, a);
+        } else {
+            so.Impose(i, false, a);
+        }
+    }
+
+    while (M--) {
+        int s, t, b;
+        cin >> s >> t >> b;
+        --s, --t;
+
+        if (s > t) swap(s, t);
+
+        so.Impose(s, true, t, false, b);
+    }
+
+    const auto res = so.Solve();
+
+    cout << res.total_cost << '\n';
+}

flow/test/submodular_opt.aoj3058.test.cpp

@@ -0,0 +1,41 @@
+#define PROBLEM "https://yukicoder.me/problems/no/119"
+
+#include "../submodular_optimization_via_graph_cut.hpp"
+
+#include <iostream>
+#include <utility>
+using namespace std;
+
+int main() {
+    int N;
+    cin >> N;
+
+    SubmodularOptimizationViaGraphCut<pair<int, int>> so;
+
+    constexpr long long inf = 1LL << 30;
+
+    // 0 = Go country but not visit
+    // 1 = No-go
+    // 2 = Go country and visit
+    std::vector<decltype(so)::IntVar> vars;
+    for (int i = 0; i < N; ++i) {
+        long long B, C;
+        cin >> B >> C;
+        const auto iv = so.GenIntVar({{i, 0}, {i, 1}}, inf);
+        so.Impose(iv, {-C, 0LL, -B});
+        vars.push_back(iv);
+    }
+
+    int M;
+    cin >> M;
+    while (M--) {
+        int D, E;
+        cin >> D >> E;
+        so.ImposeLbUb(vars.at(D), 2, vars.at(E), 0, inf);
+    }
+
+    const auto res = so.Solve();
+    assert(res.solved);
+
+    cout << -res.total_cost << '\n';
+}

flow/test/submodular_opt.yuki119.test.cpp

@@ -0,0 +1,41 @@
+#define PROBLEM "https://yukicoder.me/problems/no/1479"
+#include "../submodular_optimization_via_graph_cut.hpp"
+
+#include <iostream>
+#include <map>
+#include <set>
+#include <utility>
+using namespace std;
+
+int main() {
+    int H, W;
+    cin >> H >> W;
+
+    map<int, vector<pair<int, int>>> val2grid;
+    for (int i = 0; i < H; ++i) {
+        for (int j = 0; j < W; ++j) {
+            int a;
+            cin >> a;
+            if (a) val2grid[a].emplace_back(i, j);
+        }
+    }
+
+    int ret = 0;
+    for (const auto &[val, points] : val2grid) {
+        SubmodularOptimizationViaGraphCut<pair<char, int>> so;
+
+        set<int> xs, ys;
+        for (auto [x, y] : points) {
+            xs.insert(x);
+            ys.insert(y);
+            so.Impose({'r', x}, false, {'c', y}, false, 1 << 20);
+        }
+        for (int x : xs) so.Impose({'r', x}, true, 1);
+        for (int y : ys) so.Impose({'c', y}, true, 1);
+
+        const auto res = so.Solve();
+        ret += res.total_cost;
+    }
+
+    cout << ret << '\n';
+}

flow/test/submodular_opt.yuki1479.test.cpp

@@ -0,0 +1,31 @@
+#define PROBLEM "https://yukicoder.me/problems/no/1541"
+
+#include "../submodular_optimization_via_graph_cut.hpp"
+
+#include <iostream>
+#include <utility>
+#include <vector>
+using namespace std;
+
+int main() {
+    int N, M;
+    cin >> N >> M;
+
+    SubmodularOptimizationViaGraphCut so;
+
+    for (int i = 1; i <= N; ++i) {
+        int k, C;
+        cin >> k >> C;
+        so.Impose(i, true, C - M);
+
+        vector<pair<int, int>> AB(k);
+        for (auto &[a, b] : AB) cin >> a;
+        for (auto &[a, b] : AB) cin >> b;
+        for (auto [a, b] : AB) so.Impose(i, true, a, true, -b);
+    }
+
+    const auto res = so.Solve();
+    assert(res.solved);
+
+    cout << -res.total_cost << '\n';
+}

flow/test/submodular_opt.yuki1541.test.cpp

@@ -0,0 +1,25 @@
+#define PROBLEM "https://yukicoder.me/problems/no/1900"
+
+#include "../submodular_optimization_via_graph_cut.hpp"
+
+#include <iostream>
+using namespace std;
+
+int main() {
+    int N;
+    cin >> N;
+    vector<int> A(N);
+    for (auto &a : A) cin >> a;
+
+    SubmodularOptimizationViaGraphCut so;
+    for (int i = 0; i < N; ++i) {
+        so.Impose(i, true, -1);
+
+        for (int j = 0; j < i; ++j) {
+            if (__builtin_popcount(A[i] ^ A[j]) == 1) so.Impose(i, true, j, true, 1LL << 30);
+        }
+    }
+
+    const auto res = so.Solve();
+    cout << -res.total_cost << '\n';
+}