网络流－最大流问题 ISAP 算法解释（转自Renfei Song's Blog）
网络流－最大流问题 ISAP 算法解释
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; }
我的Template：const int maxn=1000; const int INF=0x3f3f3f3f; int s,t,first[maxn],nxt[maxn]; struct Edge{ int u,v,cap,flow; }e[maxn]; bool vis[maxn]; int q[maxn],d[maxn],p[maxn],num[maxn],cur[maxn]; void bfs()//反向bfs进行预处理。构造d[]数组; { memset(vis,false,sizeof(vis)); int head=0,tail=1; q[0]=t; d[t]=0; vis[t]=true; while(head!=tail){ int now=q[head];head++; for(int i=first[now];i;i=nxt[i]) if(!vis[e[i].u]&&e[i].cap>e[i].flow){ vis[e[i].u]=true; d[e[i].u]=d[now]+1; q[tail++]=e[i].u; } } } int Agument() { int x=t,a=INF;//从源点到汇点通过p数组追踪增广路径; while(x!=s){ a=min(a,e[p[x]].cap-e[p[x]].flow);//找到最小容量; x=e[p[x]].u; } x=t;//更新流量; while(x!=s){ e[p[x]].flow+=a; e[p[x]^1].flow-=a; x=e[p[x]].u; } return a; } int ISAP() { int flow=0;//flow记录最大流的int型变量; bfs();//宽搜构造d[]数组; memset(num,0,sizeof(num));//num即基数数组，记录在i这个距离上有几个点，进行gap优化; for(int i=1;i<=n;i++)num[d[i]]++; int x=s; memset(cur,0,sizeof(cur));//currently，即记录当前弧的数组; while(d[s]<n){ if(x==t){ flow+=Agument();//若已到汇点则增广并更新流量; x=s;//再从源点开始找新路径; } bool ok=false; for(int i=cur[x];i;i=nxt[i]) if(e[i].cap>e[i].flow&&d[x]==d[e[i].v]+1){ ok=true; p[e[i].v]=i; cur[x]=i; x=e[i].v; break; } if(!ok){//retreat操作：更新d[]数组 int mn=n-1; for(int i=first[x];i;i=nxt[i]) if(e[i].cap>e[i].flow)mn=min(mn,d[e[i].v]); if(--num[d[x]]==0)break;//gap 优化：如果x是距离汇点最后一个distance为d[x]的点，那么意味着s、t之后不再连通，则退出搜索; num[d[x]=mn+1]++;//不然再次记录; cur[x]=0; if(x!=s)x=e[p[x]].u; } } return flow; } void Link()//建立邻接表; { memset(nxt,0,sizeof(nxt)); memset(first,0,sizeof(first)); for(int a,b,c;m;m--){ scanf("%d%d%d",&a,&b,&c); e[++ecnt].u=a,e[ecnt].v=b,e[ecnt].cap=c,e[ecnt].flow=0; nxt[ecnt]=first[a],first[a]=ecnt; e[++ecnt].u=b,e[ecnt].v=a,e[ecnt].cap=0,e[ecnt].flow=0; nxt[ecnt]=first[b],first[b]=ecnt; } }
