From d05b501b1d4122f22b957e508e4a3b43cdec6446 Mon Sep 17 00:00:00 2001 From: Horace He Date: Mon, 13 Apr 2020 08:38:25 -0700 Subject: [PATCH] Updated Hungarian (#165) * Fixed description * Updated weighted matching * Added weighted matching stress tests * updated desc --- content/graph/WeightedMatching.h | 118 ++++++++---------------- stress-tests/graph/WeightedMatching.cpp | 48 ++++++++++ 2 files changed, 86 insertions(+), 80 deletions(-) create mode 100644 stress-tests/graph/WeightedMatching.cpp diff --git a/content/graph/WeightedMatching.h b/content/graph/WeightedMatching.h index 7421a53..dde02c1 100644 --- a/content/graph/WeightedMatching.h +++ b/content/graph/WeightedMatching.h @@ -1,89 +1,47 @@ /** - * Author: Stanford - * Date: Unknown - * Source: Stanford Notebook - * Description: Min cost bipartite matching. Negate costs for max cost. - * Time: O(N^3) - * Status: tested during ICPC 2015 + * Author: Benjamin Qi, Chilli + * Date: 2020-04-04 + * License: CC0 + * Source: https://github.com/bqi343/USACO/blob/master/Implementations/content/graphs%20(12)/Matching/Hungarian.h + * Description: Given array of (possibly negative) costs to complete $N$ jobs + * w/ $M$ workers $(N \le M)$, finds min cost to complete all jobs s.t. each + * worker is assigned to at most one job. Takes cost[N][M], where cost[i][j] = + * cost for i-th job to be completed by j-th worker and returns (min cost, + * match), where match[i] = worker assigned to i-th job. Negate costs for max + * cost. + * Time: O(N^2M) + * Status: Tested on kattis:cordonbleu, stress-tested */ #pragma once -typedef vector vd; -bool zero(double x) { return fabs(x) < 1e-10; } -double minCostMatching(const vector& cost, vi& L, vi& R) { - int n = sz(cost), mated = 0; - vd dist(n), u(n), v(n); - vi dad(n), seen(n); - - /// construct dual feasible solution - rep(i,0,n) { - u[i] = cost[i][0]; - rep(j,1,n) u[i] = min(u[i], cost[i][j]); - } - rep(j,0,n) { - v[j] = cost[0][j] - u[0]; - rep(i,1,n) v[j] = min(v[j], cost[i][j] - u[i]); - } - - /// find primal solution satisfying complementary slackness - L = R = vi(n, -1); - rep(i,0,n) rep(j,0,n) { - if (R[j] != -1) continue; - if (zero(cost[i][j] - u[i] - v[j])) { - L[i] = j; - R[j] = i; - mated++; - break; - } - } - - for (; mated < n; mated++) { // until solution is feasible - int s = 0; - while (L[s] != -1) s++; - fill(all(dad), -1); - fill(all(seen), 0); - rep(k,0,n) - dist[k] = cost[s][k] - u[s] - v[k]; - - int j = 0; - for (;;) { /// find closest - j = -1; - rep(k,0,n){ - if (seen[k]) continue; - if (j == -1 || dist[k] < dist[j]) j = k; +pair hungarian(const vector &a) { + if (a.empty()) return {0, {}}; + int n = sz(a) + 1, m = sz(a[0]) + 1; + vi u(n), v(m), p(m), ans(n - 1); + rep(i,1,n) { + p[0] = i; + int j0 = 0; // add "dummy" worker 0 + vi dist(m, INT_MAX), pre(m, -1); + vector done(m + 1); + do { // dijkstra + done[j0] = true; + int i0 = p[j0], j1, delta = INT_MAX; + rep(j,1,m) if (!done[j]) { + auto cur = a[i0 - 1][j - 1] - u[i0] - v[j]; + if (cur < dist[j]) dist[j] = cur, pre[j] = j0; + if (dist[j] < delta) delta = dist[j], j1 = j; } - seen[j] = 1; - int i = R[j]; - if (i == -1) break; - rep(k,0,n) { /// relax neighbors - if (seen[k]) continue; - auto new_dist = dist[j] + cost[i][k] - u[i] - v[k]; - if (dist[k] > new_dist) { - dist[k] = new_dist; - dad[k] = j; - } + rep(j,0,m) { + if (done[j]) u[p[j]] += delta, v[j] -= delta; + else dist[j] -= delta; } + j0 = j1; + } while (p[j0]); + while (j0) { // update alternating path + int j1 = pre[j0]; + p[j0] = p[j1], j0 = j1; } - - /// update dual variables - rep(k,0,n) { - if (k == j || !seen[k]) continue; - auto w = dist[k] - dist[j]; - v[k] += w, u[R[k]] -= w; - } - u[s] += dist[j]; - - /// augment along path - while (dad[j] >= 0) { - int d = dad[j]; - R[j] = R[d]; - L[R[j]] = j; - j = d; - } - R[j] = s; - L[s] = j; } - auto value = vd(1)[0]; - rep(i,0,n) value += cost[i][L[i]]; - return value; + rep(j,1,m) if (p[j]) ans[p[j] - 1] = j - 1; + return {-v[0], ans}; // min cost } diff --git a/stress-tests/graph/WeightedMatching.cpp b/stress-tests/graph/WeightedMatching.cpp new file mode 100644 index 0000000..f720ba9 --- /dev/null +++ b/stress-tests/graph/WeightedMatching.cpp @@ -0,0 +1,48 @@ +#include "../utilities/template.h" +#include "../utilities/utils.h" +#include "../utilities/random.h" + +#include "../../content/graph/WeightedMatching.h" +#include /// include-line, keep-include +#include "../../content/graph/MinCostMaxFlow.h" + +void test(int N, int mxCost, int iters) { + for (int it = 0; it < iters; it++) { + int n = randRange(1, N), m = randRange(1, N); + if (n > m) + swap(n, m); + + MCMF mcmf(n + m + 2); + int s = 0; + int t = 1; + for (int i = 0; i < n; i++) + mcmf.addEdge(s, i + 2, 1, 0); + for (int i = 0; i < m; i++) + mcmf.addEdge(2 + n + i, t, 1, 0); + + vector cost(n, vi(m)); + for (int i = 0; i < n; i++) { + for (int j = 0; j < m; j++) { + cost[i][j] = randRange(-mxCost, mxCost); + mcmf.addEdge(i + 2, 2 + n + j, 1, cost[i][j]); + } + } + mcmf.setpi(s); + auto maxflow = mcmf.maxflow(s, t); + auto matching = hungarian(cost); + assert(maxflow.first == n); + assert(maxflow.second == matching.first); + int matchSum = 0; + for (int i = 0; i < n; i++) + matchSum += cost[i][matching.second[i]]; + assert(matchSum == matching.first); + return; + } +} +signed main() { + test(25, 5, 1000); + test(100, 1000, 100); + test(100, 1, 50); + test(5, 5, 10000); + cout << "Tests passed!" << endl; +}