Hungarian.H

Go to the documentation of this file.

/*
                          Aleph_w
 
  Data structures & Algorithms
  version 2.0.0b
  https://github.com/lrleon/Aleph-w
 
  This file is part of Aleph-w library
 
  Copyright (c) 2002-2026 Leandro Rabindranath Leon
 
  Permission is hereby granted, free of charge, to any person obtaining a copy
  of this software and associated documentation files (the "Software"), to deal
  in the Software without restriction, including without limitation the rights
  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  copies of the Software, and to permit persons to whom the Software is
  furnished to do so, subject to the following conditions:
 
  The above copyright notice and this permission notice shall be included in all
  copies or substantial portions of the Software.
 
  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  SOFTWARE.
*/
 
#ifndef HUNGARIAN_H
#define HUNGARIAN_H
 
#include <cmath>
#include <limits>
#include <type_traits>
#include <utility>
#include <initializer_list>
#include <tpl_dynMat.H>
#include <tpl_array.H>
#include <htlist.H>
#include <ah-errors.H>
#include <ahFunction.H>
 
namespace Aleph {
template <typename Cost_Type = double>
struct Hungarian_Result
{
  static_assert(std::is_signed_v<std::remove_cv_t<Cost_Type>> or
                  std::is_floating_point_v<std::remove_cv_t<Cost_Type>>,
                "Cost_Type must be signed or floating-point to avoid "
                "underflow when negating costs");
  Cost_Type total_cost = Cost_Type{0};  
  Array<long> row_to_col;               
  Array<long> col_to_row;               
  size_t orig_rows = 0;                 
  size_t orig_cols = 0;                 
 
  [[nodiscard]] DynList<std::pair<size_t, size_t>> get_pairs() const
  {
    DynList<std::pair<size_t, size_t>> pairs;
    for (size_t i = 0; i < orig_rows; ++i)
      if (const long j = row_to_col[i]; j >= 0 and static_cast<size_t>(j) < orig_cols)
        pairs.append(std::make_pair(i, static_cast<size_t>(j)));
    return pairs;
  }
};
 
template <typename Cost_Type = double>
class Hungarian_Assignment
 
{
  static_assert(std::is_signed_v<std::remove_cv_t<Cost_Type>> or
                  std::is_floating_point_v<std::remove_cv_t<Cost_Type>>,
                "Cost_Type must be signed or floating-point to avoid "
                "underflow when negating costs");
  size_t n_ = 0;          // Padded square dimension
  size_t orig_rows_ = 0;  // Original row count
  size_t orig_cols_ = 0;  // Original column count
  Cost_Type total_cost_ = Cost_Type{0};
  Array<long> row_to_col_;  // row i -> column (or -1 if dummy)
  Array<long> col_to_row_;  // column j -> row (or -1 if dummy)
 
