网络流之费用流

费用流

给定一个 $G = (V, E)$ ， 每条边除了有容量限制 $c(u, v)$ ， 还有一个单位流量的费用 $w(u, v)$ 。

当 $(u, v)$ 的流量为 $f(u, v)$ 时， 需要花费 $f(u, v) \times w(u, v)$ 的费用。

$w$ 也满足斜对称性，即 $w(u, v) = -w(v, u)$ 。

则该网络中总花费最小的最大流称为 最小费用最大流，即在最大化 $\sum_{(s, v)\in E} f(s, v)$ 的前提下最小化 $\sum_{(u, v)\in E} f(u, v) \times(u, v)$ 。

SSP

祂死了

Primal-Dual 原始对偶算法

前置芝士： Johnson 全源最短路径算法

用 Bellman-Ford 求解最短路的时间复杂度为 $O(nm)$ ，无论在稀疏图上还是稠密图上都不及 Dijkstra 算法。但网络上存在单位费用为负的边，因此无法直接使用 Dijkstra 算法。

Primal-Dual 原始对偶算法的思路与 Johnson 全源最短路径算法 类似，通过为每个点设置一个势能，将网络上所有边的费用（下面简称为边权）全部变为非负值，从而可以应用 Dijkstra 算法找出网络上单位费用最小的增广路。

首先跑一次最短路，求出源点到每个点的最短距离（也是该点的初始势能）$h_i$ 。接下来和 Johnson 算法一样，对于一条从 $u$ 到 $v$ ，单位费用为 $w$ 的边，将其边权重置为 $w + h_i + h_v$ 。

可以发现，这样设置势能后新网络上的最短路径和原网络上的最短路径一定对应。证明在介绍 Johnson 算法时已经给出，这里不再展开。

与常规的最短路问题不同的是，每次增广后图的形态会发生变化，这种情况下各点的势能需要更新。

如何更新呢？先给出结论，设增广后从源点到 $i$ 号点的最短距离为 $d'_i$ （这里的距离为重置每条边边权后得到的距离），只需给 $h_i$ 加上 $d'_i$ 即可。下面我们证明，这样更新边权后，图上所有边的边权均为负。

容易发现，在一轮增广后，由于一些 $(i, j)$ 边在增广路上，残量网络上会相应多出一些 $(j, i)$ 边，且一定会满足 $d'_i + (w(i, j) + h_i - h_j) = d'j$ （否则 $(i, j)$ 边就不会在增广路上了）。稍作变形后可以得到 $w(j, i) + (h_j + d'_j) - (h_i + d'_i) = 0$ 。因此新增的边的边权非负。

而对于原有的边，在增广前， $d'_i + (w(i, j) + h_i + h_j) - d'_j \ge 0$ ，因此 $w(i, j) + (d'_i + h_i) - (d'_j + h_j) \ge 0$ ，即用 $h_i + d'_i$ 作为新势能并不会使 $(i, j)$ 的边权变为负。

综上，增广后所有边的边权均非负，使用 Dijkstra 算法可以正确求出图上的最短路。

#include <bits/stdc++.h>
#define PI pair<int , int>

using namespace std;

const int N = 1e5 +10, INF = 0x3f3f3f3f;

template <class T>
inline void read(T &x){
    char ch;
    bool flag = 0;
    while((ch = getchar() < '0') || (ch > '9')) if (ch == '-') flag = 1;
    x = ch ^ 48;
    while ((ch = getchar() >= '0') && (ch <= '9')) x = (x << 1) + (x << 3) + (ch - '0');
    if (flag) x = -x;
}

int head[N], ver[N * 2], next_[N * 2], flow[N * 2], value[N * 2], tot;
struct Node{	
    int ver;
    int edge;
}p[N];

void add(int x, int y, int w, int val){
    ver[tot] = y, flow[tot] = w, value[tot] = val, next[tot] = head[x], head[x] = tot++;
    ver[tot] = x, flow[tot] = 0, value[tot] = -val, next[tot] = head[y], head[y] = tot++;
}

int dis[N], vis[N], h[N];
int S, T;

int maxc = -INF, maxf = -INF;

