帶花樹算法

帶花樹算法

先放上大神的blog，我的認爲沒辦法比這位dalao解釋的更清楚。

帶花樹算法

在北京冬令營的時候，yby提到了「帶花樹開花」算法來解非二分圖的最大匹配。

因而，我打算看看這是個什麼玩意。其實以前，我已經對這個算法瞭解了個大概，可是。。。真的不敢去寫。

有一個叫Galil Zvi的人（應該叫計算機科學家），寫了篇論文：

Efficient Algorithms for Finding Maximal Matching in Graphs

（若是你在網上搜不到，能夠： http://builtinclz.abcz8.com/art/2012/Galil%20Zvi.pdf）

這篇論文真神啊，它解決了4個問題：

（通常圖+二分圖）的（最大匹配+最大權匹配）問題。

算法的思想、故事，請本身看論文吧。

這個論文告訴了咱們不少有趣的東西，例如：

 用Dinic實現的二分圖匹配的時間複雜度實際上是O(M*N^0.5)，這也許可以解釋爲何通常網絡流算法比Hungry要快了。

另外，帶花樹算法的正確性的證實比較困難；而其時間複雜度是能夠作到O(M*N^0.5)的，不過要詳細實現，那麼就快能到「ACM最長論文獎」了。

 

我寫了一個實例代碼：

http://builtinclz.abcz8.com/art/2012/ural1099.cpp

沒錯，這是用來解決URAL 1099 Work Schedule那題的。時間複雜度是O（N＾３）

簡述一下「帶花樹」算法吧：

它的核心思想仍是找增廣路。假設已經匹配好了一堆點，咱們從一個沒有匹配的節點s開始，使用BFS生成搜索樹。每當發現一個節點u，若是u尚未被匹配，那麼就能夠進行一次成功的增廣；不然，咱們就把節點u和它的配偶v一同接到樹上，以後把ｖ丟進隊列繼續搜索。咱們給每一個在搜索樹上的點一個類型：S或者T。當ｕ把它的配偶ｖ扔進隊列的時候，咱們把ｕ標記爲T型，ｖ標記爲S型。因而，搜索樹的樣子是這樣的：

　　　　　　　ｓ

　　　　　　 ／　　 ＼

　　　　　 ａ　　　　 ｂ

　　　　　｜　　　　｜

　　　　　ｃ　　　　ｄ

　　　　 ／　 ＼　　 ／　 ＼

　　　　 ｅ　 ｆ　　 ｕ　ｊ

　　　　｜　｜　　｜　｜

　　　　ｉ　ｊ　　ｖ　ｋ

其中，黑色豎線相連的兩個點是已經匹配好的，藍色斜線表示兩個點之間有邊，可是沒有配對。T型的用紅色，S型的用黑色。

 

這裏有個小問題：一個S型點ｄ在某一步擴展的時候發現了點ｕ，若是ｕ已經在搜索樹上了（即，出現了環），怎麼辦？

咱們規定，若是ｕ的類型是T型，就無視此次發現；（這意味着咱們找到了一個長度爲偶數的環，直接無視）

　　　　　　　ｓ

　　　　　　 ／　　 ＼

　　　　　 ａ　　　　 ｂ

　　　　　｜　　　　｜

　　　　　ｃ　　　　ｄ　　　若是連出來的邊是指向T型點的，就無視這個邊。

　　　　 ／　 ＼　　 ／　 ＼

　　　　 ｅ　 ｆ ＜ －　　 ｇ

　　　　｜　｜　　　　｜

　　　　ｉ　ｊ　　　　ｋ

不然，咱們找到了一個長度爲奇數的環，就要進行一次「縮花」的操做！所謂縮花操做，就是把這個環縮成一個點。

　　　　　　　ｓ

　　　　　　 ／　　 ＼

　　　　　 ａ　　　　 ｂ

　　　　　｜　　　　｜

　　　　　ｃ　　　　ｄ

　　　　 ／　 ＼　　 ／　 ＼

　　　　 ｅ　 ｆ　　 ｜　 ｇ

