From 7541ae58a29f18fea5612c97eeb0cf12aa689622 Mon Sep 17 00:00:00 2001 From: hitonanode <32937551+hitonanode@users.noreply.github.com> Date: Fri, 31 Mar 2023 01:03:25 +0900 Subject: [PATCH 1/3] Hungarian, prototype --- combinatorial_opt/hungarian.hpp | 106 ++++++++++++++++++++++ combinatorial_opt/test/hungarian.test.cpp | 20 ++++ 2 files changed, 126 insertions(+) create mode 100644 combinatorial_opt/hungarian.hpp create mode 100644 combinatorial_opt/test/hungarian.test.cpp diff --git a/combinatorial_opt/hungarian.hpp b/combinatorial_opt/hungarian.hpp new file mode 100644 index 00000000..254e3d09 --- /dev/null +++ b/combinatorial_opt/hungarian.hpp @@ -0,0 +1,106 @@ +#pragma once +#include +#include +#include +#include + +// Solve assignment problem by Hungarian algorithm +// dual problem: maximize sum(f) - sum(g) s.t. f_i - g_j <= C_ij +// Requires: n == m +// Ccomplexity: O(n^2 m) +// https://www.slideshare.net/joisino/ss-249394573 +// Todo: +// - generalize: n != m +// - Reduce to O(nm min(n, m)) +template std::pair> hungarian(const std::vector> &C) { + if (C.empty()) return {T(), {}}; + const int n = C.size(), m = C.front().size(); + assert(n == m); + std::vector f(n, T()), g(m, T()); + + auto chmin = [](T &x, T y) { return (x > y ? (x = y, true) : false); }; + + // Make dual variables feasible + for (int j = 0; j < m; ++j) { + g[j] = f[0] - C[0][j]; + for (int i = 0; i < n; ++i) g[j] = std::max(g[j], f[i] - C[i][j]); + } + + std::vector lmate(n, -1); + std::vector> rmate(m, std::nullopt); + + std::vector rreach(m, -1); + std::vector rprv(m, -1); + + std::vector lvisited; + + for (int i = 0; i < n; ++i) { + lvisited = {i}; + int cur = 0; + std::optional reachable_r = std::nullopt; + + while (!reachable_r.has_value()) { + + auto check_l = [&]() -> void { + int l = lvisited[cur++]; + for (int j = 0; j < m; ++j) { + if (rreach[j] == i) continue; + if (f[l] - g[j] == C[l][j]) { + rreach[j] = i; + rprv[j] = l; + if (rmate[j].has_value()) { + lvisited.push_back(rmate[j].value()); + } else { + reachable_r = j; + cur = lvisited.size(); + return; + } + } + } + }; + while (cur < int(lvisited.size())) check_l(); + + if (!reachable_r.has_value()) { + T min_diff = T(); + int min_l = -1, min_r = -1; + for (int l : lvisited) { + for (int j = 0; j < m; ++j) { + if (rreach[j] == i) continue; + T diff = C[l][j] + g[j] - f[l]; + if (min_l < 0) { + min_diff = diff; + min_l = l; + min_r = j; + } else if (chmin(min_diff, diff)) { + min_l = l; + min_r = j; + } + } + } + for (int l : lvisited) f[l] += min_diff; + for (int j = 0; j < m; ++j) { + if (rreach[j] == i) g[j] += min_diff; + } + rreach[min_r] = i; + rprv.at(min_r) = min_l; + + if (rmate[min_r].has_value()) { + lvisited.push_back(rmate[min_r].value()); + } else { + reachable_r = min_r; + } + } + } + for (int h = reachable_r.value(); h >= 0;) { + int l = rprv.at(h); + int nxth = lmate.at(l); + rmate.at(h) = l; + lmate.at(l) = h; + h = nxth; + } + } + + T sol = T(); + for (int i = 0; i < n; ++i) sol += C.at(i).at(lmate.at(i)); + return {sol, lmate}; +} diff --git a/combinatorial_opt/test/hungarian.test.cpp b/combinatorial_opt/test/hungarian.test.cpp new file mode 100644 index 00000000..fd2f46a2 --- /dev/null +++ b/combinatorial_opt/test/hungarian.test.cpp @@ -0,0 +1,20 @@ +#define PROBLEM "https://judge.yosupo.jp/problem/assignment" +#include "../hungarian.hpp" + +#include +#include +using namespace std; + +int main() { + cin.tie(nullptr), ios::sync_with_stdio(false); + int n; + cin >> n; + vector> a(n, vector(n)); + for (auto &vec : a) { + for (auto &x : vec) cin >> x; + } + auto [val, sol] = hungarian(a); + cout << val << '\n'; + for (auto x : sol) cout << x << ' '; + cout << '\n'; +} From 6f6c2cd5e86413eade58a26c9dbcca97076e2021 Mon Sep 17 00:00:00 2001 From: hitonanode <32937551+hitonanode@users.noreply.github.com> Date: Fri, 31 Mar 2023 01:11:31 +0900 Subject: [PATCH 2/3] Refactor --- combinatorial_opt/hungarian.hpp | 24 ++++++++++-------------- 1 file changed, 10 insertions(+), 14 deletions(-) diff --git a/combinatorial_opt/hungarian.hpp b/combinatorial_opt/hungarian.hpp index 254e3d09..88da71e4 100644 --- a/combinatorial_opt/hungarian.hpp +++ b/combinatorial_opt/hungarian.hpp @@ -18,8 +18,6 @@ template std::pair> hungarian(const std::vector f(n, T()), g(m, T()); - auto chmin = [](T &x, T y) { return (x > y ? (x = y, true) : false); }; - // Make dual variables feasible for (int j = 0; j < m; ++j) { g[j] = f[0] - C[0][j]; @@ -37,9 +35,9 @@ template std::pair> hungarian(const std::vector reachable_r = std::nullopt; + int reachable_r = -1; - while (!reachable_r.has_value()) { + while (reachable_r < 0) { auto check_l = [&]() -> void { int l = lvisited[cur++]; @@ -60,38 +58,36 @@ template std::pair> hungarian(const std::vector diff) { min_diff = diff; min_l = l; min_r = j; - } else if (chmin(min_diff, diff)) { - min_l = l; - min_r = j; } } } + for (int l : lvisited) f[l] += min_diff; for (int j = 0; j < m; ++j) { if (rreach[j] == i) g[j] += min_diff; } - rreach[min_r] = i; + rreach.at(min_r) = i; rprv.at(min_r) = min_l; - if (rmate[min_r].has_value()) { - lvisited.push_back(rmate[min_r].value()); + if (rmate.at(min_r).has_value()) { + lvisited.push_back(rmate.at(min_r).value()); } else { reachable_r = min_r; } } } - for (int h = reachable_r.value(); h >= 0;) { + + for (int h = reachable_r; h >= 0;) { int l = rprv.at(h); int nxth = lmate.at(l); rmate.at(h) = l; From 8a1ff1c288b758f4686e32673e8dcbaff4f76a11 Mon Sep 17 00:00:00 2001 From: hitonanode <32937551+hitonanode@users.noreply.github.com> Date: Fri, 7 Apr 2023 11:07:42 +0900 Subject: [PATCH 3/3] =?UTF-8?q?Hungarian:=20n=20!=3D=20m=20=E3=81=AB?= =?UTF-8?q?=E5=AF=BE=E5=BF=9C=EF=BC=9F?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- combinatorial_opt/hungarian.hpp | 43 +++++++++++++++---- .../test/hungarian.aoj1163.test.cpp | 28 ++++++++++++ 2 files changed, 62 insertions(+), 9 deletions(-) create mode 100644 combinatorial_opt/test/hungarian.aoj1163.test.cpp diff --git a/combinatorial_opt/hungarian.hpp b/combinatorial_opt/hungarian.hpp index 88da71e4..4b8416a3 100644 --- a/combinatorial_opt/hungarian.hpp +++ b/combinatorial_opt/hungarian.hpp @@ -6,16 +6,26 @@ // Solve assignment problem by Hungarian algorithm // dual problem: maximize sum(f) - sum(g) s.t. f_i - g_j <= C_ij -// Requires: n == m -// Ccomplexity: O(n^2 m) +// Ccomplexity: O(nm min(n, m)) // https://www.slideshare.net/joisino/ss-249394573 -// Todo: -// - generalize: n != m -// - Reduce to O(nm min(n, m)) template std::pair> hungarian(const std::vector> &C) { if (C.empty()) return {T(), {}}; + const int n = C.size(), m = C.front().size(); - assert(n == m); + + if (n < m) { + std::vector transC(m, std::vector(n)); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < m; ++j) transC.at(j).at(i) = C.at(i).at(j); + } + auto j2i = hungarian(transC); + std::vector ret(n, -1); + for (int j = 0; j < m; ++j) { + if (j2i.second.at(j) >= 0) ret.at(j2i.second.at(j)) = j; + } + return {j2i.first, ret}; + } + std::vector f(n, T()), g(m, T()); // Make dual variables feasible @@ -72,6 +82,20 @@ template std::pair> hungarian(const std::vector= m); + int l = lvisited.front(); + for (int t : lvisited) { + if (f.at(l) < f.at(t)) l = t; + } + if (lmate.at(l) >= 0) { + reachable_r = lmate.at(l); + rmate.at(lmate.at(l)) = std::nullopt; + lmate.at(l) = -1; + } + break; + } + for (int l : lvisited) f[l] += min_diff; for (int j = 0; j < m; ++j) { if (rreach[j] == i) g[j] += min_diff; @@ -88,8 +112,7 @@ template std::pair> hungarian(const std::vector= 0;) { - int l = rprv.at(h); - int nxth = lmate.at(l); + int l = rprv.at(h), nxth = lmate.at(l); rmate.at(h) = l; lmate.at(l) = h; h = nxth; @@ -97,6 +120,8 @@ template std::pair> hungarian(const std::vector= 0) sol += C.at(i).at(lmate.at(i)); + } return {sol, lmate}; } diff --git a/combinatorial_opt/test/hungarian.aoj1163.test.cpp b/combinatorial_opt/test/hungarian.aoj1163.test.cpp new file mode 100644 index 00000000..be99adf6 --- /dev/null +++ b/combinatorial_opt/test/hungarian.aoj1163.test.cpp @@ -0,0 +1,28 @@ +#define PROBLEM "https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1163" +#include "../hungarian.hpp" +#include +#include +#include +using namespace std; + +int main() { + cin.tie(nullptr), ios::sync_with_stdio(false); + + while (true) { + int m, n; + cin >> m >> n; + if (!m) break; + + vector b(m), r(n); + for (auto &x : b) cin >> x; + for (auto &x : r) cin >> x; + + vector> mat; + for (int x : b) { + mat.push_back({}); + for (int y : r) mat.back().push_back(std::gcd(x, y) > 1 ? -1 : 0); + } + + cout << -hungarian(mat).first << '\n'; + } +}