匈牙利匹配算法

18 Dec 2014

理论基础

Berge 定理:图 G 的匹配 M 是最大匹配的充要条件是 G 中不存在 M 可扩展路。

Hall 定理:设 G 是具有二划分 (X, Y) 的二部图,则 G 有饱和 X 的匹配当且仅当对 $\forall S \subseteq X, |N(S)|\ge |S|$,其中 N(S) 表示 S 的所有邻点之集。

算法步骤

输入:二部图 G=(X, Y).

输出:G 的一个最大匹配 M .

第0步:任取图 G 的一个匹配 M 设 X 中 M 非饱和点的集合为 A.

第1步:若$A=\phi$,则停止,输出当前的 M;否则,任取$x\in A$,记 S := {x},T := $\phi$,转到下一步.

第2步:若$N(S)\subseteq T$,则不存在从 x 出发的 M 增广路,令 A := A - {x},转第一步;否则,取$y\in N(S)-T$ 转下步.

第3步:若y是 M 饱和的,设 $yz\in M$,令 $S := S\bigcup {z}$,$T := T\bigcup {y}$,转第1步;否则获得一条 M 增广路 P(x, y),令 M := M$\oplus$E(P), A := A - {x, y},转第一步.

参考代码

 1 const int MAXN = 1000;
 2 int uN, vN; // u, v数目,要初始化!!!
 3 bool g[MAXN][MAXN]; // g[i][j] 表示xi与yj相连
 4 int xM[MAXN], yM[MAXN]; // 输出量
 5 bool chk[MAXN]; // 辅助量检查某轮y[v]是否被check
 6 bool SearchPath(int u){
 7     int v;
 8     for(v = 0; v < vN; v++)
 9         if(g[u][v] && !chk[v])
10         {
11             chk[v] = true;
12             if(yM[v] == -1 || SearchPath(yM[v]))
13             {
14                 yM[v] = u; xM[u] = v;
15                 return true ;
16             }
17         }
18     return false ;
19 }
20 int hungary(){
21     int u, ret = 0 ;
22     memset(xM, -1, sizeof (xM));
23     memset(yM, -1, sizeof (yM));
24     for(u = 0; u < uN; u++)
25         if(xM[u] == -1){
26             memset(chk, false, sizeof (chk));
27             if(SearchPath(u)) ret++;
28         }
29     return ret;
30 }

求赋权二部图最大匹配的 Kuhn-Munkres 算法

18 Dec 2014

最优指派问题

欲安排 n 个人员 $x_1, x_2, \cdots, x_n$ 从事 n 项工作 $y_1, y_2,\cdots, y_n$ ,已知每个人能胜任其中一项或几项工作,各人做不同工作的效率不同. 求一种工作安排,使得每个人分配一项不同于其他人的工作,且使总的工作效率最大.

图论描述

给定赋权二部图 $K_{n,n} = (X, Y) : X = \{x_1, x_2, \cdots, x_n\}, Y = \{y_1, y_2,\cdots, y_n \}$,边$x_iy_j$有权$w_{ij}$(权值为0表示不胜任工作). 求$K_{n,n}$的一个具有最大权的完美匹配.

一般地,给定一个赋权二部图G={X,Y}(未必 X = Y ),可以添加一些定点及一些权值为0的边以转化问题.

算法理论基础

  • 顶点标号与可行顶点标号

图 G 的顶点标号是从顶点到正整数集的一个映射,用l(v)表示顶点 v 的标号, w(uv)表示边(u, v)的权. 对于赋权二部图 G = {X, Y},若对每条边 e = xy,均有 $l(x)+l(y)\ge w(xy)$,则称这个标号为 G 的一个可行定点标号。

picture001

赋权二部图的可行顶点标号总是存在的. 一种平凡的可行顶点标号是:对$\forall v \in V$,

picturetes

  • 相等子图

设 G 是一个赋权二部图,l 是 G 的可行顶点标号,边 (u, v) 上的权为 w(uv),令

