123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105 |
- // Min cost bipartite matching via shortest augmenting paths
- //
- // This is an O(n^3) implementation of a shortest augmenting path
- // algorithm for finding min cost perfect matchings in dense
- // graphs. In practice, it solves 1000x1000 problems in around 1
- // second.
- //
- // cost[i][j] = cost for pairing left node i with right node j
- // Lmate[i] = index of right node that left node i pairs with
- // Rmate[j] = index of left node that right node j pairs with
- //
- // The values in cost[i][j] may be positive or negative. To perform
- // maximization, simply negate the cost[][] matrix.
- typedef vector<double> VD;
- typedef vector<VD> VVD;
- typedef vector<int> VI;
- double MinCostMatching(const VVD &cost, VI &Lmate, VI &Rmate) {
- int n = int(cost.size());
- // construct dual feasible solution
- VD u(n);
- VD v(n);
- for (int i = 0; i < n; i++) {
- u[i] = cost[i][0];
- for (int j = 1; j < n; j++) u[i] = min(u[i], cost[i][j]);
- }
- for (int j = 0; j < n; j++) {
- v[j] = cost[0][j] - u[0];
- for (int i = 1; i < n; i++) v[j] = min(v[j], cost[i][j] - u[i]);
- }
- // construct primal solution satisfying complementary slackness
- Lmate = VI(n, -1);
- Rmate = VI(n, -1);
- int mated = 0;
- for (int i = 0; i < n; i++) {
- for (int j = 0; j < n; j++) {
- if (Rmate[j] != -1) continue;
- if (fabs(cost[i][j] - u[i] - v[j]) < 1e-10) {
- Lmate[i] = j;
- Rmate[j] = i;
- mated++;
- break;
- }
- }
- }
- VD dist(n);
- VI dad(n);
- VI seen(n);
- // repeat until primal solution is feasible
- while (mated < n) {
- // find an unmatched left node
- int s = 0;
- while (Lmate[s] != -1) s++;
- // initialize Dijkstra
- fill(dad.begin(), dad.end(), -1);
- fill(seen.begin(), seen.end(), 0);
- for (int k = 0; k < n; k++)
- dist[k] = cost[s][k] - u[s] - v[k];
- int j = 0;
- while (true) {
- // find closest
- j = -1;
- for (int k = 0; k < n; k++) {
- if (seen[k]) continue;
- if (j == -1 || dist[k] < dist[j]) j = k;
- }
- seen[j] = 1;
- // termination condition
- if (Rmate[j] == -1) break;
- // relax neighbors
- const int i = Rmate[j];
- for (int k = 0; k < n; k++) {
- if (seen[k]) continue;
- const double new_dist = dist[j] + cost[i][k] - u[i] - v[k];
- if (dist[k] > new_dist) {
- dist[k] = new_dist;
- dad[k] = j;
- }
- }
- }
- // update dual variables
- for (int k = 0; k < n; k++) {
- if (k == j || !seen[k]) continue;
- const int i = Rmate[k];
- v[k] += dist[k] - dist[j];
- u[i] -= dist[k] - dist[j];
- }
- u[s] += dist[j];
- // augment along path
- while (dad[j] >= 0) {
- const int d = dad[j];
- Rmate[j] = Rmate[d];
- Lmate[Rmate[j]] = j;
- j = d;
- }
- Rmate[j] = s;
- Lmate[s] = j;
- mated++;
- }
- double value = 0;
- for (int i = 0; i < n; i++)
- value += cost[i][Lmate[i]];
- return value;
- }
|