[学习笔记] 网络最大流的HLPP算法

#define (u)的伴点集合 与(u)相隔一条边的且(u)能达到的点的集合

(0x00~ {}~Preface)

(HLPP(Highest~Label~Preflow~Push))最高标签预流推进算法是处理网络最大流里两种常用方法——增广路&预流推进中，预流推进算法的一种。据传由(tarjan)发明怎么又是他 ，并被其他科学家证明了其复杂度是紧却的(O(n^2sqrt m))。在随机数据中不逊色于普通的增广路算法，而在精心构造的数据中无法被卡，所以是一种可以替代(Dinic)的方法（随我怎么说，代码又长又难调，所以还是(Dinic)好啊( m{TAT})）

但无论怎样，(wiki)里面已经承认(HLPP)是现在最优秀的网络流算法了。

那么预流推进这个大门类里面，思想都差不多。大抵上就是我们对每个点记录超额流((Extra~Flow)) ，即允许流在非源点暂时存储，并伺机将超额流推送出去。不可推送的，就会流回源点。那么最终答案显然存储在(Extra[T])里面。

但同时这也有一个问题，就是会出现两个点相互推送不停的情况。为了防止这样，我们采用最高标号的策略，给每个点一个高度，对于一个点(u)以及它的伴点集合({v})，当且仅当(h_u = h_v + 1) 时才可以推送流。并且我们对于源点(S)，设置(h_S = N)，并对于(S)实行无限制推送。那么最后的答案就保存在(Extra[T])里面 。

但有时，我们发现有个点是”谷“，即周围点的高度都比它高，但是它有超额流。那么我们此时考虑拔高它的高度，即重贴标签((relabel))操作。

(0x01quad) 初步的算法流程

以下我们用(Extra_u)表示(u)的超额流，(h_u)表示(u)的高度，用(f_k)表示边(k)的容量。

• 首先把所有的(h_i)都置为零，并把(h_s)置为(N)（点数）。

• 将(S)的流推送到每个与(S)相邻的点，同时把他们加入一个以高度为键值得大根堆，所以每次取出的应该是高度最高的、且超额流不为零的点，并执行推送操作。

• 对于点(u)推送过程中，如果(Extra_u)减到了(0)，就立即退出（优化一）

• 对于每条出边(k)，推送的流量(F = min(f_k,Extra_u))并执行两个点（(u,v)）的超额流增减。如果(v)不在堆里面，要把(v)放到堆里面。

• 如果推送完毕(Extra[u])不为零，那么从他的伴点集合选取一个高度最小的点并记录它的高度(h_{min})，则新的(h_u = h_{min}+1)，并把(u)入堆。

好的，然后就可以撒花了……可是等等，他怎么这么慢(qaq)

接下来我们发现，重贴标签的过程似乎与(ISAP)有点点像……所以我们不妨通过一个(Gap)数组来记录”断层情况“：即如果对于一个点(u)来说，他的伴点集({v})已经不存在(h_u = h_v + 1)的点了，并且也不存在一个点(j)使得(h_j = h_u)那么这个地方就是一个断层((Gap)) ，那么也就是说，对于所有(h_i> h_u)的点来说，它们把流推送到(h_u)的高度就不能继续推送了，所以我们直接(h_i = N + 1)，让他们回流到源点。（优化二）

接下来这个优化，亲测可以提速(4000ms)，平均每个测试点提速(700) ~ (800ms)，去掉数据最小的点，每个点平均提速(1000ms)。这就是——(BFS)!

我们不妨一开始就倒着(BFS)一遍，搜出每个点离汇点的最短距离作为初始高度而不是把零作为初始高度（源点高度还是(N)。嗯，(Mr\_Spade)大佬实在太强了(qwq)

对了，代码实现方面，需要好多判断不是源点和汇点的小细节……无路赛无路赛无路赛(>\_<)！

(color{red}{C}color{cyan}{o}color{gold}{d}color{green}{e}·1)

#include <bits/stdc++.h>
//省略某些部分
#define Inf, MAXN, MAXM, to(k)
struct state{
    int num, h ;
    bool operator <(const state & now) 
    const{  return h < now.h ; }
} ; priority_queue <state> heap ; 
BFS init ;
int N, M, S, T, cnt = -1, A, B, C, D, t, min_h ;
int head[MAXN], Extra[MAXN], H[MAXN], Gap[MAXN], node ;

