POJ2175 Evacuation Plan 【 最小费用流+消圈定理】

 poj2175链接

题意：给出N个建筑物的坐标和里面的人数。给出M个防空洞的坐标和能容纳的人数。一个人的花费定义为从他的建筑物到目标防空洞的曼哈顿距离再加1。现在有人已经设计了一个避难方案。要你判断他给出的方案是否是所有人花费最少的解决方案，如果不是，请给出一个花费更小的解决方案（不必是最优）。

消圈定理

所谓消圈定理，就是在某个流  中，如果其对应的残余网络没有负圈（剩余流量为  的边视为不存在），那它一定就是当前流量下的最小费用流。反之亦然。即「 是最小费用流等价于其残余网络中没有负圈」。

证明

怎么证明？

假设任意流 ，如果在相同流量下有比他更小费用的流 。观察流 ，在流  中，除源汇外所有点的流入量等于流出量，在流  中亦然，即流  是由若干圈构成的。因为流  的费用是负的，所以在流  中，必定至少有一个负圈。

摘自日本人的书，可能有点不好理解。让我们用更直观的方式考虑。

假设一个网络是这样的，为方便起见，所有边的容量都是 。

如果流量走上路的话，其残余网络（没有特别说明的话，考虑的都是残余网络）看起来应该是这样的：

因为上路的边的流量占满了，所以现在上路只有反边。明显看到  是一个大大的负圈。我们沿着此负圈增广（每条边的流量＋1），可以得到当前流量下的最小费用流。
正确性很好证明。原流中，除了源点以及汇点每个点的入流量等于出流量，这个流一定是可行流。沿着负圈增广之后，除了源点以及汇点每个点的入流量仍然等于出流量。
就拿顶点  来说，原图中  的出入流量均为 ，增广之后，有  的流量从  流入 ，但是，却有  的流量从  流回 。
一番下来， 的出入流量仍然相等，还是 。不难证明其它顶点均如此，所以消圈后既保证了这个流是可行流，又减少了费用，所以它是正确的。
为了加深理解，对应下图，流量在圈中增广，总的流量既没有增加，也没有减少，只不过是流量从费用更少的地方流过 （），从费用大的地方退流而已（），流过的流量和退掉的流量是相等的。实质上是将从  流出的流量的方向改变，使得费用更小。
网络流的反边给了我们一个很好的反悔机制，使得我们可以对任意一个流 ，通过消负圈（可能不止一个），来得到它当前流量下的最小费用流。

可以看到，沿着负圈增广之后，已经没有负圈存在了，已经达到了当前流量下的最小费用流（也就是最小费用最大流）。所以只要有负圈，就可以增广达到更小费用。反之亦然，所以有「 是最小费用流等价于其残余网络中没有负圈」

应用

前面说了对任意一个流 ，可以通过消负圈，来得到它当前流量下的最小费用流。问题来了，怎么找负圈？残余网络可能不是联通的，因此从谁出发都不确定。既然不确定，我们添加一个超级源点 ，给每一个顶点  连一条从  的有向边，权值为 ，容量无限大（想一想为什么不能是无向边），从  出发，令 。这样就保证了图的联通性。
再将  算法执行  次， 为顶点数（含 ）。同时记录每个顶点的前驱节点 ，以及对应的边 。如果第  次仍然有顶点被更新，那么残余网络中一定有负圈。
如果没有负圈，最坏的情况是链式结构，每次只更新一个顶点，共需要  次，所以如果第  次仍然被更新，那说明被更新的顶点  的前驱中一定有负圈（它不一定在负圈上），那我们只需要不断的找它的前驱节点，每经过一个节点标记一次，那么一旦一个节点  被经过  次，节点  一定在负圈中。

但用不着这么麻烦，既然超级源点到每个顶点的边权都是 ，我们不妨直接让所有顶点的  直接为 ，这样是等效的，而且避免的加边带来的麻烦。具体操作就是让所有的 ，顶点数量  还是原来的顶点数（ 已经不存在了），然后跟之前一样用  就行了。
当然  也可以，不过条件要改成顶点入队的次数为 。

遍历负圈也是一样的方法，从刚刚找到的在负圈的节点  开始，通过  找到这条边让其流量 +1（别忘了修改反边），再走到前驱节点，如果当前节点又是  则跳出。具体细节见代码。

#include<cstring>#include<cstdio>#include<iostream>#include<algorithm>#include<vector>#include<queue>using namespace std;const int INF=0x3f3f3f3f;const int N=666;   int abs(int a){        return a>0?a:-a;}struct build{    int x,y,num;};struct edge{    int to,cap,cost;    edge(int x=0,int y=0,int z=0)    {        to=x;cap=y;cost=z;    }};vector<edge> g[N];build bu[N];int cnt[N],d[N],n,m,M,pre[N],nu[N],mp[N][N];bool vis[N];int costij(int i,int j){    return abs(bu[i].x-bu[j].x)+abs(bu[i].y-bu[j].y)+1;}void add(int from,int to,int cap1,int cap2,int cost){    if(cap1>0)    g[from].push_back(edge(to,cap1,cost));    if(cap2>0)    g[to].push_back(edge(from,cap2,-cost));}int spfa(int s){    for(int i=0;i<=M;i++)        d[i]=INF;    d[s]=0;    memset(vis,false,sizeof(vis));    memset(cnt,0,sizeof(cnt));    memset(pre,-1,sizeof(pre));    queue<int> que;    que.push(s);    vis[s]=true;    cnt[s]=1;    while(!que.empty())    {        int u=que.front();que.pop();        //cout<<g[u].size()<<"  &&&"<<endl;        for(int i=0;i<g[u].size();i++)        {            edge &e=g[u][i];            //cout<<e.cap<<"   ************"<<d[e.to]<<endl;            if(d[e.to]>d[u]+e.cost)            {                d[e.to]=d[u]+e.cost;                pre[e.to]=u;                if(!vis[e.to])                {                    que.push(e.to);                    vis[e.to]=true;                    cnt[e.to]++;                    if(cnt[e.to]>M)return e.to;                }            }        }        vis[u]=false;    }    return -1;}int main(){    while(scanf("%d%d",&n,&m)!=EOF)    {        M=n+m;        for(int i=0;i<=M;i++)            g[i].clear();        for(int i=1;i<=M;i++ )        {            scanf("%d%d%d",&bu[i].x,&bu[i].y,&bu[i].num);        }        memset(nu,0,sizeof(nu));        for(int i=1;i<=n;i++)        for(int j=1;j<=m;j++)        {            scanf("%d",&mp[i][j]);            int cost0=costij(i,j+n);            add(i,j+n,INF,mp[i][j],cost0);            nu[j]+=mp[i][j];        }        for(int j=1;j<=m;j++)        {            add(j+n,0,bu[j+n].num-nu[j],nu[j],0);        }        int u=spfa(0);        if(u==-1)        {            printf("OPTIMAL\n");        }        else        {            printf("SUBOPTIMAL\n");            memset(cnt,0,sizeof(cnt));            int v=u,now;               while(cnt[pre[v]]==0)            {                cnt[v]++;                v=pre[v];            }            now=v;            do            {                if(pre[now]<=n&&now>n)                    mp[pre[now]][now-n]++;                else if(pre[now]>n&&now<=n)                     mp[now][pre[now]-n]--;                now=pre[now];            }while(now!=v);            for(int i=1;i<=n;i++)            {                for(int j=1;j<m;j++)                    printf("%d ",mp[i][j]);                printf("%d\n",mp[i][m]);            }        }    }    return 0;}