Hungarian algorithm: Difference between revisions

Revision as of 15:09, 20 April 2025

Latest revision as of 15:38, 23 May 2025

(6 intermediate revisions by 4 users not shown)

Line 148:

*/

#includeimport <cassert>;

import std;

template <typename T, typename U>

using Pair = std::pair;

using VectorPair = std::vectorpair<T, U>;

template <typename T>

using PairVector = std::pairvector<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 <typename 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 188 ⟶ 200:

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

Vector<T> answers;

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

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

int wCur = W;

Line 194 ⟶ 206:

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

Vector<T> minTo(W + 1, inf);

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

Vector<bool> inZ(W + 1); // whether worker is in Z

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

in_ZinZ[wCur] = true;

const int j = job[wCur];

T delta = inf;

Line 204 ⟶ 216:

if (!inZ[w]) {

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

prvprev[w] = wCur;

if (ckmin(delta, minTo[w])) wNext = w;

wNext = w;

}

}

Line 223 ⟶ 236:

// update assignments along alternating path

for (int w; wCur != W; wCur = w)

job[wCur] = job[w = prvprev[wCur]];

answers.push_back(-yt[W]);

}

Line 230 ⟶ 243:

/**

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

*

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

* First job (5):

Line 243 ⟶ 258:

void sanityCheckHungarian() {

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

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

std::println(stderr, "Sanity check passed.");

}

/**

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

*/

void cordonBleu() {

int N;

Line 256 ⟶ 273:

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

Vector<Pair<int, int>> couriers(M);

for (const Pair<int, int>& b : bottles)

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

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

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

Pair<int, int> rest;

Line 275 ⟶ 292:

}

/**

* @brief Entry point into the program.

*

* @return The return code of the program.

*/

int main() {

sanityCheckHungarian();

Line 304 ⟶ 326:

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

template <typename T>

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

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 543 ⟶ 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.