最小費用流

graph/min_cost_flow.hpp · GitHub で見る

共通ページ

#include "lib/graph/min_cost_flow.hpp"

MinCostFlow<> graph(3);
graph.add_edge(0, 1, 5, -2);
graph.add_edge(1, 2, 3, 7);
auto [flow, cost] = graph.flow(0, 2, 3);

操作・型	内容・返り値	計算量
MinCostFlow<Cap,Cost,Total>(n=0)	n 頂点で構築。既定は Cap = std::int64_t、Cost = std::int64_t、Total = __int128_t	\(O(N)\)
size() / edge_count()	頂点数／追加した辺数	\(O(1)\)
reserve(m)	元の辺 m 本分の領域を確保	\(O(M)\)
add_edge(from,to,cap,cost)	容量と単位費用を指定して有向辺を追加し、0 始まりの辺番号を返す	償却 \(O(1)\)
flow(s,t,limit) / flow(s,t)	limit 以下で可能な最大の流量と、その最小費用を pair<Cap,Total> で返す。省略時の上限は Cap の最大値	初期化＋ \(O(N+Mw+M(N+M\log N)\log U)\)
slope(s,t,limit) / slope(s,t)	流量と最小費用の折れ線の頂点を vector<pair<Cap,Total>> で返す。最初は (0,0)、傾きが同じ区間は統合する	初期化＋ \(O(A(N+Mw))\)
get_edge(i)	追加時の順番で辺の情報を返す	\(O(1)\)
edges()	全ての辺の情報を vector<Edge> で返す	\(O(M)\)
reset_flow()	全ての流量を 0 に戻す。頂点と辺は残す	\(O(M)\)
potentials()	残余辺の reduced cost が非負となるポテンシャルを const vector<Cost>& で返す	初期化が必要ならその時間、準備済み: \(O(1)\)
Edge	from・to・cap・flow・cost を持つ	—（型）
total_type	費用合計の型 Total	—（型）

頂点は 0 始まりです。自己辺・多重辺・負の費用に対応します。各呼び出し前の残余グラフ全体に負閉路がないことが必要です。負の辺がある場合も、初期ポテンシャルを求めてから処理します。負閉路を含む一般の循環流は、検証用の最小費用 b-flow のように、負費用辺を上限まで流した状態から逆向きの非負費用辺へ変換して扱えます。

Cap は 64 bit 以下の標準整数型で bool を除き、容量・上限は非負です。Cost は 64 bit 以下の標準の符号付き整数型、Total は符号付き整数型です。cost に Cost の最小値は指定できません。ポテンシャル・最短距離・reduced cost とそれらの計算途中は Cost に、費用合計と積は Total に収めてください。既定の合計値は FastOutput で出力できます。

s != t が必要です。頂点数は非負、元の辺数は INT_MAX/2 以下にします。繰り返し呼ぶと現在の残余グラフから追加分を流し、返り値の流量・費用は今回の追加分です。要求した量を全て流せたかは flow == limit で確認します。辺の追加後も続けられますが、負閉路が生じないことを確認してください。

逐次最短路法では、ポテンシャル付き Dijkstra と Radix Heap を使い、reduced cost が 0 の辺はヒープを通さず処理します。flow は増加回数が多い場合に Network Simplex へ切り替え、作業量の上限に達した場合は容量スケーリングを使います。slope は全ての折れ点を求めるため、逐次最短路法を使い続けます。ポテンシャルの参照は次の変更操作まで有効です。

\(N\) 頂点・\(M\) 辺、整数幅 \(w\leq64\)、増加回数 \(A\)、容量・流量上限の最大値を \(U\geq1\) とします。初期化は非負費用からなら \(O(N+M)\)、一般には \(O(NM)\)。slope は初期化に加えて \(O(A(N+Mw))\)、flow は一定回数の逐次最短路法に加えて、切り替え時に \(O(M(N+M\log N)\log U)\) 時間です。potentials は必要なら初期化を行い、それ以外は \(O(1)\)。辺追加は償却 \(O(1)\)、get_edge・値取得は \(O(1)\)、edges・reset_flow・領域を増やす reserve は \(O(M)\)、辺追加後の隣接配列の構築は \(O(N+M)\)、空間は \(O(N+M)\)、slope の出力は追加で \(O(A)\) です。

依存ライブラリ

• Radix Heap
• D-ary Heap
• フローの残余グラフ
• 最小費用流の容量スケーリング
• 最小費用流の Network Simplex

検証コード

• test/min_cost_b_flow.test.cpp
• test/min_cost_flow.test.cpp

コード

#ifndef MIN_COST_FLOW_HPP
#define MIN_COST_FLOW_HPP

#include "../components/flow_graph.hpp"
#include "../components/min_cost_circulation.hpp"
#include "../components/min_cost_simplex.hpp"
#include "../components/RadixHeap.hpp"
#include <cstdint>
#include <type_traits>
#include <utility>

