最大流及其增广路算法入门（Emonds-Karp & Dinic & ISAP）

一.最大流的定义.

网络：我们称一张带权有向图 G = { V , E } G=\{V,E\} G={V,E}为网络，当且仅当这张有向图中有两个特殊点$s,t$，其中$s$称为原点，$t$称为汇点.

容量：在一张网络上某条边$(x,y)$的边权称之为这条边的容量，这里记为$c(x,y)$.

流函数：对于每条边有一个函数$f(x,y)\le c(x,y)$，我们称$f$为流函数.

流量守恒：对于每一个点$x=/s$且$x=/t$，流函数必须满足：
$$\sum _{(i,x)\in E}f(i,x)=\sum _{(x,j)\in E}f(x,j)$$

我们称上述流函数的性质为流量守恒.

流函数的通俗理解：给定一张水网，有$n$个蓄水池和$m$条水管（每条水管仅连接两个水池），每条水管有一个容量限制，那么某一条水管流过的水量就是这条水管的流函数.

流量守恒的通俗理解：一个蓄水池流入的水量必须与它流出的水量相等.

网络的流量：一张网络 G = { V , E , s , t } G=\{V,E,s,t\} G={V,E,s,t}（$s$为源点，$t$为汇点）的流量$f(G)$为：
$$f(G)=\sum _{(s,i)\in E}f(x,i)=\sum _{(t,i)\in E}$$

网络流量的通俗理解：从源点流到汇点的水量.

最大流：一张网络的最大流即这张网络的最大流量（显然一张网络的流量是不唯一的）.

最大流的通俗理解：源点可以流出的最多的水，汇点可以得到的最多的水.


二.增广路算法.

增广路：网络中一条增广路定义为一条从源点到汇点的路径，使得这条路径上所有边的容量都大于$0$.

在解决最大流问题的算法中，有一类算法专门利用了增广路这个东西，来计算最大流.

这一类算法有一个共同特点：每次都会找一条或多条增广路增加流量，并通过这个方法计算最大流.

不给过通常情况下这种算法并不能直接保证得到的最终流量是最大流量，这就需要利用到这一类算法另一个基本的重要思想了——建立反向边.

原网络：即最初没有经过任何操作的网络.

残余网络：当每次用一条或多条增广路流完（把增广路上的流量限制减掉这一次增加的流量）之后，剩下来的网络称为残余网络.

反向边：对于任意一个残余网络中的某条边，它的反向边本应该是不存在的.但是对于任意一条边$(u,v)$，若当前残余网络中它的容量为$c^{\prime }(u,v)$，我们会新建一条对应的反向边$(v,u)$使得$c^{\prime }(v,u)=f(u,v)$.

反向边的作用是提供了一个反悔机制，即如果之前我们往这条边流入了$v$的流量，但现在不想要了，就可以往反向边流入$v$的流量，就相当于退流了.

有了上面的这些基本概念后，我们就可以开始具体讲解最大流的几个基于增广路的算法了.


三.Emonds-Karp算法.

这是最基础的增广路最大流算法，而且也很好写…

首先，我们已经知道了上面的一些基础知识，那么接下来就只需要一个算法流程：
1.初始设置最大流为$0$，每条边的反向边容量为$0$.
2.每一次增加流量的时候，找到一条路径上边容量均大于$0$的增广路，用bfs实现.
3.找到一条可行的增广路之后，设$v$为这条增广路上的边的最小容量，把这条增广路上所有边的容量减去$v$，这些边的反向边的增广路加上$v$.最后给最大流增加$v$.
4.重复执行步骤2,3知道找不到一条合法的增广路.

那么如何找到一条边的反向边呢？我们可以在储存邻接表的时候，让边在内存池中的下标从$2$开始，并且把一条边的反向边存在这条边下标$+1$的位置，然后就可以通过直接$xor1$来得到一条边的反向边的下标了.

Emonds-Karp算法主要代码如下：

int use[N+9],pre[N+9],flow[N+9];
queue<int>q;

void Dec_path(int st,int td,int flow){
  for (int i=td;i;i=e[pre[i]^1].y)
    e[pre[i]].f-=flow,e[pre[i]^1].f+=flow; 
}

int Bfs_path(int st,int td){
  for (int i=1;i<=n;++i) use[i]=0;
  for (;!q.empty();q.pop());
  use[st]=1;flow[st]=INF;q.push(st);
  while (!q.empty()){
    int t=q.front();q.pop();
    for (int i=lin[t];i;i=e[i].next)
      if (e[i].f&&!use[e[i].y]){
        use[e[i].y]=1;
        pre[e[i].y]=i;
        flow[e[i].y]=min(flow[t],e[i].f);
        q.push(e[i].y);
        if (e[i].y==td) {Dec_path(st,td,flow[td]);return flow[td];}
	  }
  }
  return 0;
}

