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;
}