inline void Preflow_Push(){
    register int i, k ;
    for (i = 1 ; i <= N ; ++ i) 
        if(H[i] < Inf) ++ Gap[H[i]] ;
    for(k = head[S]; k != -1 ; k = E[k].next)
        if((t = E[k].f)){
            E[k].f -= t, E[k ^ 1].f += t, Extra[S] -= t, Extra[to(k)] += t ;
            if(to(k) != T && !vis[to(k)])
                heap.push((state){to(k), H[to(k)]}), vis[to(k)] = 1 ;
        }
    while(!heap.empty()){
        vis[node = heap.top().num] = 0, min_h = Inf, heap.pop() ;
        for(k = head[node] ; k != -1 ; k = E[k].next){
            if(E[k].f && H[node] == H[to(k)] + 1){
                t = min(Extra[node], E[k].f) ;
                E[k].f -= t, E[k ^ 1].f += t, Extra[node] -= t, Extra[to(k)] += t ;
                if(!vis[to(k)] && to(k) != S && to(k) != T)
                    heap.push((state){to(k), H[to(k)]}), vis[to(k)] = 1 ;
            }
            if (E[k].f) min_h = min(min_h, H[to(k)]) ;
            if (!Extra[node]) break ;
        }
        if(Extra[node]) {
            if (!--Gap[H[node]])    
                for(i = 1; i <= N ; ++ i)
                    if(i != S && i != T && H[i] > H[node] && H[i] < N + 1) H[i] = N + 1 ;
            H[node] = Inf; H[node] = min_h + 1 ; 
            heap.push((state){node, H[node]}), vis[node] = 1, ++ Gap[H[node]] ;
        }
    }
}
inline int read() {fast;}
int main(){
    Input() ;
    for (i = 1 ; i <= N ; ++ i) 
        head[i] = -1, H[i] = Inf ;
    while(M --){Add; }
    q.push(T), H[T] = 0 ;
    while(!q.empty()){
        int now = q.front() ; q.pop() ;
        for(k = head[now] ; k != -1 ; k = E[k].next)
            if (H[to(k)] > H[now] + 1)
                H[E[k].to] = H[now] + 1, q.push(E[k].to) ;
    }
    if (H[S] == 0) {cout << 0 << endl ; return 0 ;} 
    H[S] = N, Preflow_Push() ; cout << Extra[T] << endl ;
}

看起来我们加上下面这一句话的毒瘤卡常就可以有(4000ms)左右的好成绩，但事实上，这个速度依旧慢的要死。

注意！这个写法是经典写法，其时间复杂度是紧却的(oldsymbol{ m{Theta(n^2mlogn)}})的，也就是说在(oldsymbol{n=1200})时会带一个(oldsymbol{10})倍的常数

怎么优化呢——

(oldsymbol{0x02~~Optimization})

各位，你们将会见到迄今为止OI界最丧心病狂的优化（之一）……

来，我们首先思考思考普通的HLPP到底会慢在哪里：

• (STL)支持的(heap)（比如优先队列）实在是太太太…太慢了！

• 每次(Gap)优化，我们的时间复杂度是紧确的(Theta(n))。的这显然不合算，因为假设我当前的(oldsymbol{gap})（断层）正好位于倒数第一高的点和倒数第二高的点，那么也就相当于我单次会浪费(oldsymbol{Theta(n)})的时间。

事实上…普通的(HLPP)代码并不长，主要问题就是这两个。

我们考虑，如果不用堆的话怎么做呢？

呃…不用堆的意思并不是我们不从高度最大的点开始推送。这个地方需要一个(idea)——在(HLPP)中，高度函数(oldsymbol{H(x)})和点数集大小(oldsymbol{N(x)})是广义同阶的。 换句话说，我们可以考虑从高度入手。

换句话说，我们原来是通过节点编号访问节点以及其高度，现在我们如果从高度入手，再去访问节点，我们就可以做到(oldsymbol{O(n)})而不是(oldsymbol{ m{O(nlogn)}}) 。 那么由于同一高度的节点或许有很多，直接开一个(vector)。在这个地方我们用(vector)而不用二维数组建立二维关系的原因，主要是我们初始化麻烦得很，如果套用(memset)或者(fill)的话，常数之大可想而知。

