天天看点

最大流Dinic算法讲解 && ISAP 算法解释短路增值算法(MPLA)Dinic算法 网络流-最大流问题 ISAP 算法解释

转自:http://blog.csdn.net/wall_f/article/details/8207595

为了更好的介绍Dinic算法,我们先来介绍最短增广路算法。

最短增广路算法

1、顶点的层次和层次网络

顶点的层次:在残留网络中,把从源点到顶点u的最短路径长度(该长度仅仅是值路径上边的数目,与容量无关),称为顶点u的层次,记为level(u)。源点Vs的层次为0。

将残留网络中所有的顶点的层次标注出来的过程称为分层。

注意:

(1)对残留网路进行分层后,弧可能有3种可能的情况。

1、从第i层顶点指向第i+1层顶点。

2、从第i层顶点指向第i层顶点。

3、从第i层顶点指向第j层顶点(j < i)。

(2)不存在从第i层顶点指向第i+k层顶点的弧(k>=2)。

(3)并非所有的网络都能分层。

层次网络:对残留网络进行分层后,删去比汇点Vt层次更高的顶点和与汇点Vt同层的顶点(保留Vt),并删去这些顶点相关联的弧,再删去从某层顶点指向同层顶点和低层顶点的弧,所剩余的各条弧的容量与残留网络中的容量相同,这样得到的网络就是残留网络的子网络,称为层次网络,记为G''(V'',E'')。

根据层次网络定义,层次网络中任意的一条弧<u,v>,有满足level(u)+1 == level(v),这条弧也叫允许弧。直观的说,层次网络是建立在残留网络基础之上的一张“最短路径图”。从源点开始,在层次网络中沿着边不管怎么走,到达一个终点之后,经过的路径一定是终点在残留网络中的最短路径。

阻塞流:设容量网络的一个可行流为f,当该网络的层次G''中不存在增广路(即从源点Vs到汇点Vt的路径时),称该可行流f为层次网络G''的阻塞流。

2、最短路增广路径的算法思想

最短增广路的算法思想是:每次在层次网络中找一条含弧数最少的增广路进行增广。最短增广路算法的具体步骤如下:

(1)初始化容量网络和网络流。

(2)构造残留网络和层次网络,若汇点不在层次网络中,则算法结束。

(3)在层次网络中不断用BFS增广,直到层次网络中没有增广路为止;每次增广完毕,在层次网络中要去掉因改进流量而导致饱和的弧。

(4)转步骤(2)。

在最短增广路算法中,第(2)、(3)步被循环执行,将执行(2)、(3)步的一次循环称为一个阶段。在每个阶段中,首先根据残留网络建立层次网络,然后不断用BFS在层次网络中增广,直到出现阻塞流。注意每次增广后,在层次网络中要去掉因改进流量而饱和的弧。该阶段的增广完毕后,进入下一阶段。这样的不断重复,直到汇点不在层次网络中为止。汇点不在层次网络中意味着在残留网络中不存在一条从源点到汇点的路径,即没有增广路。

在程序实现的时候,并不需要真正“构造”层次网络,只需要对每个顶点标记层次,增广的时候,判断边是否满足level(v) = level(u)+1这一约束条件即可。

3、最短增广路算法复杂度分析

最短增广路的复杂度包括建立层次网络和寻找增广路两部分。

在最短增广路中,最多建立n个层次网络,每个层次网络用BFS一次遍历即可得到。一次BFS的复杂度为O(m),所以建层次图的总复杂度为O(n*m)。

每增广一次,层次网络中必定有一条边会被删除。层次网络中最多有m条边,所以认为最多可以增广m次。在最短增广路算法中,用BFS来增广,一次增广的复杂度为O(n+m),其中O(m)为BFS的花费,O(n)为修改流量的花费。所以在每一阶段寻找增广路的复杂度为O(m*(m+n)) = O(m*m)。因此n个阶段寻找增广路的算法总复杂度为O(n*m*m)。

两者之和为O(n*m*m)。

