Max Flow: Ford-Fulkerson方法

http://chuanwang66.iteye.com/blog/1446977

2. 网络流：

（1）容量限制：f(u,v)≤c(u,v)                                 ==> 单向 流速受限

（2）反对称性：f(u,v)=-f(v,u)                                ==> 在管道中向不同方向看，水流一面迎面而来 ，一面向前推我

（3）流守恒性：f(V-s-t,V)=0                                  ==> A. 中间顶点不存储流 （流进=流出 ）

                                                                                   B. 进入中间顶点的正网络流=离开该点的正网络流

注意：

       流守恒性也可以写作f(V,V-s-t)=0，即流入一个顶点的总流为0；

       在f(V-s-t,V)或f(V,V-s-t)的中，其中每一元素可以为正可以为负，不要以为不能为负（这是“残留网络”中的限制，网络流不受限）

引理26.1: G=<V,E>是流网络，f是G的一个流。

               定义f(X,Y)=∑f(x,y)，x∈X, y∈Y.

              （1）任意X⊆V,                         则 f(X,X)=0                        ==>反对称性:  f(u,v)= - f(v,u)

              （2）任意X,Y⊆V,                      则f(X,Y)=-f(Y,X)                  ==>（1）的推广

             （3）X,Y,Z⊆ V, X∩Y=∅             则f(X∪Y,Z)=f(X,Z)+f(Y,Z)

最大流最小割定理 ：

    （1）f是G的一个最大流                                     ==>达到流值|f|=f(s,V)最大

    （2）残留网络Gf不包含增广路径                         ==>不能再压入正网络流

    （3）对G的某个割(S,T) ，存在|f|=c(S,T)            ==>对最大流的限制来自最小容量的那个割

http://blog.csdn.net/smartxxyx/article/details/9293665

一、残留网络

顾名思义，残留网络是指给定网络和一个流，其对应还可以容纳的流组成的网络。具体说来，就是假定一个网络G=（V，E），其源点s，汇点t。设f为G中的一个流，对应顶点u到顶点v的流。在不超过C（u，v）的条件下（C代表边容量），从u到v之间可以压入的额外网络流量，就是边（u，v）的残余容量（residual capacity），定义如下：

r（u，v）=c（u，v）-f（u，v）

举个例子，假设（u，v）当前流量为3/4，那么就是说c（u，v）=4，f（u，v）=3，那么r（u，v）=1。

我们知道，在网络流中还有这么一条规律。从u到v已经有了3个单位流量，那么从反方向上看，也就是从v到u就有了3个单位的残留网络，这时r（v，u）=3。可以这样理解，从u到v有3个单位流量，那么从v到u就有了将这3个单位流量的压回去的能力。

我们来具体看一个例子，如下图所示一个流网络

其对应的残留网络为：

二、增广路径

在了解了残留网络后，我们来介绍增广路径。已知一个流网络G和流f，增广路径p是其残留网络Gf中从s到t的一条简单路径。形象的理解为从s到t存在一条不违反边容量的路径，向这条路径压入流量，可以增加整个网络的流值。上面的残留网络中，存在这样一条增广路径：

其可以压入4个单位的流量，压入后，我们得到一个新的流网络，其流量比原来的流网络要多4。这时我们继续在新的流网络上用同样的方法寻找增广路径，直到找不到为止。这时我们就得到了一个最大的网络流。

三、流网络的割

上面仅仅是介绍了方法，可是怎么证明当无法再寻找到增广路径时，就证明当前网络是最大流网络呢？这就需要用到最大流最小割定理。

首先介绍下，割的概念。流网络G（V，E）的割（S，T）将V划分为S和T=V-S两部分，使得s属于S，t属于T。割（S，T）的容量是指从集合S到集合T的所有边（有方向）的容量之和（不算反方向的，必须是S-àT）。如果f是一个流，则穿过割（S，T）的净流量被定义为f（S，T）（包括反向的，SàT的为正值，T—>S的负值）。将上面举的例子继续拿来，随便画一个割，如下图所示：

