SuccessiveShortestPath (最小費用流、最短路反復法)
(Src/Graph/Flow/SuccessiveShortestPath.hpp)

View this file on GitHub
Last update: 2024-01-04 02:51:07+09:00
Include: #include "Src/Graph/Flow/SuccessiveShortestPath.hpp"

概要

ライブラリの使い方

最小費用流は色々ゴニョゴニョしたいことが多いので、全メンバがpublicに置かれている。思う存分ゴニョゴニョすることができる。

テンプレート引数

template <class Cap, class Cost>

Cap: 容量の型
Cost: コストの型 (逆辺を張る都合上、符号付きの算術型しか許容しない)

コンストラクタ

(1) SuccessiveShortestPath() = default (1)
(2) SuccessiveShortestPath(usize n, usize m = usize{}) (2)

(2) $n$ は頂点数 (正整数)、 $m$ は辺数 (指定すると内部のコンテナをreserveする。普段使いでは無視して良いパラメータ)

計算量: $O(n + m)$

addEdge

u32 addEdge(u32 from, u32 to, const Cap& cap, const Cost& cost)

ネットワーク上に頂点fromから頂点toへの、容量capコストcostの辺を張る。

逆辺も張られる。

この辺に振られた一意なidが返る。逆辺のidは $\text{id} + 1$ である(使う上で覚える必要は無い)

計算量: $O(1)$

reverseId

static constexpr u32 reverseId(u32 i) noexcept

$i$ 番のidが振られた辺の逆辺のidを返す。

from

inline u32 from(u32 i) const noexcept

$i$ 番のidが振られた辺の湧き出し元の頂点番号を返す。

to

inline u32 to(u32 i) const noexcept

$i$ 番のidが振られた辺の湧き出し先の頂点番号を返す。

residual

inline const Cap& residual(u32 i) const noexcept

$i$ 番のidが振られた辺の残余容量を返す。

cost

inline const Cost& cost(u32 i) const noexcept

$i$ 番のidが振られた辺のコストを返す。

flowed

inline const Cap& flowed(u32 i) const noexcept

$i$ 番のidが振られた辺に流れた水の量を返す。

capacity

inline const Cap& capcacity(u32 i) const noexcept

$i$ 番のidが振られた辺の元々の容量を返す。

flow

bool flow(u32 s, u32 t, Cap flow)

$s$ から $t$ へ水を $\text{flow}$ だけ流す。 $\text{flow}$ 流せたら true を返す。

計算量: $O(\text{flow} \times \mid E\mid\log (\mid V\mid))$

負辺をaddEdgeしていると遅い。こういう時はbellmanford、dagdp -> updatePotential でポテンシャルを予め計算しておくと正常に動作する。

bellmanfordは $O(\midV\mid\midE\mid)$ で動作することに注意

maxflow

Cap maxflow(u32 s, u32 t)

$s$ から $t$ へ流せるだけ水を流す (最小費用最大流問題)

計算量: 最大流を $F$ として $O(F \times \mid E\mid\log (\mid V\mid))$

minCost

Cost minCost() const noexcept

現在のフローのコストを返す。

slope

std::vector<Cap> slope(u32 s, u32 t)

$s$ から $t$ への最大流を $F$ とする。 $s$ から $t$ へ水を $0$ 流した時、 $1$ 流した時、 …. 、 $F$ 流した時のフローの最小コストをstd::vector<Cost>で返す。

計算量: 最大流を $F$ として $O(F \times \mid E\mid\log (\mid V\mid))$

slopeACL

struct Line {
    Cap dn, up;
    Cost cost;
};
std::vector<Line> slopeACL(u32 s, u32 t)

atcoder::mcf_graphのslopeと似たやつを返す。最大流がでかい場合はslopeよりこっちを使った方が良さそう。

基本的にはここまでのメンバを使うことで事足りる。アドホックに色々したい、どうしても内部のメンバを呼び出したい時は下記のメンバも使えるかもしれない

edgeCost

inline Cost edgeCost(u32 i) const noexcept

$i$ 番のidが振られた辺の現在のポテンシャルを加味した辺のコストを返す。

