zoukankan      html  css  js  c++  java
  • 网络流-最大流问题 ISAP 算法解释(转自Renfei Song's Blog)

    网络流-最大流问题 ISAP 算法解释

    ISAP 是图论求最大流的算法之一,它很好的平衡了运行时间和程序复杂度之间的关系,因此非常常用。

    约定

    我们使用邻接表来表示图,表示方法可以见文章带权最短路 Dijkstra, SPFA, Bellman-Ford, ASP, Floyd-Warshall 算法分析二分图的最大匹配、完美匹配和匈牙利算法的开头(就不重复贴代码了)。在下文中,图的源点(source)表示为 s

    ,汇点(sink)表示为 t ,当前节点为 u 。建图时,需要建立双向边(设反向的边容量为 0

    )才能保证算法正确。

    引入

    求解最大流问题的一个比较容易想到的方法就是,每次在残量网络(residual network)中任意寻找一条从 s

    t 的路径,然后增广,直到不存在这样的路径为止。这就是一般增广路算法(labeling algorithm)。可以证明这种不加改进的贪婪算法是正确的。假设最大流是 f ,那么它的运行时间为O( fE) 。但是,这个运行时间并不好,因为它和最大流 f

    有关。

    人们发现,如果每次都沿着残量网络中的最短增广路增广,则运行时间可以减为 O(E2V) 

    。这就是最短增广路算法。而 ISAP 算法则是最短增广路算法的一个改进。其实,ISAP 的意思正是「改进的最短增广路」 (Improved Shortest Augmenting Path)。

    顺便说一句,上面讨论的所有算法根本上都属于增广路方法(Ford-Fulkerson method)。和它对应的就是大名鼎鼎的预流推进方法(Preflow-push method)。其中最高标号预流推进算法(Highest-label preflow-push algorithm)的复杂度可以达到 O(V2E)

    。虽然在复杂度上比增广路方法进步很多,但是预流推进算法复杂度的上界是比较紧的,因此有时差距并不会很大。

    算法解释

    概括地说,ISAP 算法就是不停地找最短增广路,找到之后增广;如果遇到死路就 retreat,直到发现s

    , t不连通,算法结束。找最短路本质上就是无权最短路径问题,因此采用 BFS 的思想。具体来说,使用一个数组d,记录每个节点到汇点t的最短距离。搜索的时候,只沿着满足d[u]=d[v]+1的边uv

    (这样的边称为允许弧)走。显然,这样走出来的一定是最短路。

    原图存在两种子图,一个是残量网络,一个是允许弧组成的图。残量网络保证可增广,允许弧保证最短路(时间界较优)。所以,在寻找增广路的过程中,一直是在残量网络中沿着允许弧寻找。因此,允许弧应该是属于残量网络的,而非原图的。换句话说,我们沿着允许弧,走的是残量网络(而非原图)中的最短路径。当我们找到沿着残量网络找到一条增广路,增广后,残量网络肯定会变化(至少少了一条边),因此决定允许弧的d

    数组要进行相应的更新(顺便提一句,Dinic 的做法就是每次增广都重新计算d数组)。然而,ISAP 「改进」的地方之一就是,其实没有必要马上更新d

    数组。这是因为,去掉一条边只可能令路径变得更长,而如果增广之前的残量网络存在另一条最短路,并且在增广后的残量网络中仍存在,那么这条路径毫无疑问是最短的。所以,ISAP 的做法是继续增广,直到遇到死路,才执行 retreat 操作。

    说到这里,大家应该都猜到了,retreat 操作的主要任务就是更新d

    数组。那么怎么更新呢?非常简单:假设是从节点u找遍了邻接边也没找到允许弧的;再设一变量m,令m等于残量网络中u的所有邻接点的d数组的最小值,然后令d[u]等于m+1即可。这是因为,进入 retreat 环节说明残量网络中ut已经不能通过(已过时)的允许弧相连,那么ut实际上在残量网络中的最短路的长是多少呢?(这正是d的定义!)显然是残量网络中u的所有邻接点和t的距离加1的最小情况。特殊情况是,残量网络中u根本没有邻接点。如果是这样,只需要把d[u]设为一个比较大的数即可,这会导致任何点到u的边被排除到残量网络以外。(严格来说只要大于等于V即可。由于最短路一定是无环的,因此任意路径长最大是V1)。修改之后,只需要把正在研究的节点u

    沿着刚才走的路退一步,然后继续搜索即可。

    讲到这里,ISAP 算法的框架内容就讲完了。对于代码本身,还有几个优化和实现的技巧需要说明。

    1. 算法执行之前需要用 BFS 初始化d

    数组,方法是从ts

    • 逆向进行。
    • 算法主体需要维护一个「当前节点」u
    • ,执行这个节点的前进、retreat 等操作。
    • 记录路径的方法非常简单,声明一个数组p,令p[i]等于增广路上到达节点i的边的序号(这样就可以找到从哪个顶点到的顶点i
    • )。需要路径的时候反向追踪一下就可以了。
    • 判断残量网络中s,t不连通的条件,就是d[s]V 。这是因为当s,t不连通时,最终残量网络中s将没有任何邻接点,对s
    • 的 retreat 将导致上面条件的成立。
    • GAP 优化。GAP 优化可以提前结束程序,很多时候提速非常明显(高达 100 倍以上)。GAP 优化是说,进入 retreat 环节后,u,t之间的连通性消失,但如果u是最后一个和t距离d[u](更新前)的点,说明此时s,t也不连通了。这是因为,虽然u,t已经不连通,但毕竟我们走的是最短路,其他点此时到t的距离一定大于d[u](更新前),因此其他点要到t,必然要经过一个和t距离为d[u]
      1. (更新前)的点。GAP 优化的实现非常简单,用一个数组记录并在适当的时候判断、跳出循环就可以了。
      2. 另一个优化,就是用一个数组保存一个点已经尝试过了哪个邻接边。寻找增广的过程实际上类似于一个 BFS 过程,因此之前处理过的邻接边是不需要重新处理的(残量网络中的边只会越来越少)。具体实现方法直接看代码就可以,非常容易理解。需要注意的一点是,下次应该从上次处理到的邻接边继续处理,而非从上次处理到的邻接边的下一条开始。

      最后说一下增广过程。增广过程非常简单,寻找增广路成功(当前节点处理到t

      )后,沿着你记录的路径走一遍,记录一路上的最小残量,然后从st

      更新流量即可。

      实现

      int source;         // 源点
      int sink;           // 汇点
      int p[max_nodes];   // 可增广路上的上一条弧的编号
      int num[max_nodes]; // 和 t 的最短距离等于 i 的节点数量
      int cur[max_nodes]; // 当前弧下标
      int d[max_nodes];   // 残量网络中节点 i 到汇点 t 的最短距离
      bool visited[max_nodes];
      
      // 预处理, 反向 BFS 构造 d 数组
      bool bfs()
      {
          memset(visited, 0, sizeof(visited));
          queue<int> Q;
          Q.push(sink);
          visited[sink] = 1;
          d[sink] = 0;
          while (!Q.empty()) {
              int u = Q.front();
              Q.pop();
              for (iterator_t ix = G[u].begin(); ix != G[u].end(); ++ix) {
                  Edge &e = edges[(*ix)^1];
                  if (!visited[e.from] && e.capacity > e.flow) {
                      visited[e.from] = true;
                      d[e.from] = d[u] + 1;
                      Q.push(e.from);
                  }
              }
          }
          return visited[source];
      }
      
      // 增广
      int augment()
      {
          int u = sink, df = __inf;
          // 从汇点到源点通过 p 追踪增广路径, df 为一路上最小的残量
          while (u != source) {
              Edge &e = edges[p[u]];
              df = min(df, e.capacity - e.flow);
              u = edges[p[u]].from;
          }
          u = sink;
          // 从汇点到源点更新流量
          while (u != source) {
              edges[p[u]].flow += df;
              edges[p[u]^1].flow -= df;
              u = edges[p[u]].from;
          }
          return df;
      }
      
      int max_flow()
      {
          int flow = 0;
          bfs();
          memset(num, 0, sizeof(num));
          for (int i = 0; i < num_nodes; i++) num[d[i]]++;
          int u = source;
          memset(cur, 0, sizeof(cur));
          while (d[source] < num_nodes) {
              if (u == sink) {
                  flow += augment();
                  u = source;
              }
              bool advanced = false;
              for (int i = cur[u]; i < G[u].size(); i++) { 
                  Edge& e = edges[G[u][i]];
                  if (e.capacity > e.flow && d[u] == d[e.to] + 1) {
                      advanced = true;
                      p[e.to] = G[u][i];
                      cur[u] = i;
                      u = e.to;
                      break;
                  }
              }
              if (!advanced) { // retreat
                  int m = num_nodes - 1;
                  for (iterator_t ix = G[u].begin(); ix != G[u].end(); ++ix)
                      if (edges[*ix].capacity > edges[*ix].flow)
                          m = min(m, d[edges[*ix].to]);
                  if (--num[d[u]] == 0) break; // gap 优化
                  num[d[u] = m+1]++;
                  cur[u] = 0;
                  if (u != source)
                      u = edges[p[u]].from;
              }
          }
          return flow;
      }


      我的Template:

      const int maxn=1000;
      const int INF=0x3f3f3f3f;
      int s,t,first[maxn],nxt[maxn];
      struct Edge{
      	int u,v,cap,flow;
      }e[maxn];
      bool vis[maxn];
      int q[maxn],d[maxn],p[maxn],num[maxn],cur[maxn];
      void bfs()//反向bfs进行预处理。构造d[]数组;
      {
          memset(vis,false,sizeof(vis));
          int head=0,tail=1;
          q[0]=t;
          d[t]=0;
          vis[t]=true;
          while(head!=tail){
              int now=q[head];head++;
              for(int i=first[now];i;i=nxt[i])
                  if(!vis[e[i].u]&&e[i].cap>e[i].flow){
                      vis[e[i].u]=true;
                      d[e[i].u]=d[now]+1;
                      q[tail++]=e[i].u;
                  }
          }
      }
      int Agument()
      {
          int x=t,a=INF;//从源点到汇点通过p数组追踪增广路径;
          while(x!=s){
              a=min(a,e[p[x]].cap-e[p[x]].flow);//找到最小容量;
              x=e[p[x]].u;
          }
          x=t;//更新流量;
          while(x!=s){
              e[p[x]].flow+=a;
              e[p[x]^1].flow-=a;
              x=e[p[x]].u;
          }
          return a;
      }
      int ISAP()
      {
          int flow=0;//flow记录最大流的int型变量;
          bfs();//宽搜构造d[]数组;
          memset(num,0,sizeof(num));//num即基数数组,记录在i这个距离上有几个点,进行gap优化;
          for(int i=1;i<=n;i++)num[d[i]]++;
          int x=s;
          memset(cur,0,sizeof(cur));//currently,即记录当前弧的数组;
          while(d[s]<n){
              if(x==t){
                  flow+=Agument();//若已到汇点则增广并更新流量;
                  x=s;//再从源点开始找新路径;
              }
              bool ok=false;
              for(int i=cur[x];i;i=nxt[i])
                  if(e[i].cap>e[i].flow&&d[x]==d[e[i].v]+1){
                      ok=true;
                      p[e[i].v]=i;
                      cur[x]=i;
                      x=e[i].v;
                      break;
                  }
              if(!ok){//retreat操作:更新d[]数组
                  int mn=n-1;
                  for(int i=first[x];i;i=nxt[i])
                      if(e[i].cap>e[i].flow)mn=min(mn,d[e[i].v]);
                  if(--num[d[x]]==0)break;//gap 优化:如果x是距离汇点最后一个distance为d[x]的点,那么意味着s、t之后不再连通,则退出搜索;
                  num[d[x]=mn+1]++;//不然再次记录;
                  cur[x]=0;
                  if(x!=s)x=e[p[x]].u;
              }
          }
          return flow;
      }
      void Link()//建立邻接表;
      {
          memset(nxt,0,sizeof(nxt));
          memset(first,0,sizeof(first));
          for(int a,b,c;m;m--){
              scanf("%d%d%d",&a,&b,&c);
              e[++ecnt].u=a,e[ecnt].v=b,e[ecnt].cap=c,e[ecnt].flow=0;
              nxt[ecnt]=first[a],first[a]=ecnt;
              e[++ecnt].u=b,e[ecnt].v=a,e[ecnt].cap=0,e[ecnt].flow=0;
              nxt[ecnt]=first[b],first[b]=ecnt;
          }
      }
      



  • 相关阅读:
    前端学习(21)~css学习:如何让一个元素水平垂直居中?
    前端学习(20)~css布局(十三)
    前端学习(19)~css3属性(十二):Flex布局图片详解
    前端学习(18)~css3属性学习(十一):动画详解
    前端学习(16)~css3属性学习(十)
    蓝桥---大臣的旅费(求树的直径)
    蓝桥---乘积最大(区间dp、数据水的话long long,暴力就能过)
    蓝桥---错误票据(注意输入)
    蓝桥---带分数(全排列)
    蓝桥---传球游戏(dp)
  • 原文地址:https://www.cnblogs.com/keshuqi/p/5957764.html
Copyright © 2011-2022 走看看