[转]带花树,Edmonds's matching algorithm,一般图最大匹配

　　看了两篇博客，觉得写得不错，便收藏之。。

　　首先是第一篇，转自某Final牛

带花树……其实这个算法很容易理解，但是实现起来非常奇葩（至少对我而言）。

除了wiki和amber的程序我找到的资料看着都不大靠谱

比如昨晚找到一篇鄙视带花树的论文，然后介绍了一种O(E)的一般图最大匹配……我以为找到了神论文，然后ACM_DIY众神纷纷表示这个是错的……于是神论文成为了”神论文“……

又比如围观nocow上带花树标程，一看……这显然是裸的匈牙利算法……货不对板啊

当然……如果二分图的匈牙利算法还不会请先围观求二分图最大匹配的匈牙利算法。

实际上任意图求最大匹配也是找增广路，但是由于奇环的出现，找增广路变得困难。

首先明确一点，增广路上是不能有重复出现的点的。

‍‍
‍

二分图中，匹配边可以看作是有向的，比如定义总是从X集指向Y集。假若定义了起点必须在X集中，那么增广路中出现该匹配边时，必然是按照这个方向的。所以一个点在增广路中的奇偶性是确定的。

而这个图中，从增广路3->1->4->5和2->4->1->6可以看出，对于有奇环的任意图，1和4这两个点在增广路中所在位置的奇偶性不再一定。于是我们考虑处理这些奇环。

定义奇环：包含2k+1个点和k条匹配边的一个环。（如果不是这样，我们找增广路不会走上去）

对于这个奇环，k条匹配覆盖了2k个点，那么显然有一个点未被覆盖。我们拿出这个点来讨论。

比如图中的1号点就是这个这个特殊的点。除了这个点以外，其它的点都被覆盖了，所以只能向外连非匹配边，而1号点可以向外连匹配边
或非匹配边。

如果1号点没有被外面的点匹配，那么无论从其它的哪个点走进来，都能以1为终点找到增广路。(要么顺时针跑到1，要么逆时针)

同理如果1号点被外面的点匹配了，那么无论从其它的哪个点走进来，都能把这个圈看成一个点，然后从1的那条匹配边穿出去。(要么顺时针，要么逆时针)

于是这个奇环就可以看成一个点，其主要特性由1号点体现(诸如和谁匹配了之流)。

这个合成点就叫做花。这个算法的思想就是不断地把奇环合成点，直至找到增广路（合成了某朵花以后就把整朵花当成一个点）。

考虑用BFS搜索增广路。

围观wiki这个图

‍

由于BFS的性质，我们找到奇环只能是和同层的点，或者下下一层的点。

然后奇环的关键点必然是这棵BFS树里深度最浅的点。然后考虑合成以后，花如何展开对应的路径，使得我们能够增广。

花套花这个东西想起来都纠结>_<。

amber的程序里面并没有把点真的合成，只是弄了一个表示集合的标号：Base，然后邻接矩阵就不用变来变去了。

对于花中连向父亲的是匹配边的点，他的增广路显然是直接顺着父亲走，而如果连向父亲的边是非匹配边的点，那么显然是往后走然后跑过红色的横插边，然后再向上跑回关键点。

注意到如果连向父子的边是匹配边的点原先是不需要Father这个域来描述的，直接用表示匹配的那个域就可以了。但是现在在花中，他的Father这个域就要起作用了，用来向后指向，然后绕过红色横插边然后再跑回关键点。