割的容量就是c(u,w)+c(v,x)=26

当前流网络的穿过割的净流量为f(u,w)+f(v,x)-f(w,v)=12+11-4=19

显然，我们有对任意一个割，穿过该割的净流量上界就是该割的容量，即不可能超过割的容量。所以网络的最大流必然无法超过网络的最小割。

流网络的流量守恒的原则，根据这个原则我们可以知道，对网络的任意割，其净流量的都是相等的。

和上面的割相比，集合S中少了u和v，从源点s到集合T的净流量都流向了u和v，而在上一个割图中，集合S到集合T的流量是等于u和v到集合T的净流量的。其中w也有流流向了u和v，而这部分流无法流向源点s，因为没有路径，所以最后这部分流量加上s到u和v的流量，在u和v之间无论如何互相传递流，最终都要流向集合T，所以这个流量值是等于s流向u和v的值的。将s比喻成一个水龙头，u和v流向别处的水流，都是来自s的，其自身不可能创造水流。所以任意割的净流量都是相等的。

万事俱备，现在来证明当残留网络Gf中不包含增广路径时，f是G的最大流。

假设Gf中不包含增广路径，即Gf不包含从s到v的路径，定义S={v:Gf中从s到v存在一条通路}，也就是Gf中s能够有通路到达的点的集合，显然这个集合不包括t，因为s到t没有通路。这时，我们令T=V-S。那么（S，T）就是一个割。如下图所示：

那么，对于顶点u属于S，v属于T，有f（u，v）=c（u，v）。否则（u，v）就存在残余流量，因而s到u加上u到v就构成了一条s到v的通路，所以v就必须属于S，矛盾。因此这时就表明当前流f是等于当前的割的容量的，因此f就是最大流。

图的匹配问题与最大流问题（三）——最大流问题Ford-Fulkerson方法Java实现

FORD-FULKERSON（G，t，s）

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 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]

如果在4行中用广度优先搜索来实现对增广路径p的计算，即找到s到t的最短增广路径，能够改进FORD-FULERSON的界，这就是Ford-Fulkerson方法的Edmonds-Karp算法

证明该算法的运行时间为O（VE*E），易知，对流增加的全部次数上界为O（VE），每次迭代时间O（E）

