SAP算法详解与模板

链接:

1. Maximum Flow: Augmenting Path Algorithms Comparison

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

2. 刘汝佳《算法艺术与信息学竞赛》 P321 ( 注: 上面的代码似乎有误,retreat()部分未回退< 详见下文or 链接1. > )

---------------------------------------------

 

关键概念与性质:

距离函数(distance function),我们说一个距离函数是有效的当且仅当满足有效条件(valid function)

(1)d(t)= 0;

(2)d(i)<= d(j)+1(如果弧rij在残余网络G(x)中);

   

性质1:

如果距离标号是有效的,那么d(i)便是残余网络中从点i到汇点的距离的下限;

证明:

复制代码
令任意的从i结点到汇点t长度为k的路径,设为i= i1-i2-i3...ik+1=t  
d(i(k)) <= d(i(k+1))+1= d(t)+1=1  
d(i(k-1)) <= d(i(k))+1=2  
d(i(k-2)) <= d(i(k-1))+1=3  
            :  
            :  
d(i(1)) <= d(i(2))+1= k  
由此可见,d(i)是i到t的距离下限     
复制代码

 

 

  

性质2:

允许弧(边):如果对于边rij>0,且d(i)= d(j)+1,那么称为允许边

允许路:一条从源点到汇点的且有允许弧组成的路径

允许路是从源点到汇点的最短增广路径。

证明:

(1)因为rij>0所以必然是一条增广路

(2)假设路径p是一条允许路包含k条弧,那么由d(i) = d(j)+1 可知,d(s)= k;

又因为d(s)是点s到汇点的距离下限,所以距离下限为k,所以p便是一条最短路。

   

性质3:

在sap算法中距离标号始终是正确,有效的。并且每次的冲标号都会是距离标号严格递增

证明:略;

 

  

伪代码:

复制代码
//algorithm shortest augment path;  
void sap()  
    {  
        x =0;  
        obtain the exact distance labels d(i);  
        i = s;  
        while (d(s)<n)  
        {  
            if (i has an admissible arc)  
            {  
                advance(i);  
                if (i == t)  
                {  
                    augment;  
                    i = s;  
                }  
            }  
            else  
                retreat(i);  
        }  
    }  
      
    //procedure advance(i);  
void advance(i)  
    {  
        let(i,j)be an admissible arc in A(t);  
        pre(j) = i;  
        i = j;  
    }  
    //procedure retreat  
void retreat()  
    {  
        d(i) = min(d(j)+1):rij>0;  
        if (i != s)  
            i = pre(i);  
    }  
复制代码

 

  

代码模板:

复制代码
    #include <iostream>  
    #include <cstring>  
    #include <cstdlib>  
      
    usingnamespace std;  
      
    constint Max =225;  
    constint oo =210000000;  
      
    int n,m,c[Max][Max],r[Max][Max],c1[Max][Max],source,sink;  
      
    //c1是c的反向网络,用于dis数组的初始化  
int dis[Max];  
      
    void initialize()// 初始化dis数组  
    {  
         int q[Max],head =0, tail =0;//BFS队列  
         memset(dis,0,sizeof(dis));  
         q[++head] = sink;  
         while(tail < head)  
         {  
             int u = q[++tail], v;  
             for(v =0; v <= sink; v++)  
             {  
                   if(!dis[v] && c1[u][v] >0)  
                   {  
                       dis[v] = dis[u] +1;  
                       q[++head] = v;  
                   }  
             }  
         }  
    }  
      
    int maxflow_sap()  
    {  
        initialize();  
        int top = source, pre[Max], flow =0, i, j, k, low[Max];  
      
        // top是当前增广路中最前面一个点。  
        memset(low,0,sizeof(low));//low数组用于保存路径中的最小容量  
while(dis[source] < n)  
        {  
            bool flag =false;  
            low[source] = oo;  
            for(i =0; i <= sink; i++)//找允许弧,根据允许弧的定义  
            {  
                if(r[top][i] >0&& dis[top] == dis[i] +1)  
                {  
                    flag =true;  
                    break;  
                }  
            }  
            if(flag)// 如果找到允许弧  
            {  
                low[i] = r[top][i];  
                if(low[i] > low[top]) low[i] = low[top];//更新low  
                pre[i] = top; top = i;  
                if(top == sink)// 如果找到一条增广路了  
                {  
                    flow += low[sink];  
                    j = top;  
                    while(j != source)// 路径回溯更新残留网络  
                    {  
                        k = pre[j];  
                        r[k][j] -= low[sink];  
                        r[j][k] += low[sink];  
                        j = k;  
                    }  
                    top = source;//从头开始再找最短路  
                    memset(low,0,sizeof(low));  
                }  
            }  
            else// 如果没有允许弧  
            {  
                int mindis =10000000;  
                for(j =0; j <= sink; j++)//找和top相邻dis最小的点  
                {  
                    if(r[top][j] >0&& mindis > dis[j] +1)  
                        mindis = dis[j] +1;  
                }  
                dis[top] = mindis;//更新top的距离值  
if(top != source) top = pre[top];// 回溯找另外的路  
            }  
        }  
        return(flow);  
    }  
