hdu 2255 奔小康赚大钱 (二分图最佳匹配,KM算法模板题)
题意:传说在遥远的地方有一个非常富裕的村落,有一天,村长决定进行制度改革:重新分配房子。 这可是一件大事,关系到人民的住房问题啊。村里共有n间房间,刚好有n家老百姓,考虑到每家都要有房住(如果有老百姓没房子住的话,容易引起不安定因素),每家必须分配到一间房子且只能得到一间房子。 另一方面,村长和另外的村领导希望得到最大的效益,这样村里的机构才会有钱.由于老百姓都比较富裕,他们都能对每一间房子在他们的经济范围内出一定的价格,比如有3间房子,一家老百姓可以对第一间出10万,对第2间出2万,对第3间出20万.(当然是在他们的经济范围内).现在这个问题就是村领导怎样分配房子才能使收入最大.(村民即使有钱购买一间房子但不一定能买到,要看村领导分配的).
思路:这是一道裸的二分图最佳匹配题目。
如果G为加权二分图,则权值和最大的完备匹配称为最佳匹配.
解决二分图最佳匹配的常用算法为KM算法,该算法的前置技能点为匈牙利算法。
参考了以下文章: 关于KM算法的详细解释 km算法 求最大权二分匹配的KM算法
先粗浅得说下我对这个算法的理解,等A掉足够多的题目之后再回来总结。
因为才刚刚学习这个算法,可能有很多没理解透彻或者理解错误的地方,还望各位多多指教。
KM算法引入了一个交顶标的概念,lx[i]表示X集合中第i个点的顶标,ly[i]表示Y集合中第i个点的顶标。至于这个顶标具体代表什么先不用考虑,只需要知道,这个概念的引入是为了解决二分图最佳匹配的问题。
如果一组xi,yi,满足lx[xi]+ly[yi]-w[xi][yi]>=0成立,那么我们称这组点为【可行的】
我们要时刻保证每组点都是可行的。
如果我们把所有满足lx[xi]+ly[yi]==w[xi][yi]的(xi,yi)导出,形成一个新的图,那么新图的最大匹配一定是原图的最佳匹配。原因也很简单:新图的权值和等于所有的顶标和,而根据【可行】点对的定义,最好的情况就是所有不等式都取了等号。
为了保证初始时lx[xi]+ly[yi]-w[xi][yi]对于任意xi,yi成立,不妨这样构造lx[],ly[]:
lx[i]=max(w[i][j]),(1=<j<=ny) ,ly[i]=0;
注意是“不妨”,也就是初始值并不是唯一的,只不过这种构造方法在满足不等式的前提下相对简单。
然后就和hungary算法基本类似,从每一个X集合中的点开始寻找增广路。
但是。。但是哪有那么好的事情。。。就那么恰好满足lx[xi]+ly[yi]=w[xi][yi]了。。?
没有满足的怎么办。
接下来就要说KM算法最关键的一步,调整顶标。
我们调整顶标的目的是让更多的点对满足式子lx[xi]+ly[yi]==w[xi][yi],从而进入相等子图,直到图中存在仅由可行边组成的完全匹配为止.
**调整的过程中仍然要保证所有点对都【可行】**也就是保证lx[xi]+ly[yi]>=w[xi][yi]一直成立。
由于当前lx[xi]+ly[yi]-w[xi][yi]是比我们预期的大,所以我们要想办法减小。
具体方法为:
方法为:将所有在增广轨中(就是在增广过程中遍历到)的X方点的标号全部减去一个常数d,所有在增广轨中的Y方点的标号全部加上一个常数d,则对于图中的任意一条边(i, j, W)(i为X方点,j为Y方点):
点对经过调整后分为四种情况:
<1>i和j都在增广轨中:此时边(i, j)的(lx[i]+ly[j])值不变,也就是这条边的可行性不变(原来是可行边则现在仍是,原来不是则现在仍不是); **<2>i在增广轨中而j不在:此时边(i, j)的(lx[i]+ly[j])的值减少了d,也就是原来这条边不是可行边(否则j就会被遍历到了),而现在可能是;** <3>j在增广轨中而i不在:此时边(i, j)的(lx[i]+ly[j])的值增加了d,也就是原来这条边不是可行边(若这条边是可行边,则在遍历到j时会紧接着执行DFS(i),此时i就会被遍历到),现在仍不是;<4>i和j都不在增广轨中:此时边(i, j)的(lx[i]+ly[j])值不变,也就是这条边的可行性不变。
简单得说:(1) (3)(4)三种情况,不会改变一组点的可行性,用白话解释就是,以前满足那个等式而进入相等子图的点对不会因为调整而从相等子图中出来,以前不满足那个等式的点不会因为调整而进入相等子图。
而对于情况(2),因为lx[xi]+ly[yi]-w[xi][yi]减少了,那么以前不满足,现在有可能满足。注意,现在还只是可能满足。
那么这个常数d到底是多少?怎么确定?
是的,接下来的着手点就在这个d上。
那么d的值应取多少?显然,整个点标不能失去可行性,也就是对于上述的第<2>类边,其lx[i]+ly[j]>=W这一性质不能被改变,故取所有第<2>类边的(lx[i]+ly[j]-W)的最小值作为d值即可。这样一方面可以保证点标的可行性,另一方面,**经过这一步后,图中至少会增加一条可行边。**
这样的最小值d,会恰好使得得到最小值d的那组点对的lx[i]+ly[j]-w的值从大于0变成0(好像是废话,我只是为了强调一下)
但是如果朴素得去找最小值更新d的话,复杂度是O(n4)(不会分析复杂度。。。mzdd)
1#define M 505
2#define inf 0x3fffffff
3bool sx[M], sy[M];
4int match[M], w[M][M], n, m, d, lx[M], ly[M];
5//n:左集元素个数; m:右集元素个数
6void init ()
7{
8 memset (w, 0, sizeof(w)); //不一定要,求最小值一般要初始化为负无穷!
9}
1bool dfs (int u)
2{
3 int v; sx[u] = true;
4 for (v = 0; v < m; v++)
5 {
6 if (!sy[v] && lx[u]+ly[v]==w[u][v])
7 {
8 sy[v] = true;
9 if (match[v] == -1 || dfs (match[v]))
10 {
11 match[v] = u;
12 return true;
13 }
14 }
15 }
16 return false;
17}
1int KM ()
2{
3 int i, j, k, sum = 0;
4 memset (ly, 0, sizeof(ly));
5 for (i = 0; i < n; i++)
6 {
7 lx[i] = -inf;
8 for (j = 0; j < m; j++)
9 if (lx[i] < w[i][j])
10 lx[i] = w[i][j];
11 }
12 memset (match, -1, sizeof(match));
13 for (i = 0; i < n; i++)
14 {
15 while (1)
16 {
17 memset (sx, false, sizeof(sx));
18 memset (sy, false, sizeof(sy));
19 if (dfs (i))
20 break;
21 d = inf;
22 for (j = 0; j < n; j++)
23 if (sx[j])
24 for (k = 0; k < m; k++)
25 if (!sy[k])
26 d = min (d, lx[j]+ly[k]-w[j][k]);
27 if (d == inf) //找不到完美匹配
28 return -1;
29 for (j = 0; j < n; j++)
30 if (sx[j])
31 lx[j] -= d;
32 for (j = 0; j < m; j++)
33 if (sy[j])
34 ly[j] += d;
35 }
36 }
37 for (i = 0; i < m; i++)
38 if (match[i] > -1)
39 sum += w[match[i]][i];
40 return sum;
41}
**O(n^3)的优化**:如果每次都花O(n^2)的时间时间去找 min(lx[i] + ly[j] - mp[i][j]);显然,总的时间复杂度是O(n^3),这里加一个slack[]数组,记录每次dfs找完美匹配时lx[i] + ly[j] - mp[i][j]的最小值,在实现多次调整时只需要在slack[]里找到调整值d就可以。
回过头来看KM算法的过程,像一种奇妙的构造法
根据给定的规则去初始化和调整顶标,就可以得到二分图的最佳匹配。
正面来思考比较难以想到这种巧妙的构造法。在理解清楚整个过程之前,对于顶标到底是什么意义,为什么要满足那个等式,都是一脸蒙逼。。
但是回过头来,就有了一种柳暗花明的感觉。
回过头来说这道题:裸的二分图最佳匹配,n<=300,存图可以用邻接矩阵。。
然后意识到km算法O(N3)的复杂度。。。似乎存图总是可以邻接矩阵?
说得差不多了,代码里还有一些细节的注释。
/* ***********************************************
Author :111qqz
Created Time :2016年06月01日 星期三 11时23分42秒
File Name :code/hdu/2255.cpp
************************************************ */
1#include <cstdio>
2#include <cstring>
3#include <iostream>
4#include <algorithm>
5#include <vector>
6#include <queue>
7#include <set>
8#include <map>
9#include <string>
10#include <cmath>
11#include <cstdlib>
12#include <ctime>
13#define fst first
14#define sec second
15#define lson l,m,rt<<1
16#define rson m+1,r,rt<<1|1
17#define ms(a,x) memset(a,x,sizeof(a))
18typedef long long LL;
19#define pi pair < int ,int >
20#define MP make_pair
1using namespace std;
2const double eps = 1E-8;
3const int dx4[4]={1,0,0,-1};
4const int dy4[4]={0,-1,1,0};
5const int inf = 0x3f3f3f3f;
6const int N=305;
7int w[N][N];
8int n;
9int link[N];
10int lx[N],ly[N];//顶标
11int slk[N]; //slk为一个优化函数,表示在DFS增广过程中,所有搜到的与该Y方点关联的边的(lx+ly-W)的最小值
12bool visx[N],visy[N];
13bool find( int u)
14{
15 visx[u] = true;
16 for ( int v = 1 ; v <= n ; v++) //这道题特殊:每两个点都有边相连。
17 { //对于有些题目中可能两个集合某些点没有边相连,我们可以看作有权值为0的边。
18 if (visy[v]) continue;
19 int tmp = lx[u] + ly[v] - w[u][v]; //要时刻保证tmp>=0,并称一样一组点标为可行的。
20// cout<<"tmp:"<<tmp<<endl;
21 if (tmp==0)
22 {
23 visy[v] = true;
24 if (link[v]==-1||find(link[v]))
25 {
26 // cout<<"link[v]:"<<link[v]<<endl;
27 link[v] = u;
28 return true;
29 }
30 }
31 else
32 if (tmp<slk[v]) slk[v] = tmp; //顺便(?)更新slk
33 }
34 return false;
35}
36int KM()
37{
1 ms(link,-1);//初始没有找到匹配,同hungary
2 ms(ly,0);
3 ms(lx,0);//因为要使得所有lx[i]+ly[j]-w[i][j]>=0,所以初始化的时候(不妨)lx[i]=max(w[i][j]),ly[i]=0;
4 for ( int i = 1 ; i <= n ; i++)
5 for ( int j = 1 ; j <= n ; j++)
6 lx[i] = max(w[i][j],lx[i]);
1 for ( int i = 1 ; i <= n ; i++)
2 {
3// ms(slk,0x3f);
4 for ( int j = 1 ; j <= n ; j++) slk[j] = inf;
5 while (true)
6 {
7 // cout<<"mio mio mio "<<endl;
8 ms(visx,false);
9 ms(visy,false);
1 if (find(i)) break;
2 int d = inf;
3 for ( int j = 1 ; j <= n ; j++)
4 {
5 if (!visy[j]&&slk[j]<d) d= slk[j]; //d取所有slk值的最小值。
6 //这样有两个作用,一个是可以保证所有所有【点标对】仍然可行。
7 //另一个是这样的d每次都可以增加一条<可行边>进入相等子图。
8 }
1 for (int j = 1 ; j <= n ; j++) if (visx[j]) lx[j]-=d;
2 for (int j = 1 ; j <= n ; j++) if (visy[j]) ly[j]+=d;else slk[j]-=d;
3 }
4 }
1 int res = 0 ;
2 for ( int i = 1 ; i <= n ; i++)
3 {
4 if (link[i]>-1)
5 {
6 res += w[link[i]][i];
7 }
8 }
9 return res;
10}
11int main()
12{
13 #ifndef ONLINE_JUDGE
14 freopen("code/in.txt","r",stdin);
15 #endif
1 while (scanf("%d",&n)!=EOF)
2 {
3 for ( int i = 1 ; i <= n ; i++)
4 for ( int j = 1 ; j <= n ; j++)
5 scanf("%d",&w[i][j]);
6 int ans = KM();
7 printf("%d\n",ans);
8 }
1 #ifndef ONLINE_JUDGE
2 fclose(stdin);
3 #endif
4 return 0;
5}