G 中以 $E_l$ 为边集的生成子图称为 G 的 l 相等子图,记为 $G_l$.

下面的例子显示了一个赋权完全二部图 $G=K_{4,4}$的平凡顶点标号,以及这种标号下的相等子图.

tex002

picture002 picture003

定理

设l是赋权二部图 G 的一个可行顶点标号. 若相等子图$G_l$有完美匹配$M^* $,则$M^* $ 是G的最大权完美匹配.

证明

由于相等子图 $ G_l $ 是G的生成子图,固 $ G_l $ 的完美匹配 $ M^* $ 也是 G 的完美匹配,而且

tex003

另一方面,对 G 的任何完美匹配 M,有:

tex004

可见 $W(M^* )\ge W(M)$ ,即 $M^* $ 是 G 的最大权完美匹配. 证毕.

算法思想

Kuhn和Munkres给出修改标号的方法,使得新的相等子图的最大匹配逐渐扩大,最终的到完美匹配的相等子图.

首先给出赋权二部图 G 的任意一个可行顶点标号(如平凡标号),然后决定相等子图 $ G_l $ ,在 $ G_l $ 中执行匈牙利算法. 若在 $ G_l $中找到完美匹配,则它就是 G 的最大完美匹配. 否则,匈牙利算法终止于 $ S\subset X, T\subset Y $ ,且 $ N_{G_{l}}=T $ 如下图所示.

picture004

设当前找到的匹配为 M. 令 对每个顶点u,修改其标号如下

tex005

可以检验$l’$仍是 G 的一个可行顶点标号.用 $l’$代替$l$,获得新的相等子图 $G’$.

参考代码

  1 /*
  2  * w 表示带权图
  3  * pop表示图一侧的顶点数
  4  *
  5  * 返回对大权匹配的权值
  6  * son_y 表示匹配方案
  7  */
  8 const int maxn = 222;
  9 const int inf = 1000000000;
 10 
 11 int mymin(int a, int b) {
 12     return a > b ? b : a;
 13 }
 14 int mymax(int a, int b) {
 15     return a > b ? a : b;
 16 }
 17 
 18 int w[maxn][maxn], x[maxn], y[maxn];
 19 int prev_x[maxn], prev_y[maxn], son_y[maxn], slack[maxn], par[maxn];
 20 int lx, ly, pop;
 21 void adjust(int v) {
 22     son_y[v] = prev_y[v];
 23     if(prev_x[son_y[v]] != -2)
 24         adjust(prev_x[son_y[v]]);
 25 }
 26 bool find(int v) {
 27     int i;
 28     for(i = 0; i < pop; i++)
 29         if(prev_y[i] == -1) {
 30             if(slack[i] > x[v] + y[i] - w[v][i]) {
 31                 slack[i] = x[v] + y[i] - w[v][i];
 32                 par[i] = v;
 33             }
 34             if(x[v] + y[i] == w[v][i]) {
 35                 prev_y[i] = v;
 36                 if(son_y[i] == -1) {
 37                     adjust(i);
 38                     return true;
 39                 }
 40                 if(prev_x[son_y[i]] != -1)
 41                     continue;
 42                 prev_x[son_y[i]] = i;
 43                 if(find(son_y[i]))
 44                     return true;
 45             }
 46         }
 47     return false;
 48 }
 49 int km() {
 50     int i, j, m;
 51     for(i = 0; i < pop; i++) {
 52         son_y[i] = -1;
 53         y[i] = 0;
 54     }
 55     for(i = 0; i < pop; i++) {
 56         x[i] = 0;
 57         for(j = 0; j < pop; j++)
 58             x[i] = mymax(x[i], w[i][j]);
 59     }
 60     bool flag;
 61     for(i = 0; i < pop; i++) {
 62         for(j = 0; j < pop; j++) {
 63             prev_x[j] = prev_y[j] = -1;
 64             slack[j] = inf;
 65         }
 66         prev_x[i] = -2;
 67         if(find(i)) continue;
 68         flag = false;
 69         while(!flag) {
 70             m = inf;
 71             for(j = 0; j < pop; j++)
 72                 if(prev_y[j] == -1)
 73                     m = mymin(m, slack[j]);
 74             for(j = 0; j < pop; j++) {
 75                 if(prev_x[j] != -1)
 76                     x[j] -= m;
 77                 if(prev_y[j] != -1)
 78                     y[j] += m;
 79                 else
 80                     slack[j] -= m;
 81             }
 82             for(j = 0; j < pop; j++)
 83                 if(prev_y[j] == -1 && !slack[j]) {
 84                     prev_y[j] = par[j];
 85                     if(son_y[j] == -1) {
 86                         adjust(j);
 87                         flag = true;
 88                         break;
 89                     }
 90                     prev_x[son_y[j]] = j;
 91                     if(find(son_y[j])) {
 92                         flag = true;
 93                         break;
 94                     }
 95                 }
 96         }
 97     }
 98     int ans = 0;
 99     for(int i = 0; i < pop; i++)
100         ans += w[son_y[i]][i];
101     return ans;
102 }