relax

bool relax(u32 i)

$i$ 番目の辺でネットワークを緩和する。緩和されて最短経路が更新されたらtrueを返す。

現在の辺容量が $0$ ならば必ず緩和は失敗する

dijkstra

bool dijkstra(u32 s, u32 t)

始点 $s$ 、終点 $t$ でネットワークにダイクストラ法を適用する。容量が正の辺のみを考慮する。(終点 $t$ に到達したらループから抜けるみたいな枝刈りは入っていない)

始点から終点へ到達可能ならtrueを返す。

dagdp

bool dagdp(u32 s, u32 t)

始点 $s$ 、終点 $t$ でネットワークにdpで最短経路を更新する。容量が正の辺のみを考慮する。(終点 $t$ に到達したらループから抜けるみたいな枝刈りは入っていない)

始点から終点へ到達可能ならtrueを返す。

当然だが、ネットワークがDAGで無い場合の挙動は知らん。

全ての逆辺の辺容量が $0$ の時以外呼び出すべきで無いと思うが..

bellmanford

bool bellmanford(u32 s, u32 t)

始点 $s$ 、終点 $t$ でネットワークにベルマンフォード法を適用する。容量が正の辺のみを考慮する。

始点から終点へ到達可能ならtrueを返す。

updatePotential

void updatePotential()

ポテンシャルを更新する。

flush

Cap flush(u32 s, u32 t, Cap up = std::numeric_limits<Cap>::max())

始点 $s$ から終点 $t$ のパスで現在発見している最短経路へ、最大upの水を流す。

流せた水の量を返す。

dijkstra(s, t)、dagdp(s, t)、bellmanford(s, t)などを呼び最短経路を求めておくこと。

emplace

void emplace(u32 from, u32 to, const Cap& cap, const Cost& cost)

コンテナに頂点fromから頂点toへの、辺容量cap、辺コストcostの辺を直接挿入する。

大体の場合ではaddEdgeメンバで事足りる。こいつを直接呼び出す必要がある状況をあまり想定できていない

Depends on

Verified with

Code

#pragma once

#include "../../Template/TypeAlias.hpp"
#include "../../Utility/U32Pair.hpp"

#include <cassert>
#include <limits>
#include <queue>
#include <type_traits>
#include <utility>
#include <vector>

namespace zawa {

template <class Cap, class Cost>
class SuccessiveShortestPath {
public:
    static_assert(std::is_signed_v<Cost>, U"Cost must be signed");
    static constexpr Cost invalid{(std::numeric_limits<Cost>::max() >> 1) - 1};
    static constexpr u32 reverseId(u32 i) noexcept {
        return i ^ 1;
    }

    struct Edge {
        u32 from{}, to{};
        Cap residual{};
        Cost cost{};
        Edge() = default;
        Edge(u32 from, u32 to, const Cap& cap, const Cost& cost)
            : from{from}, to{to}, residual{cap}, cost{cost} {}
    };

    usize n_{}, m_{};
    std::vector<Edge> edges_;
    std::vector<std::vector<u32>> g_;
    std::vector<Cost> dist_, potential_;
    std::vector<U32Pair> prev_;
    Cost mcf_{};

    constexpr usize size() const noexcept {
        return n_;
    }
    constexpr usize edgeSize() const noexcept {
        return m_;
    }
    
    SuccessiveShortestPath() = default;
    SuccessiveShortestPath(usize n, usize m = usize{}) 
        : n_{n}, m_{}, edges_{}, g_(n), dist_(n), potential_(n), prev_(n), mcf_{} {
        g_.shrink_to_fit();
        dist_.shrink_to_fit();
        potential_.shrink_to_fit();
        prev_.shrink_to_fit();
        edges_.reserve(2 * m);
    }

    void emplace(u32 from, u32 to, const Cap& cap, const Cost& cost) {
        g_[from].emplace_back(m_);
        edges_.emplace_back(from, to, cap, cost);
        m_++;
    }