　　　　｜　｜　　 ｜　｜

　　　　ｉ　ｕ ＜－＋　ｋ

這個圖縮花以後變成了５個點（一個大點，或者叫一朵花，加原來的４個點）：

縮點完成以後，還要把原來環裏面的T型點通通變成S型點，以後扔到隊列裏去。

　　＋－－－－－－－－－－－－－＋

　　｜　　　　　　　　　　　　　｜

　　｜　　　　　ｓ　　　　　　　｜

　　｜　　　　 ／　　 ＼　　　　　｜

　　｜　　　 ａ　　　　 ｂ 　　　　｜

　　｜　　　｜　　　　｜　　　　｜　　　如今是一個點了！仍是一個S點。

　　｜　　　ｃ　　　　ｄ　　　　｜

　　｜　　 ／　 ＼　　 ／　＼　　　｜

＋－｜－－ 　　 ｆ－－ｕ 　　－－－｜－－－＋

｜　｜　　　　　　　　 　　　　　｜　　　｜

｜　｜　　　　　　　　　　　　　｜　　　｜

｜　｜　　　　　　　　　　　　　｜　　　｜

｜　＋－－－－－－－－－－－－－＋　　　｜

｜　　　　　　　　　　　　　　　　　　　｜

ｅ　　　　　　　　　　　　　　　　　　　ｇ

｜　　　　　　　　　　　　　　　　　　　｜

ｉ　　　　　　　　　　　　　　　　　　　ｋ

爲何能縮成一個點呢？咱們看一個長度爲奇數的環（例如上圖中的ｓ－ｂ－ｄ－ｕ－ｆ－ｃ－ａ－），若是咱們可以給它中的任意一個點找一個出度（配偶），那麼環中的其餘點正好能夠配成對，這說明，每一個點的出度都是等效的。例如，假設咱們可以給圖中的點ｄ另找一個配偶（例如ｄ＇好了），那麼，環中的剩下６個點正好能配成３對，一個很少，一個很多（算上ｄ和ｄ＇一共４對剛恰好）。

ａ－ｓ－ｂ－ｄ －ｄ＇　　　　　　　　　ａ　ｓ－ｂ　ｄ－ｄ＇

　＼　　　　｜　　　　　　　＝＞　　　　＼　　　　　

　　ｃ －ｆ－ｕ　　　　　　　　　　　　　　ｃ　ｆ－ｕ

這就是咱們縮點的思想來源。有一個勞苦功高的計算機科學家證實了：縮點以前和縮點以後的圖是否有增廣路的狀況是相同的。

縮起來的點又叫一朵花（ｂｌｏｓｓｏｍ）．

注意到，組成一朵花的裏面可能嵌套着更小的花。

 

當咱們最終找到一條增廣路的時候，要把嵌套着的花層層展開，還原出原圖上的增廣路出來。

 

嗯，如今你對實現這個算法有任何想法嗎？

天啊，還要縮點……寫死誰。。。。。。

我一開始也是這麼想的。

我看了一眼網上某個大牛的程序，以後結合本身的想法，很努力地寫出了一個能AC的版本。

實現的要點有什麼呢？

首先，咱們不「顯式」地表示花。咱們記錄一個Nｅｘｔ數組，表示最終增廣的時候的路徑上的後繼。同時，咱們維護一個並查集，表示每一個點如今在以哪一個點爲根的花裏（一個花被縮進另外一朵花以後就不算花了）。還要記錄每一個點的標記。

主程序是一段BFS。對於每一個由ｘ發展出來的點ｙ，分４種狀況討論：

１。ｘｙ是配偶（不說夫妻，這是非二分圖。。。）或者ｘｙ如今是一個點（在一朵花裏）：直接無視

２。ｙ是T型點：直接無視

３。ｙ目前單身：太好了，進行增廣！

４。ｙ是一個S型點：縮點！縮點！

縮點的時候要進行的工做：

１。找ｘ和ｙ的LCA（的根）ｐ。找LCA能夠用各類方法。。。直接樸素也行。

