最大流 ISAP

Posted jpphy0

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了最大流 ISAP相关的知识,希望对你有一定的参考价值。

引入
求解最大流问题的一个比较容易想到的方法就是,每次在残量网络(residual network)中任意寻找一条从 s s s t t t的路径,然后增广,直到不存在这样的路径为止。这就是一般增广路算法(labeling algorithm)。可以证明这种不加改进的贪婪算法是正确的。假设最大流是 f f f,那么它的运行时间为 O (   f ⋅ ∣ E ∣ ) O\\left(\\ f\\cdot \\mid E\\mid\\right) O( fE)。但是,这个运行时间并不好,因为它和最大流 f f f有关。
人们发现,如果每次都沿着残量网络中的最短增广路增广,则运行时间可以减为$O\\left(\\mid E\\mid^2\\cdot\\mid V\\mid\\right)$。这就是最短增广路算法。而 ISAP 算法则是最短增广路算法的一个改进。其实,ISAP 的意思正是「改进的最短增广路」 (Improved Shortest Augmenting Path)。
顺便说一句,上面讨论的所有算法根本上都属于增广路方法(Ford-Fulkerson method)。和它对应的就是大名鼎鼎的预流推进方法(Preflow-push method)。其中最高标号预流推进算法(Highest-label preflow-push algorithm)的复杂度可以达到 O ( ∣ V ∣ 2 ∣ E ∣ ) O\\left(\\mid V\\mid ^2 \\sqrt{\\mid E\\mid}\\right) O(V2E )。虽然在复杂度上比增广路方法进步很多,但是预流推进算法复杂度的上界是比较紧的,因此有时差距并不会很大。
算法解释
概括地说,ISAP 算法就是不停地找最短增广路,找到之后增广;如果遇到死路就 retreat,直到发现 s s s, t t t不连通,算法结束。找最短路本质上就是无权最短路径问题,因此采用 BFS 的思想。具体来说,使用一个数组 d d d,记录每个节点到汇点 t t t的最短距离。搜索的时候,只沿着满足 d [ u ] = d [ v ] + 1 d[u]=d[v]+1 d[u]=d[v]+1的边 u → v u\\rightarrow v uv(这样的边称为允许弧)走。显然,这样走出来的一定是最短路。
原图存在两种子图,一个是残量网络,一个是允许弧组成的图。残量网络保证可增广,允许弧保证最短路(时间界较优)。所以,在寻找增广路的过程中,一直是在残量网络中沿着允许弧寻找。因此,允许弧应该是属于残量网络的,而非原图的。换句话说,我们沿着允许弧,走的是残量网络(而非原图)中的最短路径。当我们找到沿着残量网络找到一条增广路,增广后,残量网络肯定会变化(至少少了一条边),因此决定允许弧的 d d d数组要进行相应的更新(顺便提一句,Dinic 的做法就是每次增广都重新计算 d d d数组)。然而,ISAP 「改进」的地方之一就是,其实没有必要马上更新 d d d数组。这是因为,去掉一条边只可能令路径变得更长,而如果增广之前的残量网络存在另一条最短路,并且在增广后的残量网络中仍存在,那么这条路径毫无疑问是最短的。所以,ISAP 的做法是继续增广,直到遇到死路,才执行 retreat 操作。
说到这里,大家应该都猜到了,retreat 操作的主要任务就是更新 d d d数组。那么怎么更新呢?非常简单:假设是从节点 u u u找遍了邻接边也没找到允许弧的;再设一变量 m m m,令 m m m等于残量网络中 u u u的所有邻接点的 d d d数组的最小值,然后令 d [ u ] d[u] d[u]等于 m + 1 m+1 m+1即可。这是因为,进入 retreat 环节说明残量网络中 u u u t t t已经不能通过(已过时)的允许弧相连,那么 u u u t t t实际上在残量网络中的最短路的长是多少呢?(这正是 d d d的定义!)显然是残量网络中 u u u的所有邻接点和 t t t的距离加 1 1 1的最小情况。特殊情况是,残量网络中 u u u根本没有邻接点。如果是这样,只需要把 d [ u ] d[u] d[u]设为一个比较大的数即可,这会导致任何点到 u u u的边被排除到残量网络以外。(严格来说只要大于等于 ∣ V ∣ \\mid V\\mid V即可。由于最短路一定是无环的,因此任意路径长最大是 ∣ V ∣ − 1 \\mid V\\mid−1 V1)。修改之后,只需要把正在研究的节点 u u u沿着刚才走的路退一步,然后继续搜索即可。
讲到这里,ISAP 算法的框架内容就讲完了。对于代码本身,还有几个优化和实现的技巧需要说明。

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

以上是关于最大流 ISAP的主要内容,如果未能解决你的问题,请参考以下文章

最大流 ISAP

HDU-4280-Island Transport(网络流,最大流, ISAP)

模板ISAP最大流

算法学习笔记(8.1): 网络最大流算法 EK, Dinic, ISAP

圆桌问题 网络流Isap算法

CCF(引水入城:60分):最大流+ISAP算法