那么这两个问题就顺理成章地解决了。但这个地方还有一个优化，就是虽然(vector)和(list)都是线性容器，但是(list)的本质是双向链表，频繁处理插入删除操作时会具有更优秀的表现。

也就是说，原来的(Gap)数组我们可以直接用(list)做，以图更小的常数。那么这时存在一个问题，就是虽然本质上删除是容易的，但是你怎么知道要删同一高度下的哪个元素(=@__@=)？就算你知道，(list)也不知道啊2333

hhh不皮了，其实我们记录一下位置就好，即记录一下每个节点在(list)中的位置，单独开一个(Iterator)类型的(vector)记录即可。

好了，现在我们获得了(10)倍(+)的常数优势qwq，撒花花…

哦对，还有几点我debug的时候被坑死的点：

• 那个(Iterator)类型的(vector)对象是点的编号不是高度！
• 注意你的下标！下标！再说一遍，下标！因为STL自带左闭右开的性质wrnm，所以一定要注意，如果你是([1,n])选手，注意你的(assign)函数！

(color{red}{C}color{cyan}{o}color{gold}{d}color{green}{e}·2) (我觉得写的很难看但是有注释qaq)：

//writter:Orchidany(pks)
#include <bits/stdc++.h>
#pragma GCC target("avx")
#pragma GCC optimize(3)
#pragma GCC optimize("Ofast")//sb毒瘤优化

#define MAXN 10030
#define min my_min
#define BG begin()
#define gc getchar
#define rr register 
#define Iter iterator
#define INF 2147483647
#define rep(i, a, x) for(i = a ; i <= x ; ++ i)

using namespace std ;
typedef list<int> List ; int step;

struct Edge{
    int to, f, next ;
    Edge(int to,int f,int next):to(to),f(f),next(next){}//没有人发现正下方这句注释前半句和后半句都是三个音节的吗qaq
} ; vector <int> q, H, Extra, Set[MAXN], cnt ; List Gap[MAXN] ;//list，就是快（
//q:队列，H:高度，Extra：每个点的超额流，Set:…就是那个经典版HLPP里的堆，高度做第一维
int Ans, N, M, S, T, max_H, now_H ; vector <Edge> E[MAXN] ; /*vector存边（据说会快）*/vector<List::iterator> Era_pos ; //辅助定位+删除

inline void eggs() { ;}//for free~
inline int my_min(int a, int b){return a & ((a - b) >> 31) | b & ( ~ (a - b) >> 31) ;}//黑科技
inline void Add(int f, int v, int u){ E[u].push_back(Edge(v, f, E[v].size())), E[v].push_back(Edge(u, 0, E[u].size() - 1)) ; }
inline int qr(){ rr int k = 0 ; char c = gc() ; while (!isdigit(c)) c = gc() ;while (isdigit(c)) k = (k << 1) + (k << 3) + c - 48, c = gc() ; return k ; }

