1. 程式人生 > >網絡最大流·預流推進及其衍生

網絡最大流·預流推進及其衍生

assign define 節點 允許 其他 地方 class 一個 def

原文地址:http://www.orchidany.cf/2019/01/11/HLPP/#more

#define \(u\)的伴點集合 與\(u\)相隔一條邊的且\(u\)能達到的點的集合

\(0x00~ {}~Preface\)

\(HLPP(Highest~Label~Preflow~Push)\)最高標簽預流推進算法是處理網絡最大流裏兩種常用方法——增廣路&預流推進中,預流推進算法的一種。據傳由\(tarjan\)發明怎麽又是他 ,並被其他科學家證明了其復雜度是緊卻的\(O(n^2\sqrt m)\)。在隨機數據中不遜色於普通的增廣路算法,而在精心構造的數據中無法被卡,所以是一種可以替代\(Dinic\)

的方法(隨我怎麽說,代碼又長又難調,所以還是\(Dinic\)好啊\(\rm{TAT}\)

但無論怎樣,\(wiki\)裏面已經承認\(HLPP\)是現在最優秀的網絡流算法了。

那麽預流推進這個大門類裏面,思想都差不多。大抵上就是我們對每個點記錄超額流(\(Extra~Flow\)) ,即允許流在非源點暫時存儲,並伺機將超額流推送出去。不可推送的,就會流回源點。那麽最終答案顯然存儲在\(Extra[T]\)裏面。

但同時這也有一個問題,就是會出現兩個點相互推送不停的情況。為了防止這樣,我們采用最高標號的策略,給每個點一個高度,對於一個點\(u?\)以及它的伴點集合\(\{v\}?\),當且僅當\(h_u = h_v + 1?\)

時才可以推送流。並且我們對於源點\(S?\),設置\(h_S = N?\),並對於\(S?\)實行無限制推送。那麽最後的答案就保存在\(Extra[T]?\)裏面 。

但有時,我們發現有個點是”谷“,即周圍點的高度都比它高,但是它有超額流。那麽我們此時考慮拔高它的高度,即重貼標簽(\(relabel?\))操作。

\(0x01\quad\) 初步的算法流程

以下我們用\(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\)左右的好成績,但事實上,這個速度依舊慢的要死。

註意!這個寫法是經典寫法,其時間復雜度是緊卻的\(\boldsymbol{\rm{\Theta(n^2mlogn)}}\)的,也就是說在\(\boldsymbol{n=1200}\)時會帶一個\(\boldsymbol{10}\)倍的常數

怎麽優化呢——

\(\boldsymbol{0x02~~Optimization}\)

各位,你們將會見到迄今為止OI界最喪心病狂的優化(之一)……

來,我們首先思考思考普通的HLPP到底會慢在哪裏:

  • \(STL\)支持的\(heap\)(比如優先隊列)實在是太太太…太慢了!

  • 每次\(Gap\)優化,我們的時間復雜度是緊確\(\Theta(n)\)。的這顯然不合算,因為假設我當前的\(\boldsymbol{gap}\)(斷層)正好位於倒數第一高的點和倒數第二高的點,那麽也就相當於我單次會浪費\(\boldsymbol{\Theta(n)}\)的時間

事實上…普通的\(HLPP\)代碼並不長,主要問題就是這兩個。

我們考慮,如果不用堆的話怎麽做呢?

呃…不用堆的意思並不是我們不從高度最大的點開始推送。這個地方需要一個\(idea\)——在\(HLPP\)中,高度函數\(\boldsymbol{H(x)}\)和點數集大小\(\boldsymbol{N(x)}\)是廣義同階的。 換句話說,我們可以考慮從高度入手。

換句話說,我們原來是通過節點編號訪問節點以及其高度,現在我們如果從高度入手,再去訪問節點,我們就可以做到\(\boldsymbol{O(n)}\)而不是\(\boldsymbol{\rm{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='\n';
}

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)
  • 最後,人艱不拆,人艱不拆…

\(\boldsymbol{\mathfrak{writter:Orchidany(pks)}}\)

網絡最大流·預流推進及其衍生