实在是太精妙了。

  1 //Problem:http://acm.timus.ru/problem.aspx?space=1&num=1099
  2 #include <cstdio>
  3 #include <cstdlib>
  4 #include <cstring>
  5 #include <iostream>
  6 #include <algorithm>
  7 using namespace std;
  8 const int N=250;
  9 int n;
 10 int head;
 11 int tail;
 12 int Start;
 13 int Finish;
 14 int link[N];     //表示哪个点匹配了哪个点
 15 int Father[N];   //这个就是增广路的Father……但是用起来太精髓了
 16 int Base[N];     //该点属于哪朵花
 17 int Q[N];
 18 bool mark[N];
 19 bool map[N][N];
 20 bool InBlossom[N];
 21 bool in_Queue[N];
 22  
 23 void CreateGraph(){
 24     int x,y;
 25     scanf("%d",&n);
 26     while (scanf("%d%d",&x,&y)!=EOF)
 27       map[x][y]=map[y][x]=1;
 28 }
 29  
 30 void BlossomContract(int x,int y){
 31     fill(mark,mark+n+1,false);
 32     fill(InBlossom,InBlossom+n+1,false);
 33     #define pre Father[link[i]]
 34     int lca,i;
 35     for (i=x;i;i=pre) {i=Base[i]; mark[i]=true; }
 36     for (i=y;i;i=pre) {i=Base[i]; if (mark[i]) {lca=i; break;} }  //寻找lca之旅……一定要注意i=Base[i]
 37     for (i=x;Base[i]!=lca;i=pre){
 38         if (Base[pre]!=lca) Father[pre]=link[i]; //对于BFS树中的父边是匹配边的点，Father向后跳
 39         InBlossom[Base[i]]=true;
 40         InBlossom[Base[link[i]]]=true;
 41     }
 42     for (i=y;Base[i]!=lca;i=pre){
 43         if (Base[pre]!=lca) Father[pre]=link[i]; //同理
 44         InBlossom[Base[i]]=true;
 45         InBlossom[Base[link[i]]]=true;
 46     }
 47     #undef pre
 48     if (Base[x]!=lca) Father[x]=y;     //注意不能从lca这个奇环的关键点跳回来
 49     if (Base[y]!=lca) Father[y]=x;
 50     for (i=1;i<=n;i++)
 51       if (InBlossom[Base[i]]){
 52           Base[i]=lca;
 53           if (!in_Queue[i]){
 54               Q[++tail]=i;
 55               in_Queue[i]=true;     //要注意如果本来连向BFS树中父结点的边是非匹配边的点，可能是没有入队的
 56           }
 57       }
 58 }
 59  
 60 void Change(){
 61     int x,y,z;
 62     z=Finish;
 63     while (z){
 64         y=Father[z];
 65         x=link[y];
 66         link[y]=z;
 67         link[z]=y;
 68         z=x;
 69     }
 70 }
 71  
 72 void FindAugmentPath(){
 73     fill(Father,Father+n+1,0);
 74     fill(in_Queue,in_Queue+n+1,false);
 75     for (int i=1;i<=n;i++) Base[i]=i;
 76     head=0; tail=1;
 77     Q[1]=Start;
 78     in_Queue[Start]=1;
 79     while (head!=tail){
 80         int x=Q[++head];
 81         for (int y=1;y<=n;y++)
 82           if (map[x][y] && Base[x]!=Base[y] && link[x]!=y)   //无意义的边
 83             if ( Start==y || link[y] && Father[link[y]] )    //精髓地用Father表示该点是否
 84                 BlossomContract(x,y);
 85             else if (!Father[y]){
 86                 Father[y]=x;
 87                 if (link[y]){
 88                     Q[++tail]=link[y];
 89                     in_Queue[link[y]]=true;
 90                 }
 91                 else{
 92                     Finish=y;
 93                     Change();
 94                     return;
 95                 }
 96             }
 97     }
 98 }
 99  
100 void Edmonds(){
101     memset(link,0,sizeof(link));
102     for (Start=1;Start<=n;Start++)
103       if (link[Start]==0)
104         FindAugmentPath();
105 }
106  
107 void output(){
108     fill(mark,mark+n+1,false);
109     int cnt=0;
110     for (int i=1;i<=n;i++)
111       if (link[i]) cnt++;
112     printf("%d
",cnt);
113     for (int i=1;i<=n;i++)
114       if (!mark[i] && link[i]){
115           mark[i]=true;
116           mark[link[i]]=true;
117           printf("%d %d
",i,link[i]);
118       }
119 }
120  
121 int main(){
122 //    freopen("input.txt","r",stdin);
123     CreateGraph();
124     Edmonds();
125     output();
126     return 0;
127 }

　　然后还有一篇，链接请猛戳。。

在北京冬令营的时候，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