Posted on 2011-07-23 22:14
Mato_No1 閱讀(7619)
評論(9) 編輯 收藏 引用 所屬分類:
圖算法
【背景（神犇不要鄙視）】
本沙茶今年AHOI的時候，遇到裸的最佳匹配的題，竟然把KM算法搞忘了，幸虧是WJMZBMR神犇保佑，臨時亂弄一通，想起來了……這MS反映出了本沙茶以前在看某些經典算法的時候看得不深，木有理解透徹……
前幾天又遇到一道最佳匹配的題，發現KM算法竟然又忘了……米辦法，只有把這個搞死人的算法的具體過程重新看了一遍，終于懂了……
【KM算法及其具體過程】
（1）可行點標：每個點有一個標號，記lx[i]為X方點i的標號，ly[j]為Y方點j的標號。如果對于圖中的任意邊(i, j, W)都有lx[i]+ly[j]>=W，則這一組點標是可行的。特別地，對于lx[i]+ly[j]=W的邊(i, j, W)，稱為可行邊；
（2）KM算法的核心思想就是通過修改某些點的標號（但要滿足點標始終是可行的），不斷增加圖中的可行邊總數，直到圖中存在僅由可行邊組成的完全匹配為止，此時這個匹配一定是最佳的（因為由可行點標的的定義，圖中的任意一個完全匹配，其邊權總和均不大于所有點的標號之和，而僅由可行邊組成的完全匹配的邊權總和等于所有點的標號之和，故這個匹配是最佳的）。一開始，求出每個點的初始標號：lx[i]=max{e.W|e.x=i}（即每個X方點的初始標號為與這個X方點相關聯的權值最大的邊的權值），ly[j]=0（即每個Y方點的初始標號為0）。這個初始點標顯然是可行的，并且，與任意一個X方點關聯的邊中至少有一條可行邊；
（3）然后，從每個X方點開始DFS增廣。DFS增廣的過程與最大匹配的Hungary算法基本相同，只是要注意兩點：一是只找可行邊，二是要把搜索過程中遍歷到的X方點全部記下來（可以用vst搞一下），以進行后面的修改；
（4）增廣的結果有兩種：若成功（找到了增廣軌），則該點增廣完成，進入下一個點的增廣。若失敗（沒有找到增廣軌），則需要改變一些點的標號，使得圖中可行邊的數量增加。方法為：將所有在增廣軌中（就是在增廣過程中遍歷到）的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])值不變，也就是這條邊的可行性不變。
這樣，在進行了這一步修改操作后，圖中原來的可行邊仍可行，而原來不可行的邊現在則可能變為可行邊。那么d的值應取多少？顯然，整個點標不能失去可行性，也就是對于上述的第<2>類邊，其lx[i]+ly[j]>=W這一性質不能被改變，故取所有第<2>類邊的(lx[i]+ly[j]-W)的最小值作為d值即可。這樣一方面可以保證點標的可行性，另一方面，經過這一步后，圖中至少會增加一條可行邊。
（5）修改后，繼續對這個X方點DFS增廣，若還失敗則繼續修改，直到成功為止；
下面分析整個算法的時間復雜度：每次修改后，圖中至少會增加一條可行邊，故最多增廣M次、修改M次就可以找到僅由可行邊組成的完全匹配（除非圖中不存在完全匹配，這個可以通過預處理得到），故整個算法的時間復雜度為O(M * (N + 一次修改點標的時間))。而一次修改點標的時間取決于計算d值的時間，如果暴力枚舉計算，這一步的時間為O(M)，優化：可以對每個Y方點設立一個slk值，表示在DFS增廣過程中，所有搜到的與該Y方點關聯的邊的(lx+ly-W)的最小值（這樣的邊的X方點必然在增廣軌中）。每次DFS增廣前，將所有Y方點的slk值設為+∞，若增廣失敗，則取所有不在增廣軌中的Y方點的slk值的最小值為d值。這樣一次修改點標的時間降為O(N)，總時間復雜度降為O(NM)。
需要注意的一點是，在增廣過程中需要記下每個X、Y方點是否被遍歷到，即fx[i]、fy[j]。因此，在每次增廣前（不是對每個X方點增廣前）就要將所有fx和fy值清空。
代碼：
void init_d()
{
re(i, n) E[i].pre = E[i].next = i; m = n;
}
void add_edge(int a, int b, int len)
{
E[m].a = a; E[m].b = b; E[m].len = len; E[m].pre = E[a].pre; E[m].next = a; E[a].pre = m; E[E[m].pre].next = m++;
}
inline int dist(int x, int y, int x0, int y0)
{
return abs(x - x0) + abs(y - y0);
}
bool dfs(int x)
{
int y, x0; fx[x] = _FLAG;
for (int p=E[x].next; p != x; p=E[p].next) {
y = E[p].b;
if (lx[x] + ly[y] > E[p].len) {
if (lx[x] + ly[y] - E[p].len < slk[y]) slk[y] = lx[x] + ly[y] - E[p].len;
} else if (fy[y] != _FLAG) {
fy[y] = _FLAG; x0 = z[y];
if (x0 == -1 || dfs(x0)) {z[y] = x; return 1;}
}
}
return 0;
}
void solve()
{
re(i, n) {
lx[i] = -INF;
for (int p=E[i].next; p != i; p=E[p].next) if (E[p].len > lx[i]) lx[i] = E[p].len;
}
re(i, n0) {ly[i] = 0; z[i] = -1;}
int d;
re(i, n) {
re(j, n0) slk[i] = INF; _FLAG++;
while (!dfs(i)) {
d = INF; re(j, n0) if (fy[j] != _FLAG && slk[j] < d) d = slk[j];
re(j, n) if (fx[j] == _FLAG) lx[j] -= d;
re(j, n0) {slk[j] = INF; if (fy[j] == _FLAG) ly[j] += d;}
_FLAG++;
}
}
res = 0; re(i, n) res += lx[i]; re(i, n0) res += ly[i];
}