究极的最大流算法ISAP与HLPP——钱逸凡

发表于 | 分类于 | ℃热度

本文搬运自洛谷大佬钱逸凡的博客

究级的最大流算法：ISAP与HLPP

最大流的定义及基础解法这里就不再赘述了

有需要的可以看之前写的劣质文章EK与Dinic求解最大流的方法了

前言

最大流算法目前有增广路算法和预流推进算法两种，增广路算法的思想是不断地寻找增广路来增大最大流，而预流推进算法的思想是……后面再讲。

先从最熟悉的增广路算法开始讲。

ISAP(Improved Shortest Augumenting Path）：

其实Dinic已经足够高效了，但那些（闲着没事干）的计算机科学家们对Dinic的效率仍不满足，于是就有了ISAP算法。

在Dinic中我们每次增广前都进行了一次bfs来初始化每个点的深度，虽然一次标号增广了很多条路，但是我们还是很有可能要跑很多遍bfs导致效率不高。

那有没有什么办法只跑一次bfs呢?
那就是ISAP算法了！

ISAP运行过程：

1.从t到s跑一遍bfs，标记深度（为什么是从t到s呢？后面会讲）

2.从s到t跑dfs，和Dinic类似，只是当一个点跑完后，如果从上一个点传过来的flow比该点的used大（对于该点当前的深度来说，该点在该点以后的路上已经废了），则把它的深度加1，如果出现断层（某个深度没有点），结束算法

3.如果操作2没有结束算法，重复操作2

好抽象啊。

什么原理呢？

ISAP其实与Dinic差不多，但是它只跑一遍bfs，但是每个点的层数随着dfs的进行而不断提高（这样就不用反复跑bfs重新分层了），当s的深度大于n时（这就是为什么bfs要从t到s），结束算法。

先给一些数组定义：

1
2
int dep[13000],gap[13000];
//dep[i]表示节点i的深度，gap[i]表示深度为i的点的数量

bfs部分:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
void bfs(){
    memset(dep,-1,sizeof(dep));
    memset(gap,0,sizeof(gap));
    dep[t]=0;
    gap[0]=1;//
    queue<int>q;
    q.push(t);
    while(!q.empty()){
        int u=q.front();
        q.pop();
        for(int i=head[u];i;i=node[i].next){
            int v=node[i].v;
            if(dep[v]!=-1)continue;//防止重复改某个点
            q.push(v);
            dep[v]=dep[u]+1;
            gap[dep[v]]++;
        }
    }
    return;
}//初始化bfs，从t到s搜出每个点初始深度

可以看出ISAP里的bfs对 边权!=0 这一条件并没有限制，只要在后面的dfs里判断边权就好了。

接下来是dfs：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
int dfs(int u,int flow){
    if(u==t){//可以到达t
        maxflow+=flow;
        return flow;
    }
    int used=0;
    for(int i=head[u];i;i=node[i].next){
        int d=node[i].v;
        if(node[i].val&&dep[d]+1==dep[u]){
            int mi=dfs(d,min(node[i].val,flow-used));
            if(mi){
                node[i].val-=mi;
                node[i^1].val+=mi;
                used+=mi;
            }
            if(used==flow)return used;
        }
    }
    //前半段和Dinic一模一样
    //如果已经到了这里，说明该点出去的所有点都已经流过了
    //并且从前面点传过来的流量还有剩余
    //则此时，要对该点更改dep
    //使得该点与该点出去的点分隔开
    --gap[dep[u]];
    if(gap[dep[u]]==0)dep[s]=n+1;//出现断层，无法到达t了
    dep[u]++;//层++ 
    gap[dep[u]]++;//层数对应个数++
    return used; 
}

感觉和Dinic差不多。

只是多几句话来统计深度对应的点数。

这里要解释一下为什么要统计每个深度的节点数：
为了优化！

统计每个深度对应点数只为了这句话：

1
if(gap[dep[u]]==0)dep[s]=n+1;
因为我们是按照深度来往前走的，路径上的点的深度一定是连续的，而t的深度为0，如果某个深度的点不存在，那么我们就无法到达t了

此时直接结束算法可以大大节省时间。

接下来是主函数：

1
2
3
4
5
6
int ISAP(){
    maxflow=0;
    bfs();
    while(dep[s]<n)dfs(s,inf);
    return maxflow;
}

由于dfs部分（有点）相当抽象，可以结合图示理解。

以下图为例：（反向边没有画，但是有反向边）

33500

先 bfs 分层：（蓝色字体表示层数）

33501

按照深度，我们只能走S->1->5->6->t这条增广路，流量为3

在dfs过程中，从5传到6时$flow==4$而在6号点的$used==3$，所以此时不会直接返回（见上方的代码），而是把6号点的深度加1，同理s,1,5的也要加1，也就是说，跑了这条增广路后，图变成了：

33504

看见了吗？s->1->5->6->t这条路被封了，但是s->1->2->3->4->t这条路却能走了

按照深度，这次肯定是走s->1->2->3->4->t这条路。

这次的情况和上次不同，由于从2传到3时的flow==3而之后的路径的used都是3，所以2,3,4号点的深度不会改变，而1,s的情况与上次的那些点相同，深度会改变。

跑了这条路之后，图变成了：

33505

图中出现了断层（深度为4的点没有了），所以一定不存在增广路了，可以直接结束算法了。

刚才那张图如果用Dinic会怎么样？

三遍bfs，两遍dfs！

与之相比，ISAP真是太高效了！

这里顺便说一下为什么终止条件是dep[s]<n

从上面的过程可以看出：每走一遍增广路，s的层数就会加1，如果一直没有出现断层，最多跑n-dep（刚bfs完时s的深度）条增广路（因为一共就n个点）

完整代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
#include<queue>
using namespace std;
const int inf=1<<30;
int cnt=1,head[13000];
int n,m,s,t;
inline int Read(){
    int x=0;
    char c=getchar();
    while(c>'9'||c<'0')c=getchar();
    while(c>='0'&&c<='9')x=x*10+c-'0',c=getchar();
    return x;
}
struct Node{
    int v;
    int next;
    int val;
}node[250000];
inline void addedge(int u,int v,int val){
    node[++cnt].v=v;
    node[cnt].val=val;
    node[cnt].next=head[u];
    head[u]=cnt;
}
int dep[13000],gap[13000];
void bfs(){
    memset(dep,-1,sizeof(dep));
    memset(gap,0,sizeof(gap));
    dep[t]=0;
    gap[0]=1;
    queue<int>q;
    q.push(t);
    while(!q.empty()){
        int u=q.front();
        q.pop();
        for(int i=head[u];i;i=node[i].next){
            int v=node[i].v;
            if(dep[v]!=-1)continue;
            q.push(v);
            dep[v]=dep[u]+1;
            gap[dep[v]]++;
        }
    }
    return;
}
int maxflow;
int dfs(int u,int flow){
    if(u==t){
        maxflow+=flow;
        return flow;
    }
    int used=0;
    for(int i=head[u];i;i=node[i].next){
        int d=node[i].v;
        if(node[i].val&&dep[d]+1==dep[u]){
            int mi=dfs(d,min(node[i].val,flow-used));
            if(mi){
                node[i].val-=mi;
                node[i^1].val+=mi;
                used+=mi;
            }
            if(used==flow)return used;
        }
    }
    --gap[dep[u]];
    if(gap[dep[u]]==0)dep[s]=n+1;
    dep[u]++;
    gap[dep[u]]++;
    return used; 
}
int ISAP(){
    maxflow=0;
    bfs();
    while(dep[s]<n)dfs(s,inf);
    return maxflow;
}
int main(){
    int i=1;
    n=Read(),m=Read(),s=Read(),t=Read();
    int u,v,w;
    for(i=1;i<=m;i++)u=Read(),v=Read(),w=Read(),addedge(u,v,w),addedge(v,u,0);  
    printf("%d",ISAP());
    return 0;
}

当然，这东西也可以当前弧优化（原理见 Dinic ）

只需要改一点即可：

1
2
3
4
while(dep[s]<n){
    memcpy(cur,head,sizeof(head));
    dfs(s,inf); 
    }

dfs部分的修改就不放代码了，应该都会。

接下来就是相当抽象的预流推进算法了

  • 在讲解算法之前，我们先来回顾一下最大流是干嘛的：
  • 水流从一个源点s通过很多路径，经过很多点，到达汇点t，问你最多能有多少水能够到达t点。
  • 在刚面对这个问题时，你或许会有这样的思路：
  • 先从s往相邻点拼命灌水，然后让水不停地向t流，能流多少是多少。

这就是预流推进：

能推流，就推流，最终t有多少水，最大流就是多少。

预流推进算法的思想：

在讲思想之前，先引入一个概念：

余流：

  • 听名字就知道是什么意思了，即每个点当前有多少水。

预留推进算法的思想是：

  • 1.先假装s有无限多的水(余流)，从s向周围点推流（把该点的余流推给周围点，注意：推的流量不能超过边的容量也不能超过该点余流），并让周围点入队（ 注意：s和t不能入队 ）
  • 2.不断地取队首元素，对队首元素推流
  • 3.队列为空时结束算法，t点的余流即为最大流。
  • 上述思路是不是看起来很简单，也感觉是完全正确的？
  • 但是这个思路有一个问题，就是可能会出现两个点不停地来回推流的情况，一直推到TLE。
  • 怎么解决这个问题呢？
  • 给每个点一个高度，水只会从高处往低处流。在算法运行时， 不断地对有余流的点更改高度 ，直到这些点全部没有余流为止。
为什么这样就不会出现来回推流的情况了呢？

  • 当两个点开始来回推流时，它们的高度会不断上升，当它们的高度大于s时，会把余流还给s。

  • 所以在开始预流推进前要先把s的高度改为n（点数），免得一开始s周围那些点就急着把余流还给s。

先用图解来表示算法的运行过程：

以下图为例（反向边没画出来）

33590

我们用深蓝色的字表示每个点的余流，用淡蓝色的字表示每个点的高度。

刚开始肯定是：

33594

我们把s点的高度改成6（一共6个点）

33595

然后从往s拼命地灌水，水流会顺着还没满的边向周围的节点流

33598

此时1,3节点还有余流，于是更新1,3号点高度，更新为与它相邻且最低的点的高度+1。

33599

此时我们对1号节点推流。

33604

更新有余流的1,2号节点的高度

33605

为什么1的高度变成了7呢？因为1与s有反向边，因为刚才推过流，反向边有了边权，与1连通的最低的点就是s。

然后是对3节点推流（更新高度是同时进行的只是刚才为了便于理解而分开画了）

33606

然后是2号点推流，注意t点不能更新高度（否则好不容易流过去的水就会往回流了）

33607

4号点推流

33609

2号点推流

33611

3号点推流

33630

除了s和t以外，所有节点都无余流，算法结束，最大流就是t的余流。

完整代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
// luogu-judger-enable-o2
//效率相当低，吸氧气都过不了模板
#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
#include<queue>
using namespace std;
const int inf=1<<30;;
int n,m,s,t;
struct Node{
    int v;
    int val;
    int next;
}node[205202];
int top=1;
int inque[13000];
int head[13000];
int h[13000];
int e[13000];//每个点的高度，每个点的剩余（多余）流量
inline int Read(){
    int x=0;
    char c=getchar();
    while(c>'9'||c<'0')c=getchar();
    while(c>='0'&&c<='9')x=x*10+c-'0',c=getchar();
    return x;
}//读入优化 
inline int min(int a,int b){
    return a<b?a:b;
}
inline void addedge(int u,int v,int val){
    node[++top].v=v;
    node[top].val=val;
    node[top].next=head[u];
    head[u]=top;
}
inline void add(int u,int v,int val){
    addedge(u,v,val);
    addedge(v,u,0);
}
int maxflow;
inline void push_(int u,int i){
    int mi=min(e[u],node[i].val);
    node[i].val-=mi;
    node[i^1].val+=mi;
    e[u]-=mi;
    e[node[i].v]+=mi;
}//从u到v推流,i为u与v相连的那一条边的编号 
inline bool relabel(int u){
    int mh=inf;
    register int i;
    for(int i=head[u];i;i=node[i].next){
        int d=node[i].v;
        if(node[i].val){
            mh=min(mh,h[d]);
        }
    }
    if(mh==inf)return 0;//无可以流的路径
    h[u]=mh+1;//把u点拉到可流点以上（使水流能从u流向周围的点）
    return 1; 
}//判断点u是否能流向周围点
int push_relabel(){
    queue<int>q;
    h[s]=n;//把源点的高度设成n（图中的点数） 
    for(register int i=head[s];i;i=node[i].next){
        int d=node[i].v;
        int val=node[i].val;
        if(val){
            node[i].val-=val;
            node[i^1].val+=val;
            e[s]-=val;
            e[d]+=val;
            if(inque[d]==0&&d!=t&&d!=s&&relabel(d)){
                q.push(d);
                inque[d]=1;
            }
        }//给所有与s相连的点推入满流 
    }
    while(!q.empty()){
        int u=q.front();
        q.pop();
        inque[u]=0;
        for(int i=head[u];i;i=node[i].next){
            int d=node[i].v;
            if(e[u]==0)break;
            if(node[i].val==0)continue;
            if(h[u]==h[d]+1) push_(u,i);
            if(e[d]!=0&&d!=t&&d!=s){
                if(inque[d]==0&&relabel(d)){
                    q.push(d);
                    inque[d]=1;
                }
            }
        }
        if(e[u]&&inque[u]==0&&relabel(u)){
            q.push(u);
            inque[u]=1;
        }
    }
    return e[t];
}//push_relabel
int main(){
    n=Read(),m=Read(),s=Read(),t=Read();
    register int i;
    int u,v,val;
    for(i=1;i<=m;i++)u=Read(),v=Read(),val=Read(),add(u,v,val);
    printf("%d",push_relabel());
    return 0;
}

这个预流推进算法相当慢，至于P4722 【模板】最大流 加强版 / 预流推进（毒瘤题）更没希望。

那为什么这么慢，我们还要学呢？

为了下面的讲解做铺垫。

最高标号预流推进（HLPP）

算法步骤：

  • 1.先从t到s反向bfs，使每个点有一个初始高度
  • 2.从s开始向外推流，将有余流的点放入优先队列
  • 3.不断从优先队列里取出高度最高的点进行推流操作
  • 4.若推完还有余流，更新高度标号，重新放入优先队列
  • 5.当优先队列为空时结束算法，最大流即为t的余流

与基础的余流推进相比的优势：

  • 通过bfs预先处理了高度标号，并利用优先队列（闲着没事可以手写堆）使得每次推流都是高度最高的顶点，以此减少推流的次数和重标号的次数。

优化：

  • 和ISAP一样的gap优化，如果某个高度不存在，将所有比该高度高的节点标记为不可到达（使它的高度为n+1，这样就会直接向s推流了）。

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
#include<queue>
#include<vector>
using namespace std;
inline int Read(){
    int x=0;
    char c=getchar();
    while(c>'9'||c<'0')c=getchar();
    while(c>='0'&&c<='9')x=x*10+c-'0',c=getchar();
    return x;
}
const int inf=1<<30;
int top=1,head[10100];
int n,m,s,t;
int e[10100],h[10100],cnth[20100];//每个点对应的余流，高度；每个高度有多少个点 
struct cmp{
    inline bool operator () (int a,int b) const{
        return h[a]<h[b];
    }
};
struct Node{
    int v;
    int val;
    int next;
}node[400100];
inline void addedge(int u,int v,int val){
    node[++top].v=v;
    node[top].val=val;
    node[top].next=head[u];
    head[u]=top;
}
inline void add(int u,int v,int val){
    addedge(u,v,val);
    addedge(v,u,0);
}
int inque[11000];
void bfs(){
    memset(h,0x3f,sizeof(h));
    h[t]=0;
    queue<int>qu;
    qu.push(t);
    while(!qu.empty()){
        int u=qu.front();
        qu.pop();
        inque[u]=0;
        for(int i=head[u];i;i=node[i].next){
            int d=node[i].v;
            if(node[i^1].val&&h[d]>h[u]+1){//反向跑 
                h[d]=h[u]+1;
                if(inque[d]==0){
                    qu.push(d);
                    inque[d]=1;
                }
            }
        }
    }
    return;
}
priority_queue<int,vector<int>,cmp>q;
inline void push_(int u){
    for(int i=head[u];i;i=node[i].next){
        int d=node[i].v;
        if(node[i].val&&h[d]+1==h[u]){//可以推流 
            int mi=min(node[i].val,e[u]);
            node[i].val-=mi;
            node[i^1].val+=mi;
            e[u]-=mi;
            e[d]+=mi;
            if(inque[d]==0&&d!=t&&d!=s){
                q.push(d);
                inque[d]=1;
            }
            if(e[u]==0)break;//已经推完了 
        }
    }
}//推流 
inline void relabel(int u){
    h[u]=inf;
    for(int i=head[u];i;i=node[i].next){
        int d=node[i].v;
        if(node[i].val&&h[d]+1<h[u]){
            h[u]=h[d]+1;
        }
    }
}//把u的高度更改为与u相邻的最低的点的高度加1 
int hlpp(){
    register int i;
    bfs();
    if(h[s]==0x3f3f3f3f)return 0;//s与t不连通 
    h[s]=n;
    for(i=1;i<=n;i++)if(h[i]<0x3f3f3f3f)cnth[h[i]]++;//统计各个高度的点数，注意不要让下标越界
    for(i=head[s];i;i=node[i].next){
        int d=node[i].v;    
        int mi=node[i].val;
        if(mi){
            e[s]-=mi;
            e[d]+=mi;
            node[i].val-=mi;
            node[i^1].val+=mi;
            if(d!=t&&inque[d]==0&&d!=s){
            q.push(d);
            inque[d]=1; 
            }
        }
    }//从s向周围点推流
    while(!q.empty()){
        int u=q.top();
        inque[u]=0;
        q.pop();
        push_(u);
        if(e[u]){//还有余流 
            cnth[h[u]]--;
            if(cnth[h[u]]==0){
                for(int i=1;i<=n;i++){
                    if(i!=s&&i!=t&&h[i]>h[u]&&h[i]<n+1){
                        h[i]=n+1;//标记无法到达 
                    }
                }
            }//gap优化 
            relabel(u);
            cnth[h[u]]++;
            q.push(u);
            inque[u]=1;
        }
    }
    return e[t];
} 
int main(){
    n=Read(),m=Read(),s=Read(),t=Read();
    register int i;
    int u,v,val;
    for(i=1;i<=m;i++)u=Read(),v=Read(),val=Read(),add(u,v,val);
    printf("%d",hlpp());
    return 0;
}

这次终于可以过掉那道毒瘤题了。

总结：

ISAP的时间复杂度为O(n^2m),而HLPP的时间复杂度O(n^2sqrt(m)),但由于HLPP常数过大，在纯随机数据（模板题）下跑得还没有ISAP快，一般来说，用ISAP就够了。

最后感谢钱逸凡大佬的经验共享QWQ

-------------本文结束(づ￣ 3￣)づ~感谢您的阅读(*╹▽╹*)-------------