  void solve(const DynMatrix<Cost_Type> &cost)
  {
    const long n = static_cast<long>(n_);
    const auto sz = static_cast<size_t>(n + 1);
 
    // Dual variables (potentials), 1-indexed
    Array<Cost_Type> u(sz, Cost_Type{0});
    Array<Cost_Type> v(sz, Cost_Type{0});
 
    // p[j] = row matched to column j (0 = unmatched)
    Array<long> p(sz, 0L);
 
    const Cost_Type Inf = std::numeric_limits<Cost_Type>::max() / 2;
 
    // For each row, find the shortest augmenting path
    for (long i = 1; i <= n; ++i)
      {
        // "Virtual" column 0 is matched to row i
        p(0) = i;
 
        // Minimum reduced cost to reach each column
        Array<Cost_Type> dist(sz, Inf);
        dist(0) = Cost_Type{0};
 
        // way[j] = previous column on the shortest path to j
        Array<long> way(sz, 0L);
 
        // visited[j] = true if column j is in the "tree"
        Array<bool> visited(sz, false);
 
        // Dijkstra-like scan
        long j0 = 0;  // current column (start at virtual column 0)
        do
          {
            visited(j0) = true;
            const long i0 = p(j0);  // row matched to current column
            Cost_Type delta = Inf;
            long j1 = -1;
 
            for (long j = 1; j <= n; ++j)
              {
                if (visited(j))
                  continue;
 
                // Reduced cost: c[i0][j] - u[i0] - v[j]
                const Cost_Type reduced =
                  cost.read(static_cast<size_t>(i0 - 1), static_cast<size_t>(j - 1)) - u(i0) - v(j);
 
                if (reduced < dist(j))
                  {
                    dist(j) = reduced;
                    way(j) = j0;
                  }
 
                if (dist(j) < delta)
                  {
                    delta = dist(j);
                    j1 = j;
                  }
              }
 
            // Update potentials
            for (long j = 0; j <= n; ++j)
              if (visited(j))
                {
                  u(p(j)) += delta;
                  v(j) -= delta;
                }
              else
                dist(j) -= delta;
 
            j0 = j1;
        } while (p(j0) != 0);  // until we reach a free column
 
        // Augment along the path
        do
          {
            const long j1 = way(j0);
            p(j0) = p(j1);
            j0 = j1;
        } while (j0 != 0);
      }
 
    // Extract results
    total_cost_ = -v(0);  // v[0] accumulates the total cost
 
    row_to_col_ = Array<long>(orig_rows_, -1L);
    col_to_row_ = Array<long>(orig_cols_, -1L);
 
    for (long j = 1; j <= n; ++j)
      {
        const long row = p(j) - 1;  // convert to 0-based
        const long col = j - 1;
 
        if (row >= 0 and static_cast<size_t>(row) < orig_rows_ and col >= 0 and
            static_cast<size_t>(col) < orig_cols_)
          {
            row_to_col_(static_cast<size_t>(row)) = col;
            col_to_row_(static_cast<size_t>(col)) = row;
          }
      }
  }
 
public:
  explicit Hungarian_Assignment(const DynMatrix<Cost_Type> &cost)
  {
    orig_rows_ = cost.rows();
    orig_cols_ = cost.cols();
    ah_invalid_argument_if(orig_rows_ == 0 or orig_cols_ == 0)
      << "Hungarian_Assignment: cost matrix must not be empty";
 
    if constexpr (std::is_floating_point_v<Cost_Type>)
      for (size_t i = 0; i < orig_rows_; ++i)
        for (size_t j = 0; j < orig_cols_; ++j)
          ah_invalid_argument_if(not std::isfinite(cost.read(i, j)))
            << "Hungarian_Assignment: cost[" << i << "][" << j << "] is not finite";
 
    n_ = orig_rows_ > orig_cols_ ? orig_rows_ : orig_cols_;
 
    // Build padded square matrix (zeros for padding)
    DynMatrix<Cost_Type> padded(n_, n_, Cost_Type{0});
    padded.allocate();
    for (size_t i = 0; i < orig_rows_; ++i)
      for (size_t j = 0; j < orig_cols_; ++j)
        padded(i, j) = cost.read(i, j);
 
    solve(padded);
  }
 
  Hungarian_Assignment(std::initializer_list<std::initializer_list<Cost_Type>> rows)
  {
    ah_invalid_argument_if(rows.size() == 0)
      << "Hungarian_Assignment: cost matrix must not be empty";
 
    orig_rows_ = rows.size();
    orig_cols_ = rows.begin()->size();
    ah_invalid_argument_if(orig_cols_ == 0)
      << "Hungarian_Assignment: cost matrix must not have zero columns";
    n_ = orig_rows_ > orig_cols_ ? orig_rows_ : orig_cols_;
 
    DynMatrix<Cost_Type> padded(n_, n_, Cost_Type{0});
    padded.allocate();
 
    size_t i = 0;
    for (const auto &row : rows)
      {
        ah_invalid_argument_if(row.size() != orig_cols_)
          << "Hungarian_Assignment: all rows must have the same length";
        size_t j = 0;
        for (const auto &val : row)
          {
            if constexpr (std::is_floating_point_v<Cost_Type>)
              ah_invalid_argument_if(not std::isfinite(val))
                << "Hungarian_Assignment: non-finite cost value at row " << i;
            padded(i, j++) = val;
          }
        ++i;
      }
 
    solve(padded);
  }
 
  [[nodiscard]] Cost_Type get_total_cost() const noexcept
  {
    return total_cost_;
  }
 
