ISPA

来自CSDN的Rachel Zhang

4. Improved SAP 算法

本次介绍的重头戏。通常的 SAP 类算法在寻找增广路时总要先进行 BFS，BFS 的最坏情况下复杂度为 O(E)，这样使得普通 SAP 类算法最坏情况下时间复杂度达到了 O(VE2)。为了避免这种情况，Ahuja 和 Orlin 在1987年提出了Improved SAP 算法，它充分利用了距离标号的作用，每次发现顶点无出弧时不是像 Dinic 算法那样到最后进行 BFS，而是就地对顶点距离重标号，这样相当于在遍历的同时顺便构建了新的分层网络，每轮寻找之间不必再插入全图的 BFS 操作，极大提高了运行效率。国内一般把这个算法称为 SAP...显然这是不准确的，毕竟从字面意思上来看 E-K 和 Dinic 都属于 SAP，我还是习惯称为 ISAP 或改进的 SAP 算法。

与 Dinic 算法不同，ISAP 中的距离标号是每个顶点到达终点 t 的距离。同样也不需显式构造分层网络，只要保存每个顶点的距离标号即可。程序开始时用一个反向 BFS 初始化所有顶点的距离标号，之后从源点开始，进行如下三种操作：(1)当前顶点 i 为终点时增广 (2) 当前顶点有满足 dist[i] = dist[j] + 1 的出弧时前进 (3) 当前顶点无满足条件的出弧时重标号并回退一步。整个循环当源点 s 的距离标号 dist[s] >= n 时结束。对 i 点的重标号操作可概括为 dist[i] = 1 + min{dist[j] : (i,j)属于残量网络Gf}。具体算法描述如下：

algorithm Improved-Shortest-Augmenting-Path
 1 f <-- 0
 2 从终点 t 开始进行一遍反向 BFS 求得所有顶点的起始距离标号 d(i)
 3 i <-- s
 4 while d(s) < n do
 5     if i = t then    // 找到增广路
 6         Augment
 7         i <-- s      // 从源点 s 开始下次寻找
 8     if Gf 包含从 i 出发的一条允许弧 (i,j)
 9         then Advance(i)
10         else Retreat(i)    // 没有从 i 出发的允许弧则回退
11 return f

procedure Advance(i)
1 设 (i,j) 为从 i 出发的一条允许弧
2 pi(j) <-- i    // 保存一条反向路径，为回退时准备
3 i <-- j        // 前进一步，使 j 成为当前结点

procedure Retreat(i)
1 d(i) <-- 1 + min{d(j):(i,j)属于残量网络Gf}    // 称为重标号的操作
2 if i != s
3     then i <-- pi(i)    // 回退一步

procedure Augment
1 pi 中记录为当前找到的增广路 P
2 delta <-- min{rij:(i,j)属于P}
3 沿路径 P 增广 delta 的流量
4 更新残量网络 Gf

算法中的允许弧是指在残量网络中满足 dist[i] = dist[j] + 1 的弧。Retreat 过程中若从 i 出发没有弧属于残量网络 Gf 则把顶点距离重标号为 n 。

虽然 ISAP 算法时间复杂度与 Dinic 相同都是 O(V2E)，但在实际表现中要好得多。要提的一点是关于 ISAP 的一个所谓 GAP 优化。由于从 s 到 t 的一条最短路径的顶点距离标号单调递减，且相邻顶点标号差严格等于1，因此可以预见如果在当前网络中距离标号为 k (0 <= k < n) 的顶点数为 0，那么可以知道一定不存在一条从 s 到 t 的增广路径，此时可直接跳出主循环。在我的实测中，这个优化是绝对不能少的，一方面可以提高速度，另外可增强 ISAP 算法时间上的稳定性，不然某些情况下 ISAP 会出奇的费时，而且大大慢于 Dinic 算法。相对的，初始的一遍 BFS 却是可有可无，因为 ISAP 可在循环中自动建立起分层网络。实测加不加 BFS 运行时间差只有 5% 左右，代码量可节省 15~20 行。

来自CSDN的keshuqi

网络流－最大流问题 ISAP 算法解释（转自Renfei Song's Blog）

网络流－最大流问题 ISAP 算法解释
August 7, 2013 / 编程指南

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;
}

来自CSDN的TA201314

