hungarian.hpp

SECTIONGraph INCLUDEnoya/graph/hungarian.hpp

求 \(n\)≤\(m\) 矩形指派的最小或最大费用，返回匹配与势能；直接读取费用矩阵，辅助空间 O(\(n\)+\(m\))。

Complexity: Time: \(O(n^{2} m)\) for an \(n\)-by-\(m\) assignment instance. Space: \(O(n + m)\) auxiliary, excluding the input matrix.

Source: maspy

AC 记录：assignment。

跳到代码 · GitHub ↗

Implementation

ReferenceSourceExpanded

当前头文件，省略 include guard；依赖见 #include。

/// @complexity Time: O(n^2 m) for an n-by-m assignment instance.
/// Space: O(n + m) auxiliary, excluding the input matrix.

#include <cassert>
#include <tuple>
#include <utility>
#include <vector>

namespace noya {

/// @brief Assign every row to a distinct column, minimizing/maximizing cost.
/// @return (cost, match, X-potential, Y-potential). Potentials refer to C for
/// MIN=true and -C otherwise. C is unchanged; all rows have the same size,
/// n <= m, and signed costs/intermediate arithmetic must be representable.
template <typename T, bool MIN>
std::tuple<T, std::vector<int>, std::vector<T>, std::vector<T>>
hungarian(std::vector<std::vector<T>> &C) {
  int N = int(C.size());
  if (N == 0) {
    return {};
  }
  int M = int(C[0].size());
  assert(N <= M);
  for (const auto &row : C) {
    assert(int(row.size()) == M);
  }
  std::vector<int> P(M + 1), way(M + 1), used;
  std::vector<T> X(N + 1), Y(M + 1), mn(M + 1);
  std::vector<unsigned char> vis(M + 1);
  T res{};
  for (int i = 1; i <= N; i++) {
    P[0] = i;
    used.clear();
    vis.assign(M + 1, false);
    int j0 = 0;
    T last{};
    while (P[j0] != 0) {
      if (j0 != 0) {
        vis[j0] = true;
        used.push_back(j0);
      }
      int i0 = P[j0], j1 = -1;
      const auto &row = C[i0 - 1];
      T dlt{};
      for (int j = 1; j <= M; j++) {
        if (vis[j]) {
          continue;
        }
        T cur = row[j - 1];
        if constexpr (!MIN) {
          cur = -cur;
        }
        cur = cur - X[i0] - Y[j];
        // Fold the previous slack decrease into this scan. The first scan
        // initializes all slacks directly, without an infinity sentinel.
        if (j0 != 0) {
          mn[j] -= last;
        }
        if (j0 == 0 || cur < mn[j]) {
          mn[j] = cur;
          way[j] = j0;
        }
        if (j1 == -1 || mn[j] < dlt) {
          dlt = mn[j];
          j1 = j;
        }
      }
      res += dlt;
      X[i] += dlt;
      for (int j : used) {
        X[P[j]] += dlt;
        Y[j] -= dlt;
      }
      last = dlt;
      j0 = j1;
    }
    do {
      P[j0] = P[way[j0]];
      j0 = way[j0];
    } while (j0 != 0);
  }
  std::vector<int> mat(N);
  for (int j = 1; j <= M; j++) {
    if (P[j] != 0) {
      mat[P[j] - 1] = j - 1;
    }
  }
  X.erase(X.begin());
  Y.erase(Y.begin());
  if constexpr (!MIN) {
    res = -res;
  }
  return {res, std::move(mat), std::move(X), std::move(Y)};
}
} // namespace noya

 // NOYA_HUNGARIAN_HPP

#ifndef NOYA_HUNGARIAN_HPP
#define NOYA_HUNGARIAN_HPP 1

/// @complexity Time: O(n^2 m) for an n-by-m assignment instance.
/// Space: O(n + m) auxiliary, excluding the input matrix.

#include <cassert>
#include <tuple>
#include <utility>
#include <vector>