1. public class FordFulkerson  
2.     {  
3.         private int residualNetwork[][]=null;  
4.         private int flowNetwork[][]=null;  
5.           
6.         public final int N;  
7.         int parent[];  
8.         public FordFulkerson(int N)  
9.         {  
10.             this.N=N;  
11.             parent=new int[N];  
12.         }  
13.         /** 
14.          * 实现FordFulkerson方法的一种算法——edmondsKarp算法 
15.          * @param graph 
16.          * @param s 
17.          * @param t 
18.          * @return 
19.          */  
20.         public int edmondsKarpMaxFlow(int graph[][],int s,int t)  
21.         {  
22.             int length=graph.length;  
23.             int f[][]=new int[length][length];  
24.             for(int i=0;i<length;i++)  
25.             {  
26.                 Arrays.fill(f[i], 0);  
27.             }  
28.             int r[][]=residualNetwork(graph,f);  
29.             int result=augmentPath(r,s,t);  
30.               
31.             int sum=0;  
32.               
33.             while(result!=-1)  
34.             {  
35.                 int cur=t;  
36.                 while(cur!=s)  
37.                 {  
38.                     f[parent[cur]][cur]+=result;  
39.                     f[cur][parent[cur]]=-f[parent[cur]][cur];  
40.                     r[parent[cur]][cur]-=result;  
41.                     r[cur][parent[cur]]+=result;  
42.                     cur=parent[cur];  
43.                 }  
44.                   
45.                 sum+=result;  
46.                 result=augmentPath(r,s,t);  
47.             }  
48.               
49.             residualNetwork=r;  
50.             flowNetwork=f;  
51.               
52.             return sum;  
53.         }  
54.   
55.         /** 
56.          * deepCopy 
57.          * @param c 
58.          * @param f 
59.          * @return 
60.          */  
61.         private int[][] residualNetwork(int c[][],int f[][]) {  
62.             int length=c.length;  
63.             int r[][]=new int[length][length];  
64.             for(int i=0;i<length;i++)  
65.             {  
66.                 for(int j=0;j<length;j++)  
67.                 {  
68.                     r[i][j]=c[i][j]-f[i][j];  
69.                 }  
70.             }  
71.               
72.             return r;  
73.         }  
74.   
75.         /** 
76.          * 广度优先遍历，寻找增光路径，也是最短增广路径 
77.          * @param graph 
78.          * @param s 
79.          * @param t 
80.          * @return 
81.          */  
82.         public int augmentPath(int graph[][],int s,int t)  
83.         {  
84.               
85.             int maxflow=Integer.MAX_VALUE;  
86.             Arrays.fill(parent, -1);  
87.             Queue<Integer> queue=new LinkedList<Integer>();  
88.             queue.add(s);  
89.             parent[s]=s;  
90.   
91.             while(!queue.isEmpty())  
92.             {  
93.                 int p=queue.poll();  
94.                 if(p==t)  
95.                     break;  
96.                 for(int i=0;i<graph.length;i++)  
97.                 {  
98.                     if(i!=p&&parent[i]==-1&&graph[p][i]>0)  
99.                     {  
100.                         if(maxflow>graph[p][i])  
101.                             maxflow=graph[p][i];  
102.                         //flow[i]=Math.min(flow[p], graph[p][i]);  
103.                         parent[i]=p;  
104.                         queue.add(i);  
105.                     }  
106.                 }  
107.             }  
108.             if(parent[t]==-1)  
109.                 return -1;  
110.             return  maxflow;  
111.               
112.         }  
113.   
114.         public int[][] getResidualNetwork() {  
115.             return residualNetwork;  
116.         }  
117.   
118.         public int[][] getFlowNetwork() {  
119.             return flowNetwork;  
120.         }  
121.     }  

http://chuanwang66.iteye.com/blog/1446977

2. 求解最大流问题

    2.1 枚举算法

          时间复杂度： O(2^|V|·|E|)

          思路：枚举所有割(S,T)，找到容量c(S,T)最小的那个割的容量，即为最大流的流值

          评价：算法复杂度高，因此仅当顶点个数|V|较少时适用（否则整数越界）；另外，它还不能给出得到最大流的具体的网络流，只是返回了最大流值。

1. //source在最低位，对应1；sink在最高位，对应0。  
2.   
3. //sink * * * * * * source  
4. //  0   0 0 0 0 0 0   1  
5. //  0   0 0 0 0 0 1   1  
6. //  0    .  .  .  .  .  .   1  ==>每一种二进制组合对应一个割  
7. //  0   1 1 1 1 1 1   1  
8. //      |                 |  
9. //      ------ i -------  
10. //用int i枚举二进制组合  
11.   
12. int maxflow(){  
13.     int flow=INF;  
14.     int tmp;  
15.     for(int i=0;i<(1<<n-2);i++){  
16.         //枚举每一个割(S,T). S中元素对应位为1, T中元素对应为为0  
17.         s=(i<<1)|1;  
18.         tmp=0;  
19.         for(int u=0;u<n;u++)  
20.             for(int v=0;v<n;v++)  
21.                 if(  ((s>>u)&1)==1  &&  ((s>>v)&1)==0  )  
22.                     tmp+=c[u][v];  
23.         if(flow>tmp) flow=tmp;  
24.     }  
25.     return flow;  
26. } 

2.2 增广路算法

FORD-FULKERSON方法：

1. f=0  
2. while( p exists)  
3.    do f+=p  

 其中，f是流值|f|不断增加的网络流，初始时f=0，之后不断沿着增广路径压入正网络流；

            p是每次找到的增广路径

