网络流的基本算法——绝恋love枫
引言
过去听起来高深莫测的网络流算法，现在已飞入寻常百姓家了，对于每一个OIER，网络流是一个神圣的东西（个人见解），但神圣的同时，它并不是那样抽象，最形象的模型就是水流，从长江源点无限的向外流水，而大海（汇点）则在不断地‘喝水’，当然，你也可以不把它想成水，或者是其他一切可以流动的东西。而事实上，有些东西的流动比较流畅，而某些东西可能相对而言比较粘稠，流速更慢，因此，就产生了一个问题，单位时间内的总流量最多多少，这里会根据流速给定单位时间内的流量，这就是最先开启网络流之门的最大流算法，它的解决方式将在后面谈到，再想一下，如果水管是另一个物流公司所有，那么你会根据从哪里运到哪里付出一定的代价， 为了你自己的利润，显然要找一个在运的东西最多的前提下有最小费用的方案，这就引出了下一个问题，最小费用最大流。再引用某牛一句话“当然也有有钱没处花的傻子，去求最大费用最大流”，而事实上，题目会出现这个模型，为了避免你成为傻瓜，现在你要给它一个新的定义：最大收益流，这时的你，变成了物流公司的经理，而客户的路线由你规划，为了你的钱包，最大收益必不可少。
正文
第一部分.概念性问题(基本定理及定义)
对于一些网络流新手来说，有必要知道一些基本定义和基本定理，这些虽然看起来理论价值不大，但是现在的许多网络流描述需要这些专业性的词语，所以还是 有些了解为好。
首先对于图G
G 的流是一个实值函数 f， f (u, v) 表示顶点 u 到顶点 v 的流，它可以为正， 为零，也可以为负，且满足下列三个性质：
1.容量限制：对所有u, v ÎV ，要求 f (u, v) £ c(u, v) 。 反对称性：对所有u, v ÎV ，要求 f (u, v) =- f (v, u) 。
2.流守恒性：对所有u ÎV -{s, t} ，要求 å f (u, v) = 0 。
3.整个流网络 G 的流量 f = å f (s, v) 或 f= å f (u, t) 。
接下来定义各种算法中都要用到的一些东东：
1.残留网络
给定一个流网络G = (V , E) 和流 f，由 f 压得的 G 的残留网络Gf= (V , E f ) ，定义 c f (u, v) 为残留网络G f 中边 (u, v) 的容量。如果弧 (u, v) Î E 或弧 (v, u) Î E ，则 弧 (u, v) Î E f ，且 c f (u, v) = c(u, v) - f (u, v) 。
残留网络又被称为剩余图。
2.点的高度和层次，这是两个相对的概念，高度的定义为到汇点的最短路径长度，而层次则是指到源点的最短路径长度（这里的路径长度按照各个边的长度都为1算），这两个量是在最大流算法中贯穿始末的利器。
接下来引入最大流最小割定理
对了，可能有同学还不知道什么是最小割，在这里提一下
流网络 G = (V , E) 的割 (S ,T ) 将V 划分成 S 和T = V - S 两部分，使得 s Î S ，t ÎT 。定义割 (S ,T ) 的容量为 c(S ,T ),
对 于 最 小 的 c ， 它 是 最 小 割 。
3. 最 大 流 最 小 割 定 理
在 流 网 络 中，最 小 割 的 容 量 等 于 最 大 流 的 流 量 。（证 明 再 次 略 过 ）
第二部分.最大流的算法
下面步入与实际问题更加接近的算法实现部分，首先给出问题，给定一个流网络，求源到汇在单位时间内的最大流量。
最简单而效率较好的算法 是基于增广路的算法，这类算法在王欣上大牛的论文中有详细介绍，但我仍然想谈谈我的想法，希望能起到抛砖引玉的作用。基于增广路的算法主要有两种:MPLA,Dinic,SAP.其中最简单的是MPLA,最实用最简洁也是最多人用的是Dinia,SAP的范围也很广，加上GAP优化后的效率也让人咋舌，这也是最近SAP大泛滥的原因吧！个人比较喜欢Dinic，数据变态就用最高标号预流推进，SAP用的比较少，当然，用什么算法还是看你自己的感觉吧。有些人认为增广路算法格式低效，于是想出了对于每个节点操作的算法，这类算法以预留推进为顶梁柱，MPM也勉强归入这一类吧。
1.MPLA算法
即最短路径增值算法，可以有一个简单的思想，每次都找一条从源到汇的路径来增广，直到不能增广为止，之中算法的正确性是可以保证的，但效率不尽如人意，有些时候，把事情格式化反而有益，这里的MPLA就是这样，它只在层次图中找增广路，构建出层次图之后，用BFS不断增广，直到当前层次图中不再有增广路，再重新构建层次图，如果汇点不在层次图内，则源汇不再连通，最大流已经求出，否则继续执行增广，如此反复，就可以求出最大流，在程序实现时层次图不用被构建出来，只需要BFS出各点的距离标号，找路径时判断对于f(u,v)是否有d[u]+1=d[v]即可。
如果每建一次层次图成为一个阶段，则在最短路径增值算法中，最多有N个阶段，证明再次略过。
因此在整个算法中，最多有N个阶段，每个阶段构建层次图的BFS时间复杂度为O(m),建N次，因此构建层次图的总时间为O(mn),而在增广过程中，每一次增广至少删除一条边，因此增广m次，加上修改流量的时间，每一阶段的增广时间为O(m*(m+n)),共有N个阶段，所以复杂度为O(n*m*(m+n))=O(nm^2),这也是该算法的时间复杂度。
2.Dinic算法
MPLA虽然简单，但经常会点超时，我们把增广过程中的BFS改成DFS，效率会有比较大的提高么？答案是肯定的，至此我们已经得到了Dinic的算法流程，只是将MPLA的增广改为DFS，就能写出那美妙的Dinic了，同样，分析一下时间，在DFS过程中，会有前进和后退两种情况，最多前进后退N次，而增广路最多找M次，再加上N个阶段，所以Dinic的复杂度就是O(mn^2)，事实上，它也确实比MPLA快很多，简洁而比较高效，这也是许多OIER选择Dinic的理由了吧，毕竟，写它可能会节省出较长时间来完成其他题目。
View Code
1 program dinic(input,output);
2 var
3 f : array[0..1000,0..1000] of longint;
4 number : array[0..1000] of longint;
5 q : array[0..10000] of longint;
6 n,m,ans,s,t : longint;
7 procedure init;
8 var
9 x,y,c : longint;
10 i : longint;
11 begin
12 readln(m,n);
13 s:=1;
14 t:=n;
15 fillchar(f,sizeof(f),0);
16 for i:=1 to m do
17 begin
18 readln(x,y,c);
19 inc(f[x,y],c);
20 end;
21 end; { init }
22 function min(aa,bb :longint ):longint;
23 begin
24 if aa<bb then
25 exit(aa);
26 exit(bb);
27 end; { min }
28 function bfs(): boolean;
29 var
30 head,tail : longint;
31 now,i : longint;
32 begin
33 fillchar(number,sizeof(number),0);
34 head:=0;
35 tail:=1;
36 q[1]:=s;
37 number[s]:=1;
38 while head<tail do
39 begin
40 inc(head);
41 now:=q[head];
42 for i:=1 to n do
43 if f[now,i]>0 then
44 if number[i]=0 then
45 begin
46 number[i]:=number[now]+1;
47 inc(tail);
48 q[tail]:=i;
49 end;
50 end;
51 if number[t]=0 then
52 exit(false);
53 exit(true);
54 end; { bfs }
55 function dfs(now,flow :longint ):longint;
56 var
57 tmp,i : longint;
58 begin
59 if now=t then
60 exit(flow);
61 for i:=1 to n do
62 if number[i]=number[now]+1 then
63 if f[now,i]>0 then
64 begin
65 tmp:=dfs(i,min(flow,f[now,i]));
66 if tmp<>0 then
67 begin
68 inc(f[i,now],tmp);
69 dec(f[now,i],tmp);
70 exit(tmp);
71 end;
72 end;
73 exit(0);
74 end; { dfs }
75 procedure main;
76 var
77 tmp : longint;
78 begin
79 ans:=0;
80 while bfs() do
81 begin
82 tmp:=dfs(s,maxlongint>>2);
83 while tmp<>0 do
84 begin
85 inc(ans,tmp);
86 tmp:=dfs(s,maxlongint>>2);
87 end;
88 end;
89 writeln(ans);
90 end; { main }
91 begin
92 init;
93 main;
94 end.
3.SAP算法
SAP也是找最短路径来增广的算法，有这样一句话：SAP算法更易理解，实现更简单，效率更高，而也有测试表明，SAP加上重要的GAP优化后，效率仅次于最高标号预流推进算法，因此如果你想背一个模板，SAP是最佳选择。SAP在增光时充分的利用了以前的信息，当按照高度找不到增广路时，它会对节点重新标号，h[i]=min{h[j]}+1(c[i,j]>0)，这也是SAP比较核心的思想，而根据这个我们可以发现，当高度出现间隙时，一定不会存在增广路了，算法已经可以结束，因此，这里引入间隙优化(GAP),即出现间隙时结束算法。
在算法实现中，初始标号可以全部置为0，在增广过程中在逐渐提升高度，时间上可能会有常数的增加，但不改变渐进时间复杂度。同时为了简洁，SAP实现时用递归，代码不过80行左右。
View Code
1 program sap(input,output);
2 var
3 c : array[0..1000,0..1000] of longint;
4 h,vh : array[0..1000] of longint;
5 flow,n,m,ans : longint;
6 tmpflow : longint;
7 can : boolean;
8 procedure init;
9 var
10 i,j : longint;
11 xx,yy,cc : longint;
12 begin
13 fillchar(c,sizeof(c),0);
14 fillchar(h,sizeof(h),0);
15 ans:=0;
16 readln(m,n);
17 for i:=1 to m do
18 begin
19 readln(xx,yy,cc);
20 inc(c[xx,yy],cc);
21 end;
22 end; { init }
23 procedure dfs(now : longint );
24 var
25 min,tmp : longint;
26 i : longint;
27 begin
28 min:=n-1;
29 tmp:=tmpflow;
30 if now=n then
31 begin
32 can:=true;
33 inc(ans,tmpflow);
34 exit;
35 end;
36 for i:=1 to n do
37 if c[now,i]>0 then
38 begin
39 if h[i]+1=h[now] then
40 begin
41 if c[now,i]<tmpflow then
42 tmpflow:=c[now,i];
43 dfs(i);
44 if h[1]>=n then
45 exit;
46 if can then
47 break;
48 tmpflow:=tmp;
49 end;
50 if h[i]<min then
51 min:=h[i];
52 end;
53 if not can then
54 begin
55 dec(vh[h[now]]);
56 if vh[h[now]]=0 then
57 h[1]:=n;
58 h[now]:=min+1;
59 inc(vh[h[now]]);
60 end
61 else
62 begin
63 dec(c[now,i],tmpflow);
64 inc(c[i,now],tmpflow);
65 end;
66 end; { dfs }
67 begin
68 init;
69 fillchar(vh,sizeof(vh),0);
70 vh[0]:=n;
71 while h[1]<n do
72 begin
73 tmpflow:=maxlongint>>2;;
74 can:=false;
75 dfs(1);
76 end;
77 writeln(ans);
78 end.
4.MPM算法
这个算法我还没有实践过，因为它的实现过程比较繁琐，而且时间效率不高，是一个只具有理论价值的算法，这个算法每次都处理单独节点，记每个节点入流和与出流和的最小值作为thoughput(now)（定义在非源汇点），每次先从now向汇推大小为thoughput(now)的流量，在从点now向源点拉大小为thoughput(now)的流量，删除该节点，继续执行直到图中只剩下源汇。时间复杂度为O(n^3)，但时间常数较大，时间效率不高。
5.预留推进算法
以上的算法中，基本上都需要从大体上来把握全局，而预留推进算法则是将每一个顶点看作了一个战场，分别对他们进行处理，在处理过程中，存在某些时间不满足流量收支平衡，所以对预先推出的流叫做预流，下面来看算法如何将预流变成最大流的。
预留推进算法有两个主过程，push和relabel，即推进和重标号，它是在模拟水流的过程，一开始先让源的出弧全部饱和，之后随着时间的推移，不断改变顶点的高度，而又规定水流仅能从高处流向低处，所以在模拟过程中，最终会有水流入汇，而之前推出的多余的水则流回了源，那么我们每次处理的是什么节点呢？把当前节点内存有水的节点称为活跃节点，每次对活跃节点执行推流操作，直到该节点不再活跃，如果不能再推流而当前节点仍未活跃节点，就需要对它进行重新标号了，标号后再继续推流，如此重复，直到网络中不再存在活跃节点为止，这时源的流出量就是该网络的最大流。注意，对于活跃节点的定义，不包括源汇，否则你会死的很惨。
朴素的预留推进的效率还过得去，最多进行nm次饱和推进和n^2m次不饱和推进，因此总的时间复杂度为O(mn^2)
事实上，如同增广路算法引入层次图一样，定下一些规则，可以让预留推进算法有更好的时间效率，下面介绍相对而言比较好实现的FIFO预留推进算法，它用一个队列来保存活跃节点，每次从队首取出一个节点进行推进，对一个节点relabel之后把它加到队尾，如此执行，直到队列为空，这样一来，预留推进算法的时间复杂度降为O(n^3),实现的时候，可以加上同样的间隙优化，但注意，出现间隙时不要马上退出，将新标号的的高度置为n+1，继续执行程序，这样会让所有的剩水流回源，满足流量收支平衡，以便最后的统计工作。
View Code
1 program preflow(input,output);
2 var
3 f,c : array[0..2000,0..2000] of longint;
4 q,h,vh,e : array[0..2000] of longint;
5 m,n,s,t : longint;
6 flow : longint;
7 procedure init;
8 var
9 i,j : longint;
10 xx,yy,cc : longint;
11 begin
12 readln(m,n);
13 fillchar(f,sizeof(f),0);
14 fillchar(c,sizeof(c),0);
15 fillchar(e,sizeof(e),0);
16 for i:=1 to m do
17 begin
18 readln(xx,yy,cc);
19 inc(c[xx,yy],cc);
20 end;
21 s:=1;
22 t:=n;
23 end; { init }
24 procedure main;
25 var
26 i,j : longint;
27 head,tail : longint;
28 now,tmp,tmph : longint;
29 begin
30 flow:=0;
31 h[s]:=n;
32 head:=0;
33 tail:=0;
34 for i:=1 to n do
35 begin
36 e[i]:=c[s,i];
37 f[s,i]:=c[s,i];
38 f[i,s]:=-f[s,i];
39 if (e[i]>0)and(i<>t) then
40 begin
41 inc(tail);
42 q[tail]:=i;
43 inc(vh[h[i]]);
44 end;
45 end;
46 while head<tail do
47 begin
48 inc(head);
49 now:=q[head];
50 for i:=1 to n do
51 if (c[now,i]>f[now,i])and(h[now]=h[i]+1)and(e[now]>0) then
52 begin
53 tmp:=c[now,i]-f[now,i];
54 if tmp>e[now] then
55 tmp:=e[now];
56 inc(f[now,i],tmp);
57 dec(f[i,now],tmp);
58 dec(e[now],tmp);
59 inc(e[i],tmp);
60 if (e[i]=tmp)and(i<>s)and(i<>t) then
61 begin
62 inc(tail);
63 q[tail]:=i;
64 end;
65 end;
66 if (e[now]>0)and(now<>s)and(now<>t) then
67 begin
68 tmph:=h[now];
69 dec(vh[tmph]);
70 h[now]:=$FFFF;
71 for i:=1 to n do
72 if (c[now,i]>f[now,i])and(h[now]>h[i]+1) then
73 h[now]:=h[i]+1;
74 inc(tail);
75 q[tail]:=now;
76 inc(vh[h[now]]);
77 if vh[tmph]=0 then
78 for i:=1 to n do
79 if (h[i]>tmph)and(h[i]<n) then
80 begin
81 dec(vh[h[i]]);
82 h[i]:=n;
83 inc(vh[n]);
84 end;
85 end;
86 end;
87 flow:=0;
88 for i:=1 to n do
89 inc(flow,f[s,i]);
90 end; { main }
91 begin
92 init;
93 main;
94 writeln(flow);
95 end.
下面介绍最后一个，也是编程难度最大，时间表现不同凡响的算法，最高标号预流推进，它的思想是既然水是从高处向低处流的，那么如果从低处开始会做许多重复工作，不如从最高点开始流，留一次就解决问题。再直观一些，引用黑书上的话“让少数的节点聚集大量的盈余，然后通过对这些节点的检查把非饱和推进变成一串连续的饱和推进”。在程序现实现时，用一个表list来储存所有的活跃节点，其中list(h)存储高的为h的活跃节点，同时记录一个level,为最高标号，每次查找时依次从level,level-1……查找，直到找到节点为止，这时从表内删掉这个节点，对它进行Push,Relabel操作，直到该节点不再活跃，继续进行，直到表内不在存在活跃节点。
它的复杂度为O(n^2*m^(1/2)),时间效率很优秀（当然，如果你刻意构造卡预留推进的数据，它比MPLA还慢也是有可能的）。
View Code
1 program hign_node_flow(input,output);
2 var
3 c : array[0..1000,0..1000] of longint; {保存原图}
4 f : array[0..1000,0..1000] of longint; {保存当前的预流图}
5 h : array[0..1000] of longint; {保存各个节点当前高度}
6 vh : array[0..1000] of longint; {保存各个高度节点的数量}
7 e : array[0..1100] of longint; {保存各个节点的盈余}
8 level : longint; {当前所有活跃节点的最高高度}
9 l : array[0..1000,0..1000] of longint; {保存活跃节点的表,l[i,0]表示高度为i的活跃节点数，这也是不能用vh数组的原因}
10 n,m,s,t : longint; {节点数，边数，源，汇}
11 listsum : longint; {记录当前在表内的元素个数}
12 flow : longint; {记录流量}
13 inlist : array[0..1000] of boolean; {节点是否在表内}
14 q : array[0..10000] of longint; {用于BFS扩展的队列}
15 procedure init;
16 var
17 i,xx,yy,cc : longint;
18 begin
19 readln(m,n);
20 fillchar(f,sizeof(f),0);
21 fillchar(c,sizeof(c),0);
22 fillchar(e,sizeof(e),0);
23 fillchar(h,sizeof(h),0);
24 fillchar(vh,sizeof(vh),0);
25 for i:=1 to m do
26 begin
27 readln(xx,yy,cc);
28 inc(c[xx,yy],cc);{注意某些情况下有重边，这样处理比较保险}
29 end;
30 s:=1;
31 t:=n;
32 end; { init }
33 procedure insect(now :longint ); {在活跃节点表内插入节点now}
34 begin
35 inlist[now]:=true; {标记now节点在表内}
36 inc(listsum); {表中元素增加1}
37 inc(l[h[now],0]); {高度为h[now]的活跃节点数增加1}
38 l[h[now],l[h[now],0]]:=now; {表中高度为h[now]的第l[h[now],0]个活跃节点为now}
39 if h[now]>level then {更新活跃节点最高高度}
40 level:=h[now];
41 end; { insect }
42 procedure bfs(); {利用BFS(反向的)，求的各个节点的高度}
43 var
44 head,tail,i : longint;
45 begin
46 head:=0;
47 tail:=1;
48 q[1]:=t;
49 h[t]:=1; {汇点的高度为1}
50 while head<tail do
51 begin
52 inc(head);
53 for i:=1 to n do
54 if c[i,q[head]]>0 then {存在边}
55 if h[i]=0 then {i节点高度没有求出}
56 begin
57 h[i]:=h[q[head]]+1; {求的节点i的高度}
58 inc(tail);
59 q[tail]:=i;
60 end;
61 end;
62 end; { bfs }
63 procedure previous(); {预流推进的预处理}
64 var
65 i : longint;
66 begin
67 for i:=1 to n do
68 begin
69 e[i]:=c[s,i]; {让源点的出弧饱和，则弧的指向点的盈余要改变}
70 f[s,i]:=c[s,i]; {源点出弧饱和}
71 f[i,s]:=-f[s,i]; {反向弧的处理}
72 if (e[i]>0)and(i<>t)and(not inlist[i]) then {节点i成为活跃节点,且不是汇点，没有在表内(其实也不可能在表内)}
73 insect(i);
74 end;
75 h[1]:=n;
76 for i:=1 to n-1 do
77 inc(vh[h[i]]);
78 end; { previous }
79 functio