int Max_flow(int st,int td){
  int res=0;
  while (2333){
    int t=Bfs_path(st,td);
    if (!t) return res;
    res+=t;
  }
}

可以证明这个算法的时间复杂度是$O(n{m}^2)$的，至于证明…等我学了再写…


四.Dinic算法.

Dinic算法是基于Emonds-Karp算法的一个改进算法.

Dinic算法的优秀之处在于，每一次在找增广路的时候，先把残余网络（将流量限制为$0$的边看做不存在）划分成一张分层图（暂时舍弃一些边），然后每次把这一张分层图上的增广路全部找完之后再重新增广，这使得这个算法可以一次找到多条增广路.

那么Dinic的算法的流程如下：
1.初始设置最大流为$0$，每条边的反向边容量为$0$.
2.用bfs将残余网络划分成一张分层图，并且同时通过是否划分到了汇点判定是否还有增广路，若没有划分到汇点说明得到最大流，直接退出.
3.在分层图上dfs，每次dfs得到源点最多可以通过这个分层图流出多少流量，并依次处理好路径上的边的流量.最后加入最大流中.
4.重复2,3直到退出.

详细讲一下dfs的实现.dfs有三个参数$dfs(k,t,f)$，表示当前点为$k$，汇点为$t$，当前点最多可以流出$f$的流量.枚举$k$可以直接到达的所有点，每次枚举到一个点$y$的时候，执行$dfs(y,t,\min (f,f(k,y)))$，并且返回得到$y$流出的流量大小$maxf$，当前的$f$减去$maxf$，并将当前函数的返回值加上$maxf$，同时处理边$(k,y)$的流量变化即可.

不过注意到这个过程中，如果一个点$k$已经被dfs了一次并且运行到了某条边$(k,y)$可用流量为$0$了，那么就直接退出.然后第二次再次dfs这个点的时候，容易发现这次运行肯定是从上次遍历到的边$(k,y)$开始的，这个时候我们就发现了一个可以优化的地方.

具体怎么优化呢？考虑记录一个数组$cur[k]$表示第$k$个点当前遍历到的边的编号，每一次dfs点$k$的时候就直接从$cur[k]$开始，并且每一次遍历完一条边后都更新$cur[k]$（可以用地址实现），这样可以做到不重复遍历任意一条边.这个优化我们称之为当前弧优化.

不过注意，当前弧优化在每一次得到分层图之后都要重新初始化一遍.

当前弧优化后的Dinic主要代码如下：

int dis[N+9];
queue<int>q;

bool Bfs_div(int st,int td){
  for (int i=1;i<=n;++i) dis[i]=0;
  dis[st]=1;q.push(st);
  while (!q.empty()){
    int t=q.front();q.pop();
    for (int i=lin[t];i;i=e[i].next)
      if (e[i].f&&!dis[e[i].y])
        dis[e[i].y]=dis[t]+1,q.push(e[i].y);
  }
  return dis[td];
}

int cur[N+9];

int Dfs_path(int k,int td,int flow){
  if (k==td) return flow;
  int res=0;
  for (int &i=cur[k];i;i=e[i].next)
    if (e[i].f&&dis[k]+1==dis[e[i].y]){
      int t=Dfs_path(e[i].y,td,min(flow,e[i].f));
      flow-=t;res+=t;
      e[i].f-=t;e[i^1].f+=t;
      if (!flow) return res;
    }
  return res;
}

int Max_flow(int st,int td){
  int res=0;
  for (;Bfs_div(st,td);res+=Dfs_path(st,td,INF))
    for (int i=1;i<=n;++i) cur[i]=lin[i];
  return res;
}

关于时间复杂度，容易发现每一次分层都是$O(n+m)$的，每一次dfs至多是$O(nm)$的，最多只会操作$n$次，所以是$O(n^{2}m)$.

其实我根本不会证明…

还有一点，用Dinic跑二分图匹配复杂度被证明了不到$O(n^{2}m)$，而是$O(m\sqrt{n})$.

你觉得我会证明吗…


五.ISAP算法.

ISAP算法是在Dinic算法的基础上进一步优化.

容易发现的是，Dinic算法每一次都要重新分层得到一张新的分层图，然而事实上每次新的分层图与旧的分层图区别并不大，这是不是显得有些浪费呢？

