网络流算法

     将每条有向边想象成传输物质的管道。每个管道都有一个固定的容量，可以看作是物质能流经该管道的最大速度（譬如可以想象为水流和河槽）。顶点是管道间的连接点，除了源点（S，Source）和汇点（T，Target）以外，物质只流经这些顶点。而不聚集在顶点中。

     注：下文提到的数字，基本都可加单位 “单位流量”来理解。

一、网络流基础

     网络流中的最大流研究的问题是：在不违背容量限制的条件下，把物质从源点传输到汇点的最大速率是多少？上边的图示中，还有些需要说明，设G = （V，E）是一个流网络，其容量函数为c（c(u, v)，Capicity）。s 为源点，t 为汇点。G 的流是一个函数 f （Flow）。上图中每条边标记的是流网络的容量 c(u, v)，右图中，G中的流 |f| = 19。图中只显示正网络流。如果 f（u，v） > 0，则标记（u，v）为 f（u，v）/ c（u，v）（斜杠仅仅用来区分流和容量，不表示相除）。如果 f（u，v）<= 0，边（u，v）只标记它的容量。譬如（s，v₁）这条边，最大流量限制是16，但实际流过11，所以表示为11 / 16。

    网络流算法要基于三种思想：残留网络（Residual Network），增广路径（Augmenting Path）和割（Cut）。

     和开头同一个例子，上图，图（b）就是残留网络，仔细观察一下，应该就可以明白，譬如边（s，v1），从 s 到 v1 流过11，还剩5，所以残余容量是5，当然，流量是守恒的，所以反向从 v1 到 s 流过11。剩余的容量 + 反向平衡的流量共同构成了残留网络。对于名词“残留容量（Residual Capacity）”的定义：在不超过容量c（u，v）的条件下，从 u 到 v 之间可以压入的额外网络流量，就是（u，v）的残留容量（Residual Capacity），公式定义是：c_f（u，v） = c（u，v） – f（u，v）。而残留网络G_f = （V，E_f）。

     图（b），阴影覆盖的边为增广路径 p ，其残留容量为 c_f（p） = c_f（v₂，v₃） = 4。看图应该就可以大概理解什么是增广路了。看图可以发现 s 到 v₂还可以通过 5，v₂到v₃还可以通过 4，v₃到 t 还可以流 5，按照常识，也可得出，这条路径还可以再流过 4 的结论。而从 s 输送 4 到 t 的这条路就是增广路。"增广路径 p”的定义：p 为残留网络G_f中从 s 到 t 的一条简单路径。在不违反边的容量限制条件下，增广路径上的每条边（u，v）可以容纳从 u 到 v 的某额外正网络流。Starvae师兄在讲解增广路时，有一段通俗的描述：

     假如有这么一条路，这条路从源点开始一直一段一段的连到了汇点，并且，这条路上的每一段都满足流量<容量，注意，是严格的<,而不是<=。那么，我们一定能找到这条路上的每一段的(容量-流量)的值当中的最小值delta。我们把这条路上每一段的流量都加上这个delta，一定可以保证这个流依然是可行流。这样我们就得到了一个更大的流，他的流量是之前的流量+delta，而这条路就叫做增广路。

另外，这里得出的，计算增广路流量的公式非常重要：c_f（p） = min { c_f（u，v）：（u，v）在 p 上}，因为最大流算法的核心基本也就是寻找增广路了。

