KM算法

　　昨晚纠结了半天是看KM算法还是复习软工，最后顶不住挂课的压力复习软工，今天考试还算可以，希望不要挂掉。。。上午上课+考试，感觉比连打几场比赛还累。T_T

下午队内DIY，中途停电。。。-_-///然后回来开始研究KM。。。

找了很多资料，下面这几个算是让我理解一点的吧，感谢大牛！

http://cuitianyi.com/blog/%E6%B1%82%E6%9C%80%E5%A4%A7%E6%9D%83%E4%BA%8C%E5%88%86%E5%8C%B9%E9%85%8D%E7%9A%84km%E7%AE%97%E6%B3%95/

http://www.cppblog.com/MatoNo1/archive/2012/04/26/151724.html

http://blog.csdn.net/liguanxing/article/details/5665646

 

KM算法过程：

1、初始化可行性顶标lx[], ly[]。

2、用类似Hungry算法的思想求完全匹配。

3、找不到则调整lx[], ly[]的值，然后回到1。

4、如果二分图已经是完全匹配的则退出，ans = ∑(lx[i] + ly[i]);

写点自己的理解：

顶标：对一个点i，lx[i]或则ly[i]，表示在i点上做的一个标志，当然，这个标志是有规律的，见下面；

关于可行性顶标：已知i, j的权值为w[i][j]，lx[i] + ly[j] >= w[i][j]。则lx[i],ly[j]这一组顶标是可行的；

完全匹配：在二分图X集合或者Y集合中，X集合中的点都对应匹配或者Y集合中的点都对应匹配；

如何实现二分图的完全匹配：如果所有满足lx[i] + ly[j] == w[i][j]的边组导出一个完全匹配图来，则这个图是二分图的完全匹配；

具体实现：（摘自）

当前的顶标的导出子图并不一定存在完美匹配。这时，可以用某种方法对顶标进行调整。调整的方法是：根据最后一次不成功的寻找交错路的DFS，取所有i被访问到而j没被访问到的边(i,j)的lx[i]+ly[j]-w[i][j]的最小值d。将交错树中的所有左端点的顶标减小d，右端点的顶标增加d。经过这样的调整以后：原本在导出子图里面的边，两边的顶标都变了，不等式的等号仍然成立，仍然在导出子图里面；原本不在导出子图里面的边，它的左端点的顶标减小了，右端点的顶标没有变，而且由于d的定义，不等式仍然成立，所以他就可能进入了导出子图里。

初始时随便指定一个可行顶标，比如说lx[i]=max{w[i][j]|j是右边的点}，ly[i]=0。然后对每个顶点进行类似Hungary算法的find过程，如果某次find没有成功，则按照这次find访问到的点对可行顶标进行上述调整。这样就可以逐步找到完美匹配了。

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就可以。

 算法模板：

//这里写的是一个求最小权匹配的，把原图的权值转成负数了
//POJ 2195

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>

#define REP(i, n)   for(i = 0; i < n; ++i)
#define FOR(i, l, h)    for(i = l; i <= h; ++i)
#define FORD(i, h, l)   for(i = h; i >= l; --i)
#define CL(arr, val)    memset(arr, val, sizeof(arr))

using namespace std;

const int N = 1024;
const int inf = ~0u>>2;

int x[N], y[N];
int lx[N];
int ly[N];
int mp[N][N];
int num[N][N];
char c[N][N];
int linky[N];
int slack[N];
bool visx[N];
bool visy[N];
int n, cnt;

int Abs(int x) {
    return x < 0 ? -x : x;
}

bool dfs(int x) {
    int y, d;
    visx[x] = true;
    REP(y, cnt) {
        d = lx[x] + ly[y] - mp[x][y];
        if(d == 0 && !visy[y]) {
            visy[y] = true;
            if(linky[y] == -1 || dfs(linky[y])) {
                linky[y] = x;
                return true;
            }
        } else slack[y] = min(slack[y], d);
    }
    return false;
}

void KM() {
    int i, j, k, ans = 0, d;
    REP(i, cnt) {
        ly[i] = 0; linky[i] = -1; lx[i] = -inf;
        REP(j, cnt) {
            lx[i] = max(lx[i], mp[i][j]);
        }
    }
    REP(k, cnt) {
        CL(visx, false); 
        CL(visy, false);
        REP(i, cnt)    slack[i] = inf;
        while(!dfs(k)) {
            d = inf;
            REP(i, cnt) if(!visy[i]) d = min(d, slack[i]);
            REP(i, cnt) {
                if(visx[i]) {visx[i] = false; lx[i] -= d;}
                if(visy[i]) {visy[i] = false; ly[i] += d;}
            }
        }
    }
    REP(i, cnt)    ans += lx[i] + ly[i];
    ans /= -2;
    cout << ans << endl;
}

void build(int n, int m) {
    int i, j;
    cnt = 0;
    CL(c, 0);
    REP(i, n) {
        scanf("%s", c[i]);
        REP(j, m) {
            if(c[i][j] == 'H') {
                x[cnt] = i; y[cnt++] = j;
            }
        }
    }
    int l = cnt;
    REP(i, n) {
        REP(j, m) {
            if(c[i][j] == 'm') {
                x[cnt] = i; y[cnt++] = j;
            }
        }
    }
    REP(i, cnt) REP(j, cnt) mp[i][j] = -inf;
    REP(i, l) {
        FOR(j, l, cnt - 1) {
            mp[j][i] = mp[i][j] = -(Abs(x[i] - x[j]) + Abs(y[i] - y[j]));
        }
    }
}

int main() {
    //freopen("data.in", "r", stdin);

    int n, m;
    while(~scanf("%d%d", &n, &m)) {
        if(n + m == 0)  break;
        build(n, m);
        KM();
    }
    return 0;
}