Flow
(graph/flow/flow.hpp)
- View this file on GitHub
- Last update: 2026-08-04 16:49:58+09:00
- Include:
#include "graph/flow/flow.hpp"
Overview
graph/flow/flow.hpp includes flow-network algorithms. Flow networks are
directed: an edge u -> v only sends flow from u to v.
For an undirected capacity between u and v, MaxFlow provides
add_undirected_edge(u, v, cap), which stores one shared-capacity residual
pair. Other directed-flow classes represent it with two directed edges.
Included Headers
| Header | Graph orientation | Contents |
|---|---|---|
graph/flow/bounded_flow.hpp |
Directed flow network | Feasible flow with lower/upper bounds, balances, and negative flow intervals. |
graph/flow/bounded_min_cost_flow.hpp |
Directed flow network | Minimum-cost feasible flow with lower/upper bounds, balances, and negative flow intervals. |
graph/flow/gomory_hu.hpp |
Undirected capacitated graph | Gomory-Hu cut tree and pairwise minimum-cut queries. |
graph/flow/max_flow.hpp |
Directed or shared-capacity undirected network | $O(N^2 \sqrt M)$ highest-label preflow-push, flow-limited Dinic, and minimum cut. |
graph/flow/min_cost_flow.hpp |
Directed flow network | Minimum-cost flow with potentials. |
Complexity
This header is an include bundle and provides no runtime operation by itself. See the included algorithm pages for public interfaces and complexities.
Depends on
- Bounded Flow
(graph/flow/bounded_flow.hpp)
- Bounded Min Cost Flow
(graph/flow/bounded_min_cost_flow.hpp)
- Gomory-Hu Tree
(graph/flow/gomory_hu.hpp)
- Max Flow
(graph/flow/max_flow.hpp)
- Min Cost Flow
(graph/flow/min_cost_flow.hpp)
Required by
Verified with
- verify/graph/cow_game.test.cpp
- verify/graph/flow/flow_algorithms.test.cpp
- verify/graph/graph_algorithms.test.cpp
- verify/graph/range_edge_graph.test.cpp
Code
#ifndef M1UNE_FLOW_FLOW_HPP
#define M1UNE_FLOW_FLOW_HPP 1
#include "bounded_flow.hpp"
#include "bounded_min_cost_flow.hpp"
#include "gomory_hu.hpp"
#include "max_flow.hpp"
#include "min_cost_flow.hpp"
#endif // M1UNE_FLOW_FLOW_HPP#line 1 "graph/flow/flow.hpp"
#line 1 "graph/flow/bounded_flow.hpp"
#include <cassert>
#include <optional>
#include <vector>
#line 1 "graph/flow/max_flow.hpp"
#include <algorithm>
#line 6 "graph/flow/max_flow.hpp"
#include <cstddef>
#include <limits>
#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 9 "graph/flow/bounded_flow.hpp"
namespace m1une {
namespace flow {
template <class Cap>
struct BoundedFlow {
struct Edge {
int from;
int to;
Cap lower;
Cap upper;
};
struct ResultEdge {
int from;
int to;
Cap lower;
Cap upper;
Cap flow;
};
struct Result {
std::vector<ResultEdge> edges;
std::vector<Cap> balance;
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:
int _n;
std::vector<Edge> _edges;
std::vector<Cap> _balance;
public:
BoundedFlow() : BoundedFlow(0) {}
explicit BoundedFlow(int n) : _n(n), _balance(n, Cap(0)) {
assert(0 <= n);
}
int size() const {
return _n;
}
int edge_count() const {
return int(_edges.size());
}
int add_edge(int from, int to, Cap lower, Cap upper) {
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});
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> feasible_flow() const {
return feasible_flow(_balance);
}
std::optional<Result> feasible_flow(const std::vector<Cap>& balance) const {
assert(int(balance.size()) == _n);
int ss = _n, tt = _n + 1;
MaxFlow<Cap> mf(_n + 2);
std::vector<int> edge_ids;
edge_ids.reserve(_edges.size());
std::vector<Cap> need = balance;
for (const auto& e : _edges) {
edge_ids.push_back(mf.add_edge(e.from, e.to, e.upper - e.lower));
need[e.from] -= e.lower;
need[e.to] += e.lower;
}
Cap positive_sum = Cap(0), negative_sum = Cap(0);
for (int v = 0; v < _n; v++) {
if (need[v] > Cap(0)) {
positive_sum += need[v];
mf.add_edge(ss, v, need[v]);
} else if (need[v] < Cap(0)) {
negative_sum += -need[v];
mf.add_edge(v, tt, -need[v]);
}
}
if (positive_sum != negative_sum) return std::nullopt;
if (mf.max_flow(ss, tt) != positive_sum) return std::nullopt;
Result result;
result.balance = balance;
result.edges.reserve(_edges.size());
for (int i = 0; i < int(_edges.size()); i++) {
auto used = mf.get_edge(edge_ids[i]).flow;
const auto& e = _edges[i];
result.edges.push_back(ResultEdge{e.from, e.to, e.lower, e.upper, e.lower + used});
}
return result;
}
std::optional<Result> feasible_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 feasible_flow(balance);
}
};
template <class Cap>
using BFlow = BoundedFlow<Cap>;
} // namespace flow
} // namespace m1une
#line 1 "graph/flow/bounded_min_cost_flow.hpp"
#line 7 "graph/flow/bounded_min_cost_flow.hpp"
#include <cmath>
#include <functional>
#line 11 "graph/flow/bounded_min_cost_flow.hpp"
#include <queue>
#include <utility>
#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/gomory_hu.hpp"
#line 9 "graph/flow/gomory_hu.hpp"
namespace m1une {
namespace flow {
template <class Cap>
struct GomoryHu {
struct Edge {
int u;
int v;
Cap cap;
};
private:
struct FlowEdge {
int to;
int rev;
Cap cap;
Cap initial_cap;
};
int _n;
bool _built = false;
std::vector<Edge> _edges;
std::vector<Edge> _tree_edges;
std::vector<int> _parent;
std::vector<Cap> _cut_value;
std::vector<std::vector<std::pair<int, Cap>>> _tree;
std::vector<std::vector<int>> _up;
std::vector<std::vector<Cap>> _minimum;
std::vector<int> _depth;
std::vector<std::vector<FlowEdge>> _graph;
std::vector<Cap> _excess;
std::vector<int> _height;
std::vector<int> _height_count;
std::vector<int> _current;
std::vector<bool> _active;
std::vector<std::vector<int>> _buckets;
std::vector<int> _queue;
int _highest;
long long _work;
long long _work_limit;
void add_flow_edge(int u, int v, Cap cap) {
if (u == v || cap == Cap(0)) return;
int ui = int(_graph[u].size());
int vi = int(_graph[v].size());
_graph[u].push_back(FlowEdge{v, vi, cap, cap});
_graph[v].push_back(FlowEdge{u, ui, cap, cap});
}
void reset_flow() {
for (auto& edges : _graph) {
for (auto& edge : edges) edge.cap = edge.initial_cap;
}
}
void activate(int v, int s, int t) {
int dead = 2 * _n;
if (v == s || v == t || _active[v] || _excess[v] == Cap(0) || _height[v] >= dead) return;
_active[v] = true;
_buckets[_height[v]].push_back(v);
_highest = std::max(_highest, _height[v]);
}
void rebuild_buckets(int s, int t) {
for (auto& bucket : _buckets) bucket.clear();
std::fill(_active.begin(), _active.end(), false);
_highest = -1;
for (int v = 0; v < _n; v++) activate(v, s, t);
}
void global_relabel(int s, int t) {
int dead = 2 * _n;
int unreachable = _n + 1;
std::fill(_height.begin(), _height.end(), unreachable);
std::fill(_height_count.begin(), _height_count.end(), 0);
std::fill(_current.begin(), _current.end(), 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& edge : _graph[v]) {
const FlowEdge& reverse = _graph[edge.to][edge.rev];
if (reverse.cap == Cap(0) || _height[edge.to] != unreachable) continue;
_height[edge.to] = _height[v] + 1;
_queue[tail++] = edge.to;
}
}
for (int v = 0; v < _n; v++) {
_height[v] = std::min(_height[v], dead);
_height_count[_height[v]]++;
}
rebuild_buckets(s, t);
_work = 0;
}
void push(int v, FlowEdge& edge, int s, int t) {
if (edge.cap == Cap(0) || _height[v] != _height[edge.to] + 1) return;
Cap sent = std::min(_excess[v], edge.cap);
if (sent == Cap(0)) return;
bool was_zero = _excess[edge.to] == Cap(0);
edge.cap -= sent;
_graph[edge.to][edge.rev].cap += sent;
_excess[v] -= sent;
_excess[edge.to] += sent;
if (was_zero) activate(edge.to, s, t);
}
void gap(int height, int s, int t) {
int unreachable = _n + 1;
for (int v = 0; v < _n; v++) {
if (v == s || v == t || _height[v] <= height || _height[v] >= _n) continue;
_height_count[_height[v]]--;
_height[v] = unreachable;
_height_count[_height[v]]++;
_current[v] = 0;
}
rebuild_buckets(s, t);
}
bool relabel(int v, int s, int t) {
int dead = 2 * _n;
int old_height = _height[v];
int new_height = dead;
_work += int(_graph[v].size());
for (const auto& edge : _graph[v]) {
if (edge.cap != Cap(0)) new_height = std::min(new_height, _height[edge.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, s, t);
return true;
}
return false;
}
void discharge(int v, int s, int t) {
while (_excess[v] != Cap(0) && _height[v] < 2 * _n) {
if (_current[v] == int(_graph[v].size())) {
if (relabel(v, s, t)) return;
continue;
}
FlowEdge& edge = _graph[v][_current[v]];
_work++;
if (edge.cap != Cap(0) && _height[v] == _height[edge.to] + 1) {
push(v, edge, s, t);
} else {
_current[v]++;
}
}
activate(v, s, t);
}
Cap max_flow(int s, int t) {
reset_flow();
std::fill(_excess.begin(), _excess.end(), Cap(0));
for (auto& edge : _graph[s]) {
Cap sent = edge.cap;
if (sent == Cap(0)) continue;
edge.cap = Cap(0);
_graph[edge.to][edge.rev].cap += sent;
_excess[edge.to] += sent;
}
global_relabel(s, t);
while (_highest >= 0) {
if (_buckets[_highest].empty()) {
_highest--;
continue;
}
int v = _buckets[_highest].back();
_buckets[_highest].pop_back();
if (!_active[v] || _height[v] != _highest) continue;
_active[v] = false;
discharge(v, s, t);
if (_work >= _work_limit) global_relabel(s, t);
}
return _excess[t];
}
std::vector<bool> source_side(int s) {
std::vector<bool> visited(_n, false);
int head = 0;
int tail = 0;
visited[s] = true;
_queue[tail++] = s;
while (head < tail) {
int v = _queue[head++];
for (const auto& edge : _graph[v]) {
if (edge.cap == Cap(0) || visited[edge.to]) continue;
visited[edge.to] = true;
_queue[tail++] = edge.to;
}
}
return visited;
}
void build_query_table() {
int log = 1;
while ((1LL << log) <= std::max(1, _n)) log++;
const Cap infinity = std::numeric_limits<Cap>::max();
_up.assign(log, std::vector<int>(_n, 0));
_minimum.assign(log, std::vector<Cap>(_n, infinity));
_depth.assign(_n, 0);
if (_n == 0) return;
std::vector<int> order;
order.reserve(_n);
order.push_back(0);
for (int i = 0; i < int(order.size()); i++) {
int v = order[i];
for (auto [to, cap] : _tree[v]) {
if (to == _up[0][v] && v != 0) continue;
_up[0][to] = v;
_minimum[0][to] = cap;
_depth[to] = _depth[v] + 1;
order.push_back(to);
}
}
for (int k = 1; k < log; k++) {
for (int v = 0; v < _n; v++) {
int middle = _up[k - 1][v];
_up[k][v] = _up[k - 1][middle];
_minimum[k][v] = std::min(_minimum[k - 1][v], _minimum[k - 1][middle]);
}
}
}
public:
GomoryHu() : GomoryHu(0) {}
explicit GomoryHu(int n) : _n(n) {
assert(0 <= n);
}
int size() const {
return _n;
}
int edge_count() const {
return int(_edges.size());
}
int add_edge(int u, int v, Cap cap) {
assert(0 <= u && u < _n);
assert(0 <= v && v < _n);
assert(Cap(0) <= cap);
_built = false;
int id = int(_edges.size());
_edges.push_back(Edge{u, v, cap});
return id;
}
void build() {
std::vector<Edge> flow_edges;
flow_edges.reserve(_edges.size());
for (auto edge : _edges) {
if (edge.u == edge.v || edge.cap == Cap(0)) continue;
if (edge.u > edge.v) std::swap(edge.u, edge.v);
flow_edges.push_back(edge);
}
std::sort(flow_edges.begin(), flow_edges.end(), [](const Edge& lhs, const Edge& rhs) {
return std::pair<int, int>(lhs.u, lhs.v) < std::pair<int, int>(rhs.u, rhs.v);
});
int unique_edges = 0;
for (const auto& edge : flow_edges) {
if (unique_edges > 0 && flow_edges[unique_edges - 1].u == edge.u &&
flow_edges[unique_edges - 1].v == edge.v) {
flow_edges[unique_edges - 1].cap += edge.cap;
} else {
flow_edges[unique_edges++] = edge;
}
}
flow_edges.resize(unique_edges);
_graph.assign(_n, {});
std::vector<int> degree(_n, 0);
for (const auto& edge : flow_edges) {
degree[edge.u]++;
degree[edge.v]++;
}
for (int v = 0; v < _n; v++) _graph[v].reserve(degree[v]);
for (const auto& edge : flow_edges) add_flow_edge(edge.u, edge.v, edge.cap);
_excess.resize(_n);
_height.resize(_n);
_height_count.resize(2 * _n + 1);
_current.resize(_n);
_active.resize(_n);
_buckets.resize(2 * _n + 1);
_queue.resize(_n);
long long arc_count = 0;
for (const auto& edges : _graph) arc_count += int(edges.size());
_work_limit = std::max(1LL, 4 * arc_count + _n);
_parent.assign(_n, 0);
_cut_value.assign(_n, std::numeric_limits<Cap>::max());
for (int s = 1; s < _n; s++) {
int t = _parent[s];
Cap flow = max_flow(s, t);
std::vector<bool> cut = source_side(s);
for (int v = s + 1; v < _n; v++) {
if (_parent[v] == t && cut[v]) _parent[v] = s;
}
if (cut[_parent[t]]) {
_parent[s] = _parent[t];
_parent[t] = s;
_cut_value[s] = _cut_value[t];
_cut_value[t] = flow;
} else {
_cut_value[s] = flow;
}
}
_tree.assign(_n, {});
_tree_edges.clear();
if (_n > 0) _tree_edges.reserve(_n - 1);
for (int v = 1; v < _n; v++) {
int p = _parent[v];
Cap cap = _cut_value[v];
_tree_edges.push_back(Edge{v, p, cap});
_tree[v].emplace_back(p, cap);
_tree[p].emplace_back(v, cap);
}
build_query_table();
_built = true;
}
const std::vector<Edge>& tree_edges() const {
assert(_built);
return _tree_edges;
}
const std::vector<int>& parent() const {
assert(_built);
return _parent;
}
const std::vector<Cap>& cut_values() const {
assert(_built);
return _cut_value;
}
Cap min_cut(int u, int v) const {
assert(_built);
assert(0 <= u && u < _n);
assert(0 <= v && v < _n);
assert(u != v);
Cap result = std::numeric_limits<Cap>::max();
if (_depth[u] < _depth[v]) std::swap(u, v);
int difference = _depth[u] - _depth[v];
for (int k = 0; difference > 0; k++, difference >>= 1) {
if ((difference & 1) == 0) continue;
result = std::min(result, _minimum[k][u]);
u = _up[k][u];
}
if (u == v) return result;
for (int k = int(_up.size()) - 1; k >= 0; k--) {
if (_up[k][u] == _up[k][v]) continue;
result = std::min(result, _minimum[k][u]);
result = std::min(result, _minimum[k][v]);
u = _up[k][u];
v = _up[k][v];
}
result = std::min(result, _minimum[0][u]);
result = std::min(result, _minimum[0][v]);
return result;
}
};
} // namespace flow
} // namespace m1une
#line 1 "graph/flow/min_cost_flow.hpp"
#line 5 "graph/flow/min_cost_flow.hpp"
#include <array>
#include <bit>
#line 12 "graph/flow/min_cost_flow.hpp"
#include <type_traits>
#line 15 "graph/flow/min_cost_flow.hpp"
#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
#line 9 "graph/flow/flow.hpp"