Задача о назначениях. Решение с помощью min-cost-flow

Задача имеет две эквивалентные постановки:

Дана квадратная матрица A[1..N,1..N]. Нужно выбрать в ней N элементов так, чтобы в каждой строке и столбце был выбран ровно один элемент, а сумма значений этих элементов была наименьшей.
Имеется N заказов и N станков. Про каждый заказ известна стоимость его изготовления на каждом станке. На каждом станке можно выполнять только один заказ. Требуется распределить все заказы по станкам так, чтобы минимизировать суммарную стоимость.

Здесь мы рассмотрим решение задачи на основе алгоритма нахождения потока минимальной стоимости (min-cost-flow), решив задачу о назначениях за O (N⁵).

Описание

Построим двудольную сеть: имеется исток S, сток T, в первой доле находятся N вершин (соответствующие строкам матрицы или заказам), во второй - тоже N вершин (соответствующие столбцам матрицы или станкам). Между каждой вершиной i первой доли и каждой вершиной j второй доли проведём ребро с пропускной способностью 1 и стоимостью A_ij. От истока S проведём рёбра ко всем вершинам i первой доли с пропускной способностью 1 и стоимостью 0. От каждой вершины второй доли j к стоку T проведём ребро с пропускной способностью 1 и стоимостью 0.

Найдём в полученной сети максимальный поток минимальной стоимости. Очевидно, величина потока будет равна N. Далее, очевидно, что для каждой вершины i из первой доли найдётся ровно одна вершина j из второй доли, такая, что поток F_ij = 1. Наконец, очевидно, это взаимно однозначное соответствие между вершинами первой доли и вершинами второй доли является решением задачи (поскольку найденный поток имеет минимальную стоимость, то сумма стоимостей выбранных рёбер будет наименьшей из возможных, что и является критерием оптимальности).

Асимптотика этого решения задачи о назначениях зависит от того, каким алгоритмом производится поиск максимального потока минимальной стоимости. Асимптотика составит O (N³) при использовании алгоритма Дейкстры или O (N⁴) при использовании алгоритма Форда-Беллмана.

Реализация

Приведённая здесь реализация длинноватая, возможно, её можно значительно сократить.

typedef vector<int> vint;
typedef vector<vint> vvint;

const int INF = 1000*1000*1000;


int main()
{
	int n;
	vvint a (n, vint (n));
	... чтение a ...

	int m = n * 2 + 2;
	vvint f (m, vint (m));
	int s = m-2, t = m-1;
	int cost = 0;
	for (;;)
	{
		vector<int> dist (m, INF);
		vector<int> p (m);
		vector<int> type (m, 2);
		deque<int> q;
		dist[s] = 0;
		p[s] = -1;
		type[s] = 1;
		q.push_back (s);
		for (; !q.empty(); )
		{
			int v = q.front(); q.pop_front();
			type[v] = 0;
			if (v == s)
			{
				for (int i=0; i<n; ++i)
					if (f[s][i] == 0)
					{
						dist[i] = 0;
						p[i] = s;
						type[i] = 1;
						q.push_back (i);
					}
			}
			else
			{
				if (v < n)
				{
					for (int j=n; j<n+n; ++j)
						if (f[v][j] < 1 && dist[j] > dist[v] + a[v][j-n])
						{
							dist[j] = dist[v] + a[v][j-n];
							p[j] = v;
							if (type[j] == 0)
								q.push_front (j);
							else if (type[j] == 2)
								q.push_back (j);
							type[j] = 1;
						}
				}
				else
				{
					for (int j=0; j<n; ++j)
						if (f[v][j] < 0 && dist[j] > dist[v] - a[j][v-n])
						{
							dist[j] = dist[v] - a[j][v-n];
							p[j] = v;
							if (type[j] == 0)
								q.push_front (j);
							else if (type[j] == 2)
								q.push_back (j);
							type[j] = 1;
						}
				}
			}
		}

		int curcost = INF;
		for (int i=n; i<n+n; ++i)
			if (f[i][t] == 0 && dist[i] < curcost)
			{
				curcost = dist[i];
				p[t] = i;
			}
		if (curcost == INF) break;
		cost += curcost;
		for (int cur=t; cur!=-1; cur=p[cur])
		{
			int prev = p[cur];
			if (prev!=-1)
				f[cur][prev] = - (f[prev][cur] = 1);
		}

	}

	printf ("%d\n", cost);
	for (int i=0; i<n; ++i)
		for (int j=0; j<n; ++j)
			if (f[i][j+n] == 1)
				printf ("%d ", j+1);

}