Flujo de costo mínimo

En esta sección voy a recorrer brevemente mi comprensión del flujo de costo mínimo. Recomiendo mucho leer Min Cost Flow de CP-Algorithms para entender primero la idea de la solución. Además, conviene revisar el tutorial de TopCoder para una explicación más detallada.

Desigualdad triangular

En teoría de grafos, la desigualdad triangular (triangle inequality) afirma que si hay una arista $u \rightarrow v$ de peso $w$ , y $d_k$ es el camino más corto hasta el nodo $k$ (para alguna definición razonable de camino más corto), entonces $d_v - d_u \leq w$ .

Recursos
Fuente	Recurso	Notas
Wikipedia	Triangle Inequality	Principalmente en geometría, pero transmite la idea.

Algoritmo de Johnson

La idea principal del algoritmo de Johnson es que si todos los pesos de arista son positivos, entonces ejecutar Dijkstra desde cada nodo daría un algoritmo $\mathcal{O}(VE \log E)$ . Si hay aristas negativas, el algoritmo de Johnson define una función de potencial $\pi$ tal que para toda arista $u,v,w$ se cumple $w>=\pi(u)-\pi(v)$ . Entonces cada peso de arista se puede transformar en $w \rightarrow w + \pi(v)-\pi(u)$ , lo que resulta en peso positivo. Esta condición coincide con la desigualdad triangular, así que podemos elegir un nodo arbitrariamente y ejecutar un algoritmo de camino más corto $\mathcal{O}(VE)$ para determinar esta función.

Recursos
Fuente	Recurso	Notas
GFG	Johnson's Algorithm for All-pairs shortest paths
Wikipedia	Johnson's Algorithm

Flujo de costo mínimo

La idea general del flujo de costo mínimo es empujar flujo repetidamente a lo largo del camino más corto. Como los grafos de flujo tienen aristas negativas, cada paso de forma naive tomaría $\mathcal{O}(VE)$ de tiempo. Para acelerarlo, podemos usar la misma función de potencial del algoritmo de Johnson para emplear Dijkstra en este proceso. En este caso debemos usar la distancia desde $S$ , el nodo fuente, como función $\pi$ . En cada paso, ejecutamos Dijkstra usando la función $\pi$ , actualizamos $\pi$ para que coincida con las distancias actuales, y luego empujamos flujo a lo largo del camino más corto, invirtiendo aristas según haga falta. Cuando se alcanza el flujo pedido o el sumidero es inalcanzable, terminamos.

Recursos
Fuente	Recurso	Notas
CP-Algorithms	Minimum-cost flow	Nota: no usa la solución óptima, pero explica bien el concepto.
TopCoder	Minimum Cost Flow Algorithms

Implementación

Dicho esto, aquí está mi implementación.


template <int MN, int MM>
struct MCF  // MN = nodes, MM = edges [assume edges one-directional]
{
  public:
	int N, M, S, T;
	int flow[MM * 2], cap[MM * 2], hd[MN], nx[MM * 2], to[MM * 2], cost[MM * 2];
	int pi[MN], p[MN], d[MN];
	int vis[MN];
	void init(int n, int s, int t) {
		N = n, S = s, T = t;
		memset(hd, -1, sizeof hd);
	}
	void adde1(int a, int b, int f, int c) {
		nx[M] = hd[a], hd[a] = M;
		to[M] = b, cost[M] = c, cap[M] = f;
		M++;
	}
	void adde(int a, int b, int f, int c) {
		adde1(a, b, f, c);
		adde1(b, a, 0, -c);
	}
	void setpi() {
		std::queue<int> q;
		memset(pi, 0x3e, sizeof pi);
		memset(vis, 0, sizeof vis);
		q.push(S);
		pi[S] = 0;
		for (int n; !q.empty();) {
			n = q.front();
			q.pop();
			for (int id = hd[n], x; ~id; id = nx[id]) {
				if (cap[id] - flow[id] <= 0) continue;
				x = to[id];
				if (ckmin(pi[x], pi[n] + cost[id])) assert(++vis[x] <= N), q.push(x);
			}
		}
	}
	struct state {
	  public:
		int n, d;
		bool operator>(state o) const { return d > o.d; }
	};
	void dijk() {
		std::priority_queue<state, std::vector<state>, std::greater<state>> q;
		memset(p, -1, N * sizeof p[0]);
		memset(vis, 0, N * sizeof vis[0]);
		memset(d, 0x3e, N * sizeof d[0]);
 
		d[S] = 0;
		q.push({S, 0});
		for (int n; !q.empty();) {
			n = q.top().n;
			q.pop();
			if (vis[n]) continue;
			vis[n] = 1;
			for (int id = hd[n], x, w; ~id; id = nx[id]) {
				if (cap[id] - flow[id] <= 0) continue;
				x = to[id];
				w = cost[id] + pi[n] - pi[x];
				if (ckmin(d[x], w + d[n])) p[x] = id, q.push({x, d[x]});
			}
		}
	}
	int mincost(int F) {
		setpi();
		int C = 0;
		while (F > 0) {
			dijk();
			if (d[T] == INF) return INF;
			int c = d[T] + pi[T] - pi[S], f = F;
			for (int x = T; x != S; x = to[p[x] ^ 1]) ckmin(f, cap[p[x]] - flow[p[x]]);
			C += c * f;
			for (int x = T; x != S; x = to[p[x] ^ 1]) {
				flow[p[x]] += f;
				flow[p[x] ^ 1] -= f;
			}
			F -= f;
			for (int i = 0; i < N; ++i) pi[i] += d[i];
		}
		return C;
	}
};

