引入

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

人们发现,如果每次都沿着残量网络中的最短增广路增广,则运行时间可以减为 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∣‾‾‾√)。虽然在复杂度上比增广路方法进步很多,但是预流推进算法复杂度的上界是比较紧的,因此有时差距并不会很大。

算法解释

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

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

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

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

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

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


书上的模版:

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <vector>
#include <queue>
using namespace std;
#define N 1010
#define INF 0x3f3f3f3f

struct Edge {
    int from, to, cap, flow;
    Edge() {}
    Edge(int from, int to, int cap, int flow): from(from), to(to), cap(cap), flow(flow) {}
};

struct ISAP {
    int p[N], num[N], cur[N], d[N];
    int t, s, n, m;
    bool vis[N];

    vector<int> G[N];
    vector<Edge> edges;

    void init(int n, int m) {
        this->n = n; this->m = m;
        for (int i = 0; i <= n; i++) {
            G[i].clear();
            d[i] = INF;
        }
        edges.clear();
    }

    void AddEdge(int from, int to, int cap) {
        edges.push_back(Edge(from, to, cap, 0));
        edges.push_back(Edge(to, from, 0, 0));
        int m = edges.size();
        G[from].push_back(m - 2);
        G[to].push_back(m - 1);
    }

    bool BFS() {
        memset(vis, 0, sizeof(vis));

        queue<int> Q;
        d[t] = 0;
        vis[t] = 1;
        Q.push(t);

        while (!Q.empty()) {
            int u = Q.front();
            Q.pop();

            for (int i = 0; i < G[u].size(); i++) {
                Edge &e = edges[G[u][i] ^ 1];
                if (!vis[e.from] && e.cap > e.flow) {
                    vis[e.from] = true;
                    d[e.from] = d[u] + 1;
                    Q.push(e.from);
                }
            }
        }
        return vis[s];
    }

    int Augment() {
        int u = t, flow = INF;
        while (u != s) {
            Edge &e = edges[p[u]];
            flow = min(flow, e.cap - e.flow);
            u = edges[p[u]].from;
        }

        u = t;
        while (u != s) {
            edges[p[u]].flow += flow;
            edges[p[u] ^ 1].flow -= flow;
            u = edges[p[u]].from;
        }
        return flow;
    }

    int Maxflow(int s, int t) {
        this->s = s; this->t = t;
        int flow = 0;
        BFS();
        if (d[s] >= n)
            return 0;

        memset(num, 0, sizeof(num));
        memset(cur, 0, sizeof(cur));
        for (int i = 0; i < n; i++)
            if (d[i] < INF)
                num[d[i]]++;
        int u = s;

        while (d[s] < n) {
            if (u == t) {
                flow += Augment();
                u = s;
            }
            bool ok = false;
            for (int i = cur[u]; i < G[u].size(); i++) {
                Edge &e = edges[G[u][i]];
                if (e.cap > e.flow && d[u] == d[e.to] + 1) {
                    ok = true;
                    p[e.to] = G[u][i]; 
                    cur[u] = i;
                    u = e.to;
                    break;
                }
            }

            if (!ok) {
                int Min = n - 1;
                for (int i = 0; i < G[u].size(); i++) {
                    Edge &e = edges[G[u][i]];
                    if (e.cap > e.flow)
                        Min = min(Min, d[e.to]);
                }
                if (--num[d[u]] == 0)
                    break;
                num[d[u] = Min + 1]++;
                cur[u] = 0;
                if (u != s)
                    u = edges[p[u]].from;
            }
        }
        return flow;
    }
};

ISAP isap;

int main() {
    return 0;
}

第二种,改进的

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <vector>
#include <queue>
using namespace std;
const int MAXNODE = 210;
const int MAXEDGE = 100010;
typedef int Type;
const Type INF = 0x3f3f3f3f;

struct Edge {
    int u, v, next;
    Type cap, flow;
    Edge() {}
    Edge(int u, int v, Type cap, Type flow, int next): u(u), v(v), cap(cap), flow(flow), next(next) {}
};

struct ISAP {
    int n, m, s, t;
    Edge edges[MAXEDGE];
    int head[MAXNODE], p[MAXNODE], num[MAXNODE], cur[MAXNODE], d[MAXNODE];
    bool vis[MAXNODE];

    void init(int n) {
        this->n = n;
        memset(head, -1, sizeof(head));
        m = 0;
    }

    void AddEdge(int u, int v, Type cap) {
        edges[m] = Edge(u, v, cap, 0, head[u]);
        head[u] = m++;
        edges[m] = Edge(v, u, 0, 0, head[v]);
        head[v] = m++;
    }

    bool BFS() {
        memset(vis, 0, sizeof(vis));
        queue<int> Q;
        for (int i = 0; i < n; i++)
            d[i] = INF;
        d[t] = 0;
        vis[t] = 1;
        Q.push(t);

        while (!Q.empty()) {
            int u = Q.front(); Q.pop();

            for (int i = head[u]; ~i; i = edges[i].next) {
                Edge &e = edges[i ^ 1];
                if (!vis[e.u] && e.cap > e.flow) {
                    vis[e.u] = true;
                    d[e.u] = d[u] + 1;
                    Q.push(e.u);
                }
            }
        }
        return vis[s];
    }

    Type Augment() {
        int u = t;
        Type flow = INF;
        while (u != s) {
            Edge &e = edges[p[u]];
            flow = min(flow, e.cap - e.flow);
            u = edges[p[u]].u;
        }

        u = t;
        while (u != s) {
            edges[p[u]].flow += flow;
            edges[p[u] ^ 1].flow -= flow;
            u = edges[p[u]].u;
        }
        return flow;
    }

    Type Maxflow(int s, int t) {
        this->s = s; this->t = t;
        Type flow = 0;
        BFS();
        //如果s-->t走不通
        if (d[s] >= n)
            return 0;

        memset(num, 0, sizeof(num));
        for (int i = 0; i < n; i++)
            cur[i] = head[i];

        for (int i = 0; i < n; i++)
            if (d[i] < INF) 
                num[d[i]]++;

        int u = s;
        while (d[s] < n) {
            if (u == t) {
                flow += Augment();
                u = s;
            }

            bool ok = false;//纪录是否找到了下一个点
            for (int i = cur[u]; ~i; i = edges[i].next) {
                Edge &e = edges[i];
                if (e.cap > e.flow && d[u] == d[e.v] + 1) {
                    ok = true;
                    p[e.v] = i;//点v由第i条边增广得到
                    cur[u] = i;//尝试到第i条边
                    u = e.v;
                    break;
                }
            }

            //如果没找到下一个点,表示u到t的最短路要变长了,或者没路可走了
            if (!ok) {
                //找寻u到下一个点的最短路
                int Min = n - 1;
                for (int i = head[u]; ~i; i = edges[i].next) {
                    Edge &e = edges[i];
                    if (e.cap > e.flow)
                        Min = min(Min, d[e.v]);
                }
                if (--num[d[u]] == 0)//GAP优化
                    break;
                num[d[u] = Min + 1]++;
                cur[u] = head[u];
                //返回前一个点,因为该点的最短距离已经变了
                if (u != s)
                    u = edges[p[u]].u;
            }
        }
        return flow;
    }
}isap;

int main() {
    return 0;
}