namespace noya {

/// @brief Assign every row to a distinct column, minimizing/maximizing cost.
/// @return (cost, match, X-potential, Y-potential). Potentials refer to C for
/// MIN=true and -C otherwise. C is unchanged; all rows have the same size,
/// n <= m, and signed costs/intermediate arithmetic must be representable.
template <typename T, bool MIN>
std::tuple<T, std::vector<int>, std::vector<T>, std::vector<T>>
hungarian(std::vector<std::vector<T>> &C) {
  int N = int(C.size());
  if (N == 0) {
    return {};
  }
  int M = int(C[0].size());
  assert(N <= M);
  for (const auto &row : C) {
    assert(int(row.size()) == M);
  }
  std::vector<int> P(M + 1), way(M + 1), used;
  std::vector<T> X(N + 1), Y(M + 1), mn(M + 1);
  std::vector<unsigned char> vis(M + 1);
  T res{};
  for (int i = 1; i <= N; i++) {
    P[0] = i;
    used.clear();
    vis.assign(M + 1, false);
    int j0 = 0;
    T last{};
    while (P[j0] != 0) {
      if (j0 != 0) {
        vis[j0] = true;
        used.push_back(j0);
      }
      int i0 = P[j0], j1 = -1;
      const auto &row = C[i0 - 1];
      T dlt{};
      for (int j = 1; j <= M; j++) {
        if (vis[j]) {
          continue;
        }
        T cur = row[j - 1];
        if constexpr (!MIN) {
          cur = -cur;
        }
        cur = cur - X[i0] - Y[j];
        // Fold the previous slack decrease into this scan. The first scan
        // initializes all slacks directly, without an infinity sentinel.
        if (j0 != 0) {
          mn[j] -= last;
        }
        if (j0 == 0 || cur < mn[j]) {
          mn[j] = cur;
          way[j] = j0;
        }
        if (j1 == -1 || mn[j] < dlt) {
          dlt = mn[j];
          j1 = j;
        }
      }
      res += dlt;
      X[i] += dlt;
      for (int j : used) {
        X[P[j]] += dlt;
        Y[j] -= dlt;
      }
      last = dlt;
      j0 = j1;
    }
    do {
      P[j0] = P[way[j0]];
      j0 = way[j0];
    } while (j0 != 0);
  }
  std::vector<int> mat(N);
  for (int j = 1; j <= M; j++) {
    if (P[j] != 0) {
      mat[P[j] - 1] = j - 1;
    }
  }
  X.erase(X.begin());
  Y.erase(Y.begin());
  if constexpr (!MIN) {
    res = -res;
  }
  return {res, std::move(mat), std::move(X), std::move(Y)};
}
} // namespace noya

#endif // NOYA_HUNGARIAN_HPP

/// @complexity Time: O(n^2 m) for an n-by-m assignment instance.
/// Space: O(n + m) auxiliary, excluding the input matrix.

#include <cassert>
#include <tuple>
#include <utility>
#include <vector>

namespace noya {

/// @brief Assign every row to a distinct column, minimizing/maximizing cost.
/// @return (cost, match, X-potential, Y-potential). Potentials refer to C for
/// MIN=true and -C otherwise. C is unchanged; all rows have the same size,
/// n <= m, and signed costs/intermediate arithmetic must be representable.
template <typename T, bool MIN>
std::tuple<T, std::vector<int>, std::vector<T>, std::vector<T>>
hungarian(std::vector<std::vector<T>> &C) {
  int N = int(C.size());
  if (N == 0) {
    return {};
  }
  int M = int(C[0].size());
  assert(N <= M);
  for (const auto &row : C) {
    assert(int(row.size()) == M);
  }
  std::vector<int> P(M + 1), way(M + 1), used;
  std::vector<T> X(N + 1), Y(M + 1), mn(M + 1);
  std::vector<unsigned char> vis(M + 1);
  T res{};
  for (int i = 1; i <= N; i++) {
    P[0] = i;
    used.clear();
    vis.assign(M + 1, false);
    int j0 = 0;
    T last{};
    while (P[j0] != 0) {
      if (j0 != 0) {
        vis[j0] = true;
        used.push_back(j0);
      }
      int i0 = P[j0], j1 = -1;
      const auto &row = C[i0 - 1];
      T dlt{};
      for (int j = 1; j <= M; j++) {
        if (vis[j]) {
          continue;
        }
        T cur = row[j - 1];
        if constexpr (!MIN) {
          cur = -cur;
        }
        cur = cur - X[i0] - Y[j];
        // Fold the previous slack decrease into this scan. The first scan
        // initializes all slacks directly, without an infinity sentinel.
        if (j0 != 0) {
          mn[j] -= last;
        }
        if (j0 == 0 || cur < mn[j]) {
          mn[j] = cur;
          way[j] = j0;
        }
        if (j1 == -1 || mn[j] < dlt) {
          dlt = mn[j];
          j1 = j;
        }
      }
      res += dlt;
      X[i] += dlt;
      for (int j : used) {
        X[P[j]] += dlt;
        Y[j] -= dlt;
      }
      last = dlt;
      j0 = j1;
    }
    do {
      P[j0] = P[way[j0]];
      j0 = way[j0];
    } while (j0 != 0);
  }
  std::vector<int> mat(N);
  for (int j = 1; j <= M; j++) {
    if (P[j] != 0) {
      mat[P[j] - 1] = j - 1;
    }
  }
  X.erase(X.begin());
  Y.erase(Y.begin());
  if constexpr (!MIN) {
    res = -res;
  }
  return {res, std::move(mat), std::move(X), std::move(Y)};
}
} // namespace noya

 // NOYA_HUNGARIAN_HPP