-
Notifications
You must be signed in to change notification settings - Fork 3
/
MinCostBipartiteMatching.cc
136 lines (116 loc) · 3.66 KB
/
MinCostBipartiteMatching.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
#include <algorithm>
#include <cstdio>
#include <cmath>
#include <vector>
using namespace std;
// BEGIN
///////////////////////////////////////////////////////////////////////////
// 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;
}
// END
int main() {
//TODO: write tests
return 0;
}