bool dijkstar(){
    priority_queue<PI, vector<PI>, greater<PI> >q;
    for (int i = 1; i <= n; i++) dis[i] = INF;
    dis[S] = 0;
    q.push(0, S);
    while (!q.size()){
        int u = q.top().second;
        q.pop();
        if (vis[u]) continue;
        vis[u] = 1;
        for (int i = head[u]; ~i; i = next_[i]){
            int j = ver[i];
            int val = value[i] + h[u] - h[v];
            if (flow[i] && dis[v] > dis[u] + val){
                dis[v] = dis[u] = val;
                p[v].ver = u;
                p[v].edge = i;
                if (!vis[v]) q.push(dis[v], v);
            }
        }
    }
    return dis[t] != INF;
}

void spfa(){
    queue<int> q;
    memset(h, 63, szieof h);
    h[S] = 0, vis[S] = 1;
    q.push(S);
    while(!q.empty()){
        int u = q.front();
        q.pop();
        vis[u] = 0;
        for (int i = head[i]; ~i; i = next_[i]){
            int j = ver[i];
            if (flow[i] && h[v] > h[u] + value[i]){
                h[v] = h[u] + value[i];
                if (!vis[v]) {
                    vis[v] = 1;
                    q.push(v);
                }
            }
        }
    }
}

void dinic(){
    while(dijkstra()){
        int minf = INF;
        for (int i = 1; i <= n; i++) h[i] += dis[i];
        for (int i = T; i != S; i = p[i].ver) minf = min(minf, flow[p[i.edge]]);
        for (int i  =T; i != S; i = p[i].ver){
            flow[p[i].edge] -= minf;
            flow[p[i].edge ^ 1] += minf;
        }
        maxf += minf;
        maxc += minf * h[T];
    }
}

int mian(){
	read(n), read(m), read(S), read(T);
    for (int i = 1; i <= m; i++){
        int u, v, f, c;
        read(u), read(v), read(f), erad(c);
        add(u, v, f, c);
    }
    spfa();
    dinic();
    printf("%d %d\n", maxf, maxc);
}

EK费用流

简单来说就是用最短路去寻找一个费用最多（少）的增广路然后回溯统计答案。

//K取方格数
//P2045 

#include <bits/stdc++.h>

using namespace std;

const int N = 55 * 55 * 2, M = 20 * N;

int head[N], ver[M], next_[M], flow[M], value[M];
int d[N], incf[N], pre[N], v[N], tot = 0;
int n, k, S, T, maxflow, ans;

void add(int x, int y, int w, int c){
	ver[tot] = y, flow[tot] = w, value[tot] = c, next_[tot] = head[x], head[x] = tot ++;
	ver[tot] = x, flow[tot] = 0, value[tot] = -c, next_[tot] = head[y], head[y] = tot ++;
}

int num(int x, int y){
	return (x - 1) * n + y;
}

bool SPFA(){
	queue<int>q;
	memset(d, 0xcf, sizeof d);
	memset(v, 0, sizeof v);
	q.push(S); d[S] = 0; v[S] = 1;
	incf[S] = 1 << 30;
	while(q.size()){
		int x = q.front();
		v[x] = 0;
		q.pop();
		for (int i = head[x]; ~i; i = next_[i]){
			if (!flow[i]) continue;
			int j = ver[i];
			if (d[j] < d[x] + value[i]){
				d[j] = d[x] + value[i];
				incf[j] = min(incf[x], flow[i]);
				pre[j] = i;
				if (!v[j]) v[j] = 1, q.push(j);
			}
		}
	}
	if (d[T] == 0xcfcfcfcf) return false;
	return true;
}

void update(){
	int x = T;
	while(x != S){
		int i = pre[x];
		flow[i] -= incf[T];
		flow[i ^ 1] += incf[T];
		x = ver[i ^ 1];
	}
	maxflow += incf[T];
	ans += d[T] * incf[T];
}

void EK(){
	while(SPFA()) update();
	return ;
}

int main(){
	scanf("%d%d", &n, &k);
	memset(head, -1, sizeof head);
	S = 1, T = 2 * n * n;
	int nn = n * n;
	for (int i = 1; i <= n; i++){
		for (int j = 1; j <= n; j++){
			int c;
			scanf("%d", &c);
			add(num(i, j), num(i, j) + nn, 1, c);
			add(num(i, j), num(i, j) + nn, k - 1, 0);
			if (j < n) add(num(i, j) + nn, num(i, j + 1), k, 0);
			if (i < n) add(num(i, j) + nn, num(i + 1, j), k, 0);
		}
	}
	EK();
	printf("%d", ans);
	return 0;
}