diff --git a/combinatorial_opt/matroid_intersection.hpp b/combinatorial_opt/matroid_intersection.hpp index 99ee9e67..5778c097 100644 --- a/combinatorial_opt/matroid_intersection.hpp +++ b/combinatorial_opt/matroid_intersection.hpp @@ -3,7 +3,7 @@ #include #include -// 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 } } -// (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 diff --git a/combinatorial_opt/matroid_intersection_dijkstra.hpp b/combinatorial_opt/matroid_intersection_dijkstra.hpp new file mode 100644 index 00000000..e4d76f57 --- /dev/null +++ b/combinatorial_opt/matroid_intersection_dijkstra.hpp @@ -0,0 +1,115 @@ +#pragma once + +#include +#include +#include + +// 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 +bool augment_matroid_intersection_dijkstra( + Matroid1 &m1, // Matroid, size n, updated + Matroid2 &m2, // Matroid, size n, updated + std::vector &I, // Size k maximum weight intersection, size n, updated + const std::vector &weight, // Weights of elements, size n + std::vector &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 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; // (sum of weights, num. of vertices) + std::vector dp(gt + 1, {-1, -1}); + std::vector prv(gt + 1, -1); // prv[i] >= 0 => i is reachable (i != gs) + + using Tpl = std::pair; + std::priority_queue, std::greater> pq; // (dist, len, now) + std::vector> 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; +} diff --git a/combinatorial_opt/matroid_intersection_dijkstra.md b/combinatorial_opt/matroid_intersection_dijkstra.md new file mode 100644 index 00000000..4d8a08a2 --- /dev/null +++ b/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 weight(n); +vector potential(n); +vector 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. diff --git a/combinatorial_opt/test/matroid_intersection_dijkstra.aoj1605.test.cpp b/combinatorial_opt/test/matroid_intersection_dijkstra.aoj1605.test.cpp new file mode 100644 index 00000000..2f18e812 --- /dev/null +++ b/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 +#include +#include +#include +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> partition(2); + vector R{K, N - 1 - K}; + vector> edges; + vector 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 potential(weight.size()); + vector 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'; + } + } +} diff --git a/combinatorial_opt/test/matroid_intersection_dijkstra.aoj_grl_2_b.test.cpp b/combinatorial_opt/test/matroid_intersection_dijkstra.aoj_grl_2_b.test.cpp new file mode 100644 index 00000000..03897da2 --- /dev/null +++ b/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 +#include +#include +using namespace std; + +int main() { + int V, E, r; + cin >> V >> E >> r; + vector> partition(V); + vector R(V, 1); + R.at(r) = 0; + vector> edges; + vector 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 potential(weight.size()); + vector 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'; +}