mcmf4.cpp

/* Min cost max flow (Edmonds-Karp relabelling + fast heap Dijkstra)
 * Based on code by Frank Chu and Igor Naverniouk
 * (http://shygypsy.com/tools/mcmf4.cpp)
 *
 * COMPLEXITY:
 *      - Worst case: O(min(m*log(m)*flow, n*m*log(m)*fcost))
 * FIELD TESTING:
 *      - Valladolid 10594: Data Flow
 * REFERENCE:
 *      Edmonds, J., Karp, R.  "Theoretical Improvements in Algorithmic
 *          Efficieincy for Network Flow Problems".
 *      This is a slight improvement of Frank Chu's implementation.
 **/

#define Inf (LLONG_MAX/2)
#define BUBL { \
    t = q[i]; q[i] = q[j]; q[j] = t; \
    t = inq[q[i]]; inq[q[i]] = inq[q[j]]; inq[q[j]] = t; }
#define Pot(u,v) (d[u] + pi[u] - pi[v])
struct MinCostMaxFlow {
	typedef long long LL;
	int n, qs;
	vector<vector<LL> > cap, cost, fnet;
	vector<vector<int> > adj;
	vector<LL> d, pi;
	vector<int> deg, par, q, inq;
	
	// n = number of vertices
	MinCostMaxFlow(int n): n(n), qs(0), deg(n+1), par(n+1), d(n+1), q(n+1), inq(n+1), pi(n+1), cap(n+1, vector<LL>(n+1)), cost(cap), fnet(cap), adj(n+1, vector<int>(n+1)) {}
	
	// call to add a directed edge. vertices are 0-based
	// ALL COSTS MUST BE NON-NEGATIVE
	void AddEdge(int from, int to, LL cap_, LL cost_) {
		cap[from][to] = cap_; cost[from][to] = cost_;
	}
	
	bool dijkstra( int s, int t ) {
		fill(d.begin(), d.end(), 0x3f3f3f3f3f3f3f3fLL);
		fill(par.begin(), par.end(), -1);
		fill(inq.begin(), inq.end(), -1);
		d[s] = qs = 0;
		inq[q[qs++] = s] = 0;
		par[s] = n;
		while( qs )	{
			int u = q[0]; inq[u] = -1;
			q[0] = q[--qs];
			if( qs ) inq[q[0]] = 0;
			for( int i = 0, j = 2*i + 1, t; j < qs; i = j, j = 2*i + 1 ) {
				if( j + 1 < qs && d[q[j + 1]] < d[q[j]] ) j++;
				if( d[q[j]] >= d[q[i]] ) break;
				BUBL;
			}
			for( int k = 0, v = adj[u][k]; k < deg[u]; v = adj[u][++k] ) {
				if( fnet[v][u] && d[v] > Pot(u,v) - cost[v][u] ) 
					d[v] = Pot(u,v) - cost[v][par[v] = u];			
				if( fnet[u][v] < cap[u][v] && d[v] > Pot(u,v) + cost[u][v] ) 
					d[v] = Pot(u,v) + cost[par[v] = u][v];
				if( par[v] == u ) {
					if( inq[v] < 0 ) { inq[q[qs] = v] = qs; qs++; }
					for( int i = inq[v], j = ( i - 1 )/2, t;
						 d[q[i]] < d[q[j]]; i = j, j = ( i - 1 )/2 )
						 BUBL;
				}
			}
		}
		for( int i = 0; i < n; i++ ) if( pi[i] < Inf ) pi[i] += d[i];
		return par[t] >= 0;
	}

	// Returns: (flow, total cost) between source s and sink t
	// Call this once only. fnet[i][j] contains the flow from i to j. Careful, fnet[i][j] and fnet[j][i] could both be positive.
	pair<LL, LL> mcmf4(int s, int t) {
		for( int i = 0; i < n; i++ )
			for( int j = 0; j < n; j++ )
				if( cap[i][j] || cap[j][i] ) adj[i][deg[i]++] = j;
		LL flow = 0; LL fcost = 0;
		while( dijkstra( s, t ) ) {
			LL bot = LLONG_MAX;
			for( int v = t, u = par[v]; v != s; u = par[v = u] )
				bot = min(bot, fnet[v][u] ? fnet[v][u] : ( cap[u][v] - fnet[u][v] ));
			for( int v = t, u = par[v]; v != s; u = par[v = u] )
				if( fnet[v][u] ) { fnet[v][u] -= bot; fcost -= bot * cost[v][u]; }
				else { fnet[u][v] += bot; fcost += bot * cost[u][v]; }
			flow += bot;
		}
		return make_pair(flow, fcost);
	}
};