inline void Init_label(){//等价于一开始的那个BFS，只执行一次
    rr int i, h = 0, t = 0 ;q.clear(), q.resize(N) ; 
    H.assign(N + 1, N + 1) ; H[T] = 0 ; q[t ++] = T ;//从T（高度小的）向前标号
    while (h < t){//队列……BFS……真熟悉啊……嗝……
        rr int now = q[h] ; ++ h ;
        for (vector <Edge> :: Iter k = E[now].begin() ; k != E[now].end() ; ++ k)
            if (H[k->to] == N + 1 && E[k->to][k->next].f) H[k->to] = H[now] + 1, ++ cnt[H[k->to]], q[t ++] = k->to ;
    }
    rep(i, 0, N + 1) Set[i].clear(), Gap[i].clear() ;//还是清空一下比较好吧
    rep(i, 0, N) 
        if (H[i]  <N + 1)  
            Era_pos[i] = Gap[H[i]].insert(Gap[H[i]].BG, i), //疑似insert函数的返回值是一个指针qaq
            (Extra[i]>0) ? Set[H[i]].push_back(i) : eggs() ;//这个彩蛋（eggs）是因为，三目运算符":"两边类型需要形同…
    max_H = now_H = H[q[-- t]] ; //更新，BFS的性质，最后一个元素一定高度最大（除了源点）
}
inline void Push(int x, Edge &e){//单独写出来的push函数，好像很方便？
    rr int now_flow = min(Extra[x], e.f) ;
    Extra[x] -= now_flow, e.f -= now_flow, Extra[e.to] += now_flow, E[e.to][e.next].f += now_flow ;
    if (Extra[e.to] > 0 && Extra[e.to] <= now_flow) Set[H[e.to]].push_back(e.to) ;  // push it into "heap"
}
inline void _Push(int x){
    rr int i, x_h = N, t = H[x] ; 
    for (vector <Edge> :: Iter k = E[x].BG ; k != E[x].end() ; ++ k)
        if (k->f > 0)//如果可以流
            if (H[k->to] == H[x] - 1) { Push(x, *k) ; if (!Extra[x]) return ;} else x_h = min(x_h, H[k->to] + 1) ;
    if (cnt[H[x]] <= 1){//如果出现断层了
        for(i = t ; i <= max_H ; Gap[i].clear(), ++ i)//这个gap的for肯定比O(n)优秀
            for(List::Iter k = Gap[i].BG ; k != Gap[i].end() ; ++ k) cnt[H[*k]] --, H[*k] = N ; 
        max_H = t - 1 ; /*断层以上的高度都没用了*/return ;
    }
    -- cnt[t], Era_pos[x] = Gap[t].erase(Era_pos[x]) ; H[x] = x_h ; if (x_h == N) return ; //重贴标签操作，为当前点删除原来的高度
    ++ cnt[x_h], Era_pos[x] = Gap[x_h].insert(Gap[x_h].begin(), x), max_H = max(now_H = x_h, max_H), Set[x_h].push_back(x) ;//增添新的高度
}
inline int HLPP(){
    rr int i, now ; H.assign(N, 0) ; H[S] = N ; Era_pos.resize(N)  ;
    rep(i, 0, N - 1) if (i != S) Era_pos[i] = Gap[H[i]].insert(Gap[H[i]].BG, i) ; 
    cnt.assign(N, 0), cnt[0] = N - 1 ; Extra.assign(N, 0), Extra[S] = INF, Extra[T] =- INF ;
    rep(i, 0, E[S].size() - 1) Push(S, E[S][i]) ;  //下面源点要单独拿出来推送，因为源点推送时高度差不需要=1.
    Init_label() ; //初始化（BFS）
    while (now_H >= 0) //正式开始HLPP（泪目）
        if (Set[now_H].empty()) now_H -- ; //高度递减，实现一个堆的效果
        else now = Set[now_H].back(), Set[now_H].pop_back(), _Push(now) ;
    return Extra[T] + INF ;
}
int main(){
    N = qr(),; rr int i ;//下面的++N是为了日后好操作qaq
    rep(i, 1, M) Add(qr(), qr(), qr()) ; ++ N, Ans = HLPP() ; cout << Ans << endl ; return 0 ; 
}

下面是个(fread)卡常版本(qaq)

#include <bits/stdc++.h>
#pragma GCC target("avx")
#pragma GCC optimize(3)
#pragma GCC optimize("Ofast")

#define MAXN 1202
#define min my_min
#define BG begin()
#define rr register
#define swap my_swap 
#define Iter iterator
#define INF 2147483647
#define rep(i, a, x) for(i = a ; i <= x ; ++ i)

using namespace std ;
typedef list<int> List ; int step;
const int ch_top=4e7+3;
char ch[ch_top],*now_r=ch-1,*now_w=ch-1;

inline int read(){
    while(*++now_r<'0');
    register int x=*now_r-'0';
    while(*++now_r>='0')x=x*10+*now_r-'0';
    return x;
}
inline void write(int x){
    static char st[20];static int top;
    while(st[++top]='0'+x%10,x/=10);
    while(*++now_w=st[top],--top);
    *++now_w='
';
}

struct Edge{
    int to, f, next ;
    Edge(int to,int f,int next):to(to),f(f),next(next){}
} ; vector <int> q, H, Extra, Set[MAXN], cnt ; List Gap[MAXN] ;
int Ans, N, M, S, T, max_H, now_H ; vector <Edge> E[MAXN] ; vector<List::iterator> Era_pos ; 

