Minimum Cost Flow
(flow/min_cost_flow.hpp)

• View this file on GitHub
• Last update: 2024-01-08 01:08:59+09:00
• Include: #include "flow/min_cost_flow.hpp"

Description

フローネットワークの最小費用流を求める．この実装では，primal-dual 法を用いている．

Operations

• MinCostFlow(int V)
  • グラフを $V$ 頂点で初期化する
  • 時間計算量: $O(V)$
• void add_edge(int u, int v, Cap cap, Cost cost)
  • 容量 $cap$，コスト $cost$ の辺 $(u, v)$ を追加する
  • 時間計算量: $O(1)$
• void add_edge(int u, int v, Cap lb, Cap ub, Cost cost)
  • 最小流量 $lb$, 容量 $ub$，コスト $cost$ の辺 $(u, v)$ を追加する
  • 時間計算量: $O(1)$
• Cost min_cost_flow(int s, int t, Cap f, bool arbitrary)
  • 始点 $s$ から終点 $t$ への流量 $f$ の最小費用流を求める．arbitrary == true の場合，流量は $f$ 以下の任意の値とする．
  • 時間計算量: $O(fE\log V)$

Note

このライブラリがそのまま使える場合は，すべての辺のコストが非負である普通の最小費用流のとき．以下，いろいろな状況での使い方を説明する．

• 負辺がある場合
  • ポテンシャルの初期値の計算に，負辺があっても動作する最短路アルゴリズムを用いる必要がある．Bellman-Ford algorithm を用いることができる．また，グラフがDAGである場合は，トポロジカルソートしてDPすることができる．calculate_initial_potential()という private メソッドを用意しているのでその中を自分で書き換える．
  • 蟻本に載っているテク（超頂点を作って頑張る）を用いて負辺を除去することもできる．
• 負閉路がある場合
  • Bellman-Ford algorithm で負閉路を見つけてそこに流せるだけ流しておけば良い．

Verified with

• test/aoj/GRL_6_B.test.cpp

Code

#pragma once
#include <algorithm>
#include <functional>
#include <limits>
#include <queue>
#include <utility>
#include <vector>

template <typename Cap, typename Cost>
class MinCostFlow {
   public:
    MinCostFlow() = default;
    explicit MinCostFlow(int V) : V(V), G(V), add(0) {}

    void add_edge(int u, int v, Cap cap, Cost cost) {
        G[u].emplace_back(v, cap, cost, (int)G[v].size());
        G[v].emplace_back(u, 0, -cost, (int)G[u].size() - 1);
    }

    void add_edge(int u, int v, Cap lb, Cap ub, Cost cost) {
        add_edge(u, v, ub - lb, cost);
        add_edge(u, v, lb, cost - M);
        add += M * lb;
    }

    Cost min_cost_flow(int s, int t, Cap f, bool arbitrary = false) {
        Cost ret = add;
        std::vector<Cost> dist(V);
        std::vector<int> prevv(V), preve(V);
        using P = std::pair<Cost, int>;
        std::priority_queue<P, std::vector<P>, std::greater<P>> pq;

        auto h = calculate_initial_potential(s);

        while (f > 0) {
            // update h using dijkstra
            std::ranges::fill(dist, INF);
            dist[s] = 0;
            pq.emplace(0, s);
            while (!pq.empty()) {
                auto [d, v] = pq.top();
                pq.pop();
                if (dist[v] < d) continue;
                for (int i = 0; i < (int)G[v].size(); ++i) {
                    Edge& e = G[v][i];
                    Cost ndist = dist[v] + e.cost + h[v] - h[e.to];
                    if (e.cap > 0 && dist[e.to] > ndist) {
                        dist[e.to] = ndist;
                        prevv[e.to] = v;
                        preve[e.to] = i;
                        pq.emplace(dist[e.to], e.to);
                    }
                }
            }

            if (!arbitrary && dist[t] == INF) return -1;
            for (int v = 0; v < V; ++v) h[v] += dist[v];

            if (arbitrary && h[t] >= 0) break;

            Cap d = f;
            for (int v = t; v != s; v = prevv[v]) {
                d = std::min(d, G[prevv[v]][preve[v]].cap);
            }
            f -= d;
            ret += d * h[t];
            for (int v = t; v != s; v = prevv[v]) {
                Edge& e = G[prevv[v]][preve[v]];
                e.cap -= d;
                G[v][e.rev].cap += d;
            }
        }
        return ret;
    }

   private:
    struct Edge {
        int to;
        Cap cap;
        Cost cost;
        int rev;
        Edge(int to, Cap cap, Cost cost, int rev)
            : to(to), cap(cap), cost(cost), rev(rev) {}
    };