２。在Nｅｘｔ數組中把ｘ和ｙ接起來（表示它們造成環了！）

３。從ｘ、ｙ分別走到ｐ，修改並查集使得它們都變成一家人，同時沿路把Nｅｘｔ數組接起來。

 

Nｅｘｔ數組很奇妙。每時每刻，它實際造成了若干個掛在一塊兒的雙向鏈表來表示一朵花內部的走法。

　　　　　－－－－

　　　　／　　　　＼ ＜－－＋

　　　　｜　　　　｜　　　｜

　　　　｜　　　　｜ ＜－－＋

　　　　 ｖ　　　　 ｖ

　　　－－－－－－－－－－

　　／　　　　　　　　　　＼

＋－　　　　　　　　　　　　－－＋

｜　　　　　　　　　　　　　　　｜

｜　　　　　　　　　　　　　　　｜

＋－－－－ ＞ｓ　　 ＜－－－－－－＋　　　　　

 

 

 

有權圖的最大匹配怎麼作？

看論文吧。。。用相似KM的方法，不過，是給每一個花再來一個權值。真的很複雜。。。

有一我的寫了代碼，好像是GPL許可證。。。你最好想辦法搜到它的網站來看看版權的問題；總之，我先貼出來：

http://builtinclz.abcz8.com/art/2012/mwmatching.py、

 

筆者的補充

一開始總以爲這東西麼。。。沒有什麼用。。。由於到了考場上確定有不少特殊的性質（好比沒有奇環），或者數據規模很小（暴力枚舉全部點貪心）。這樣能夠用更少的時間獲得一個比較讓人滿意的分數。可是就有這麼倆集訓隊爺非得在論文裏討論通常圖的最大匹配。。。這就使省選出這類題目的可能性無限變大了。考慮到今年出了圓方樹，因此仍是去學習了一下。。。

筆者認爲理解了上面dalao的講解以後，對於帶花樹算法的實質內容以及實現形式都應該有了一個比較深的認識。可是代碼確實不是很好打以致於筆者也是各類蒐羅才寫出了一份看上去比較優美的。而後加上了一些筆者本身的看法，也就是註釋，但願可以幫到各位理解代碼吧。。