1. FORD-FULKERSON(G,s,t)  
2.     for each edge(u,v)∈E[G]  
3.          do  f[u,v] <-      0  
4.                f[v,u] <-      0  
5.     while there exists a path p from s to t in the residual network Gf  
6.              do cf(p) <- min{cf(u,v): (u,v) is in p}  
7.                   for  each edge(u,v) in p  
8.                         do f[u,v] <-      f[u,v] + cf(p)      //“沿着”增广路径，有向边上的流“增加”增广路径残留容量cf(p)  
9.                              f[v,u] <-      - f[u,v]               //“逆着”增广路径，有向边上的流“减少”增广路径残留容量cf(p)  

B. BFS搜索增广路径 （即“Edmonds-Karp算法”）

    时间复杂度：O(|V|·|E|^2)      ==>每次搜索增广路径需要O(|E|)；可以增广O(|V||E|)次，原因见下面评价

    思路：每次BFS得到一条最短的增广路径增广

    评价：经证明(independently proved by Dinic in 1970, and by Edmonds & Karp in 1972)，若每次从最短增广路径开始增广，则最多只用进行O(|V||E|)次增广，而不是DFS搜索增广路径的O(c)次。其中“最短增广路径”指的是：当每条边长度都为1时的最短路径，显然，最朴素的想法就是用BFS来找这种最短路径。

BFS搜索增广路径——Edmonds-karp算法：

1. const int N=210;//边  
2. const int INF=0x7FFFFFFF;  
3. int n,m;  
4. int map[N][N];  //残留图，初始状态标识每条有向边的容量的图也可以看做是一个残留图  
5. int pi[N];  //BFS的前驱图; pi[start]=0；pi[i]=-1，如果i没有在BFS中被扩展过  
6. int flow_in[N]; //流入i的最大流量是flow_in[i]  
7. int start,end;  
8. queue<int> q;  
9.   
10. int bfs(){  
11.     int i,t;  
12.     while(!q.empty()) q.pop();  
13.     memset(pi,-1,sizeof(pi));  
14.   
15.     pi[start]=0;  
16.     flow_in[start]=INF; //key!  
17.   
18.     q.push(start);  
19.     while(!q.empty()){  
20.         t=q.front();  
21.         q.pop();  
22.   
23.         if(t==end) break;  
24.   
25.         for(i=1;i<=m;i++){  
26.             if(pi[i]==-1 && map[t][i]){  
27.                 flow_in[i]=flow_in[t]<map[t][i]?flow_in[t]:map[t][i];  
28.                 q.push(i);  
29.                 pi[i]=t;  
30.             }  
31.         }  
32.     }  
33.     if(pi[end]==-1) return -1;  //不存在增广路径  
34.     else return flow_in[m];         //还存在增光路径，不保证flow_in[m]达到最大可能值  
35. }  
36.   
37. int Edmonds_Karp(){  
38.     int max_flow_in=0;  //流f的流值|f|  
39.     int cf_p;       //增广路径的残留容量Cf(p)  
40.     int now,pre;  
41.   
42.     while((cf_p=bfs())!=-1){  
43.         //1. 流值|f|增加本次增广路径的残留容量cf_p  
44.         max_flow_in+=cf_p;   
45.   
46.         //2. 顺着和逆着增广路径压入网络流  
47.         now=end;  
48.         while(now!=start){  
49.             pre=pi[now];  
50.             map[pre][now]-=cf_p;    //更新正向边的实际容量  
51.             map[now][pre]+=cf_p;    //添加反向边  
52.             now=pre;  
53.         }  
54.     }  
55.     return max_flow_in;  
56. } 

A. DFS搜索增广路径

    时间复杂度：O(|E|·c) , c是最大流 ==>每次搜索增广路径需要O(|E|)；可以增广O(c)次，∵每次可以增广大小为1的流

    思路：每次DFS得到一条盲目的增广路径增广

    评价：经证明，若每次从最短增广路径开始增广，则最多只用进行O(|V||E|)次增广，而不是DFS搜索增广路径的O(c)次。从这个意义上讲，DFS搜索增广路径是盲目的。