以上介绍最短增广路算法只是为下面介绍Dinic算法而提供给大家一个铺垫,有了以上的基础,接下来我们来介绍Dinic算法,Dinic其实是最短增广路算法的优化版。

连续最短增广路算法----Dinic算法

1、Dinic算法思路

Dinic算法的思想也是分阶段地在层次网络中增广。它与最短增广路算法不同之处是:最短增广路每个阶段执行完一次BFS增广后,要重新启动BFS从源点Vs开始寻找另一条增广路;而在Dinic算法中,只需一次DFS过程就可以实现多次增广,这是Dinic算法的巧妙之处。Dinic算法具体步骤如下:

(1)初始化容量网络和网络流。

(2)构造残留网络和层次网络,若汇点不再层次网络中,则算法结束。

(3)在层次网络中用一次DFS过程进行增广,DFS执行完毕,该阶段的增广也执行完毕。

(4)转步骤(2)。

在Dinic的算法步骤中,只有第(3)步与最短增广路相同。在下面实例中,将会发现DFS过程将会使算法的效率有非常大的提高。

Dinic算法复杂度分析

与最短增广路算法一样,Dinic算法最多被分为n个阶段,每个阶段包括建层次网络和寻找增广路两部分,其中建立层次网络的复杂度仍是O(n*m)。

现在来分析DFS过程的总复杂度。在每一阶段,将DFS分成两部分分析。

(1)修改增广路的流量并后退的花费。在每一阶段,最多增广m次,每次修改流量的费用为O(n)。而一次增广后在增广路中后退的费用也为O(n)。所以在每一阶段中,修改增广路以及后退的复杂度为O(m*n)。

(2)DFS遍历时的前进与后退。在DFS遍历时,如果当前路径的最后一个顶点能够继续扩展,则一定是沿着第i层的顶点指向第i+1层顶点的边向汇点前进了一步。因为增广路经长度最长为n,所以最坏的情况下前进n步就会遇到汇点。在前进的过程中,可能会遇到没有边能够沿着继续前进的情况,这时将路径中的最后一个点在层次图中删除。

注意到每后退一次必定会删除一个顶点,所以后退的次数最多为n次。在每一阶段中,后退的复杂度为O(n)。

假设在最坏的情况下,所有的点最后均被退了回来,一共共后退了n次,这也就意味着,有n次的前进被“无情”地退了回来,这n次前进操作都没有起到“寻找增广路”的作用。除去这n次前进和n次后退,其余的前进都对最后找到增广路做了贡献。增广路最多找到m次。每次最多前进n个点。所以所有前进操作最多为n+m*n次,复杂度为O(n*m)。

于是得到,在每一阶段中,DFS遍历时前进与后退的花费为O(m*n)。

综合以上两点,一次DFS的复杂度为O(n*m)。因此,Dinic算法的总复杂度即O(m*n*n)。

下面给出 Dinic具体实现的代码。

预处理:

[cpp]  view plain  copy

  1. #include <iostream>  
  2. #include <cstdlib>  
  3. #include <cstdio>  
  4. #include <cstring>  
  5. #include <string>  
  6. #include <algorithm>  
  7. #include <queue>  
  8. using namespace std;  
  9. const int MAXN = 210;  
  10. const int MAXM = 210*210;  
  11. const int INF = 0x3f3f3f3f;  
  12. struct Edge  
  13. {  
  14.     int v, f;  
  15.     int next;  
  16. }edge[MAXM];  
  17. int n, m;  
  18. int cnt;  
  19. int first[MAXN], level[MAXN];  
  20. int q[MAXN];  
  21. void init()  
  22. {  
  23.     cnt = 0;  
  24.     memset(first, -1, sizeof(first));  
  25. }  
  26. void read_graph(int u, int v, int f)  
  27. {  
  28.     edge[cnt].v = v, edge[cnt].f = f;  
  29.     edge[cnt].next = first[u], first[u] = cnt++;  
  30.     edge[cnt].v = u, edge[cnt].f = 0;  //增加一条反向弧,容量为0   
  31.     edge[cnt].next = first[v], first[v] = cnt++;  
  32. }  

BFS构建层次网络:

[cpp]  view plain  copy

  1. int bfs(int s, int t) //构建层次网络   
  2. {  
  3.     memset(level, 0, sizeof(level));  
  4.     level[s] = 1;  
  5.     int front = 0, rear = 1;  
  6.     q[front] = s;  
  7.     while(front < rear)  
  8.     {  
  9.         int x = q[front++];  
  10.         if(x == t) return 1;  
  11.         for(int e = first[x]; e != -1; e = edge[e].next)  
  12.         {  
  13.             int v = edge[e].v, f = edge[e].f;  
  14.             if(!level[v] && f)  
  15.             {  
  16.                 level[v] = level[x] + 1;  
  17.                 q[rear++] = v;  
  18.             }  
  19.         }  
  20.     }  
  21.     return 0;  
  22. }  

DFS找寻增广路:

[cpp]  view plain  copy

  1. int dfs(int u, int maxf, int t)  
  2. {  
  3.     if(u == t) return maxf;  
  4.     int ret = 0;  
  5.     for(int e = first[u]; e != -1; e = edge[e].next)  
  6.     {  
  7.         int v = edge[e].v, f = edge[e].f;  
  8.         if(level[u] + 1 == level[v] && f)  
  9.         {  
  10.             int Min = min(maxf-ret, f);  
  11.             f = dfs(v, Min, t);  
  12.             edge[e].f -= f;  
  13.             edge[e^1].f += f;  
  14.             ret += f;  
  15.             if(ret == maxf) return ret;  
  16.         }  
  17.     }  
  18.     return ret;  
  19. }  

Dinic:

[cpp]  view plain  copy

  1. int Dinic(int s, int t) //Dinic  
  2. {  
  3.     int ans = 0;  
  4.     while(bfs(s, t)) ans += dfs(s, INF, t);  
  5.     return ans;  
  6. }  

层次图

最大流Dinic算法讲解 &amp;&amp; ISAP 算法解释短路增值算法(MPLA)Dinic算法 网络流-最大流问题 ISAP 算法解释
最大流Dinic算法讲解 &amp;&amp; ISAP 算法解释短路增值算法(MPLA)Dinic算法 网络流-最大流问题 ISAP 算法解释

短路增值算法(MPLA)

算法流程

最大流Dinic算法讲解 &amp;&amp; ISAP 算法解释短路增值算法(MPLA)Dinic算法 网络流-最大流问题 ISAP 算法解释

汇点不在层次图内意味着在剩余图中不存在一条从源点到汇点的路径,即没有增广路。

在程序实现的时候,层次图并不用被“建”出来,我们只需对每个顶点标记层次,增广的时候,判断边是否满足level(u) +1= level(v)这一约束即可。

Dinic算法

Dinic算法的思想也是分阶段地在层次图中增广。

它与最短路径增值算法不同之处是:在Dinic算法中,我们用一个dfs过程代替多次bfs来寻找阻塞流。下面给出其算法步骤:

算法流程

最大流Dinic算法讲解 &amp;&amp; ISAP 算法解释短路增值算法(MPLA)Dinic算法 网络流-最大流问题 ISAP 算法解释

增广过程图解

最大流Dinic算法讲解 &amp;&amp; ISAP 算法解释短路增值算法(MPLA)Dinic算法 网络流-最大流问题 ISAP 算法解释
最大流Dinic算法讲解 &amp;&amp; ISAP 算法解释短路增值算法(MPLA)Dinic算法 网络流-最大流问题 ISAP 算法解释

伪代码描述

最大流Dinic算法讲解 &amp;&amp; ISAP 算法解释短路增值算法(MPLA)Dinic算法 网络流-最大流问题 ISAP 算法解释

在程序里,p表示找到的增广路径,p.top为路径中的最后一个顶点。一开始,p中只有源点。

整个While循环分为2个操作。如果p的最后一个顶点为汇点,也就是说找到了增广路,那么对p增广,注意到增广后一定有一条或多条p中的边被删除了。这时,我们使增广路径后退至p中从源点可到达的最后一个顶点。