关于日期的计算

29 Nov 2014

给定一个日期,问这个日期是星期几。

罗马教皇格里高利13世在1582年组织了一批天文学家,根据哥白尼日心说计算出来的数据,对儒略历作了修改。将1582年10月5日到14日之间的10天宣布撤销,继10月4日之后为10月15日。后来人们将这一新的历法称为“格里高利历”,也就是今天世界上通用的历法。不同国家取消旧历法启用新历法的年代不同,导致蔡勒公式的不同版本。

第一个方法书计算这个日期与今天的距离X,假设今天是星期y,那么给定日期是星期((y - X) % 7 + 7) % 7(给定日期是今天之前),或者星期 (y + X) % 7 + 1(给定日期是未来的日期)。

第二个方法是直接使用蔡勒公式:

Week = (Day + 2 X Month + 3 X (Month + 1) / 5) + Year + Year / 4 - Year / 100 + Year / 400) % 7

当日期在1752年9月3日之前时:

Week = (Day + 2 X Month + 3 X (Month + 1) / 5) + Year + Year / 4 + 5) % 7
 1 int watDay(int d, int m, int y) {
 2     int ret;
 3     if (m == 1 || m == 2) {
 4         m += 12;
 5         y--;
 6     }
 7     if ((y < 1752) || (y == 1752 && m < 9) || (y == 1752 && m == 9 && d == 3))
 8         ret = (d + 2 * m + 3 * (m + 1) / 5 + y + y / 4 + 5) % 7;
 9     else
10         ret = (d + 2 * m + 3 * (m + 1) / 5 + y + y / 4 - y / 100 + y / 400) % 7;
11     return ret;
12 }

日期相隔天数计算

给定两个日期A,B,求A,B间隔多少天。

计算公元元年到A和B分别有多少天,然后两个值相减即可。下面的代码不考虑公元前以及历史上的日历调整。

 1 const int days = 365;
 2 const int s[] = {0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
 3 
 4 bool Isleep(int y) {
 5     if (y % 400 == 0 || y % 100 && y % 4 == 0) return 1;
 6     return 0;
 7 }
 8 int leap(int y) {
 9     if (!y) return 0;
10     return y / 4 - y / 100 + y / 400;
11 }
12 int calc(int day, int mon, int year) {
13     int res = (year - 1) * days + leap(year - 1);
14     for (int i = 1; i < mon; i++) {
15         res += s[i];
16     }
17     if (Isleep(year) && mon > 2) res++;
18     res += day;
19     return res;
20 }
21 int countDay(int da, int ma, int ya, int db, int mb, int yb) {
22     return abs(calc(da, ma, ya) - calc(db, mb, yb));
23 }