安塔的ISAP和HLPP学习笔记
木木！
2019-09-25 16:50:16
# 前置知识
关于这些知识，可能会提及，但是不会详细介绍。请确认你已掌握以下知识：
+ [Dinic算法](https://www.cnblogs.com/linzhengmin/p/9313216.html)
+ 网络流的基本用语，如残量网络
+ [在`priority_queue`里面使用自定义`cmp`函数](https://blog.csdn.net/qq_30339595/article/details/79566788)、
# ISAP
`ISAP`，全称`Improved Shortest Augment Path`。这里提到它`Improved`，指它不需要在残量网络上重复执行`bfs`。它的基本思路是，一次`BFS`确定标号之后边流边修改。
由于要修改，而对于一个点，我们可以很容易地枚举它的出边，但是除非存反图，没有办法枚举入边。所以标号不再是到`s`的最短距离，而是到`t`的最短距离。这样可以很方便地进行修改。
首先做一次类似`Dinic`的计算。如果传进来的`flow`还有剩余，就需要提高标号以便于继续找增广路。由于标号的范围可以认为是`1-n`的，所以可以直接给标号`+1`。
按照这个思路，我们写出来的`isap`类似于一个`while`套`Dinic`。大概长这样：
（递归版）
```cpp
int isap(int s,int flow)
{
if(s==t)
{
return flow;
}
int flowed = 0;
while(flow && tag[s]<n)
{
for(int p=beg[s]; p; p=nxt[p])
{
if(ci[p] && tag[ed[p]]==tag[s]-1)
{
int f = isap(ed[p],min(ci[p],flow));
flow -= f;
flowed += f;
ci[p] -= f;
ci[p^1] += f;
if(!flow)
{
return flowed;
}
}
}
++tag[s];
}
return flowed;
}
```
这个就是`isap`的基本骨架了。当然，`Dinic`的运行速度有一半都是它的优化贡献的。`isap`同样也需要有很多优化。
## 1.优化-当前弧优化
这里的`ISAP`同样可以加当前弧优化，和`Dinic`的一样。
这里要注意一点：在`++tag[s]`的时候，必须将当前弧还原（`cnt[s]=beg[s]`）。
## 2.优化-gap优化
> 如果标号出现断层，说明当前已经不存在一个能够从s到t的增广路了。此时可以直接结束算法。
——这句话是对非递归版`ISAP`说的。
对于递归版`ISAP`，同样可以维护一个标号数量数组`tagn[]`，每次更改标号时候，如果对应标号被改没了，可以直接`break`。但是要注意，不能直接结束算法，因为前面可能还有点能补上来。
由于比较抽象，附代码。
```cpp
int isap(int s,int flow)
{
if(s==t)
{
return flow;
}
if(!tag[s] || tag[s]>n) //小优化
{
return 0;
}
int flowed = 0;
while(flow && tag[s]<n)
{
for(int p=cnt[s]; p; p=nxt[p])
{
if(ci[p] && tag[ed[p]]==tag[s]-1)
{
int f = isap(ed[p],min(ci[p],flow));
flow -= f;
flowed += f;
ci[p] -= f;
ci[p^1] += f;
if(!flow)
{
cnt[s] = p;
return flowed;
}
}
}
--tagn[tag[s]];
if(!tagn[tag[s]]) //直接退出
{
++tag[s];
++tagn[tag[s]];
return flowed;
}
++tag[s];
++tagn[tag[s]];
cnt[s] = beg[s];
}
return flowed;
}
```
## 3.修正-vis优化
如果你直接拿这份代码去交[P3376模板](https://www.luogu.org/problem/P3376)，你会发现你`TLE`了一个点。
因为访问一个点之后，它可能没有任何高度合适的点，所以它就会将`tag`自增，然后再访问下一个点。下一个点可能也没有任何点，因此它也会将`tag`自增。由此，我们发现这个过程可以一直重复。更重要的是，其中的某一个点可能在`tag`自增后比中间还在递归栈内的某一个点的`tag`大，因此可能会形成环。
由于我们可能会重复访问一个点多达$\Theta(n)$次，这个增广路的长度最大是$\Theta(n^2)$的，因此你就会华丽丽地`TLE`。
修正方法是，增加一个`vis`数组，以避免重复访问某一个点导致增广路过长。
然后就是最终版代码：
```cpp
#include <queue>
#include <cstdio>
using namespace std;
int beg[10005];
int ed[200005];
int ci[200005];
int nxt[200005];
int top = 1;
void addedge(int a,int b,int c)
{
++top;
ed[top] = b;
ci[top] = c;
nxt[top] = beg[a];
beg[a] = top;
}
int n,t;
int tag[10005];
int cnt[10005];
int vis[10005];
int tagn[20005];
void rbfs()
{
for(int i=1; i<=n; ++i)
{
cnt[i] = beg[i];
}
tag[t] = 1;
tagn[1] = 1;
queue<int> q;
q.push(t);
while(!q.empty())
{
int th = q.front();
q.pop();
for(int p=beg[th]; p; p=nxt[p])
{
if(!ci[p] && !tag[ed[p]])
{
tag[ed[p]] = tag[th]+1;
++tagn[tag[ed[p]]];
q.push(ed[p]);
}
}
}
}
int isap(int s,int flow)
{
if(s==t)
{
return flow;
}
if(!tag[s] || tag[s]>n)
{
return 0;
}
vis[s] = 1;
int flowed = 0;
while(flow && tag[s]<n)
{
for(int p=cnt[s]; p; p=nxt[p])
{
if(ci[p] && !vis[ed[p]] && tag[ed[p]]==tag[s]-1)
{
int f = isap(ed[p],min(ci[p],flow));
flow -= f;
flowed += f;
ci[p] -= f;
ci[p^1] += f;
if(!flow)
{
cnt[s] = p;
vis[s] = 0;
return flowed;
}
}
}
--tagn[tag[s]];
if(!tagn[tag[s]])
{
return flowed;
}
++tag[s];
++tagn[tag[s]];
cnt[s] = beg[s];
}
vis[s] = 0;
return flowed;
}
int main()
{
int m,s;
scanf("%d%d%d%d",&n,&m,&s,&t);
for(int i=1; i<=m; ++i)
{
int a,b,c;
scanf("%d%d%d",&a,&b,&c);
addedge(a,b,c);
addedge(b,a,0);
}
rbfs();
printf("%d",isap(s,0x7f7f7f7f));
}
```
# HLPP
这里讲解`HLPP`的方式和途径和别的博客不一样。
显然，我们并没有将`ISAP`的所有优化搞完。如果同一个点被调用多次，我们可以将这些调用合并起来。即，`isap(s,3)`连着`isap(s,5)`等价于`isap(s,8)`。
合并调用，类似于线段树的合并操作，可以用`lazy_algorithm`的思想。我们开一个队列，存下那些调用，而不是每一次都着急地运行这些调用。
当然，我们也不能和上面一样使用`isap`的返回值了。不过这不是大问题。
```cpp
void hlpp(int s)
{
flow[s] = 0x7f7f7f7f7f7f7f7f;
queue<int> q;
q.push(s);
while(!q.empty())
{
int th = q.top();
q.pop();
while(flow[th] && tag[th]<n)
{
for(int p=beg[th]; p&&flow[th]; p=nxt[p])
{
if(ci[p] && tag[ed[p]]==tag[th]-1)
{
int f = min((long long)ci[p],flow[th]);
flow[ed[p]] += f;
ci[p] -= f;
flow[th] -= f;
ci[p^1] += f;
if(ed[p]!=t)
{
pq.push(ed[p]);
}
}
}
if(flow[th])
{
++tag[th];
}
}
}
}
```
`flow[]`数组就相当于调用`isap`时的参数`flow`。这里使用的是没有加任何优化的`isap`形式，得到了`hlpp`的骨架。当然，这时候它还不能被称作`HLPP`，因为缺了个最关键的优化。
我们约定，对于一个点做一次类似于`isap`的操作称为“推流”，而在那个类似`isap`的过程中给一个点的`flow[]`增加称为“推向”
## 1.优化-最高标号优化
> 预流推进+最高标号优化=最高标号预流推进（HLPP）
由于我们应用了`lazy_algorithm`，我们可以自由地决定先对哪些点推流。显然，如果可以的话，我们希望对标号尽量高的点推流。因为这样的话，可以尽量快地让一个点达成`tag[s]=n+1`。（“迫害原则”）
当然，你可以决定先对`flow[]`大的点推流，也可以决定先对`tag[]`小的点推流。你也可以对各种方案都试试，然后选择效果最好的。
Cheriyan和Maheshwari在1988年证明了先对`tag[]`大的点推流的复杂度为$O(n^2\sqrt{m})$，所以，我们只讲先对`tag[]`大的点推流的做法。
另外一个理由就是，先对`tag[]`大的点推流，就是先对远离`t`的点推流。这样，就像挤牙膏一样，把流都往`t`推。
至于实现原理，只需要使用`priority_queue`就好了。
## 2.优化-inq优化
类似`SLF`优化的`SPFA`，一个点只需要进队一次，所以可以记录`inq`数组，防止重复进队。
## 3.优化-当前弧优化
与`ISAP`的当前弧优化差不多。由于我们实现的是加`lazy_algorithm`的`ISAP`，当前弧优化可以完美兼容。
## 4.优化-gap优化
这里，出现断层仍旧不可以直接结束代码。但是，如果出现断层，可以确定`tag`比断层大的所有点都不能推流到`t`了。所以，可以直接将他们的`tag[]`变成`n+1`，来表示这个点已经废了。
附HLPP代码：
```cpp
#include <queue>
#include <cstdio>
using namespace std;
int beg[10005];
int ed[240005];
int nxt[240005];
int ci[240005];
int top = 1;
void addedge(int a,int b,int c)
{
++top;
ed[top] = b;
ci[top] = c;
nxt[top] = beg[a];
beg[a] = top;
}
int n,t;
int tag[10005];
int inq[10005];
int cnt[10005];
int tagn[10005];
long long flow[10005];
void rbfs()
{
for(int i=1; i<=n; ++i)
{
cnt[i] = beg[i];
}
tag[t] = 1;
queue<int> q;
q.push(t);
tagn[1] = 1;
while(!q.empty())
{
int th = q.front();
q.pop();
for(int p=beg[th]; p; p=nxt[p])
{
if(!tag[ed[p]])
{
tag[ed[p]] = tag[th]+1;
++tagn[tag[ed[p]]];
q.push(ed[p]);
}
}
}
}
class cmp
{
public:
bool operator() (int a,int b)
{
return tag[a]<tag[b];
}
};
void hlpp(int s)
{
flow[s] = 0x7f7f7f7f7f7f7f7f;
priority_queue<int,vector<int>,cmp> pq;
inq[s] = 1;
pq.push(s);
while(!pq.empty())
{
int th = pq.top();
inq[th] = 0;
pq.pop();
while(flow[th] && tag[th]<n)
{
for(int p=cnt[th]; p; p=nxt[p])
{
if(ci[p] && tag[ed[p]]==tag[th]-1)
{
int f = min((long long)ci[p],flow[th]);
flow[ed[p]] += f;
ci[p] -= f;
flow[th] -= f;
ci[p^1] += f;
if(ed[p]!=t && !inq[ed[p]])
{
inq[ed[p]] = 1;
pq.push(ed[p]);
}
if(!flow[th])
{
cnt[th] = p;
break;
}
}
}
if(flow[th])
{
--tagn[tag[th]];
if(!tagn[tag[th]])
{
for(int i=1; i<=n; ++i)
{
if(tag[i]>=tag[th])
{
tag[i] = n+1;
}
}
break;
}
++tag[th];
++tagn[tag[th]];
cnt[th] = beg[th];
}
}
}
}
int main()
{
int m,s;
scanf("%d%d%d%d",&n,&m,&s,&t);
for(int i=1; i<=m; ++i)
{
int a,b,c;
scanf("%d%d%d",&a,&b,&c);
addedge(a,b,c);
addedge(b,a,0);
}
rbfs();
hlpp(s);
printf("%lld\n",flow[t]);
}
```
# 总结
ISAP和HLPP在[P3376](https://www.luogu.org/problem/P3376)下跑速差不多，`HLPP`加那三个优化会稍快一些。`ISAP`比`HLPP`好写一些。在[P4722](https://www.luogu.org/problem/P4722)，`HLPP`几乎是卡着最坏复杂度过的，`ISAP`根本过不了。
具体而言，在网络流算法里面，`HLPP`几乎肯定是首选。如果不想码太多代码，可以选择`ISAP`。如果还不想码太多代码，可以选择`Dinic`。一般而言，`Dinic`已经足够。
# 参考资料
感谢以下作者的无私付出：
+ [究级的最大流算法：ISAP与HLPP -钱逸凡](https://www.luogu.org/blog/ONE-PIECE/jiu-ji-di-zui-tai-liu-suan-fa-isap-yu-hlpp)
+ [最大流算法-ISAP -permui](https://www.cnblogs.com/owenyu/p/6852664.html)