    u32 addEdge(u32 from, u32 to, const Cap& cap, const Cost& cost) {
        assert(from < size());
        assert(to < size());
        u32 res{static_cast<u32>(m_)};
        emplace(from, to, cap, cost);
        emplace(to, from, Cap{}, -cost);
        return res;
    }

    inline u32 from(u32 i) const noexcept {
        return edges_[i].from;
    }
    inline u32 to(u32 i) const noexcept {
        return edges_[i].to;
    }
    inline const Cap& residual(u32 i) const noexcept {
        return edges_[i].residual;
    }
    inline const Cost& cost(u32 i) const noexcept {
        return edges_[i].cost;
    }
    inline const Cap& flowed(u32 i) const noexcept {
        assert(i < edgeSize());
        return residual(i ^ 1);
    }
    inline const Cap& capcacity(u32 i) const noexcept {
        assert(i < edgeSize());
        return residual(i) + flowed(i);
    }

    inline Cost edgeCost(u32 i) const noexcept {
        return potential_[from(i)] + cost(i) - potential_[to(i)];
    }
    bool relax(u32 i) {
        if (residual(i) > 0 and dist_[to(i)] > dist_[from(i)] + edgeCost(i)) {
            dist_[to(i)] = dist_[from(i)] + edgeCost(i);
            prev_[to(i)] = U32Pair{from(i), i};
            return true;
        }
        return false;
    }

    bool dijkstra(u32 s, u32 t) {
        assert(s < size());
        assert(t < size());
        std::fill(dist_.begin(), dist_.end(), invalid);
        dist_[s] = Cost{};
        using qt = std::pair<Cost, u32>;
        std::priority_queue<qt, std::vector<qt>, std::greater<qt>> que;
        que.emplace(dist_[s], s);
        while (que.size()) {
            auto [d, v]{que.top()};
            que.pop();
            if (dist_[v] < d) continue;
            for (u32 i : g_[v]) {
                if (relax(i)) {
                    que.emplace(dist_[to(i)], to(i));
                }
            }
        }
        return dist_[t] < invalid;
    }

    bool dagdp(u32 s, u32 t) {
        assert(s < size());
        assert(t < size());
        std::fill(dist_.begin(), dist_.end(), invalid);
        dist_[s] = Cost{};
        std::vector<u32> ord;
        ord.reserve(size());
        std::vector<bool> vis(size());
        vis[s] = true;
        ord.push_back(s);
        for (u32 i{} ; i < ord.size() ; i++) {
            u32 v{ord[i]};
            for (auto e : g_[v]) {
                if (!vis[to(e)]) {
                    ord.push_back(to(e));
                    vis[to(e)] = true;
                }
                relax(e);
            }
        }
        return dist_[t] < invalid;
    }

    bool bellmanford(u32 s, u32 t) {
        assert(s < size());
        assert(t < size());
        std::fill(dist_.begin(), dist_.end(), invalid);
        dist_[s] = Cost{};
        for (u32 i{} ; i + 1 < size() ; i++) {
            for (u32 j{} ; j < edgeSize() ; j++) {
                relax(j);
            }
        }
        return dist_[t] < invalid;
    }

    void updatePotential() {
        for (u32 v{} ; v < size() ; v++) {
            potential_[v] = potential_[v] + (dist_[v] == invalid ? Cost{} : dist_[v]);
        }
    }

    Cap flush(u32 s, u32 t, Cap up = std::numeric_limits<Cap>::max()) {
        for (u32 v{t} ; v != s ; v = prev_[v].first()) {
            up = std::min(up, residual(prev_[v].second()));
        }
        for (u32 v{t} ; v != s ; v = prev_[v].first()) {
            edges_[prev_[v].second()].residual -= up;
            edges_[prev_[v].second() ^ 1].residual += up;
        }
        return up;
    }

    bool flow(u32 s, u32 t, Cap flow) {
        assert(s < size());
        assert(t < size());
        while (flow > 0 and dijkstra(s, t)) {
            updatePotential();
            Cap water{flush(s, t, flow)};
            mcf_ += potential_[t] * water;
            flow -= water;
        }
        return flow == 0;
    }

