最大流（二）——SAP算法

找了半天，“百度百科”上这篇文章还能比较直观地帮我理解SAP算法:
 
求最大流有一种经典的算法，就是每次找增广路时用BFS找，保证找到的增广路是弧数最少的(即边权都为1时的最短路径)，也就是所谓的Edmonds-Karp算法。可以证明的是在使用最短路增广时增广过程不超过|V|*|E|次，每次BFS的时间都是O(|E|)，所以Edmonds-Karp的时间复杂度就是O(|V|*|E|^2)。
 
如果能让每次寻找增广路时的时间复杂度降下来，那么就能提高算法效率了，使用距离标号d（有时候称为“高度h”）的最短增广路算法就是这样的。所谓距离标号 ，就是某个点到汇点的最少的弧的数量（即边权值为1时某个点到汇点的最短路径长度） 。设点i的标号为d[i]，那么如果将满足d[i]=d[j]+1的弧(i,j)叫做允许弧 ，且增广时只走允许弧，那么就可以达到“怎么走都是最短路”的效果 。每个点的初始标号可以在一开始用一次从汇点沿所有反向边的BFS求出，实践中可以初始设全部点的距离标号为0，问题就是如何在增广过程中维护这个距离标号。

 

 

Sam评注：
(1).
    外国大牛的文章http://community.topcoder.com/tc?module=Static&d1=tutorials&d2=maxFlowRevisited#1 中有这样一句话非常重要：“We can distance function exactly if each i in V, d[i] equals the length of the shortest path from i to t(sink) in the residual network.”
    上面这句话意思是说：只有当d[i]恰好等于i到汇点t的最短路径（边权都为1）时，才称d为“距离函数”。而每次用SAP寻找增广路径初始，d[i] (i∈V)都被初始化为0。因此一开始，只有汇点的距离函数d[t]=0满足这个条件；对其他顶点，d[i]只能称为“i的距离函数值的下界/估计”，因为它<真正的d[i]。于是，外国大牛又冒出这样一句话：
“It is also easy to prove that if d(s)≥|V|, then the residual network contains no path from the source to the sink.”
(2).
    寻找一条最短（边最少）增广路径中用到的“递归机制”，实际上保证了最短路径的层次图是从汇点向源点扩展的。而每次用SAP寻找一条增广路径开始时，d[t]=0已经是真正的最短路径，因此如果有最短增广路径，完成SAP搜索后d[i]=真正的最短路径。

void  find_path_sap(int  cur){
           ...
           find_path_sap(i);
           ...
}

 

维护距离标号的方法是这样的：当找增广路过程中发现某点出发没有允许弧时，将这个点的距离标号设为由它出发的所有弧的终点的距离标号的最小值加一。这种维护距离标号的方法的正确性我就不证了。由于距离标号的存在，由于“怎么走都是最短路”，所以就可以采用DFS找增广路，用一个栈保存当前路径的弧即可。当某个点的距离标号被改变时，栈中指向它的那条弧肯定已经不是允许弧了，所以就让它出栈，并继续用栈顶的弧的端点增广。为了使每次找增广路的时间变成均摊O(V)，还有一个重要的优化是对于每个点保存“当前弧”：初始时当前弧是邻接表的第一条弧；在邻接表中查找时从当前弧开始查找，找到了一条允许弧，就把这条弧设为当前弧；改变距离标号时，把当前弧重新设为邻接表的第一条弧，还有一种在常数上有所优化的写法是改变距离标号时把当前弧设为那条提供了最小标号的弧。当前弧的写法之所以正确就在于任何时候我们都能保证在邻接表中当前弧的前面肯定不存在允许弧。
 
还有一个常数优化是在每次找到路径并增广完毕之后不要将路径中所有的顶点退栈，而是只将瓶颈边以及之后的边退栈，这是借鉴了Dinic算法的思想。注意任何时候待增广的“当前点”都应该是栈顶的点的终点。这的确只是一个常数优化，由于当前边结构的存在，我们肯定可以在O(n)的时间内复原路径中瓶颈边之前的所有边。
 
几个重要优化：
 
1.邻接表优化：
如果顶点多的话，往往N^2存不下，这时候就要存边：
存每条边的出发点，终止点和价值，然后排序一下，再记录每个出发点的位置。以后要调用从出发点出发的边时候，只需要从记录的位置开始找即可（其实可以用链表）。优点是时间加快空间节省，缺点是编程复杂度将变大，所以在题目允许的情况下，建议使用邻接矩阵。
 
2.GAP优化：
如果一次重标号时，出现距离断层，则可以证明ST无可行流，此时则可以直接退出算法。
 
3.当前弧优化：
为了使每次找增广路的时间变成均摊O(V)，还有一个重要的优化是对于每个点保存“当前弧”：初始时当前弧是邻接表的第一条弧；在邻接表中查找时从当前弧开始查找，找到了一条允许弧，就把这条弧设为当前弧；改变距离标号时，把当前弧重新设为邻接表的第一条弧。
 
学过之后又看了算法速度的比较，发现如果写好的话SAP的速度不会输给HLPP。

 