最大流最小割定理：一个流是最大流，当且仅当它的残留网络不包含增广路径。

    流网络的割（S，T）将V划分为 S 和 T = V – S两部分，使得 s ∈ S，t ∈ T。如果 f 是一个流，则穿过割 （S，T）的净流被定义为 f (S，T）。割（S，T）的容量为 c（S，T）。一个网络的最小割就是网络中所有割中最小容量的割。

     上图中，流网络的割（S = {s，v₁，v₂}，T = {v₃，v₄，t}），S中的顶点是黑色，T中的顶点是白色。穿越（S，T）的净流量为：f（v₁，v₃） + f（v₂，v₃） + f（v₂，v₄） = 12 + （-4） + 11 = 19；容量为：c（v₁，v₃） + c（v₂，v₄） = 12 + 14 = 26。可以看出，割的净流由双向的网络流组成。而割的容量仅由 S 到 T 的边计算而得。

     还有一个很有重要的知识：反向边。如图（a），1是源点，4是汇点。很明显如果第一次迭代找到的增广路是1→2→4和1→3→4，则最大流是2，但是如果第一次迭代找到的增广路是1→2→3→4，那么流量只有1，如图（b）是残留网络，这时候，反向边就起作用了，由于反向边的原因，第二次迭代的时候，又找到一条增广路1→2→3→4。这样下来，总的流量还是2。还是Starvae师兄的解释：

      当我们第二次的增广路走 3→2 这条反向边的时候，就相当于把 2→3 这条正向边已经是用了的流量给“退”了回去，不走 2→3 这条路，而改走从 2 点出发的其他的路也就是 2→4。（有人问如果这里没有 2→4 怎么办，这时假如没有 2→4 这条路的话，最终这条增广路也不会存在，因为他根本不能走到汇点）同时本来在 3→4 上的流量由 1→3→4 这条路来“接管”。而最终 2→3 这条路正向流量1，反向流量1，等于没有流量。反向边的作用就是给程序一个可以后悔的机会。—— 一语中的啊。（这句话我加的^^）

二、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 p from s to t in the residual network G_f

5      do c_f(p) ← min {c_f(u, v) : (u, v) is in p}

6         for each edge (u, v) in p

7             do f[u, v] ← f[u, v] + c_f(p)

8                f[v, u] ← -f[u, v]

     从伪代码可以看出，Ford-Fulkerson的核心过程就是：第4~8行的 while 循环反复找出 G_f 中的增广路径 p，并把沿 p 的流 f 加上其残留容量c_f（p）。当不再有增广路径时，流 f 就是一个最大流。完整图示：

   （f）最后在 while 循环测试的残留网络，它没有增光路径，所以（e）中显示的流就是最大流。

三、压入重标记push_relabel算法O(VE)

Push-relabel用到一个很有趣的概念一Preflow(前置流)他允许流进的量比流出的量还要多，有水来就先流进来，流不出去再说。

Push-relabel算法的流程如下：  

1 )我们先假定一个高度函数 h ( u ) ，他代表u点的高度，只有

h ( u )比较高的点才能够将水流到h ( u )比较低的点；  

2 ) 在程序一开始的时候，让source node的高度是n ( node数) ，  其它点的高度都是 0 ，这样source node才有足够的高度可以流往其它地方 ；  

3 ) 然后， 让source node往其它所有跟他直接相邻的node ，都流水管宽度的水量 ( 流过去之后当然要计算剩下的网络情形,即计算 residual  edges(残留网络)

4 )对所有active node( 目前有水的 node )做 relabel（重标记）的动作：  在当某个node明明有水，但是他所连出去的所有对象的 h ( u ) 都比他还高， 则让他的h ( u ) 增加为至少有一条水管可以流出去的量，也就是让这个有水的 active node的高度变成比他连往的”高度最小的 n o d e ”+l ， ( 流过去之后还是要计算剩下的网络情形

5 )对所有可以做push（压入）动作的node做push的动作。  所谓的Push动作是指 ： 当某个node有水，并且他有可以流出去的边， 且他刚好比可以流出去的那个点高度高一点点 ( 高度恰好比他高 1 ) ，那就把某个node 的水流过去，要流多少呢?以下两者取 min。  流出去的水管的量( 也就是说，这个active node的水量很多，  足够把这条水管塞满( 饱和) ，这个时候就叫做saturating Push)（饱和压入）某个 n o d e 现在的水量( 这个 n o d e的水量不足以把流出去的这个水管填满，称作non saturating Push（不饱和压入）  

6 )重复Relabel和Push的工作，一直到没有active node为止，此时从source node所流出的总流量( P r e f l o w) ，就是这个图的最大流量。  

Push-relabel  algorithm 提供了最大流另一方向的思考，且就效率而言，Push -Relabel的复杂度为 o (vE )

这里以POJ 1459为例

[cpp]  view plain  copy