1. int n;  //有向边数  
2. int m;  //顶点数  
3. int f[N][N];    //网络流  
4. int cf[N][N];   //残留网络  
5. bool vis[N];    //DFS中是否访问过  
6.     
7. const int inf = 1<<29;    
8.     
9. int inline min(int x,int y){return x>y?y:x;}    
10.     
11. bool dfs(int s,int t,int &cf_path)    
12. {    
13.     vis[s]=true;    
14.     if(s==t)  
15.         return true;   
16.   
17.     for(int i=1;i<=m;++i)    
18.     {      
19.         if(vis[i]==false && cf[s][i]>0)    
20.         {    
21.             int temp = cf_path;   
22.   
23.             cf_path=min(cf_path,cf[s][i]);  //(a)式  
24.             if(dfs(i,t,cf_path))    
25.             {    
26.                 f[s][i]+=cf_path;   
27.   
28.                 cf[s][i]-=cf_path;    
29.                 cf[i][s]+=cf_path;  //允许“反悔”  
30.   
31.                 return true;    
32.             }  
33.   
34.             //沿着<s,i>向下深搜失败，回溯时恢复增广路径的残留容量cf[path];  
35.             //实际上是在撤销(a)式的影响!  
36.             cf_path = temp;  
37.         }    
38.     }    
39.     return false;    
40. }    
41.     
42. int find(int s=1,int t=m)    //这实际上是DFS的外围框架（回忆DFS算法实现有内外两层）  
43. {    
44.     int cf_path = inf;      //本次可增广的流量，在DFS过程中被限制得越来越小  
45.     memset(vis,0,sizeof(vis));    
46.   
47.     dfs(s,t,cf_path);  
48.   
49.     if(cf_path>=inf)  
50.         return 0;    
51.     else   
52.         return cf_path;    
53. }    
54.     
55. int max_flow()    
56. {    
57.     int ret=0,i;    
58.     memset(f,0,sizeof(f));    
59.     while(find());   
60.   
61.     for(i=1;i<=m;++i)    
62.         ret+=f[1][i];   //网络流的值=从源出发的流量和f(s,V)=流入汇点的流量和f(V,t)  
63.     return ret;    
64. } 

SAP算法：