在前面的SAP程序加入一些注释，观察增广过程：

/*
测试用例： 
6 5
1 2 2
1 3 5
2 4 1
3 4 6
3 5 1
4 5 5
*/
#include<stdio.h>  
#include<string.h>  
#include <iostream>
using namespace std;

#define INF 1<<29;  
  
struct sap  
{  
    int s,t;    //s:源点 t:汇点  
    int m;      //顶点数  
    int cf[305][305];   //残留容量  
  
    int flow;       //当前最大流  
    int cf_path;    //本次可增广的流量（本次增广路径的残留容量）
    bool flag;      //本次增广路径是否找到  
    int vh[1000];   //各高度的顶点个数  
    int h[1000];    //各顶点的高度，又称为“距离标号”——到汇点距离的下界（边权值都为1）
  
    sap(){  
        flow =cf_path =s =t =0;  
        flag=false;  
        memset(cf,0,sizeof(cf));  
        memset(vh,0,sizeof(vh));  
        memset(h,0,sizeof(h));  
    }  
  
    /*find_path_sap(int cur):增广路径已经推进到顶点cur，尝试向下一个顶点推进
       主要思路：
             递归地尝试沿着当前点cur的所有允许弧推进增广路径。
             如果失败，将当前点cur的高度（距离函数下界）更新到minH(minH:残留网络中所有cur指向节点的最小高度)
    */
    void find_path_sap(int cur){  
        //1. 检查是否走到增广路径末端  
        //*a. 终止条件：cur是本次增广路径末端节点(汇点)  
        if(cur==t){  
            flow+=cf_path;  
            flag=true;  
            return;  
        }  
  
        //2. 尝试递归/回溯 找增广路径  
        //如果cur仅比邻接点高1，则尝试最大限度向它送出流  
        int i;  
        int minH=m-1;  
        int tmp_cf_path=cf_path;  
        for(i=1;i<=m;i++){  
            if(cf[cur][i]){  
                if(h[i]+1==h[cur]){  
                    if(cf[cur][i]<cf_path)  
                        cf_path=cf[cur][i];  
                    find_path_sap(i);  
                      
                    //*b. 终止条件: 如果h[1]增长到m，表示通过本次"find_path_sap(i)"递归搜索，不能找到一条增广路径  
                    if(h[1]>=m)    
                        return;  
  
                    //*c. 终止条件: 通过本次"find_path_sap(i)"递归搜索，找到了一条增广路径  
                    if(flag)  
                        break;  
  
                    //通过本次"find_path_sap(i)"递归没有找到一条增广路径，回溯时需要恢复状态  
                    cf_path=tmp_cf_path;   
                }  
                if(h[i]<minH)  
                    minH=h[i];  
            }  
        }  
        //以上for()执行完，minH为cur的邻接顶点中高度最小的点的高度  
  
        //3. 检查回溯结果  
        //   成功：逆着增广路径压入增广的流  
        //   失败：顶点cur高度增加1  
        if(flag){  
            //上面代码中的某次递归成功找到一条增广路径(从cur找下去，找到一条增广路径)  
            cf[cur][i]-=cf_path;  
            cf[i][cur]+=cf_path;  
        }else{  
            vh[h[cur]]--;  
            if(vh[h[cur]]==0)  
                h[1]=m;     //作用到终止条件"*b"  

            h[cur]=minH+1;  
            vh[h[cur]]++;  
        }  
    }  
      
    //Ford-Fulkerson算法框架  
    int solve(){  
        vh[0]=m;  
        flow=0;

        int cnt=1;
		//h[1]保持<m，一旦增长到m，就不再存在任何增广路径
        while(h[1]<m)  
        {  
            flag=false;
            cf_path=INF;
            find_path_sap(s);
            
            cout<<"第"<<cnt++<<"次: "<<flow<<endl;
            for(int i=1;i<=m;i++){
                for(int j=1;j<=m;j++){
                    if(cf[i][j]>0){
                        cout<<"cf["<<i<<","<<j<<"]= "<<cf[i][j]<<endl;
                    }
                }
            }
            
            for(int i=1;i<=m;i++){
                cout<<"h["<<i<<"]"<<"\t";
            }
            cout<<endl;
            for(int i=1;i<=m;i++){
                cout<<h[i]<<"\t";
            }
            cout<<endl;
            cout<<"-----------------------"<<endl;
        }  
        return flow;  
    }  
      
    void addEdge(int x,int y,int c){  
            cf[x][y]+=c;  
    }  
      
};  
  
int main(){  
    int n;	//有向边数
	int m;	//顶点数
    while(scanf("%d %d",&n,&m)!=-1)  
    {  
        sap nt= sap();  
        nt.s=1;  
        nt.t=m;  
        nt.m=m;  
        for(int i=1;i<=n;i++)  
        {  
            int x,y,c;  
            scanf("%d %d %d",&x,&y,&c);  
            nt.addEdge(x,y,c);  
        }  
        printf("%d\n",nt.solve());  
    }  
    return 0;  
}