【资料】ISAP写法--最大流

来源：互联网发布：linux tail 行编辑：程序博客网时间：2024/05/04 19:48

转载大神的博客：点击打开链接

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

约定

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

引入

求解最大流问题的一个比较容易想到的方法就是，每次在残量网络（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 算法的框架内容就讲完了。对于代码本身，还有几个优化和实现的技巧需要说明。

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

最后说一下增广过程。增广过程非常简单，寻找增广路成功（当前节点处理到 t ）后，沿着你记录的路径走一遍，记录一路上的最小残量，然后从 s 到 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;}

自己的代码：

#include<cstdlib>#include<cstdio>#include<algorithm>#include<cstring>#include<iostream>#include<cmath>#include<vector>#include<queue>using namespace std;struct edge{int from,to,cap,flow;};vector<edge>edges;//edges[e]&edges[e^1]are couples;vector<int>G[50];int map[50][50],max_node,s,t;//s:source,t:terminal;int d[50],p[50],num[50],cur[50],mark_m[50],mark_d[50];//p[] is  to remark the former augmenting edge//cur:the current node's subscript;bool bfs(){    bool visit[max_node];    memset(visit,false,sizeof(visit));    queue<int>Q;    Q.push(t);visit[t]=1;    d[t]=0;    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(!visit[e.from]&&e.cap>e.flow)            {                visit[e.from]=true;                d[e.from]=d[u]+1;                Q.push(e.from);            }        }    }    return visit[s];}int augment(){    int x=t,df=0X3f3f3f3f;//df is the minimal flow of the road    while(x!=s)    {        edge &e= edges[p[x]];        df=min (df,e.cap-e.flow);        x=edges[p[x]].from;    }    x=t;    while(x!=s)    {        edge &e= edges[p[x]];        e.flow+=df;        edges[p[x]^1].flow-=df;        x=edges[p[x]].from;    }    return df;}int max_flow(){    int flow =0;    bfs();    memset(num,0,sizeof(num));    memset(cur,0,sizeof(cur));    for(int i=0;i<max_node;i++)    num[d[i]]++;    int x=s;    while(d[s]<max_node)    {        if(x==t)        {            flow+=augment();            x=s;        }        int ok=0;        for(int i=cur[x];i<G[x].size();i++)        {            edge &e =edges[G[x][i]];            if(e.cap>e.flow&&d[x]==d[e.to]+1)            {                ok=1;                p[e.to]=G[x][i];                cur[x]=i;                x = e.to ; break;            }        }        if(!ok)        {            int m=max_node-1;            for(int i=0;i<G[x].size();i++)            {               edge &e =edges[G[x][i]];               if(e.cap>e.flow) m = min(m,d[e.to]);            }            if(--num[d[x]]==0) break;            num[d[x]=m+1]++;            cur[x]=0;            if(x!=s) x=edges[p[x]].from;        }    }    return flow;}void addedges(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);}int main(){    int human,cat,dog,relation,i,j,k;    while(scanf("%d%d%d",&human,&cat,&dog)!=EOF)    {        memset(map,0,sizeof(map));        memset(mark_m,0,sizeof(mark_m));        memset(mark_d,0,sizeof(mark_d));        scanf("%d",&relation);        while(relation--)        {            int a,b;            scanf("%d%d",&a,&b);            if(a<=human) mark_m[a]++;            if(a>cat+human) mark_d[a]++;            if(b<=human) mark_m[b]++;            if(b>cat+human) mark_d[b]++;            map[a][b]=map[b][a]=1;        }        max_node=human+cat+dog+2;        s=0,t=max_node-1;        for(i=0;i<50;i++)         G[i].clear();        for( i=1;i<=human;i++)        addedges(0,i,cat-mark_m[i]);        for( i=human+cat+1;i<=human+cat+dog;i++)        addedges(i,human+cat+dog+1,cat-mark_d[i]);        for(i=1;i<=human;i++)        {            for(j=human+1;j<=human+cat;j++)            {                if(map[i][j]==1)                continue;                for(k=human+cat+1;k<=human+cat+dog;k++)                {                    if(map[j][k]==1) continue;                    addedges(i,j,1);                    addedges(j,k,1);                }            }        }        printf("%d\n",max_flow());    }    return 0;}

0 0