ISAP學習筆記

ISAP是OI中求最大流的常用方法之一。相對於Dinic，ISAP的速度提升了很多，但編碼複雜度也上升了不少。

約定

採用鄰接表儲存圖，對於每條弧，增加一條容量為0的逆向邊。
d陣列代表每個點到終點的距離。

引入

求解最大流有一種基礎方法，每次在殘量網路中dfs找到任意路徑增廣，直到不存在這樣的道路為止。這就是最短路增廣演算法。但是這樣的方法效率太低。
所以我們比較容易想到的一種方法是每次沿最短路增廣。這就是最短路增廣演算法。而ISAP則是針對最短路增廣的一種改進(Improved Shortest Augmenting Path)。

演算法概述

ISAP基於這樣一個事實，每次增廣之後，殘量網路中的最短路不會變短，但如果殘量網路中存在另一條最短路，那麼可以繼續增廣，直到遇到死路為止，所以每次增廣後再次bfs找最短路不是必須的(而Dinic則需要每次增廣後bfs)。

那麼遇到死路後怎麼辦呢?對於節點u，如果它周圍的任意節點v都不滿足d[v]=d[u]+1,那麼記min為d[v]的最小值,d[u]=d[v]+1。並退回上次經過的節點繼續尋找增廣路。

但是對於一種特殊情況，如果d[u]=k的節點已完全消失，那麼可以直接終止程式，因為此時d值>k的點以無法經過d值等於k的點到達t(d[t]=0),也就是說s必定無法到達k，這就是大名鼎鼎的gap優化。

說到這裡，演算法流程已大致清楚了，但還有一些程式實現的細節需要注意。

實現

#include <bits/stdc++.h>
using namespace std;

const int maxn = 10005, maxm = 100005;
const int Inf = 1 << 30;

struct edge {
    int from, to, next;
    int cap, flow;
}e[maxm * 2];
int h[maxn], tot;
int n, m, s, t, flow;

int cur[maxn]; //當前弧編號
int p[maxn]; //增廣路上一條邊編號(第i點的入邊)
 

int num[maxn]; //和t最短路距離為i的點數量
int d[maxn]; //殘量網路中到t的最短距離

inline void add(int u, int v, int w)
{
    e[tot++] = (edge) {u, v, h[u], w, 0};
    h[u] = tot - 1;
}

inline int gi()
{
    char c = getchar();
    while(c < '0' || c > '9') c = getchar();
    int sum = 0;
    while('0' <= c && c <= '9') sum = sum * 10 + c - 48, c = getchar();
    return sum;
}

queue<int> que;
//預處理,反向bfs構造d陣列
bool bfs()
{
    memset(d, 63, sizeof(int) * (n + 1));
    d[t] = 0;
    que.push(t);
    while(!que.empty()) {
        int u = que.front(); que.pop();
        ++num[d[u]];
        for(int v, i = h[u]; ~i; i = e[i].next) {
            v = e[i].to;
            if(e[i].cap == 0 && d[v] > d[u] + 1) {
                d[v] = d[u] + 1;
                que.push(v);
            }
        }
    }
    return d[s] < n;
}

//增廣
int augment()
{
    int u = t, Min = Inf;
    //從匯點到源點通過p追蹤增廣路徑,Min為一路上的最小流量
    while(u != s) {
        edge E = e[p[u]];
        Min = min(Min, E.cap - E.flow);
        u = e[p[u]].from;
    }
    u = t;
    while(u != s) {
        e[p[u]].flow += Min;
        e[p[u] ^ 1].flow -= Min;
        u = e[p[u]].from;
    }
    return Min;
}

int main()
{
    n = gi(); m = gi(); s = gi(); t = gi();

    memset(h, -1, sizeof(int) * (n + 1));
    for(int u, v, w, i = 1; i <= m; ++i) {
        u = gi(); v = gi(); w = gi();
        add(u, v, w); add(v, u, 0);
    }

    if(!bfs()) printf("0\n"), exit(0);

    int u = s;
    memcpy(cur, h, sizeof(int) * (n + 1));
    while(d[s] < n) {
        if(u == t) {
            flow += augment();
            u = s;
        }
        bool advanced = false;
        for(int i = cur[u]; ~i; i = e[i].next) {
            int v = e[i].to;
            if(e[i].cap > e[i].flow && d[u] == d[v] + 1) {
                advanced = true;
                p[v] = i;  cur[u] = i; u = v;
                break;
            }
        }
        if(!advanced) { //retreat操作
            int Min = n - 1;
            for(int i = h[u]; ~i; i = e[i].next)
                if(e[i].cap > e[i].flow) Min = min(Min, d[e[i].to]);
            if(--num[d[u]] == 0) break; //gap優化
            num[d[u] = Min + 1]++;
            cur[u] = h[u];
            if(u != s) u = e[p[u]].from;
        }
    }

    printf("%d\n", flow);

    return 0;
}