先来学算导。。
26.1流网络
①与第二版不同，算导第三版修改了关于流的约定，约定所有流均非负；并删除了流的第二条性质（反向对称）。同时放弃了给源点以定义，允许源点入度不为0，相应的，也允许汇点出度不为0.
②一些奇怪的性质包括：流守恒，所有点必在s~>t上，两点之间最多只有一条边。
③在网络流中遇到了大量加边加点的技巧：多个源汇点的话可以添加一个超级源汇点；为反向平行边添加虚拟节点，将一条边裂解成两个。
④网络中的流为一个以0为下界的凸集。
26.2Ford-Fulkerson方法
①换取等价增加：在残留网络中寻找从s->t的简单路径。（在流网络中，任意从s到t的简单路径都是一个流。）
②等价关系：
最大流<=>在残留网络中无可增广路<=>最小割。
③EK的时间复杂度证明：
每个可增广路至少有一条关键边；
若约定每条边权重为1，则EK时在残留图中任意节点到s点的最短路径必不下降。
④26.2-11无向图的边联通性是指使图变为非连通图所需要删除的最少边数k。
即，最小割！
·转换成网络流问题：O(V^2E)
⑤26.2-13在一个满足一个最小的前提小，找另一个最小？在保证一个不变的前提下扩大另一个的区别！
假定我们希望找到一个流网络G的所有最小切割中包含边的条数最少的切割，这里假定G的所有容量都是整数值。说明如何修改Ｇ的容量来创建一个新的流网络G',使得G‘中的任意一个最小切割是G中包含边的条数最少的最小切割。
以Ec(u,v)+1代c(u,v)即可。
乘E不会改变相对大小，但却会放大相对差距；+1会使得边数不同则割的容量不同。
ISAP算法：
①一个几乎完全类似的性质：每次增广之后，任意节点到汇点的距离也必然不下降。
推论：若一条路径在某次增广之前是最短路径，则在其增广之后只有两种可能，要么不存在了，要么依然为最短路径（因为不会有更短的路径被诞生）。
所以呢？！我们大可为每个节点保存一个下界函数d(i),作为i到t的最短距离的下界；在每次增广时，用接近汇点优先的策略寻找是否存在符合d的最短路。
对于每个节点，我们还需要保存它当前处理到哪条边了，因为d(i)是时间t的不下降函数，所以若在当前边之前的所有边在之前的遍历中已经不可能有允许拓展，则在当前状态中亦然不可能有允许拓展——之所以这样做，实际上是保证了其上界时间复杂度。
②如果在某一个节点上出现了当前d(i)过大的情况，就执行retreated，寻找最小后继，更新d(i)
——这样之所以可行是因为其更新的先后顺序决定的，若d(i)被过大，则必然是由于更新出该值的后继节点被更新导致。
——retreated操作保证了d(i)的紧密递增，一个点的d(i)不可能比其后继节点的d大2,即其必是从s开始的一棵下降连续或断裂上升的树；而不可能出现断裂下降的情况。
③GAP优化：
我们保存d(i)的每个值在当前G中的数量，若在某次retreated之后，若在残留图G中不存在与原先的d(i)相等的值，则必然意味着s与t不相连了，即我们找到了最小割、最大流。

GAP优化的必要性：如果不使用GAP优化的话，注意在最小割被找出之后实际上我们需要再遍历一整遍s所在的割并依次retreated每个点使得s最终被retreated，而这样做几乎相当于遍历一遍邻接表，这显然不是我们想要的，所以我们需要使用GAP优化快速判断出是否s与t已被割开。
④时间复杂度分析：
每条边至多成为关键边|V|/2次，因为关键边是会消失的边，在它消失之后即其完全成为了其原先的反向边，此时若其再成为关键路径，就意味着由d(u)=d(v)+1,变成了d'(u)=d'(v)-1,即d'(u)<=d(u)-2,因此增广次数至多为|E||V|/2次。
下面我们来讨论增广时的开销，首先注意到寻找增广路的过程是需要把点集至多遍历一遍的，这样做的花费是|V|.
然而需要注意到我们还进行了遍历邻接表的操作，我们试图均摊这样做的开销：对于每个节点v，其最多能直接抵达|V|个节点，那么出边也最多会被遍历|V|次，因为d(i)的取值范围是[1,|V|],而对于每个会出现的d(i),都需要遍历且仅遍历一遍其所有出边。
综上，ISAP算法让我们看到了下界函数在使不下降操作中的（贪心）作用，这不仅仅是对于网络流的ISAP算法，对于其他题来说也是很有用处的。

