通过 USACO 4.2.1 Ditch 学习一下最大流算法 。可惜它给的测试数据几乎没有任何杀伤力,后面测试时我们采用 DD_engi 写的程序生成的加强版数据。
总体上来说,最大流算法分为两大类:增广路 (Augmenting Path) 和预流推进重标号 (Push Relabel) 。也有算法同时借鉴了两者的长处,如 Improved SAP 。本篇主要介绍增广路类算法,思想、复杂度及实际运行效率比较,并试图从中选择一种兼顾代码复杂度和运行效率的较好方案。以下我们将会看到,有时理论分析的时间复杂度并不能很好的反映一种算法的实际效率。
所有增广路算法的基础都是 Ford - Fulkerson 方法。称之为方法而不是算法是因为 Ford - Fulkerson 只提供了一类思想,在此之上的具体操作可有不同的实现方案。
给定一个有向网络 G(V,E) 以及源点 s 终点 t ,FF 方法描述如下:
假设有向网络 G 中边 (i,j) 的容量为 c(i,j) ,当前流量为 f(i,j) ,则此边的剩余流量即为 r(i,j) = c(i,j) - f(i,j) ,其反向边的剩余流量为 r(j,i) = f(i,j) 。有向网中所有剩余流量 r(i,j) > 0 的边构成残量网络 Gf ,增广路径p即是残量网络中从源点 s 到终点 t 的路径。
沿路径 p 增广流量 f 的操作基本都是相同的,各算法的区别就在于寻找增广路径 p 的方法不同。例如可以寻找从 s 到 t 的最短路径,或者流量最大的路径。
Shortest Augmenting Path (SAP) 是每次寻找最短增广路的一类算法,Edmonds - Karp 算法以及后来著名的 Dinic 算法都属于此。SAP 类算法可统一描述如下:
在无权边的有向图中寻找最短路,最简单的方法就是广度优先搜索 (BFS),E-K 算法就直接来源于此。每次用一遍 BFS 寻找从源点 s 到终点 t 的最短路作为增广路径,然后增广流量 f 并修改残量网络,直到不存在新的增广路径。
E-K 算法的时间复杂度为 O(VE2),由于 BFS 要搜索全部小于最短距离的分支路径之后才能找到终点,因此可以想象频繁的 BFS 效率是比较低的。实践中此算法使用的机会较少。
BFS 寻找终点太慢,而 DFS 又不能保证找到最短路径。1970年 Dinic 提出一种思想,结合了 BFS 与 DFS 的优势,采用构造分层网络的方法可以较快找到最短增广路,此算法又称为阻塞流算法 (Blocking Flow Algorithm)。
首先定义分层网络 AN(f)。在残量网络中从源点 s 起始进行 BFS,这样每个顶点在 BFS 树中会得到一个距源点 s 的距离 d,如 d(s) = 0,直接从 s 出发可到达的点距离为 1,下一层距离为2 ... 。称所有具有相同距离的顶点位于同一层,在分层网络中,只保留满足条件 d(i) + 1 = d(j) 的边,这样在分层网络中的任意路径就成为到达此顶点的最短路径。
Dinic 算法每次用一遍 BFS 构建分层网络 AN(f),然后在 AN(f) 中一遍 DFS 找到所有到终点 t 的路径增广;之后重新构造 AN(f),若终点 t 不在 AN(f) 中则算法结束。DFS 部分算法可描述如下:
实际代码中不必真的用一个图来存储分层网络,只需保存每个顶点的距离标号并在 DFS 时判断 dist[i] + 1 = dist[j] 即可。Dinic 的时间复杂度为 O(V2E)。由于较少的代码量和不错的运行效率,Dinic 在实践中比较常用。具体代码可参考 DD_engi 网络流算法评测包中的标程,这几天 dinic 算法的实现一共看过和比较过将近 10 个版本了,DD 写的那个在效率上是数一数二的,逻辑上也比较清晰。
本次介绍的重头戏。通常的 SAP 类算法在寻找增广路时总要先进行 BFS,BFS 的最坏情况下复杂度为 O(E),这样使得普通 SAP 类算法最坏情况下时间复杂度达到了 O(VE2)。为了避免这种情况,Ahuja 和 Orlin 在1987年提出了Improved SAP 算法,它充分利用了距离标号的作用,每次发现顶点无出弧时不是像 Dinic 算法那样到最后进行 BFS,而是就地对顶点距离重标号,这样相当于在遍历的同时顺便构建了新的分层网络,每轮寻找之间不必再插入全图的 BFS 操作,极大提高了运行效率。国内一般把这个算法称为 SAP...显然这是不准确的,毕竟从字面意思上来看 E-K 和 Dinic 都属于 SAP,我还是习惯称为 ISAP 或改进的 SAP 算法。
与 Dinic 算法不同,ISAP 中的距离标号是每个顶点到达终点 t 的距离。同样也不需显式构造分层网络,只要保存每个顶点的距离标号即可。程序开始时用一个反向 BFS 初始化所有顶点的距离标号,之后从源点开始,进行如下三种操作:(1)当前顶点 i 为终点时增广 (2) 当前顶点有满足 dist[i] = dist[j] + 1 的出弧时前进 (3) 当前顶点无满足条件的出弧时重标号并回退一步。整个循环当源点 s 的距离标号 dist[s] >= n 时结束。对 i 点的重标号操作可概括为 dist[i] = 1 + min{dist[j] : (i,j)属于残量网络Gf}。具体算法描述如下:
算法中的允许弧是指在残量网络中满足 dist[i] = dist[j] + 1 的弧。Retreat 过程中若从 i 出发没有弧属于残量网络 Gf 则把顶点距离重标号为 n 。
虽然 ISAP 算法时间复杂度与 Dinic 相同都是 O(V2E),但在实际表现中要好得多。要提的一点是关于 ISAP 的一个所谓 GAP 优化。由于从 s 到 t 的一条最短路径的顶点距离标号单调递减,且相邻顶点标号差严格等于1,因此可以预见如果在当前网络中距离标号为 k (0 <= k < n) 的顶点数为 0,那么可以知道一定不存在一条从 s 到 t 的增广路径,此时可直接跳出主循环。在我的实测中,这个优化是绝对不能少的,一方面可以提高速度,另外可增强 ISAP 算法时间上的稳定性,不然某些情况下 ISAP 会出奇的费时,而且大大慢于 Dinic 算法。相对的,初始的一遍 BFS 却是可有可无,因为 ISAP 可在循环中自动建立起分层网络。实测加不加 BFS 运行时间差只有 5% 左右,代码量可节省 15~20 行。
1972年还是那个 E-K 组合提出的另一种最大流算法。每次寻找增广路径时不找最短路径,而找容量最大的。可以预见,此方法与 SAP 类算法相比可更快逼近最大流,从而降低增广操作的次数。实际算法也很简单,只用把前面 E-K 算法的 BFS 部分替换为一个类 Dijkstra 算法即可。USACO 4.2 节的说明详细介绍了此算法,这里就不详述了。
时间复杂度方面。BFS 是 O(E),简单 Dijkstra 是 O(V2),因此效果可想而知。但提到 Dijkstra 就不能不提那个 Heap 优化,虽然 USACO 的算法例子中没有用 Heap ,我自己还是实现了一个加 Heap 的版本,毕竟 STL 的优先队列太好用了不加白不加啊。效果也是非常明显的,但比起 Dinic 或 ISAP 仍然存在海量差距,这里就不再详细介绍了。
不知道怎么翻比较好,索性就这么放着吧。叫什么的都有,容量缩放算法、容量变尺度算法等,反正就那个意思。类似于二分查找的思想,寻找增广路时不必非要局限于寻找最大容量,而是找到一个可接受的较大值即可,一方面有效降低寻找增广路时的复杂度,另一方面增广操作次数也不会增加太多。时间复杂度 O(E2logU) 实际效率嘛大约稍好于最前面 BFS 的 E-K 算法,稀疏图时表现较优,但仍然不敌 Dinic 与 ISAP。
重头戏之二,虽然引用比较多,哎~
首先引用此篇强文 《Maximum Flow: Augmenting Path Algorithms Comparison》
对以上算法在稀疏图、中等稠密图及稠密图上分别进行了对比测试。直接看结果吧:
稀疏图:
ISAP 轻松拿下第一的位置,图中最左边的 SAP 应该指的是 E-K 算法,这里没有比较 Dinic 算法是个小遗憾吧,他把 Dinic 与 SAP 归为一类了。最大流量路径的简单 Dijkstra 实现实在是太失败了 - -,好在 Heap 优化后还比较能接受……可以看到 Scaling 算法也有不错的表现。
中等稠密图:
ISAP 依然领先。最大流量算法依然不太好过……几个 Scaling 类算法仍然可接受。
稠密图:
ISAP……你无敌了!这次可以看出 BFS 的 Scaling 比 DFS 实现好得多,而且几乎与 Improved Scaling 不相上下,代码量却差不少。看来除 ISAP 外 BFS Scaling 也是个不错的选择。
最后补个我自己实测的图,比较算法有很多是 DD 网络流算法评测包里的标程,评测系统用的 Cena,评测数据为 DD ditch 数据生成程序生成的加强版数据:
我一共写了 7 个版本的 ISAP ,Dinic 算法也写了几个递归版的但效率都不高,只放上来一个算了。从上图来看似乎 ISAP 全面超越了号称最大流最快速算法的 HLPP,但在另外一台机器上测试结果与此却不大相同,有时 ISAP 与 HLPP 基本持平,有时又稍稍慢一些。在这种差距非常小的情况下似乎编译器的效果也比较明显。那个 HLPP 是用 PASCAL 写的,我不知在 Win32 平台下编译代码效率如何,至少我的几个 ISAP 用 VC2008 + SP1 编译比用 g++ 要强不少,也可能是参数设置的问题。
不过这些都是小事,关键是见证了 ISAP 的实际效率。从上面可以看出不加 GAP 优化的 ISAP 有几个测试点干脆无法通过,而不加 BFS 却无甚大影响,递归与非递归相差在 7% 左右的样子。综合以上表现,推荐采用 ISAP 不加 BFS,非递归 + GAP 优化的写法,Ditch 这道题一共也才 80 行左右代码。要提的一点是 GAP 优化用递归来表现的话不如 while 循环来得直接。另外,看起来 HLPP 表现确实很优秀,有机会也好好研究一下吧,预流推进重标号算法也是一大类呢……
最后附上一个 ISAP + GAP + BFS 的非递归版本代码,网络采用邻接表 + 反向弧指针:
2009年4月13日 02:55
兄弟 你这篇文章是我看到的网上关于最大流的最好的文章了 我受益匪浅。。。
我找了很久 大多数都是垃圾
2009年5月08日 05:09
写得不错,赞一个~
2009年5月24日 04:29
终于找到了,好东西要顶!
2009年7月25日 21:07
基本上是Network Flows: Thoery, Algorithms, and Applications前几章的小结。不过写的的确不错!
感兴趣的再去看看原书英文版吧!
2009年10月01日 18:35
为什么在pku的oj上提交时re呢
2009年10月02日 07:03
@ally: 嗯,确实没有仔细检查呢
这段代码本来是写给 USACO 上的,只是为反映算法而已,在 POJ 上变成了多组测试数据,而两组数据间少了初始化代码,所以最后就产生了链表终点没有 NULL 指针终结的错误;而且多了几组很阴险的数据 ,比如自己到自己的边啦等需要多检查一下。
在外层 while 循环中添加对 Net dist numbs 三组的 memset
读入数据后添加一个 if(u == v)continue;
即可
2009年10月02日 17:38
谢谢!我本来是用pascal的,现在改用c,所以想找些程序看看,你写的不错哦。。。。
2009年10月10日 00:27
excess scaling....
O(n*m+n^2*log(U))
2009年10月10日 09:32
@DarkRaven: 多谢指正!写这篇时只是为了做题在网上临时找了些资料(貌似写的比较清楚的很少)和自己的一点试验总结,并没有找到论文原文,所以理论错误应该是有不少吧。不过本来就不是很正式的介绍,而且也不想再修改,希望读到的只看个大概不要以此为准就好 = =
研究的话 Ahuja 和 Orlin 的文章就有必要仔细看下了 呵呵
2010年1月05日 20:34
文章写得太好了 DD 网络流算法评测包能否发给我一份?谢谢 wangchaohui@hotmail.com
2010年4月15日 03:02
请问下,PKU3281(http://acm.pku.edu.cn/JudgeOnline/problem?id=3281)能不能用isap写?
2010年4月15日 03:09
大牛帮忙看看这个模板错哪了,谢谢了!
#include <iostream>
#include <queue>
#define msize 1024 //最大顶点数目
#define inf 0x7fffffff
using namespace std;
int d[msize]; //标号
int r[msize][msize]; //残留网络,初始为原图
int num[msize]; //num[i]表示标号为i的顶点数有多少
int pre[msize];
int n,m,s,t; //n个顶点,m条边,从源点s到汇点t
void ini_d() //BFS计算标号,汇点t标号为0
{
int k;
queue<int>Q;
memset(d,-1,sizeof(d));
memset(num,0,sizeof(num));
Q.push(t);
d[t]=0;
num[0]=1;
while (!Q.empty())
{
k=Q.front();
Q.pop();
for (int i=0;i<n;i++)
{
if (d[i]==-1&&r[i][k]>0)
{
d[i]=d[k]+1;
Q.push(i);
num[d[i]]++;
}
}
}
}
int findAlowArc(int i) //从i出发寻找允许弧
{
int j;
for (j=0;j<n;j++)
if (r[i][j]>0&&d[i]==d[j]+1) return j;
return -1;
}
int reLable(int i) //重新标号
{
int mm=inf;
for (int j=0;j<n;j++)
if (r[i][j]>0) mm=mm<d[j]+1? mm : d[j]+1;
return mm==inf ? n : mm;
}
int maxFlow(int s,int t) //从源点s出发的最大流
{
ini_d();
int flow=0,i=s,j,k;
int delta; //增量
pre[s]=-1;
while (d[s]<n)
{
j=findAlowArc(i);
if (j>=0)
{
pre[j]=i;
i=j;
if (i==t) //更新残留网络
{
delta=inf;
for (k=t;k!=s;k=pre[k])
delta = delta<r[pre[k]][k] ? delta : r[pre[k]][k];
for (k=t;k!=s;k=pre[k])
{
r[pre[k]][k] -= delta;
r[k][pre[k]] += delta;
}
flow += delta;
i=s;
}
}
else
{
int x=reLable(i); //重新标号
num[x]++;
num[d[i]]--;
if (num[d[i]]==0) return flow; //间隙优化
d[i]=x;
if (i!=s) i=pre[i];
}
}
return flow;
}
我PKU3281的代码是这样的:
#include <iostream>
#include <queue>
#define msize 500 //最大顶点数目
#define inf 0x7fffffff
using namespace std;
int d[msize]; //标号
int r[msize][msize]; //残留网络,初始为原图
int num[msize]; //num[i]表示标号为i的顶点数有多少
int pre[msize];
int n,m,s,t; //n个顶点,m条边,从源点s到汇点t
void ini_d() //BFS计算标号,汇点t标号为0
{
int k;
queue<int>Q;
memset(d,-1,sizeof(d));
memset(num,0,sizeof(num));
Q.push(t);
d[t]=0;
num[0]=1;
while (!Q.empty())
{
k=Q.front();
Q.pop();
for (int i=0;i<n;i++)
{
if (d[i]==-1&&r[i][k]>0)
{
d[i]=d[k]+1;
Q.push(i);
num[d[i]]++;
}
}
}
}
int findAlowArc(int i) //从i出发寻找允许弧
{
int j;
for (j=0;j<n;j++)
if (r[i][j]>0&&d[i]==d[j]+1) return j;
return -1;
}
int reLable(int i) //重新标号
{
int mm=inf;
for (int j=0;j<n;j++)
if (r[i][j]>0) mm=mm<d[j]+1? mm : d[j]+1;
return mm==inf ? n : mm;
}
int maxFlow() //从源点s出发的最大流
{
ini_d();
int flow=0,i=s,j,k;
int delta; //增量
pre[s]=-1;
while (d[s]<n)
{
j=findAlowArc(i);
if (j>=0)
{
pre[j]=i;
i=j;
if (i==t) //更新残留网络
{
delta=inf;
for (k=t;k!=s;k=pre[k])
delta = delta<r[pre[k]][k] ? delta : r[pre[k]][k];
for (k=t;k!=s;k=pre[k])
{
r[pre[k]][k] -= delta;
r[k][pre[k]] += delta;
}
flow += delta;
i=s;
}
}
else
{
int x=reLable(i); //重新标号
num[x]++;
num[d[i]]--;
if (num[d[i]]==0) return flow; //间隙优化
d[i]=x;
if (i!=s) i=pre[i];
}
}
return flow;
}
int main()
{
int D,F,N,ff,dd,i,j,k;
while(scanf("%d%d%d",&N,&F,&D)!=EOF)
{
s=0;
t=F+N+N+D+1;
n=t+1;
memset(r,0,sizeof(r));
for(i=1;i<=N;i++)
{
scanf("%d%d",&ff,&dd);
for(j=0;j<ff;j++)
{
scanf("%d",&k);
r[k][i+F]=1;
}
for(j=1;j<=dd;j++)
{
scanf("%d",&k);
r[F+N+i][F+N+N+k]=1;
}
}
for(i=1;i<=F;i++)
r[s][i]=1;
for(i=1;i<=D;i++)
r[F+N+N+i][t]=1;
for(i=1;i<=N;i++)
r[i+F][i+F+N]=1;
printf("%d\n",maxFlow());
}
return 0;
}
下面是一组数据:
30 40 60
6 6 7 23 24 29 33 40 7 8 20 25 41 48
3 1 19 25 28 58
5 5 11 12 19 24 36 7 10 14 17 24
4 6 9 13 32 40 3 5 29 32 41 45
3 2 18 20 31 28 52
3 5 10 17 35 12 19 40 53 55
3 8 2 8 31 3 8 9 18 21 42 47 55
5 7 8 11 14 29 40 4 7 13 17 28 45 49
5 2 11 18 28 32 35 27 29
3 8 7 38 40 4 11 17 26 28 30 38 55
3 2 5 9 25 12 33
4 5 2 13 23 24 3 12 34 52 59
4 3 1 21 30 35 6 23 40
2 9 30 31 18 23 27 31 34 38 51 57 58
6 6 8 13 18 28 30 32 1 30 41 50 57 59
4 5 5 19 28 32 38 42 43 44 59
3 3 4 21 32 13 16 33
4 6 12 15 20 28 5 6 21 38 43 57
2 2 6 17 36 48
3 3 29 32 35 2 6 54
1 4 18 11 19 48 55
2 3 2 35 8 10 57
2 10 4 29 3 4 5 13 17 20 31 40 41 59
3 3 7 12 17 35 45 57
3 2 20 21 24 2 6
5 5 5 16 18 29 35 25 30 36 39 41
4 5 5 27 33 35 33 35 45 49 55
0 6 7 18 29 34 46 60
3 11 27 31 39 11 14 23 25 29 32 36 40 44 48 56
1 5 21 5 17 25 47 49
output: 29
但我的输出时 26
2010年4月19日 00:49
@hi_boy2010: 额,应该是你的模板有问题。这道题我照搬你的 main() 函数,只是把最大流部分换成自己的程序就 AC 了。(16ms,为了适应性用的 vector 可能效率有影响)
你的模板在 1273 的 Ditch 也是 WA 好像
我一般都用邻接表存图,用矩阵的话存在一个问题,如果两点间双向都有一定流量限制的边,用矩阵就没法处理,而且效率也低一些。3281这题应该不用考虑这问题。你的模板我大致看了下也没发现什么,可能有比较隐蔽的bug吧。。。
2010年4月19日 02:02
嗯,谢谢指点。我找到错了。memset(d,-1,sizeof(d));d应该初始化为>=n 。相应的将if (d[i]==-1&&r[i][k]>0)改下就可以了。但就是不知道是什么原因....
2010年7月18日 19:12
我怎么利用此模板实现计算一个点到多个点的最大流量?不是多个收点,而是计算一个点到其余多个点的最大流量,有点类似计算一个点到多个点的最短路。如果每次调用此程序,文件读入的时间太大了。请指教下
2010年7月18日 19:14
能否将文件读入写成个函数,只调用一次,后面直接调用图结构就可以了,怎么做呢,我是新手,方便的话请email联系
2010年8月02日 00:45
@旷野蓝风: 写成函数的话当然可以啊,读文件就是 main() 函数里 while 部分。具体怎么读得看数据格式,然后计算之前把建好的图先复制一下,因为计算的时候会改变原图数据。代码很久没看了,要写得先复习一下呵呵。email 怎么联系?
2010年11月07日 04:23
你好,我是一个算法初学者,苦于一个问题:测试数据!你也知道OJ上给的例子都不是很庞大的数据,不能AC都不知道自己错在哪里,所以看你文章开头说的DD_engi数据?能说的详细一些吗?想问一下你平时的测试数据怎么得到的?
2010年11月08日 00:16
@huihui: 抱歉,隔的时间有点长,那个DD_engi评测包已经忘了丢到哪了(大概是随着上一块硬盘一起废掉了)。思路是写一段程序,随机生成一个比较大的图(具体规模、稀疏或稠密或有哪些特点可以用不同的方法来实现)作为输入数据,保存在一个文本文件里如in.txt,然后用一个标准程序的输出out.txt作为检测标准,然后依次执行各个算法的标准代码,判断结果是否一致同时统计执行时间并为各参测程序打分排序。后面这部分工作可用评测系统CENA完成。当然前面数据也可以生成多组输入一并检测
我并不是专业参加算法竞赛的,只是平时有兴趣自己做一些题目而已,所以某些经验可能不是很有代表性。我的测试数据,一开始是自己思考各种一般或极端情况构造数据,然后到讨论区看别人一些有代表性的测试数据,偶尔也会写程序产生多组一起检测。这点USACO比较好测试数据都能看到,但那边题目比较少。如果上面测试都通过仍然AC不掉就只能检查自己程序是不是有潜在错误,这种情况比较少,其他也没有什么很好的办法。
2010年11月08日 00:18
刚才评论忘记登录了 = =
2010年11月08日 04:06
我也是算法业余爱好者,自己做些题,希望以后能多多交流~
2010年12月22日 05:10
能否给我发一份测试数据?我写了个,还不知对错。谢谢。我的邮箱是roosephu@gmail.com
2013年8月14日 18:01
少年你已经逆天了。。。。。。
2021年1月28日 22:23
感谢博主的分享,请问最大容量路径算法的类 Dijkstra 算法是不是将容量取为赋值,计算路径中容量最小值求得容量最大路径,方便给一下示例程序吗?感谢!
2022年8月22日 02:19
There are several public and private schools in Punjab that are associated with this board; all of these schools operate under its supervision and are also controlled by it. Many students from this state took the class 11th exams this year as well, and they are all currently waiting for the Punjab +1 PSEB 11th Question Paper 2023 of these exams. Punjab Plus One Previous Paper 2023 Punjab Board will as soon as possible announce the Punjab +1 Important Question Paper 2023 for class 11th via their official website. Numerous pupils took the class 11th exams that were administered by this board in the previous academic year, which was done in the past. In this test, girls performed better than guys.
2024年1月17日 21:34
오피타임's serene environment is my escape from daily stress.
2024年1月18日 04:49
제주안마's ability to relieve stress, ease pain, and enhance mental focus makes it a top choice for those looking to unwind on Jeju Island.