inline void eggs() { ;}//for free~
inline int my_min(int a, int b){return a & ((a - b) >> 31) | b & ( ~ (a - b) >> 31) ;}
inline void Add(int f, int v, int u){ E[u].push_back(Edge(v, f, E[v].size())), E[v].push_back(Edge(u, 0, E[u].size() - 1)) ; }

inline void Init_label(){
    rr int i, h = 0, t = 0 ;q.clear(), q.resize(N) ; 
    rr int qaq = N + 1 ; H.assign(qaq, qaq) ; H[T] = 0 ; q[t ++] = T ;
    while (h < t){
        rr int now = q[h], qwq = H[now] + 1 ; ++ h ; 
        for (vector <Edge> :: Iter k = E[now].begin() ; k != E[now].end() ; ++ k)
            if (H[k->to] == qaq && E[k->to][k->next].f) H[k->to] = qwq, ++ cnt[H[k->to]], q[t ++] = k->to ;
    }
    rep(i, 0, N - 1) Set[i].clear(), Gap[i].clear() ;
    rep(i, 0, N - 1) if (H[i] < N)  Era_pos[i] = Gap[H[i]].insert(Gap[H[i]].BG, i), (Extra[i] > 0) ? Set[H[i]].push_back(i) : eggs() ;
    max_H = now_H = H[q[-- t]] ; 
}
inline void Push(int x, Edge &e){
    rr int now_flow = min(Extra[x], e.f) ;
    Extra[x] -= now_flow, e.f -= now_flow, Extra[e.to] += now_flow, E[e.to][e.next].f += now_flow ;
    if (Extra[e.to] > 0 && Extra[e.to] <= now_flow) Set[H[e.to]].push_back(e.to) ;  // push it into heap
}
inline void _Push(int x){
    rr int i, x_h = N, t = H[x] ; 
    for (vector <Edge> :: Iter k = E[x].BG ; k != E[x].end() ; ++ k)
        if (k->f > 0)
            if (H[k->to] == H[x] - 1) { Push(x, *k) ; if (!Extra[x]) return ;} else x_h = min(x_h, H[k->to] + 1) ;
    if (cnt[H[x]] <= 1){
        for(i = t ; i <= max_H ; Gap[i].clear(), ++ i)
            for(List::Iter k = Gap[i].BG ; k != Gap[i].end() ; ++ k) cnt[H[*k]] --, H[*k] = N ; max_H = t - 1 ; return ;
    }
    -- cnt[t], Era_pos[x] = Gap[t].erase(Era_pos[x]) ; H[x] = x_h ; if (x_h == N) return ; 
    ++ cnt[x_h], Era_pos[x] = Gap[x_h].insert(Gap[x_h].begin(), x), max_H = max(now_H = x_h, max_H), Set[x_h].push_back(x) ;
}
int HLPP(){
    rr int i, now ; H.assign(N, 0) ; H[S] = N ; cnt.assign(N, 0) ; Era_pos.resize(N) ;
    rep(i, 0, N - 1) if (i != S) Era_pos[i] = Gap[H[i]].insert(Gap[H[i]].BG, i) ; cnt[0] = N - 1 ;
    Extra.assign(N, 0), Extra[S] = INF, Extra[T] = -INF ; rep(i, 0, E[S].size() - 1) Push(S, E[S][i]) ;  Init_label() ; 
    while (now_H >= 0) if (Set[now_H].empty()) now_H -- ; else now = Set[now_H].back(), Set[now_H].pop_back(), _Push(now) ;return Extra[T] + INF ;
}
int main(){
    fread(ch,1,ch_top,stdin); N = read(), M = read(), S = read(), T = read() ; rr int i ;
    rep(i, 1, M) Add(read(), read(), read()) ; ++ N, Ans = HLPP() ; write(Ans) ; fwrite(ch,1,now_w-ch,stdout) ;
}

撒fa~

(0x03~~)后记

• 这道题的经典版本好几个月之前我写了一天……然后今天又翻出来，发现了巨佬KevinYu抛了一个玉，我就打算优化一波……毒瘤啊，什么(vector)存边、(list)我都是第一次用呜呜……
• 不得不说…某些大佬的码风真是不可看啊…都写题解了怎么还这么…这么…（虽然自己的也不咋地qaq）
• 最后，人艰不拆，人艰不拆…

(oldsymbol{mathfrak{writter:Orchidany(pks)}})