网络流-消圈定理

想法

如何判定在当前流量下，一个费用流是否是最小费用流？

这个问题等同于，是否存在一种方案，在不改变总流量的情况下改变一些边的流量，最终减小费用。网络流的消圈定理就对此作出了解答。

消圈定理

可行流 \(f\) 是当前流量下的最小费用流 \(\Leftrightarrow\) 残余网络不存在负环

把单点看成一个长度为 0 的环。

设 \(f'\) 是一个流量不变，费用更小的可行流，那么 \(f'-f\) 的每个点都满足流入流量都等于流出流量。这是因为 \(f,f'\) 的除源点汇点的其他点都分别满足这个条件，而总流量不变，所以源汇的流量也平衡。

这时 \(f'-f\) 是由一些环组成的。若不存在环，那么不可能满足流量平衡；若除了环还存在一些链，那么链头和链尾也就不满足要求。

由于 \(w(f')<w(f)\) ，那么 \(w(f'-f)<0\) ，所以这些环中一定存在一个负环。

应用

有了这个定理，我们就能判定当前流量下费用流是否是最小费用流了——若有负环就不是，否则是。

这也给我们提供了一个得到最小费用流的方法，即不断在负环上增广，直到找不到负环 。

直接做的复杂度上界为 \(O(nm^2cw)\) ，其中 \(c,w\) 分别为最大容量和最大费用，整个过程可能执行 \(mcu\) 次，一次最坏找到负环为 \(O(nm)\) 。

代码

POJ2175 学会了这个方法就很简单了。

用spfa来判断负环，用队列实现 600ms 左右，换成栈变成 300ms 左右。然而dfs TLE 了。讲道理，栈不就是dfs吗？！

#include<cstdio>
#include<cctype>
#include<climits>
#include<cstring>
#include<algorithm>
using namespace std;
inline char nchar() {
	static const int bufl=1<<20;
	static char buf[bufl],*a,*b;
	return a==b && (b=(a=buf)+fread(buf,1,bufl,stdin),a==b)?EOF:*a++;
}
inline int read() {
	int x=0,f=1;
	char c=nchar();
	for (;!isdigit(c);c=nchar()) if (c=='-') f=-1;
	for (;isdigit(c);c=nchar()) x=x*10+c-'0';
	return x*f;
}
const int maxn=205;
const int maxe=maxn*maxn*2;
int n,m,bt[maxn],ct[maxn];
struct D {
	int x,y,d;
	inline int operator - (const D &b) const {
		return abs(x-b.x)+abs(y-b.y)+1;
	}
} b[maxn],c[maxn];
namespace graph {
	struct edge {
		int u,v,c,w,nxt;
	} e[maxe];
	int h[maxn],tot=1,pre[maxn],d[maxn],cnt[maxn],que[maxe],ql,qr; // pre e
	bool inq[maxn],vis[maxn];
	inline int add(int x,int y,int c,int w) {
		e[++tot]=(edge){x,y,c,w,h[x]};
		return h[x]=tot;
	}
	void print(int x) {
		for (int i=x;;i=e[pre[i]].u) {
			const int &k=pre[i];
			--e[k].c,++e[k^1].c;
			if (e[pre[i]].u==x) break;
		}
		puts("SUBOPTIMAL");
		for (int i=1;i<=n;++i) {
			static int a[maxn];
			for (int j=h[i],v=e[j].v;j;j=e[j].nxt,v=e[j].v) a[v-n]=e[j^1].c;
			for (int j=1;j<=m;++j) printf("%d%c",a[j]," \n"[j==m]);
		}
	}
	inline int find(int x) {
		for (int i=x;;i=e[pre[i]].u) if (!vis[i]) vis[i]=true; else return i;
	}
	inline void cancle(int s,int t) {
		memset(d,0x3f,sizeof d);
		d[t]=0;
		que[ql=qr=1]=t;
		while (ql<=qr) {
			int x=que[qr--];
			for (int i=h[x],v=e[i].v;i;i=e[i].nxt,v=e[i].v) if (e[i].c && d[v]>d[x]+e[i].w) {
				d[v]=d[x]+e[i].w;
				pre[v]=i;
				if (inq[v]) continue;
				inq[v]=true;
				if (++cnt[que[++qr]=v]>t+1) {
					int p=find(v);
					print(p);
					exit(0);
				}
			}
			inq[x]=false;
		}
	}
}
using graph::add;
int main() {
#ifndef ONLINE_JUDGE
	freopen("test.in","r",stdin);
#endif
	n=read(),m=read();
	int S=0,T=n+m+1;
	for (int i=1;i<=n;++i) {
		b[i].x=read(),b[i].y=read();
		bt[i]=add(S,i,b[i].d=read(),0);
		add(i,S,0,0);
	}
	for (int i=1;i<=m;++i) {
		c[i].x=read(),c[i].y=read();
		ct[i]=add(n+i,T,c[i].d=read(),0);
		add(T,n+i,0,0);
	}
	for (int i=1;i<=n;++i) for (int j=1;j<=m;++j) {
		int d=b[i]-c[j],x=read();
		add(i,n+j,INT_MAX-x,d);
		add(n+j,i,x,-d);
		graph::e[bt[i]].c-=x,graph::e[bt[i]^1].c+=x;
		graph::e[ct[j]].c-=x,graph::e[ct[j]^1].c+=x;
	}
	graph::cancle(S,T);
	puts("OPTIMAL");
	return 0;
}