ISAP就是在Dinic的基础上优化了bfs分层这个过程.ISAP在每一个dfs函数的结尾，有一个重新分层的过程，也就是说，若当前节点的流还没流完，那么就修改当前点的深度.

事实上，为了方便重新分层，ISAP改变了深度的意义——从源点出发经过的点数变成了到达汇点经过的点数.这样，每一次若当前点的流没流完，就把当前点的深度$+1$，相当于让它距离源点更近，距离汇点更远.

同时，ISAP流完不是靠一次dfs出的流量为$0$来判断的了，而是看有一个深度是不是没点了，没点就说明流完了.为了实现这个过程，ISAP使用了一个$gap$数组，$gap[i]$表示深度为$i$的点数量，每次修改一个点深度的时候改变一下$gap$数组，并判定原来的深度是不是为$0$即可.

那么ISAP就讲完了，注意Dinic的当前弧优化也可以用到ISAP上.

ISAP算法的主要代码如下：

int dis[N+9],gap[N+9];
queue<int>q;

void Bfs_div(int st,int td){
  gap[dis[td]=1]=1;q.push(td);
  while (!q.empty()){
    int t=q.front();q.pop();
    for (int i=lin[t];i;i=e[i].next)
      if (!dis[e[i].y]){
        ++gap[dis[e[i].y]=dis[t]+1];
        q.push(e[i].y);
	  }
  }
}

int cur[N+9];

int Dfs_path(int k,int st,int td,int flow){
  if (k==td) return flow;
  int res=0;
  for (int &i=cur[k];i;i=e[i].next)
    if (dis[k]==dis[e[i].y]+1){
      int t=Dfs_path(e[i].y,st,td,min(flow,e[i].f));
      flow-=t;res+=t;
      e[i].f-=t;e[i^1].f+=t;
      if (!flow) return res;
	}
  if (!(--gap[dis[k]])) dis[st]=n+1;
  ++gap[++dis[k]];
  cur[k]=lin[k];
  return res;
}

int Max_flow(int st,int td){
  int res=0;
  Bfs_div(st,td);      //其实呢，ISAP把一开始bfs分层的过程去掉也是没关系的...
  for (int i=1;i<=n;++i) cur[i]=lin[i];
  for (;dis[st]<=n;res+=Dfs_path(st,st,td,INF));
  return res;
}

ISAP的时间复杂度与Dinic其实是一样的，都是$O(n^{2}m)$.不过ISAP跑二分图匹配的时候不会有$O(m\sqrt{n})$的复杂度保证，也就是说用ISAP跑二分图匹配时间复杂度也是$O(n^{2}m)$的，不过也没什么关系网络流复杂度本来就是O(玄学).

证明是不存在的，这辈子都不存在的…


六.例题与代码.

题目：luogu3376.

Emonds-Karp算法代码如下：

#include<bits/stdc++.h>
  using namespace std;
 
#define Abigail inline void
typedef long long LL;

const int N=10000,M=100000,INF=(1<<30)-1;

int n,m,st,td;
struct side{
  int y,next,f;
}e[M*2+9];
int lin[N+9],cs=1;

void Ins(int x,int y,int v){e[++cs].y=y;e[cs].f=v;e[cs].next=lin[x];lin[x]=cs;}

int use[N+9],pre[N+9],flow[N+9];
queue<int>q;

void Dec_path(int st,int td,int flow){
  for (int i=td;i;i=e[pre[i]^1].y)
    e[pre[i]].f-=flow,e[pre[i]^1].f+=flow; 
}

int Bfs_path(int st,int td){
  for (int i=1;i<=n;++i) use[i]=0;
  for (;!q.empty();q.pop());
  use[st]=1;flow[st]=INF;q.push(st);
  while (!q.empty()){
    int t=q.front();q.pop();
    for (int i=lin[t];i;i=e[i].next)
      if (e[i].f&&!use[e[i].y]){
        use[e[i].y]=1;
        pre[e[i].y]=i;
        flow[e[i].y]=min(flow[t],e[i].f);
        q.push(e[i].y);
        if (e[i].y==td) {Dec_path(st,td,flow[td]);return flow[td];}
	  }
  }
  return 0;
}

int Max_flow(int st,int td){
  int res=0;
  while (2333){
    int t=Bfs_path(st,td);
    if (!t) return res;
    res+=t;
  }
}