    Cap maxflow(u32 s, u32 t) {
        assert(s < size());
        assert(t < size());
        Cap flow{};
        while (dijkstra(s, t)) {
            updatePotential();
            Cap water{flush(s, t)};
            mcf_ += potential_[t] * water;
            flow += water;
        }
        return flow;
    }

    std::vector<Cost> slope(u32 s, u32 t) {
        assert(s < size());
        assert(t < size());
        Cap flow{};
        std::vector<Cost> res;
        while (dijkstra(s, t)) {
            updatePotential();
            Cap water{flush(s, t)};
            for (u32 i{} ; i < water ; i++) {
                res.emplace_back(mcf_);
                mcf_ += potential_[t];
                flow++;
            }
        }
        res.emplace_back(mcf_);
        return res;
    }

    struct Line {
        Cap dn{}, up{};
        Cost slope{};
        Line() = default;
        Line(Cap dn, Cap up, Cost slope) : dn{dn}, up{up}, slope{slope} {}
    };
    std::vector<Line> slopeACL(u32 s, u32 t) {
        assert(s < size());
        assert(t < size()); 
        Cap flow{};
        std::vector<Line> res;
        while (dijkstra(s, t)) {
            updatePotential();
            Cap water{flush(s, t)};
            mcf_ += potential_[t] * water;
            res.emplace_back(flow, flow + water, potential_[t]);
            flow += water;
        }
        return res;
    }

    Cost minCost() const noexcept {
        return mcf_;
    }
};

} // namespace zawa

#line 2 "Src/Graph/Flow/SuccessiveShortestPath.hpp"

#line 2 "Src/Template/TypeAlias.hpp"

#include <cstdint>
#include <cstddef>

namespace zawa {

using i16 = std::int16_t;
using i32 = std::int32_t;
using i64 = std::int64_t;
using i128 = __int128_t;

using u8 = std::uint8_t;
using u16 = std::uint16_t;
using u32 = std::uint32_t;
using u64 = std::uint64_t;

using usize = std::size_t;

} // namespace zawa
#line 2 "Src/Utility/U32Pair.hpp"

#line 4 "Src/Utility/U32Pair.hpp"

#include <functional>
#include <iostream>

namespace zawa {

class U32Pair {
private:
    static constexpr u32 SHIFT{32};
    static constexpr u32 MASK{static_cast<u32>((1LL << SHIFT) - 1)};
    u64 value_{};
public:
    constexpr U32Pair() {}
    constexpr U32Pair(u32 first, u32 second) {
        value_ = (static_cast<u64>(first) << SHIFT) | second;
    }
    constexpr u32 first() const noexcept {
        return static_cast<u32>(value_ >> SHIFT);
    }
    constexpr u32 second() const noexcept {
        return static_cast<u32>(value_ & MASK);
    }
    constexpr u64 combined() const noexcept {
        return value_;
    }
    constexpr U32Pair& operator=(const U32Pair& rhs) {
        value_ = rhs.value_;
        return *this;
    }
    friend constexpr bool operator==(const U32Pair& lhs, const U32Pair& rhs) {
        return lhs.value_ == rhs.value_;
    }
    friend constexpr bool operator!=(const U32Pair& lhs, const U32Pair& rhs) {
        return lhs.value_ != rhs.value_;
    }
    friend constexpr bool operator<(const U32Pair& lhs, const U32Pair& rhs) {
        return lhs.value_ < rhs.value_;
    }
    friend constexpr bool operator<=(const U32Pair& lhs, const U32Pair& rhs) {
        return lhs.value_ <= rhs.value_;
    }
    friend constexpr bool operator>(const U32Pair& lhs, const U32Pair& rhs) {
        return lhs.value_ > rhs.value_;
    }
    friend constexpr bool operator>=(const U32Pair& lhs, const U32Pair& rhs) {
        return lhs.value_ >= rhs.value_;
    }
    friend std::ostream& operator<<(std::ostream& os, const U32Pair& pair) {
        os << '(' << pair.first() << ',' << pair.second() << ')';
        return os;
    }
};

struct U32PairHash {
    usize operator()(const U32Pair& pair) const noexcept {
        return std::hash<u64>{}(pair.combined());
    }
};

} // namespace zawa
#line 5 "Src/Graph/Flow/SuccessiveShortestPath.hpp"

