sap算法详解与模板

关键概念与性质：

距离函数（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);      }