general_weighted_matching.hpp¶
求一般图的最大权匹配；适合非二分图中边有收益、每点至多匹配一次的问题。
Complexity: Time: O(n m log n). Space: O(n+m).
AC 记录：general_weighted_matching。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(n m log n). Space: O(n+m).
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <functional>
#include <queue>
#include <tuple>
#include <utility>
#include <vector>
namespace noya {
namespace general_weighted_matching_internal {
using namespace std;
template <typename CostType, typename TotalCostType>
class gabow_maximum_weighted_matching {
public:
using cost_t = CostType;
using tcost_t = TotalCostType;
private:
enum Label { kSeparated = -2, kInner = -1, kFree = 0, kOuter = 1 };
static constexpr cost_t Inf = cost_t(1) << (sizeof(cost_t) * 8 - 2);
private:
template <typename T>
class BinaryHeap {
public:
struct Node {
#ifdef FAST_COMPARE
bool operator < (const Node& rhs) const { return val < rhs.val; }
#else
bool operator < (const Node& rhs) const {
if (val < rhs.val) { return true; }
if (rhs.val < val) { return false; }
return id < rhs.id;
}
#endif
T val; int id;
};
BinaryHeap() {}
BinaryHeap(int N) : sz_(0), nd(N + 1), idx(N, 0) {}
int size() const { return sz_; }
bool empty() const { return sz_ == 0; }
void clear() { while (sz_ > 0) idx[nd[sz_--].id] = 0; }
T min() const { return nd[1].val; }
int argmin() const { return nd[1].id; } // argmin ?
T get_val(int id) const { return nd[idx[id]].val; }
void pop() { if (sz_ > 0) pop(1); }
void erase(int id) { if (idx[id]) pop(idx[id]); }
bool has(int id) const { return idx[id] != 0; }
void update(int id, T v) {
if (!has(id)) return push(id, v);
bool up = (v < nd[idx[id]].val);
nd[idx[id]].val = v;
if (up) up_heap(idx[id]);
else down_heap(idx[id]);
}
void decrease_key(int id, T v) {
if (!has(id)) return push(id, v);
if (v < nd[idx[id]].val) nd[idx[id]].val = v, up_heap(idx[id]);
}
void push(int id, T v) {
// assert(!has(id));
idx[id] = ++sz_; nd[sz_] = {v, id};
up_heap(sz_);
}
private:
void pop(int pos) {
idx[nd[pos].id] = 0;
if (pos == sz_) { --sz_; return; }
bool up = (nd[sz_].val < nd[pos].val);
nd[pos] = nd[sz_--]; idx[nd[pos].id] = pos;
if (up) up_heap(pos);
else down_heap(pos);
}
void swap_node(int a, int b) {
swap(nd[a], nd[b]); idx[nd[a].id] = a; idx[nd[b].id] = b;
}
void down_heap(int pos) {
for (int k = pos, nk = k; 2 * k <= sz_; k = nk) {
if (nd[2 * k] < nd[nk]) nk = 2 * k;
if (2 * k + 1 <= sz_ && nd[2 * k + 1] < nd[nk]) nk = 2 * k + 1;
if (nk == k) break;
swap_node(k, nk);
}
}
void up_heap(int pos) {
for (int k = pos; k > 1 && nd[k] < nd[k >> 1]; k >>= 1) swap_node(k, k >> 1);
}
int sz_;
vector<Node> nd;
vector<int> idx;
};
template <typename Key>
class PairingHeaps {
private:
struct Node {
Node() : pre(-1) {} // "pre < 0" means the node is unused.
Node(Key v) : key(v), v1(0), nxt(0), pre(0) {}
Key key; int v1, nxt, pre;
};
public:
PairingHeaps(int H, int N) : hp(H), nd(N) {
// It consists of `H` Pairing heaps.
// Each heap-node ID can appear at most 1 time(s) among heaps
// and should be in [1, N).
}
void clear(int h) { if (hp[h]) clear_rec(hp[h]), hp[h] = 0; }
void clear_all() {
for (size_t i = 0; i < hp.size(); ++i) hp[i] = 0;
for (size_t i = 0; i < nd.size(); ++i) nd[i] = Node();
}
bool empty(int h) const { return !hp[h]; }
bool used(int v) const { return nd[v].pre >= 0; }
Key min(int h) const { return nd[hp[h]].key; }
int argmin(int h) const { return hp[h]; }
void pop(int h) {
// assert(!empty(h));
erase(h, hp[h]);
}
void push(int h, int v, Key key) {
// assert(!used(v));
nd[v] = Node(key);
hp[h] = merge(hp[h], v);
}
void erase(int h, int v) {
if (!used(v)) return;
int w = two_pass_pairing(nd[v].v1);
if (!nd[v].pre) hp[h] = w;
else {
cut(v);
hp[h] = merge(hp[h], w);
}
nd[v].pre = -1;
}
void decrease_key(int h, int v, Key key) {
if (!used(v)) return push(h, v, key);
if (!nd[v].pre) nd[v].key = key;
else {
cut(v); nd[v].key = key;
hp[h] = merge(hp[h], v);
}
}
private:
void clear_rec(int v) {
for (; v; v = nd[v].nxt) {
if (nd[v].v1) clear_rec(nd[v].v1);
nd[v].pre = -1;
}
}
inline void cut(int v) {
auto& n = nd[v]; int pv = n.pre, nv = n.nxt;
auto& pn = nd[pv];
if (pn.v1 == v) pn.v1 = nv;
else pn.nxt = nv;
nd[nv].pre = pv;
n.nxt = n.pre = 0;
}
int merge(int l, int r) {
if (!l) return r;
if (!r) return l;
if (nd[l].key > nd[r].key) swap(l, r);
int lc = nd[r].nxt = nd[l].v1;
nd[l].v1 = nd[lc].pre = r;
return nd[r].pre = l;
}
int two_pass_pairing(int rt) {
if (!rt) return 0;
int a = rt; rt = 0;
while (a) {
int b = nd[a].nxt, na = 0;
nd[a].pre = nd[a].nxt = 0;
if (b) na = nd[b].nxt, nd[b].pre = nd[b].nxt = 0;
a = merge(a, b);
nd[a].nxt = rt; rt = a; a = na;
}
int s = nd[rt].nxt; nd[rt].nxt = 0;
while (s) {
int t = nd[s].nxt; nd[s].nxt = 0;
rt = merge(rt, s);
s = t;
}
return rt;
}
private:
vector<int> hp;
vector<Node> nd;
};
template <typename T>
struct PriorityQueue : public priority_queue< T, vector<T>, greater<T>> {
PriorityQueue() {}
PriorityQueue(int N) { this->c.reserve(N);}
T min() const { return this->top(); }
void clear() { this->c.clear(); }
};
template <typename T>
struct Queue {
Queue() {}
Queue(int N) : qh(0), qt(0), dat(N) {}
T operator [] (int i) const { return dat[i]; }
void enqueue(int u) { dat[qt++] = u; }
int dequeue() { return dat[qh++]; }
bool empty() const { return qh == qt; }
void clear() { qh = qt = 0; }
int size() const { return qt; }
int qh, qt;
vector<T> dat;
};
public:
struct InputEdge { int u1, to; cost_t cst; };
private:
template <typename T> using ModifiableHeap = BinaryHeap<T>;
template <typename T> using ModifiableHeaps = PairingHeaps<T>;
template <typename T> using FastHeap = PriorityQueue<T>;
struct Edge { int to; cost_t cst; };
struct Link { int u1, to; };
struct Node {
struct NodeLink { int b, v; };
Node() {}
Node(int u) : fa(0), siz(1) { lnk[0] = lnk[1] = {u, u}; }
int next_v() const { return lnk[0].v; }
int next_b() const { return lnk[0].b; }
int prev_v() const { return lnk[1].v; }
int prev_b() const { return lnk[1].b; }
int fa, siz;
NodeLink lnk[2];
};
struct Event {
Event() {}
Event(cost_t tm, int id) : tm(tm), id(id) {}
#ifdef FAST_COMPARE
bool operator < (const Event& rhs) const { return tm < rhs.tm; }
#else
bool operator < (const Event& rhs) const {
if (tm < rhs.tm) { return true; }
if (rhs.tm < tm) { return false; }
return id < rhs.id;
}
#endif
bool operator > (const Event& rhs) const { return rhs.operator<(*this); }
cost_t tm; int id;
};
struct EdgeEvent {
EdgeEvent() {}
EdgeEvent(cost_t tm, int u1, int to) : tm(tm), u1(u1), to(to) {}
#ifdef FAST_COMPARE
bool operator < (const EdgeEvent& rhs) const { return tm < rhs.tm; }
#else
bool operator < (const EdgeEvent& rhs) const {
if (tm < rhs.tm) {
return true;
}
if (tm > rhs.tm) {
return false;
}
return make_pair(u1, to) < make_pair(rhs.u1, rhs.to);
}
#endif
bool operator > (const EdgeEvent& rhs) const { return rhs.operator<(*this); }
cost_t tm; int u1, to;
};
public:
gabow_maximum_weighted_matching(int N, const vector<InputEdge>& in)
: N(N), B((N - 1) / 2), S(N + B + 1), ofs(N + 2), es(in.size() * 2),
h2(S), hps(S, S), h3(es.size()), h4(S) {
for (auto& e : in) ofs[e.u1 + 1]++, ofs[e.to + 1]++;
for (int i = 1; i <= N + 1; ++i) ofs[i] += ofs[i - 1];
for (auto& e : in) {
es[ofs[e.u1]++] = {e.to, e.cst * 2};
es[ofs[e.to]++] = {e.u1, e.cst * 2};
}
for (int i = N + 1; i > 0; --i) ofs[i] = ofs[i - 1];
ofs[0] = 0;
}
tcost_t maximum_weighted_matching(vector<pair<int, int>> &mat, bool ini=false) {
initialize();
set_potential();
if (ini) find_maximal_matching();
for (int u = 1; u <= N; ++u) if (!mt[u]) do_edmonds_search(u);
tcost_t ret = compute_optimal_value();
mat.clear();
for (int u = 1; u <= N; ++u) if (mt[u] > u) {
mat.push_back({u, mt[u]});
}
return ret;
}
private:
tcost_t compute_optimal_value() const {
tcost_t ret = 0;
for (int u = 1; u <= N; ++u) if (mt[u] > u) {
cost_t mx = 0;
for (int eid = ofs[u]; eid < ofs[u + 1]; ++eid) {
if (es[eid].to == mt[u]) mx = max(mx, es[eid].cst);
}
ret += mx;
}
return ret >> 1;
}
inline tcost_t reduced_cost(int u, int v, const Edge& e) const {
return tcost_t(pot[u]) + pot[v] - e.cst;
}
void rematch(int v, int w) {
int t = mt[v]; mt[v] = w;
if (mt[t] != v) return;
if (lnk[v].to == sf1[lnk[v].to]) {
mt[t] = lnk[v].u1;
rematch(mt[t], t);
} else {
int x = lnk[v].u1, y = lnk[v].to;
rematch(x, y); rematch(y, x);
}
}
void fix_mate_and_base(int b) {
if (b <= N) return;
int bv = bas[b], mv = nd[bv].lnk[0].v, bmv = nd[bv].lnk[0].b;
int d = (nd[bmv].lnk[1].v == mt[mv]) ? 0 : 1;
while (1) {
int mv = nd[bv].lnk[d].v, bmv = nd[bv].lnk[d].b;
if (nd[bmv].lnk[1 ^ d].v != mt[mv]) break;
fix_mate_and_base(bv); fix_mate_and_base(bmv);
bv = nd[bmv].lnk[d].b;
}
fix_mate_and_base(bas[b] = bv);
mt[b] = mt[bv];
}
void reset_time() {
tc1 = 0; ev1 = {Inf, 0};
}
void reset_blossom(int b) {
lab[b] = kFree; lnk[b].u1 = 0; sl[b] = Inf; lz[b] = 0;
}
void reset_all() {
lab[0] = kFree; lnk[0].u1 = 0;
for (int v = 1; v <= N; ++v) { // should be optimized for sparse graphs.
if (lab[v] == kOuter) pot[v] -= tc1;
else {
int bv = sf1[v];
pot[v] += lz[bv];
if (lab[bv] == kInner) pot[v] += tc1 - tc[bv];
}
reset_blossom(v);
}
for (int b = N + 1, r = B - ubi; r > 0 && b < S; ++b) if (bas[b] != b) {
if (sf1[b] == b) {
fix_mate_and_base(b);
if (lab[b] == kOuter) pot[b] += (tc1 - tc[b]) << 1;
else if (lab[b] == kInner) fix_blossom_potential<kInner>(b);
else fix_blossom_potential<kFree>(b);
}
hps.clear(b);
reset_blossom(b); --r;
}
que.clear();
reset_time(); h2.clear();
h3.clear(); h4.clear();
}
void do_edmonds_search(int rt) {
if (pot[rt] == 0) return;
link_blossom(sf1[rt], {0, 0});
push_outer_and_fix_potentials(sf1[rt], 0);
for (bool aug = false; !aug; ) {
aug = augment(rt);
if (aug) break;
aug = adjust_dual_variables(rt);
}
reset_all();
}
template <Label Lab>
inline cost_t fix_blossom_potential(int b) {
// Return the amount.
// (If v is an atom, the pot[v] will not be changed.)
cost_t d = lz[b]; lz[b] = 0;
if (Lab == kInner) {
cost_t dt = tc1 - tc[b];
if (b > N) pot[b] -= dt << 1;
d += dt;
}
return d;
}
template <Label Lab>
inline void update_heap2(int x, int y, int by, cost_t t) {
if (t >= sl[y]) return;
sl[y] = t; bf[y] = x;
if (y == by) {
if (Lab != kInner) h2.decrease_key(y, EdgeEvent(t + lz[y], x, y));
} else {
int gy = grp[y];
if (gy != y) {
if (t >= sl[gy]) return;
sl[gy] = t;
}
hps.decrease_key(by, gy, EdgeEvent(t, x, y));
if (Lab == kInner) return;
EdgeEvent m = hps.min(by);
h2.decrease_key(by, EdgeEvent(m.tm + lz[by], m.u1, m.to));
}
}
void activate_heap2_node(int b) {
if (b <= N) {
if (sl[b] < Inf) h2.push(b, EdgeEvent(sl[b] + lz[b], bf[b], b));
} else {
if (hps.empty(b)) return;
EdgeEvent m = hps.min(b);
h2.push(b, EdgeEvent(m.tm + lz[b], m.u1, m.to));
}
}
void swap_blossom(int a, int b) {
// Assume that `b` is a maximal blossom.
swap(bas[a], bas[b]); if (bas[a] == a) bas[a] = b;
swap(son[a], son[b]); if (son[a] == a) son[a] = b;
swap(lnk[a], lnk[b]);
swap(mt[a], mt[b]);
swap(pot[a], pot[b]); swap(lz[a], lz[b]);
swap(tc[a], tc[b]);
for (int d = 0; d < 2; ++d) nd[nd[a].lnk[d].b].lnk[1 ^ d].b = b;
swap(nd[a], nd[b]);
}
void set_surface_and_group(int b, int sf, int g) {
sf1[b] = sf, grp[b] = g;
if (b <= N) return;
for (int bb = bas[b]; sf1[bb] != sf; bb = nd[bb].next_b()) {
set_surface_and_group(bb, sf, g);
}
}
void merge_smaller_blossoms(int bid) {
int lb = bid, ms = 1;
for (int bta = bas[bid], b = bta; ;) {
if (nd[b].siz > ms) ms = nd[b].siz, lb = b;
if ((b = nd[b].next_b()) == bta) break;
}
for (int bta = bas[bid], b = bta; ;) {
if (b != lb) set_surface_and_group(b, lb, b);
if ((b = nd[b].next_b()) == bta) break;
}
grp[lb] = lb;
if (ms > 1) {
sf1[bid] = son[bid] = lb;
swap_blossom(lb, bid);
} else son[bid] = 0;
}
void contract(int x, int y, int eid) {
int bx = sf1[x], by = sf1[y]; assert(bx != by);
const int h = -(eid + 1);
lnk[sf1[mt[bx]]].u1 = lnk[sf1[mt[by]]].u1 = h;
int lca = -1;
while (1) {
if (mt[by] != 0) swap(bx, by);
bx = lca = sf1[lnk[bx].u1];
if (lnk[sf1[mt[bx]]].u1 == h) break;
lnk[sf1[mt[bx]]].u1 = h;
}
const int bid = ub[--ubi]; assert(ubi >= 0);
int tsz = 0;
for (int d = 0; d < 2; ++d) {
for (int bv = sf1[x]; bv != lca; ) {
int mv = mt[bv], bmv = sf1[mv], v = mt[mv];
int f = lnk[v].u1, t = lnk[v].to;
tsz += nd[bv].siz + nd[bmv].siz;
lnk[mv] = {x, y};
if (bv > N) pot[bv] += (tc1 - tc[bv]) << 1;
if (bmv > N) h4.erase(bmv);
push_outer_and_fix_potentials(bmv, fix_blossom_potential<kInner>(bmv));
nd[bv].lnk[d] = {bmv, mv};
nd[bmv].lnk[1 ^ d] = {bv, v}; nd[bmv].lnk[d] = {bv = sf1[f], f};
nd[bv].lnk[1 ^ d] = {bmv, t};
}
nd[sf1[x]].lnk[1 ^ d] = {sf1[y], y};
swap(x, y);
}
if (lca > N) pot[lca] += (tc1 - tc[lca]) << 1;
nd[bid].siz = tsz + nd[lca].siz;
bas[bid] = lca; lnk[bid] = lnk[lca]; mt[bid] = mt[lca];
lab[bid] = kOuter;
sf1[bid] = bid; tc[bid] = tc1;
pot[bid] = 0; lz[bid] = 0;
merge_smaller_blossoms(bid); // O(n log n) time / Edmonds search
}
void link_blossom(int v, Link l) {
lnk[v] = {l.u1, l.to};
if (v <= N) return;
int b = bas[v]; link_blossom(b, l);
int pb = nd[b].prev_b();
l = {nd[pb].next_v(), nd[b].prev_v()};
for (int bv = b; ; ) {
int bw = nd[bv].next_b();
if (bw == b) break;
link_blossom(bw, l);
Link nl = {nd[bw].prev_v(), nd[bv].next_v()};
bv = nd[bw].next_b();
link_blossom(bv, nl);
}
}
void push_outer_and_fix_potentials(int v, cost_t d) {
lab[v] = kOuter;
if (v > N) {
for (int b = bas[v]; lab[b] != kOuter; b = nd[b].next_b()) {
push_outer_and_fix_potentials(b, d);
}
} else {
pot[v] += tc1 + d;
if (pot[v] < ev1.tm) ev1 = {pot[v], v};
que.enqueue(v);
}
}
bool grow(int x, int y) {
int by = sf1[y];
bool vis = (lab[by] != kFree);
if (!vis) link_blossom(by, {0, 0});
lab[by] = kInner; tc[by] = tc1; h2.erase(by);
if (y != by) h4.update(by, tc1 + (pot[by] >> 1));
int z = mt[by];
if (z == 0) {
rematch(x, y); rematch(y, x);
return true;
}
int bz = sf1[z];
if (!vis) link_blossom(bz, {x, y});
else lnk[bz] = lnk[z] = {x, y};
push_outer_and_fix_potentials(bz, fix_blossom_potential<kFree>(bz));
tc[bz] = tc1; h2.erase(bz);
return false;
}
void free_blossom(int bid) {
ub[ubi++] = bid;
bas[bid] = bid;
}
int recalculate_minimum_slack(int b, int g) {
// Return the destination of the best edge of blossom `g`.
if (b <= N) {
if (sl[b] >= sl[g]) return 0;
sl[g] = sl[b]; bf[g] = bf[b];
return b;
}
int v = 0;
for (int bta = bas[b], bb = bta; ; ) {
int w = recalculate_minimum_slack(bb, g);
if (w != 0) v = w;
if ((bb = nd[bb].next_b()) == bta) break;
}
return v;
}
void construct_smaller_components(int b, int sf, int g) {
sf1[b] = sf, grp[b] = g; // `grp[b] = g` is unneeded.
if (b <= N) return;
for (int bb = bas[b]; sf1[bb] != sf; bb = nd[bb].next_b()) {
if (bb == son[b]) {
construct_smaller_components(bb, sf, g);
} else {
set_surface_and_group(bb, sf, bb);
int to = 0;
if (bb > N) sl[bb] = Inf, to = recalculate_minimum_slack(bb, bb);
else if (sl[bb] < Inf) to = bb;
if (to > 0) hps.push(sf, bb, EdgeEvent(sl[bb], bf[bb], to));
}
}
}
void move_to_largest_blossom(int bid) {
const int h = son[bid];
cost_t d = (tc1 - tc[bid]) + lz[bid]; lz[bid] = 0;
for (int bta = bas[bid], b = bta; ;) {
tc[b] = tc1;
lz[b] = d;
if (b != h) construct_smaller_components(b, b, b), hps.erase(bid, b);
if ((b = nd[b].next_b()) == bta) break;
}
if (h > 0) swap_blossom(h, bid), bid = h;
free_blossom(bid);
}
void expand(int bid) {
int mv = mt[bas[bid]];
move_to_largest_blossom(bid); // O(n log n) time / Edmonds search
Link ol = lnk[mv];
int ob = sf1[mt[mv]], rt = sf1[ol.to];
int d = (mt[rt] == nd[rt].lnk[0].v) ? 1 : 0;
for (int b = nd[ob].lnk[d ^ 1].b; b != rt; ) {
lab[b] = kSeparated; activate_heap2_node(b); b = nd[b].lnk[d ^ 1].b;
lab[b] = kSeparated; activate_heap2_node(b); b = nd[b].lnk[d ^ 1].b;
}
for (int b = ob; ; b = nd[b].lnk[d].b) {
lab[b] = kInner;
int nb = nd[b].lnk[d].b;
if (b == rt) lnk[mt[b]] = ol;
else lnk[mt[b]] = {nd[b].lnk[d].v, nd[nb].lnk[d ^ 1].v};
lnk[sf1[mt[b]]] = lnk[mt[b]]; // fix tree links
if (b > N) {
if (pot[b] == 0) expand(b);
else h4.push(b, tc1 + (pot[b] >> 1));
}
if (b == rt) break;
push_outer_and_fix_potentials(nb, fix_blossom_potential<kInner>(b = nb));
}
}
bool augment(int rt) {
// Return true if an augmenting path is found.
while (!que.empty()) {
int x = que.dequeue(), bx = sf1[x];
if (pot[x] == tc1) {
if (x != rt) rematch(x, 0);
return true;
}
for (int eid = ofs[x]; eid < ofs[x + 1]; ++eid) {
auto& e = es[eid]; int y = e.to, by = sf1[y];
if (bx == by) continue;
Label l = lab[by];
if (l == kOuter) {
cost_t t = reduced_cost(x, y, e) >> 1; // < 2 * Inf
if (t == tc1) {
contract(x, y, eid); bx = sf1[x];
} else if (t < ev1.tm) {
h3.emplace(t, x, eid);
}
} else {
tcost_t t = reduced_cost(x, y, e); // < 3 * Inf
if (t >= Inf) continue;
if (l != kInner) {
if (cost_t(t) + lz[by] == tc1) {
if (grow(x, y)) return true;
} else update_heap2<kFree>(x, y, by, t);
} else {
if (mt[x] != y) update_heap2<kInner>(x, y, by, t);
}
}
}
}
return false;
}
bool adjust_dual_variables(int rt) {
// delta1 : rematch
cost_t t1 = ev1.tm;
// delta2 : grow
cost_t t2 = Inf;
if (!h2.empty()) t2 = h2.min().tm;
// delta3 : contract : O(m log n) time / Edmonds search [ bottleneck (?) ]
cost_t t3 = Inf;
while (!h3.empty()) {
EdgeEvent e = h3.min();
int x = e.u1, y = es[e.to].to; // e.to is some edge id.
if (sf1[x] != sf1[y]) {
t3 = e.tm;
break;
} else h3.pop();
}
// delta4 : expand
cost_t t4 = Inf;
if (!h4.empty()) t4 = h4.min();
// -- events --
cost_t tn = min(min(t1, t2), min(t3, t4));
assert(tc1 <= tn && tn < Inf);
tc1 = tn;
if (tc1 == ev1.tm) {
int x = ev1.id;
if (x != rt) rematch(x, 0);
return true;
}
while (!h2.empty() && h2.min().tm == tc1) {
int x = h2.min().u1, y = h2.min().to;
if (grow(x, y)) return true; // `grow` function will call `h2.erase(by)`.
}
while (!h3.empty() && h3.min().tm == tc1) {
int x = h3.min().u1, eid = h3.min().to;
int y = es[eid].to; h3.pop();
if (sf1[x] == sf1[y]) continue;
contract(x, y, eid);
}
while (!h4.empty() && h4.min() == tc1) {
int b = h4.argmin(); h4.pop();
expand(b);
}
return false;
}
private:
void initialize() {
que = Queue<int>(N);
mt.assign(S, 0);
lnk.assign(S, {0, 0});
lab.assign(S, kFree);
bas.resize(S); for (int u = 1; u < S; ++u) bas[u] = u;
sf1.resize(S); for (int u = 1; u < S; ++u) sf1[u] = u;
pot.resize(S);
nd.resize(S); for (int b = 1; b < S; ++b) nd[b] = Node(b);
ub.resize(B); for (int i = 0; i < B; ++i) ub[i] = N + B - i;
ubi = B;
// for O(nm log n) implementation
reset_time();
tc.resize(S);
sl.resize(S); for (int i = 0; i < S; ++i) sl[i] = Inf;
bf.assign(S, 0);
son.assign(S, 0);
lz.assign(S, 0);
grp.resize(S); for (int i = 0; i < S; ++i) grp[i] = i;
}
void set_potential() {
for (int u = 1; u <= N; ++u) {
cost_t mx = 0;
for (int eid = ofs[u]; eid < ofs[u + 1]; ++eid) {
mx = max(mx, es[eid].cst);
}
pot[u] = mx >> 1;
}
}
void find_maximal_matching() {
// Find a maximal matching naively.
for (int u = 1; u <= N; ++u) if (!mt[u]) {
for (int eid = ofs[u]; eid < ofs[u + 1]; ++eid) {
auto& e = es[eid]; int v = e.to;
if (mt[v] > 0 || reduced_cost(u, v, e) > 0) continue;
mt[u] = v; mt[v] = u;
break;
}
}
}
private:
const int N, B, S; // N = |V|, B = (|V| - 1) / 2, S = N + B + 1
vector<int> ofs;
vector<Edge> es;
Queue<int> que;
vector<int> mt, sf1, bas;
vector<Link> lnk;
vector<Label> lab;
vector<cost_t> pot;
vector<int> ub; int ubi;
vector<Node> nd;
// for O(nm log n) implementation
vector<int> son, grp;
vector<cost_t> tc, lz, sl;
vector<int> bf;
cost_t tc1;
Event ev1;
ModifiableHeap<EdgeEvent> h2;
ModifiableHeaps<EdgeEvent> hps;
FastHeap<EdgeEvent> h3;
ModifiableHeap<cost_t> h4;
};
} // namespace general_weighted_matching_internal
template <class TotalCost> struct weighted_matching_result {
TotalCost wt{};
std::vector<std::pair<int, int>> es;
};
/// @brief Compute a maximum-weight matching in a general undirected graph.
/// A primal-dual Edmonds search grows an alternating forest through zero
/// reduced-cost edges. Odd cycles are contracted into blossoms; dual changes
/// schedule grow, contract, and expand events. Gabow's heavy-blossom grouping
/// and meldable heaps make all slack/event maintenance near-linear per search.
template <class Cost, class TotalCost = Cost>
weighted_matching_result<TotalCost> general_weighted_matching(
int n1,
const std::vector<std::tuple<int, int, Cost>> &es1) {
using solver = general_weighted_matching_internal::
gabow_maximum_weighted_matching<Cost, TotalCost>;
std::vector<typename solver::InputEdge> es;
es.reserve(es1.size());
for (auto [a1, b1, wt] : es1) {
assert(0 <= a1 && a1 < n1);
assert(0 <= b1 && b1 < n1);
assert(a1 != b1 && wt >= Cost{});
es.push_back({a1 + 1, b1 + 1, wt});
}
solver mat(n1, es);
weighted_matching_result<TotalCost> res;
res.wt = mat.maximum_weighted_matching(res.es);
for (auto &[a1, b1] : res.es) {
--a1;
--b1;
}
return res;
}
} // namespace noya
#ifndef NOYA_GENERAL_WEIGHTED_MATCHING_HPP
#define NOYA_GENERAL_WEIGHTED_MATCHING_HPP 1
/// @complexity Time: O(n m log n). Space: O(n+m).
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <functional>
#include <queue>
#include <tuple>
#include <utility>
#include <vector>
namespace noya {
namespace general_weighted_matching_internal {
using namespace std;
template <typename CostType, typename TotalCostType>
class gabow_maximum_weighted_matching {
public:
using cost_t = CostType;
using tcost_t = TotalCostType;
private:
enum Label { kSeparated = -2, kInner = -1, kFree = 0, kOuter = 1 };
static constexpr cost_t Inf = cost_t(1) << (sizeof(cost_t) * 8 - 2);
private:
template <typename T>
class BinaryHeap {
public:
struct Node {
#ifdef FAST_COMPARE
bool operator < (const Node& rhs) const { return val < rhs.val; }
#else
bool operator < (const Node& rhs) const {
if (val < rhs.val) { return true; }
if (rhs.val < val) { return false; }
return id < rhs.id;
}
#endif
T val; int id;
};
BinaryHeap() {}
BinaryHeap(int N) : sz_(0), nd(N + 1), idx(N, 0) {}
int size() const { return sz_; }
bool empty() const { return sz_ == 0; }
void clear() { while (sz_ > 0) idx[nd[sz_--].id] = 0; }
T min() const { return nd[1].val; }
int argmin() const { return nd[1].id; } // argmin ?
T get_val(int id) const { return nd[idx[id]].val; }
void pop() { if (sz_ > 0) pop(1); }
void erase(int id) { if (idx[id]) pop(idx[id]); }
bool has(int id) const { return idx[id] != 0; }
void update(int id, T v) {
if (!has(id)) return push(id, v);
bool up = (v < nd[idx[id]].val);
nd[idx[id]].val = v;
if (up) up_heap(idx[id]);
else down_heap(idx[id]);
}
void decrease_key(int id, T v) {
if (!has(id)) return push(id, v);
if (v < nd[idx[id]].val) nd[idx[id]].val = v, up_heap(idx[id]);
}
void push(int id, T v) {
// assert(!has(id));
idx[id] = ++sz_; nd[sz_] = {v, id};
up_heap(sz_);
}
private:
void pop(int pos) {
idx[nd[pos].id] = 0;
if (pos == sz_) { --sz_; return; }
bool up = (nd[sz_].val < nd[pos].val);
nd[pos] = nd[sz_--]; idx[nd[pos].id] = pos;
if (up) up_heap(pos);
else down_heap(pos);
}
void swap_node(int a, int b) {
swap(nd[a], nd[b]); idx[nd[a].id] = a; idx[nd[b].id] = b;
}
void down_heap(int pos) {
for (int k = pos, nk = k; 2 * k <= sz_; k = nk) {
if (nd[2 * k] < nd[nk]) nk = 2 * k;
if (2 * k + 1 <= sz_ && nd[2 * k + 1] < nd[nk]) nk = 2 * k + 1;
if (nk == k) break;
swap_node(k, nk);
}
}
void up_heap(int pos) {
for (int k = pos; k > 1 && nd[k] < nd[k >> 1]; k >>= 1) swap_node(k, k >> 1);
}
int sz_;
vector<Node> nd;
vector<int> idx;
};
template <typename Key>
class PairingHeaps {
private:
struct Node {
Node() : pre(-1) {} // "pre < 0" means the node is unused.
Node(Key v) : key(v), v1(0), nxt(0), pre(0) {}
Key key; int v1, nxt, pre;
};
public:
PairingHeaps(int H, int N) : hp(H), nd(N) {
// It consists of `H` Pairing heaps.
// Each heap-node ID can appear at most 1 time(s) among heaps
// and should be in [1, N).
}
void clear(int h) { if (hp[h]) clear_rec(hp[h]), hp[h] = 0; }
void clear_all() {
for (size_t i = 0; i < hp.size(); ++i) hp[i] = 0;
for (size_t i = 0; i < nd.size(); ++i) nd[i] = Node();
}
bool empty(int h) const { return !hp[h]; }
bool used(int v) const { return nd[v].pre >= 0; }
Key min(int h) const { return nd[hp[h]].key; }
int argmin(int h) const { return hp[h]; }
void pop(int h) {
// assert(!empty(h));
erase(h, hp[h]);
}
void push(int h, int v, Key key) {
// assert(!used(v));
nd[v] = Node(key);
hp[h] = merge(hp[h], v);
}
void erase(int h, int v) {
if (!used(v)) return;
int w = two_pass_pairing(nd[v].v1);
if (!nd[v].pre) hp[h] = w;
else {
cut(v);
hp[h] = merge(hp[h], w);
}
nd[v].pre = -1;
}
void decrease_key(int h, int v, Key key) {
if (!used(v)) return push(h, v, key);
if (!nd[v].pre) nd[v].key = key;
else {
cut(v); nd[v].key = key;
hp[h] = merge(hp[h], v);
}
}
private:
void clear_rec(int v) {
for (; v; v = nd[v].nxt) {
if (nd[v].v1) clear_rec(nd[v].v1);
nd[v].pre = -1;
}
}
inline void cut(int v) {
auto& n = nd[v]; int pv = n.pre, nv = n.nxt;
auto& pn = nd[pv];
if (pn.v1 == v) pn.v1 = nv;
else pn.nxt = nv;
nd[nv].pre = pv;
n.nxt = n.pre = 0;
}
int merge(int l, int r) {
if (!l) return r;
if (!r) return l;
if (nd[l].key > nd[r].key) swap(l, r);
int lc = nd[r].nxt = nd[l].v1;
nd[l].v1 = nd[lc].pre = r;
return nd[r].pre = l;
}
int two_pass_pairing(int rt) {
if (!rt) return 0;
int a = rt; rt = 0;
while (a) {
int b = nd[a].nxt, na = 0;
nd[a].pre = nd[a].nxt = 0;
if (b) na = nd[b].nxt, nd[b].pre = nd[b].nxt = 0;
a = merge(a, b);
nd[a].nxt = rt; rt = a; a = na;
}
int s = nd[rt].nxt; nd[rt].nxt = 0;
while (s) {
int t = nd[s].nxt; nd[s].nxt = 0;
rt = merge(rt, s);
s = t;
}
return rt;
}
private:
vector<int> hp;
vector<Node> nd;
};
template <typename T>
struct PriorityQueue : public priority_queue< T, vector<T>, greater<T>> {
PriorityQueue() {}
PriorityQueue(int N) { this->c.reserve(N);}
T min() const { return this->top(); }
void clear() { this->c.clear(); }
};
template <typename T>
struct Queue {
Queue() {}
Queue(int N) : qh(0), qt(0), dat(N) {}
T operator [] (int i) const { return dat[i]; }
void enqueue(int u) { dat[qt++] = u; }
int dequeue() { return dat[qh++]; }
bool empty() const { return qh == qt; }
void clear() { qh = qt = 0; }
int size() const { return qt; }
int qh, qt;
vector<T> dat;
};
public:
struct InputEdge { int u1, to; cost_t cst; };
private:
template <typename T> using ModifiableHeap = BinaryHeap<T>;
template <typename T> using ModifiableHeaps = PairingHeaps<T>;
template <typename T> using FastHeap = PriorityQueue<T>;
struct Edge { int to; cost_t cst; };
struct Link { int u1, to; };
struct Node {
struct NodeLink { int b, v; };
Node() {}
Node(int u) : fa(0), siz(1) { lnk[0] = lnk[1] = {u, u}; }
int next_v() const { return lnk[0].v; }
int next_b() const { return lnk[0].b; }
int prev_v() const { return lnk[1].v; }
int prev_b() const { return lnk[1].b; }
int fa, siz;
NodeLink lnk[2];
};
struct Event {
Event() {}
Event(cost_t tm, int id) : tm(tm), id(id) {}
#ifdef FAST_COMPARE
bool operator < (const Event& rhs) const { return tm < rhs.tm; }
#else
bool operator < (const Event& rhs) const {
if (tm < rhs.tm) { return true; }
if (rhs.tm < tm) { return false; }
return id < rhs.id;
}
#endif
bool operator > (const Event& rhs) const { return rhs.operator<(*this); }
cost_t tm; int id;
};
struct EdgeEvent {
EdgeEvent() {}
EdgeEvent(cost_t tm, int u1, int to) : tm(tm), u1(u1), to(to) {}
#ifdef FAST_COMPARE
bool operator < (const EdgeEvent& rhs) const { return tm < rhs.tm; }
#else
bool operator < (const EdgeEvent& rhs) const {
if (tm < rhs.tm) {
return true;
}
if (tm > rhs.tm) {
return false;
}
return make_pair(u1, to) < make_pair(rhs.u1, rhs.to);
}
#endif
bool operator > (const EdgeEvent& rhs) const { return rhs.operator<(*this); }
cost_t tm; int u1, to;
};
public:
gabow_maximum_weighted_matching(int N, const vector<InputEdge>& in)
: N(N), B((N - 1) / 2), S(N + B + 1), ofs(N + 2), es(in.size() * 2),
h2(S), hps(S, S), h3(es.size()), h4(S) {
for (auto& e : in) ofs[e.u1 + 1]++, ofs[e.to + 1]++;
for (int i = 1; i <= N + 1; ++i) ofs[i] += ofs[i - 1];
for (auto& e : in) {
es[ofs[e.u1]++] = {e.to, e.cst * 2};
es[ofs[e.to]++] = {e.u1, e.cst * 2};
}
for (int i = N + 1; i > 0; --i) ofs[i] = ofs[i - 1];
ofs[0] = 0;
}
tcost_t maximum_weighted_matching(vector<pair<int, int>> &mat, bool ini=false) {
initialize();
set_potential();
if (ini) find_maximal_matching();
for (int u = 1; u <= N; ++u) if (!mt[u]) do_edmonds_search(u);
tcost_t ret = compute_optimal_value();
mat.clear();
for (int u = 1; u <= N; ++u) if (mt[u] > u) {
mat.push_back({u, mt[u]});
}
return ret;
}
private:
tcost_t compute_optimal_value() const {
tcost_t ret = 0;
for (int u = 1; u <= N; ++u) if (mt[u] > u) {
cost_t mx = 0;
for (int eid = ofs[u]; eid < ofs[u + 1]; ++eid) {
if (es[eid].to == mt[u]) mx = max(mx, es[eid].cst);
}
ret += mx;
}
return ret >> 1;
}
inline tcost_t reduced_cost(int u, int v, const Edge& e) const {
return tcost_t(pot[u]) + pot[v] - e.cst;
}
void rematch(int v, int w) {
int t = mt[v]; mt[v] = w;
if (mt[t] != v) return;
if (lnk[v].to == sf1[lnk[v].to]) {
mt[t] = lnk[v].u1;
rematch(mt[t], t);
} else {
int x = lnk[v].u1, y = lnk[v].to;
rematch(x, y); rematch(y, x);
}
}
void fix_mate_and_base(int b) {
if (b <= N) return;
int bv = bas[b], mv = nd[bv].lnk[0].v, bmv = nd[bv].lnk[0].b;
int d = (nd[bmv].lnk[1].v == mt[mv]) ? 0 : 1;
while (1) {
int mv = nd[bv].lnk[d].v, bmv = nd[bv].lnk[d].b;
if (nd[bmv].lnk[1 ^ d].v != mt[mv]) break;
fix_mate_and_base(bv); fix_mate_and_base(bmv);
bv = nd[bmv].lnk[d].b;
}
fix_mate_and_base(bas[b] = bv);
mt[b] = mt[bv];
}
void reset_time() {
tc1 = 0; ev1 = {Inf, 0};
}
void reset_blossom(int b) {
lab[b] = kFree; lnk[b].u1 = 0; sl[b] = Inf; lz[b] = 0;
}
void reset_all() {
lab[0] = kFree; lnk[0].u1 = 0;
for (int v = 1; v <= N; ++v) { // should be optimized for sparse graphs.
if (lab[v] == kOuter) pot[v] -= tc1;
else {
int bv = sf1[v];
pot[v] += lz[bv];
if (lab[bv] == kInner) pot[v] += tc1 - tc[bv];
}
reset_blossom(v);
}
for (int b = N + 1, r = B - ubi; r > 0 && b < S; ++b) if (bas[b] != b) {
if (sf1[b] == b) {
fix_mate_and_base(b);
if (lab[b] == kOuter) pot[b] += (tc1 - tc[b]) << 1;
else if (lab[b] == kInner) fix_blossom_potential<kInner>(b);
else fix_blossom_potential<kFree>(b);
}
hps.clear(b);
reset_blossom(b); --r;
}
que.clear();
reset_time(); h2.clear();
h3.clear(); h4.clear();
}
void do_edmonds_search(int rt) {
if (pot[rt] == 0) return;
link_blossom(sf1[rt], {0, 0});
push_outer_and_fix_potentials(sf1[rt], 0);
for (bool aug = false; !aug; ) {
aug = augment(rt);
if (aug) break;
aug = adjust_dual_variables(rt);
}
reset_all();
}
template <Label Lab>
inline cost_t fix_blossom_potential(int b) {
// Return the amount.
// (If v is an atom, the pot[v] will not be changed.)
cost_t d = lz[b]; lz[b] = 0;
if (Lab == kInner) {
cost_t dt = tc1 - tc[b];
if (b > N) pot[b] -= dt << 1;
d += dt;
}
return d;
}
template <Label Lab>
inline void update_heap2(int x, int y, int by, cost_t t) {
if (t >= sl[y]) return;
sl[y] = t; bf[y] = x;
if (y == by) {
if (Lab != kInner) h2.decrease_key(y, EdgeEvent(t + lz[y], x, y));
} else {
int gy = grp[y];
if (gy != y) {
if (t >= sl[gy]) return;
sl[gy] = t;
}
hps.decrease_key(by, gy, EdgeEvent(t, x, y));
if (Lab == kInner) return;
EdgeEvent m = hps.min(by);
h2.decrease_key(by, EdgeEvent(m.tm + lz[by], m.u1, m.to));
}
}
void activate_heap2_node(int b) {
if (b <= N) {
if (sl[b] < Inf) h2.push(b, EdgeEvent(sl[b] + lz[b], bf[b], b));
} else {
if (hps.empty(b)) return;
EdgeEvent m = hps.min(b);
h2.push(b, EdgeEvent(m.tm + lz[b], m.u1, m.to));
}
}
void swap_blossom(int a, int b) {
// Assume that `b` is a maximal blossom.
swap(bas[a], bas[b]); if (bas[a] == a) bas[a] = b;
swap(son[a], son[b]); if (son[a] == a) son[a] = b;
swap(lnk[a], lnk[b]);
swap(mt[a], mt[b]);
swap(pot[a], pot[b]); swap(lz[a], lz[b]);
swap(tc[a], tc[b]);
for (int d = 0; d < 2; ++d) nd[nd[a].lnk[d].b].lnk[1 ^ d].b = b;
swap(nd[a], nd[b]);
}
void set_surface_and_group(int b, int sf, int g) {
sf1[b] = sf, grp[b] = g;
if (b <= N) return;
for (int bb = bas[b]; sf1[bb] != sf; bb = nd[bb].next_b()) {
set_surface_and_group(bb, sf, g);
}
}
void merge_smaller_blossoms(int bid) {
int lb = bid, ms = 1;
for (int bta = bas[bid], b = bta; ;) {
if (nd[b].siz > ms) ms = nd[b].siz, lb = b;
if ((b = nd[b].next_b()) == bta) break;
}
for (int bta = bas[bid], b = bta; ;) {
if (b != lb) set_surface_and_group(b, lb, b);
if ((b = nd[b].next_b()) == bta) break;
}
grp[lb] = lb;
if (ms > 1) {
sf1[bid] = son[bid] = lb;
swap_blossom(lb, bid);
} else son[bid] = 0;
}
void contract(int x, int y, int eid) {
int bx = sf1[x], by = sf1[y]; assert(bx != by);
const int h = -(eid + 1);
lnk[sf1[mt[bx]]].u1 = lnk[sf1[mt[by]]].u1 = h;
int lca = -1;
while (1) {
if (mt[by] != 0) swap(bx, by);
bx = lca = sf1[lnk[bx].u1];
if (lnk[sf1[mt[bx]]].u1 == h) break;
lnk[sf1[mt[bx]]].u1 = h;
}
const int bid = ub[--ubi]; assert(ubi >= 0);
int tsz = 0;
for (int d = 0; d < 2; ++d) {
for (int bv = sf1[x]; bv != lca; ) {
int mv = mt[bv], bmv = sf1[mv], v = mt[mv];
int f = lnk[v].u1, t = lnk[v].to;
tsz += nd[bv].siz + nd[bmv].siz;
lnk[mv] = {x, y};
if (bv > N) pot[bv] += (tc1 - tc[bv]) << 1;
if (bmv > N) h4.erase(bmv);
push_outer_and_fix_potentials(bmv, fix_blossom_potential<kInner>(bmv));
nd[bv].lnk[d] = {bmv, mv};
nd[bmv].lnk[1 ^ d] = {bv, v}; nd[bmv].lnk[d] = {bv = sf1[f], f};
nd[bv].lnk[1 ^ d] = {bmv, t};
}
nd[sf1[x]].lnk[1 ^ d] = {sf1[y], y};
swap(x, y);
}
if (lca > N) pot[lca] += (tc1 - tc[lca]) << 1;
nd[bid].siz = tsz + nd[lca].siz;
bas[bid] = lca; lnk[bid] = lnk[lca]; mt[bid] = mt[lca];
lab[bid] = kOuter;
sf1[bid] = bid; tc[bid] = tc1;
pot[bid] = 0; lz[bid] = 0;
merge_smaller_blossoms(bid); // O(n log n) time / Edmonds search
}
void link_blossom(int v, Link l) {
lnk[v] = {l.u1, l.to};
if (v <= N) return;
int b = bas[v]; link_blossom(b, l);
int pb = nd[b].prev_b();
l = {nd[pb].next_v(), nd[b].prev_v()};
for (int bv = b; ; ) {
int bw = nd[bv].next_b();
if (bw == b) break;
link_blossom(bw, l);
Link nl = {nd[bw].prev_v(), nd[bv].next_v()};
bv = nd[bw].next_b();
link_blossom(bv, nl);
}
}
void push_outer_and_fix_potentials(int v, cost_t d) {
lab[v] = kOuter;
if (v > N) {
for (int b = bas[v]; lab[b] != kOuter; b = nd[b].next_b()) {
push_outer_and_fix_potentials(b, d);
}
} else {
pot[v] += tc1 + d;
if (pot[v] < ev1.tm) ev1 = {pot[v], v};
que.enqueue(v);
}
}
bool grow(int x, int y) {
int by = sf1[y];
bool vis = (lab[by] != kFree);
if (!vis) link_blossom(by, {0, 0});
lab[by] = kInner; tc[by] = tc1; h2.erase(by);
if (y != by) h4.update(by, tc1 + (pot[by] >> 1));
int z = mt[by];
if (z == 0) {
rematch(x, y); rematch(y, x);
return true;
}
int bz = sf1[z];
if (!vis) link_blossom(bz, {x, y});
else lnk[bz] = lnk[z] = {x, y};
push_outer_and_fix_potentials(bz, fix_blossom_potential<kFree>(bz));
tc[bz] = tc1; h2.erase(bz);
return false;
}
void free_blossom(int bid) {
ub[ubi++] = bid;
bas[bid] = bid;
}
int recalculate_minimum_slack(int b, int g) {
// Return the destination of the best edge of blossom `g`.
if (b <= N) {
if (sl[b] >= sl[g]) return 0;
sl[g] = sl[b]; bf[g] = bf[b];
return b;
}
int v = 0;
for (int bta = bas[b], bb = bta; ; ) {
int w = recalculate_minimum_slack(bb, g);
if (w != 0) v = w;
if ((bb = nd[bb].next_b()) == bta) break;
}
return v;
}
void construct_smaller_components(int b, int sf, int g) {
sf1[b] = sf, grp[b] = g; // `grp[b] = g` is unneeded.
if (b <= N) return;
for (int bb = bas[b]; sf1[bb] != sf; bb = nd[bb].next_b()) {
if (bb == son[b]) {
construct_smaller_components(bb, sf, g);
} else {
set_surface_and_group(bb, sf, bb);
int to = 0;
if (bb > N) sl[bb] = Inf, to = recalculate_minimum_slack(bb, bb);
else if (sl[bb] < Inf) to = bb;
if (to > 0) hps.push(sf, bb, EdgeEvent(sl[bb], bf[bb], to));
}
}
}
void move_to_largest_blossom(int bid) {
const int h = son[bid];
cost_t d = (tc1 - tc[bid]) + lz[bid]; lz[bid] = 0;
for (int bta = bas[bid], b = bta; ;) {
tc[b] = tc1;
lz[b] = d;
if (b != h) construct_smaller_components(b, b, b), hps.erase(bid, b);
if ((b = nd[b].next_b()) == bta) break;
}
if (h > 0) swap_blossom(h, bid), bid = h;
free_blossom(bid);
}
void expand(int bid) {
int mv = mt[bas[bid]];
move_to_largest_blossom(bid); // O(n log n) time / Edmonds search
Link ol = lnk[mv];
int ob = sf1[mt[mv]], rt = sf1[ol.to];
int d = (mt[rt] == nd[rt].lnk[0].v) ? 1 : 0;
for (int b = nd[ob].lnk[d ^ 1].b; b != rt; ) {
lab[b] = kSeparated; activate_heap2_node(b); b = nd[b].lnk[d ^ 1].b;
lab[b] = kSeparated; activate_heap2_node(b); b = nd[b].lnk[d ^ 1].b;
}
for (int b = ob; ; b = nd[b].lnk[d].b) {
lab[b] = kInner;
int nb = nd[b].lnk[d].b;
if (b == rt) lnk[mt[b]] = ol;
else lnk[mt[b]] = {nd[b].lnk[d].v, nd[nb].lnk[d ^ 1].v};
lnk[sf1[mt[b]]] = lnk[mt[b]]; // fix tree links
if (b > N) {
if (pot[b] == 0) expand(b);
else h4.push(b, tc1 + (pot[b] >> 1));
}
if (b == rt) break;
push_outer_and_fix_potentials(nb, fix_blossom_potential<kInner>(b = nb));
}
}
bool augment(int rt) {
// Return true if an augmenting path is found.
while (!que.empty()) {
int x = que.dequeue(), bx = sf1[x];
if (pot[x] == tc1) {
if (x != rt) rematch(x, 0);
return true;
}
for (int eid = ofs[x]; eid < ofs[x + 1]; ++eid) {
auto& e = es[eid]; int y = e.to, by = sf1[y];
if (bx == by) continue;
Label l = lab[by];
if (l == kOuter) {
cost_t t = reduced_cost(x, y, e) >> 1; // < 2 * Inf
if (t == tc1) {
contract(x, y, eid); bx = sf1[x];
} else if (t < ev1.tm) {
h3.emplace(t, x, eid);
}
} else {
tcost_t t = reduced_cost(x, y, e); // < 3 * Inf
if (t >= Inf) continue;
if (l != kInner) {
if (cost_t(t) + lz[by] == tc1) {
if (grow(x, y)) return true;
} else update_heap2<kFree>(x, y, by, t);
} else {
if (mt[x] != y) update_heap2<kInner>(x, y, by, t);
}
}
}
}
return false;
}
bool adjust_dual_variables(int rt) {
// delta1 : rematch
cost_t t1 = ev1.tm;
// delta2 : grow
cost_t t2 = Inf;
if (!h2.empty()) t2 = h2.min().tm;
// delta3 : contract : O(m log n) time / Edmonds search [ bottleneck (?) ]
cost_t t3 = Inf;
while (!h3.empty()) {
EdgeEvent e = h3.min();
int x = e.u1, y = es[e.to].to; // e.to is some edge id.
if (sf1[x] != sf1[y]) {
t3 = e.tm;
break;
} else h3.pop();
}
// delta4 : expand
cost_t t4 = Inf;
if (!h4.empty()) t4 = h4.min();
// -- events --
cost_t tn = min(min(t1, t2), min(t3, t4));
assert(tc1 <= tn && tn < Inf);
tc1 = tn;
if (tc1 == ev1.tm) {
int x = ev1.id;
if (x != rt) rematch(x, 0);
return true;
}
while (!h2.empty() && h2.min().tm == tc1) {
int x = h2.min().u1, y = h2.min().to;
if (grow(x, y)) return true; // `grow` function will call `h2.erase(by)`.
}
while (!h3.empty() && h3.min().tm == tc1) {
int x = h3.min().u1, eid = h3.min().to;
int y = es[eid].to; h3.pop();
if (sf1[x] == sf1[y]) continue;
contract(x, y, eid);
}
while (!h4.empty() && h4.min() == tc1) {
int b = h4.argmin(); h4.pop();
expand(b);
}
return false;
}
private:
void initialize() {
que = Queue<int>(N);
mt.assign(S, 0);
lnk.assign(S, {0, 0});
lab.assign(S, kFree);
bas.resize(S); for (int u = 1; u < S; ++u) bas[u] = u;
sf1.resize(S); for (int u = 1; u < S; ++u) sf1[u] = u;
pot.resize(S);
nd.resize(S); for (int b = 1; b < S; ++b) nd[b] = Node(b);
ub.resize(B); for (int i = 0; i < B; ++i) ub[i] = N + B - i;
ubi = B;
// for O(nm log n) implementation
reset_time();
tc.resize(S);
sl.resize(S); for (int i = 0; i < S; ++i) sl[i] = Inf;
bf.assign(S, 0);
son.assign(S, 0);
lz.assign(S, 0);
grp.resize(S); for (int i = 0; i < S; ++i) grp[i] = i;
}
void set_potential() {
for (int u = 1; u <= N; ++u) {
cost_t mx = 0;
for (int eid = ofs[u]; eid < ofs[u + 1]; ++eid) {
mx = max(mx, es[eid].cst);
}
pot[u] = mx >> 1;
}
}
void find_maximal_matching() {
// Find a maximal matching naively.
for (int u = 1; u <= N; ++u) if (!mt[u]) {
for (int eid = ofs[u]; eid < ofs[u + 1]; ++eid) {
auto& e = es[eid]; int v = e.to;
if (mt[v] > 0 || reduced_cost(u, v, e) > 0) continue;
mt[u] = v; mt[v] = u;
break;
}
}
}
private:
const int N, B, S; // N = |V|, B = (|V| - 1) / 2, S = N + B + 1
vector<int> ofs;
vector<Edge> es;
Queue<int> que;
vector<int> mt, sf1, bas;
vector<Link> lnk;
vector<Label> lab;
vector<cost_t> pot;
vector<int> ub; int ubi;
vector<Node> nd;
// for O(nm log n) implementation
vector<int> son, grp;
vector<cost_t> tc, lz, sl;
vector<int> bf;
cost_t tc1;
Event ev1;
ModifiableHeap<EdgeEvent> h2;
ModifiableHeaps<EdgeEvent> hps;
FastHeap<EdgeEvent> h3;
ModifiableHeap<cost_t> h4;
};
} // namespace general_weighted_matching_internal
template <class TotalCost> struct weighted_matching_result {
TotalCost wt{};
std::vector<std::pair<int, int>> es;
};
/// @brief Compute a maximum-weight matching in a general undirected graph.
/// A primal-dual Edmonds search grows an alternating forest through zero
/// reduced-cost edges. Odd cycles are contracted into blossoms; dual changes
/// schedule grow, contract, and expand events. Gabow's heavy-blossom grouping
/// and meldable heaps make all slack/event maintenance near-linear per search.
template <class Cost, class TotalCost = Cost>
weighted_matching_result<TotalCost> general_weighted_matching(
int n1,
const std::vector<std::tuple<int, int, Cost>> &es1) {
using solver = general_weighted_matching_internal::
gabow_maximum_weighted_matching<Cost, TotalCost>;
std::vector<typename solver::InputEdge> es;
es.reserve(es1.size());
for (auto [a1, b1, wt] : es1) {
assert(0 <= a1 && a1 < n1);
assert(0 <= b1 && b1 < n1);
assert(a1 != b1 && wt >= Cost{});
es.push_back({a1 + 1, b1 + 1, wt});
}
solver mat(n1, es);
weighted_matching_result<TotalCost> res;
res.wt = mat.maximum_weighted_matching(res.es);
for (auto &[a1, b1] : res.es) {
--a1;
--b1;
}
return res;
}
} // namespace noya
#endif // NOYA_GENERAL_WEIGHTED_MATCHING_HPP
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <functional>
#include <queue>
#include <tuple>
#include <utility>
#include <vector>
/// @complexity Time: O(n m log n). Space: O(n+m).
namespace noya {
namespace general_weighted_matching_internal {
using namespace std;
template <typename CostType, typename TotalCostType>
class gabow_maximum_weighted_matching {
public:
using cost_t = CostType;
using tcost_t = TotalCostType;
private:
enum Label { kSeparated = -2, kInner = -1, kFree = 0, kOuter = 1 };
static constexpr cost_t Inf = cost_t(1) << (sizeof(cost_t) * 8 - 2);
private:
template <typename T>
class BinaryHeap {
public:
struct Node {
#ifdef FAST_COMPARE
bool operator < (const Node& rhs) const { return val < rhs.val; }
#else
bool operator < (const Node& rhs) const {
if (val < rhs.val) { return true; }
if (rhs.val < val) { return false; }
return id < rhs.id;
}
#endif
T val; int id;
};
BinaryHeap() {}
BinaryHeap(int N) : sz_(0), nd(N + 1), idx(N, 0) {}
int size() const { return sz_; }
bool empty() const { return sz_ == 0; }
void clear() { while (sz_ > 0) idx[nd[sz_--].id] = 0; }
T min() const { return nd[1].val; }
int argmin() const { return nd[1].id; } // argmin ?
T get_val(int id) const { return nd[idx[id]].val; }
void pop() { if (sz_ > 0) pop(1); }
void erase(int id) { if (idx[id]) pop(idx[id]); }
bool has(int id) const { return idx[id] != 0; }
void update(int id, T v) {
if (!has(id)) return push(id, v);
bool up = (v < nd[idx[id]].val);
nd[idx[id]].val = v;
if (up) up_heap(idx[id]);
else down_heap(idx[id]);
}
void decrease_key(int id, T v) {
if (!has(id)) return push(id, v);
if (v < nd[idx[id]].val) nd[idx[id]].val = v, up_heap(idx[id]);
}
void push(int id, T v) {
// assert(!has(id));
idx[id] = ++sz_; nd[sz_] = {v, id};
up_heap(sz_);
}
private:
void pop(int pos) {
idx[nd[pos].id] = 0;
if (pos == sz_) { --sz_; return; }
bool up = (nd[sz_].val < nd[pos].val);
nd[pos] = nd[sz_--]; idx[nd[pos].id] = pos;
if (up) up_heap(pos);
else down_heap(pos);
}
void swap_node(int a, int b) {
swap(nd[a], nd[b]); idx[nd[a].id] = a; idx[nd[b].id] = b;
}
void down_heap(int pos) {
for (int k = pos, nk = k; 2 * k <= sz_; k = nk) {
if (nd[2 * k] < nd[nk]) nk = 2 * k;
if (2 * k + 1 <= sz_ && nd[2 * k + 1] < nd[nk]) nk = 2 * k + 1;
if (nk == k) break;
swap_node(k, nk);
}
}
void up_heap(int pos) {
for (int k = pos; k > 1 && nd[k] < nd[k >> 1]; k >>= 1) swap_node(k, k >> 1);
}
int sz_;
vector<Node> nd;
vector<int> idx;
};
template <typename Key>
class PairingHeaps {
private:
struct Node {
Node() : pre(-1) {} // "pre < 0" means the node is unused.
Node(Key v) : key(v), v1(0), nxt(0), pre(0) {}
Key key; int v1, nxt, pre;
};
public:
PairingHeaps(int H, int N) : hp(H), nd(N) {
// It consists of `H` Pairing heaps.
// Each heap-node ID can appear at most 1 time(s) among heaps
// and should be in [1, N).
}
void clear(int h) { if (hp[h]) clear_rec(hp[h]), hp[h] = 0; }
void clear_all() {
for (size_t i = 0; i < hp.size(); ++i) hp[i] = 0;
for (size_t i = 0; i < nd.size(); ++i) nd[i] = Node();
}
bool empty(int h) const { return !hp[h]; }
bool used(int v) const { return nd[v].pre >= 0; }
Key min(int h) const { return nd[hp[h]].key; }
int argmin(int h) const { return hp[h]; }
void pop(int h) {
// assert(!empty(h));
erase(h, hp[h]);
}
void push(int h, int v, Key key) {
// assert(!used(v));
nd[v] = Node(key);
hp[h] = merge(hp[h], v);
}
void erase(int h, int v) {
if (!used(v)) return;
int w = two_pass_pairing(nd[v].v1);
if (!nd[v].pre) hp[h] = w;
else {
cut(v);
hp[h] = merge(hp[h], w);
}
nd[v].pre = -1;
}
void decrease_key(int h, int v, Key key) {
if (!used(v)) return push(h, v, key);
if (!nd[v].pre) nd[v].key = key;
else {
cut(v); nd[v].key = key;
hp[h] = merge(hp[h], v);
}
}
private:
void clear_rec(int v) {
for (; v; v = nd[v].nxt) {
if (nd[v].v1) clear_rec(nd[v].v1);
nd[v].pre = -1;
}
}
inline void cut(int v) {
auto& n = nd[v]; int pv = n.pre, nv = n.nxt;
auto& pn = nd[pv];
if (pn.v1 == v) pn.v1 = nv;
else pn.nxt = nv;
nd[nv].pre = pv;
n.nxt = n.pre = 0;
}
int merge(int l, int r) {
if (!l) return r;
if (!r) return l;
if (nd[l].key > nd[r].key) swap(l, r);
int lc = nd[r].nxt = nd[l].v1;
nd[l].v1 = nd[lc].pre = r;
return nd[r].pre = l;
}
int two_pass_pairing(int rt) {
if (!rt) return 0;
int a = rt; rt = 0;
while (a) {
int b = nd[a].nxt, na = 0;
nd[a].pre = nd[a].nxt = 0;
if (b) na = nd[b].nxt, nd[b].pre = nd[b].nxt = 0;
a = merge(a, b);
nd[a].nxt = rt; rt = a; a = na;
}
int s = nd[rt].nxt; nd[rt].nxt = 0;
while (s) {
int t = nd[s].nxt; nd[s].nxt = 0;
rt = merge(rt, s);
s = t;
}
return rt;
}
private:
vector<int> hp;
vector<Node> nd;
};
template <typename T>
struct PriorityQueue : public priority_queue< T, vector<T>, greater<T>> {
PriorityQueue() {}
PriorityQueue(int N) { this->c.reserve(N);}
T min() const { return this->top(); }
void clear() { this->c.clear(); }
};
template <typename T>
struct Queue {
Queue() {}
Queue(int N) : qh(0), qt(0), dat(N) {}
T operator [] (int i) const { return dat[i]; }
void enqueue(int u) { dat[qt++] = u; }
int dequeue() { return dat[qh++]; }
bool empty() const { return qh == qt; }
void clear() { qh = qt = 0; }
int size() const { return qt; }
int qh, qt;
vector<T> dat;
};
public:
struct InputEdge { int u1, to; cost_t cst; };
private:
template <typename T> using ModifiableHeap = BinaryHeap<T>;
template <typename T> using ModifiableHeaps = PairingHeaps<T>;
template <typename T> using FastHeap = PriorityQueue<T>;
struct Edge { int to; cost_t cst; };
struct Link { int u1, to; };
struct Node {
struct NodeLink { int b, v; };
Node() {}
Node(int u) : fa(0), siz(1) { lnk[0] = lnk[1] = {u, u}; }
int next_v() const { return lnk[0].v; }
int next_b() const { return lnk[0].b; }
int prev_v() const { return lnk[1].v; }
int prev_b() const { return lnk[1].b; }
int fa, siz;
NodeLink lnk[2];
};
struct Event {
Event() {}
Event(cost_t tm, int id) : tm(tm), id(id) {}
#ifdef FAST_COMPARE
bool operator < (const Event& rhs) const { return tm < rhs.tm; }
#else
bool operator < (const Event& rhs) const {
if (tm < rhs.tm) { return true; }
if (rhs.tm < tm) { return false; }
return id < rhs.id;
}
#endif
bool operator > (const Event& rhs) const { return rhs.operator<(*this); }
cost_t tm; int id;
};
struct EdgeEvent {
EdgeEvent() {}
EdgeEvent(cost_t tm, int u1, int to) : tm(tm), u1(u1), to(to) {}
#ifdef FAST_COMPARE
bool operator < (const EdgeEvent& rhs) const { return tm < rhs.tm; }
#else
bool operator < (const EdgeEvent& rhs) const {
if (tm < rhs.tm) {
return true;
}
if (tm > rhs.tm) {
return false;
}
return make_pair(u1, to) < make_pair(rhs.u1, rhs.to);
}
#endif
bool operator > (const EdgeEvent& rhs) const { return rhs.operator<(*this); }
cost_t tm; int u1, to;
};
public:
gabow_maximum_weighted_matching(int N, const vector<InputEdge>& in)
: N(N), B((N - 1) / 2), S(N + B + 1), ofs(N + 2), es(in.size() * 2),
h2(S), hps(S, S), h3(es.size()), h4(S) {
for (auto& e : in) ofs[e.u1 + 1]++, ofs[e.to + 1]++;
for (int i = 1; i <= N + 1; ++i) ofs[i] += ofs[i - 1];
for (auto& e : in) {
es[ofs[e.u1]++] = {e.to, e.cst * 2};
es[ofs[e.to]++] = {e.u1, e.cst * 2};
}
for (int i = N + 1; i > 0; --i) ofs[i] = ofs[i - 1];
ofs[0] = 0;
}
tcost_t maximum_weighted_matching(vector<pair<int, int>> &mat, bool ini=false) {
initialize();
set_potential();
if (ini) find_maximal_matching();
for (int u = 1; u <= N; ++u) if (!mt[u]) do_edmonds_search(u);
tcost_t ret = compute_optimal_value();
mat.clear();
for (int u = 1; u <= N; ++u) if (mt[u] > u) {
mat.push_back({u, mt[u]});
}
return ret;
}
private:
tcost_t compute_optimal_value() const {
tcost_t ret = 0;
for (int u = 1; u <= N; ++u) if (mt[u] > u) {
cost_t mx = 0;
for (int eid = ofs[u]; eid < ofs[u + 1]; ++eid) {
if (es[eid].to == mt[u]) mx = max(mx, es[eid].cst);
}
ret += mx;
}
return ret >> 1;
}
inline tcost_t reduced_cost(int u, int v, const Edge& e) const {
return tcost_t(pot[u]) + pot[v] - e.cst;
}
void rematch(int v, int w) {
int t = mt[v]; mt[v] = w;
if (mt[t] != v) return;
if (lnk[v].to == sf1[lnk[v].to]) {
mt[t] = lnk[v].u1;
rematch(mt[t], t);
} else {
int x = lnk[v].u1, y = lnk[v].to;
rematch(x, y); rematch(y, x);
}
}
void fix_mate_and_base(int b) {
if (b <= N) return;
int bv = bas[b], mv = nd[bv].lnk[0].v, bmv = nd[bv].lnk[0].b;
int d = (nd[bmv].lnk[1].v == mt[mv]) ? 0 : 1;
while (1) {
int mv = nd[bv].lnk[d].v, bmv = nd[bv].lnk[d].b;
if (nd[bmv].lnk[1 ^ d].v != mt[mv]) break;
fix_mate_and_base(bv); fix_mate_and_base(bmv);
bv = nd[bmv].lnk[d].b;
}
fix_mate_and_base(bas[b] = bv);
mt[b] = mt[bv];
}
void reset_time() {
tc1 = 0; ev1 = {Inf, 0};
}
void reset_blossom(int b) {
lab[b] = kFree; lnk[b].u1 = 0; sl[b] = Inf; lz[b] = 0;
}
void reset_all() {
lab[0] = kFree; lnk[0].u1 = 0;
for (int v = 1; v <= N; ++v) { // should be optimized for sparse graphs.
if (lab[v] == kOuter) pot[v] -= tc1;
else {
int bv = sf1[v];
pot[v] += lz[bv];
if (lab[bv] == kInner) pot[v] += tc1 - tc[bv];
}
reset_blossom(v);
}
for (int b = N + 1, r = B - ubi; r > 0 && b < S; ++b) if (bas[b] != b) {
if (sf1[b] == b) {
fix_mate_and_base(b);
if (lab[b] == kOuter) pot[b] += (tc1 - tc[b]) << 1;
else if (lab[b] == kInner) fix_blossom_potential<kInner>(b);
else fix_blossom_potential<kFree>(b);
}
hps.clear(b);
reset_blossom(b); --r;
}
que.clear();
reset_time(); h2.clear();
h3.clear(); h4.clear();
}
void do_edmonds_search(int rt) {
if (pot[rt] == 0) return;
link_blossom(sf1[rt], {0, 0});
push_outer_and_fix_potentials(sf1[rt], 0);
for (bool aug = false; !aug; ) {
aug = augment(rt);
if (aug) break;
aug = adjust_dual_variables(rt);
}
reset_all();
}
template <Label Lab>
inline cost_t fix_blossom_potential(int b) {
// Return the amount.
// (If v is an atom, the pot[v] will not be changed.)
cost_t d = lz[b]; lz[b] = 0;
if (Lab == kInner) {
cost_t dt = tc1 - tc[b];
if (b > N) pot[b] -= dt << 1;
d += dt;
}
return d;
}
template <Label Lab>
inline void update_heap2(int x, int y, int by, cost_t t) {
if (t >= sl[y]) return;
sl[y] = t; bf[y] = x;
if (y == by) {
if (Lab != kInner) h2.decrease_key(y, EdgeEvent(t + lz[y], x, y));
} else {
int gy = grp[y];
if (gy != y) {
if (t >= sl[gy]) return;
sl[gy] = t;
}
hps.decrease_key(by, gy, EdgeEvent(t, x, y));
if (Lab == kInner) return;
EdgeEvent m = hps.min(by);
h2.decrease_key(by, EdgeEvent(m.tm + lz[by], m.u1, m.to));
}
}
void activate_heap2_node(int b) {
if (b <= N) {
if (sl[b] < Inf) h2.push(b, EdgeEvent(sl[b] + lz[b], bf[b], b));
} else {
if (hps.empty(b)) return;
EdgeEvent m = hps.min(b);
h2.push(b, EdgeEvent(m.tm + lz[b], m.u1, m.to));
}
}
void swap_blossom(int a, int b) {
// Assume that `b` is a maximal blossom.
swap(bas[a], bas[b]); if (bas[a] == a) bas[a] = b;
swap(son[a], son[b]); if (son[a] == a) son[a] = b;
swap(lnk[a], lnk[b]);
swap(mt[a], mt[b]);
swap(pot[a], pot[b]); swap(lz[a], lz[b]);
swap(tc[a], tc[b]);
for (int d = 0; d < 2; ++d) nd[nd[a].lnk[d].b].lnk[1 ^ d].b = b;
swap(nd[a], nd[b]);
}
void set_surface_and_group(int b, int sf, int g) {
sf1[b] = sf, grp[b] = g;
if (b <= N) return;
for (int bb = bas[b]; sf1[bb] != sf; bb = nd[bb].next_b()) {
set_surface_and_group(bb, sf, g);
}
}
void merge_smaller_blossoms(int bid) {
int lb = bid, ms = 1;
for (int bta = bas[bid], b = bta; ;) {
if (nd[b].siz > ms) ms = nd[b].siz, lb = b;
if ((b = nd[b].next_b()) == bta) break;
}
for (int bta = bas[bid], b = bta; ;) {
if (b != lb) set_surface_and_group(b, lb, b);
if ((b = nd[b].next_b()) == bta) break;
}
grp[lb] = lb;
if (ms > 1) {
sf1[bid] = son[bid] = lb;
swap_blossom(lb, bid);
} else son[bid] = 0;
}
void contract(int x, int y, int eid) {
int bx = sf1[x], by = sf1[y]; assert(bx != by);
const int h = -(eid + 1);
lnk[sf1[mt[bx]]].u1 = lnk[sf1[mt[by]]].u1 = h;
int lca = -1;
while (1) {
if (mt[by] != 0) swap(bx, by);
bx = lca = sf1[lnk[bx].u1];
if (lnk[sf1[mt[bx]]].u1 == h) break;
lnk[sf1[mt[bx]]].u1 = h;
}
const int bid = ub[--ubi]; assert(ubi >= 0);
int tsz = 0;
for (int d = 0; d < 2; ++d) {
for (int bv = sf1[x]; bv != lca; ) {
int mv = mt[bv], bmv = sf1[mv], v = mt[mv];
int f = lnk[v].u1, t = lnk[v].to;
tsz += nd[bv].siz + nd[bmv].siz;
lnk[mv] = {x, y};
if (bv > N) pot[bv] += (tc1 - tc[bv]) << 1;
if (bmv > N) h4.erase(bmv);
push_outer_and_fix_potentials(bmv, fix_blossom_potential<kInner>(bmv));
nd[bv].lnk[d] = {bmv, mv};
nd[bmv].lnk[1 ^ d] = {bv, v}; nd[bmv].lnk[d] = {bv = sf1[f], f};
nd[bv].lnk[1 ^ d] = {bmv, t};
}
nd[sf1[x]].lnk[1 ^ d] = {sf1[y], y};
swap(x, y);
}
if (lca > N) pot[lca] += (tc1 - tc[lca]) << 1;
nd[bid].siz = tsz + nd[lca].siz;
bas[bid] = lca; lnk[bid] = lnk[lca]; mt[bid] = mt[lca];
lab[bid] = kOuter;
sf1[bid] = bid; tc[bid] = tc1;
pot[bid] = 0; lz[bid] = 0;
merge_smaller_blossoms(bid); // O(n log n) time / Edmonds search
}
void link_blossom(int v, Link l) {
lnk[v] = {l.u1, l.to};
if (v <= N) return;
int b = bas[v]; link_blossom(b, l);
int pb = nd[b].prev_b();
l = {nd[pb].next_v(), nd[b].prev_v()};
for (int bv = b; ; ) {
int bw = nd[bv].next_b();
if (bw == b) break;
link_blossom(bw, l);
Link nl = {nd[bw].prev_v(), nd[bv].next_v()};
bv = nd[bw].next_b();
link_blossom(bv, nl);
}
}
void push_outer_and_fix_potentials(int v, cost_t d) {
lab[v] = kOuter;
if (v > N) {
for (int b = bas[v]; lab[b] != kOuter; b = nd[b].next_b()) {
push_outer_and_fix_potentials(b, d);
}
} else {
pot[v] += tc1 + d;
if (pot[v] < ev1.tm) ev1 = {pot[v], v};
que.enqueue(v);
}
}
bool grow(int x, int y) {
int by = sf1[y];
bool vis = (lab[by] != kFree);
if (!vis) link_blossom(by, {0, 0});
lab[by] = kInner; tc[by] = tc1; h2.erase(by);
if (y != by) h4.update(by, tc1 + (pot[by] >> 1));
int z = mt[by];
if (z == 0) {
rematch(x, y); rematch(y, x);
return true;
}
int bz = sf1[z];
if (!vis) link_blossom(bz, {x, y});
else lnk[bz] = lnk[z] = {x, y};
push_outer_and_fix_potentials(bz, fix_blossom_potential<kFree>(bz));
tc[bz] = tc1; h2.erase(bz);
return false;
}
void free_blossom(int bid) {
ub[ubi++] = bid;
bas[bid] = bid;
}
int recalculate_minimum_slack(int b, int g) {
// Return the destination of the best edge of blossom `g`.
if (b <= N) {
if (sl[b] >= sl[g]) return 0;
sl[g] = sl[b]; bf[g] = bf[b];
return b;
}
int v = 0;
for (int bta = bas[b], bb = bta; ; ) {
int w = recalculate_minimum_slack(bb, g);
if (w != 0) v = w;
if ((bb = nd[bb].next_b()) == bta) break;
}
return v;
}
void construct_smaller_components(int b, int sf, int g) {
sf1[b] = sf, grp[b] = g; // `grp[b] = g` is unneeded.
if (b <= N) return;
for (int bb = bas[b]; sf1[bb] != sf; bb = nd[bb].next_b()) {
if (bb == son[b]) {
construct_smaller_components(bb, sf, g);
} else {
set_surface_and_group(bb, sf, bb);
int to = 0;
if (bb > N) sl[bb] = Inf, to = recalculate_minimum_slack(bb, bb);
else if (sl[bb] < Inf) to = bb;
if (to > 0) hps.push(sf, bb, EdgeEvent(sl[bb], bf[bb], to));
}
}
}
void move_to_largest_blossom(int bid) {
const int h = son[bid];
cost_t d = (tc1 - tc[bid]) + lz[bid]; lz[bid] = 0;
for (int bta = bas[bid], b = bta; ;) {
tc[b] = tc1;
lz[b] = d;
if (b != h) construct_smaller_components(b, b, b), hps.erase(bid, b);
if ((b = nd[b].next_b()) == bta) break;
}
if (h > 0) swap_blossom(h, bid), bid = h;
free_blossom(bid);
}
void expand(int bid) {
int mv = mt[bas[bid]];
move_to_largest_blossom(bid); // O(n log n) time / Edmonds search
Link ol = lnk[mv];
int ob = sf1[mt[mv]], rt = sf1[ol.to];
int d = (mt[rt] == nd[rt].lnk[0].v) ? 1 : 0;
for (int b = nd[ob].lnk[d ^ 1].b; b != rt; ) {
lab[b] = kSeparated; activate_heap2_node(b); b = nd[b].lnk[d ^ 1].b;
lab[b] = kSeparated; activate_heap2_node(b); b = nd[b].lnk[d ^ 1].b;
}
for (int b = ob; ; b = nd[b].lnk[d].b) {
lab[b] = kInner;
int nb = nd[b].lnk[d].b;
if (b == rt) lnk[mt[b]] = ol;
else lnk[mt[b]] = {nd[b].lnk[d].v, nd[nb].lnk[d ^ 1].v};
lnk[sf1[mt[b]]] = lnk[mt[b]]; // fix tree links
if (b > N) {
if (pot[b] == 0) expand(b);
else h4.push(b, tc1 + (pot[b] >> 1));
}
if (b == rt) break;
push_outer_and_fix_potentials(nb, fix_blossom_potential<kInner>(b = nb));
}
}
bool augment(int rt) {
// Return true if an augmenting path is found.
while (!que.empty()) {
int x = que.dequeue(), bx = sf1[x];
if (pot[x] == tc1) {
if (x != rt) rematch(x, 0);
return true;
}
for (int eid = ofs[x]; eid < ofs[x + 1]; ++eid) {
auto& e = es[eid]; int y = e.to, by = sf1[y];
if (bx == by) continue;
Label l = lab[by];
if (l == kOuter) {
cost_t t = reduced_cost(x, y, e) >> 1; // < 2 * Inf
if (t == tc1) {
contract(x, y, eid); bx = sf1[x];
} else if (t < ev1.tm) {
h3.emplace(t, x, eid);
}
} else {
tcost_t t = reduced_cost(x, y, e); // < 3 * Inf
if (t >= Inf) continue;
if (l != kInner) {
if (cost_t(t) + lz[by] == tc1) {
if (grow(x, y)) return true;
} else update_heap2<kFree>(x, y, by, t);
} else {
if (mt[x] != y) update_heap2<kInner>(x, y, by, t);
}
}
}
}
return false;
}
bool adjust_dual_variables(int rt) {
// delta1 : rematch
cost_t t1 = ev1.tm;
// delta2 : grow
cost_t t2 = Inf;
if (!h2.empty()) t2 = h2.min().tm;
// delta3 : contract : O(m log n) time / Edmonds search [ bottleneck (?) ]
cost_t t3 = Inf;
while (!h3.empty()) {
EdgeEvent e = h3.min();
int x = e.u1, y = es[e.to].to; // e.to is some edge id.
if (sf1[x] != sf1[y]) {
t3 = e.tm;
break;
} else h3.pop();
}
// delta4 : expand
cost_t t4 = Inf;
if (!h4.empty()) t4 = h4.min();
// -- events --
cost_t tn = min(min(t1, t2), min(t3, t4));
assert(tc1 <= tn && tn < Inf);
tc1 = tn;
if (tc1 == ev1.tm) {
int x = ev1.id;
if (x != rt) rematch(x, 0);
return true;
}
while (!h2.empty() && h2.min().tm == tc1) {
int x = h2.min().u1, y = h2.min().to;
if (grow(x, y)) return true; // `grow` function will call `h2.erase(by)`.
}
while (!h3.empty() && h3.min().tm == tc1) {
int x = h3.min().u1, eid = h3.min().to;
int y = es[eid].to; h3.pop();
if (sf1[x] == sf1[y]) continue;
contract(x, y, eid);
}
while (!h4.empty() && h4.min() == tc1) {
int b = h4.argmin(); h4.pop();
expand(b);
}
return false;
}
private:
void initialize() {
que = Queue<int>(N);
mt.assign(S, 0);
lnk.assign(S, {0, 0});
lab.assign(S, kFree);
bas.resize(S); for (int u = 1; u < S; ++u) bas[u] = u;
sf1.resize(S); for (int u = 1; u < S; ++u) sf1[u] = u;
pot.resize(S);
nd.resize(S); for (int b = 1; b < S; ++b) nd[b] = Node(b);
ub.resize(B); for (int i = 0; i < B; ++i) ub[i] = N + B - i;
ubi = B;
// for O(nm log n) implementation
reset_time();
tc.resize(S);
sl.resize(S); for (int i = 0; i < S; ++i) sl[i] = Inf;
bf.assign(S, 0);
son.assign(S, 0);
lz.assign(S, 0);
grp.resize(S); for (int i = 0; i < S; ++i) grp[i] = i;
}
void set_potential() {
for (int u = 1; u <= N; ++u) {
cost_t mx = 0;
for (int eid = ofs[u]; eid < ofs[u + 1]; ++eid) {
mx = max(mx, es[eid].cst);
}
pot[u] = mx >> 1;
}
}
void find_maximal_matching() {
// Find a maximal matching naively.
for (int u = 1; u <= N; ++u) if (!mt[u]) {
for (int eid = ofs[u]; eid < ofs[u + 1]; ++eid) {
auto& e = es[eid]; int v = e.to;
if (mt[v] > 0 || reduced_cost(u, v, e) > 0) continue;
mt[u] = v; mt[v] = u;
break;
}
}
}
private:
const int N, B, S; // N = |V|, B = (|V| - 1) / 2, S = N + B + 1
vector<int> ofs;
vector<Edge> es;
Queue<int> que;
vector<int> mt, sf1, bas;
vector<Link> lnk;
vector<Label> lab;
vector<cost_t> pot;
vector<int> ub; int ubi;
vector<Node> nd;
// for O(nm log n) implementation
vector<int> son, grp;
vector<cost_t> tc, lz, sl;
vector<int> bf;
cost_t tc1;
Event ev1;
ModifiableHeap<EdgeEvent> h2;
ModifiableHeaps<EdgeEvent> hps;
FastHeap<EdgeEvent> h3;
ModifiableHeap<cost_t> h4;
};
} // namespace general_weighted_matching_internal
template <class TotalCost> struct weighted_matching_result {
TotalCost wt{};
std::vector<std::pair<int, int>> es;
};
/// @brief Compute a maximum-weight matching in a general undirected graph.
/// A primal-dual Edmonds search grows an alternating forest through zero
/// reduced-cost edges. Odd cycles are contracted into blossoms; dual changes
/// schedule grow, contract, and expand events. Gabow's heavy-blossom grouping
/// and meldable heaps make all slack/event maintenance near-linear per search.
template <class Cost, class TotalCost = Cost>
weighted_matching_result<TotalCost> general_weighted_matching(
int n1,
const std::vector<std::tuple<int, int, Cost>> &es1) {
using solver = general_weighted_matching_internal::
gabow_maximum_weighted_matching<Cost, TotalCost>;
std::vector<typename solver::InputEdge> es;
es.reserve(es1.size());
for (auto [a1, b1, wt] : es1) {
assert(0 <= a1 && a1 < n1);
assert(0 <= b1 && b1 < n1);
assert(a1 != b1 && wt >= Cost{});
es.push_back({a1 + 1, b1 + 1, wt});
}
solver mat(n1, es);
weighted_matching_result<TotalCost> res;
res.wt = mat.maximum_weighted_matching(res.es);
for (auto &[a1, b1] : res.es) {
--a1;
--b1;
}
return res;
}
} // namespace noya