Min Cost Flow
(graph/flow/min_cost_flow.hpp)
- View this file on GitHub
- Last update: 2026-08-04 16:49:58+09:00
- Include:
#include "graph/flow/min_cost_flow.hpp"
Overview
MinCostFlow<Cap, Cost> sends flow from a source s to a sink t while
minimizing total cost. Each edge has a capacity and a cost per unit of flow.
flow uses an adaptive one-shot solver. Small flows and graphs whose terminal
capacity can be carried by only a few arcs use successive shortest augmenting
paths. When a non-negative-cost instance necessarily uses many terminal arcs,
flow instead solves an exact residual flow with network simplex. The simplex
path has a pivot threshold and falls back to a polynomial capacity-scaling
solver. If the terminal-capacity upper bound is infeasible because of an
internal bottleneck, a maximum-flow pass finds the exact sendable value first.
When that upper bound saturates both terminal cuts, those forced terminal arcs
are contracted directly into vertex balances before the exact solve.
slope always uses successive shortest augmenting paths because it must retain
every flow-cost breakpoint. Its shortest-path implementation uses potentials,
early sink termination, and a radix heap for fixed-width integer costs. A
binary heap is used for other cost types. Initial Bellman-Ford work is skipped
on the first call when all added costs are non-negative.
Graph Orientation
Directed flow network. An edge added by add_edge(from, to, cap, cost) can
send flow only from from to to. For an undirected capacity, add both
directions with the desired costs.
The graph is stateful. Running flow or slope changes residual capacities and
stores the chosen flow. Use get_edge or edges afterward to inspect the
result.
How to Use It
Create MinCostFlow<Cap, Cost> mcf(n), add directed edges with capacity and
cost, and call mcf.flow(s, t, flow_limit).
The returned pair is {sent_flow, minimum_cost}. If the requested
flow_limit cannot be fully sent, sent_flow will be smaller.
Use slope(s, t, flow_limit) when you need the minimum cost for every
breakpoint of the amount of flow. The returned vector starts with {0, 0} and
then adds one entry after each augmentation.
When the edge count is known, call reserve_edges(m) before adding edges. It
reserves edge metadata and adjacency capacity with average-degree headroom.
If every endpoint degree is known, reserve_edges(m, degrees) reserves exact
adjacency capacities. Each original edge contributes one to both endpoint
degrees, and a self-loop contributes two to its endpoint.
Edge Fields
| Field | Type | Meaning |
|---|---|---|
from |
int |
Original edge source. |
to |
int |
Original edge destination. |
cap |
Cap |
Original capacity currently assigned to this edge. |
flow |
Cap |
Flow currently sent through this edge. |
cost |
Cost |
Cost per unit of flow on this edge. |
Methods
| Method | Signature | Description | Complexity |
|---|---|---|---|
| Constructor | MinCostFlow() |
Creates an empty flow graph. | $O(1)$ |
| Constructor | explicit MinCostFlow(int n) |
Creates a graph with n vertices. |
$O(N)$ |
size |
int size() const |
Returns the number of vertices. | $O(1)$ |
edge_count |
int edge_count() const |
Returns the number of original edges. | $O(1)$ |
reserve_edges |
void reserve_edges(int edge_count) |
Reserves edge metadata and average-degree adjacency headroom. | $O(N + M)$ when reallocation occurs |
reserve_edges |
void reserve_edges(int edge_count, const std::vector<int>& degrees) |
Reserves edge metadata and exact residual adjacency capacities. | $O(N + M)$ when reallocation occurs |
add_edge |
int add_edge(int from, int to, Cap cap, Cost cost) |
Adds a directed edge and returns its edge id. | Amortized $O(1)$ |
get_edge |
Edge get_edge(int i) const |
Returns the current state of original edge i. |
$O(1)$ |
edges |
std::vector<Edge> edges() const |
Returns all original edges with current flow. | $O(M)$ |
flow |
std::pair<Cap, Cost> flow(int s, int t) |
Sends as much flow as possible with the adaptive one-shot solver. | See below |
flow |
std::pair<Cap, Cost> flow(int s, int t, Cap flow_limit) |
Sends at most flow_limit flow with the adaptive one-shot solver. |
See below |
slope |
std::vector<std::pair<Cap, Cost>> slope(int s, int t) |
Returns flow-cost breakpoints using successive shortest paths. | See below |
slope |
std::vector<std::pair<Cap, Cost>> slope(int s, int t, Cap flow_limit) |
Returns flow-cost breakpoints up to flow_limit. |
See below |
Time Complexity
Let $F$ be the number of shortest-path augmentations, $U$ the maximum residual capacity, and $W$ the number of bits in the integer distance key.
The successive-shortest-path path performs initial Bellman-Ford in $O(NM)$
when negative residual costs need it. Each binary-heap shortest path takes
$O(M\log N)$. For fixed-width integer Cost, the radix-heap bound is
$O((N+M)W)$ per augmentation instead. Thus slope, and flow when it selects
this path, take
with the binary heap, or $O(NM + F(N+M)W)$ with the radix heap. The $O(NM)$ term is skipped on a fresh graph whose original costs are all non-negative.
The network-simplex path is guarded by at most a constant number of simplex attempts, a polynomial capacity-scaling fallback, and, when necessary, one maximum-flow computation. Its conservative worst-case bound is
\[O\left( N^2M + (N+M)^2 + M\log U\,(M+N\log N) \right).\]If the requested value or the terminal-capacity upper bound is feasible, the $O(N^2M)$ maximum-flow term is skipped. The dispatch examines only the two terminal adjacency lists: the simplex path is considered only for signed integer capacities, signed costs, non-negative original costs, at least 64 edges, and a requested value that requires at least eight arcs at both terminals. These are performance heuristics and do not affect correctness.
Notes
Costs may be negative, but the residual graph must not contain a reachable
negative-cost cycle. If such a cycle exists, the minimum cost is not
well-defined. Adding any negative-cost edge keeps flow on the
successive-shortest-path path.
Calling flow or slope a second time continues from the current residual
state and sends additional flow.
Example
#include "graph/flow/min_cost_flow.hpp"
#include <iostream>
int main() {
m1une::flow::MinCostFlow<long long, long long> mcf(4);
mcf.add_edge(0, 1, 2, 1);
mcf.add_edge(0, 2, 1, 2);
mcf.add_edge(1, 2, 1, 0);
mcf.add_edge(1, 3, 1, 3);
mcf.add_edge(2, 3, 2, 1);
auto [flow, cost] = mcf.flow(0, 3, 2);
std::cout << flow << " " << cost << "\n"; // 2 5
}
Depends on
Required by
Verified with
- verify/graph/cow_game.test.cpp
- verify/graph/flow/flow_algorithms.test.cpp
- verify/graph/flow/min_cost_flow.test.cpp
- verify/graph/graph_algorithms.test.cpp
- verify/graph/range_edge_graph.test.cpp
Code
#ifndef M1UNE_FLOW_MIN_COST_FLOW_HPP
#define M1UNE_FLOW_MIN_COST_FLOW_HPP 1
#include <algorithm>
#include <array>
#include <bit>
#include <cassert>
#include <cstddef>
#include <functional>
#include <limits>
#include <optional>
#include <type_traits>
#include <utility>
#include <vector>
#include "bounded_min_cost_flow.hpp"
#include "max_flow.hpp"
namespace m1une {
namespace flow {
template <class Cap, class Cost>
struct MinCostFlow {
struct Edge {
int from;
int to;
Cap cap;
Cap flow;
Cost cost;
};
private:
struct InternalEdge {
int to;
int rev;
Cap cap;
Cost cost;
};
int _n;
std::vector<std::pair<int, int>> _pos;
std::vector<std::vector<InternalEdge>> _g;
bool _has_negative_cost;
bool _has_flow;
template <class Key>
struct RadixHeap {
using Unsigned = std::make_unsigned_t<Key>;
static constexpr int bits = std::numeric_limits<Unsigned>::digits;
std::array<std::vector<std::pair<Unsigned, int>>, bits + 1> bucket;
Unsigned last = 0;
std::size_t count = 0;
static int index(Unsigned first, Unsigned second) {
return int(std::bit_width(first ^ second));
}
void clear() {
for (auto& values : bucket) values.clear();
last = 0;
count = 0;
}
bool empty() const {
return count == 0;
}
void push(Key key, int vertex) {
Unsigned value = static_cast<Unsigned>(key);
assert(last <= value);
bucket[index(value, last)].emplace_back(value, vertex);
count++;
}
std::pair<Key, int> pop() {
if (bucket[0].empty()) {
int i = 1;
while (bucket[i].empty()) i++;
last = bucket[i][0].first;
for (const auto& value : bucket[i]) {
last = std::min(last, value.first);
}
for (const auto& value : bucket[i]) {
bucket[index(value.first, last)].push_back(value);
}
bucket[i].clear();
}
auto [key, vertex] = bucket[0].back();
bucket[0].pop_back();
count--;
return {static_cast<Key>(key), vertex};
}
};
template <class Key>
struct BinaryHeap {
using Value = std::pair<Key, int>;
std::vector<Value> heap;
void clear() {
heap.clear();
}
bool empty() const {
return heap.empty();
}
void push(Key key, int vertex) {
heap.emplace_back(key, vertex);
std::push_heap(heap.begin(), heap.end(), std::greater<Value>());
}
Value pop() {
std::pop_heap(heap.begin(), heap.end(), std::greater<Value>());
Value result = heap.back();
heap.pop_back();
return result;
}
};
template <
class Key,
bool UseRadix =
std::numeric_limits<Key>::is_integer && sizeof(Key) <= 8
>
struct HeapSelector {
using Type = BinaryHeap<Key>;
};
template <class Key>
struct HeapSelector<Key, true> {
using Type = RadixHeap<Key>;
};
bool use_network_simplex(int s, int t, Cap flow_limit) const {
if (_has_negative_cost) return false;
if (_pos.size() < 64) return false;
auto add_saturated = [](Cap first, Cap second) {
const Cap maximum = std::numeric_limits<Cap>::max();
return maximum - first < second ? maximum : first + second;
};
struct TerminalCapacity {
Cap total = Cap(0);
std::array<Cap, 7> largest{};
};
auto add_capacity = [&](TerminalCapacity& terminal, Cap cap) {
terminal.total = add_saturated(terminal.total, cap);
for (Cap& current : terminal.largest) {
if (cap <= current) break;
std::swap(cap, current);
}
};
TerminalCapacity source;
for (const auto& e : _g[s]) {
if (e.to == s) continue;
add_capacity(source, e.cap);
}
TerminalCapacity sink;
for (const auto& e : _g[t]) {
if (e.to == t) continue;
Cap cap = _g[e.to][e.rev].cap;
add_capacity(sink, cap);
}
Cap target = std::min(
flow_limit,
std::min(source.total, sink.total)
);
if (target == Cap(0)) return false;
auto requires_eight_arcs = [&](const TerminalCapacity& terminal) {
Cap sum = Cap(0);
for (Cap cap : terminal.largest) {
sum = add_saturated(sum, cap);
}
return sum < target;
};
return requires_eight_arcs(source) && requires_eight_arcs(sink);
}
std::pair<Cap, Cost> network_simplex_flow(
int s,
int t,
Cap flow_limit
) {
struct ResidualArc {
int edge;
bool reverse;
};
using Solver = BoundedMinCostFlow<Cap, Cost, Cost>;
std::vector<ResidualArc> arcs;
arcs.reserve(2 * _pos.size());
for (int i = 0; i < int(_pos.size()); i++) {
auto [from, idx] = _pos[i];
const auto& e = _g[from][idx];
const auto& reverse = _g[e.to][e.rev];
if (e.cap != Cap(0)) {
arcs.push_back(ResidualArc{i, false});
}
if (reverse.cap != Cap(0)) {
arcs.push_back(ResidualArc{i, true});
}
}
auto add_saturated = [](Cap first, Cap second, bool& exact) {
const Cap maximum = std::numeric_limits<Cap>::max();
if (maximum - first < second) {
exact = false;
return maximum;
}
return first + second;
};
bool source_capacity_exact = true;
Cap source_capacity = Cap(0);
for (const auto& e : _g[s]) {
if (e.to == s) continue;
source_capacity = add_saturated(
source_capacity,
e.cap,
source_capacity_exact
);
}
bool sink_capacity_exact = true;
Cap sink_capacity = Cap(0);
for (const auto& e : _g[t]) {
if (e.to == t) continue;
sink_capacity = add_saturated(
sink_capacity,
_g[e.to][e.rev].cap,
sink_capacity_exact
);
}
Cap target = std::min(
flow_limit,
std::min(source_capacity, sink_capacity)
);
if (target == Cap(0)) return {Cap(0), Cost(0)};
struct ArcData {
int from;
int to;
Cap cap;
Cost cost;
};
auto arc_data = [&](const ResidualArc& arc) {
auto [from, idx] = _pos[arc.edge];
const auto& e = _g[from][idx];
const auto& reverse = _g[e.to][e.rev];
return arc.reverse
? ArcData{e.to, from, reverse.cap, reverse.cost}
: ArcData{from, e.to, e.cap, e.cost};
};
auto apply_flow = [&](const ResidualArc& arc, Cap amount) {
auto [from, idx] = _pos[arc.edge];
auto& e = _g[from][idx];
auto& reverse = _g[e.to][e.rev];
if (arc.reverse) {
reverse.cap -= amount;
e.cap += amount;
} else {
e.cap -= amount;
reverse.cap += amount;
}
};
bool target_infeasible = false;
if (
source_capacity_exact && sink_capacity_exact &&
target == source_capacity && target == sink_capacity
) {
Solver terminal_solver(_n);
terminal_solver.reserve_edges(int(arcs.size()));
std::vector<Cap> balance(_n, Cap(0));
std::vector<int> internal_arcs;
std::vector<int> fixed_arcs;
internal_arcs.reserve(arcs.size());
fixed_arcs.reserve(_g[s].size() + _g[t].size());
Cost fixed_cost = Cost(0);
for (int i = 0; i < int(arcs.size()); i++) {
ArcData data = arc_data(arcs[i]);
if (data.from == s) {
if (data.to == s) continue;
fixed_arcs.push_back(i);
fixed_cost += Cost(data.cap) * data.cost;
if (data.to != t) balance[data.to] += data.cap;
} else if (data.to == t) {
if (data.from == t) continue;
fixed_arcs.push_back(i);
fixed_cost += Cost(data.cap) * data.cost;
balance[data.from] -= data.cap;
} else if (data.to != s && data.from != t) {
terminal_solver.add_edge(
data.from,
data.to,
Cap(0),
data.cap,
data.cost
);
internal_arcs.push_back(i);
}
}
auto terminal_result = terminal_solver.min_cost_flow(balance);
if (terminal_result) {
for (int i : fixed_arcs) {
apply_flow(arcs[i], arc_data(arcs[i]).cap);
}
for (int i = 0; i < int(internal_arcs.size()); i++) {
apply_flow(
arcs[internal_arcs[i]],
terminal_result->flow(i)
);
}
_has_flow = true;
return {target, fixed_cost + terminal_result->cost};
}
target_infeasible = true;
}
Solver solver(_n);
solver.reserve_edges(int(arcs.size()));
for (const auto& arc : arcs) {
ArcData data = arc_data(arc);
solver.add_edge(
data.from,
data.to,
Cap(0),
data.cap,
data.cost
);
}
Cap sent = target;
std::optional<typename Solver::Result> result;
if (
!target_infeasible &&
target != std::numeric_limits<Cap>::max()
) {
result = solver.min_cost_st_flow(s, t, target);
}
if (!result) {
MaxFlow<Cap> feasible(_n);
feasible.reserve_edges(int(arcs.size()));
for (const auto& arc : arcs) {
auto [from, idx] = _pos[arc.edge];
const auto& e = _g[from][idx];
const auto& reverse = _g[e.to][e.rev];
if (arc.reverse) {
feasible.add_edge(e.to, from, reverse.cap);
} else {
feasible.add_edge(from, e.to, e.cap);
}
}
sent = feasible.max_flow(s, t, target);
if (sent == Cap(0)) return {Cap(0), Cost(0)};
result = solver.min_cost_st_flow(s, t, sent);
}
assert(result.has_value());
for (int i = 0; i < int(arcs.size()); i++) {
auto [from, idx] = _pos[arcs[i].edge];
auto& e = _g[from][idx];
auto& reverse = _g[e.to][e.rev];
Cap amount = result->flow(i);
if (arcs[i].reverse) {
reverse.cap -= amount;
e.cap += amount;
} else {
e.cap -= amount;
reverse.cap += amount;
}
}
_has_flow = true;
return {sent, result->cost};
}
void init_potential(int s, std::vector<Cost>& potential, Cost cost_inf) const {
if (!_has_negative_cost && !_has_flow) {
potential.assign(_n, Cost(0));
return;
}
potential.assign(_n, cost_inf);
potential[s] = Cost(0);
for (int iter = 0; iter < _n - 1; iter++) {
bool updated = false;
for (int v = 0; v < _n; v++) {
if (potential[v] == cost_inf) continue;
for (const auto& e : _g[v]) {
if (e.cap == Cap(0)) continue;
Cost nd = potential[v] + e.cost;
if (nd < potential[e.to]) {
potential[e.to] = nd;
updated = true;
}
}
}
if (!updated) break;
}
for (int v = 0; v < _n; v++) {
if (potential[v] == cost_inf) potential[v] = Cost(0);
}
}
public:
MinCostFlow() : MinCostFlow(0) {}
explicit MinCostFlow(int n)
: _n(n), _g(n), _has_negative_cost(false), _has_flow(false) {
assert(0 <= n);
}
int size() const {
return _n;
}
int edge_count() const {
return int(_pos.size());
}
void reserve_edges(int edge_count) {
assert(0 <= edge_count);
_pos.reserve(edge_count);
if (_n == 0 || edge_count == 0 ||
2 * std::size_t(edge_count) < std::size_t(_n)) {
return;
}
const std::size_t average_degree =
(3 * std::size_t(edge_count) + std::size_t(_n) - 1)
/ std::size_t(_n);
for (auto& edges : _g) edges.reserve(average_degree);
}
void reserve_edges(int edge_count, const std::vector<int>& degrees) {
assert(0 <= edge_count);
assert(int(degrees.size()) == _n);
_pos.reserve(edge_count);
for (int v = 0; v < _n; v++) {
assert(0 <= degrees[v]);
_g[v].reserve(degrees[v]);
}
}
int add_edge(int from, int to, Cap cap, Cost cost) {
assert(0 <= from && from < _n);
assert(0 <= to && to < _n);
assert(Cap(0) <= cap);
_has_negative_cost = _has_negative_cost || cost < Cost(0);
int id = int(_pos.size());
int from_id = int(_g[from].size());
int to_id = int(_g[to].size());
if (from == to) to_id++;
_pos.emplace_back(from, from_id);
_g[from].push_back(InternalEdge{to, to_id, cap, cost});
_g[to].push_back(InternalEdge{from, from_id, Cap(0), -cost});
return id;
}
Edge get_edge(int i) const {
assert(0 <= i && i < int(_pos.size()));
auto [from, idx] = _pos[i];
const auto& e = _g[from][idx];
const auto& re = _g[e.to][e.rev];
return Edge{from, e.to, e.cap + re.cap, re.cap, e.cost};
}
std::vector<Edge> edges() const {
std::vector<Edge> result;
result.reserve(_pos.size());
for (int i = 0; i < int(_pos.size()); i++) result.push_back(get_edge(i));
return result;
}
std::pair<Cap, Cost> flow(int s, int t) {
return flow(s, t, std::numeric_limits<Cap>::max());
}
std::pair<Cap, Cost> flow(int s, int t, Cap flow_limit) {
assert(0 <= s && s < _n);
assert(0 <= t && t < _n);
assert(s != t);
assert(Cap(0) <= flow_limit);
if (flow_limit == Cap(0)) return {Cap(0), Cost(0)};
if constexpr (
std::numeric_limits<Cap>::is_integer &&
std::numeric_limits<Cap>::is_signed &&
std::numeric_limits<Cost>::is_signed
) {
if (use_network_simplex(s, t, flow_limit)) {
return network_simplex_flow(s, t, flow_limit);
}
}
auto result = slope(s, t, flow_limit);
return result.back();
}
std::vector<std::pair<Cap, Cost>> slope(int s, int t) {
return slope(s, t, std::numeric_limits<Cap>::max());
}
std::vector<std::pair<Cap, Cost>> slope(int s, int t, Cap flow_limit) {
assert(0 <= s && s < _n);
assert(0 <= t && t < _n);
assert(s != t);
assert(Cap(0) <= flow_limit);
const Cost cost_inf = std::numeric_limits<Cost>::max() / Cost(4);
std::vector<Cost> potential, dist(_n);
std::vector<int> prev_v(_n), prev_e(_n);
std::vector<int> settled;
settled.reserve(_n);
typename HeapSelector<Cost>::Type que;
init_potential(s, potential, cost_inf);
std::vector<std::pair<Cap, Cost>> result;
result.emplace_back(Cap(0), Cost(0));
Cap flow = 0;
Cost cost = 0;
while (flow < flow_limit) {
std::fill(dist.begin(), dist.end(), cost_inf);
dist[s] = Cost(0);
settled.clear();
que.clear();
que.push(Cost(0), s);
while (!que.empty()) {
auto [d, v] = que.pop();
if (dist[v] != d) continue;
settled.push_back(v);
if (v == t) break;
for (int i = 0; i < int(_g[v].size()); i++) {
const auto& e = _g[v][i];
if (e.cap == Cap(0)) continue;
Cost nd = d + e.cost + potential[v] - potential[e.to];
if (nd >= dist[e.to]) continue;
dist[e.to] = nd;
prev_v[e.to] = v;
prev_e[e.to] = i;
que.push(nd, e.to);
}
}
if (dist[t] == cost_inf) break;
for (int v : settled) {
potential[v] += dist[v] - dist[t];
}
Cap add = flow_limit - flow;
for (int v = t; v != s; v = prev_v[v]) {
add = std::min(add, _g[prev_v[v]][prev_e[v]].cap);
}
Cost path_cost = potential[t] - potential[s];
for (int v = t; v != s; v = prev_v[v]) {
auto& e = _g[prev_v[v]][prev_e[v]];
e.cap -= add;
_g[e.to][e.rev].cap += add;
}
flow += add;
cost += Cost(add) * path_cost;
result.emplace_back(flow, cost);
}
_has_flow = _has_flow || flow != Cap(0);
return result;
}
};
} // namespace flow
} // namespace m1une
#endif // M1UNE_FLOW_MIN_COST_FLOW_HPP#line 1 "graph/flow/min_cost_flow.hpp"
#include <algorithm>
#include <array>
#include <bit>
#include <cassert>
#include <cstddef>
#include <functional>
#include <limits>
#include <optional>
#include <type_traits>
#include <utility>
#include <vector>
#line 1 "graph/flow/bounded_min_cost_flow.hpp"
#line 7 "graph/flow/bounded_min_cost_flow.hpp"
#include <cmath>
#line 11 "graph/flow/bounded_min_cost_flow.hpp"
#include <queue>
#line 14 "graph/flow/bounded_min_cost_flow.hpp"
namespace m1une {
namespace flow {
template <
class Cap,
class Cost,
class TotalCost = Cost,
std::size_t PivotLimitFactor = 8
>
struct BoundedMinCostFlow {
static_assert(std::numeric_limits<Cap>::is_integer);
static_assert(std::numeric_limits<Cap>::is_signed);
static_assert(std::numeric_limits<Cost>::is_specialized);
static_assert(std::numeric_limits<Cost>::is_signed);
struct Edge {
int from;
int to;
Cap lower;
Cap upper;
Cost cost;
};
struct ResultEdge {
int from;
int to;
Cap lower;
Cap upper;
Cap flow;
Cost cost;
};
struct Result {
std::vector<ResultEdge> edges;
std::vector<Cap> balance;
std::vector<Cost> potential;
TotalCost cost;
ResultEdge get_edge(int i) const {
assert(0 <= i && i < int(edges.size()));
return edges[i];
}
Cap flow(int i) const {
assert(0 <= i && i < int(edges.size()));
return edges[i].flow;
}
};
private:
struct NetworkEdge {
int to;
Cap cap;
Cost cost;
};
struct NetworkSimplexSolver {
enum class Status {
optimal,
infeasible,
pivot_limit_reached,
};
struct Parent {
int vertex;
int edge;
Cap up;
Cap down;
};
int n;
std::vector<NetworkEdge> edges;
std::vector<Cap> excess;
std::vector<Cost> potential;
std::size_t pivot_count = 0;
NetworkSimplexSolver(int vertex_count, const std::vector<Cap>& balance)
: n(vertex_count), excess(balance) {}
void reserve_edges(int edge_count) {
edges.reserve(2 * (edge_count + n));
}
int add_edge(int from, int to, Cap lower, Cap upper, Cost cost) {
int id = int(edges.size()) / 2;
edges.push_back(NetworkEdge{to, upper - lower, cost});
edges.push_back(NetworkEdge{from, Cap(0), -cost});
excess[from] -= lower;
excess[to] += lower;
return id;
}
Status solve(std::size_t pivot_limit) {
pivot_count = 0;
const int original_edge_count = int(edges.size());
potential.assign(n + 1, Cost(0));
Cost artificial_cost = Cost(1);
for (int edge = 0; edge < original_edge_count; edge += 2) {
artificial_cost += edges[edge].cost < Cost(0)
? -edges[edge].cost : edges[edge].cost;
}
std::vector<Parent> parent(n);
edges.reserve(original_edge_count + 2 * n);
for (int vertex = 0; vertex < n; vertex++) {
if (excess[vertex] >= Cap(0)) {
edges.push_back(NetworkEdge{n, Cap(0), artificial_cost});
edges.push_back(NetworkEdge{vertex, excess[vertex], -artificial_cost});
potential[vertex] = -artificial_cost;
} else {
edges.push_back(NetworkEdge{n, -excess[vertex], -artificial_cost});
edges.push_back(NetworkEdge{vertex, Cap(0), artificial_cost});
potential[vertex] = artificial_cost;
}
int edge = int(edges.size()) - 2;
parent[vertex] = Parent{
n, edge, edges[edge].cap, edges[edge ^ 1].cap
};
}
std::vector<int> depth(n + 1, 1);
depth[n] = 0;
std::vector<int> next(2 * (n + 1));
std::vector<int> previous(2 * (n + 1));
auto connect = [&](int first, int second) {
next[first] = second;
previous[second] = first;
};
for (int vertex = 0; vertex <= n; vertex++) {
connect(2 * vertex, 2 * vertex + 1);
}
for (int vertex = 0; vertex < n; vertex++) {
connect(2 * vertex + 1, next[2 * n]);
connect(2 * n, 2 * vertex);
}
auto push_flow = [&](int entering_edge) {
const int first = edges[entering_edge ^ 1].to;
const int second = edges[entering_edge].to;
const Cost cycle_cost =
edges[entering_edge].cost
+ potential[first] - potential[second];
Cap amount = edges[entering_edge].cap;
bool leave_first_side = true;
int leaving_vertex = second;
int first_ancestor = first;
int second_ancestor = second;
auto move_first_up = [&] {
if (parent[first_ancestor].down < amount) {
amount = parent[first_ancestor].down;
leaving_vertex = first_ancestor;
leave_first_side = true;
}
first_ancestor = parent[first_ancestor].vertex;
};
auto move_second_up = [&] {
if (parent[second_ancestor].up <= amount) {
amount = parent[second_ancestor].up;
leaving_vertex = second_ancestor;
leave_first_side = false;
}
second_ancestor = parent[second_ancestor].vertex;
};
if (depth[first_ancestor] >= depth[second_ancestor]) {
int difference = depth[first_ancestor] - depth[second_ancestor];
for (int i = 0; i < difference; i++) move_first_up();
} else {
int difference = depth[second_ancestor] - depth[first_ancestor];
for (int i = 0; i < difference; i++) move_second_up();
}
while (first_ancestor != second_ancestor) {
move_first_up();
move_second_up();
}
const int ancestor = first_ancestor;
if (amount != Cap(0)) {
int vertex = first;
while (vertex != ancestor) {
parent[vertex].up += amount;
parent[vertex].down -= amount;
vertex = parent[vertex].vertex;
}
vertex = second;
while (vertex != ancestor) {
parent[vertex].up -= amount;
parent[vertex].down += amount;
vertex = parent[vertex].vertex;
}
}
int vertex = first;
int new_parent = second;
std::pair<Cap, Cap> parent_capacities{
edges[entering_edge].cap - amount,
edges[entering_edge ^ 1].cap + amount
};
Cost potential_difference = -cycle_cost;
if (!leave_first_side) {
std::swap(vertex, new_parent);
std::swap(parent_capacities.first, parent_capacities.second);
potential_difference = -potential_difference;
}
int parent_edge = entering_edge ^ (leave_first_side ? 0 : 1);
while (new_parent != leaving_vertex) {
int new_depth = depth[new_parent];
int tour_index = 2 * vertex;
while (tour_index != 2 * vertex + 1) {
if ((tour_index & 1) == 0) {
new_depth++;
potential[tour_index / 2] += potential_difference;
depth[tour_index / 2] = new_depth;
} else {
new_depth--;
}
tour_index = next[tour_index];
}
connect(previous[2 * vertex], next[2 * vertex + 1]);
connect(2 * vertex + 1, next[2 * new_parent]);
connect(2 * new_parent, 2 * vertex);
std::swap(parent[vertex].edge, parent_edge);
parent_edge ^= 1;
std::swap(parent[vertex].up, parent_capacities.first);
std::swap(parent[vertex].down, parent_capacities.second);
std::swap(parent_capacities.first, parent_capacities.second);
int old_parent = parent[vertex].vertex;
parent[vertex].vertex = new_parent;
new_parent = vertex;
vertex = old_parent;
}
edges[parent_edge].cap = parent_capacities.first;
edges[parent_edge ^ 1].cap = parent_capacities.second;
};
bool pivot_limit_reached = false;
auto pivot = [&](int entering_edge) {
if (pivot_count == pivot_limit) {
pivot_limit_reached = true;
return false;
}
push_flow(entering_edge);
pivot_count++;
return true;
};
const int candidate_limit = std::max(
int(0.2 * std::sqrt(double(original_edge_count))), 10
);
const int minor_limit = std::max(candidate_limit / 10, 3);
std::vector<int> candidates;
candidates.reserve(candidate_limit);
auto minor_pivot = [&] {
Cost best_cost = Cost(0);
int best_edge = -1;
int index = 0;
while (index < int(candidates.size())) {
int edge = candidates[index];
if (edges[edge].cap == Cap(0)) {
candidates[index] = candidates.back();
candidates.pop_back();
continue;
}
Cost reduced_cost =
edges[edge].cost
+ potential[edges[edge ^ 1].to]
- potential[edges[edge].to];
if (reduced_cost >= Cost(0)) {
candidates[index] = candidates.back();
candidates.pop_back();
continue;
}
if (reduced_cost < best_cost) {
best_cost = reduced_cost;
best_edge = edge;
}
index++;
}
if (best_edge == -1) return false;
return pivot(best_edge);
};
int edge = 0;
while (true) {
for (int iteration = 0; iteration < minor_limit; iteration++) {
if (!minor_pivot()) break;
}
if (pivot_limit_reached) return Status::pivot_limit_reached;
Cost best_cost = Cost(0);
int best_edge = -1;
candidates.clear();
for (int scanned = 0; scanned < int(edges.size()); scanned++) {
if (edges[edge].cap != Cap(0)) {
Cost reduced_cost =
edges[edge].cost
+ potential[edges[edge ^ 1].to]
- potential[edges[edge].to];
if (reduced_cost < Cost(0)) {
if (reduced_cost < best_cost) {
best_cost = reduced_cost;
best_edge = edge;
}
candidates.push_back(edge);
if (int(candidates.size()) == candidate_limit) break;
}
}
edge++;
if (edge == int(edges.size())) edge = 0;
}
if (candidates.empty()) break;
if (!pivot(best_edge)) return Status::pivot_limit_reached;
}
for (int vertex = 0; vertex < n; vertex++) {
edges[parent[vertex].edge].cap = parent[vertex].up;
edges[parent[vertex].edge ^ 1].cap = parent[vertex].down;
}
bool feasible = true;
for (int vertex = 0; vertex < n; vertex++) {
int artificial_edge = original_edge_count + 2 * vertex;
if (
(excess[vertex] >= Cap(0)
&& edges[artificial_edge ^ 1].cap != Cap(0))
|| (excess[vertex] < Cap(0)
&& edges[artificial_edge].cap != Cap(0))
) {
feasible = false;
break;
}
}
potential.pop_back();
return feasible ? Status::optimal : Status::infeasible;
}
Cap edge_flow(int edge_id, Cap lower) const {
return lower + edges[2 * edge_id + 1].cap;
}
};
struct ScalingEdge {
int to;
int reverse;
Cap cap;
Cap flow;
Cost cost;
};
struct ScalingSolver {
int n;
std::vector<std::vector<ScalingEdge>> graph;
std::vector<std::pair<int, int>> positions;
std::vector<Cap> excess;
std::vector<Cost> potential;
std::vector<Cost> distance;
std::vector<int> parent_vertex;
std::vector<int> parent_edge;
std::vector<int> excess_vertices;
std::vector<int> deficit_vertices;
Cost farthest = Cost(0);
ScalingSolver(int vertex_count, const std::vector<Cap>& balance)
: n(vertex_count), graph(vertex_count), excess(balance),
potential(vertex_count, Cost(0)) {}
void reserve_edges(int edge_count) {
positions.reserve(edge_count);
}
int add_edge(int from, int to, Cap lower, Cap upper, Cost cost) {
int id = int(positions.size());
int from_edge = int(graph[from].size());
int to_edge = int(graph[to].size());
if (from == to) to_edge++;
positions.emplace_back(from, from_edge);
graph[from].push_back(ScalingEdge{
to, to_edge, upper, Cap(0), cost
});
graph[to].push_back(ScalingEdge{
from, from_edge, -lower, Cap(0), -cost
});
return id;
}
Cap residual_capacity(int from, int edge_id) const {
const auto& edge = graph[from][edge_id];
return edge.cap - edge.flow;
}
Cost residual_cost(int from, const ScalingEdge& edge) const {
return edge.cost + potential[from] - potential[edge.to];
}
void push(int from, int edge_id, Cap amount) {
auto& edge = graph[from][edge_id];
edge.flow += amount;
graph[edge.to][edge.reverse].flow -= amount;
}
void saturate_negative(Cap delta) {
excess_vertices.clear();
deficit_vertices.clear();
for (int from = 0; from < n; from++) {
for (
int edge_id = 0;
edge_id < int(graph[from].size());
edge_id++
) {
const auto& edge = graph[from][edge_id];
Cap residual = edge.cap - edge.flow;
residual -= residual % delta;
if (
residual_cost(from, edge) < Cost(0)
|| residual < Cap(0)
) {
int to = edge.to;
push(from, edge_id, residual);
excess[from] -= residual;
excess[to] += residual;
}
}
}
for (int vertex = 0; vertex < n; vertex++) {
if (excess[vertex] > Cap(0)) {
excess_vertices.push_back(vertex);
} else if (excess[vertex] < Cap(0)) {
deficit_vertices.push_back(vertex);
}
}
}
bool dual(Cap delta) {
excess_vertices.erase(
std::remove_if(
excess_vertices.begin(), excess_vertices.end(),
[&](int vertex) { return excess[vertex] < delta; }
),
excess_vertices.end()
);
deficit_vertices.erase(
std::remove_if(
deficit_vertices.begin(), deficit_vertices.end(),
[&](int vertex) { return excess[vertex] > -delta; }
),
deficit_vertices.end()
);
const Cost unreachable = std::numeric_limits<Cost>::max();
distance.assign(n, unreachable);
parent_vertex.assign(n, -1);
parent_edge.assign(n, -1);
using QueueEntry = std::pair<Cost, int>;
std::priority_queue<
QueueEntry,
std::vector<QueueEntry>,
std::greater<QueueEntry>
> queue;
for (int vertex : excess_vertices) {
distance[vertex] = Cost(0);
queue.emplace(Cost(0), vertex);
}
farthest = Cost(0);
int reached_deficits = 0;
while (!queue.empty()) {
auto [current_distance, from] = queue.top();
queue.pop();
if (distance[from] != current_distance) continue;
farthest = current_distance;
if (excess[from] <= -delta) reached_deficits++;
if (reached_deficits >= int(deficit_vertices.size())) break;
for (
int edge_id = 0;
edge_id < int(graph[from].size());
edge_id++
) {
const auto& edge = graph[from][edge_id];
if (edge.cap - edge.flow < delta) continue;
Cost next_distance =
current_distance + residual_cost(from, edge);
if (next_distance >= distance[edge.to]) continue;
distance[edge.to] = next_distance;
parent_vertex[edge.to] = from;
parent_edge[edge.to] = edge_id;
queue.emplace(next_distance, edge.to);
}
}
for (int vertex = 0; vertex < n; vertex++) {
potential[vertex] += std::min(distance[vertex], farthest);
}
return reached_deficits > 0;
}
void primal(Cap delta) {
for (int sink : deficit_vertices) {
if (distance[sink] > farthest) continue;
Cap amount = -excess[sink];
int root = sink;
while (parent_edge[root] != -1) {
int from = parent_vertex[root];
amount = std::min(
amount,
residual_capacity(from, parent_edge[root])
);
root = from;
}
amount = std::min(amount, excess[root]);
amount -= amount % delta;
if (amount <= Cap(0)) continue;
int vertex = sink;
while (parent_edge[vertex] != -1) {
int from = parent_vertex[vertex];
int edge_id = parent_edge[vertex];
push(from, edge_id, amount);
if (residual_capacity(from, edge_id) == Cap(0)) {
parent_edge[vertex] = -1;
}
vertex = from;
}
excess[sink] += amount;
excess[root] -= amount;
}
}
bool solve() {
Cap scale_bound = Cap(1);
for (Cap value : excess) {
scale_bound = std::max(scale_bound, value);
scale_bound = std::max(scale_bound, -value);
}
for (const auto& edges : graph) {
for (const auto& edge : edges) {
Cap residual = edge.cap - edge.flow;
scale_bound = std::max(scale_bound, residual);
scale_bound = std::max(scale_bound, -residual);
}
}
Cap delta = Cap(1);
while (delta <= scale_bound / Cap(2)) delta *= Cap(2);
while (true) {
saturate_negative(delta);
while (dual(delta)) primal(delta);
if (delta == Cap(1)) break;
delta /= Cap(2);
}
return excess_vertices.empty() && deficit_vertices.empty();
}
Cap edge_flow(int edge_id, Cap) const {
auto [from, index] = positions[edge_id];
return graph[from][index].flow;
}
};
int _n;
std::vector<Edge> _edges;
std::vector<Cap> _balance;
template <class Solver>
Result make_result(
const std::vector<Cap>& balance,
const Solver& solver,
std::vector<Cost> potential
) const {
Result result;
result.balance = balance;
result.cost = TotalCost(0);
result.edges.reserve(_edges.size());
for (int i = 0; i < int(_edges.size()); i++) {
const auto& edge = _edges[i];
Cap flow = solver.edge_flow(i, edge.lower);
result.cost += TotalCost(flow) * TotalCost(edge.cost);
result.edges.push_back(ResultEdge{
edge.from,
edge.to,
edge.lower,
edge.upper,
flow,
edge.cost
});
}
result.potential = std::move(potential);
return result;
}
std::vector<Cost> residual_potential(
const std::vector<ResultEdge>& edges
) const {
std::vector<Cost> potential(_n, Cost(0));
bool updated = false;
for (int iteration = 0; iteration < _n; iteration++) {
updated = false;
for (const ResultEdge& edge : edges) {
if (
edge.flow < edge.upper
&& potential[edge.to] > potential[edge.from] + edge.cost
) {
potential[edge.to] = potential[edge.from] + edge.cost;
updated = true;
}
if (
edge.lower < edge.flow
&& potential[edge.from] > potential[edge.to] - edge.cost
) {
potential[edge.from] = potential[edge.to] - edge.cost;
updated = true;
}
}
if (!updated) break;
}
assert(!updated);
return potential;
}
std::optional<Result> polynomial_min_cost_flow_impl(
const std::vector<Cap>& balance
) const {
ScalingSolver solver(_n, balance);
solver.reserve_edges(int(_edges.size()));
for (const auto& edge : _edges) {
solver.add_edge(
edge.from,
edge.to,
edge.lower,
edge.upper,
edge.cost
);
}
if (!solver.solve()) return std::nullopt;
Result result = make_result(balance, solver, {});
result.potential = residual_potential(result.edges);
return result;
}
public:
BoundedMinCostFlow() : BoundedMinCostFlow(0) {}
explicit BoundedMinCostFlow(int n) : _n(n), _balance(n, Cap(0)) {
assert(0 <= n);
}
int size() const {
return _n;
}
int edge_count() const {
return int(_edges.size());
}
void reserve_edges(int edge_count) {
assert(0 <= edge_count);
_edges.reserve(edge_count);
}
int add_edge(int from, int to, Cap lower, Cap upper, Cost cost) {
assert(0 <= from && from < _n);
assert(0 <= to && to < _n);
assert(lower <= upper);
int id = int(_edges.size());
_edges.push_back(Edge{from, to, lower, upper, cost});
return id;
}
Edge get_edge(int i) const {
assert(0 <= i && i < int(_edges.size()));
return _edges[i];
}
std::vector<Edge> edges() const {
return _edges;
}
void set_balance(int v, Cap b) {
assert(0 <= v && v < _n);
_balance[v] = b;
}
void add_balance(int v, Cap b) {
assert(0 <= v && v < _n);
_balance[v] += b;
}
void add_supply(int v, Cap supply) {
assert(Cap(0) <= supply);
add_balance(v, supply);
}
void add_demand(int v, Cap demand) {
assert(Cap(0) <= demand);
add_balance(v, -demand);
}
Cap balance(int v) const {
assert(0 <= v && v < _n);
return _balance[v];
}
const std::vector<Cap>& balances() const {
return _balance;
}
std::optional<Result> min_cost_flow() const {
return min_cost_flow(_balance);
}
std::optional<Result> min_cost_flow(const std::vector<Cap>& balance) const {
assert(int(balance.size()) == _n);
Cap balance_sum = Cap(0);
for (Cap value : balance) balance_sum += value;
if (balance_sum != Cap(0)) return std::nullopt;
NetworkSimplexSolver solver(_n, balance);
solver.reserve_edges(int(_edges.size()));
for (const auto& edge : _edges) {
solver.add_edge(edge.from, edge.to, edge.lower, edge.upper, edge.cost);
}
const std::size_t graph_size =
std::size_t(_n) + _edges.size() + 1;
std::size_t pivot_limit = 0;
if constexpr (PivotLimitFactor != 0) {
const std::size_t maximum =
std::numeric_limits<std::size_t>::max();
pivot_limit = graph_size > maximum / PivotLimitFactor
? maximum : PivotLimitFactor * graph_size;
}
auto status = solver.solve(pivot_limit);
if (status == NetworkSimplexSolver::Status::infeasible) {
return std::nullopt;
}
if (status == NetworkSimplexSolver::Status::pivot_limit_reached) {
return polynomial_min_cost_flow_impl(balance);
}
return make_result(balance, solver, std::move(solver.potential));
}
std::optional<Result> min_cost_flow_polynomial() const {
return min_cost_flow_polynomial(_balance);
}
std::optional<Result> min_cost_flow_polynomial(
const std::vector<Cap>& balance
) const {
assert(int(balance.size()) == _n);
Cap balance_sum = Cap(0);
for (Cap value : balance) balance_sum += value;
if (balance_sum != Cap(0)) return std::nullopt;
return polynomial_min_cost_flow_impl(balance);
}
std::optional<Result> min_cost_st_flow(int s, int t, Cap flow_value) const {
assert(0 <= s && s < _n);
assert(0 <= t && t < _n);
assert(s != t);
std::vector<Cap> balance = _balance;
balance[s] += flow_value;
balance[t] -= flow_value;
return min_cost_flow(balance);
}
std::optional<Result> min_cost_st_flow_polynomial(
int s,
int t,
Cap flow_value
) const {
assert(0 <= s && s < _n);
assert(0 <= t && t < _n);
assert(s != t);
std::vector<Cap> balance = _balance;
balance[s] += flow_value;
balance[t] -= flow_value;
return min_cost_flow_polynomial(balance);
}
};
template <
class Cap,
class Cost,
class TotalCost = Cost,
std::size_t PivotLimitFactor = 8
>
using BMinCostFlow = BoundedMinCostFlow<
Cap,
Cost,
TotalCost,
PivotLimitFactor
>;
} // namespace flow
} // namespace m1une
#line 1 "graph/flow/max_flow.hpp"
#line 9 "graph/flow/max_flow.hpp"
namespace m1une {
namespace flow {
template <class Cap>
struct MaxFlow {
struct Edge {
int from;
int to;
Cap cap;
Cap flow;
};
private:
struct InternalEdge {
int to;
int rev;
Cap cap;
};
struct Position {
int from;
int edge;
};
int _n;
std::vector<Position> _pos;
std::vector<std::vector<InternalEdge>> _g;
Cap highest_label_preflow_push(int s, int t) {
const int dead = 2 * _n;
const int unreachable = _n + 1;
std::vector<Cap> excess(_n, Cap(0));
std::vector<int> state(8 * std::size_t(_n) + 2);
int* height = state.data();
int* height_count = height + _n;
int* current = height_count + dead + 1;
int* queue = current + _n;
int* next = queue + _n;
int* bucket_head = next + _n;
std::vector<char> active(_n, false);
int highest = -1;
long long work = 0;
const long long arc_count =
2LL * static_cast<long long>(_pos.size());
const long long work_limit = std::max(1LL, 4 * arc_count + _n);
auto activate = [&](int v) {
if (v == s || v == t || active[v] || excess[v] == Cap(0) ||
height[v] >= dead) {
return;
}
active[v] = true;
next[v] = bucket_head[height[v]];
bucket_head[height[v]] = v;
highest = std::max(highest, height[v]);
};
auto rebuild_buckets = [&]() {
std::fill(bucket_head, bucket_head + dead + 1, -1);
std::fill(active.begin(), active.end(), false);
highest = -1;
for (int v = 0; v < _n; v++) activate(v);
};
auto global_relabel = [&]() {
std::fill(height, height + _n, unreachable);
std::fill(height_count, height_count + dead + 1, 0);
std::fill(current, current + _n, 0);
int head = 0;
int tail = 0;
height[t] = 0;
height[s] = _n;
queue[tail++] = t;
while (head != tail) {
int v = queue[head++];
for (const auto& e : _g[v]) {
if (e.to == s || height[e.to] != unreachable) continue;
const auto& reverse = _g[e.to][e.rev];
if (reverse.cap == Cap(0)) continue;
height[e.to] = height[v] + 1;
queue[tail++] = e.to;
}
}
for (int v = 0; v < _n; v++) height_count[height[v]]++;
rebuild_buckets();
work = 0;
};
auto gap = [&](int empty_height) {
for (int v = 0; v < _n; v++) {
if (v == s || v == t || height[v] <= empty_height ||
height[v] >= _n) {
continue;
}
height_count[height[v]]--;
height[v] = unreachable;
height_count[height[v]]++;
current[v] = 0;
}
rebuild_buckets();
};
auto relabel = [&](int v) -> bool {
int old_height = height[v];
int new_height = dead;
work += int(_g[v].size());
for (const auto& e : _g[v]) {
if (e.cap != Cap(0)) {
new_height = std::min(new_height, height[e.to] + 1);
}
}
height_count[old_height]--;
height[v] = std::min(new_height, dead);
height_count[height[v]]++;
current[v] = 0;
if (old_height < _n && height_count[old_height] == 0) {
gap(old_height);
return true;
}
return false;
};
auto push = [&](int v, InternalEdge& e) {
Cap sent = std::min(excess[v], e.cap);
bool was_zero = excess[e.to] == Cap(0);
e.cap -= sent;
_g[e.to][e.rev].cap += sent;
excess[v] -= sent;
excess[e.to] += sent;
if (was_zero) activate(e.to);
};
auto discharge = [&](int v) {
while (excess[v] != Cap(0) && height[v] < dead) {
if (current[v] == int(_g[v].size())) {
if (relabel(v)) return;
continue;
}
auto& e = _g[v][current[v]];
work++;
if (e.cap != Cap(0) && height[v] == height[e.to] + 1) {
push(v, e);
} else {
current[v]++;
}
}
activate(v);
};
for (auto& e : _g[s]) {
if (e.to == s || e.cap == Cap(0)) continue;
Cap sent = e.cap;
e.cap = Cap(0);
_g[e.to][e.rev].cap += sent;
excess[e.to] += sent;
}
global_relabel();
while (highest >= 0) {
if (bucket_head[highest] == -1) {
highest--;
continue;
}
int v = bucket_head[highest];
bucket_head[highest] = next[v];
if (!active[v] || height[v] != highest) continue;
active[v] = false;
discharge(v);
if (work >= work_limit) global_relabel();
}
return excess[t];
}
public:
MaxFlow() : MaxFlow(0) {}
explicit MaxFlow(int n) : _n(n), _g(n) {
assert(0 <= n);
}
int size() const {
return _n;
}
int edge_count() const {
return int(_pos.size());
}
void reserve_edges(int edge_count) {
assert(0 <= edge_count);
_pos.reserve(edge_count);
if (_n == 0 || edge_count == 0 ||
2 * std::size_t(edge_count) < std::size_t(_n)) {
return;
}
const std::size_t average_degree =
(3 * std::size_t(edge_count) + std::size_t(_n) - 1)
/ std::size_t(_n);
for (auto& edges : _g) edges.reserve(average_degree);
}
void reserve_edges(int edge_count, const std::vector<int>& degrees) {
assert(0 <= edge_count);
assert(int(degrees.size()) == _n);
_pos.reserve(edge_count);
for (int v = 0; v < _n; v++) {
assert(0 <= degrees[v]);
_g[v].reserve(degrees[v]);
}
}
int add_edge(int from, int to, Cap cap) {
assert(0 <= from && from < _n);
assert(0 <= to && to < _n);
assert(Cap(0) <= cap);
int id = int(_pos.size());
int from_id = int(_g[from].size());
int to_id = int(_g[to].size());
if (from == to) to_id++;
_pos.push_back(Position{from, from_id});
_g[from].push_back(InternalEdge{to, to_id, cap});
_g[to].push_back(InternalEdge{from, from_id, Cap(0)});
return id;
}
int add_undirected_edge(int first, int second, Cap cap) {
static_assert(std::numeric_limits<Cap>::is_signed);
assert(0 <= first && first < _n);
assert(0 <= second && second < _n);
assert(Cap(0) <= cap);
assert(cap <= std::numeric_limits<Cap>::max() / Cap(2));
int id = int(_pos.size());
int first_id = int(_g[first].size());
int second_id = int(_g[second].size());
if (first == second) second_id++;
_pos.push_back(Position{first, ~first_id});
_g[first].push_back(InternalEdge{second, second_id, cap});
_g[second].push_back(InternalEdge{first, first_id, cap});
return id;
}
Edge get_edge(int i) const {
assert(0 <= i && i < int(_pos.size()));
const auto& position = _pos[i];
int from = position.from;
bool undirected = position.edge < 0;
int idx = undirected ? ~position.edge : position.edge;
const auto& e = _g[from][idx];
const auto& re = _g[e.to][e.rev];
if (undirected) {
return Edge{
from,
e.to,
(e.cap + re.cap) / Cap(2),
(re.cap - e.cap) / Cap(2)
};
}
return Edge{from, e.to, e.cap + re.cap, re.cap};
}
std::vector<Edge> edges() const {
std::vector<Edge> result;
result.reserve(_pos.size());
for (int i = 0; i < int(_pos.size()); i++) result.push_back(get_edge(i));
return result;
}
void change_edge(int i, Cap new_cap, Cap new_flow) {
assert(0 <= i && i < int(_pos.size()));
assert(Cap(0) <= new_cap);
auto& position = _pos[i];
int from = position.from;
bool undirected = position.edge < 0;
int idx = undirected ? ~position.edge : position.edge;
auto& e = _g[from][idx];
auto& re = _g[e.to][e.rev];
if (undirected) {
assert(new_cap <= std::numeric_limits<Cap>::max() / Cap(2));
assert(-new_cap <= new_flow && new_flow <= new_cap);
e.cap = new_cap - new_flow;
re.cap = new_cap + new_flow;
} else {
assert(Cap(0) <= new_flow && new_flow <= new_cap);
e.cap = new_cap - new_flow;
re.cap = new_flow;
}
}
Cap max_flow(int s, int t) {
assert(0 <= s && s < _n);
assert(0 <= t && t < _n);
assert(s != t);
return highest_label_preflow_push(s, t);
}
Cap max_flow_push_relabel(int s, int t) {
assert(0 <= s && s < _n);
assert(0 <= t && t < _n);
assert(s != t);
return highest_label_preflow_push(s, t);
}
Cap max_flow_dinic(int s, int t) {
return max_flow(s, t, std::numeric_limits<Cap>::max());
}
Cap max_flow(int s, int t, Cap flow_limit) {
assert(0 <= s && s < _n);
assert(0 <= t && t < _n);
assert(s != t);
std::vector<int> work(3 * std::size_t(_n));
int* level = work.data();
int* iter = level + _n;
int* queue = iter + _n;
auto bfs = [&]() -> bool {
std::fill(level, level + _n, -1);
int head = 0;
int tail = 0;
level[s] = 0;
queue[tail++] = s;
while (head != tail) {
int v = queue[head++];
for (const auto& e : _g[v]) {
if (level[e.to] != -1 || e.cap == Cap(0)) continue;
level[e.to] = level[v] + 1;
if (e.to == t) return true;
queue[tail++] = e.to;
}
}
return level[t] != -1;
};
auto dfs = [&](auto&& self, int v, Cap up) -> Cap {
if (v == s) return up;
Cap result = Cap(0);
const int current_level = level[v];
auto& edges = _g[v];
const int edge_count = int(edges.size());
for (int& i = iter[v]; i < edge_count; i++) {
auto& e = edges[i];
if (level[e.to] + 1 != current_level) continue;
auto& reverse = _g[e.to][e.rev];
if (reverse.cap == Cap(0)) continue;
Cap d = self(
self,
e.to,
std::min(up - result, reverse.cap)
);
if (d == Cap(0)) continue;
e.cap += d;
reverse.cap -= d;
result += d;
if (result == up) return result;
}
level[v] = _n;
return result;
};
Cap flow = 0;
while (flow < flow_limit && bfs()) {
std::fill(iter, iter + _n, 0);
flow += dfs(dfs, t, flow_limit - flow);
}
return flow;
}
std::vector<bool> min_cut(int s) const {
assert(0 <= s && s < _n);
std::vector<bool> visited(_n, false);
std::vector<int> queue(_n);
int head = 0;
int tail = 0;
visited[s] = true;
queue[tail++] = s;
while (head != tail) {
int v = queue[head++];
for (const auto& e : _g[v]) {
if (e.cap == Cap(0) || visited[e.to]) continue;
visited[e.to] = true;
queue[tail++] = e.to;
}
}
return visited;
}
};
} // namespace flow
} // namespace m1une
#line 18 "graph/flow/min_cost_flow.hpp"
namespace m1une {
namespace flow {
template <class Cap, class Cost>
struct MinCostFlow {
struct Edge {
int from;
int to;
Cap cap;
Cap flow;
Cost cost;
};
private:
struct InternalEdge {
int to;
int rev;
Cap cap;
Cost cost;
};
int _n;
std::vector<std::pair<int, int>> _pos;
std::vector<std::vector<InternalEdge>> _g;
bool _has_negative_cost;
bool _has_flow;
template <class Key>
struct RadixHeap {
using Unsigned = std::make_unsigned_t<Key>;
static constexpr int bits = std::numeric_limits<Unsigned>::digits;
std::array<std::vector<std::pair<Unsigned, int>>, bits + 1> bucket;
Unsigned last = 0;
std::size_t count = 0;
static int index(Unsigned first, Unsigned second) {
return int(std::bit_width(first ^ second));
}
void clear() {
for (auto& values : bucket) values.clear();
last = 0;
count = 0;
}
bool empty() const {
return count == 0;
}
void push(Key key, int vertex) {
Unsigned value = static_cast<Unsigned>(key);
assert(last <= value);
bucket[index(value, last)].emplace_back(value, vertex);
count++;
}
std::pair<Key, int> pop() {
if (bucket[0].empty()) {
int i = 1;
while (bucket[i].empty()) i++;
last = bucket[i][0].first;
for (const auto& value : bucket[i]) {
last = std::min(last, value.first);
}
for (const auto& value : bucket[i]) {
bucket[index(value.first, last)].push_back(value);
}
bucket[i].clear();
}
auto [key, vertex] = bucket[0].back();
bucket[0].pop_back();
count--;
return {static_cast<Key>(key), vertex};
}
};
template <class Key>
struct BinaryHeap {
using Value = std::pair<Key, int>;
std::vector<Value> heap;
void clear() {
heap.clear();
}
bool empty() const {
return heap.empty();
}
void push(Key key, int vertex) {
heap.emplace_back(key, vertex);
std::push_heap(heap.begin(), heap.end(), std::greater<Value>());
}
Value pop() {
std::pop_heap(heap.begin(), heap.end(), std::greater<Value>());
Value result = heap.back();
heap.pop_back();
return result;
}
};
template <
class Key,
bool UseRadix =
std::numeric_limits<Key>::is_integer && sizeof(Key) <= 8
>
struct HeapSelector {
using Type = BinaryHeap<Key>;
};
template <class Key>
struct HeapSelector<Key, true> {
using Type = RadixHeap<Key>;
};
bool use_network_simplex(int s, int t, Cap flow_limit) const {
if (_has_negative_cost) return false;
if (_pos.size() < 64) return false;
auto add_saturated = [](Cap first, Cap second) {
const Cap maximum = std::numeric_limits<Cap>::max();
return maximum - first < second ? maximum : first + second;
};
struct TerminalCapacity {
Cap total = Cap(0);
std::array<Cap, 7> largest{};
};
auto add_capacity = [&](TerminalCapacity& terminal, Cap cap) {
terminal.total = add_saturated(terminal.total, cap);
for (Cap& current : terminal.largest) {
if (cap <= current) break;
std::swap(cap, current);
}
};
TerminalCapacity source;
for (const auto& e : _g[s]) {
if (e.to == s) continue;
add_capacity(source, e.cap);
}
TerminalCapacity sink;
for (const auto& e : _g[t]) {
if (e.to == t) continue;
Cap cap = _g[e.to][e.rev].cap;
add_capacity(sink, cap);
}
Cap target = std::min(
flow_limit,
std::min(source.total, sink.total)
);
if (target == Cap(0)) return false;
auto requires_eight_arcs = [&](const TerminalCapacity& terminal) {
Cap sum = Cap(0);
for (Cap cap : terminal.largest) {
sum = add_saturated(sum, cap);
}
return sum < target;
};
return requires_eight_arcs(source) && requires_eight_arcs(sink);
}
std::pair<Cap, Cost> network_simplex_flow(
int s,
int t,
Cap flow_limit
) {
struct ResidualArc {
int edge;
bool reverse;
};
using Solver = BoundedMinCostFlow<Cap, Cost, Cost>;
std::vector<ResidualArc> arcs;
arcs.reserve(2 * _pos.size());
for (int i = 0; i < int(_pos.size()); i++) {
auto [from, idx] = _pos[i];
const auto& e = _g[from][idx];
const auto& reverse = _g[e.to][e.rev];
if (e.cap != Cap(0)) {
arcs.push_back(ResidualArc{i, false});
}
if (reverse.cap != Cap(0)) {
arcs.push_back(ResidualArc{i, true});
}
}
auto add_saturated = [](Cap first, Cap second, bool& exact) {
const Cap maximum = std::numeric_limits<Cap>::max();
if (maximum - first < second) {
exact = false;
return maximum;
}
return first + second;
};
bool source_capacity_exact = true;
Cap source_capacity = Cap(0);
for (const auto& e : _g[s]) {
if (e.to == s) continue;
source_capacity = add_saturated(
source_capacity,
e.cap,
source_capacity_exact
);
}
bool sink_capacity_exact = true;
Cap sink_capacity = Cap(0);
for (const auto& e : _g[t]) {
if (e.to == t) continue;
sink_capacity = add_saturated(
sink_capacity,
_g[e.to][e.rev].cap,
sink_capacity_exact
);
}
Cap target = std::min(
flow_limit,
std::min(source_capacity, sink_capacity)
);
if (target == Cap(0)) return {Cap(0), Cost(0)};
struct ArcData {
int from;
int to;
Cap cap;
Cost cost;
};
auto arc_data = [&](const ResidualArc& arc) {
auto [from, idx] = _pos[arc.edge];
const auto& e = _g[from][idx];
const auto& reverse = _g[e.to][e.rev];
return arc.reverse
? ArcData{e.to, from, reverse.cap, reverse.cost}
: ArcData{from, e.to, e.cap, e.cost};
};
auto apply_flow = [&](const ResidualArc& arc, Cap amount) {
auto [from, idx] = _pos[arc.edge];
auto& e = _g[from][idx];
auto& reverse = _g[e.to][e.rev];
if (arc.reverse) {
reverse.cap -= amount;
e.cap += amount;
} else {
e.cap -= amount;
reverse.cap += amount;
}
};
bool target_infeasible = false;
if (
source_capacity_exact && sink_capacity_exact &&
target == source_capacity && target == sink_capacity
) {
Solver terminal_solver(_n);
terminal_solver.reserve_edges(int(arcs.size()));
std::vector<Cap> balance(_n, Cap(0));
std::vector<int> internal_arcs;
std::vector<int> fixed_arcs;
internal_arcs.reserve(arcs.size());
fixed_arcs.reserve(_g[s].size() + _g[t].size());
Cost fixed_cost = Cost(0);
for (int i = 0; i < int(arcs.size()); i++) {
ArcData data = arc_data(arcs[i]);
if (data.from == s) {
if (data.to == s) continue;
fixed_arcs.push_back(i);
fixed_cost += Cost(data.cap) * data.cost;
if (data.to != t) balance[data.to] += data.cap;
} else if (data.to == t) {
if (data.from == t) continue;
fixed_arcs.push_back(i);
fixed_cost += Cost(data.cap) * data.cost;
balance[data.from] -= data.cap;
} else if (data.to != s && data.from != t) {
terminal_solver.add_edge(
data.from,
data.to,
Cap(0),
data.cap,
data.cost
);
internal_arcs.push_back(i);
}
}
auto terminal_result = terminal_solver.min_cost_flow(balance);
if (terminal_result) {
for (int i : fixed_arcs) {
apply_flow(arcs[i], arc_data(arcs[i]).cap);
}
for (int i = 0; i < int(internal_arcs.size()); i++) {
apply_flow(
arcs[internal_arcs[i]],
terminal_result->flow(i)
);
}
_has_flow = true;
return {target, fixed_cost + terminal_result->cost};
}
target_infeasible = true;
}
Solver solver(_n);
solver.reserve_edges(int(arcs.size()));
for (const auto& arc : arcs) {
ArcData data = arc_data(arc);
solver.add_edge(
data.from,
data.to,
Cap(0),
data.cap,
data.cost
);
}
Cap sent = target;
std::optional<typename Solver::Result> result;
if (
!target_infeasible &&
target != std::numeric_limits<Cap>::max()
) {
result = solver.min_cost_st_flow(s, t, target);
}
if (!result) {
MaxFlow<Cap> feasible(_n);
feasible.reserve_edges(int(arcs.size()));
for (const auto& arc : arcs) {
auto [from, idx] = _pos[arc.edge];
const auto& e = _g[from][idx];
const auto& reverse = _g[e.to][e.rev];
if (arc.reverse) {
feasible.add_edge(e.to, from, reverse.cap);
} else {
feasible.add_edge(from, e.to, e.cap);
}
}
sent = feasible.max_flow(s, t, target);
if (sent == Cap(0)) return {Cap(0), Cost(0)};
result = solver.min_cost_st_flow(s, t, sent);
}
assert(result.has_value());
for (int i = 0; i < int(arcs.size()); i++) {
auto [from, idx] = _pos[arcs[i].edge];
auto& e = _g[from][idx];
auto& reverse = _g[e.to][e.rev];
Cap amount = result->flow(i);
if (arcs[i].reverse) {
reverse.cap -= amount;
e.cap += amount;
} else {
e.cap -= amount;
reverse.cap += amount;
}
}
_has_flow = true;
return {sent, result->cost};
}
void init_potential(int s, std::vector<Cost>& potential, Cost cost_inf) const {
if (!_has_negative_cost && !_has_flow) {
potential.assign(_n, Cost(0));
return;
}
potential.assign(_n, cost_inf);
potential[s] = Cost(0);
for (int iter = 0; iter < _n - 1; iter++) {
bool updated = false;
for (int v = 0; v < _n; v++) {
if (potential[v] == cost_inf) continue;
for (const auto& e : _g[v]) {
if (e.cap == Cap(0)) continue;
Cost nd = potential[v] + e.cost;
if (nd < potential[e.to]) {
potential[e.to] = nd;
updated = true;
}
}
}
if (!updated) break;
}
for (int v = 0; v < _n; v++) {
if (potential[v] == cost_inf) potential[v] = Cost(0);
}
}
public:
MinCostFlow() : MinCostFlow(0) {}
explicit MinCostFlow(int n)
: _n(n), _g(n), _has_negative_cost(false), _has_flow(false) {
assert(0 <= n);
}
int size() const {
return _n;
}
int edge_count() const {
return int(_pos.size());
}
void reserve_edges(int edge_count) {
assert(0 <= edge_count);
_pos.reserve(edge_count);
if (_n == 0 || edge_count == 0 ||
2 * std::size_t(edge_count) < std::size_t(_n)) {
return;
}
const std::size_t average_degree =
(3 * std::size_t(edge_count) + std::size_t(_n) - 1)
/ std::size_t(_n);
for (auto& edges : _g) edges.reserve(average_degree);
}
void reserve_edges(int edge_count, const std::vector<int>& degrees) {
assert(0 <= edge_count);
assert(int(degrees.size()) == _n);
_pos.reserve(edge_count);
for (int v = 0; v < _n; v++) {
assert(0 <= degrees[v]);
_g[v].reserve(degrees[v]);
}
}
int add_edge(int from, int to, Cap cap, Cost cost) {
assert(0 <= from && from < _n);
assert(0 <= to && to < _n);
assert(Cap(0) <= cap);
_has_negative_cost = _has_negative_cost || cost < Cost(0);
int id = int(_pos.size());
int from_id = int(_g[from].size());
int to_id = int(_g[to].size());
if (from == to) to_id++;
_pos.emplace_back(from, from_id);
_g[from].push_back(InternalEdge{to, to_id, cap, cost});
_g[to].push_back(InternalEdge{from, from_id, Cap(0), -cost});
return id;
}
Edge get_edge(int i) const {
assert(0 <= i && i < int(_pos.size()));
auto [from, idx] = _pos[i];
const auto& e = _g[from][idx];
const auto& re = _g[e.to][e.rev];
return Edge{from, e.to, e.cap + re.cap, re.cap, e.cost};
}
std::vector<Edge> edges() const {
std::vector<Edge> result;
result.reserve(_pos.size());
for (int i = 0; i < int(_pos.size()); i++) result.push_back(get_edge(i));
return result;
}
std::pair<Cap, Cost> flow(int s, int t) {
return flow(s, t, std::numeric_limits<Cap>::max());
}
std::pair<Cap, Cost> flow(int s, int t, Cap flow_limit) {
assert(0 <= s && s < _n);
assert(0 <= t && t < _n);
assert(s != t);
assert(Cap(0) <= flow_limit);
if (flow_limit == Cap(0)) return {Cap(0), Cost(0)};
if constexpr (
std::numeric_limits<Cap>::is_integer &&
std::numeric_limits<Cap>::is_signed &&
std::numeric_limits<Cost>::is_signed
) {
if (use_network_simplex(s, t, flow_limit)) {
return network_simplex_flow(s, t, flow_limit);
}
}
auto result = slope(s, t, flow_limit);
return result.back();
}
std::vector<std::pair<Cap, Cost>> slope(int s, int t) {
return slope(s, t, std::numeric_limits<Cap>::max());
}
std::vector<std::pair<Cap, Cost>> slope(int s, int t, Cap flow_limit) {
assert(0 <= s && s < _n);
assert(0 <= t && t < _n);
assert(s != t);
assert(Cap(0) <= flow_limit);
const Cost cost_inf = std::numeric_limits<Cost>::max() / Cost(4);
std::vector<Cost> potential, dist(_n);
std::vector<int> prev_v(_n), prev_e(_n);
std::vector<int> settled;
settled.reserve(_n);
typename HeapSelector<Cost>::Type que;
init_potential(s, potential, cost_inf);
std::vector<std::pair<Cap, Cost>> result;
result.emplace_back(Cap(0), Cost(0));
Cap flow = 0;
Cost cost = 0;
while (flow < flow_limit) {
std::fill(dist.begin(), dist.end(), cost_inf);
dist[s] = Cost(0);
settled.clear();
que.clear();
que.push(Cost(0), s);
while (!que.empty()) {
auto [d, v] = que.pop();
if (dist[v] != d) continue;
settled.push_back(v);
if (v == t) break;
for (int i = 0; i < int(_g[v].size()); i++) {
const auto& e = _g[v][i];
if (e.cap == Cap(0)) continue;
Cost nd = d + e.cost + potential[v] - potential[e.to];
if (nd >= dist[e.to]) continue;
dist[e.to] = nd;
prev_v[e.to] = v;
prev_e[e.to] = i;
que.push(nd, e.to);
}
}
if (dist[t] == cost_inf) break;
for (int v : settled) {
potential[v] += dist[v] - dist[t];
}
Cap add = flow_limit - flow;
for (int v = t; v != s; v = prev_v[v]) {
add = std::min(add, _g[prev_v[v]][prev_e[v]].cap);
}
Cost path_cost = potential[t] - potential[s];
for (int v = t; v != s; v = prev_v[v]) {
auto& e = _g[prev_v[v]][prev_e[v]];
e.cap -= add;
_g[e.to][e.rev].cap += add;
}
flow += add;
cost += Cost(add) * path_cost;
result.emplace_back(flow, cost);
}
_has_flow = _has_flow || flow != Cap(0);
return result;
}
};
} // namespace flow
} // namespace m1une