1. struct sap    
2. {    
3.     int s,t;    //s:源点 t:汇点    
4.     int m;      //顶点数    
5.     int cf[305][305];   //残留容量    
6.     
7.     int flow;       //当前最大流    
8.     int cf_path;    //本次可增广的流量（本次增广路径的残留容量）  
9.     bool flag;      //本次增广路径是否找到    
10.     int vh[1000];   //各高度的顶点个数    
11.     int h[1000];    //各顶点的高度，又称为“距离标号”——到汇点距离的下界（边权值都为1）  
12.     
13.     sap(){    
14.         flow =cf_path =s =t =0;    
15.         flag=false;    
16.         memset(cf,0,sizeof(cf));    
17.         memset(vh,0,sizeof(vh));    
18.         memset(h,0,sizeof(h));    
19.     }    
20.     
21.     /*find_path_sap(int cur):增广路径已经推进到顶点cur，尝试向下一个顶点推进 
22.        主要思路： 
23.              递归地尝试沿着当前点cur的所有允许弧推进增广路径。 
24.              如果失败，将当前点cur的高度（距离函数下界）更新到minH(minH:残留网络中所有cur指向节点的最小高度) 
25.     */  
26.     void find_path_sap(int cur){    
27.         //1. 检查是否走到增广路径末端    
28.         //*a. 终止条件：cur是本次增广路径末端节点(汇点)    
29.         if(cur==t){    
30.             flow+=cf_path;    
31.             flag=true;    
32.             return;    
33.         }    
34.     
35.         //2. 尝试递归/回溯 找增广路径    
36.         //如果cur仅比邻接点高1，则尝试最大限度向它送出流    
37.         int i;    
38.         int minH=m-1;    
39.         int tmp_cf_path=cf_path;    
40.         for(i=1;i<=m;i++){    
41.             if(cf[cur][i]){    
42.                 if(h[i]+1==h[cur]){    
43.                     if(cf[cur][i]<cf_path)    
44.                         cf_path=cf[cur][i];    
45.                     find_path_sap(i);    
46.                         
47.                     //*b. 终止条件: 如果h[1]增长到m，表示通过本次"find_path_sap(i)"递归搜索，不能找到一条增广路径    
48.                     if(h[1]>=m)      
49.                         return;    
50.     
51.                     //*c. 终止条件: 通过本次"find_path_sap(i)"递归搜索，找到了一条增广路径    
52.                     if(flag)    
53.                         break;    
54.     
55.                     //通过本次"find_path_sap(i)"递归没有找到一条增广路径，回溯时需要恢复状态    
56.                     cf_path=tmp_cf_path;     
57.                 }    
58.                 if(h[i]<minH)    
59.                     minH=h[i];    
60.             }    
61.         }    
62.         //以上for()执行完，minH为cur的邻接顶点中高度最小的点的高度    
63.     
64.         //3. 检查回溯结果    
65.         //   成功：逆着增广路径压入增广的流    
66.         //   失败：顶点cur高度增加1    
67.         if(flag){    
68.             //上面代码中的某次递归成功找到一条增广路径(从cur找下去，找到一条增广路径)    
69.             cf[cur][i]-=cf_path;    
70.             cf[i][cur]+=cf_path;    
71.         }else{    
72.             vh[h[cur]]--;    
73.             if(vh[h[cur]]==0)    
74.                 h[1]=m;     //作用到终止条件"*b"    
75.   
76.             h[cur]=minH+1;    
77.             vh[h[cur]]++;    
78.         }    
79.     }    
80.         
81.     //Ford-Fulkerson算法框架    
82.     int solve(){    
83.         vh[0]=m;    
84.         flow=0;  
85.   
86.         //h[1]保持<m，一旦增长到m，就不再存在任何增广路径  
87.         while(h[1]<m)    
88.         {    
89.             flag=false;  
90.             cf_path=INF;  
91.             find_path_sap(s);    
92.         }    
93.         return flow;    
94.     }    
95.         
96.     void addEdge(int x,int y,int c){    
97.             cf[x][y]+=c;    
98.     }    
99.         
100. };    
101.     
102. int main(){    
103.     int n;  //有向边数  
104.     int m;  //顶点数  
105.     while(scanf("%d %d",&n,&m)!=-1)    
106.     {    
107.         sap nt= sap();    
108.         nt.s=1;    
109.         nt.t=m;    
110.         nt.m=m;    
111.         for(int i=1;i<=n;i++)    
112.         {    
113.             int x,y,c;    
114.             scanf("%d %d %d",&x,&y,&c);    
115.             nt.addEdge(x,y,c);    
116.         }    
117.         printf("%d\n",nt.solve());    
118.     }    
119.     return 0;    
120. }   

C. SAP (shortest augmenting path)

    时间复杂度：O(|V|^2·|E|)      ==>每次搜索增广路径需要O(|V|)；可以增广O(|V||E|)次，原因见下面评价

    思路：每次"沿着从源到汇的允许弧（必然是源到汇的最短路径）"得到一条最短的增广路径增广。

             允许弧：如果一条弧(u,v)满足h[u]=h[v]+1 ，即前点仅比后点高1

    评价：经证明(independently proved by Dinic in 1970, and byEdmonds & Karp in 1972)，若每次从最短增广路径开始增广，则最多只用进行O(|V||E|)次增广，而不是DFS搜索增广路径的O(c)次。本质类似"BFS搜索增广路径"，只不过找最短路径的方法不同。再次注意，其中“最短增广路径”指的是：当每条边长度都为1时的最短路径

Learn more:

1. Dinic算法

2. 预流推进算法

Also check http://community.topcoder.com/tc?module=Static&d1=tutorials&d2=maxFlowRevisited