代码 来自csdn的一样菜

#include"stdio.h"
#include"string.h"
#include"queue"
#include"stack"
#include"iostream"
#include"stdlib.h"
#define inf 99999999
#define M 50000
using namespace std;
struct st
{
    int u,v,w,next;
}edge[M];
int t,head[M],cur[M],q[M],gap[M],dis[M];
void init()
{
    t=0;
    memset(head,-1,sizeof(head));
}
void add(int u,int v,int w)
{
    edge[t].u=u;
    edge[t].v=v;
    edge[t].w=w;
    edge[t].next=head[u];
    head[u]=t++;

    edge[t].u=v;
    edge[t].v=u;
    edge[t].w=0;
    edge[t].next=head[v];
    head[v]=t++;
}
void bfs(int start,int endl)//建立到汇点的距离层次图存在dis[]数组中
{
    int rear=0,i,j;
    memset(dis,-1,sizeof(dis));
    memset(gap,0,sizeof(gap));//gap[x]记录dis[i]=x出现了多少次
    dis[endl]=0;
    gap[dis[endl]]=1;
    q[rear++]=endl;
    for(i=0;i<rear;i++)
    {
        for(j=head[q[i]];j!=-1;j=edge[j].next)
        {
            int v=edge[j].v;
            if(dis[v]==-1)
            {
                ++gap[dis[v]=dis[q[i]]+1];
                q[rear++]=v;
            }
        }
    }
}
int SAP(int start,int endl,int n)
{
    int ans=0;
    bfs(start,endl);
    int cur[M];//代替head数组
    memcpy(cur,head,sizeof(head));
    int stack[M],top=0;//建立手工栈
    int u=start,i;
    while(dis[start]<n)
    {
        if(u==endl)//当搜到终点时即找到从原点到汇点的增光路，正常处理即可
        {
            int mini=inf,tep;
            for(i=0;i<top;i++)
            {
                if(mini>edge[stack[i]].w)
                {
                    mini=edge[stack[i]].w;
                    tep=i;
                }
            }
            for(i=0;i<top;i++)
            {
                edge[stack[i]].w-=mini;
                edge[stack[i]^1].w+=mini;
            }
            ans+=mini;
            top=tep;
            u=edge[stack[top]].u;//此时的u为变容量为0的u
        }
        if(dis[u]&&gap[dis[u]-1]==0)//出现了断层，没有增广路
            break;
        for(i=cur[u];i!=-1;i=edge[i].next)//遍历与u相连的未遍历的节点
        {
            int v=edge[i].v;
            if(dis[v]!=-1)
            {
                if(edge[i].w&&dis[u]==dis[v]+1)//层次关系找到允许路径
                    break;
            }
        }
        if(i!=-1)//找到允许弧
        {
            cur[u]=i;
            stack[top++]=i;
            u=edge[i].v;
        }
        else//无允许的路径，修改标号 当前点的标号比与之相连的点中最小的多1
        {
            int mini=n;
            for(i=head[u];i!=-1;i=edge[i].next)
            {
                if(edge[i].w==0)continue;
                int v=edge[i].v;
                if(mini>dis[v])//找到与u相连的v中dep[v]最小的点
                {
                    mini=dis[v];
                    cur[u]=i;//最小标号就是最新的允许弧
                }
            }
            --gap[dis[u]];//dep[u] 的个数变化了 所以修改gap
            ++gap[dis[u]=mini+1];//将dep[u]设为min(dep[v]) + 1, 同时修改相应的gap[]
            if(u!=start)//该点非源点&&以u开始的允许弧不存在，退点
                u=edge[stack[--top]].u;
        }
    }
    return ans;
}
int main()
{
    int n,m;
    while(scanf("%d%d",&m,&n)!=-1)
    {
        init();
        while(m--)
        {
            int a,b,c;
            scanf("%d%d%d",&a,&b,&c);
            add(a,b,c);
        }
        int ans=SAP(1,n,n);
        printf("%d\n",ans);
    }
}