Weighted matroid intersection (using Dijkstra's algorithm)

combinatorial_opt/matroid_intersection.hpp

@@ -3,7 +3,7 @@
 #include <cassert>
 #include <vector>
 
-// Find augmenting path of matroid intersection.
+// Find minimum weight augmenting path of matroid intersection.
 // m1, m2: matroids
 // I: independent set (will be updated if augmenting path is found)
 //
@@ -46,7 +46,7 @@ bool matroid_intersection_augment(Matroid1 &m1, Matroid2 &m2, std::vector<bool>
     }
 }
 
-// (Min weight) matroid intersection solver
+// Minimum weight matroid intersection solver
 // Algorithm based on http://dopal.cs.uec.ac.jp/okamotoy/lect/2015/matroid/
 // Complexity: O(Cn^2 + n^3) (C : circuit query, non-weighted)
 template <class Matroid1, class Matroid2, class T = int>

combinatorial_opt/matroid_intersection_dijkstra.hpp

@@ -0,0 +1,115 @@
+#pragma once
+
+#include <cassert>
+#include <queue>
+#include <vector>
+
+// Find maximum weight size k + 1 intersection of two matroids using Dijkstra's algorithm
+// Return `true` iff larger intersection is found.
+// Complexity: O(Cn + nk log n) (C: circuit query)
+template <class Matroid1, class Matroid2, class T = int>
+bool augment_matroid_intersection_dijkstra(
+    Matroid1 &m1,                 // Matroid, size n, updated
+    Matroid2 &m2,                 // Matroid, size n, updated
+    std::vector<bool> &I,         // Size k maximum weight intersection, size n, updated
+    const std::vector<T> &weight, // Weights of elements, size n
+    std::vector<T> &potential     // Potential, size n, updated
+) {
+    const int n = I.size();
+
+    assert((int)m1.size() == n);
+    assert((int)m2.size() == n);
+    assert((int)weight.size() == n);
+    assert((int)potential.size() == n);
+
+    m1.set(I);
+    m2.set(I);
+
+    {
+        int arghi = -1;
+        for (int e = 0; e < n; ++e) {
+            if (I.at(e)) continue;
+            if (arghi < 0 or weight.at(arghi) < weight.at(e)) arghi = e;
+        }
+        if (arghi < 0) return false;
+        if (m1.circuit(arghi).empty() and m2.circuit(arghi).empty()) {
+            I.at(arghi) = true;
+            return true;
+        }
+    }
+
+    auto l = [&](int e) -> T { return e < n ? (I.at(e) ? weight.at(e) : -weight.at(e)) : T(); };
+    auto pot = [&](int e) -> T { return e < n ? potential.at(e) : T(); };
+    auto edge_len = [&](int s, int t) -> T { return l(t) - pot(t) + pot(s); };
+
+    const int gs = n, gt = n + 1;
+    std::vector<int> on_set;
+    for (int e = 0; e < n; ++e) {
+        if (I.at(e)) on_set.push_back(e);
+    }
+
+    // Find minimum weight (& minimum num. of vertices) gs-gt path
+    using Dist = std::pair<T, int>; // (sum of weights, num. of vertices)
+    std::vector<Dist> dp(gt + 1, {-1, -1});
+    std::vector<int> prv(gt + 1, -1); // prv[i] >= 0 => i is reachable (i != gs)
+
+    using Tpl = std::pair<Dist, int>;
+    std::priority_queue<Tpl, std::vector<Tpl>, std::greater<Tpl>> pq; // (dist, len, now)
+    std::vector<std::vector<int>> to(dp.size());
+
+    for (int e = 0; e < n; ++e) {
+        if (I.at(e)) continue;
+
+        const auto c1 = m1.circuit(e), c2 = m2.circuit(e);
+
+        if (c1.empty()) {
+            to.at(e).push_back(gt);
+            for (int f : on_set) to.at(e).push_back(f);
+        }
+        for (int f : c1) {
+            if (f != e) to.at(e).push_back(f);
+        }
+
+        if (c2.empty()) {
+            dp.at(e) = Dist{edge_len(gs, e), 1};
+            prv.at(e) = gs;
+            pq.emplace(dp.at(e), e);
+        }
+        for (int f : c2) {
+            if (f != e) to.at(f).push_back(e);
+        }
+    }
+
+    while (!pq.empty()) {
+        const auto [dnow, now] = pq.top();
+        pq.pop();
+        if (prv.at(now) >= 0 and dp.at(now) < dnow) continue;
+
+        for (int nxt : to.at(now)) {
+            const auto w = edge_len(now, nxt);
+            // if (w < T() and now < n and nxt < n) assert(false); // for debug
+
+            Dist dnxt(dnow.first + w, dnow.second + 1);
+
+            if (prv.at(nxt) < 0 or dnxt < dp.at(nxt)) {
+                dp.at(nxt) = dnxt;
+                prv.at(nxt) = now;
+                if (nxt != gt) pq.emplace(dnxt, nxt);
+            }
+        }
+    }
+
+    if (prv.at(gt) < 0) return false;
+
+    for (int e = 0; e < n; ++e) {
+        auto [dist, len] = dp.at(e);
+        if (len >= 0) potential.at(e) += dist;
+    }
+
+    for (int e = prv.at(gt); e != gs; e = prv.at(e)) {
+        potential.at(e) -= l(e);
+        I.at(e) = !I.at(e);
+    }
+
+    return true;
+}

combinatorial_opt/matroid_intersection_dijkstra.md

