Hungarian algorithm: Difference between revisions

Revision as of 00:11, 28 March 2025

Latest revision as of 15:38, 23 May 2025

(8 intermediate revisions by 4 users not shown)

Line 148:

*/

#includeimport <cassert>;

#include <iostream>

import std;

#include <limits>

#include <vector>

template <typename T, typename U>

using namespace std;

using Pair = std::pair<T, U>;

template <typename T>

using Vector = std::vector<T>;

template <typename T>

using NumericLimits = std::numeric_limits<T>;

/**

* @brief Checks if b < a

*

* Sets a = min(a, b)

* @param a The first parameter to check

* @param b The second parameter to check

* @tparam The type to perform the check on

* @return true if b < a

*/

template <classtypename T>
constexpr bool ckmin(T & a, const T & b) {
return b < a ? a = b, 1true : 0false;
}

/**

* @brief Performs the Hungarian algorithm.

*

* Given J jobs and W workers (J <= W), computes the minimum cost to assign each

* prefix of jobs to distinct workers.

Line 172 ⟶ 188:

* to assign the first (j+1) jobs to distinct workers

*/

template <typename T>

template <class T> vectorVector<T> hungarian(const vectorVector<vectorVector<T>>& &C) {

const int J = (int)size(C), W = (int)size(C[0]);

const int J = static_cast<int>(C.size());

const int W = static_cast<int>(C[0].size());

assert(J <= W);

// job[w] = job assigned to w-th worker, or -1 if no job assigned

// note: a W-th worker was added for convenience

vectorVector<int> job(W + 1, -1);

vectorVector<T> ys(J),;
Vector<T> yt(W + 1); // potentials

// -yt[W] will equal the sum of all deltas

vectorVector<T> answers;

const T inf = numeric_limitsNumericLimits<T>::max();

for (int j_curjCur = 0; j_curjCur < J; ++j_curjCur) { // assign j_curjCur-th job

int w_curwCur = W;

job[w_curwCur] = j_curjCur;

// min reduced cost over edges from Z to worker w

vectorVector<T> min_tominTo(W + 1, inf);

vectorVector<int> prvprev(W + 1, -1); // previous worker on alternating path

vectorVector<bool> in_ZinZ(W + 1); // whether worker is in Z

while (job[w_curwCur] != -1) { // runs at most j_curjCur + 1 times

in_ZinZ[w_curwCur] = true;

const int j = job[w_curwCur];

T delta = inf;

int w_nextwNext;

for (int w = 0; w < W; ++w) {

if (!in_ZinZ[w]) {

if (ckmin(min_tominTo[w], C[j][w] - ys[j] - yt[w]))

prvprev[w] = w_curwCur;

if (ckmin(delta, min_tominTo[w])) w_next = w;

wNext = w;

}

}

// delta will always be nonnegative,

// except possibly during the first time this loop runs

// if any entries of C[j_curjCur] are negative

for (int w = 0; w <= W; ++w) {

if (in_ZinZ[w]) ys[job[w]] += delta, yt[w] -= delta;{

else min_to ys[job[w]] -+= delta;

yt[w] -= delta;

} else {

minTo[w] -= delta;

}

}

w_curwCur = w_nextwNext;

}

// update assignments along alternating path

for (int w; w_curwCur != W; w_curwCur = w) job[w_cur] = job[w = prv[w_cur]];

job[wCur] = job[w = prev[wCur]];

answers.push_back(-yt[W]);

}

Line 218 ⟶ 243:

/**

* @brief Performs a sanity check for the Hungarian algorithm.

*

* Sanity check: https://en.wikipedia.org/wiki/Hungarian_algorithm#Example

* First job (5):

Line 229 ⟶ 256:

* wash windows: Bob -> 3

*/

void sanity_check_hungariansanityCheckHungarian() {

vectorVector<vectorVector<int>> costs{{8, 5, 9}, {4, 2, 4}, {7, 3, 8}};

assert((hungarian(costs) == vectorVector<int>{5, 9, 15}));

cerr <<std::println(stderr, "Sanity check passed.\n");

}

/**

// * @brief solves https://open.kattis.com/problems/cordonbleu

void cordon_bleu() {

*/

int N, M;

void cordon_bleucordonBleu() {

cin >> N >> M;

vector<pair<int, int>> B(N), C(M);

int N, M;

vector<pair<int, int>> bottles(N), couriers(M);

for (auto &b std:: bottles) cin >> b.firstN >> b.secondM;

Vector<Pair<int, int>> B(N);

for (auto &c : couriers) cin >> c.first >> c.second;

pairVector<Pair<int, int>> restC(M);

vectorVector<pairPair<int, int>> bottles(N), couriers(M);

cin >> rest.first >> rest.second;

vectorVector<vectorPair<int>> costs(N, vector<int>> couriers(N + M - 1));

autofor dist(const = [&](pairPair<int, int>& x, pair<int, int>b: ybottles) {

returnstd::cin abs(x.first>> - yb.first) + abs(x.second ->> yb.second);

for (const Pair<int, int>& c: couriers)

for (auto &c : couriers) std::cin >> c.first >> c.second;

Pair<int, int> rest;

std::cin >> rest.first >> rest.second;

Vector<Vector<int>> costs(N, Vector<int>(N + M - 1));

auto dist = [&](const Pair<int, int>& x, const Pair<int, int>& y) -> int {

return std::abs(x.first - y.first) + std::abs(x.second - y.second);

};

for (int b = 0; b < N; ++b) {

for (int c = 0; c < M; ++c) { // courier -> bottle -> restaurant

costs[b][c] = dist(couriers[c], bottles[b]) + dist(bottles[b], rest);

for (int _ = 0; _ < N - 1; ++_) { // restaurant -> bottle -> restaurant

dist(couriers[c], bottles[b]) + dist(bottles[b], rest);

}

for (int _ = 0; _ < N - 1; ++_) { // restaurant -> bottle -> restaurant

costs[b][_ + M] = 2 * dist(bottles[b], rest);

}

}

cout << std::println(hungarian(costs).back() << "\n");

}

/**

* @brief Entry point into the program.

*

* @return The return code of the program.

*/

int main() {

sanity_check_hungariansanityCheckHungarian();

cordon_bleucordonBleu();

}

</syntaxhighlight>

Line 290 ⟶ 326:

<syntaxhighlight lang="cpp" line="1">

template <typename T>

template <class T> vectorVector<T> hungarian(const vectorVector<vectorVector<T>>& &C) {

const int J = (int)size(C), W = (int)size(C[0]);

const int J = (static_cast<int)size>(C), W = (int).size(C[0]));

const int W = static_cast<int>(C[0].size());

assert(J <= W);

// job[w] = job assigned to w-th worker, or -1 if no job assigned

// note: a W-th worker was added for convenience

vectorVector<int> job(W + 1, -1);

vectorVector<T> h(W); // Johnson potentials

vectorVector<T> answers;

T ans_curansCur = 0;

const T inf = numeric_limitsNumericLimits<T>::max();

// assign j_curjCur-th job using Dijkstra with potentials

for (int j_curjCur = 0; j_curjCur < J; ++j_curjCur) {

int w_curwCur = W; // unvisited worker with minimum distance

job[w_curwCur] = j_curjCur;

vectorVector<T> dist(W + 1, inf); // Johnson-reduced distances

dist[W] = 0;

vectorVector<bool> vis(W + 1); // whether visited yet

vectorVector<int> prvprev(W + 1, -1); // previous worker on shortest path

while (job[w_curwCur] != -1) { // Dijkstra step: pop min worker from heap

T min_distminDist = inf;

vis[w_curwCur] = true;

int w_nextwNext = -1; // next unvisited worker with minimum distance

// consider extending shortest path by w_curwCur -> job[w_curwCur] -> w

for (int w = 0; w < W; ++w) {

if (!vis[w]) {

// sum of reduced edge weights w_curwCur -> job[w_curwCur] -> w

T edge = C[job[w_curwCur]][w] - h[w];

if (w_curwCur != W) {

edge -= C[job[w_curwCur]][w_curwCur] - h[w_curwCur];

assert(edge >= 0); // consequence of Johnson potentials

}

if (ckmin(dist[w], dist[w_curwCur] + edge)) prv[w] = w_cur;

if (ckmin(min_dist, dist prev[w])) w_next = wwCur;

if (ckmin(minDist, dist[w]))

wNext = w;

}

}

w_curwCur = w_nextwNext;

}

for (int w = 0; w < W; ++w) { // update potentials

ckmin(dist[w], dist[w_curwCur]);

h[w] += dist[w];

}

ans_cur += h[w_curwCur];

for (int w; w_curwCur != W; w_curwCur = w) job[w_cur] = job[w = prv[w_cur]];

answers.push_back(ans_cur) job[wCur] = job[w = prev[wCur]];

answers.push_back(ansCur);

}

return answers;

Line 529 ⟶ 570:

If the number of starred zeros is {{mvar|n}} (or in the general case <math>min(n,m)</math>, where {{mvar|n}} is the number of people and {{mvar|m}} is the number of jobs), the algorithm terminates. See the [[#Result|Result]] subsection below on how to interpret the results.

Otherwise, find the lowest uncovered value. Subtract this from every unmarkeduncovered element and add it to every element covered by two lines. Go back to step 4.

This is equivalent to subtracting a number from all rows which are not covered and adding the same number to all columns which are covered. These operations do not change optimal assignments.