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

最大流的定义及基础解法这里就不再赘述了 有需要的可以看之前写的~~劣质~~文章[EK](https://www.luogu.org/blog/ONE-PIECE/wang-lao-liu-di-zong-jie)与[Dinic](https://www.luogu.org/blog/ONE-PIECE/wang-lao-liu-jiang-xie-zhi-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），结束算法。 先给一些数组定义： ``` int dep[13000],gap[13000]; //dep[i]表示节点i的深度，gap[i]表示深度为i的点的数量 ``` bfs部分 ``` 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： ``` 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差不多。 只是多几句话来统计深度对应的点数。 ### 这里要解释一下为什么要统计每个深度的点数： ### 为了优化！ ？ 统计每个深度对应点数只为了这句话： ``` if(gap[dep[u]]==0)dep[s]=n+1; ``` __ 因为我们是按照深度来往前走的，路径上的点的深度一定是连续的，而t的深度为0，如果某个深度的点不存在，那么我们就无法到达t了 __ 此时直接结束算法可以大大节省时间。 接下来是主函数： ``` int ISAP(){ maxflow=0; bfs(); while(dep[s]<n)dfs(s,inf); return maxflow; } ``` ## 由于dfs部分~~（有点）~~相当抽象，可以结合图示理解。 以下图为例：（反向边没有画，但是有反向边）  先bfs分层：（蓝色字体表示层数）  按照深度，我们只能走S->1->5->6->t这条增广路，流量为3 在dfs过程中，从5传到6时flow==4而在6号点的used==3，所以此时不会直接返回（见上方的代码），而是把6号点的深度加1，同理s,1,5的也要加1，也就是说，跑了这条增广路后，图变成了：  看见了吗？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的情况与上次的那些点相同，深度会改变。 跑了这条路之后，图变成了：  图中出现了断层（深度为4的点没有了），所以一定不存在增广路了，可以直接结束算法了。 # 刚才那张图如果用Dinic会怎么样？ ## 三遍bfs，两遍dfs！ ## 与之相比，ISAP真是太高效了！ 这里顺便说一下为什么终止条件是dep[s]<n 从上面的过程可以看出：每走一遍增广路，s的层数就会加1，如果一直没有出现断层，最多跑n-dep（刚bfs完时s的深度）条增广路（因为一共就n个点） 完整代码： ``` #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](https://www.luogu.org/blog/ONE-PIECE/wang-lao-liu-jiang-xie-zhi-dinic) ） 只需要改一点即可： ``` 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。** ------------ # 先用图解来表示算法的运行过程： 以下图为例（反向边没画出来）  我们用深蓝色的字表示每个点的余流，用淡蓝色的字表示每个点的高度。 刚开始肯定是：  我们把s点的高度改成6（一共6个点）  然后从往s拼命地灌水，水流会顺着还没满的边向周围的节点流  此时1,3节点还有余流，于是更新1,3号点高度，更新为与它相邻且最低的点的高度+1。  此时我们会对1号节点推流。  更新有余流的1,2号节点的高度  为什么1的高度变成了7呢？因为1与s有反向边，因为刚才推过流，反向边有了边权，与1连通的最低的点就是s。 然后是对3节点推流（更新高度是同时进行的只是刚才为了便于理解而分开画了）  然后是2号点推流，注意t点不能更新高度（否则好不容易流过去的水就会往回流了）  4号点推流  2号点推流  1号点推流  3号点推流  除了s和t以外，所有节点都无余流，算法结束，最大流就是t的余流。 完整代码: ``` // 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 【模板】最大流 加强版 / 预流推进（毒瘤题）](https://www.luogu.org/problemnew/show/P4722)更没希望。 那为什么这么慢，我们还要学呢？ 为了下面的讲解做铺垫。 ------------ ------------ # 最高标号预流推进（HLPP） ## 算法步骤： 1.先从t到s反向bfs，使每个点有一个初始高度 2.从s开始向外推流，将有余流的点放入优先队列 3.不断从优先队列里取出高度最高的点进行推流操作 4.若推完还有余流，更新高度标号，重新放入优先队列 5.当优先队列为空时结束算法，最大流即为t的余流 ## 与基础的余流推进相比的优势： 通过bfs预先处理了高度标号，并利用优先队列~~（闲着没事可以手写堆）~~使得每次推流都是高度最高的顶点，以此减少推流的次数和重标号的次数。 ## 优化： **和ISAP一样的gap优化**，如果某个高度不存在，将所有比该高度高的节点标记为不可到达（使它的高度为n+1，这样就会直接向s推流了）。 代码： ``` #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; } ``` 这次终于可以过掉那道[毒瘤题](https://www.luogu.org/problemnew/show/P4722)了 ------------ # 总结： ISAP的时间复杂度为O(n^2m),而HLPP的时间复杂度O(n^2sqrt(m)),但由于**HLPP常数过大**，在**纯随机数据**（模板题）下跑得还没有ISAP快，一般来说，用ISAP就够了。 证据： [ISAP模板题记录](https://www.luogu.org/record/show?rid=10909476) [HLPP模板题记录](https://www.luogu.org/record/show?rid=10894011) ~~当然也有可能只是我写的HLPP常数太大了。~~ ------------ 终于写完了。 # 求赞