@@ -0,0 +1,35 @@
+---
+title: Weighted matroid intersection using Dijkstra's algorithm
+documentation_of: ./matroid_intersection_dijkstra.hpp
+---
+
+重み付きマトロイド交叉（交差）問題 (matroid intersection) のソルバ．二つのマトロイド $M\_{1} = (E, \mathcal{I}\_{1}), M_{2} = (E, \mathcal{I}\_{2})$ とそのサイズ $k$ 最大重み共通独立集合 $I\_k$ ，および $I\_k$ の導出に使われたポテンシャルを入力とし，サイズ $k + 1$ の最大重み共通独立集合を求めるか，これが存在しないことを示す．引数として与えたポテンシャルは適宜更新される．
+
+計算量は $n = \|E\|$，マトロイドクラスのサーキットクエリ一回あたりの計算量を $c$ として $O(nc + nk \log n)$ である．ポテンシャルの利用により内部で解く最短経路問題の辺重みを全て非負にでき，高速な Dijkstra 法を利用できる．ただしコンテストでの実用上は Bellman–Ford 法ベースの `matroid_intersection_augment()` の方が高速に動作するようである．
+
+## 使用方法
+
+```cpp
+UserDefinedMatroid m1, m2;
+vector<T> weight(n);
+vector<T> potential(n);
+vector<bool> I(n);
+
+assert(m1.size() == n);
+assert(m2.size() == n);
+
+while (augment_matroid_intersection_dijkstra(m1, m2, I, weight, potential)) {
+    continue;
+}
+```
+
+## 問題例
+
+- [AOJ GRL_2_B: 最小全域有向木](https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=GRL_2_B) 分割マトロイドとグラフマトロイドの重み付き交差でも解ける．
+- [AOJ 1605: Bridge Construction Planning （橋の建造計画）](https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1605&lang=ja) 分割マトロイドとグラフマトロイドの重み付き交差．
+- [2021 ICPC Asia Taiwan Online Programming Contest I. ICPC Kingdom - Codeforces](http://codeforces.com/gym/103373/problem/I) 横断マトロイドとグラフマトロイドの重み付き交差（横断マトロイドを考える代わりに辺を仮想的に増やし，分割マトロイドと解釈することも可能）．
+
+## 文献・リンク集
+
+- [1] C. Brezovec, G. Cornuéjols, and F. Glover, "Two algorithms for weighted matroid intersection,"
+  Mathematical Programming, 36(1), 39-53, 1986.

combinatorial_opt/test/matroid_intersection_dijkstra.aoj1605.test.cpp

@@ -0,0 +1,46 @@
+#define PROBLEM "https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1605"
+#include "../matroid_intersection_dijkstra.hpp"
+#include "../matroids/graphic_matroid.hpp"
+#include "../matroids/partition_matroid.hpp"
+#include <iostream>
+#include <numeric>
+#include <utility>
+#include <vector>
+using namespace std;
+
+int main() {
+    cin.tie(nullptr);
+    ios::sync_with_stdio(false);
+    while (true) {
+        int N, M, K;
+        cin >> N >> M >> K;
+        if (N == 0) break;
+        vector<vector<int>> partition(2);
+        vector<int> R{K, N - 1 - K};
+        vector<pair<int, int>> edges;
+        vector<int> weight;
+        for (int e = 0; e < M; ++e) {
+            int u, v, w;
+            char l;
+            cin >> u >> v >> w >> l;
+            --u, --v;
+            partition[l == 'B'].push_back(e);
+            edges.emplace_back(u, v);
+            weight.push_back(-w);
+        }
+        PartitionMatroid M1(edges.size(), partition, R);
+        GraphicMatroid M2(N, edges);
+
+        vector<int> potential(weight.size());
+        vector<bool> ret(weight.size());
+        while (augment_matroid_intersection_dijkstra(M1, M2, ret, weight, potential)) continue;
+        int ne = accumulate(ret.begin(), ret.end(), 0);
+        if (ne < N - 1) {
+            cout << "-1\n";
+        } else {
+            int sum = 0;
+            for (int e = 0; e < M; ++e) sum -= ret.at(e) * weight.at(e);
+            cout << sum << '\n';
+        }
+    }
+}

combinatorial_opt/test/matroid_intersection_dijkstra.aoj_grl_2_b.test.cpp

@@ -0,0 +1,35 @@
+#define PROBLEM "https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=GRL_2_B"
+#include "../matroid_intersection_dijkstra.hpp"
+#include "../matroids/graphic_matroid.hpp"
+#include "../matroids/partition_matroid.hpp"
+#include <iostream>
+#include <utility>
+#include <vector>
+using namespace std;
+
+int main() {
+    int V, E, r;
+    cin >> V >> E >> r;
+    vector<vector<int>> partition(V);
+    vector<int> R(V, 1);
+    R.at(r) = 0;
+    vector<pair<int, int>> edges;
+    vector<int> weight;
+    for (int e = 0; e < E; ++e) {
+        int s, t, w;
+        cin >> s >> t >> w;
+        partition.at(t).push_back(e);
+        edges.emplace_back(s, t);
+        weight.emplace_back(-w);
+    }
+    PartitionMatroid M1(E, partition, R);
+    GraphicMatroid M2(V, edges);
+
+    vector<int> potential(weight.size());
+    vector<bool> sol(weight.size());
+    while (augment_matroid_intersection_dijkstra(M1, M2, sol, weight, potential)) continue;
+
+    int ret = 0;
+    for (int e = 0; e < E; ++e) ret -= sol.at(e) * weight.at(e);
+    cout << ret << '\n';
+}