复制代码

 

 

  

运用gap优化:

即当标号中出现了不连续标号的情况时,即可以证明已经不存在新的增广流,此时的流量即为最大流。

简单证明下:

假设不存在标号为k的结点,那么这时候可以将所有的结点分成两部分,一部分为d(i)>k,另一部分为d(i)<k

如此分成两份,因为性质2可知,允许路为最短的增广路,又因为不存在从>k部分到<k部分的增广流,那么有最

大流最小割定理可知此时便是最大流。

  

优化代码:要注意在标号的时候不能直接把所有的初始为0,而应该为-1,否则会在gap优化的时候出现问题,不满足递增的性质,切记!

 

  

复制代码
    #include <iostream>  
    #include <cstring>  
    #include <cstdlib>  
      
    usingnamespace std;  
      
    constint Max =225;  
    constint oo =210000000;  
      
    int n,m,c[Max][Max],r[Max][Max],c1[Max][Max],source,sink;  
    int gap[Max];//用gap来记录当前标号为i的个数,用于gap优化;  
      
    //c1是c的反向网络,用于dis数组的初始化  
int dis[Max];  
      
    void initialize()// 初始化dis数组  
    {  
        memset(gap,0,sizeof(gap));  
        int q[Max],head =0, tail =0;//BFS队列  
        memset(dis,0xff,sizeof(dis));  
        q[++head] = sink;  
        dis[sink] =0;  
        gap[0]++;  
        while(tail < head)  
        {  
            int u = q[++tail], v;  
            for(v =0; v <= sink; v++)  
            {  
                if(!dis[v] && c1[u][v] >0)  
                {  
                    dis[v] = dis[u] +1;  
                    gap[dis[v]]++;  
                    q[++head] = v;  
                }  
            }  
         }  
    }  
      
    int maxflow_sap()  
    {  
        initialize();  
        int top = source, pre[Max], flow =0, i, j, k, low[Max];  
      
        // top是当前增广路中最前面一个点。  
        memset(low,0,sizeof(low));//low数组用于保存路径中的最小容量  
while(dis[source] < n)  
        {  
            bool flag =false;  
            low[source] = oo;  
            for(i =0; i <= sink; i++)//找允许弧,根据允许弧的定义  
            {  
                if(r[top][i] >0&& dis[top] == dis[i] +1&&dis[i]>=0)  
                {  
                    flag =true;  
                    break;  
                }  
            }  
            if(flag)// 如果找到允许弧  
            {  
                low[i] = r[top][i];  
                if(low[i] > low[top])   
                    low[i] = low[top];//更新low  
                pre[i] = top; top = i;  
                if(top == sink)// 如果找到一条增广路了  
                {  
                    flow += low[sink];  
                    j = top;  
                    while(j != source)// 路径回溯更新残留网络  
                    {  
                        k = pre[j];  
                        r[k][j] -= low[sink];  
                        r[j][k] += low[sink];  
                        j = k;  
                    }  
                    top = source;//从头开始再找最短路  
                    memset(low,0,sizeof(low));  
                }  
            }  
            else// 如果没有允许弧  
            {  
                int mindis =10000000;  
                for(j =0; j <= sink; j++)//找和top相邻dis最小的点  
                {  
                    if(r[top][j] >0&& mindis > dis[j] +1&&dis[j]>=0)  
                        mindis = dis[j] +1;  
                }  
                gap[dis[top]]--;  
                if (gap[dis[top]] ==0)  
                    break;  
                gap[mindis]++;  
                dis[top] = mindis;//更新top的距离值  
if(top != source) top = pre[top];// 回溯找另外的路  
            }  
        }  
        return(flow);  
    }  