#include <cassert>
#include <limits>
#include <queue>
#include <type_traits>
#include <utility>
#include <vector>

namespace zawa {

template <class Cap, class Cost>
class SuccessiveShortestPath {
public:
    static_assert(std::is_signed_v<Cost>, U"Cost must be signed");
    static constexpr Cost invalid{(std::numeric_limits<Cost>::max() >> 1) - 1};
    static constexpr u32 reverseId(u32 i) noexcept {
        return i ^ 1;
    }

    struct Edge {
        u32 from{}, to{};
        Cap residual{};
        Cost cost{};
        Edge() = default;
        Edge(u32 from, u32 to, const Cap& cap, const Cost& cost)
            : from{from}, to{to}, residual{cap}, cost{cost} {}
    };

    usize n_{}, m_{};
    std::vector<Edge> edges_;
    std::vector<std::vector<u32>> g_;
    std::vector<Cost> dist_, potential_;
    std::vector<U32Pair> prev_;
    Cost mcf_{};

    constexpr usize size() const noexcept {
        return n_;
    }
    constexpr usize edgeSize() const noexcept {
        return m_;
    }
    
    SuccessiveShortestPath() = default;
    SuccessiveShortestPath(usize n, usize m = usize{}) 
        : n_{n}, m_{}, edges_{}, g_(n), dist_(n), potential_(n), prev_(n), mcf_{} {
        g_.shrink_to_fit();
        dist_.shrink_to_fit();
        potential_.shrink_to_fit();
        prev_.shrink_to_fit();
        edges_.reserve(2 * m);
    }

    void emplace(u32 from, u32 to, const Cap& cap, const Cost& cost) {
        g_[from].emplace_back(m_);
        edges_.emplace_back(from, to, cap, cost);
        m_++;
    }

    u32 addEdge(u32 from, u32 to, const Cap& cap, const Cost& cost) {
        assert(from < size());
        assert(to < size());
        u32 res{static_cast<u32>(m_)};
        emplace(from, to, cap, cost);
        emplace(to, from, Cap{}, -cost);
        return res;
    }

    inline u32 from(u32 i) const noexcept {
        return edges_[i].from;
    }
    inline u32 to(u32 i) const noexcept {
        return edges_[i].to;
    }
    inline const Cap& residual(u32 i) const noexcept {
        return edges_[i].residual;
    }
    inline const Cost& cost(u32 i) const noexcept {
        return edges_[i].cost;
    }
    inline const Cap& flowed(u32 i) const noexcept {
        assert(i < edgeSize());
        return residual(i ^ 1);
    }
    inline const Cap& capcacity(u32 i) const noexcept {
        assert(i < edgeSize());
        return residual(i) + flowed(i);
    }

    inline Cost edgeCost(u32 i) const noexcept {
        return potential_[from(i)] + cost(i) - potential_[to(i)];
    }
    bool relax(u32 i) {
        if (residual(i) > 0 and dist_[to(i)] > dist_[from(i)] + edgeCost(i)) {
            dist_[to(i)] = dist_[from(i)] + edgeCost(i);
            prev_[to(i)] = U32Pair{from(i), i};
            return true;
        }
        return false;
    }

    bool dijkstra(u32 s, u32 t) {
        assert(s < size());
        assert(t < size());
        std::fill(dist_.begin(), dist_.end(), invalid);
        dist_[s] = Cost{};
        using qt = std::pair<Cost, u32>;
        std::priority_queue<qt, std::vector<qt>, std::greater<qt>> que;
        que.emplace(dist_[s], s);
        while (que.size()) {
            auto [d, v]{que.top()};
            que.pop();
            if (dist_[v] < d) continue;
            for (u32 i : g_[v]) {
                if (relax(i)) {
                    que.emplace(dist_[to(i)], to(i));
                }
            }
        }
        return dist_[t] < invalid;
    }