1. #define MIN INT_MIN  
2. #define MAX INT_MAX  
3. #define N 110  
4. int min(int a,int b){return a>b?b:a;}  
5. int c[N][N];//残留容量  
6. int ef[N];//顶点余流  
7. int h[N];//顶点高度  
8. int n;  
9. int push_relabel(int s,int t){  
10.     int i,j;  
11.     int ans = 0;  
12.     memset(h,0,sizeof(h));  
13.     h[s] = t+1;//源点初始高度  
14.     memset(ef,0,sizeof(ef));  
15.     ef[s] = MAX;//源点初始余流  
16.     queue<int> qq;  
17.     qq.push(s);  
18.     while(!qq.empty()){  
19.         int u = qq.front();  
20.         qq.pop();  
21.         for(i=0;i<=t;i++){  
22.             int p;  
23.             int v = i;  
24.             if(c[u][v]<ef[u])p = c[u][v];  
25.             else p = ef[u];  
26.             if(p>0 && (u==s || h[u] == h[v] +1)){  
27.                 c[u][v] -= p;  
28.                 c[v][u] += p;  
29.                 if(v==t)ans+=p;//如果到达了汇点，就将流值加入到最大流中  
30.                 ef[u] -= p;  
31.                 ef[v] += p;  
32.                 if(v!=s && v!=t)qq.push(v);//只有既不是源点也不是汇点才进队  
33.             }  
34.         }  
35.         //如果不是源点且仍有余流，则重标记高度再进队。  
36.         //这里只是简单的将高度增加了一个单位，也可以像上面所说的一样赋值为最低的相邻顶点的高度高一个单位  
37.         if(u!= s && u!=t && ef[u]>0) {  
38.             h[u]++;  
39.             qq.push(u);  
40.         }  
41.     }  
42.     return ans;  
43. }  
44. int main(){  
45.     int np,nc,m;  
46.     while(scanf("%d%d%d%d",&n,&np,&nc,&m) != -1){  
47.         int s = n,t = n+1;  
48.         int i,j;  
49.         memset(c,0,sizeof(c));  
50.         char ss[30];  
51.         for(i=0;i<m;i++){  
52.             int u,v,w;  
53.             scanf("%s",ss);  
54.             sscanf(ss,"(%d,%d)%d",&u,&v,&w);  
55.             c[u][v] += w;  
56.         }  
57.         for(i=0;i<np;i++){  
58.             int u,w;  
59.             scanf("%s",ss);  
60.             sscanf(ss,"(%d)%d",&u,&w);  
61.             c[s][u] += w;  
62.         }  
63.         for(i=0;i<nc;i++){  
64.             int v,w;  
65.             scanf("%s",ss);  
66.             sscanf(ss,"(%d)%d",&v,&w);  
67.             c[v][t] += w;  
68.         }  
69.         printf("%d\n",push_relabel(s,t));  
70.     }  
71.     return 0;  
72. }  

四、Edmonds-Karp（EK）算法  O(V*E*E)

     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 效率是比较低的。实践中此算法使用的机会较少。

这里以 POJ 1273 为例，这里可以作为EK模板，也是我的第一道网络流，mark一下

[cpp]  view plain  copy

1. #define MIN INT_MIN  
2. #define MAX INT_MAX  
3. #define N 204  
4.   
5. int c[N][N];//边容量  
6. int f[N][N];//边实际流量  
7. int pre[N];//记录增广路径  
8. int res[N];//残余网络  
9. queue<int> qq;  
10. void init(){  
11.     while(!qq.empty())qq.pop();  
12.     memset(c,0,sizeof(c));  
13.     memset(f,0,sizeof(f));  
14. }  
15. int EK(int s,int t){  
16.     int i,j;  
17.     int ans=0;  
18.     while(1){  
19.         memset(res,0,sizeof(res));  
20.         res[s] = MAX;//源点的残留网络要置为无限大！否则下面找增广路出错  
21.         pre[s] = -1;  
22.         qq.push(s);  
23.         //bfs找增广路径  
24.         while(!qq.empty()){  
25.             int x = qq.front();  
26.             qq.pop();  
27.             for(i=1;i<=t;i++){  
28.                 if(!res[i] && f[x][i] < c[x][i]){  
29.                     qq.push(i);  
30.                     pre[i] = x;  
31.                     res[i] = min(c[x][i] - f[x][i], res[x]);//这里类似dp，如果有增广路，那么res[t]就是增广路的最小权  
32.                 }  
33.             }  
34.         }  
35.         if(res[t]==0)break;//找不到增广路就退出  
36.         int k = t;  
37.         while(pre[k]!=-1){  
38.             f[pre[k]][k] += res[t];//正向边加上新的流量  
39.             f[k][pre[k]] -= res[t];//反向边要减去新的流量，反向边的作用是给程序一个后悔的机会  
40.             k = pre[k];  
41.         }  
42.         ans += res[t];  
43.     }  
44.     return ans;  
45. }  
46. int main(){  
47.     int n,m;  
48.     while(scanf("%d%d",&n,&m) != -1){  
49.         int i,j;  
50.         init();  
51.         while(n--){  
52.             int a,b,v;  
53.             scanf("%d%d%d",&a,&b,&v);  
54.             c[a][b]+=v;  
55.         }  
56.         printf("%d\n",EK(1,m));  
57.     }  
58.     return 0;  
59. }  

