min_cost_b_flow.hpp¶
用容量缩放处理点供需、边流量上下界与最小/最大费用；可独立指定总费用类型，并用 Dijkstra 生成规范化势和最优性证书。
Complexity: Time: \(O((V + E)^{2} \log (V + E + 1) \log (U + 1))\) conservative scaling bound for a fixed scaling factor; get_potential: \(O((V + E) \log (V + E + 1))\). \(U\) bounds absolute initial residual capacities/balances. Space: \(O(V + E)\).
Source: Library Checker
AC 记录：min_cost_b_flow。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
// SPDX-License-Identifier: Apache-2.0
// Adapted from Library Checker, with changes for libra.
// Sources and license: third_party/README.md.
/// @complexity Time: O((V+E)^2 log(V+E+1) log(U+1)) conservative scaling bound
/// for a fixed scaling factor; get_potential: O((V+E) log(V+E+1)).
/// U bounds absolute initial residual capacities/balances. Space: O(V+E).
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <functional>
#include <limits>
#include <numeric>
#include <queue>
#include <tuple>
#include <utility>
#include <vector>
namespace noya {
/// @brief Minimum- or maximum-cost feasible b-flow with lower and upper bounds.
/// Lower bounds are represented as reverse residual capacity and vertex
/// supplies are maintained as excesses. Capacity scaling repeatedly saturates
/// negative reduced-cost residual edges, then a multi-source Dijkstra dual
/// step updates potentials and a primal step routes multiples of the scale from
/// excess vertices to deficits. When the scale reaches one, residual reduced
/// costs certify optimality; a remaining excess certifies infeasibility.
enum Objective {
MINIMIZE = 1,
MAXIMIZE = -1,
};
enum class Status {
OPTIMAL,
INFEASIBLE,
};
template <class Flow, class Cost, Objective objective = Objective::MINIMIZE,
Flow SCALING_FACTOR = 2>
class MinCostFlow {
static_assert(SCALING_FACTOR >= 2);
static_assert(std::numeric_limits<Flow>::is_integer &&
std::numeric_limits<Flow>::is_signed);
using V_id = uint32_t;
using E_id = uint32_t;
class Edge {
friend class MinCostFlow;
V_id src, dst;
Flow flw, cap;
Cost cst;
E_id rev;
public:
Edge() = default;
Edge(const V_id src, const V_id dst, const Flow cap, const Cost cst,
const E_id rev)
: src(src), dst(dst), flw(0), cap(cap), cst(cst), rev(rev) {}
[[nodiscard]] Flow residual_cap() const { return cap - flw; }
};
public:
class EdgePtr {
friend class MinCostFlow;
const MinCostFlow *obj;
V_id v;
E_id e;
EdgePtr(const MinCostFlow *const obj, const V_id v, const E_id e)
: obj(obj), v(v), e(e) {}
[[nodiscard]] const Edge &edge() const { return obj->g[v][e]; }
[[nodiscard]] const Edge &rev() const {
const Edge &e = edge();
return obj->g[e.dst][e.rev];
}
public:
EdgePtr() = default;
[[nodiscard]] V_id src() const { return v; }
[[nodiscard]] V_id dst() const { return edge().dst; }
[[nodiscard]] Flow flow() const { return edge().flw; }
[[nodiscard]] Flow lower() const { return -rev().cap; }
[[nodiscard]] Flow upper() const { return edge().cap; }
[[nodiscard]] Cost cost() const { return edge().cst; }
[[nodiscard]] Cost gain() const { return -edge().cst; }
};
private:
V_id n;
std::vector<std::vector<Edge>> g;
std::vector<Flow> b;
public:
MinCostFlow() : n(0) {}
V_id add_vertex() {
++n;
g.resize(n);
b.resize(n);
return n - 1;
}
std::vector<V_id> add_vertices(const size_t siz) {
std::vector<V_id> ret(siz);
std::iota(std::begin(ret), std::end(ret), n);
n += siz;
g.resize(n);
b.resize(n);
return ret;
}
EdgePtr add_edge(const V_id src, const V_id dst, const Flow lo, const Flow hi,
const Cost cst) {
assert(src < n && dst < n);
const E_id e = g[src].size(), re = src == dst ? e + 1 : g[dst].size();
assert(lo <= hi);
g[src].emplace_back(Edge{src, dst, hi, cst * objective, re});
g[dst].emplace_back(Edge{dst, src, -lo, -cst * objective, e});
return EdgePtr{this, src, e};
}
void add_supply(const V_id v, const Flow amt) { b[v] += amt; }
void add_demand(const V_id v, const Flow amt) { b[v] -= amt; }
private:
// Variables used in calculation
const Cost inf = std::numeric_limits<Cost>::max();
Cost far;
std::vector<Cost> pot;
std::vector<Cost> dis;
std::vector<Edge *> fa; // out-forrest.
using Item = std::pair<Cost, V_id>;
struct Heap : std::priority_queue<Item, std::vector<Item>, std::greater<>> {
void clear() { this->c.clear(); }
};
Heap pq; // Empty between searches; retain the allocation for the next one.
std::vector<V_id> ev, dv;
void enqueue(Cost distance, V_id vertex) { pq.emplace(distance, vertex); }
void push(Edge &e, const Flow amt) {
e.flw += amt;
g[e.dst][e.rev].flw -= amt;
}
Cost residual_cost(const V_id src, const V_id dst, const Edge &e) {
return e.cst + pot[src] - pot[dst];
}
bool dual(const Flow dlt) {
ev.erase(std::remove_if(std::begin(ev), std::end(ev),
[&](const V_id v) { return b[v] < dlt; }),
std::end(ev));
dv.erase(std::remove_if(std::begin(dv), std::end(dv),
[&](const V_id v) { return b[v] > -dlt; }),
std::end(dv));
if (ev.empty() || dv.empty()) return false;
dis.assign(n, inf);
fa.assign(n, nullptr);
for (const auto v : ev)
enqueue(dis[v] = 0, v);
far = 0;
std::size_t nd = 0;
while (!pq.empty()) {
auto [d, u] = pq.top();
pq.pop();
if (dis[u] < d)
continue;
far = d;
if (b[u] <= -dlt)
++nd;
if (nd >= dv.size())
break;
for (auto &e : g[u]) {
if (e.residual_cap() < dlt)
continue;
const auto v = e.dst;
const auto vtx = d + residual_cost(u, v, e);
if (vtx >= dis[v])
continue;
enqueue(dis[v] = vtx, v);
fa[v] = &e;
}
}
pq.clear();
for (V_id v = 0; v < n; ++v) {
pot[v] += std::min(dis[v], far);
}
return nd > 0;
}
void primal(const Flow dlt) {
for (const auto t : dv) {
if (dis[t] > far)
continue;
Flow f = -b[t];
V_id v;
for (v = t; fa[v] != nullptr; v = fa[v]->src) {
f = std::min(f, fa[v]->residual_cap());
}
f = std::min(f, b[v]);
f -= f % dlt;
if (f <= 0)
continue;
for (v = t; fa[v] != nullptr;) {
auto &e = *fa[v];
push(e, f);
V_id u = fa[v]->src;
if (e.residual_cap() <= 0)
fa[v] = nullptr;
v = u;
}
b[t] += f;
b[v] -= f;
}
}
void saturate_negative(const Flow dlt) {
ev.clear();
dv.clear();
for (auto &es : g)
for (auto &e : es) {
Flow ca1 = e.residual_cap();
ca1 -= ca1 % dlt;
if (ca1 == 0) continue;
const Cost cs1 = residual_cost(e.src, e.dst, e);
if (cs1 < 0 || ca1 < 0) {
push(e, ca1);
if (e.src != e.dst) {
b[e.src] -= ca1;
b[e.dst] += ca1;
}
}
}
for (V_id v = 0; v < n; ++v)
if (b[v] != 0) {
(b[v] > 0 ? ev : dv).emplace_back(v);
}
}
public:
/// @brief Solve using signed integral flows. Flow must hold residuals,
/// excesses and their negations; Cost must hold reduced costs/distances.
/// Choose Total wide enough for all objective products and partial sums,
/// e.g. solve<__int128>() when edge costs fit int64 but their total does not.
template <class Total = Cost> std::pair<Status, Total> solve() {
pot.resize(n);
Flow in1 = 1;
for (const auto t : b)
in1 = std::max({in1, t, -t});
for (const auto &es : g)
for (const auto &e : es)
in1 = std::max({in1, e.residual_cap(), -e.residual_cap()});
Flow dlt = 1;
while (dlt <= in1 / SCALING_FACTOR)
dlt *= SCALING_FACTOR;
for (; dlt; dlt /= SCALING_FACTOR) {
saturate_negative(dlt);
while (dual(dlt))
primal(dlt);
}
Total val = get_result_value<Total>();
if (ev.empty() && dv.empty()) {
return {Status::OPTIMAL, val / Total(objective)};
} else {
return {Status::INFEASIBLE, val / Total(objective)};
}
}
std::vector<Cost> get_potential() {
assert(pot.size() == n); // Call after solve().
if (n == 0) return pot;
auto [low, high] = std::minmax_element(pot.begin(), pot.end());
if (*low == *high) {
std::fill(pot.begin(), pot.end(), Cost(0));
return pot;
}
// Canonical shortest distances from a new zero-cost source to every
// vertex. Reweight by the feasible pot so residual edges are nonnegative.
Cost top = *high;
std::vector<Cost> distance(n);
for (V_id v = 0; v < n; ++v) enqueue(distance[v] = top - pot[v], v);
while (!pq.empty()) {
auto [d, u] = pq.top();
pq.pop();
if (distance[u] < d) continue;
for (const auto &e : g[u]) if (e.residual_cap() > 0) {
Cost cost = residual_cost(u, e.dst, e);
assert(cost >= 0);
if (d >= distance[e.dst] || cost >= distance[e.dst] - d) continue;
distance[e.dst] = d + cost;
enqueue(distance[e.dst], e.dst);
}
}
for (V_id v = 0; v < n; ++v) pot[v] -= top - distance[v];
return pot;
}
// Internal minimization objective: negated user objective for MAXIMIZE.
template <class T> T get_result_value() {
T val = 0;
for (const auto &es : g)
for (std::size_t i = 0; i < es.size(); ++i) {
const auto &e = es[i];
// The two residual directions have the same flow*cost product.
// Count one, avoiding an unnecessary doubled intermediate total.
if (e.src < e.dst || (e.src == e.dst && i < e.rev))
val += T(e.flw) * T(e.cst);
}
return val;
}
std::vector<size_t> get_cut() {
std::vector<size_t> res;
if (ev.empty())
return res;
for (size_t v = 0; v < n; ++v) {
if (dv.empty() || (dis[v] < inf))
res.emplace_back(v);
}
return res;
}
};
} // namespace noya
// NOYA_MIN_COST_B_FLOW_HPP
#ifndef NOYA_MIN_COST_B_FLOW_HPP
#define NOYA_MIN_COST_B_FLOW_HPP 1
// SPDX-License-Identifier: Apache-2.0
// Adapted from Library Checker, with changes for libra.
// Sources and license: third_party/README.md.
/// @complexity Time: O((V+E)^2 log(V+E+1) log(U+1)) conservative scaling bound
/// for a fixed scaling factor; get_potential: O((V+E) log(V+E+1)).
/// U bounds absolute initial residual capacities/balances. Space: O(V+E).
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <functional>
#include <limits>
#include <numeric>
#include <queue>
#include <tuple>
#include <utility>
#include <vector>
namespace noya {
/// @brief Minimum- or maximum-cost feasible b-flow with lower and upper bounds.
/// Lower bounds are represented as reverse residual capacity and vertex
/// supplies are maintained as excesses. Capacity scaling repeatedly saturates
/// negative reduced-cost residual edges, then a multi-source Dijkstra dual
/// step updates potentials and a primal step routes multiples of the scale from
/// excess vertices to deficits. When the scale reaches one, residual reduced
/// costs certify optimality; a remaining excess certifies infeasibility.
enum Objective {
MINIMIZE = 1,
MAXIMIZE = -1,
};
enum class Status {
OPTIMAL,
INFEASIBLE,
};
template <class Flow, class Cost, Objective objective = Objective::MINIMIZE,
Flow SCALING_FACTOR = 2>
class MinCostFlow {
static_assert(SCALING_FACTOR >= 2);
static_assert(std::numeric_limits<Flow>::is_integer &&
std::numeric_limits<Flow>::is_signed);
using V_id = uint32_t;
using E_id = uint32_t;
class Edge {
friend class MinCostFlow;
V_id src, dst;
Flow flw, cap;
Cost cst;
E_id rev;
public:
Edge() = default;
Edge(const V_id src, const V_id dst, const Flow cap, const Cost cst,
const E_id rev)
: src(src), dst(dst), flw(0), cap(cap), cst(cst), rev(rev) {}
[[nodiscard]] Flow residual_cap() const { return cap - flw; }
};
public:
class EdgePtr {
friend class MinCostFlow;
const MinCostFlow *obj;
V_id v;
E_id e;
EdgePtr(const MinCostFlow *const obj, const V_id v, const E_id e)
: obj(obj), v(v), e(e) {}
[[nodiscard]] const Edge &edge() const { return obj->g[v][e]; }
[[nodiscard]] const Edge &rev() const {
const Edge &e = edge();
return obj->g[e.dst][e.rev];
}
public:
EdgePtr() = default;
[[nodiscard]] V_id src() const { return v; }
[[nodiscard]] V_id dst() const { return edge().dst; }
[[nodiscard]] Flow flow() const { return edge().flw; }
[[nodiscard]] Flow lower() const { return -rev().cap; }
[[nodiscard]] Flow upper() const { return edge().cap; }
[[nodiscard]] Cost cost() const { return edge().cst; }
[[nodiscard]] Cost gain() const { return -edge().cst; }
};
private:
V_id n;
std::vector<std::vector<Edge>> g;
std::vector<Flow> b;
public:
MinCostFlow() : n(0) {}
V_id add_vertex() {
++n;
g.resize(n);
b.resize(n);
return n - 1;
}
std::vector<V_id> add_vertices(const size_t siz) {
std::vector<V_id> ret(siz);
std::iota(std::begin(ret), std::end(ret), n);
n += siz;
g.resize(n);
b.resize(n);
return ret;
}
EdgePtr add_edge(const V_id src, const V_id dst, const Flow lo, const Flow hi,
const Cost cst) {
assert(src < n && dst < n);
const E_id e = g[src].size(), re = src == dst ? e + 1 : g[dst].size();
assert(lo <= hi);
g[src].emplace_back(Edge{src, dst, hi, cst * objective, re});
g[dst].emplace_back(Edge{dst, src, -lo, -cst * objective, e});
return EdgePtr{this, src, e};
}
void add_supply(const V_id v, const Flow amt) { b[v] += amt; }
void add_demand(const V_id v, const Flow amt) { b[v] -= amt; }
private:
// Variables used in calculation
const Cost inf = std::numeric_limits<Cost>::max();
Cost far;
std::vector<Cost> pot;
std::vector<Cost> dis;
std::vector<Edge *> fa; // out-forrest.
using Item = std::pair<Cost, V_id>;
struct Heap : std::priority_queue<Item, std::vector<Item>, std::greater<>> {
void clear() { this->c.clear(); }
};
Heap pq; // Empty between searches; retain the allocation for the next one.
std::vector<V_id> ev, dv;
void enqueue(Cost distance, V_id vertex) { pq.emplace(distance, vertex); }
void push(Edge &e, const Flow amt) {
e.flw += amt;
g[e.dst][e.rev].flw -= amt;
}
Cost residual_cost(const V_id src, const V_id dst, const Edge &e) {
return e.cst + pot[src] - pot[dst];
}
bool dual(const Flow dlt) {
ev.erase(std::remove_if(std::begin(ev), std::end(ev),
[&](const V_id v) { return b[v] < dlt; }),
std::end(ev));
dv.erase(std::remove_if(std::begin(dv), std::end(dv),
[&](const V_id v) { return b[v] > -dlt; }),
std::end(dv));
if (ev.empty() || dv.empty()) return false;
dis.assign(n, inf);
fa.assign(n, nullptr);
for (const auto v : ev)
enqueue(dis[v] = 0, v);
far = 0;
std::size_t nd = 0;
while (!pq.empty()) {
auto [d, u] = pq.top();
pq.pop();
if (dis[u] < d)
continue;
far = d;
if (b[u] <= -dlt)
++nd;
if (nd >= dv.size())
break;
for (auto &e : g[u]) {
if (e.residual_cap() < dlt)
continue;
const auto v = e.dst;
const auto vtx = d + residual_cost(u, v, e);
if (vtx >= dis[v])
continue;
enqueue(dis[v] = vtx, v);
fa[v] = &e;
}
}
pq.clear();
for (V_id v = 0; v < n; ++v) {
pot[v] += std::min(dis[v], far);
}
return nd > 0;
}
void primal(const Flow dlt) {
for (const auto t : dv) {
if (dis[t] > far)
continue;
Flow f = -b[t];
V_id v;
for (v = t; fa[v] != nullptr; v = fa[v]->src) {
f = std::min(f, fa[v]->residual_cap());
}
f = std::min(f, b[v]);
f -= f % dlt;
if (f <= 0)
continue;
for (v = t; fa[v] != nullptr;) {
auto &e = *fa[v];
push(e, f);
V_id u = fa[v]->src;
if (e.residual_cap() <= 0)
fa[v] = nullptr;
v = u;
}
b[t] += f;
b[v] -= f;
}
}
void saturate_negative(const Flow dlt) {
ev.clear();
dv.clear();
for (auto &es : g)
for (auto &e : es) {
Flow ca1 = e.residual_cap();
ca1 -= ca1 % dlt;
if (ca1 == 0) continue;
const Cost cs1 = residual_cost(e.src, e.dst, e);
if (cs1 < 0 || ca1 < 0) {
push(e, ca1);
if (e.src != e.dst) {
b[e.src] -= ca1;
b[e.dst] += ca1;
}
}
}
for (V_id v = 0; v < n; ++v)
if (b[v] != 0) {
(b[v] > 0 ? ev : dv).emplace_back(v);
}
}
public:
/// @brief Solve using signed integral flows. Flow must hold residuals,
/// excesses and their negations; Cost must hold reduced costs/distances.
/// Choose Total wide enough for all objective products and partial sums,
/// e.g. solve<__int128>() when edge costs fit int64 but their total does not.
template <class Total = Cost> std::pair<Status, Total> solve() {
pot.resize(n);
Flow in1 = 1;
for (const auto t : b)
in1 = std::max({in1, t, -t});
for (const auto &es : g)
for (const auto &e : es)
in1 = std::max({in1, e.residual_cap(), -e.residual_cap()});
Flow dlt = 1;
while (dlt <= in1 / SCALING_FACTOR)
dlt *= SCALING_FACTOR;
for (; dlt; dlt /= SCALING_FACTOR) {
saturate_negative(dlt);
while (dual(dlt))
primal(dlt);
}
Total val = get_result_value<Total>();
if (ev.empty() && dv.empty()) {
return {Status::OPTIMAL, val / Total(objective)};
} else {
return {Status::INFEASIBLE, val / Total(objective)};
}
}
std::vector<Cost> get_potential() {
assert(pot.size() == n); // Call after solve().
if (n == 0) return pot;
auto [low, high] = std::minmax_element(pot.begin(), pot.end());
if (*low == *high) {
std::fill(pot.begin(), pot.end(), Cost(0));
return pot;
}
// Canonical shortest distances from a new zero-cost source to every
// vertex. Reweight by the feasible pot so residual edges are nonnegative.
Cost top = *high;
std::vector<Cost> distance(n);
for (V_id v = 0; v < n; ++v) enqueue(distance[v] = top - pot[v], v);
while (!pq.empty()) {
auto [d, u] = pq.top();
pq.pop();
if (distance[u] < d) continue;
for (const auto &e : g[u]) if (e.residual_cap() > 0) {
Cost cost = residual_cost(u, e.dst, e);
assert(cost >= 0);
if (d >= distance[e.dst] || cost >= distance[e.dst] - d) continue;
distance[e.dst] = d + cost;
enqueue(distance[e.dst], e.dst);
}
}
for (V_id v = 0; v < n; ++v) pot[v] -= top - distance[v];
return pot;
}
// Internal minimization objective: negated user objective for MAXIMIZE.
template <class T> T get_result_value() {
T val = 0;
for (const auto &es : g)
for (std::size_t i = 0; i < es.size(); ++i) {
const auto &e = es[i];
// The two residual directions have the same flow*cost product.
// Count one, avoiding an unnecessary doubled intermediate total.
if (e.src < e.dst || (e.src == e.dst && i < e.rev))
val += T(e.flw) * T(e.cst);
}
return val;
}
std::vector<size_t> get_cut() {
std::vector<size_t> res;
if (ev.empty())
return res;
for (size_t v = 0; v < n; ++v) {
if (dv.empty() || (dis[v] < inf))
res.emplace_back(v);
}
return res;
}
};
} // namespace noya
#endif // NOYA_MIN_COST_B_FLOW_HPP
// SPDX-License-Identifier: Apache-2.0
// Adapted from Library Checker, with changes for libra.
// Sources and license: third_party/README.md.
/// @complexity Time: O((V+E)^2 log(V+E+1) log(U+1)) conservative scaling bound
/// for a fixed scaling factor; get_potential: O((V+E) log(V+E+1)).
/// U bounds absolute initial residual capacities/balances. Space: O(V+E).
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <functional>
#include <limits>
#include <numeric>
#include <queue>
#include <tuple>
#include <utility>
#include <vector>
namespace noya {
/// @brief Minimum- or maximum-cost feasible b-flow with lower and upper bounds.
/// Lower bounds are represented as reverse residual capacity and vertex
/// supplies are maintained as excesses. Capacity scaling repeatedly saturates
/// negative reduced-cost residual edges, then a multi-source Dijkstra dual
/// step updates potentials and a primal step routes multiples of the scale from
/// excess vertices to deficits. When the scale reaches one, residual reduced
/// costs certify optimality; a remaining excess certifies infeasibility.
enum Objective {
MINIMIZE = 1,
MAXIMIZE = -1,
};
enum class Status {
OPTIMAL,
INFEASIBLE,
};
template <class Flow, class Cost, Objective objective = Objective::MINIMIZE,
Flow SCALING_FACTOR = 2>
class MinCostFlow {
static_assert(SCALING_FACTOR >= 2);
static_assert(std::numeric_limits<Flow>::is_integer &&
std::numeric_limits<Flow>::is_signed);
using V_id = uint32_t;
using E_id = uint32_t;
class Edge {
friend class MinCostFlow;
V_id src, dst;
Flow flw, cap;
Cost cst;
E_id rev;
public:
Edge() = default;
Edge(const V_id src, const V_id dst, const Flow cap, const Cost cst,
const E_id rev)
: src(src), dst(dst), flw(0), cap(cap), cst(cst), rev(rev) {}
[[nodiscard]] Flow residual_cap() const { return cap - flw; }
};
public:
class EdgePtr {
friend class MinCostFlow;
const MinCostFlow *obj;
V_id v;
E_id e;
EdgePtr(const MinCostFlow *const obj, const V_id v, const E_id e)
: obj(obj), v(v), e(e) {}
[[nodiscard]] const Edge &edge() const { return obj->g[v][e]; }
[[nodiscard]] const Edge &rev() const {
const Edge &e = edge();
return obj->g[e.dst][e.rev];
}
public:
EdgePtr() = default;
[[nodiscard]] V_id src() const { return v; }
[[nodiscard]] V_id dst() const { return edge().dst; }
[[nodiscard]] Flow flow() const { return edge().flw; }
[[nodiscard]] Flow lower() const { return -rev().cap; }
[[nodiscard]] Flow upper() const { return edge().cap; }
[[nodiscard]] Cost cost() const { return edge().cst; }
[[nodiscard]] Cost gain() const { return -edge().cst; }
};
private:
V_id n;
std::vector<std::vector<Edge>> g;
std::vector<Flow> b;
public:
MinCostFlow() : n(0) {}
V_id add_vertex() {
++n;
g.resize(n);
b.resize(n);
return n - 1;
}
std::vector<V_id> add_vertices(const size_t siz) {
std::vector<V_id> ret(siz);
std::iota(std::begin(ret), std::end(ret), n);
n += siz;
g.resize(n);
b.resize(n);
return ret;
}
EdgePtr add_edge(const V_id src, const V_id dst, const Flow lo, const Flow hi,
const Cost cst) {
assert(src < n && dst < n);
const E_id e = g[src].size(), re = src == dst ? e + 1 : g[dst].size();
assert(lo <= hi);
g[src].emplace_back(Edge{src, dst, hi, cst * objective, re});
g[dst].emplace_back(Edge{dst, src, -lo, -cst * objective, e});
return EdgePtr{this, src, e};
}
void add_supply(const V_id v, const Flow amt) { b[v] += amt; }
void add_demand(const V_id v, const Flow amt) { b[v] -= amt; }
private:
// Variables used in calculation
const Cost inf = std::numeric_limits<Cost>::max();
Cost far;
std::vector<Cost> pot;
std::vector<Cost> dis;
std::vector<Edge *> fa; // out-forrest.
using Item = std::pair<Cost, V_id>;
struct Heap : std::priority_queue<Item, std::vector<Item>, std::greater<>> {
void clear() { this->c.clear(); }
};
Heap pq; // Empty between searches; retain the allocation for the next one.
std::vector<V_id> ev, dv;
void enqueue(Cost distance, V_id vertex) { pq.emplace(distance, vertex); }
void push(Edge &e, const Flow amt) {
e.flw += amt;
g[e.dst][e.rev].flw -= amt;
}
Cost residual_cost(const V_id src, const V_id dst, const Edge &e) {
return e.cst + pot[src] - pot[dst];
}
bool dual(const Flow dlt) {
ev.erase(std::remove_if(std::begin(ev), std::end(ev),
[&](const V_id v) { return b[v] < dlt; }),
std::end(ev));
dv.erase(std::remove_if(std::begin(dv), std::end(dv),
[&](const V_id v) { return b[v] > -dlt; }),
std::end(dv));
if (ev.empty() || dv.empty()) return false;
dis.assign(n, inf);
fa.assign(n, nullptr);
for (const auto v : ev)
enqueue(dis[v] = 0, v);
far = 0;
std::size_t nd = 0;
while (!pq.empty()) {
auto [d, u] = pq.top();
pq.pop();
if (dis[u] < d)
continue;
far = d;
if (b[u] <= -dlt)
++nd;
if (nd >= dv.size())
break;
for (auto &e : g[u]) {
if (e.residual_cap() < dlt)
continue;
const auto v = e.dst;
const auto vtx = d + residual_cost(u, v, e);
if (vtx >= dis[v])
continue;
enqueue(dis[v] = vtx, v);
fa[v] = &e;
}
}
pq.clear();
for (V_id v = 0; v < n; ++v) {
pot[v] += std::min(dis[v], far);
}
return nd > 0;
}
void primal(const Flow dlt) {
for (const auto t : dv) {
if (dis[t] > far)
continue;
Flow f = -b[t];
V_id v;
for (v = t; fa[v] != nullptr; v = fa[v]->src) {
f = std::min(f, fa[v]->residual_cap());
}
f = std::min(f, b[v]);
f -= f % dlt;
if (f <= 0)
continue;
for (v = t; fa[v] != nullptr;) {
auto &e = *fa[v];
push(e, f);
V_id u = fa[v]->src;
if (e.residual_cap() <= 0)
fa[v] = nullptr;
v = u;
}
b[t] += f;
b[v] -= f;
}
}
void saturate_negative(const Flow dlt) {
ev.clear();
dv.clear();
for (auto &es : g)
for (auto &e : es) {
Flow ca1 = e.residual_cap();
ca1 -= ca1 % dlt;
if (ca1 == 0) continue;
const Cost cs1 = residual_cost(e.src, e.dst, e);
if (cs1 < 0 || ca1 < 0) {
push(e, ca1);
if (e.src != e.dst) {
b[e.src] -= ca1;
b[e.dst] += ca1;
}
}
}
for (V_id v = 0; v < n; ++v)
if (b[v] != 0) {
(b[v] > 0 ? ev : dv).emplace_back(v);
}
}
public:
/// @brief Solve using signed integral flows. Flow must hold residuals,
/// excesses and their negations; Cost must hold reduced costs/distances.
/// Choose Total wide enough for all objective products and partial sums,
/// e.g. solve<__int128>() when edge costs fit int64 but their total does not.
template <class Total = Cost> std::pair<Status, Total> solve() {
pot.resize(n);
Flow in1 = 1;
for (const auto t : b)
in1 = std::max({in1, t, -t});
for (const auto &es : g)
for (const auto &e : es)
in1 = std::max({in1, e.residual_cap(), -e.residual_cap()});
Flow dlt = 1;
while (dlt <= in1 / SCALING_FACTOR)
dlt *= SCALING_FACTOR;
for (; dlt; dlt /= SCALING_FACTOR) {
saturate_negative(dlt);
while (dual(dlt))
primal(dlt);
}
Total val = get_result_value<Total>();
if (ev.empty() && dv.empty()) {
return {Status::OPTIMAL, val / Total(objective)};
} else {
return {Status::INFEASIBLE, val / Total(objective)};
}
}
std::vector<Cost> get_potential() {
assert(pot.size() == n); // Call after solve().
if (n == 0) return pot;
auto [low, high] = std::minmax_element(pot.begin(), pot.end());
if (*low == *high) {
std::fill(pot.begin(), pot.end(), Cost(0));
return pot;
}
// Canonical shortest distances from a new zero-cost source to every
// vertex. Reweight by the feasible pot so residual edges are nonnegative.
Cost top = *high;
std::vector<Cost> distance(n);
for (V_id v = 0; v < n; ++v) enqueue(distance[v] = top - pot[v], v);
while (!pq.empty()) {
auto [d, u] = pq.top();
pq.pop();
if (distance[u] < d) continue;
for (const auto &e : g[u]) if (e.residual_cap() > 0) {
Cost cost = residual_cost(u, e.dst, e);
assert(cost >= 0);
if (d >= distance[e.dst] || cost >= distance[e.dst] - d) continue;
distance[e.dst] = d + cost;
enqueue(distance[e.dst], e.dst);
}
}
for (V_id v = 0; v < n; ++v) pot[v] -= top - distance[v];
return pot;
}
// Internal minimization objective: negated user objective for MAXIMIZE.
template <class T> T get_result_value() {
T val = 0;
for (const auto &es : g)
for (std::size_t i = 0; i < es.size(); ++i) {
const auto &e = es[i];
// The two residual directions have the same flow*cost product.
// Count one, avoiding an unnecessary doubled intermediate total.
if (e.src < e.dst || (e.src == e.dst && i < e.rev))
val += T(e.flw) * T(e.cst);
}
return val;
}
std::vector<size_t> get_cut() {
std::vector<size_t> res;
if (ev.empty())
return res;
for (size_t v = 0; v < n; ++v) {
if (dv.empty() || (dis[v] < inf))
res.emplace_back(v);
}
return res;
}
};
} // namespace noya
// NOYA_MIN_COST_B_FLOW_HPP