最大流问题:如下图所示,有一个管道系统,源点是1,汇点是4,顶点之间的边表示有一条容量为权值的管道。问汇点可以接受的最大流量。这就是一个最大流问题。
流网络:假设G(V,E) 是一个有限的有向图,它的每条边(u,v)∈E都有一个非负值实数的容量c(u,v) 。如果(u,v)不属于E,我们假设c(u,v) = 0 。我们区别两个顶点:一个源点s和一个汇点t。一道网络流是一个对于所有结点u和v都有以下特性的实数函数f:V×V→R
容量限制 (Capacity Constraints):
|
f(u,v)≤c(u,v)一条边的流不能超过它的容量。
|
斜对称 (Skew Symmetry):
|
f(u,v)=-f(v,u)由
u到
v的净流必须是由
v到
u的净流的相反。(既然要看网络流,这是一定要知道的)
|
流守恒 (Flow Conservation):
|
除非
u=
s或
u=
t,否则Σ(w∈V)f(u,w)=0 一结点的净流是零,除了“制造”流的源点和“消耗”流的汇点。
|
注意f(u,v) 是由u到v的净流。如果该图代表一个实质的网络,由u到v有 4 单位的实际流及由v到u有 3 单位的实际流,那么f(u,v) = 1 及f(v,u) = − 1。
最大流的流值等于最小割(截集)的容量。
Ford-Fulkerson方法
Ford-Fulkerson方法,该方法也称作“扩充路径方法”,该方法是大量算法的基础,有多种实现方法。Ford-Fulkerson算法在实际中并不常用,但是它提供了一种思想:先找到一条从源点到汇点的增广路径,这条路径的容量是其中容量最小的边的容量。然后,通过不断找增广路,一步步扩大流量,当找不到增广路时,就得到最大流了(最大流最小割定理)。可以看看Ford-Fulkerson算法的伪代码,
FORD-FULKERSON(G,s,t)
1 for each edge(u,v)∈ E[G]
2 do f[u,v]=0
3. f[v,u]=0
4. while there exists a path from s to t in the residual network Gf
5. do cf(p) ← min {cf(u, v) : (u, v) is in p}
6. for each edge (u, v) in p
7. do f[u, v] ← f[u, v] + cf(p)
8. f[v, u] ← -f[u, v]
网络流算法要基于三种思想:残留网络(Residual Network),增广路径(Augmenting Path)和割(Cut)。这里不再一一说明。
从伪代码可以看出,Ford-Fulkerson的核心过程就是:第4~8行的 while 循环反复找出 Gf 中的增广路径 p,并把沿 p 的流 f 加上其残留容量cf(p)。当不再有增广路径时,流 f 就是一个最大流。而之所以称其为方法而不是算法,就是由于找到增广路径的方法不同,其复杂度也不同,它只是提供一种解出最大流的一种思想而不是直接求出最大流。
Edmonds-Karp(EK)算法
EK算法基于Ford-Fulkerson算法,唯一的区别是将第 4 行用BFS(广度优先搜索)来实现对增广路径 p 的计算。EK算法伪代码基本和上边的Ford-Fulkerson算法一样。类似用DFS实现的还有Dinic算法。它们都属于SAP(Shortest Augmenting Path)算法,从英文即可看出,它们每次都在寻找最短增广路。对于EK算法,每次用一遍 BFS 寻找从源点 s 到终点 t 的最短路作为增广路径,然后增广流量 f 并修改残量网络,直到不存在新的增广路径。E-K 算法的时间复杂度为 O(VE^2),适用于稀疏边,由于 BFS 要搜索全部小于最短距离的分支路径之后才能找到终点,因此频繁的 BFS 效率是比较低的。
代码模板如下:
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <queue>
#include <algorithm>
#include <iostream>
using namespace std;
int n,m; //n个顶点,m条边
int c[205][205]; //c[i][j]=k表示从i到j有容量为k的边
int pre[105]; //记录增广路径
bool bfs(int s,int t) //寻找增广路径
{
queue<int> q;
bool v[205];
memset(v,0,sizeof(v));
q.push(s);
v[s]=1; pre[s]=s;
while (!q.empty())
{
int w=q.front();
q.pop();
for (int i=1;i<=n;i++)
{
if ((c[w][i])&&(!v[i]))
{
pre[i]=w;
v[i]=1;
if (i==t) return 1; //找到增广路径
q.push(i);
}
}
}
return 0;
}
int main()
{
while (scanf("%d%d",&m,&n)!=EOF)
{
int ans=0;
memset(c,0,sizeof(c));
for (int i=0;i<m;i++)
{
int j,k,l;
scanf("%d%d%d",&j,&k,&l);
c[j][k]+=l; //有重边的话就把容量增加
}
while (bfs(1,n))
{
int d=c[pre[n]][n];
for (int i=n;i!=1;i=pre[i])
{
d=min(d,c[pre[i]][i]); //找到可以增加的流
}
for (int i=n;i!=1;i=pre[i])
{
c[pre[i]][i]-=d;
c[i][pre[i]]+=d; //建立反向边
}
ans+=d;
}
printf("%d\n",ans);
}
return 0;
}
这里有个建立反向边的思想,简单的来讲反向边就是可以给程序一个后悔的机会。我在网上看到一个解释不错:http://www.wutianqi.com/?p=3107,我把它先粘在这里。
下面我们来考虑如何求最大流。
首先来看一下基本的网络流最大流模型。
有n个点,有m条有向边,有一个点很特殊,只出不进,叫做源点,通常规定为1号点。另一个点也很特殊,只进不出,叫做汇点,通常规定为n号点。每条有向边上有两个量,容量和流量,从i到j的容量通常用c[I,j]表示,流量则通常是f[I,j]。通常可以把这些边想象成道路,流量就是这条道路的车流量,容量就是道路可承受的最大的车流量。很显然的,流量<=容量。而对于每个不是源点和汇点的点来说,可以类比的想象成没有存储功能的货物的中转站,所有”进入”他们的流量和等于所有从他本身”出去”的流量。
把源点比作工厂的话,问题就是求从工厂最大可以发出多少货物,是不至于超过道路的容量限制,也就是,最大流。
比如这个图。每条边旁边的数字表示它的容量。
首先,假如所有边上的流量都没有超过容量(不大于容量),那么就把这一组流量,或者说,这个流,称为一个可行流。一个最简单的例子就是,零流,即所有的流量都是0的流。
我们就从这个零流开始考虑,假如有这么一条路,这条路从源点开始一直一段一段的连到了汇点,并且,这条路上的每一段都满足流量<容量,注意,是严格的<,而不是<=。那么,我们一定能找到这条路上的每一段的(容量-流量)的值当中的最小值 delta。我们把这条路上每一段的流量都加上这个delta,一定可以保证这个流依然是可行流,这是显然的。
这样我们就得到了一个更大的流,他的流量是之前的流量+delta,而这条路就叫做增广路。
我们不断地从起点开始寻找增广路,每次都对其进行增广,直到源点和汇点不连通,也就是找不到增广路为止。当找不到增广路的时候,当前的流量就是最大流,这个结论非常重要。
寻找增广路的时候我们可以简单的从源点开始做bfs,并不断修改这条路上的delta量,直到找到源点或者找不到增广路。
这里要先补充一点,在程序实现的时候,我们通常只是用一个c数组来记录容量,而不记录流量,当流量+1的时候,我们可以通过容量-1来实现,以方便程序的实现。
先来看看BFS部分的代码(C/C++实现):
// 用BFS来判断从结点s到t的路径上是否还有delta
// 即判断s,t之间是否还有增广路径,若有,返回1
bool BFS(int s, int t)
{
queue<int> que;
memset(pre, -1, sizeof(pre));
memset(vis, false, sizeof(vis));
pre[s] = s;
vis[s] = true;
que.push(s);
int p;
while(!que.empty())
{
p = que.front();
que.pop();
for(int i=1; i<=M; ++i)
{
if(r[p][i]>0 && !vis[i])
{
pre[i] = p;
vis[i] = true;
if(i == t) // 存在增广路径
return true;
que.push(i);
}
}
}
return false;
}
但事实上并没有这么简单,上面所说的增广路还不完整,比如说下面这个网络流模型。
如果顶点多的话,往往N^2存不下,这时候就要存边:
存每条边的出发点,终止点和价值,然后排序一下,再记录每个出发点的位置。以后要调用从出发点出发的边时候,只需要从记录的位置开始找即可(其实可以用链表)。优点是时间加快空间节省,缺点是编程复杂度将变大,所以在题目允许的情况下,建议使用邻接矩阵。
如果一次重标号时,出现距离断层,则可以证明ST无可行流,此时则可以直接退出算法。
为了使每次找增广路的时间变成均摊O(V),还有一个重要的优化是对于每个点保存“当前弧”:初始时当前弧是邻接表的第一条弧;在邻接表中查找时从当前弧开始查找,找到了一条允许弧,就把这条弧设为当前弧;改变距离标号时,把当前弧重新设为邻接表的第一条弧。
SAP简化的描述是:程序开始时用一个反向 BFS 初始化所有顶点的距离标号,之后从源点开始,进行如下三种操作:(1)当前顶点 i 为终点时增广 (2) 当前顶点有满足 dist[i] = dist[j] + 1 的出弧时前进 (3) 当前顶点无满足条件的出弧时重标号并回退一步。整个循环当源点 s 的距离标号 dist[s] >= n 时结束。对 i 点的重标号操作可概括为 dist[i] = 1 + min{dist[j] : (i,j)属于残量网络Gf}。
实际应用中可以不用bfs初始化dist[i]而是把所有dist[i]设为0,这样可以减少编程复杂度同时并不影响时间复杂度。
代码模板如下:
我们第一次找到了1-2-3-4这条增广路,这条路上的delta值显然是1。于是我们修改后得到了下面这个流。(图中的数字是容量)
这时候(1,2)和(3,4)边上的流量都等于容量了,我们再也找不到其他的增广路了,当前的流量是1。
但这个答案明显不是最大流,因为我们可以同时走1-2-4和1-3-4,这样可以得到流量为2的流。
那么我们刚刚的算法问题在哪里呢?问题就在于我们没有给程序一个”后悔”的机会,应该有一个不走(2-3-4)而改走(2-4)的机制。那么如何解决这个问题呢?回溯搜索吗?那么我们的效率就上升到指数级了。
而这个算法神奇的利用了一个叫做反向边的概念来解决这个问题。即每条边(I,j)都有一条反向边(j,i),反向边也同样有它的容量。
我们直接来看它是如何解决的:
在第一次找到增广路之后,在把路上每一段的容量减少delta的同时,也把每一段上的反方向的容量增加delta。即在Dec(c[x,y],delta)的同时,inc(c[y,x],delta)
我们来看刚才的例子,在找到1-2-3-4这条增广路之后,把容量修改成如下
这时再找增广路的时候,就会找到1-3-2-4这条可增广量,即delta值为1的可增广路。将这条路增广之后,得到了最大流2。
那么,这么做为什么会是对的呢?我来通俗的解释一下吧。
事实上,当我们第二次的增广路走3-2这条反向边的时候,就相当于把2-3这条正向边已经是用了的流量给”退”了回去,不走2-3这条路,而改走从2点出发的其他的路也就是2-4。(有人问如果这里没有2-4怎么办,这时假如没有2-4这条路的话,最终这条增广路也不会存在,因为他根本不能走到汇点)同时本来在3-4上的流量由1-3-4这条路来”接管”。而最终2-3这条路正向流量1,反向流量1,等于没有流量。
这就是这个算法的精华部分,利用反向边,使程序有了一个后悔和改正的机会。而这个算法和我刚才给出的代码相比只多了一句话而已。
Improved SAP(ISAP)算法
SAP算法(by dd_engi):求最大流有一种经典的算法,就是每次找增广路时用BFS找,保证找到的增广路是弧数最少的,也就是所谓的 Edmonds-Karp 算法。可以证明的是在使用最短路增广时增广过程不超过 V * E次,每次 BFS 的时间都是O(E),所以 Edmonds-Karp 的时间复杂度就是O(V * E^2)。
如果能让每次寻找增广路时的时间复杂度降下来,那么就能提高算法效率了,使用距离标号的最短增广路算法就是这样的。所谓距离标号,就是某个点到汇点的最少的弧的数量(另外一种距离标号是从源点到该点的最少的弧的数量,本质上没什么区别)。设点 i 的标号为D[i],那么如果将满足D[i] = D[j] + 1的弧(i,j))叫做允许弧,且增广时只走允许弧,那么就可以达到“怎么走都是最短路”的效果。每个点的初始标号可以在一开始用一次从汇点沿所有反向边的BFS求出,实践中可以初始设全部点的距离标号为0,问题就是如何在增广过程中维护这个距离标号。
维护距离标号的方法是这样的:当找增广路过程中发现某点出发没有允许弧时,将这个点的距离标号设为由它出发的所有弧的终点的距离标号的最小值加一。这种维护距离标号的方法的正确性我就不证了。由于距离标号的存在,由于“怎么走都是最短路”,所以就可以采用DFS找增广路,用一个栈保存当前路径的弧即可。当某个点的距离标号被改变时,栈中指向它的那条弧肯定已经不是允许弧了,所以就让它出栈,并继续用栈顶的弧的端点增广。为了使每次找增广路的时间变成均摊O(V),还有一个重要的优化是对于每个点保存“当前弧”:初始时当前弧是邻接表的第一条弧;在邻接表中查找时从当前弧开始查找,找到了一条允许弧,就把这条弧设为当前弧;改变距离标号时,把当前弧重新设为邻接表的第一条弧,还有一种在常数上有所优化的写法是改变距离标号时把当前弧设为那条提供了最小标号的弧。当前弧的写法之所以正确就在于任何时候我们都能保证在邻接表中当前弧的前面肯定不存在允许弧。
还有一个常数优化是在每次找到路径并增广完毕之后不要将路径中所有的顶点退栈,而是只将瓶颈边以及之后的边退栈,这是借鉴了Dinic算法的思想。注意任何时候待增广的“当前点”都应该是栈顶的点的终点。这的确只是一个常数优化,由于当前边结构的存在,我们肯定可以在O(n)的时间内复原路径中瓶颈边之前的所有边。
优化:
1.邻接表优化:如果顶点多的话,往往N^2存不下,这时候就要存边:存每条边的出发点,终止点和价值,然后排序一下,再记录每个出发点的位置。以后要调用从出发点出发的边时候,只需要从记录的位置开始找即可(其实可以用链表)。优点是时间加快空间节省,缺点是编程复杂度将变大,所以在题目允许的情况下,建议使用邻接矩阵。
2.GAP优化:如果一次重标号时,出现距离断层,则可以证明ST无可行流,此时则可以直接退出算法/
3.当前弧优化:为了使每次找增广路的时间变成均摊O(V),还有一个重要的优化是对于每个点保存“当前弧”:初始时当前弧是邻接表的第一条弧;在邻接表中查找时从当前弧开始查找,找到了一条允许弧,就把这条弧设为当前弧;改变距离标号时,把当前弧重新设为邻接表的第一条弧。
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <queue>
#include <algorithm>
#include <iostream>
using namespace std;
int n,m; //n个顶点,m条边
int c[205][205]; //c[i][j]=k表示从i到j有容量为k的边
int pre[205]; //记录增广路径
int cur[205]; //当前弧优化
int gap[205]; //gap优化
int d[205]; //距离标号
/*void bfs(int s,int t) //bfs预处理d[i],可以去掉
{
queue <int> q;
q.push(t);
while (!q.empty())
{
int x=q.front();
q.pop();
for (int i=1;i<n;i++)
{
if ((!d[i])&&(c[i][x]))
{
d[i]=d[x]+1;
q.push(i);
}
}
}
}*/
int augment(int s,int t) //增广
{
int x=t,m=INT_MAX;
while (x!=s)
{
m=min(m,c[pre[x]][x]);
x=pre[x];
}
x=t;
while (x!=s)
{
c[pre[x]][x]-=m;
c[x][pre[x]]+=m;
x=pre[x];
}
return m;
}
int Isap(int s,int t)
{
int ans=0;
memset(d,0,sizeof(d));
memset(gap,0,sizeof(gap));
//bfs(1,n);
gap[0]=n; //将所有d[i]设为0
for (int i=1;i<n;i++) { cur[i]=1;}
pre[s]=s;
int x=s;
while (d[s]<n)
{
if (x==t)
{
ans+=augment(1,n);
x=s;
}
bool find=0;
for (int i=cur[x];i<=n;i++)
{
if ((c[x][i])&&(d[x]==d[i]+1))
{
find=1;
pre[i]=x;
cur[x]=i;
x=i;
break;
}
}
if (!find)
{
if (--gap[d[x]]==0) break; //gap优化
int m=n; //retreat
for (int i=1;i<=n;i++)
if (c[x][i]) m=min(m,d[i]);
d[x]=m+1;
gap[d[x]]++;
cur[x]=1;
if (x!=s) x=pre[x];
}
}
return ans;
}
int main()
{
while (scanf("%d%d",&m,&n)!=EOF)
{
memset(c,0,sizeof(c));
for (int i=0;i<m;i++)
{
int j,k,l;
scanf("%d%d%d",&j,&k,&l);
c[j][k]+=l; //有重边的话就把容量增加
}
printf("%d\n",Isap(1,n));
}
return 0;
}