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

原文地址:https://www.cnblogs.com/longdouhzt/p/2166187.html