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,r,c1,source,sink；
//c1是c的反向网络，用于dis数组的初始化
int dis;
void initialize()// 初始化dis数组
{
int q,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 && c1 >0)
{
dis = dis +1;
q[++head] = v;
}
}
}
}
int maxflow_sap()
{
initialize();
int top = source, pre, flow =0, i, j, k, low;
// top是当前增广路中最前面一个点。
memset(low,0,sizeof(low));//low数组用于保存路径中的最小容量
while(dis < n)
{
bool flag =false;
low = oo;
for(i =0; i <= sink; i++)//找允许弧，根据允许弧的定义
{
if(r >0&& dis == dis +1)
{
flag =true;
break;
}
}
if(flag)// 如果找到允许弧
{
low = r;
if(low > low) low = low;//更新low
pre = top; top = i;
if(top == sink)// 如果找到一条增广路了
{
flow += low;
j = top;
while(j != source)// 路径回溯更新残留网络
{
k = pre;
r -= low;
r += low;
j = k;
}
top = source;//从头开始再找最短路
memset(low,0,sizeof(low));
}
}
else// 如果没有允许弧
{
int mindis =10000000;
for(j =0; j <= sink; j++)//找和top相邻dis最小的点
{
if(r >0&& mindis > dis +1)
mindis = dis +1;
}
dis = mindis;//更新top的距离值
if(top != source) top = pre;// 回溯找另外的路
}
}
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,r,c1,source,sink;
int gap;//用gap来记录当前标号为i的个数，用于gap优化；
//c1是c的反向网络，用于dis数组的初始化
int dis;
void initialize()// 初始化dis数组
{
memset(gap,0,sizeof(gap));
int q,head =0, tail =0;//BFS队列
memset(dis,0xff,sizeof(dis));
q[++head] = sink;
dis =0;
gap++;
while(tail < head)
{
int u = q[++tail], v;
for(v =0; v <= sink; v++)
{
if(!dis && c1 >0)
{
dis = dis +1;
gap]++;
q[++head] = v;
}
}
}
}
int maxflow_sap()
{
initialize();
int top = source, pre, flow =0, i, j, k, low;
// top是当前增广路中最前面一个点。
memset(low,0,sizeof(low));//low数组用于保存路径中的最小容量
while(dis < n)
{
bool flag =false;
low = oo;
for(i =0; i <= sink; i++)//找允许弧，根据允许弧的定义
{
if(r >0&& dis == dis +1&&dis>=0)
{
flag =true;
break;
}
}
if(flag)// 如果找到允许弧
{
low = r;
if(low > low)
low = low;//更新low
pre = top; top = i;
if(top == sink)// 如果找到一条增广路了
{
flow += low;
j = top;
while(j != source)// 路径回溯更新残留网络
{
k = pre;
r -= low;
r += low;
j = k;
}
top = source;//从头开始再找最短路
memset(low,0,sizeof(low));
}
}
else// 如果没有允许弧
{
int mindis =10000000;
for(j =0; j <= sink; j++)//找和top相邻dis最小的点
{
if(r >0&& mindis > dis +1&&dis>=0)
mindis = dis +1;
}
gap]--;
if (gap] ==0)
break;
gap++;
dis = mindis;//更新top的距离值
if(top != source) top = pre;// 回溯找另外的路
}
}
return(flow);
}
注意：在运用sap的时候必须要时刻的保证标号的两个性质，因此，不能再初始标号的时候将全部初始为0层，因为有些点是不具有层数的，或者说是层数是无穷大的，不可达的。
网上找的比较清晰地模板sap，有各种优化：
#include <stdio.h>
#include <string.h>
#define INF 2100000000
#define MAXN 301
int SAP(int map[],int v_count,int s,int t) //邻接矩阵，节点总数，始点，汇点
{
int i;
int cur_flow,max_flow,cur,min_label,temp; //当前流，最大流，当前节点，最小标号，临时变量
char flag; //标志当前是否有可行流
int cur_arc,label,neck; //当前弧，标号，瓶颈边的入点（姑且这么叫吧）
int label_count,back_up,pre; //标号为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=v_count; //全部初始化为距离为0
neck=s;
max_flow=0;
cur=s;
cur_flow=INF;
//循环代替递归
while(label<v_count)
{
back_up=cur_flow;
flag=0;
//选择允许路径（此处还可用邻接表优化）
for(i=cur_arc;i<v_count;i++) //当前弧优化
{
if(map!=0&&label==label+1) //找到允许路径
{
flag=1;
cur_arc=i; //更新当前弧
if(map<cur_flow) //更新当前流
{
cur_flow=map;
neck=cur; //瓶颈为当前节点
}
else
{
neck=neck; //瓶颈相对前驱节点不变
}
pre=cur; //记录前驱
cur=i;
if(i==t) //找到可行流
{
max_flow+=cur_flow; //更新最大流
//修改残量网络
while(cur!=s)
{
if(map]!=INF)map]-=cur_flow;
back_up -= cur_flow;
if(map]!=INF)map]+=cur_flow;
cur=pre;
}
//优化，瓶颈之后的节点出栈
cur=neck;
cur_flow=back_up;
}
break;
}
}
if(flag)continue;
min_label=v_count-1; //初始化min_label为节点总数-1
//找到相邻的标号最小的节点
for(i=0;i<v_count;i++)
{
if(map!=0&&label<min_label)
{
min_label=label;
temp=i;
}
}
cur_arc=temp; //记录当前弧，下次从提供最小标号的节点开始搜索
label_count]--; //修改标号纪录
if(label_count]==0)break; //GAP优化
label=min_label+1; //修改当前节点标号
label_count]++; //修改标号记录
if(cur!=s)
{
//从栈中弹出一个节点
cur=pre;
cur_flow=back_up;
}
}
return(max_flow);
}
【转自】
http://blog.iyunv.com/liguanxing/article/details/5783804
页:
[1]