注册 登录  
 加关注
查看详情
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

z55250825

一只蒟蒻

 
 
 

日志

 
 

【算法】【网络流】SAP  

2013-03-23 17:22:29|  分类: 算法 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

注明这是个转载...来自百度百科的SAP,只是想方便查找...

SAP算法

 

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)的时间内复原路径中瓶颈边之前的所有边。

 

 

代码:

var

  flag:boolean;

  jl,min,flow,aug,j,m,n,tmp,a,b,c,i:longint;

  his,pre,dis,vh,di:array[0..1024] of longint;

  map:array[1..1024,1..1024] of longint;

begin

  assign(input,'ditch.in');reset(input);

  assign(output,'ditch.out');rewrite(output);

  readln(m,n);

  for i:=1 to m do

   begin

     readln(a,b,c);

     inc(map[a,b],c);

   end;

 

vh[0]:=n;

 

for i:=1 to n do di[i]:=1;{当前弧初始为第一条弧}

 

i:=1;{S开始搜}

aug:=maxlongint;

while dis[1]<n do

 begin

  his[i]:=aug;{标记,以便以后返回这个值}

  flag:=false;{这个是是否有找到允许弧的标志}

  for j:=di[i] to n do{从当前弧开始搜}

     begin

       if (map[i,j]>0)and(dis[j]+1=dis[i]) then{找到允许弧}

       begin

         flag:=true;

         di[i]:=j;{标记当前弧}

         if map[i,j]<aug then aug:=map[i,j];

         pre[j]:=i;{记录前驱}

         i:=j;

         if i=n then{找到增广路}

           begin

            inc(flow,aug);

             while i<>1 do

               begin

                tmp:=i;

                i:=pre[i];

                dec(map[i,tmp],aug);

                inc(map[tmp,i],aug);

              end;

            aug:=maxlongint;

          end;

        break;{找到允许弧则退出查找}

      end;

    end;

 if flag then continue;{找到允许弧}

 min:=n-1;{没有允许弧了,需要重标号}

 for j:=1 to n do

   begin

     if (map[i,j]>0)and(dis[j]<min) then

          begin

             jl:=j;

             min:=dis[j];

          end;

   end;

 di[i]:=jl;

 dec(vh[dis[i]]);{GAP优化}

 if vh[dis[i]]=0 then break;

 dis[i]:=min+1;

 inc(vh[dis[i]]);

 if i<>1 then

    begin

        i:=pre[i];{返回上一层}

       aug:=his[i];{知道之前记录这个值的用处了吧}

    end;

 end;

 write(flow);

 close(input);

 close(output);

end.

 

优化:

 

1.邻接表优化:

如果顶点多的话,往往N^2存不下,这时候就要存边:

存每条边的出发点,终止点和价值,然后排序一下,再记录每个出发点的位置。以后要调用从出发点出发的边时候,只需要从记录的位置开始找即可(其实可以用链表)。优点是时间加快空间节省,缺点是编程复杂度将变大,所以在题目允许的情况下,建议使用邻接矩阵。

 

2.GAP优化:

如果一次重标号时,出现距离断层,则可以证明ST无可行流,此时则可以直接退出算法。

 

3.当前弧优化:

为了使每次找增广路的时间变成均摊O(V),还有一个重要的优化是对于每个点保存“当前弧”:初始时当前弧是邻接表的第一条弧;在邻接表中查找时从当前弧开始查找,找到了一条允许弧,就把这条弧设为当前弧;改变距离标号时,把当前弧重新设为邻接表的第一条弧。

 

 

sap算法心得

 

网络最大流算法是网络流算法的基础,实现方法很多,但时间复杂度与编程复杂

度难于兼顾。一方面,诸如预流推进等高效算法难于编写调试,有时实际效果欠佳

(参见dd_engi的评测);另一方面,基于增广路的算法大多时间效率不高。于是许多人选择了相对简单的Dinic算法。事实上,SAP算法更易于理解,时间效率更高,编程更简单,思路更清晰。

 

(名词解释)SAPShortest Augmenting Paths: 最短增广路

 

算法基本框架:

· 定义距离标号为到汇点距离的下界。

· 在初始距离标号的基础上,不断沿可行弧找增广路增广,一般采用深度优先。可行弧定义为:{(i,j)|h[i]=h[j]+1}

· 遍历当前节点完成后,为了使下次再来的时候有路可走(当然要满足距离标号的性质:不超过真实距离),重标号当前节点为:min{h[j]|(i,j)}+1