复制代码

 

 

  

  

注意:在运用sap的时候必须要时刻的保证标号的两个性质,因此,不能再初始标号的时候将全部初始为0层,因为有些点是不具有层数的,或者说是层数是无穷大的,不可达的。

  

网上找的比较清晰地模板sap,有各种优化:

 

复制代码
    #include <stdio.h>  
    #include <string.h>  
    #define INF 2100000000  
    #define MAXN 301  
      
    int SAP(int map[][MAXN],int v_count,int s,int t)      //邻接矩阵,节点总数,始点,汇点  
    {  
        int i;  
        int cur_flow,max_flow,cur,min_label,temp;         //当前流,最大流,当前节点,最小标号,临时变量  
char flag;                                        //标志当前是否有可行流  
int cur_arc[MAXN],label[MAXN],neck[MAXN];         //当前弧,标号,瓶颈边的入点(姑且这么叫吧)  
int label_count[MAXN],back_up[MAXN],pre[MAXN];    //标号为i节点的数量,cur_flow的纪录,当前流路径中前驱  
      
        //初始化  
        memset(label,0,MAXN*sizeof(int));  
        memset(label_count,0,MAXN*sizeof(int));  
      
        memset(cur_arc,0,MAXN*sizeof(int));  
        label_count[0]=v_count;                           //全部初始化为距离为0  
      
        neck[s]=s;  
        max_flow=0;  
        cur=s;  
        cur_flow=INF;  
      
        //循环代替递归  
while(label[s]<v_count)  
        {  
            back_up[cur]=cur_flow;  
            flag=0;  
      
            //选择允许路径(此处还可用邻接表优化)  
for(i=cur_arc[cur];i<v_count;i++)    //当前弧优化  
            {  
               if(map[cur][i]!=0&&label[cur]==label[i]+1)    //找到允许路径  
               {  
                   flag=1;  
                   cur_arc[cur]=i;    //更新当前弧  
if(map[cur][i]<cur_flow)    //更新当前流  
                   {  
                       cur_flow=map[cur][i];  
                       neck[i]=cur;     //瓶颈为当前节点  
                   }  
                   else  
                   {  
                       neck[i]=neck[cur];     //瓶颈相对前驱节点不变  
                   }  
                   pre[i]=cur;    //记录前驱  
                   cur=i;  
                   if(i==t)    //找到可行流  
                   {  
                       max_flow+=cur_flow;    //更新最大流  
      
                       //修改残量网络  
while(cur!=s)  
                       {  
                           if(map[pre[cur]][cur]!=INF)map[pre[cur]][cur]-=cur_flow;  
                           back_up[cur] -= cur_flow;  
                           if(map[cur][pre[cur]]!=INF)map[cur][pre[cur]]+=cur_flow;  
                           cur=pre[cur];  
                       }  
      
                       //优化,瓶颈之后的节点出栈  
                       cur=neck[t];  
                       cur_flow=back_up[cur];   
                   }  
                   break;  
               }  
            }  
            if(flag)continue;  
      
            min_label=v_count-1;    //初始化min_label为节点总数-1  
      
            //找到相邻的标号最小的节点     
for(i=0;i<v_count;i++)  
            {  
                if(map[cur][i]!=0&&label[i]<min_label)  
                {  
                    min_label=label[i];  
                    temp=i;  
                }  
            }  
            cur_arc[cur]=temp;    //记录当前弧,下次从提供最小标号的节点开始搜索  
            label_count[label[cur]]--;    //修改标号纪录  
if(label_count[label[cur]]==0)break;    //GAP优化  
            label[cur]=min_label+1;    //修改当前节点标号  
            label_count[label[cur]]++;     //修改标号记录  
if(cur!=s)  
            {  
               //从栈中弹出一个节点  
               cur=pre[cur];  
               cur_flow=back_up[cur];  
            }  
        }  
        return(max_flow);  
    }  
复制代码

 

 

【转自】

http://blog.csdn.net/liguanxing/article/details/5783804

 


评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值