También conviene ver la implementación de Benq y la implementación de KACTL (que son bastante mejores que la mía).

Problemas

Fuente	Nombre	Dificultad	Tags	Solución
CF	Build String	Normal	MCF	—
CF	Four Melodies	Normal	MCF	—
CF	April Fools' Problem (medium)	Normal		—
CF	Team Building	Difícil	Bitmasks, MCF	Solución
CF	For the Emperor!	Difícil		—
CF	Cow and Exercise	Muy difícil		—
CF	One Billion Shades of Grey	Muy difícil		Solución

Aplicaciones

Problema de asignación

Hecho	Fuente	Nombre	Dificultad	Tags	Solución
	CSES	★ Task Assignment	Fácil		en el módulo

El problema de asignación se entiende mejor como un problema de flujo máximo de costo mínimo. En nuestra formulación, asignaremos $W$ trabajadores a $J$ trabajos, $J \ge W$ , y el costo de asignar el $i$ -ésimo trabajo al $j$ -ésimo trabajador es $C_{i,j}$ .

Crearemos una red de flujo con un nodo fuente $S$ , nodos de trabajo $J_i$ , nodos de trabajador $W_j$ y un nodo sumidero $E$ .

Tendremos las siguientes aristas:

$(S, J_i)$ con $(capacity, cost) = (1, 0)$
$(J_i, W_j)$ con $(1, C_{i,j})$
$(W_j, E)$ con $(1, 0)$

La respuesta será el costo mínimo del flujo máximo del grafo. Un ejemplo:

Hungarian Algorithm Flows Diagram

Para resolverlo, usaremos el algoritmo de flujo máximo de costo mínimo detallado arriba. En cada iteración del algoritmo, asignaremos un trabajo a un trabajador. En cada iteración, primero asignamos el trabajo a un trabajador auxiliar y luego intentamos hallar un camino aumentante de costo mínimo. Para optimizar esto, solo ejecutaremos Dijkstra sobre los nodos de trabajadores. Para los detalles, consultar el código.

Como hay $J$ iteraciones (una por cada trabajo) y cada iteración toma $O(W^2)$ de tiempo, la complejidad temporal total es $O(J W^2)$ .

Implementación

De Wikipedia (que a su vez copia de e-maxx):


#include <bits/stdc++.h>
using namespace std;
 
int ckmin(int &a, int b) { return a > b ? ((a = b), true) : false; }
 
/**
 * @return the jobs of each worker in the optimal assignment,
 * or -1 if the worker is not assigned
 */
template <class T> vector<int> hungarian(const vector<vector<T>> &C) {
	int J = C.size();
	int W = C[0].size();
	assert(J <= W);
 
	// job[w] = trabajo asignado al w-ésimo trabajador, o -1 si no hay trabajo
	// nota: se agregó un W-ésimo trabajador por conveniencia
	vector<int> job(W + 1, -1);
	vector<T> h(W);  // potenciales de Johnson
 
	const T inf = numeric_limits<T>::max();
	// asignar el j_cur-ésimo trabajo usando Dijkstra con potenciales
	for (int j_cur = 0; j_cur < J; j_cur++) {
		int w_cur = W;  // trabajador no visitado con distancia mínima
		job[w_cur] = j_cur;
 
		vector<T> dist(W + 1, inf);  // distancias reducidas de Johnson
		dist[W] = 0;
		vector<bool> vis(W + 1);     // si ya se visitó
		vector<int> prv(W + 1, -1);  // trabajador anterior en el camino más corto
		while (job[w_cur] != -1) {   // paso de Dijkstra: extraer el trabajador mín.
			T min_dist = inf;
			vis[w_cur] = true;
			int w_next = -1;  // siguiente trabajador no visitado con dist. mínima
 
			// considerar extender el camino más corto por w_cur -> job[w_cur] -> w
			for (int w = 0; w < W; w++) {
				if (!vis[w]) {
					// suma de pesos reducidos w_cur -> job[w_cur] -> w
					T edge = C[job[w_cur]][w] - h[w];
					if (w_cur != W) {
						edge -= C[job[w_cur]][w_cur] - h[w_cur];
						assert(edge >= 0);
					}
					if (ckmin(dist[w], dist[w_cur] + edge)) { prv[w] = w_cur; }
					if (ckmin(min_dist, dist[w])) { w_next = w; }
				}
			}
			w_cur = w_next;
		}
 
		for (int w = 0; w < W; w++) {  // actualizar potenciales
			ckmin(dist[w], dist[w_cur]);
			h[w] += dist[w];
		}
 
		while (w_cur != W) {  // actualizar la asignación de trabajos
			job[w_cur] = job[prv[w_cur]];
			w_cur = prv[w_cur];
		}
	}
 
	return job;
}
 
int main() {
	int n;
	cin >> n;
	vector<vector<int>> c(n, vector<int>(n));
 
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++) { cin >> c[i][j]; }
	}
 
	vector<int> mat = hungarian(c);
	int cost = 0;
	for (int i = 0; i < n; i++) { cost += c[mat[i]][i]; }
	cout << cost << endl;
	for (int i = 0; i < n; i++) { cout << mat[i] + 1 << ' ' << i + 1 << endl; }
}

Recursos
Fuente	Recurso	Notas
Wikipedia	Hungarian Algorithm
YouTube - Algorithms Thread	Hungarian for non Hungarians
Topcoder	Assignment Problem and Hungarian Algorithm

Hecho	Fuente	Nombre	Dificultad	Tags	Solución
	CF	Vasya and Endless Credits	Fácil		—
	Kattis	Cordon Bleu	Fácil		—