[可能有用科技]从Push-Relable到HLPP
前言
应该没有毒瘤出题人会去卡 \(\textsf{Dinic}\) 和 \(\textsf{ISAP}\) 吧。
最大流
首先明确我们要干什么 —— 求最大流。
然后用那个最常见的比喻 : 源点是自来水厂向外输水,源点输出的所有水都需要汇入污水处理厂,但是水管是有流量上限的。
那么有一个非常 \(\textsf{naive}\) 的想法就是从源点狂暴加大水压疯狂灌水,反正源点水量无限,灌水灌到汇点无法得到更多流量。这时候相当于没有从源点到汇点的有空余流量的路径了,也就是没有 增广路,根据增广路定理,这时候的流量就是最大流了。
然后有两个人说,那么我直接从源点往外沿着增广路走出去不就行了吗,每次增广一条路径到汇点,最后肯定能求出最大流。
这两个人就是 \(\textsf{L.R.Ford,Jr.}\) 和 \(\textsf{D.R.Fulkerson}\),这就是他们提出的 \(\textsf{Ford–Fulkerson}\) 算法(更正确地说 : \(\textsf{Ford–Fulkerson}\) 方法),而 \(\textsf{Edmonds-Karp,Dinic,ISAP}\) 都是其具体的实现方式。
而这个方法仍然是根据路径进行增广,而且需要全程保证 流守恒性,于是有人觉得这实在是太逊了。
预流推进
首先提出两个新思路 :
- 放弃 “增广路”,直接对点维护
- 求最大流过程中暂时忽略流守恒
流量守恒去他*
考虑对于每个点维护流,在过程中忽略流守恒性,暴力灌满边然后让结点维护流量。套用上面的比喻,在每家修建一个无限大的储水池,那么就可以对每个点维护超额流量然后把流量传递给别的结点了。
具体过程如下 :
- 把与源点连通的点加满。
- 每次更新能把超额流量传出的点直到除了源汇点外全部点都没有超额流量
- 此时汇点超额流量就是最大流
但是这个过程是极不完善的,由于流量从正边流过去后可以从反边流回来,直接照着这个过程实现就会回流。
在学 \(\textsf{Dinic}\) 的时候就知道,为了防止回流,我们为每个点分配一个顶标,那么这里仿照这种方式,从汇点到源点分配高度,对于一个点 \(u\),只能将超额流量分给满足 \(h(u) = h(v) + 1\) 且连通的点 \(v\)。
这时我们就完成了一个合理的 推流(\(\bf{push}\)) ,但是不管是 \(\textsf{Dinic}\) 还是 \(\textsf{ISAP}\) 都需要调整顶标才能求出全部的增广路,这里我们为了使得超额流量流出足够多,需要一个 重标号(\(\bf{relable}\)) 的过程。
对于一个点 \(u\),如果当前标号下已经尽可能多地推流但是仍有剩余的超额流量,那么就从所有与其联通的点中找出一个高度最小的点 \(v\) 然后令 \(h(u) \leftarrow h(v) + 1\)。
最后为了不让源点流出的流量直接回到源点,我们将源点的初始高度设定为 \(n\)。
可以发现最后多余的无从推流的超额流量会沿反边回到源点,最后的图形态也是流量平衡的。
而这就是 \(\textsf{Push-Relable}\) 预流推进,复杂度 \(\mathcal{O} (n^2m)\) 和常用的 \(\textsf{Dinic}\) 一样。
重新整理一下思路 :
- 把和源点连通的点超额流量加满。
- 对于每个流量溢出的点 \(u\),向满足 \(h(u) = h(v) + 1\) 且连通的点 \(v\) 推流。
- 如果对于 \(u\) 所有可推流的点都推流后仍然溢出,更新 \(u\) 的高度为与 \(u\) 连通所有点中高度最小者的高度加一。
- 流量平衡后,汇点的超额流量就是最大流。
码量和其他(除了 HLPP) 的最大流差不多,常数略大。
但是比较好理解也好写。
Code :
int head[N],ecnt = -1;
struct Edge {
int nxt,to,cap;
}e[M << 1];
inline void AddEdge(int st,int ed,int cap) {
e[++ecnt] = (Edge) {head[st],ed,cap},head[st] = ecnt;
e[++ecnt] = (Edge) {head[ed],st,0},head[ed] = ecnt;
}
int n,m,s,t;
int h[N];
ll ex[N];
inline void init() {
repg(i,s) {
int v = e[i].to;
ex[v] += e[i].cap,ex[s] -= e[i].cap;
e[i ^ 1].cap = e[i].cap,e[i].cap = 0;
}
h[s] = n;
}
inline bool push(int ed) {
const int u = e[ed ^ 1].to,v = e[ed].to;
int f = std::min(ex[u],(ll)e[ed].cap);
ex[u] -= f,ex[v] += f;
e[ed].cap -= f,e[ed ^ 1].cap += f;
return ex[u];
}
inline void relable(int u) {
h[u] = INF;
repg(i,u) if(e[i].cap)
h[u] = std::min(h[u],h[e[i].to]);
++h[u];
}
inline ll PushRelable() {
bool flag = 1;
init();
while(flag) {
flag = 0;
rep(u,1,n) if(ex[u] > 0 && u != s && u != t) {
flag = 1;
bool re = 1;
repg(i,u) {
if(e[i].cap && h[u] == h[e[i].to] + 1)
re &= push(i);
if(!re) break;
}
if(re) relable(u);
}
}
return ex[t];
}
HLPP
\(\mathsf{HLPP}\),全称 最高标号预流推进。
通过上面的过程我们发现,先把标号高的结点的超额流量送出是更好的早死早超生,那么维护一个大根堆维护 \(h(u)\) 最大的可更新的 \(u\),每次选择堆顶元素推流。
然后为了不带 \(\log\),直接开 \(2n - 1\) 个桶维护高度。
但是由上面可以发现即使最后超额流量不流回源点是不对最大流产生影响的,于是只开 \(n\) 个桶即可。
似乎 HLPP 也有当前弧优化,但是这个咱不会。
可以看看 神 \(\mathrm{{\color{black}M}{\color{red}in\_25}}\) 的实现 : \(\texttt{Link}\)
复杂度 \(\mathcal{O} (n^2\sqrt{m})\)。
怎么说吧……就是在普通版的最大流也挺快的。
模板 : \(\texttt{Link}\)
Code :
int head[N],ecnt = -1;
struct Edge {
int nxt,to,cap;
}e[M << 1];
inline void AddEdge(int st,int ed,int cap) {
e[++ecnt] = (Edge) {head[st],ed,cap},head[st] = ecnt;
e[++ecnt] = (Edge) {head[ed],st,0},head[ed] = ecnt;
}
int n,m,s,t;
int h[N],ex[N];
int gap[N];
std::stack<int> buc[N];
int mxht;// max height
inline bool push(int u) {
bool pre = (u == s);
repg(i,u) {
const int v = e[i].to;
if(!e[i].cap || (!pre && h[u] != h[v] + 1))
continue;
int tf = pre ? e[i].cap : std::min(e[i].cap,ex[u]);
if (v != s && v != t && !ex[v])
buc[h[v]].push(v),mxht = std::max(mxht,h[v]);
ex[u] -= tf,ex[v] += tf;
e[i].cap -= tf,e[i ^ 1].cap += tf;
if(!ex[u]) return 0;
}
return 1;
}
inline void relable(int u) {
h[u] = INF;
repg(i,u) if(e[i].cap)
h[u] = std::min(h[u],h[e[i].to]);
if(++h[u] < n) {
buc[h[u]].push(u);
mxht = std::max(mxht,h[u]);
++gap[h[u]];
}
}
bool bfs() {
mems(h,63),h[t] = 0;
std::queue<int> q;q.push(t);
while(!q.empty()) {
int u = q.front();q.pop();
repg(i,u) {
int v = e[i].to;
if(e[i ^ 1].cap && h[v] > h[u] + 1) {
h[v] = h[u] + 1;
q.push(v);
}
}
}
return h[s] != 0x3f3f3f3f;
}
inline int GetHmax() {
while(buc[mxht].empty() && mxht > -1)
--mxht;
return (~mxht) ? buc[mxht].top() : 0;
}
inline int HLPP() {
if(!bfs()) return 0;
mems(gap,0);
rep(i,1,n) if(h[i] != 0x3f3f3f3f)
++gap[h[i]];
h[s] = n,push(s);
int u;
while((u = GetHmax())) {
buc[mxht].pop();
if(push(u)) {
if(!--gap[h[u]]) {
rep(i,1,n) if(i != s && i != t && h[i] > h[u] && h[i] < n + 1)
h[i] = n + 1;
}
relable(u);
}
}
return ex[t];
}