USACO 421 Ditch 网络最大流问题算法小结.docx
《USACO 421 Ditch 网络最大流问题算法小结.docx》由会员分享，可在线阅读，更多相关《USACO 421 Ditch 网络最大流问题算法小结.docx（11页珍藏版）》请在冰豆网上搜索。
USACO421Ditch网络最大流问题算法小结
USACO4.2.1Ditch网络最大流问题算法小结
通过USACO4.2.1Ditch学习一下最大流算法。
可惜它给的测试数据几乎没有任何杀伤力，后面测试时我们采用DD_engi写的程序生成的加强版数据。
总体上来说，最大流算法分为两大类：
增广路（AugmentingPath）和预流推进重标号（PushRelabel）。
也有算法同时借鉴了两者的长处，如ImprovedSAP。
本篇主要介绍增广路类算法，思想、复杂度及实际运行效率比较，并试图从中选择一种兼顾代码复杂度和运行效率的较好方案。
以下我们将会看到，有时理论分析的时间复杂度并不能很好的反映一种算法的实际效率。
1.Ford-Fulkerson方法
所有增广路算法的基础都是Ford-Fulkerson方法。
称之为方法而不是算法是因为Ford-Fulkerson只提供了一类思想，在此之上的具体操作可有不同的实现方案。
给定一个有向网络G（V,E）以及源点s终点t，FF方法描述如下：
Ford-Fulkerson方法（G,s,t）
1将各边上流量f初始化为0
2while存在一条增广路径p
3do沿路径p增广流量f
4returnf
假设有向网络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的最短路径，或者流量最大的路径。
2.Edmonds-Karp算法
ShortestAugmentingPath（SAP）是每次寻找最短增广路的一类算法，Edmonds-Karp算法以及后来著名的Dinic算法都属于此。
SAP类算法可统一描述如下：
ShortestAugmentingPath
1x<--0
2while在残量网络Gx中存在增广路s~>t
3do找一条最短的增广路径P
4delta<--min{rij:
（i,j）属于P}
5沿P增广delta大小的流量
6更新残量网络Gx
7returnx
在无权边的有向图中寻找最短路，最简单的方法就是广度优先搜索（BFS），E-K算法就直接来源于此。
每次用一遍BFS寻找从源点s到终点t的最短路作为增广路径，然后增广流量f并修改残量网络，直到不存在新的增广路径。
E-K算法的时间复杂度为O（VE2），由于BFS要搜索全部小于最短距离的分支路径之后才能找到终点，因此可以想象频繁的BFS效率是比较低的。
实践中此算法使用的机会较少。
3.Dinic算法
BFS寻找终点太慢，而DFS又不能保证找到最短路径。
1970年Dinic提出一种思想，结合了BFS与DFS的优势，采用构造分层网络的方法可以较快找到最短增广路，此算法又称为阻塞流算法（BlockingFlowAlgorithm）。
首先定义分层网络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部分算法可描述如下：
1p<--s
2whiles的出度>0do
3u<--p.top
4ifu!
=tthen
5ifu的出度>0then
6设（u,v）为AN（f）中一条边
7p<--p,v
8else
9从p和AN（f）中删除点u以及和u连接的所有边
10else
11沿p增广
12令p.top为从s沿p可到达的最后顶点
13endwhile
实际代码中不必真的用一个图来存储分层网络，只需保存每个顶点的距离标号并在DFS时判断dist[i]+1=dist[j]即可。
Dinic的时间复杂度为O（V2E）。
由于较少的代码量和不错的运行效率，Dinic在实践中比较常用。
具体代码可参考DD_engi网络流算法评测包中的标程，这几天dinic算法的实现一共看过和比较过将近10个版本了，DD写的那个在效率上是数一数二的，逻辑上也比较清晰。
4.ImprovedSAP算法
本次介绍的重头戏。
通常的SAP类算法在寻找增广路时总要先进行BFS，BFS的最坏情况下复杂度为O（E），这样使得普通SAP类算法最坏情况下时间复杂度达到了O（VE2）。
为了避免这种情况，Ahuja和Orlin在1987年提出了ImprovedSAP算法，它充分利用了距离标号的作用，每次发现顶点无出弧时不是像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}。
具体算法描述如下：
algorithmImproved-Shortest-Augmenting-Path
1f<--0
2从终点t开始进行一遍反向BFS求得所有顶点的起始距离标号d（i）
3i<--s
4whiled（s）<ndo
5ifi=tthen//找到增广路
6Augment
7i<--s//从源点s开始下次寻找
8ifGf包含从i出发的一条允许弧（i,j）
9thenAdvance（i）
10elseRetreat（i）//没有从i出发的允许弧则回退
11returnfprocedureAdvance（i）
1设（i,j）为从i出发的一条允许弧
2pi（j）<--i//保存一条反向路径，为回退时准备
3i<--j//前进一步，使j成为当前结点procedureRetreat（i）
1d（i）<--1+min{d（j）:
（i,j）属于残量网络Gf}//称为重标号的操作
2ifi!
=s
3theni<--pi（i）//回退一步procedureAugment
1pi中记录为当前找到的增广路P
2delta<--min{rij:
（i,j）属于P}
3沿路径P增广delta的流量
4更新残量网络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行。
5.最大容量路径算法（MaximumCapacityPathAlgorithm）
1972年还是那个E-K组合提出的另一种最大流算法。
每次寻找增广路径时不找最短路径，而找容量最大的。
可以预见，此方法与SAP类算法相比可更快逼近最大流，从而降低增广操作的次数。
实际算法也很简单，只用把前面E-K算法的BFS部分替换为一个类Dijkstra算法即可。
USACO4.2节的说明详细介绍了此算法，这里就不详述了。
时间复杂度方面。
BFS是O（E），简单Dijkstra是O（V2），因此效果可想而知。
但提到Dijkstra就不能不提那个Heap优化，虽然USACO的算法例子中没有用Heap，我自己还是实现了一个加Heap的版本，毕竟STL的优先队列太好用了不加白不加啊。
效果也是非常明显的，但比起Dinic或ISAP仍然存在海量差距，这里就不再详细介绍了。
6.CapacityScalingAlgorithm
不知道怎么翻比较好，索性就这么放着吧。
叫什么的都有，容量缩放算法、容量变尺度算法等，反正就那个意思。
类似于二分查找的思想，寻找增广路时不必非要局限于寻找最大容量，而是找到一个可接受的较大值即可，一方面有效降低寻找增广路时的复杂度，另一方面增广操作次数也不会增加太多。
时间复杂度O（E2logU）实际效率嘛大约稍好于最前面BFS的E-K算法，稀疏图时表现较优，但仍然不敌Dinic与ISAP。
7.算法效率实测！
重头戏之二，虽然引用比较多，哎~
首先引用此篇强文《MaximumFlow:
AugmentingPathAlgorithmsComparison》
对以上算法在稀疏图、中等稠密图及稠密图上分别进行了对比测试。
直接看结果吧：
稀疏图：
ISAP轻松拿下第一的位置，图中最左边的SAP应该指的是E-K算法，这里没有比较Dinic算法是个小遗憾吧，他把Dinic与SAP归为一类了。
最大流量路径的简单Dijkstra实现实在是太失败了--，好在Heap优化后还比较能接受……可以看到Scaling算法也有不错的表现。
中等稠密图：
ISAP依然领先。
最大流量算法依然不太好过……几个Scaling类算法仍然可接受。
稠密图：
ISAP……你无敌了！
这次可以看出BFS的Scaling比DFS实现好得多，而且几乎与ImprovedScaling不相上下，代码量却差不少。
看来除ISAP外BFSScaling也是个不错的选择。
最后补个我自己实测的图，比较算法有很多是DD网络流算法评测包里的标程，评测系统用的Cena，评测数据为DDditch数据生成程序生成的加强版数据：
我一共写了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的非递归版本代码，网络采用邻接表+反向弧指针：
#include<cstdio>
#include<memory>
usingnamespacestd;constintmaxnode=1024;
constintinfinity=2100000000;structedge{
intver;//vertex
intcap;//capacity
intflow;//currentflowinthisarc
edge*next;//nextarc
edge*rev;//reversearc
edge（）{}
edge（intVertex,intCapacity,edge*Next）
:
ver（Vertex）,cap（Capacity）,flow（0）,next（Next）,rev（（edge*）NULL）{}
void*operatornew（size_t,void*p）{
returnp;
}
}*Net[maxnode];
intdist[maxnode]={0},numbs[maxnode]={0},src,des,n;voidrev_BFS（）{
intQ[maxnode],head=0,tail=0;
for（inti=1;i<=n;++i）{
dist[i]=maxnode;
numbs[i]=0;
}Q[tail++]=des;
dist[des]=0;
numbs[0]=1;
while（head!
=tail）{
intv=Q[head++];
for（edge*e=Net[v];e;e=e->next）{
if（e->rev->cap==0||dist[e->ver]<maxnode）continue;
dist[e->ver]=dist[v]+1;
++numbs[dist[e->ver]];
Q[tail++]=e->ver;
}
}
}intmaxflow（）{
intu,totalflow=0;
edge*CurEdge[maxnode],*revpath[maxnode];
for（inti=1;i<=n;++i）CurEdge[i]=Net[i];
u=src;
while（dist[src]<n）{
if（u==des）{//findanaugmentingpath
intaugflow=infinity;
for（inti=src;i!
=des;i=CurEdge[i]->ver）
augflow=min（augflow,CurEdge[i]->cap）;
for（inti=src;i!
=des;i=CurEdge[i]->ver）{
CurEdge[i]->cap-=augflow;
CurEdge[i]->rev->cap+=augflow;
CurEdge[i]->flow+=augflow;
CurEdge[i]->rev->flow-=augflow;
}
totalflow+=augflow;
u=src;
}edge*e;
for（e=CurEdge[u];e;e=e->next）
if（e->cap>0&&dist[u]==dist[e->ver]+1）break;
if（e）{//findanadmissiblearc,thenAdvance
CurEdge[u]=e;
revpath[e->ver]=e->rev;
u=e->ver;
}else{//noadmissiblearc,thenrelabelthisvertex
if（0==（--numbs[dist[u]]））break;//GAPcut,Important!
CurEdge[u]=Net[u];
intmindist=n;
for（edge*te=Net[u];te;te=te->next）
if（te->cap>0）mindist=min（mindist,dist[te->ver]）;
dist[u]=mindist+1;
++numbs[dist[u]];
if（u!
=src）
u=revpath[u]->ver;//Backtrack
}
}
returntotalflow;
}intmain（）{
intm,u,v,w;
freopen（"ditch.in","r",stdin）;
freopen（"ditch.out","w",stdout）;
while（scanf（"%d%d",&m,&n）!
=EOF）{//POJ1273needthiswhileloop
edge*buffer=newedge[2*m];
edge*data=buffer;
src=1;des=n;
while（m--）{
scanf（"%d%d%d",&u,&v,&w）;
Net[u]=new（（void*）data++）edge（v,w,Net[u]）;
Net[v]=new（（void*）data++）edge（u,0,Net[v]）;
Net[u]->rev=Net[v];
Net[v]->rev=Net[u];
}
rev_BFS（）;
printf（"%d\n",maxflow（））;
delete[]buffer;
}
return0;
}
升级会员 