如果p的最后一个顶点不为汇点,那么观察最后那个的顶点u 。若在层次图中存在从u连出的一条边,比如(u,v),我们就将顶点v放入路径p中,继续dfs遍历;否则,点u对之后的dfs遍历就没有用了,我们将点u以及层次图中连到u的所有边删除,并且在p中后退一个点。

Dfs过程将会不断重复这2个操作,直到从源点连出的边全部被删除为止。

转自:http://www.renfei.org/blog/isap.html

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

August 7, 2013 /  编程指南

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

约定

我们使用邻接表来表示图,表示方法可以见文章带权最短路 Dijkstra, SPFA, Bellman-Ford, ASP, Floyd-Warshall 算法分析或二分图的最大匹配、完美匹配和匈牙利算法的开头(就不重复贴代码了)。在下文中,图的源点(source)表示为  s s ,汇点(sink)表示为  t t ,当前节点为  u u 。建图时,需要建立双向边(设反向的边容量为  0 0)才能保证算法正确。

引入

求解最大流问题的一个比较容易想到的方法就是,每次在残量网络(residual network)中任意寻找一条从  s s 到  t t 的路径,然后增广,直到不存在这样的路径为止。这就是一般增广路算法(labeling algorithm)。可以证明这种不加改进的贪婪算法是正确的。假设最大流是  f f,那么它的运行时间为 O( f⋅∣E∣) O( f⋅∣E∣) 。但是,这个运行时间并不好,因为它和最大流  f f 有关。

人们发现,如果每次都沿着残量网络中的最短增广路增广,则运行时间可以减为  O(∣E∣2⋅∣V∣)  O(∣E∣2⋅∣V∣)  。这就是最短增广路算法。而 ISAP 算法则是最短增广路算法的一个改进。其实,ISAP 的意思正是「改进的最短增广路」 (Improved Shortest Augmenting Path)。

顺便说一句,上面讨论的所有算法根本上都属于增广路方法(Ford-Fulkerson method)。和它对应的就是大名鼎鼎的预流推进方法(Preflow-push method)。其中最高标号预流推进算法(Highest-label preflow-push algorithm)的复杂度可以达到  O(∣V∣2∣E∣−−−√) O(∣V∣2∣E∣)。虽然在复杂度上比增广路方法进步很多,但是预流推进算法复杂度的上界是比较紧的,因此有时差距并不会很大。

算法解释

概括地说,ISAP 算法就是不停地找最短增广路,找到之后增广;如果遇到死路就 retreat,直到发现 s s,  t t不连通,算法结束。找最短路本质上就是无权最短路径问题,因此采用 BFS 的思想。具体来说,使用一个数组 d d,记录每个节点到汇点 t t的最短距离。搜索的时候,只沿着满足 d[u]=d[v]+1 d[u]=d[v]+1的边 u→v u→v(这样的边称为允许弧)走。显然,这样走出来的一定是最短路。

原图存在两种子图,一个是残量网络,一个是允许弧组成的图。残量网络保证可增广,允许弧保证最短路(时间界较优)。所以,在寻找增广路的过程中,一直是在残量网络中沿着允许弧寻找。因此,允许弧应该是属于残量网络的,而非原图的。换句话说,我们沿着允许弧,走的是残量网络(而非原图)中的最短路径。当我们找到沿着残量网络找到一条增广路,增广后,残量网络肯定会变化(至少少了一条边),因此决定允许弧的 d d数组要进行相应的更新(顺便提一句,Dinic 的做法就是每次增广都重新计算 d d数组)。然而,ISAP 「改进」的地方之一就是,其实没有必要马上更新 d d数组。这是因为,去掉一条边只可能令路径变得更长,而如果增广之前的残量网络存在另一条最短路,并且在增广后的残量网络中仍存在,那么这条路径毫无疑问是最短的。所以,ISAP 的做法是继续增广,直到遇到死路,才执行 retreat 操作。

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

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

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

最后说一下增广过程。增广过程非常简单,寻找增广路成功(当前节点处理到 t t)后,沿着你记录的路径走一遍,记录一路上的最小残量,然后从 s s到 t t更新流量即可。

实现

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;
}