    static constexpr Cost INF = std::numeric_limits<Cost>::max() / 2;
    // large constant used for minimum flow requirement for edges
    static constexpr Cost M = INF / 1e9;

    int V;
    std::vector<std::vector<Edge>> G;
    Cost add;

    std::vector<Cost> calculate_initial_potential(int s) {
        std::vector<Cost> h(V);
        // if all costs are nonnegative, then do nothing
        return h;

        // if there is a negative edge,
        // use Bellman-Ford or topological sort and a DP (for DAG)
        // std::fill(h.begin(), h.end(), INF);
        // h[s] = 0;
        // for (int i = 0; i < V - 1; ++i) {
        //     for (int v = 0; v < V; ++v) {
        //         for (auto& e : G[v]) {
        //             if (e.cap > 0 && h[v] != INF && h[e.to] > h[v] + e.cost)
        //             {
        //                 h[e.to] = h[v] + e.cost;
        //             }
        //         }
        //     }
        // }

        // return h;
    }
};

#line 2 "flow/min_cost_flow.hpp"
#include <algorithm>
#include <functional>
#include <limits>
#include <queue>
#include <utility>
#include <vector>

template <typename Cap, typename Cost>
class MinCostFlow {
   public:
    MinCostFlow() = default;
    explicit MinCostFlow(int V) : V(V), G(V), add(0) {}

    void add_edge(int u, int v, Cap cap, Cost cost) {
        G[u].emplace_back(v, cap, cost, (int)G[v].size());
        G[v].emplace_back(u, 0, -cost, (int)G[u].size() - 1);
    }

    void add_edge(int u, int v, Cap lb, Cap ub, Cost cost) {
        add_edge(u, v, ub - lb, cost);
        add_edge(u, v, lb, cost - M);
        add += M * lb;
    }

    Cost min_cost_flow(int s, int t, Cap f, bool arbitrary = false) {
        Cost ret = add;
        std::vector<Cost> dist(V);
        std::vector<int> prevv(V), preve(V);
        using P = std::pair<Cost, int>;
        std::priority_queue<P, std::vector<P>, std::greater<P>> pq;

        auto h = calculate_initial_potential(s);

        while (f > 0) {
            // update h using dijkstra
            std::ranges::fill(dist, INF);
            dist[s] = 0;
            pq.emplace(0, s);
            while (!pq.empty()) {
                auto [d, v] = pq.top();
                pq.pop();
                if (dist[v] < d) continue;
                for (int i = 0; i < (int)G[v].size(); ++i) {
                    Edge& e = G[v][i];
                    Cost ndist = dist[v] + e.cost + h[v] - h[e.to];
                    if (e.cap > 0 && dist[e.to] > ndist) {
                        dist[e.to] = ndist;
                        prevv[e.to] = v;
                        preve[e.to] = i;
                        pq.emplace(dist[e.to], e.to);
                    }
                }
            }

            if (!arbitrary && dist[t] == INF) return -1;
            for (int v = 0; v < V; ++v) h[v] += dist[v];

            if (arbitrary && h[t] >= 0) break;

            Cap d = f;
            for (int v = t; v != s; v = prevv[v]) {
                d = std::min(d, G[prevv[v]][preve[v]].cap);
            }
            f -= d;
            ret += d * h[t];
            for (int v = t; v != s; v = prevv[v]) {
                Edge& e = G[prevv[v]][preve[v]];
                e.cap -= d;
                G[v][e.rev].cap += d;
            }
        }
        return ret;
    }

   private:
    struct Edge {
        int to;
        Cap cap;
        Cost cost;
        int rev;
        Edge(int to, Cap cap, Cost cost, int rev)
            : to(to), cap(cap), cost(cost), rev(rev) {}
    };

    static constexpr Cost INF = std::numeric_limits<Cost>::max() / 2;
    // large constant used for minimum flow requirement for edges
    static constexpr Cost M = INF / 1e9;

    int V;
    std::vector<std::vector<Edge>> G;
    Cost add;

    std::vector<Cost> calculate_initial_potential(int s) {
        std::vector<Cost> h(V);
        // if all costs are nonnegative, then do nothing
        return h;

        // if there is a negative edge,
        // use Bellman-Ford or topological sort and a DP (for DAG)
        // std::fill(h.begin(), h.end(), INF);
        // h[s] = 0;
        // for (int i = 0; i < V - 1; ++i) {
        //     for (int v = 0; v < V; ++v) {
        //         for (auto& e : G[v]) {
        //             if (e.cap > 0 && h[v] != INF && h[e.to] > h[v] + e.cost)
        //             {
        //                 h[e.to] = h[v] + e.cost;
        //             }
        //         }
        //     }
        // }

        // return h;
    }
};

Back to top page