Weighted matroid intersection using Dijkstra's algorithm
(combinatorial_opt/matroid_intersection_dijkstra.hpp)

• View this file on GitHub
• Last update: 2023-07-23 18:47:47+09:00
• Include: #include "combinatorial_opt/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 法を利用できる．

使用方法

UserDefinedMatroid m1, m2;
vector<T> weight(n);
vector<T> potential(n + 2);  // 2 個の補助頂点のポテンシャルの情報も持つ
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: 最小全域有向木 分割マトロイドとグラフマトロイドの重み付き交差でも解ける．
• Deltix Round, Summer 2021 H. DIY Tree - Codeforces グラフマトロイドと分割マトロイドの最小重み共通独立集合問題に帰着される．
• AOJ 1605: Bridge Construction Planning （橋の建造計画） 分割マトロイドとグラフマトロイドの重み付き交差．
• 2021 ICPC Asia Taiwan Online Programming Contest I. ICPC Kingdom - Codeforces 横断マトロイドとグラフマトロイドの重み付き交差（横断マトロイドを考える代わりに辺を仮想的に増やし，分割マトロイドと解釈することも可能）．

文献・リンク集

• [1] C. Brezovec, G. Cornuéjols, and F. Glover, “Two algorithms for weighted matroid intersection,” Mathematical Programming, 36(1), 39-53, 1986.

Verified with

• combinatorial_opt/test/matroid_intersection_dijkstra.aoj1605.test.cpp
• combinatorial_opt/test/matroid_intersection_dijkstra.aoj_grl_2_b.test.cpp

Code

#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 + 2, updated
) {
    const int n = I.size();

    assert((int)m1.size() == n);
    assert((int)m2.size() == n);
    assert((int)weight.size() == n);
    assert(potential.empty() or ((int)potential.size() == n) or ((int)potential.size() == n + 2));

    m1.set(I);
    m2.set(I);

    potential.resize(n + 2);

    auto l = [&](int e) -> T { return e < n ? (I.at(e) ? weight.at(e) : -weight.at(e)) : T(); };
    auto edge_len = [&](int s, int t) -> T { return l(t) - potential.at(t) + potential.at(s); };

    if (true) { // 自明な追加が可能かチェック（省略してもアルゴリズムは正当）
        int max_elem = -1;
        for (int e = 0; e < n; ++e) {
            if (!I.at(e) and (max_elem < 0 or weight.at(max_elem) < weight.at(e))) max_elem = e;
        }
        if (max_elem < 0) return false;
        for (int e = 0; e < n; ++e) {
            if (!I.at(e) and weight.at(e) == weight.at(max_elem) and m1.circuit(e).empty() and
                m2.circuit(e).empty()) {
                potential.at(e) -= l(e);
                I.at(e) = true;
                return true;
            }
        }
    }

    // Find minimum length (& minimum num. of vertices) gs-gt path
    const int gs = n, gt = n + 1;
    std::vector<std::vector<int>> to(gt + 1);

    bool has_gs_edge = false, has_gt_edge = false;

    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);
            if (!has_gt_edge) {
                has_gt_edge = true;
                potential.at(gt) = potential.at(e);
            }
            if (T el = edge_len(e, gt); el < T()) potential.at(gt) += el;
        }
        for (int f : c1) {
            if (f != e) to.at(e).push_back(f);
        }

        if (c2.empty()) {
            to.at(gs).push_back(e);
            if (!has_gs_edge) {
                has_gs_edge = true;
                potential.at(gs) = potential.at(e) - l(e);
            }
            if (T el = edge_len(gs, e); el < T()) potential.at(gs) -= el;
        }
        for (int f : c2) {
            if (f != e) to.at(f).push_back(e);
        }
    }

    if (const T e0 = potential.at(gs); e0 != T()) {
        for (auto &p : potential) p -= e0;
    }

    if (!has_gs_edge or !has_gt_edge) return false;

    std::vector<bool> potential_fixed(gt + 1);

    T potential_add_unfixed_es = T();

    auto fix_potential = [&](int e) -> void {
        assert(!potential_fixed.at(e));
        potential_fixed.at(e) = true;
        potential.at(e) += potential_add_unfixed_es;
    };

    std::priority_queue<std::pair<T, int>, std::vector<std::pair<T, int>>, std::greater<>> pq;
    std::vector<T> dijkstra(gt + 1);
    std::vector<int> prv(gt + 1, -1);

    pq.emplace(T(), gs);

    while (!pq.empty()) {
        const int e = pq.top().second;
        pq.pop();
        if (potential_fixed.at(e)) continue;
        if (e != gs) potential_add_unfixed_es = edge_len(prv.at(e), e);

        std::vector<std::pair<int, int>> push_cands;

        auto rec = [&](auto &&self, int cur) -> bool {
            if (cur == gt) return true;
            fix_potential(cur);

            for (int nxt : to.at(cur)) {
                if (potential_fixed.at(nxt)) continue;

                const T len = edge_len(cur, nxt) - potential_add_unfixed_es;
                // if (len < T()) std::cerr << cur << ' ' << nxt << ' ' << len << std::endl;
                assert(len >= T());

                if (len == T()) {
                    prv.at(nxt) = cur;
                    if (self(self, nxt)) return true;
                } else {
                    if (prv.at(nxt) == -1 or potential_add_unfixed_es + len < dijkstra.at(nxt)) {
                        dijkstra.at(nxt) = potential_add_unfixed_es + len;
                        prv.at(nxt) = cur;
                        push_cands.emplace_back(nxt, cur);
                    }
                }
            }
            return false;
        };
        if (rec(rec, e)) break;

        for (auto [nxt, now] : push_cands) {
            if (prv.at(nxt) == now) pq.emplace(dijkstra.at(nxt), nxt);
        }
    }

    for (int e = 0; e < gt + 1; ++e) {
        if (!potential_fixed.at(e)) fix_potential(e);
    }

    if (prv.at(gt) < 0) return false;

    prv.assign(gt + 1, -1);
    std::queue<int> q;
    q.push(gs);

    for (int now = q.front(); now != gt; now = q.front()) {
        q.pop();
        for (int nxt : to.at(now)) {
            if (prv.at(nxt) == -1 and edge_len(now, nxt) == T()) {
                prv.at(nxt) = now;
                q.push(nxt);
            }
        }
    }

    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;
}

