带花树算法

带花树算法

先放上大神的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;
}