#include <iostream>
#include <cstdlib>
#include <cmath>
#include <string>
#include <cstring>
#include <cstdio>
#include <algorithm>
#include <queue>
#include <set>
#include <map>
#define re register
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
#define MAXN 5001
#define MAXM 1000001
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
#define ms(arr) memset(arr, 0, sizeof(arr))
const int inf = 0x3f3f3f3f;
int f[MAXN],head[MAXN],next[MAXN],num,match[MAXN],t,vis[MAXN],mark[MAXN],n;
//f是並查集，match是匹配數組，mark表示S/T點。
int x[MAXN],y[MAXN],L;
struct po
{
    int from,nxt,to;
}edge[MAXM];
queue<int> q;
inline int read()
{
    int x=0,c=1;
    char ch=' ';
    while((ch>'9'||ch<'0')&&ch!='-')ch=getchar();
    while(ch=='-') c*=-1,ch=getchar();
    while(ch<='9'&&ch>='0')x=x*10+ch-'0',ch=getchar();
    return x*c;
}
int find(int x)
{
    return f[x]==x?f[x]:f[x]=find(f[x]);
}
inline void add_edge(int from,int to)
{
    edge[++num].nxt=head[from];
    edge[num].to=to;
    head[from]=num;
}
inline void add(int from,int to)
{
    add_edge(from,to);
    add_edge(to,from);
}
inline int lca(int x,int y)//樸素找lca，雖然筆者本身是死記硬背的。。貌似只會樹剖找lca。。
{
    t++;
    while(1){
        if(x){
            x=find(x);//點要對應到相應的花上。
            if(vis[x]==t) return x;
            vis[x]=t;
            if(match[x]) x=next[match[x]];
            else x=0;
        }
        swap(x,y);
    }
}
inline void flower(int a,int p)
{
    while(a!=p){
        int b=match[a],c=next[b];
        if(find(c)!=p) next[c]=b;
        //next數組是用來標記花中的路徑，結合match就實現了雙向鏈表的功能。
        //咱們能夠知道一個奇環裏的全部點只要向外面連邊就有可能與外面匹配，因此全部奇環中的點都         
　　　　　　要變成S型點扔到隊列中。又由於環內部確定是匹配飽和的，因此這一堆點最多隻能匹配出來一             
　　　　　　個點，因此每次匹配到一個點就能夠跳出循環了。
        if(mark[b]==2) q.push(b),mark[b]=1;
        if(mark[c]==2) q.push(c),mark[c]=1;
        f[find(a)]=find(b);f[find(b)]=find(c);
        a=c;
    }
}
inline void work(int s)
{
    for(re int i=1;i<=n;i++)
        next[i]=mark[i]=vis[i]=0,f[i]=i;//每一個階段都須要從新標記
    while(!q.empty()) q.pop();
    mark[s]=1;q.push(s);
    while(!q.empty()){
        if(match[s]) return;
        int x=q.front();q.pop();//隊列中的點均爲S型
        for(re int i=head[x];i;i=edge[i].nxt){
            int y=edge[i].to;
            if(match[x]==y) continue;//x與y是配偶，忽略掉
            if(find(x)==find(y)) continue;//x與y在一朵花裏，忽略掉
            if(mark[y]==2) continue;//y是T形點，忽略掉。
            if(mark[y]==1) {//y是s型點，須要縮點
                int r=lca(x,y);
                if(find(x)!=r) next[x]=y;//x和r不在同一個花朵， next標記花朵內的路徑
                if(find(y)!=r) next[y]=x;//同上
                flower(x,r);flower(y,r);//將r--x--y--r縮成一個點
            }
            else if(!match[y]){//y沒有匹配過，能夠增廣
                next[y]=x;
                for(re int u=y;u;){//增廣操做
                    int v=next[u],w=match[v];
                    match[v]=u;match[u]=v;u=w;
                }
                break;
            } else {//y點已經匹配過，將其配偶做爲S型點拉進來，y點自己爲T型點，繼續搜索。
                next[y]=x;
                mark[match[y]]=1;
                q.push(match[y]);
                mark[y]=2;
            }
        }
    }
}
inline void out(){
    for(re int i=1;i<=n;i++)
            if(!match[i]){
                printf("NO\n");
                return;
            }
        printf("YES\n");
}
int main() 
{
    //freopen("date.in","r",stdin);
    while(cin>>n){
        memset(match,0,sizeof(match));
        memset(head,0,sizeof(head));num=0;t=0;
        for(re int i=1;i<=n;i++)
            x[i]=read(),y[i]=read();
        L=read();
        for(re int i=1;i<=n;i++)
            for(re int j=i+1;j<=n;j++){
                if(abs(x[i]-x[j])+abs(y[i]-y[j])<=L){
                    add(i,j);
                }
            }
        for(re int i=1;i<=n;i++)
            if(!match[i]) work(i);
        out();
    }
    return 0;
}

 

 而後還有一個就是帶權帶花樹的代碼以及實現。。筆者找到了一份比價簡單明白的代碼，雖然本身尚未研究，不過給你們放上來，但願對你們有所幫助。