#line 2 "combinatorial_opt/matroid_intersection_dijkstra.hpp"

#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 + 2, updated
) {
    const int n = I.size();

    assert((int)m1.size() == n);
    assert((int)m2.size() == n);
    assert((int)weight.size() == n);
    assert(potential.empty() or ((int)potential.size() == n) or ((int)potential.size() == n + 2));

    m1.set(I);
    m2.set(I);

    potential.resize(n + 2);

    auto l = [&](int e) -> T { return e < n ? (I.at(e) ? weight.at(e) : -weight.at(e)) : T(); };
    auto edge_len = [&](int s, int t) -> T { return l(t) - potential.at(t) + potential.at(s); };

    if (true) { // 自明な追加が可能かチェック（省略してもアルゴリズムは正当）
        int max_elem = -1;
        for (int e = 0; e < n; ++e) {
            if (!I.at(e) and (max_elem < 0 or weight.at(max_elem) < weight.at(e))) max_elem = e;
        }
        if (max_elem < 0) return false;
        for (int e = 0; e < n; ++e) {
            if (!I.at(e) and weight.at(e) == weight.at(max_elem) and m1.circuit(e).empty() and
                m2.circuit(e).empty()) {
                potential.at(e) -= l(e);
                I.at(e) = true;
                return true;
            }
        }
    }

    // Find minimum length (& minimum num. of vertices) gs-gt path
    const int gs = n, gt = n + 1;
    std::vector<std::vector<int>> to(gt + 1);

    bool has_gs_edge = false, has_gt_edge = false;

    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);
            if (!has_gt_edge) {
                has_gt_edge = true;
                potential.at(gt) = potential.at(e);
            }
            if (T el = edge_len(e, gt); el < T()) potential.at(gt) += el;
        }
        for (int f : c1) {
            if (f != e) to.at(e).push_back(f);
        }

        if (c2.empty()) {
            to.at(gs).push_back(e);
            if (!has_gs_edge) {
                has_gs_edge = true;
                potential.at(gs) = potential.at(e) - l(e);
            }
            if (T el = edge_len(gs, e); el < T()) potential.at(gs) -= el;
        }
        for (int f : c2) {
            if (f != e) to.at(f).push_back(e);
        }
    }

    if (const T e0 = potential.at(gs); e0 != T()) {
        for (auto &p : potential) p -= e0;
    }

    if (!has_gs_edge or !has_gt_edge) return false;

    std::vector<bool> potential_fixed(gt + 1);

    T potential_add_unfixed_es = T();

    auto fix_potential = [&](int e) -> void {
        assert(!potential_fixed.at(e));
        potential_fixed.at(e) = true;
        potential.at(e) += potential_add_unfixed_es;
    };

    std::priority_queue<std::pair<T, int>, std::vector<std::pair<T, int>>, std::greater<>> pq;
    std::vector<T> dijkstra(gt + 1);
    std::vector<int> prv(gt + 1, -1);

    pq.emplace(T(), gs);

    while (!pq.empty()) {
        const int e = pq.top().second;
        pq.pop();
        if (potential_fixed.at(e)) continue;
        if (e != gs) potential_add_unfixed_es = edge_len(prv.at(e), e);

        std::vector<std::pair<int, int>> push_cands;

        auto rec = [&](auto &&self, int cur) -> bool {
            if (cur == gt) return true;
            fix_potential(cur);

            for (int nxt : to.at(cur)) {
                if (potential_fixed.at(nxt)) continue;

                const T len = edge_len(cur, nxt) - potential_add_unfixed_es;
                // if (len < T()) std::cerr << cur << ' ' << nxt << ' ' << len << std::endl;
                assert(len >= T());

                if (len == T()) {
                    prv.at(nxt) = cur;
                    if (self(self, nxt)) return true;
                } else {
                    if (prv.at(nxt) == -1 or potential_add_unfixed_es + len < dijkstra.at(nxt)) {
                        dijkstra.at(nxt) = potential_add_unfixed_es + len;
                        prv.at(nxt) = cur;
                        push_cands.emplace_back(nxt, cur);
                    }
                }
            }
            return false;
        };
        if (rec(rec, e)) break;

        for (auto [nxt, now] : push_cands) {
            if (prv.at(nxt) == now) pq.emplace(dijkstra.at(nxt), nxt);
        }
    }

    for (int e = 0; e < gt + 1; ++e) {
        if (!potential_fixed.at(e)) fix_potential(e);
    }

    if (prv.at(gt) < 0) return false;

    prv.assign(gt + 1, -1);
    std::queue<int> q;
    q.push(gs);

    for (int now = q.front(); now != gt; now = q.front()) {
        q.pop();
        for (int nxt : to.at(now)) {
            if (prv.at(nxt) == -1 and edge_len(now, nxt) == T()) {
                prv.at(nxt) = now;
                q.push(nxt);
            }
        }
    }

    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;
}

Back to top page