  [[nodiscard]] long get_assignment(const size_t row) const
  {
    ah_out_of_range_error_if(row >= orig_rows_) << "Hungarian_Assignment::get_assignment: row "
                                                << row << " out of range [0, " << orig_rows_ << ")";
    return row_to_col_[row];
  }
 
  [[nodiscard]] DynList<std::pair<size_t, size_t>> get_assignments() const
  {
    DynList<std::pair<size_t, size_t>> pairs;
    for (size_t i = 0; i < orig_rows_; ++i)
      if (const long j = row_to_col_[i]; j >= 0 and static_cast<size_t>(j) < orig_cols_)
        pairs.append(std::make_pair(i, static_cast<size_t>(j)));
    return pairs;
  }
 
  [[nodiscard]] const Array<long> &get_row_assignments() const noexcept
  {
    return row_to_col_;
  }
 
  [[nodiscard]] const Array<long> &get_col_assignments() const noexcept
  {
    return col_to_row_;
  }
 
  [[nodiscard]] Array<long> extract_row_assignments() && noexcept
  {
    return std::move(row_to_col_);
  }
 
  [[nodiscard]] Array<long> extract_col_assignments() && noexcept
  {
    return std::move(col_to_row_);
  }
 
  [[nodiscard]] size_t dimension() const noexcept
  {
    return n_;
  }
 
  [[nodiscard]] size_t rows() const noexcept
  {
    return orig_rows_;
  }
 
  [[nodiscard]] size_t cols() const noexcept
  {
    return orig_cols_;
  }
};
 
template <typename Cost_Type>
[[nodiscard]] Hungarian_Result<Cost_Type> hungarian_assignment(const DynMatrix<Cost_Type> &cost)
{
  static_assert(std::is_signed_v<std::remove_cv_t<Cost_Type>> or
                  std::is_floating_point_v<std::remove_cv_t<Cost_Type>>,
                "Cost_Type must be signed or floating-point to avoid "
                "underflow when negating costs");
 
  Hungarian_Assignment<Cost_Type> ha(cost);
  Hungarian_Result<Cost_Type> result;
  // Read scalar fields before moving from ha.
  result.total_cost = ha.get_total_cost();
  result.orig_rows = ha.rows();
  result.orig_cols = ha.cols();
  result.row_to_col = std::move(ha).extract_row_assignments();
  result.col_to_row = std::move(ha).extract_col_assignments();
  return result;
}
 
template <typename Cost_Type>
[[nodiscard]] Hungarian_Result<Cost_Type> hungarian_max_assignment(const DynMatrix<Cost_Type> &cost)
{
  static_assert(std::is_signed_v<std::remove_cv_t<Cost_Type>> or
                  std::is_floating_point_v<std::remove_cv_t<Cost_Type>>,
                "Cost_Type must be signed or floating-point to avoid "
                "underflow when negating costs");
 
  // Use a promoted signed type for negation to avoid overflow when
  // Cost_Type is a signed integer and a cell contains its minimum value.
  using Common = std::common_type_t<Cost_Type, long long>;
  using Promoted =
    std::conditional_t<std::is_floating_point_v<Common>, Common, std::make_signed_t<Common>>;
 
  const size_t rows = cost.rows();
  const size_t cols = cost.cols();
  DynMatrix<Promoted> negated(rows, cols);
  negated.allocate();
  for (size_t i = 0; i < rows; ++i)
    for (size_t j = 0; j < cols; ++j)
      {
        if constexpr (std::is_integral_v<Promoted> and std::is_signed_v<Promoted>)
          {
            auto v = static_cast<Promoted>(cost.read(i, j));
            ah_overflow_error_if(v == std::numeric_limits<Promoted>::min())
              << "Cannot negate minimum integer value";
            negated(i, j) = -v;
          }
        else
          negated(i, j) = -static_cast<Promoted>(cost.read(i, j));
      }
 
  auto inner = hungarian_assignment(negated);
 
  Hungarian_Result<Cost_Type> result;
  result.total_cost = static_cast<Cost_Type>(-inner.total_cost);
  result.row_to_col = std::move(inner.row_to_col);
  result.col_to_row = std::move(inner.col_to_row);
  result.orig_rows = inner.orig_rows;
  result.orig_cols = inner.orig_cols;
  return result;
}
}  // end namespace Aleph
 
#endif  // HUNGARIAN_H