网络最大流, 算法效率军备竞赛
网络最大流·Advanced
历史背景
对于网络流的初步, 半年前我写过一篇文章总结, Dinic 先生的两个算法. 在后来的日子里, 陆续有别人学习了最大流算法, 在议论着算法效率问题, 其中胜出者是@sxy, 他用 Dinic 某个神奇优化以最慢的一个点的时间为
我只有
55ms 哪个点
我的意思是总时间
然后年少轻狂的我拿着我的算法交了一发加强版, 只有
本以为这个事就这么过去了, 但是随着sxy用新学的 ISAP 以
但是光是用阴招算什么好汉, 我思前想后, 要想超越sxy, 抢走他的最优解比登天还难, 谁都知道模板题的常数已经卷到什么程度了, 所以我下定决心要过掉加强版, 这起码比优化到最优解简单.
ISAP
简介
ISAP Algorithm, Improved Shortest Augmenting Paths Algorithm, 是一种网络最大流算法, 翻译成中文是:
由于 Dinic 在算法竞赛中完全够用, 被卡机率非常小, 如果被卡了, ISAP 也能解决这些丧心病狂的题目, 所以论实用性, 学习 ISAP 的性价比还是很高的.
古人云: "知己知彼, 百战不殆", 我特意在sxy的博客中学习了 ISAP 算法, 这里挂一下友链.
从 Dinic 说起
先回顾一下 Dinic 的过程, 每次从
ISAP 算法就是从分层图的思想中优化过来的, 将多次分层优化成了单次分层.
至于这样为什么正确, 其实要从分层的意义开始分析. Dinic 算法之所以先分层再推流, 是为了让流单向流动, 防止两个点来回踢皮球, 或者环形流的出现, 也就是说我们只需要给每个点分一个先后, 或者说海拔高度, 然后让它们的流自行流入 T 即可.
那为什么要多次分层呢? 因为每次推流完成之后, 会有一些边的余量变成
前面提出一个问题, 为什么深度只能变深而不能变浅. 考虑每个点的深度是如何被更新的, 由于 BFS 的特性, 每个点一定是被起点最浅的入边更新的, 如果这条边不连通了, 那么这个点的深度就要被别的边更新, 而别的入边起点的深度一定不会更浅, 所以一个点的深度只能增加.
还有一个显而易见的结论, 如果规定流只能从浅到深单向移动, 则一个点的流只能直接流向比它深度大
接下来思考为什么可以不用重新搜索就维护每个点的深度, 从一个点的深度增加的条件考虑, 两种原因:
- 情况一: 所有入边均不再连通
由于前面证明了本轮推流的每条入边的起点的深度都是这个点的深度减
- 情况二: 所有入边的起点的深度增加
这种情况可以分析更新这些起点的深度增加的原因, 并且追溯到所有情况二的边界, 这些点深度增加的原因是情况一, 也就是所有入边都被榨干.
情况一和情况二, 总结来说就是一个点的深度增加的充要条件, 是在当前分层情况的可行边种,
一个点的深度增加后, 它的入边的起点也会变成其它深度更大的点 (深度比它的深度小
算法核心
前置结论证明完了, 接下来就是 ISAP 的精髓.
由于无法快速判断
这时, 推流的规则变成了将一个点的流推向深度比这个点小
证明了每次只有特定节点的深度会改变, 加上我们可以快速判断每个点的深度是否增加, 就可以动态地维护每个点的深度, 避免了每次 DFS 前都要跑 BFS.
判断了深度是否增加, 那如果增加, 增加多少呢?
事实上, 我们可以增加
Gap
Gap, 缺口, 豁口, 裂口.
建立一个数组,
我们把这种判断称作: "Gap 优化".
复杂度分析
在
因此, ISAP 的复杂度是
实现
实现起来非常简单, 细节也非常少, 只打了半小时就一遍过了, 但是跑得巨慢, 可能是评测机的锅 (评测机: MMP)
简单啊概括一下, 一遍 BFS 十分清新, 就是一遍 BFS.
DFS 的流程也比较简单, 符合条件的边就访问, 能推流就推流, 真正推到 T 点的流量会占用容量, 其余的多余流量回收, 推到别的出边.
由于一次只推
由于答案一般较小, 这种反复不会影响效率表现.
因为单次 DFS 源点流量没有爆 unsigned, 所以汇点实际的到的流量一定也不会爆 unsigned, 为了效率, 我们就用 unsigned 进行中间运算, 只有 unsigned long long 求每次 DFS 的返回值总和.
unsigned Hd(0), Tl(0), Gap[10005], m, n, Cnt(0), C, D, t;
unsigned long long Ans;
struct Node;
struct Edge {
Node *To;
Edge *Nxt;
unsigned Contain;
}E[10005], *CntE(E);
struct Node {
Edge *Fst;
unsigned Dep;
}N[205], *Que[205], *A, *B, *S, *T;
unsigned DFS(Node *x, unsigned Come) {
if(x == T) return Come;
register Edge *Sid(x->Fst);
register unsigned Go(0), Tmp;
while (Sid) {
if((Sid->To->Dep + 1 == x->Dep) && (Sid->Contain)) {
Tmp = DFS(Sid->To, min(Come, Sid->Contain));
Come -= Tmp;
Go += Tmp;
Sid->Contain -= Tmp;
(Sid + 1)->Contain += Tmp;
if(!Come) return Go;
}
Sid = Sid->Nxt;
}
if(!(--Gap[(x->Dep)++])) S->Dep = n + 1;
++Gap[x->Dep];
return Go;
}
signed main() {
n = RD(), m = RD(), S = N + RD(), T = N + RD();
for (register unsigned i(1); i <= m; ++i) {
A = N + RD(), B = N + RD(), C = RD();
(++CntE)->Nxt = A->Fst;
A->Fst = CntE;
CntE->To = B;
CntE->Contain = C;
(++CntE)->Nxt = B->Fst;
B->Fst = CntE;
CntE->To = A;
}
T->Dep = 1, Que[++Tl] = T;
register Node *x;
while(Hd < Tl) {
x = Que[++Hd];
register Edge *Sid(x->Fst);
while (Sid) {
if(!(Sid->To->Dep)) {
++Gap[Sid->To->Dep = x->Dep + 1];
Que[++Tl] = Sid->To;
}
Sid = Sid->Nxt;
}
}
while (S->Dep <= n) Ans += DFS(S, 0xffffffff);
printf("%llu\n", Ans);
return Wild_Donkey;
}当然, 和sxy的 ISAP 一样, 我的 ISAP 同样在加强版中得到了
HLPP
简介
要想通过加强版, 我们这些凡人都需要用到一个算法, HLPP.
HLPP Algorithm, 缩写源自于 Highest label Preflow Push Algirithm (最高标号预流推进算法)
在 Wikipedia 中, 预流推进算法的名称是: Push–relabel Maximum Flow Algorithm. 分析百科文本中对时间复杂度
所以根据这么详细的名字, 我们就可以顾名思义一下了:
-
通过给节点编号来确定单向流的方向, 这一点和 ISAP/Dinic 一样.
-
每次选择最高标号的节点来推流, 可以让本次推流中, 当前推流的节点在推流之后不会接受到任何流
所以这应该更加类似于在 ISAP 的基础上引入了推流节点的选择, 用 BFS 代替了 DFS. 让程序从跟着流从
算法思路
仍然是从
然后将
这时类似于 BFS, 每次弹出堆顶一个点, 这个点一定是之前有流量推过来的点, 用它的流量往高度比它小
-
推完所有合法的出边后该点仍有流量没有推走, 这时增加它的高度, 重新入队, 下一轮会尝试将流推到更高的点中. 因为既然这个点弹出前是标号最大, 那么它增加后再次入队, 这期间没有其它点加入, 它一定还是堆顶. 增加的数量也没必要非得是
1 , 因为高度增加是为了流向其它点, 所以取出边中, 本次高度不足以流向的点的高度最小值再加上1 作为本次高度增加后的值. 这样既不会漏掉任何可行的出边, 又不会出现大量无意义重复运算. -
推完了, 不用进行任何操作, 继续取堆顶推流即可.
如果推到了
最后的答案就是
由于根号的复杂度过于玄学, 所以我对此非常疑惑, 不过时间复杂度一定和每个点的入队次数有关, 也和优先队列有关, 所以暂不证明.
Gap
Gap, 缺口, 豁口, 裂口.
和 ISAP 不同的是, HLPP 的高度不是分层图的深度, 所以 Gap 优化不能和 ISAP 一样.
在 HLPP 中, 高度在堆顶之上的点, 除了
实现方面每次一个
实现
代码细节较多, 码量没有提高太多, 由于写代码时对算法理解程度不高, 花了不少时间才过了题, 但是总体逻辑比较清晰.
单调队列方面比较绕, 因为避免将整个结构体作为队列元素, 我们只维护指针即可, 但是指针的运算符没法或是说, 我不会单独定义, 所以我开了一个新的结构体 Que, 只有一个成员变量, Node 指针 Que 的运算符定义为按
因为 HLPP 是基于 ISAP 修改得到的, 所以代码会有雷同, 但不同的是: 它可以 AC 加强版 (虽然很慢).
和 ISAP 一样, 为了方便 BFS, 我们将
unsigned Hd(0), Tl(0), Gap[1205], m, n, Cnt(0), C, D, t, Tmp(0);
struct Node;
struct Edge {
Node *To;
Edge *Nxt;
unsigned Contain;
}E[240005], *CntE(E - 1);
struct Node {
Edge *Fst;
unsigned Dep, Contain;
}N[1205], *Qu[1205], *A, *B, *S, *T;
struct Que {
Node *P;
inline const char operator<(const Que &x) const {
return this->P->Dep < x.P->Dep;
}
};
priority_queue <Que> Q;
signed main() {
n = RD(), m = RD(), S = N + RD(), T = N + RD();
for (register unsigned i(1); i <= m; ++i) {
A = N + RD(), B = N + RD(), C = RD();
if(A == B) continue;
(++CntE)->Nxt = A->Fst;
A->Fst = CntE;
CntE->To = B;
CntE->Contain = C;
(++CntE)->Nxt = B->Fst;
B->Fst = CntE;
CntE->To = A;
}
T->Dep = 1, Qu[++Tl] = T;
register Node *x;
while(Hd < Tl) {
x = Qu[++Hd];
register Edge *Sid(x->Fst);
while (Sid) {
if((!(Sid->To->Dep)) && (!(Sid->Contain))) {
++Gap[Sid->To->Dep = x->Dep + 1];
Qu[++Tl] = Sid->To;
}
Sid = Sid->Nxt;
}
}
--Gap[S->Dep];
++Gap[S->Dep = n + 1];
register Que Pu;
register Edge *Sid(S->Fst);
while (Sid) {
if(Sid->Contain) {
if(Sid->To != T && (!(Sid->To->Contain))) {
Pu.P = Sid->To;
Q.push(Pu);
}
Sid->To->Contain += Sid->Contain;
(Sid + 1)->Contain = Sid->Contain;
Sid->Contain = 0;
}
Sid = Sid->Nxt;
}
while(Q.size()) {
x = (Q.top()).P, Q.pop();
register unsigned Real;
Sid = x->Fst;
Tmp = 0x3f3f3f3f;
while(Sid) {
if(Sid->Contain) {
if(Sid->To->Dep + 1 == x->Dep) {
Real = min(x->Contain, Sid->Contain);
if(!Real) {Sid = Sid->Nxt; continue;}
x->Contain -= Real;
Sid->Contain -= Real;
E[(Sid - E) ^ 1].Contain += Real;
if(Sid->To != S && Sid->To != T && (!(Sid->To->Contain))) {
Pu.P = Sid->To, Q.push(Pu);
}
Sid->To->Contain += Real;
if(!(x->Contain)) break;
} else Tmp = min(Tmp, Sid->To->Dep);
}
Sid = Sid->Nxt;
}
if(x->Contain) {
if(!(--Gap[x->Dep])) {
for (register unsigned i(1); i <= n; ++i) {
if(N + i != S && N + i != T && N[i].Dep > x->Dep) {
N[i].Dep = n + 2;
}
}
}
++Gap[x->Dep = Tmp + 1];
Pu.P = x;
Q.push(Pu);
}
}
printf("%u\n", T->Contain);
return Wild_Donkey;
}