Edge disjoint minimum spanning forests

combinatorial_opt/edge_disjoint_min_spanning_forests.hpp

@@ -0,0 +1,121 @@
+#pragma once
+#include <algorithm>
+#include <cassert>
+#include <numeric>
+#include <utility>
+#include <vector>
+
+// Max size min weight two spanning forests
+// Complexity: O(NM + M log M)
+// Reference: https://www.slideshare.net/tmaehara/ss-17402143
+// Verified:
+// - https://www.codechef.com/submit/HAMEL
+// - https://community.topcoder.com/stat?c=problem_statement&pm=14909&rd=17198
+template <class Label, class W = int>
+std::pair<std::vector<bool>, std::vector<bool>>
+edge_disjoint_min_spanning_forests(const std::vector<std::pair<Label, Label>> &edges,
+                                   const std::vector<W> &ws = {}) {
+    assert(ws.empty() or ws.size() == edges.size());
+    const int M = edges.size();
+    std::vector<Label> lbl(M * 2);
+    for (int e = 0; e < M; e++) lbl[e * 2] = edges[e].first, lbl[e * 2 + 1] = edges[e].second;
+    std::sort(lbl.begin(), lbl.end());
+    lbl.erase(std::unique(lbl.begin(), lbl.end()), lbl.end());
+    const int N = lbl.size();
+
+    std::vector<std::pair<int, int>> uvs(M);
+    for (int e = 0; e < M; e++) {
+        int u = std::lower_bound(lbl.begin(), lbl.end(), edges[e].first) - lbl.begin();
+        int v = std::lower_bound(lbl.begin(), lbl.end(), edges[e].second) - lbl.begin();
+        uvs[e] = {u, v};
+    }
+
+    std::vector<int> es(M);
+    std::iota(es.begin(), es.end(), 0);
+    if (!ws.empty()) std::sort(es.begin(), es.end(), [&](int i, int j) { return ws[i] < ws[j]; });
+
+    std::vector<std::vector<bool>> I(2, std::vector<bool>(M));
+    std::vector<std::vector<std::pair<int, int>>> to(N);
+
+    int nb_accepted_edges = 0;
+    auto accept_edge = [&](int e) {
+        nb_accepted_edges++;
+        int u = uvs[e].first, v = uvs[e].second;
+        to[u].emplace_back(v, e);
+        to[v].emplace_back(u, e);
+    };
+
+    auto dfs = [&](int head, const std::vector<bool> &I) -> std::vector<std::pair<int, int>> {
+        std::vector<int> st{head};
+        std::vector<std::pair<int, int>> prv(N, {-1, -1});
+        prv[head] = {head, -1};
+        while (!st.empty()) {
+            int now = st.back();
+            st.pop_back();
+            for (auto p : to[now]) {
+                int nxt = p.first, e = p.second;
+                if (!I[e] or prv[nxt].first >= 0) continue;
+                prv[nxt] = {now, e}, st.push_back(nxt);
+            }
+        }
+        return prv;
+    };
+
+    std::vector<int> prveid(M, -1), visited(N);
+    std::vector<std::vector<int>> L(2);
+    std::vector<std::vector<std::pair<int, int>>> prv(2);
+    for (const int e : es) {
+        if (nb_accepted_edges > 2 * (N - 1)) break;
+        const int u = uvs[e].first, v = uvs[e].second;
+
+        bool found = false;
+
+        for (int d = 0; d < 2; d++) {
+            prv[d] = dfs(u, I[d]);
+            if (prv[d][v].first < 0) {
+                accept_edge(e);
+                I[d][e] = 1;
+                found = true;
+                break;
+            }
+        }
+        if (found) continue;
+
+        visited.assign(N, 0);
+        visited[u] = 1;
+        L[0] = {e}, L[1] = {};
+
+        int ehead = -1;
+        prveid[e] = -1;
+        for (int i = 0;; i ^= 1) {
+            if (L[i].empty()) break;
+            L[i ^ 1].clear();
+            while (!L[i].empty()) {
+                const int exy = L[i].back();
+                L[i].pop_back();
+                int x = uvs[exy].first, y = uvs[exy].second;
+                if (prv[i][x].first < 0 or prv[i][y].first < 0) {
+                    ehead = exy;
+                    break;
+                }
+                if (!visited[x]) std::swap(x, y);
+                while (!visited[y]) {
+                    int nxty = prv[i][y].first, nxte = prv[i][y].second;
+                    L[i ^ 1].push_back(nxte);
+                    visited[y] = 1;
+                    y = nxty;
+                    prveid[nxte] = exy;
+                }
+            }
+            if (ehead >= 0) {
+                accept_edge(e);
+                int c = I[0][ehead];
+                for (; ehead >= 0; ehead = prveid[ehead], c ^= 1) {
+                    I[c][ehead] = 1, I[c ^ 1][ehead] = 0;
+                }
+                break;
+            }
+        }
+    }
+    return {I[0], I[1]};
+}