    bool dagdp(u32 s, u32 t) {
        assert(s < size());
        assert(t < size());
        std::fill(dist_.begin(), dist_.end(), invalid);
        dist_[s] = Cost{};
        std::vector<u32> ord;
        ord.reserve(size());
        std::vector<bool> vis(size());
        vis[s] = true;
        ord.push_back(s);
        for (u32 i{} ; i < ord.size() ; i++) {
            u32 v{ord[i]};
            for (auto e : g_[v]) {
                if (!vis[to(e)]) {
                    ord.push_back(to(e));
                    vis[to(e)] = true;
                }
                relax(e);
            }
        }
        return dist_[t] < invalid;
    }

    bool bellmanford(u32 s, u32 t) {
        assert(s < size());
        assert(t < size());
        std::fill(dist_.begin(), dist_.end(), invalid);
        dist_[s] = Cost{};
        for (u32 i{} ; i + 1 < size() ; i++) {
            for (u32 j{} ; j < edgeSize() ; j++) {
                relax(j);
            }
        }
        return dist_[t] < invalid;
    }

    void updatePotential() {
        for (u32 v{} ; v < size() ; v++) {
            potential_[v] = potential_[v] + (dist_[v] == invalid ? Cost{} : dist_[v]);
        }
    }

    Cap flush(u32 s, u32 t, Cap up = std::numeric_limits<Cap>::max()) {
        for (u32 v{t} ; v != s ; v = prev_[v].first()) {
            up = std::min(up, residual(prev_[v].second()));
        }
        for (u32 v{t} ; v != s ; v = prev_[v].first()) {
            edges_[prev_[v].second()].residual -= up;
            edges_[prev_[v].second() ^ 1].residual += up;
        }
        return up;
    }

    bool flow(u32 s, u32 t, Cap flow) {
        assert(s < size());
        assert(t < size());
        while (flow > 0 and dijkstra(s, t)) {
            updatePotential();
            Cap water{flush(s, t, flow)};
            mcf_ += potential_[t] * water;
            flow -= water;
        }
        return flow == 0;
    }

    Cap maxflow(u32 s, u32 t) {
        assert(s < size());
        assert(t < size());
        Cap flow{};
        while (dijkstra(s, t)) {
            updatePotential();
            Cap water{flush(s, t)};
            mcf_ += potential_[t] * water;
            flow += water;
        }
        return flow;
    }

    std::vector<Cost> slope(u32 s, u32 t) {
        assert(s < size());
        assert(t < size());
        Cap flow{};
        std::vector<Cost> res;
        while (dijkstra(s, t)) {
            updatePotential();
            Cap water{flush(s, t)};
            for (u32 i{} ; i < water ; i++) {
                res.emplace_back(mcf_);
                mcf_ += potential_[t];
                flow++;
            }
        }
        res.emplace_back(mcf_);
        return res;
    }

    struct Line {
        Cap dn{}, up{};
        Cost slope{};
        Line() = default;
        Line(Cap dn, Cap up, Cost slope) : dn{dn}, up{up}, slope{slope} {}
    };
    std::vector<Line> slopeACL(u32 s, u32 t) {
        assert(s < size());
        assert(t < size()); 
        Cap flow{};
        std::vector<Line> res;
        while (dijkstra(s, t)) {
            updatePotential();
            Cap water{flush(s, t)};
            mcf_ += potential_[t] * water;
            res.emplace_back(flow, flow + water, potential_[t]);
            flow += water;
        }
        return res;
    }

    Cost minCost() const noexcept {
        return mcf_;
    }
};

} // namespace zawa

SuccessiveShortestPath (最小費用流、最短路反復法) (Src/Graph/Flow/SuccessiveShortestPath.hpp)

概要

ライブラリの使い方

テンプレート引数

コンストラクタ

addEdge

reverseId

from

to

residual

cost

flowed

capacity

flow

maxflow

minCost

slope

slopeACL

edgeCost

relax

dijkstra

dagdp

bellmanford

updatePotential

flush

emplace

Depends on

Verified with

Code

SuccessiveShortestPath (最小費用流、最短路反復法)
(Src/Graph/Flow/SuccessiveShortestPath.hpp)