This documentation is automatically generated by online-judge-tools/verification-helper
View the Project on GitHub ei1333/library
#include "graph/flow/hungarian.hpp"
二部グラフの最小重み最大マッチングをハンガリアン法により求める. 特に $N=M$ (左側の頂点数と右側の頂点数が一致する) の場合, 割当問題とも言われる.
hungarian(A)
A
first
second
#include "../../math/matrix/matrix.hpp" /** * @brief Hungarian(二部グラフの最小重み最大マッチング) * @docs docs/hungarian.md */ template< typename T > pair< T, vector< int > > hungarian(Matrix< T > &A) { const T infty = numeric_limits< T >::max(); const int N = (int) A.height(); const int M = (int) A.width(); vector< int > P(M), way(M); vector< T > U(N, 0), V(M, 0), minV; vector< bool > used; for(int i = 1; i < N; i++) { P[0] = i; minV.assign(M, infty); used.assign(M, false); int j0 = 0; while(P[j0] != 0) { int i0 = P[j0], j1 = 0; used[j0] = true; T delta = infty; for(int j = 1; j < M; j++) { if(used[j]) continue; T curr = A[i0][j] - U[i0] - V[j]; if(curr < minV[j]) minV[j] = curr, way[j] = j0; if(minV[j] < delta) delta = minV[j], j1 = j; } for(int j = 0; j < M; j++) { if(used[j]) U[P[j]] += delta, V[j] -= delta; else minV[j] -= delta; } j0 = j1; } do { P[j0] = P[way[j0]]; j0 = way[j0]; } while(j0 != 0); } return {-V[0], P}; }
#line 1 "math/matrix/matrix.hpp" template< class T > struct Matrix { vector< vector< T > > A; Matrix() {} Matrix(size_t n, size_t m) : A(n, vector< T >(m, 0)) {} Matrix(size_t n) : A(n, vector< T >(n, 0)) {}; size_t size() const { if(A.empty()) return 0; assert(A.size() == A[0].size()); return A.size(); } size_t height() const { return (A.size()); } size_t width() const { return (A[0].size()); } inline const vector< T > &operator[](int k) const { return (A.at(k)); } inline vector< T > &operator[](int k) { return (A.at(k)); } static Matrix I(size_t n) { Matrix mat(n); for(int i = 0; i < n; i++) mat[i][i] = 1; return (mat); } Matrix &operator+=(const Matrix &B) { size_t n = height(), m = width(); assert(n == B.height() && m == B.width()); for(int i = 0; i < n; i++) for(int j = 0; j < m; j++) (*this)[i][j] += B[i][j]; return (*this); } Matrix &operator-=(const Matrix &B) { size_t n = height(), m = width(); assert(n == B.height() && m == B.width()); for(int i = 0; i < n; i++) for(int j = 0; j < m; j++) (*this)[i][j] -= B[i][j]; return (*this); } Matrix &operator*=(const Matrix &B) { size_t n = height(), m = B.width(), p = width(); assert(p == B.height()); vector< vector< T > > C(n, vector< T >(m, 0)); for(int i = 0; i < n; i++) for(int j = 0; j < m; j++) for(int k = 0; k < p; k++) C[i][j] = (C[i][j] + (*this)[i][k] * B[k][j]); A.swap(C); return (*this); } Matrix &operator^=(long long k) { Matrix B = Matrix::I(height()); while(k > 0) { if(k & 1) B *= *this; *this *= *this; k >>= 1LL; } A.swap(B.A); return (*this); } Matrix operator+(const Matrix &B) const { return (Matrix(*this) += B); } Matrix operator-(const Matrix &B) const { return (Matrix(*this) -= B); } Matrix operator*(const Matrix &B) const { return (Matrix(*this) *= B); } Matrix operator^(const long long k) const { return (Matrix(*this) ^= k); } friend ostream &operator<<(ostream &os, Matrix &p) { size_t n = p.height(), m = p.width(); for(int i = 0; i < n; i++) { os << "["; for(int j = 0; j < m; j++) { os << p[i][j] << (j + 1 == m ? "]\n" : ","); } } return (os); } T determinant() { Matrix B(*this); assert(width() == height()); T ret = 1; for(int i = 0; i < width(); i++) { int idx = -1; for(int j = i; j < width(); j++) { if(B[j][i] != 0) idx = j; } if(idx == -1) return (0); if(i != idx) { ret *= -1; swap(B[i], B[idx]); } ret *= B[i][i]; T vv = B[i][i]; for(int j = 0; j < width(); j++) { B[i][j] /= vv; } for(int j = i + 1; j < width(); j++) { T a = B[j][i]; for(int k = 0; k < width(); k++) { B[j][k] -= B[i][k] * a; } } } return (ret); } }; #line 2 "graph/flow/hungarian.hpp" /** * @brief Hungarian(二部グラフの最小重み最大マッチング) * @docs docs/hungarian.md */ template< typename T > pair< T, vector< int > > hungarian(Matrix< T > &A) { const T infty = numeric_limits< T >::max(); const int N = (int) A.height(); const int M = (int) A.width(); vector< int > P(M), way(M); vector< T > U(N, 0), V(M, 0), minV; vector< bool > used; for(int i = 1; i < N; i++) { P[0] = i; minV.assign(M, infty); used.assign(M, false); int j0 = 0; while(P[j0] != 0) { int i0 = P[j0], j1 = 0; used[j0] = true; T delta = infty; for(int j = 1; j < M; j++) { if(used[j]) continue; T curr = A[i0][j] - U[i0] - V[j]; if(curr < minV[j]) minV[j] = curr, way[j] = j0; if(minV[j] < delta) delta = minV[j], j1 = j; } for(int j = 0; j < M; j++) { if(used[j]) U[P[j]] += delta, V[j] -= delta; else minV[j] -= delta; } j0 = j1; } do { P[j0] = P[way[j0]]; j0 = way[j0]; } while(j0 != 0); } return {-V[0], P}; }