四、Improved SAP（ISAP）算法

      ISAP字面意思是改良的最短增广路算法。关于ISAP，一位叫 DD_engi 的神牛讲非常清楚，引用一下：

     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)，还有一个重要的优化是对于每个点保存“当前弧”：初始时当前弧是邻接表的第一条弧；在邻接表中查找时从当前弧开始查找，找到了一条允许弧，就把这条弧设为当前弧；改变距离标号时，把当前弧重新设为邻接表的第一条弧。

另外，ISAP简化的描述是：程序开始时用一个反向 BFS 初始化所有顶点的距离标号，之后从源点开始，进行如下三种操作：(1)当前顶点 i 为终点时增广 (2) 当前顶点有满足 dist[i] = dist[j] + 1 的出弧时前进 (3) 当前顶点无满足条件的出弧时重标号并回退一步。整个循环当源点 s 的距离标号 dist[s] >= n 时结束。对 i 点的重标号操作可概括为 dist[i] = 1 + min{dist[j] : (i,j)属于残量网络G_f}。

借用庄神的模板 http://www.zlinkin.com/?p=34 ，用它来过POJ 3469简直无敌！比Dinic快好几倍

[cpp]  view plain  copy

1. const int MAXN=20010;  
2. const int MAXM=500010;  
3. int n,m;//n为点数 m为边数  
4. int h[MAXN];  
5. int gap[MAXN];  
6.  int p[MAXN],ecnt;  
7. int source,sink;  
8. struct edge{  
9.     int v;//边的下一点  
10.     int next;//下一条边的编号  
11.     int val;//边权值  
12. }e[MAXM];  
13.   
14. inline void init(){memset(p,-1,sizeof(p));eid=0;}  
15.    
16. //有向  
17. inline void insert1(int from,int to,int val){  
18.     e[ecnt].v=to;  
19.     e[ecnt].val=val;  
20.     e[ecnt].next=p[from];  
21.     p[from]=eid++;  
22.    
23.     swap(from,to);  
24.    
25.     e[ecnt].v=to;  
26.     e[ecnt].val=0;  
27.     e[ecnt].next=p[from];  
28.     p[from]=eid++;  
29. }  
30.    
31. //无向  
32. inline void insert2(int from,int to,int val){  
33.     e[ecnt].v=to;  
34.     e[ecnt].val=val;  
35.     e[ecnt].next=p[from];  
36.     p[from]=eid++;  
37.    
38.     swap(from,to);  
39.    
40.     e[ecnt].v=to;  
41.     e[ecnt].val=val;  
42.     e[ecnt].next=p[from];  
43.     p[from]=eid++;  
44. }  
45.    
46.   
47. inline int dfs(int pos,int cost){  
48.     if (pos==sink){  
49.         return cost;  
50.     }  
51.    
52.     int j,minh=n-1,lv=cost,d;  
53.    
54.     for (j=p[pos];j!=-1;j=e[j].next){  
55.         int v=e[j].v,val=e[j].val;  
56.         if(val>0){  
57.             if (h[v]+1==h[pos]){  
58.                 if (lv<e[j].val) d=lv;  
59.                 else d=e[j].val;  
60.    
61.                 d=dfs(v,d);  
62.                 e[j].val-=d;  
63.                 e[j^1].val+=d;  
64.                 lv-=d;  
65.                 if (h[source]>=n) return cost-lv;  
66.                 if (lv==0) break;  
67.             }  
68.    
69.             if (h[v]<minh)   minh=h[v];  
70.         }  
71.     }  
72.    
73.     if (lv==cost){  
74.         --gap[h[pos]];  
75.         if (gap[h[pos]]==0) h[source]=n;  
76.         h[pos]=minh+1;  
77.         ++gap[h[pos]];  
78.     }  
79.    
80.     return cost-lv;  
81.    
82. }  
83.    
84. int sap(int st,int ed){  
85.    
86.     source=st;  
87.     sink=ed;  
88.     int ans=0;  
89.     memset(gap,0,sizeof(gap));  
90.     memset(h,0,sizeof(h));  
91.    
92.     gap[st]=n;  
93.    
94.     while (h[st]<n){  
95.         ans+=dfs(st,INT_MAX);  
96.     }  
97.    
98.     return ans;  
99. }  

对于EK算法与ISAP算法的区别：

     EK算法每次都要重新寻找增广路，寻找过程只受残余网络的影响，如果改变残余网络，则增广路的寻找也会随之改变；SAP算法预处理出了增广路的寻找大致路径，若中途改变残余网络，则此算法将重新进行。EK处理在运算过程中需要不断加边的最大流比SAP更有优势。