template<class Cap = std::int64_t, class Cost = std::int64_t, class Total = __int128_t>
class MinCostFlow {
	static_assert(std::is_integral_v<Cap> && !std::is_same_v<Cap, bool> && sizeof(Cap) <= 8);
	static_assert(std::is_integral_v<Cost> && std::is_signed_v<Cost> && sizeof(Cost) <= 8);
	static_assert(std::numeric_limits<Total>::is_integer && std::numeric_limits<Total>::is_signed);
	flow_detail::Graph<Cap> graph_;
	std::vector<Cost> cost_, potential_;
	bool valid_ = false;
	void vertex([[maybe_unused]] int v) const { assert(0 <= v && v < graph_.n); }
	int index(int i) const { assert(0 <= i && std::size_t(i) < graph_.arcs.size() / 2); return 2 * i; }
	void initialize() {
		if (valid_) return;
		potential_.assign(graph_.n, Cost(0));
		bool changed = false;
		for (int round = 0; round < graph_.n; ++round) {
			changed = false;
			for (int i = 0; i < int(cost_.size()); ++i) {
				if (!graph_.arcs[i].cap) continue;
				int from = graph_.arcs[i ^ 1].to, to = graph_.arcs[i].to;
				Cost next = potential_[from] + cost_[i];
				if (next < potential_[to]) potential_[to] = next, changed = true;
			}
			if (!changed) break;
		}
		assert(!changed);
		valid_ = true;
	}
	template<class Report> std::pair<Cap, Total> run(int s, int t, Cap limit, Report report, int rounds = -1) {
		vertex(s), vertex(t);
		assert(s != t && limit >= 0);
		graph_.build();
		initialize();
		auto& g = graph_;
		const Cost inf = std::numeric_limits<Cost>::max();
		std::vector<Cost> distance(g.n);
		std::vector<int> parent(g.n), zero, settled;
		std::vector<unsigned char> visited(g.n);
		zero.reserve(g.n), settled.reserve(g.n);
		RadixHeap<Cost> heap;
		Cap sent = 0;
		Total total = 0;
		while (sent < limit && rounds != 0) {
			std::fill(distance.begin(), distance.end(), inf);
			std::fill(visited.begin(), visited.end(), 0);
			zero.clear(), settled.clear(), heap.clear();
			distance[s] = 0;
			zero.push_back(s);
			while (!zero.empty() || !heap.empty()) {
				int v;
				if (!zero.empty()) v = zero.back(), zero.pop_back();
				else {
					auto [d, u] = heap.pop_value();
					if (d != distance[u]) continue;
					v = u;
				}
				if (visited[v]) continue;
				visited[v] = 1;
				settled.push_back(v);
				if (v == t) break;
				for (int j = g.start[v]; j < g.start[v + 1]; ++j) {
					int i = g.order[j];
					const auto& e = g.arcs[i];
					if (!e.cap || visited[e.to]) continue;
					Cost reduced = cost_[i] + potential_[v] - potential_[e.to];
					assert(reduced >= 0);
					if (reduced >= inf - distance[v]) continue;
					Cost next = distance[v] + reduced;
					if (next >= distance[e.to]) continue;
					distance[e.to] = next, parent[e.to] = i;
					if (!reduced) zero.push_back(e.to);
					else heap.push({next, e.to});
				}
			}
			if (!visited[t]) break;
			for (int v : settled) potential_[v] -= distance[t] - distance[v];
			Cap amount = limit - sent;
			for (int v = t; v != s; v = g.arcs[parent[v] ^ 1].to) amount = std::min(amount, g.arcs[parent[v]].cap);
			for (int v = t; v != s; v = g.arcs[parent[v] ^ 1].to) {
				int i = parent[v];
				g.arcs[i].cap -= amount, g.arcs[i ^ 1].cap += amount;
			}
			Cost price = potential_[t] - potential_[s];
			sent += amount, total += Total(amount) * Total(price);
			report(sent, total, price);
			if (rounds > 0) --rounds;
		}
		return {sent, total};
	}

  public:
	using total_type = Total;
	struct Edge { int from, to; Cap cap, flow; Cost cost; };
	explicit MinCostFlow(int n = 0) : graph_(n) {}
	int size() const { return graph_.n; }
	int edge_count() const { return int(cost_.size() / 2); }
	void reserve(std::size_t m) { graph_.reserve(m); cost_.reserve(2 * m); }
	int add_edge(int from, int to, Cap cap, Cost cost) {
		assert(cost > std::numeric_limits<Cost>::min());
		int id = graph_.add_edge(from, to, cap);
		cost_.push_back(cost), cost_.push_back(-cost);
		valid_ = false;
		return id;
	}
	Edge get_edge(int i) const {
		int j = index(i);
		const auto& a = graph_.arcs[j];
		const auto& b = graph_.arcs[j ^ 1];
		return {b.to, a.to, Cap(a.cap + b.cap), b.cap, cost_[j]};
	}
	std::vector<Edge> edges() const {
		std::vector<Edge> result;
		result.reserve(edge_count());
		for (int i = 0; i < edge_count(); ++i) result.push_back(get_edge(i));
		return result;
	}
	void reset_flow() { graph_.reset_flow(); valid_ = false; }
	const std::vector<Cost>& potentials() { initialize(); return potential_; }
	std::pair<Cap, Total> flow(int s, int t, Cap limit = std::numeric_limits<Cap>::max()) {
		int count = 0, budget = 8;
		auto result = run(s, t, limit, [&](Cap, Total, Cost) { ++count; }, budget);
		if (count == budget && result.first < limit) {
			auto extra = flow_detail::simplex_flow<Cap, Cost, Total>(graph_, cost_, s, t, limit - result.first);
			if (!extra) extra = flow_detail::scaling_flow<Cap, Cost, Total>(graph_, cost_, s, t, limit - result.first);
			result.first += extra->first, result.second += extra->second;
			valid_ = false;
		}
		return result;
	}
	std::vector<std::pair<Cap, Total>> slope(int s, int t, Cap limit = std::numeric_limits<Cap>::max()) {
		std::vector<std::pair<Cap, Total>> result{{Cap(0), Total(0)}};
		Cost previous = 0;
		bool first = true;
		run(s, t, limit, [&](Cap sent, Total total, Cost price) {
			if (!first && previous == price) result.pop_back();
			result.emplace_back(sent, total);
			previous = price, first = false;
		});
		return result;
	}
};

#endif

Back to top page