转：网络流---isap

引入

求解最大流问题的一个比较容易想到的方法就是，每次在残量网络（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更新流量即可。

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

模板：

#include <queue> 
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
const int MAXNODE = 210;
const int MAXEDGE = 10010;
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[MAXNODE];
    int head[MAXEDGE], 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!=-1; 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!=-1; 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;   // 此n-1表示距离； 
                for(int i=head[u]; i!=-1; i=edges[i].next){
                    Edge &e=edges[i];
                if(e.cap>e.flow)
                    Min=min(Min, d[e.v]);
                }
                if(--num[d[u]]==0)
                    break;
                num[d[u]=Min+1]++;
                cur[u]=head[u];
                //返回前一个点 ; 
                if(u!=s)
                    u=edges[p[u]].u;
            }
        }
        return flow; 
    }
}isap;
int main()
{
    //isap.init(100);
    //isap.Maxflow(0, 99);
    return 0;
}