· 重标号后当前节点处理完毕。当源点被重标号后,检查其距离标号,当大于顶点数时,图中已不存在可增广路,此时算法结束;否则再次从源点遍历。

· 理论复杂度为:O(n2m)

 

我的心得:

· 理论上初始标号要用反向BFS求得,实践中可以全部设为0,可以证明:这样做

不改变渐进时间复杂度。

· 理论上可以写出许多子程序并迭代实现,但判断琐碎,没有层次,实践中用递归简单明了,增加的常数复杂度比起SAP速度微乎其微却极大降低编程复杂度,易

于编写调试。

· GAP优化★ (性价比极高,推荐!详见程序“//GAP”部分)

注意到我们在某次增广后,最大流可能已经求出,因此算法做了许多无用功。可

以发现,距离标号是单调增的。这启示我们如果标号中存在“间隙”,则图中不

会再有可增广路,于是算法提前终止。实践中我们使用数组vh[i]记录标号为i

顶点个数,若重标号使得vh中原标号项变为0,则停止算法。

 

附效率测试与比较:

算法高标推进SAP(GAP) Dinic SAP(NOGAP) BFS+EK

时间(s) 5.67 6.02 11.35 23.90 55.77

 

:高标推进程序在180行左右。Dinic用的是dd_engiBFS的标程,带GAP优化的SAP参见程序。

 

附源程序与注释(注意以下程序未用BFS,只用到递归,实现简单):

 

program ditch; //USACOditch为例

var

c:array[0..1000,0..1000] of cardinal; //邻接矩阵

h,vh:array[0..1000] of cardinal;

n,augc,rd1,rd2,m,i:cardinal; //augc为增广路容量

flow:cardinal=0;

found:boolean; //记录是否已达汇点

procedure aug(const m:cardinal);

var

i,augco,minh:cardinal;

begin

minh:=n-1;

augco:=augc;

if m=n then begin

found:=true;

inc(flow,augc);

exit;

end;

for i:=1 to n do

if c[m,i]>0 then begin

if h[i]+1=h[m] then begin

if c[m,i]<augc then augc:=c[m,i];

aug(i);

if h[1]>=n then exit; //GAP

if found then break;

augc:=augco;

end;

if h[i]<minh then minh:=h[i];

end;

if not found then begin //重标号

dec(vh[h[m]]); //GAP

if vh[h[m]]=0 then h[1]:=n; //GAP

h[m]:=minh+1;

inc(vh[h[m]]); //GAP

end else begin //修改残量

dec(c[m,i],augc);

inc(c[i,m],augc);

end;

end;

begin

assign(input,'ditch.in');

assign(output,'ditch.out');

reset(input);

rewrite(output);

readln(m,n);

fillchar(c,sizeof(c),0);

for i:=1 to m do begin

readln(rd1,rd2,augc);

inc(c[rd1,rd2],augc);

end;

fillchar(h,sizeof(h),0);

fillchar(vh,sizeof(vh),0);

vh[0]:=n;

while h[1]<n do begin

augc:=$FFFFFFFF;

found:=false;

aug(1);

end;

writeln(flow);

close(input);

close(output);

end.

 

Reference:

 

1.几种高级网络流算法的实测, dd_engi

2.Network Flows: Theory, Algorithms, and Applications, Ahuja, R.K.

3.http://www.gnocuil.cn/blog/2007/11/网络流的sap算法, gnocuil.

 

 

c++语言版本:(pku 1273

 

#include<iostream>

using namespace std;

 

#define maxn 9999999

#define N 110

 

int n,augc,m,i,flow;

int h[N],vh[N];

int c[N][N];

bool found;

 

void aug(int m)

{

int i,augco,minh;

minh=n-1;

augco=augc;

if(n==m) { found=true; flow+=augc; return; }

for(i=1;i<=n;i++)

{

   if(c[m][i]>0){

    if(h[i]+1==h[m]){

     if(c[m][i]<augc) augc=c[m][i];

     aug(i);

     if(h[1]>=n) return;//GAP

     if(found) break;

     augc=augco;

    }

    if(h[i]<minh) minh=h[i];

   }

}

if(!found){

   vh[h[m]]--;

   if(vh[h[m]]==0) h[1]=n;

   h[m]=minh+1;

   vh[h[m]]++;

}

else { c[m][i]-=augc;c[i][m]+=augc; }

}

 

int main()

{

int x,y,z;

while(scanf("%d%d",&m,&n)!=EOF)

{

   memset(c,0,sizeof(c));

   for(i=1;i<=m;i++)

   {

    scanf("%d%d%d",&x,&y,&z);

    c[x][y]+=z;

   }

   memset(h,0,sizeof(h));

   memset(vh,0,sizeof(vh));

   vh[0]=n;flow=0;

   while(h[1]<n)

   {

    augc=maxn;

    found=false;

    aug(1);

   }

   printf("%d\n",flow);

}

return 0;

}

 

## bfs改进版 pascal ##

 

program profit;

const

name='profit';

type

toetype=^etype;

etype=record

t,c:longint;

p,next:toetype;

end;

 

var

e:array[-55000..5500] of toetype;

h:array[-55000..5500] of longint;

vh:array[0..70000] of longint;

n,m,i,rd,rd1,rd2:longint;

per:longint=0;

 

procedure addedge(const sn,tn,cn:longint);

var

newe:toetype;

begin

new(newe);

with newe^ do begin

t:=tn;

c:=cn;

next:=e[sn];

new(p);

with p^ do begin

   t:=sn;

   c:=0;

   next:=e[tn];

   p:=newe;

end;

e[tn]:=p;

end;

e[sn]:=newe;

end;

 

function aug(const u,t:longint):longint;

var

minh,l,d:longint;

ie:toetype;

begin

if u=0 then begin aug:=t;exit;end;

l:=t;

minh:=n+m+1;

ie:=e[u];

while ie<>nil do begin

if ie^.c>0 then begin

   if h[ie^.t]+1=h[u] then begin

    if l<ie^.c then d:=l else d:=ie^.c;

    d:=aug(ie^.t,d);

    dec(ie^.c,d);

    inc(ie^.p^.c,d);

    dec(l,d);

    if h[n+1]>=n+m+2 then begin aug:=t-l;exit;end;

    if l=0 then break;

   end;

   if h[ie^.t]<minh then minh:=h[ie^.t];

end;

ie:=ie^.next;

end;

if l=t then begin

dec(vh[h[u]]);

if vh[h[u]]=0 then h[n+1]:=n+m+2;

h[u]:=minh+1;

inc(vh[h[u]]);

end;

aug:=t-l;

end;

 

begin

assign(input,name+'.in');

assign(output,name+'.out');

reset(input);

rewrite(output);

 

readln(n,m);

fillchar(e,sizeof(e),0);

for i:=1 to n do begin

read(rd);

addedge(n+1,i,rd);

end;

readln;

for i:=1 to m do begin

readln(rd1,rd2,rd);

inc(per,rd);

addedge(rd1,-i,maxlongint);

addedge(rd2,-i,maxlongint);

addedge(-i,0,rd);

end;

fillchar(h,sizeof(h),0);

fillchar(vh,sizeof(vh),0);

vh[0]:=n+m+2;

while h[n+1]<n+m+2 do

dec(per,aug(n+1,maxlongint));

writeln(per);

 

close(input);

close(output);

end.

 

==这个是SAP==

 

program ditch;

Const Maxn = 3000;

Inf=MaxLongint;

Fin ='Ditch.in';

Fout = 'Ditch.out';

Var Cf,g:Array[0..Maxn,0..Maxn]Of Longint;

    d,pre,list,Min,Cur,Numb:Array[0..Maxn]Of Longint;

    n,m,MaxFlow:Longint;

 

Procedure Init;

Var i,a,b,c,Open,Closed:Longint;

Begin

    Assign(Input,Fin);Reset(Input);

    Readln(m,n);

    For i:=1 To m Do Begin

      Readln(a,b,c);

      Inc(Cf[a,b],c);

      Inc(G[a,0]);G[a,G[a,0]]:=b;

      Inc(G[b,0]);G[b,G[b,0]]:=a;

    End;

 

    Fillchar(D,Sizeof(D),$7f);

    D[n]:=0;List[1]:=n;Open:=1;Closed:=0;

    While Open>Closed Do Begin

      Inc(Closed);a:=List[Closed];

      For i:=1 To G[a,0] Do

        If D[a]+1<D[G[a,i]] Then Begin

          Inc(Open);List[Open]:=G[a,i];D[G[a,i]]:=D[a]+1;

        End;

    End;

    For i:=1 To n Do Inc(Numb[D[i]]);

    Close(Input);

End;

 

Procedure Modify;

Var i:Longint;

Begin

    i:=n;

    While Pre[i]<>0 Do Begin

      Inc(Cf[i,Pre[i]],Min[n]);

      Dec(Cf[Pre[i],i],Min[n]);

      i:=Pre[i];

    End;

    Inc(MaxFlow,Min[n]);

End;

 

Procedure Print;

Begin

    Assign(Output,Fout);Rewrite(Output);

    Writeln(MaxFlow);

    Close(Output);

    Halt;

End;

 

Procedure Main;

Var i,j,k:Longint;

      Flag:Boolean;

Begin

    i:=1;MaxFlow:=0;Min[1]:=Inf;

 

    While D[1]<n Do Begin

      Flag:=False;k:=Inf;

      While Cur[i]<=G[i,0] Do Begin  

        If (Cf[i,G[i,cur[i]]] > 0)And(D[G[i,cur[i]]]+1=D[i]) Then Begin

          Flag:=True;

          Pre[G[i,cur[i]]]:=i;

          If Cf[i,G[i,cur[i]]]<Min[i] Then

            Min[G[i,cur[i]]]:=Cf[i,G[i,cur[i]]]

          Else Min[G[i,cur[i]]]:=Min[i];

          i:=G[i,cur[i]];

          If i=n Then Begin

            Modify;

            i:=1;

          End;

        End

        Else Inc(Cur[i]);

      End;

 

      If Not Flag Then Begin

        Dec(Numb[D[i]]);

        If Numb[D[i]]=0 Then Print;

        For j:=1 To G[i,0] Do If (Cf[i,G[i,j]]>0) And (D[G[i,j]]<k) Then

          k:=D[G[i,j]];

        If k<Inf Then

          D[i]:=k+1

        Else

          D[i]:=n;

        Cur[i]:=1;

        Inc(Numb[D[i]]);

        If i<>1 Then i:=Pre[i];

      End;

    End;

End;

 

begin

Init;

Main;

Print;

end.

 

## c++ 语言版 sap ##

 

//sap 网络最大流算法

#include<iostream>

using namespace std;

 

const int N = 3000;

const int inf=0x7fffffff;

 

int cf[N][N],g[N][N];

int d[N],pre[N],list[N],minn[N],cur[N],numb[N];

int n,m,maxflow,i,j,x,y,z,open,closed;

 

void init()

{

memset(cf,0,sizeof(cf));

memset(g,0,sizeof(g));

memset(d,inf,sizeof(d));

memset(list,0,sizeof(list));

memset(pre,0,sizeof(pre));

memset(numb,0,sizeof(numb));

memset(minn,0,sizeof(minn));

memset(cur,0,sizeof(cur));

for(i=1;i<=m;i++)

{

   scanf("%d%d%d",&x,&y,&z);

   cf[x][y]+=z;

   g[x][0]++;g[x][g[x][0]]=y;

   g[y][0]++;g[y][g[y][0]]=x;

}

 

d[n]=0;list[1]=n;open=1;closed=0;

while(open>closed)

{

   closed++;x=list[closed];

   for(i=1;i<=g[x][0];i++)

    if(d[x]+1<d[g[x][i]]){

     open++;

     list[open]=g[x][i];

     d[g[x][i]]=d[x]+1;

    }

}

for(i=1;i<=n;i++) numb[d[i]]++;

}

 

void modify()

{

int i=n;

while(pre[i]!=0)

{

   cf[i][pre[i]]+=minn[n];

   cf[pre[i]][i]-=minn[n];

   i=pre[i];

}

maxflow+=minn[n];

}

 

void sap()

{

int i,j,k;bool flag;

i=1;maxflow=0;minn[1]=inf;

 

while(d[1]<n)

{

   flag=false;k=inf;

   while(cur[i]<=g[i][0]) {

    if(cf[i][g[i][cur[i]]]>0 && d[g[i][cur[i]]]+1==d[i]){

     flag=true;

     pre[g[i][cur[i]]]=i;

     if(cf[i][g[i][cur[i]]]<minn[i]) minn[g[i][cur[i]]]=cf[i][g[i][cur[i]]];

     else minn[g[i][cur[i]]]=minn[i];

     i=g[i][cur[i]];

     if(i==n) { modify(); i=1; }

    }

    else cur[i]++;

   }

 

   if(!flag){

    numb[d[i]]--;

    if(numb[d[i]]==0) return;

    for(j=1;j<=g[i][0];j++)

     if(cf[i][g[i][j]]>0 && d[g[i][j]]<k) k=d[g[i][j]];

    if(k<inf) d[i]=k+1; else d[i]=n;

    cur[i]=1;

    numb[d[i]]++;

    if(i!=1) i=pre[i]; 

   }

}

}

 

int main()

{

while(scanf("%d%d",&m,&n)!=EOF)

{

   init();

   sap();

   printf("%d\n",maxflow);

}

return 0;

}

  评论这张
 
阅读(124)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2018