Abigail into(){
  scanf("%d%d%d%d",&n,&m,&st,&td);
  int x,y,v;
  for (int i=1;i<=m;++i){
  	scanf("%d%d%d",&x,&y,&v);
  	Ins(x,y,v);Ins(y,x,0);
  }
}

Abigail outo(){
  printf("%d\n",Max_flow(st,td));
}

int main(){
  into();
  outo();
  return 0;
}

Dinic算法代码如下：

#include<bits/stdc++.h>
  using namespace std;
 
#define Abigail inline void
typedef long long LL;

const int N=10000,M=100000,INF=(1<<30)-1;

int n,m,st,td;
struct side{
  int y,next,f;
}e[M*2+9];
int lin[N+9],cs=1;

void Ins(int x,int y,int v){e[++cs].y=y;e[cs].f=v;e[cs].next=lin[x];lin[x]=cs;}

int dis[N+9];
queue<int>q;

bool Bfs_div(int st,int td){
  for (int i=1;i<=n;++i) dis[i]=0;
  dis[st]=1;q.push(st);
  while (!q.empty()){
    int t=q.front();q.pop();
    for (int i=lin[t];i;i=e[i].next)
      if (e[i].f&&!dis[e[i].y])
        dis[e[i].y]=dis[t]+1,q.push(e[i].y);
  }
  return dis[td];
}

int cur[N+9];

int Dfs_path(int k,int td,int flow){
  if (k==td) return flow;
  int res=0;
  for (int &i=cur[k];i;i=e[i].next)
    if (e[i].f&&dis[k]+1==dis[e[i].y]){
      int t=Dfs_path(e[i].y,td,min(flow,e[i].f));
      flow-=t;res+=t;
      e[i].f-=t;e[i^1].f+=t;
      if (!flow) return res;
    }
  return res;
}

int Max_flow(int st,int td){
  int res=0;
  for (;Bfs_div(st,td);res+=Dfs_path(st,td,INF))
    for (int i=1;i<=n;++i) cur[i]=lin[i];
  return res;
}

Abigail into(){
  scanf("%d%d%d%d",&n,&m,&st,&td);
  int x,y,v;
  for (int i=1;i<=m;++i){
  	scanf("%d%d%d",&x,&y,&v);
  	Ins(x,y,v);Ins(y,x,0);
  }
}

Abigail outo(){
  printf("%d\n",Max_flow(st,td));
}

int main(){
  into();
  outo();
  return 0;
}

ISAP算法代码如下：

#include<bits/stdc++.h>
  using namespace std;
 
#define Abigail inline void
typedef long long LL;

const int N=10000,M=100000,INF=(1<<30)-1;

int n,m,st,td;
struct side{
  int y,next,f;
}e[M*2+9];
int lin[N+9],cs=1;

void Ins(int x,int y,int v){e[++cs].y=y;e[cs].f=v;e[cs].next=lin[x];lin[x]=cs;}

int dis[N+9],gap[N+9];
queue<int>q;

void Bfs_div(int st,int td){
  gap[dis[td]=1]=1;q.push(td);
  while (!q.empty()){
  	int t=q.front();q.pop();
    for (int i=lin[t];i;i=e[i].next)
      if (!dis[e[i].y])
        ++gap[dis[e[i].y]=dis[t]+1],q.push(e[i].y);
  }
} 

int cur[N+9];

int Dfs_path(int k,int st,int td,int flow){
  if (k==td) return flow;
  int res=0;
  for (int &i=cur[k];i;i=e[i].next)
    if (e[i].f&&dis[e[i].y]+1==dis[k]){
      int t=Dfs_path(e[i].y,st,td,min(flow,e[i].f));
      flow-=t;res+=t;
      e[i].f-=t;e[i^1].f+=t;
      if (!flow) return res;
    }
  if (!(--gap[dis[k]])) dis[st]=n+1;
  ++gap[++dis[k]];
  cur[k]=lin[k];
  return res; 
}

int Max_flow(int st,int td){
  int res=0;
  Bfs_div(st,td);
  for (int i=1;i<=n;++i) cur[i]=lin[i];
  for (;dis[st]<=n;res+=Dfs_path(st,st,td,INF));
  return res;
}

Abigail into(){
  scanf("%d%d%d%d",&n,&m,&st,&td);
  int x,y,v;
  for (int i=1;i<=m;++i){
  	scanf("%d%d%d",&x,&y,&v);
  	Ins(x,y,v);Ins(y,x,0);
  }
}

Abigail outo(){
  printf("%d\n",Max_flow(st,td));
}

int main(){
  into();
  outo();
  return 0;
}