#include<bits/stdc++.h>
using namespace std;
//from vfleaking
//本身進行一些進行一些小修改 
#define INF INT_MAX
#define MAXN 400
struct edge{
    int u,v,w;
    edge(){}
    edge(int u,int v,int w):u(u),v(v),w(w){}
};
int n,n_x;
edge g[MAXN*2+1][MAXN*2+1];
int lab[MAXN*2+1];
int match[MAXN*2+1],slack[MAXN*2+1],st[MAXN*2+1],pa[MAXN*2+1];
int flower_from[MAXN*2+1][MAXN+1],S[MAXN*2+1],vis[MAXN*2+1];
vector<int> flower[MAXN*2+1];
queue<int> q;
inline int e_delta(const edge &e){ // does not work inside blossoms
    return lab[e.u]+lab[e.v]-g[e.u][e.v].w*2;
}
inline void update_slack(int u,int x){
    if(!slack[x]||e_delta(g[u][x])<e_delta(g[slack[x]][x]))slack[x]=u;
}
inline void set_slack(int x){
    slack[x]=0;
    for(int u=1;u<=n;++u)
        if(g[u][x].w>0&&st[u]!=x&&S[st[u]]==0)update_slack(u,x);
}
void q_push(int x){
    if(x<=n)q.push(x);
    else for(size_t i=0;i<flower[x].size();i++)q_push(flower[x][i]);
}
inline void set_st(int x,int b){
    st[x]=b;
    if(x>n)for(size_t i=0;i<flower[x].size();++i)
            set_st(flower[x][i],b);
}
inline int get_pr(int b,int xr){
    int pr=find(flower[b].begin(),flower[b].end(),xr)-flower[b].begin();
    if(pr%2==1){//檢查他在前一層圖是奇點還是偶點
        reverse(flower[b].begin()+1,flower[b].end());
        return (int)flower[b].size()-pr;
    }else return pr;
}
inline void set_match(int u,int v){
    match[u]=g[u][v].v;
    if(u>n){
        edge e=g[u][v];
        int xr=flower_from[u][e.u],pr=get_pr(u,xr);
        for(int i=0;i<pr;++i)set_match(flower[u][i],flower[u][i^1]);
        set_match(xr,v);
        rotate(flower[u].begin(),flower[u].begin()+pr,flower[u].end());
    }
}
inline void augment(int u,int v){
    for(;;){
        int xnv=st[match[u]];
        set_match(u,v);
        if(!xnv)return;
        set_match(xnv,st[pa[xnv]]);
        u=st[pa[xnv]],v=xnv;
    }
}
inline int get_lca(int u,int v){
    static int t=0;
    for(++t;u||v;swap(u,v)){
        if(u==0)continue;
        if(vis[u]==t)return u;
        vis[u]=t;//這種方法能夠不用清空v陣列 
        u=st[match[u]];
        if(u)u=st[pa[u]];
    }
    return 0;
}
inline void add_blossom(int u,int lca,int v){
    int b=n+1;
    while(b<=n_x&&st[b])++b;
    if(b>n_x)++n_x;
    lab[b]=0,S[b]=0;
    match[b]=match[lca];
    flower[b].clear();
    flower[b].push_back(lca);
    for(int x=u,y;x!=lca;x=st[pa[y]])
        flower[b].push_back(x),flower[b].push_back(y=st[match[x]]),q_push(y);
    reverse(flower[b].begin()+1,flower[b].end());
    for(int x=v,y;x!=lca;x=st[pa[y]])
        flower[b].push_back(x),flower[b].push_back(y=st[match[x]]),q_push(y);
    set_st(b,b);
    for(int x=1;x<=n_x;++x)g[b][x].w=g[x][b].w=0;
    for(int x=1;x<=n;++x)flower_from[b][x]=0;
    for(size_t i=0;i<flower[b].size();++i){
        int xs=flower[b][i];
        for(int x=1;x<=n_x;++x)
            if(g[b][x].w==0||e_delta(g[xs][x])<e_delta(g[b][x]))
                g[b][x]=g[xs][x],g[x][b]=g[x][xs];
        for(int x=1;x<=n;++x)
            if(flower_from[xs][x])flower_from[b][x]=xs;
    }
    set_slack(b);
}
inline void expand_blossom(int b){ // S[b] == 1
    for(size_t i=0;i<flower[b].size();++i)
        set_st(flower[b][i],flower[b][i]);
    int xr=flower_from[b][g[b][pa[b]].u],pr=get_pr(b,xr);
    for(int i=0;i<pr;i+=2){
        int xs=flower[b][i],xns=flower[b][i+1];
        pa[xs]=g[xns][xs].u;
        S[xs]=1,S[xns]=0;
        slack[xs]=0,set_slack(xns);
        q_push(xns);
    }
    S[xr]=1,pa[xr]=pa[b];
    for(size_t i=pr+1;i<flower[b].size();++i){
        int xs=flower[b][i];
        S[xs]=-1,set_slack(xs);
    }
    st[b]=0;
}
inline bool on_found_edge(const edge &e){
    int u=st[e.u],v=st[e.v];
    if(S[v]==-1){
        pa[v]=e.u,S[v]=1;
        int nu=st[match[v]];
        slack[v]=slack[nu]=0;
        S[nu]=0,q_push(nu);
    }else if(S[v]==0){
        int lca=get_lca(u,v);
        if(!lca)return augment(u,v),augment(v,u),true;
        else add_blossom(u,lca,v);
    }
    return false;
}
inline bool matching(){
    memset(S+1,-1,sizeof(int)*n_x);
    memset(slack+1,0,sizeof(int)*n_x);
    q=queue<int>();
    for(int x=1;x<=n_x;++x)
        if(st[x]==x&&!match[x])pa[x]=0,S[x]=0,q_push(x);
    if(q.empty())return false;
    for(;;){
        while(q.size()){
            int u=q.front();q.pop();
            if(S[st[u]]==1)continue;
            for(int v=1;v<=n;++v)
                if(g[u][v].w>0&&st[u]!=st[v]){
                    if(e_delta(g[u][v])==0){
                        if(on_found_edge(g[u][v]))return true;
                    }else update_slack(u,st[v]);
                }
        }
        int d=INF;
        for(int b=n+1;b<=n_x;++b)
            if(st[b]==b&&S[b]==1)d=min(d,lab[b]/2);
        for(int x=1;x<=n_x;++x)
            if(st[x]==x&&slack[x]){
                if(S[x]==-1)d=min(d,e_delta(g[slack[x]][x]));
                else if(S[x]==0)d=min(d,e_delta(g[slack[x]][x])/2);
            }
        for(int u=1;u<=n;++u){
            if(S[st[u]]==0){
                if(lab[u]<=d)return 0;
                lab[u]-=d;
            }else if(S[st[u]]==1)lab[u]+=d;
        }
        for(int b=n+1;b<=n_x;++b)
            if(st[b]==b){
                if(S[st[b]]==0)lab[b]+=d*2;
                else if(S[st[b]]==1)lab[b]-=d*2;
            }
        q=queue<int>();
        for(int x=1;x<=n_x;++x)
            if(st[x]==x&&slack[x]&&st[slack[x]]!=x&&e_delta(g[slack[x]][x])==0)
                if(on_found_edge(g[slack[x]][x]))return true;
        for(int b=n+1;b<=n_x;++b)
            if(st[b]==b&&S[b]==1&&lab[b]==0)expand_blossom(b);
    }
    return false;
}
inline pair<long long,int> weight_blossom(){
    memset(match+1,0,sizeof(int)*n);
    n_x=n;
    int n_matches=0;
    long long tot_weight=0;
    for(int u=0;u<=n;++u)st[u]=u,flower[u].clear();
    int w_max=0;
    for(int u=1;u<=n;++u)
        for(int v=1;v<=n;++v){
            flower_from[u][v]=(u==v?u:0);
            w_max=max(w_max,g[u][v].w);
        }
    for(int u=1;u<=n;++u)lab[u]=w_max;
    while(matching())++n_matches;
    for(int u=1;u<=n;++u)
        if(match[u]&&match[u]<u)
            tot_weight+=g[u][match[u]].w;
    return make_pair(tot_weight,n_matches);
}
inline void init_weight_graph(){
    for(int u=1;u<=n;++u)
        for(int v=1;v<=n;++v)
            g[u][v]=edge(u,v,0);
}
int main(){
    int m;
    scanf("%d%d",&n,&m);
    init_weight_graph();
    for(int i=0;i<m;++i){
        int u,v,w;
        scanf("%d%d%d",&u,&v,&w);
        g[u][v].w=g[v][u].w=w;
    }
    printf("%lld\n",weight_blossom().first);
    for(int u=1;u<=n;++u)printf("%d ",match[u]);puts("");
    return 0;
}