Content deleted Content added
mark last step of matrix interpretation as unclear |
add O(n^3) Hungarian algorithm + code |
||
Line 105:
===Proof that {{mvar|y}} remains a potential===
To show that {{mvar|y}} remains a potential after being adjusted, it suffices to show that no edge has its total potential increased beyond its cost. This is already established for edges in {{mvar|M}} by the preceding paragraph, so consider an arbitrary edge {{mvar|uv}} from {{mvar|S}} to {{mvar|T}}. If <math>y(u)</math> is increased by {{math|Δ}}, then either <math>v \in Z \cap T</math>, in which case <math>y(v)</math> is decreased by {{math|Δ}}, leaving the total potential of the edge unchanged, or <math>v \in T \setminus Z</math>, in which case the definition of {{math|Δ}} guarantees that <math>y(u)+y(v)+\Delta \leq c(u,v)</math>. Thus {{mvar|y}} remains a potential.
==The algorithm in <math>O(n^3)</math> time==
Suppose there are <math>J</math> jobs and <math>W</math> workers (<math>J\le W</math>). We describe how to compute for each prefix of jobs the minimum total cost to assign each of these jobs to distinct workers. Specifically, we add the <math>j</math>th job and update the total cost in time <math>O(jW)</math>, yielding an overall time complexity of <math>O\left(\sum_{j=1}^Jj\cdot W\right)=O(J^2W)</math>. Note that this is much better than <math>O(W^3)</math> when the number of jobs is small relative to the number of workers.
===Adding the <math>j</math>th job in <math>O(jW)</math> time===
We use the same notation to the previous section, though we modify their definitions as necessary. Let <math>S_j</math> denote the set of the first <math>j</math> jobs and <math>T</math> denote the set of all workers. Before the <math>j</math>th step of the algorithm, we assume that we have a matching that matches all jobs in <math>S_{j-1}</math> and potentials such that the potentials of all unmatched workers are equal and at least the maximum potential over all matched workers. During the <math>j</math>th step, we add the <math>j</math>th job to <math>S_{j-1}</math> to form <math>S_j</math> and initialize <math>Z=\{j\}</math>. At all times, every vertex in <math>Z</math> will be reachable from the <math>j</math>th job in <math>G_y</math>.
While <math>Z</math> does not contain a worker that has not been assigned a job, let
:<math>\Delta := \min\{c(j, w)-y(j)-y(w):j\in Z\cap S_j, w\in T\setminus Z\}</math>
and <math>w_{\text{next}}</math> denote any <math>w</math> at which the minimum is attained. After adjusting the potentials in the way described in the previous section, there is now a tight edge from <math>Z</math> to <math>w_{\text{next}}</math>.
* If <math>w_{\text{next}}</math> is unmatched, then we have an [[Matching (graph theory)|augmenting path]] in the subgraph of tight edges from <math>j</math> to <math>w_{\text{next}}</math>. After toggling the matching along this path, we have now matched the first <math>j</math> jobs, and this procedure terminates.
* Otherwise, we add <math>w_{\text{next} }</math> and the job matched with it to <math>Z</math>.
Recomputing <math>\Delta</math> and <math>w_{\text{next}}</math> after adjusting potentials can be done in <math>O(W)</math> time. Case 1 can occur at most <math>j-1</math> times before case 2 occurs and the procedure terminates.
===Implementation in C++===
For convenience of implementation, the code below adds an additional worker <math>w_{W}</math> such that <math>y(w_{W})</math> stores the negation of the sum of all <math>\Delta</math> computed so far. After the <math>j</math>th job is added and the matching updated, the cost of the current matching equals the sum of all <math>\Delta</math> computed so far, or <math>-y(w_{W})</math>.<syntaxhighlight lang="cpp" line="1">
/**
* Solution to https://open.kattis.com/problems/cordonbleu using Hungarian
* algorithm.
* Based on https://github.com/mpfeifer1/Kattis/blob/master/cordonbleu.cpp
*/
#include <cassert>
#include <iostream>
#include <limits>
#include <vector>
using namespace std;
/**
* Sets a = min(a, b)
* @return true if b < a
*/
template <class T> bool ckmin(T &a, const T &b) { return b < a ? a = b, 1 : 0; }
/**
* Given J jobs and W workers (J <= W), computes the minimum cost to assign each
* prefix of jobs to distinct workers.
*
* @tparam T a type large enough to represent integers on the order of J *
* max(|C|)
* @param C a matrix of dimensions JxW such that C[j][w] = cost to assign j-th
* job to w-th worker (possibly negative)
*
* @return a vector of length J, with the j-th entry equaling the minimum cost
* to assign the first (j+1) jobs to distinct workers
*/
template <class T> vector<T> hungarian(const vector<vector<T>> &C) {
const int J = (int)size(C), W = (int)size(C[0]);
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
vector<int> job(W + 1, -1);
vector<T> ys(J), yt(W + 1); // potentials
// -yt[W] will equal the sum of all deltas
vector<T> answers;
const T inf = numeric_limits<T>::max();
for (int j_cur = 0; j_cur < J; ++j_cur) { // assign j_cur-th job
int w_cur = W;
job[w_cur] = j_cur;
// minimum reduced cost over edges entering worker w
vector<T> dist(W + 1, inf);
vector<int> prv(W + 1, -1); // for restoring alternating path
vector<bool> vis(W + 1); // whether worker is in Z
while (job[w_cur] != -1) { // runs at most j_cur + 1 times
vis[w_cur] = true;
const int j = job[w_cur];
T delta = inf;
int w_next;
for (int w = 0; w < W; ++w) {
if (!vis[w]) {
if (ckmin(dist[w], C[j][w] - ys[j] - yt[w])) prv[w] = w_cur;
if (ckmin(delta, dist[w])) w_next = w;
}
}
// delta will always be non-negative,
// except during the first iteration if one of C[j_cur][*] is
// negative
for (int w = 0; w <= W; ++w) {
if (vis[w]) ys[job[w]] += delta, yt[w] -= delta;
else dist[w] -= delta;
}
w_cur = w_next;
}
// update assignments along alternating path
for (int w; w_cur != -1; w_cur = w) job[w_cur] = job[w = prv[w_cur]];
answers.push_back(-yt[W]);
}
return answers;
}
/**
* Sanity check: https://en.wikipedia.org/wiki/Hungarian_algorithm#Example
* First job (5):
* clean bathroom: Bob -> 5
* First + second jobs (9):
* clean bathroom: Bob -> 5
* sweep floors: Alice -> 4
* First + second + third jobs (15):
* clean bathroom: Alice -> 8
* sweep floors: Dora -> 4
* wash windows: Bob -> 3
*/
void sanity_check_hungarian() {
vector<vector<int>> costs{{8, 5, 9}, {4, 2, 4}, {7, 3, 8}};
assert((hungarian(costs) == vector<int>{5, 9, 15}));
cerr << "Sanity check passed.\n";
}
// solves https://open.kattis.com/problems/cordonbleu
void cordon_bleu() {
int N, M;
cin >> N >> M;
vector<pair<int, int>> B(N), C(M);
vector<pair<int, int>> bottles(N), couriers(M);
for (auto &b : bottles) cin >> b.first >> b.second;
for (auto &c : couriers) cin >> c.first >> c.second;
pair<int, int> rest;
cin >> rest.first >> rest.second;
vector<vector<int>> costs(N, vector<int>(N + M - 1));
auto dist = [&](pair<int, int> x, pair<int, int> y) {
return abs(x.first - y.first) + 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
costs[b][_ + M] = 2 * dist(bottles[b], rest);
}
}
cout << hungarian(costs).back() << "\n";
}
int main() {
sanity_check_hungarian();
cordon_bleu();
}
</syntaxhighlight>
==Matrix interpretation==
| |||