combinatorial_opt/edge_disjoint_min_spanning_forests.md

@@ -0,0 +1,42 @@
+---
+title: Edge disjoint minimum spanning forests （無向グラフにおける最大辺数最小重みの辺素な二つの全域森）
+documentation_of: ./edge_disjoint_min_spanning_forests.hpp
+---
+
+$n$ 頂点 $m$ 辺の重み付き無向グラフを入力として，辺を共有しない二つの全域森であって辺数の和が最大となるもののうち重み最小のものを一つ構築する．
+本実装の計算量は重み付きのとき $O(nm + m \log m)$，重みなしのとき $O(nm)$．
+
+## アルゴリズムの概要
+
+アルゴリズムの大枠は [Matroid union （マトロイドの合併） | cplib-cpp](https://hitonanode.github.io/cplib-cpp/combinatorial_opt/matroid_union.hpp) と同様．
+
+## 使用方法
+
+```cpp
+vector<pint> edges;
+vector<int> weights;
+while (M--) {
+    int u, v, w;
+    cin >> u >> v >> w;
+    edges.emplace_back(u, v);
+    weights.push_back(w);
+}
+vector<bool> I1, I2;
+tie(I1, I2) = edge_disjoint_min_spanning_forests(edges); // 重みなし
+tie(I1, I2) = edge_disjoint_min_spanning_forests(edges, w); // 重み付き
+```
+
+返り値の `I1`, `I2` はともに長さ $m$ のベクトル．`I1[e]` の値が `true` のとき，$e$ 本目の辺が一つ目の全域森に含まれていることを示す（`I2` の解釈も同様）．
+
+## 問題例
+
+- [Hamel Paths \| CodeChef](https://www.codechef.com/problems/HAMEL) $n \le 64, m \le n(n - 1) / 2$, 最大 64 ケース．
+- [2018 TCO Round 3A 1000 ColoringEdgesDiv1](https://community.topcoder.com/stat?c=problem_statement&pm=14909&rd=17198) $n \le 1000, m = 3n/2$
+
+## 参考文献・リンク
+
+- [様々な全域木問題](https://www.slideshare.net/tmaehara/ss-17402143)
+- [1] J. Roskind & R. E. Tarjan, 
+  "A note on finding minimum-cost edge-disjoint spanning trees,"
+  Mathematics of Operations Research, 10(4), 701-708, 1985.
+  - 一般に $k$ 個の全域森にグラフを分割する $O(m \log m + k^2 n^2)$ のアルゴリズムが存在するらしい．

combinatorial_opt/matroid_union.md

@@ -28,7 +28,7 @@ $\|E\| = n$ として，$e = 1, \dots, n$ に $s$, $t$ を加えた $n + 2$ 頂
 
 - [Hamel Paths \| CodeChef](https://www.codechef.com/problems/HAMEL) 無向グラフから2つの全域木を作る．
 - [SRM 685 Div.1 450 FoxAirline2](https://community.topcoder.com/stat?c=problem_statement&pm=14194&rd=16689) 無向グラフから2つの全域木を作れるか判定する．
-- [2018 TCO Round 3A 1000 ColoringEdgesDiv1](https://community.topcoder.com/stat?c=problem_statement&pm=14909&rd=17198) 各頂点の次数が3の単純無向グラフ（$n \le 1000, m = 3n/2$）の辺を2つの全域森に分割する方法を構成する．最初に乱択で辺を追加して全域森2つの初期解を構成し，その後全ての辺を割り当て終えるまで `augment_union_matroid()` を真面目に使用するというヒューリスティックによって最大 200 ms 程度の実行時間で用意された全ケースに通る（ただし，後述するように無向グラフの全域森への分割にはより高速なアルゴリズムが存在する）．
+- [2018 TCO Round 3A 1000 ColoringEdgesDiv1](https://community.topcoder.com/stat?c=problem_statement&pm=14909&rd=17198) 各頂点の次数が3の単純無向グラフ（$n \le 1000, m = 3n/2$）の辺を2つの全域森に分割する方法を構成する．最初に乱択で辺を追加して全域森2つの初期解を構成し，その後全ての辺を割り当て終えるまで `augment_union_matroid()` を真面目に使用するというヒューリスティックによって [最大 200 ms 程度の実行時間で用意された全ケースに通る](https://vjudge.net/solution/32788901)（ただし，後述するように無向グラフの全域森への分割にはより高速なアルゴリズムが存在する：`edge_disjoint_min_spanning_forests.hpp` 参照）．
 
 ## 参考文献・リンク