「POJ2175」Evacuation Plan

题目

点这里看题目。

分析

显而易见的方法：直接建图跑一个最小费用最大流，然后比较自己得到的结果和给出的方案。

但是这里有 \(O(n^2)\) 条边，再加上流量可以被构造得很大，因此这种方法不出所料地超时了。

注意到，题目给出的方案一定是一个最大流的方案，但不一定是最小费用的方案。换句话说，题目其实给出了一个残量网络。根据完全不众所周知的消圈原理，如果这个残量网络上存在负环，那么这个流就一定不是最小费用流。另外，如果搜索到了一个负环，我们也可以在这个负环上调整 1 的流量，从而得到一种更优的方案，因此也可以很容易地构造出方案。

理论的复杂度应为 \(O(n^2m)\)，勉强能够过去。如果使用 DFS 版本的 SPFA 并加上若干优化，那么就能跑得飞快。

补充一下消圈算法：

消圈算法就是几乎直接由消圈原理得出的算法。

我们首先找出一个合法的最大流方案，而后不断地调整负环，使得残量网络上再也找不出负环，就得到了最小费用最大流。

不同于常用的 SSP 方法导出的算法，无论原网络上有无负环，这个算法都可以正常地运行。

据传，精细实现后的消圈算法，复杂度可以降到 \(O(nm^2\log n)\)，这个算法是多项式的！

小结：仔细观察给出信息的含义，并充分利用它。

代码

#include <queue>
#include <cstdio>

#define rep( i, a, b ) for( int i = (a) ; i <= (b) ; i ++ )
#define per( i, a, b ) for( int i = (a) ; i >= (b) ; i -- )

const int INF = 0x3f3f3f3f;
const int MAXN = 1005, MAXE = 1e5 + 5;

template<typename _T>
void read( _T &x ) {
    x = 0; char s = getchar(); bool f = false;
    while( s < '0' || '9' < s ) { f = s == '-', s = getchar(); }
    while( '0' <= s && s <= '9' ) { x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar(); }
    if( f ) x = -x;
}

template<typename _T>
void write( _T x ) {
    if( x < 0 ) putchar( '-' ), x = -x;
    if( 9 < x ) write( x / 10 );
    putchar( x % 10 + '0' );
}

template<typename _T>
_T ABS( const _T a ) {
    return a < 0 ? -a : a;
}

struct Edge {
    int to, nxt, c, w;
}Graph[MAXE << 1];

int cir[MAXN], tot;
int stk[MAXN], top;

int dist[MAXN], dep[MAXN];
bool vis[MAXN];

int head[MAXN];
int ntot, cnt = 1;

int E[105][105];
int X[MAXN], Y[MAXN], B[MAXN];
int P[MAXN], Q[MAXN], C[MAXN];

int N, M;

inline void AddEdge( const int from, const int to, const int C, const int W ) {
    Graph[++ cnt].to = to, Graph[cnt].nxt = head[from];
    Graph[cnt].c = C, Graph[cnt].w = W, head[from] = cnt;
}

bool DFS( const int u ) {
    vis[u] = true;
    stk[dep[u] = ++ top] = u;
    for( int i = head[u], v ; i ; i = Graph[i].nxt )
        if( Graph[i].c && dist[v = Graph[i].to] > dist[u] + Graph[i].w ) {
            dist[v] = dist[u] + Graph[i].w;
            if( ! vis[v] ) {
                if( DFS( v ) )
                    return true;
            }
            else {
                rep( j, dep[v], dep[u] )
                    cir[++ tot] = stk[j];
                return true;
            }
        }
    vis[u] = false, top --;
    return false;
}

int main() {
    read( N ), read( M ), ntot = N + M;
    const int s = ++ ntot, t = ++ ntot;
    rep( i, 1, N ) 
        read( X[i] ), read( Y[i] ), read( B[i] );
    rep( i, 1, M ) 
        read( P[i] ), read( Q[i] ), read( C[i] );
    rep( i, 1, N ) rep( j, 1, M ) {
        read( E[i][j] );
        int wei = ABS( X[i] - P[j] ) + ABS( Y[i] - Q[j] ) + 1;
        AddEdge( i, j + N, INF, wei ), AddEdge( j + N, i, E[i][j], - wei );
    }
    rep( i, 1, N ) {
        int su = 0;
        rep( j, 1, M ) su += E[i][j];
        AddEdge( s, i, B[i] - su, 0 );
        AddEdge( i, s, su, 0 );
    }
    rep( j, 1, M ) {
        int su = 0;
        rep( i, 1, N ) su += E[i][j];
        AddEdge( j + N, t, C[j] - su, 0 );
        AddEdge( t, j + N, su, 0 );
    }
    rep( i, 1, ntot ) {
        rep( j, 1, ntot ) dist[j] = 0, vis[j] = false;
        tot = top = dist[i] = 0; 
        if( DFS( i ) ) break;
    }
    if( tot == 0 ) puts( "OPTIMAL" );
    else {
        puts( "SUBOPTIMAL" );
        cir[++ tot] = cir[1];
        for( int i = 1 ; i < tot ; i ++ ) {
            int u = cir[i], v = cir[i + 1];
            if( u <= N && N < v && v <= N + M ) E[u][v - N] ++;
            if( v <= N && N < u && u <= N + M ) E[v][u - N] --;
        }
        rep( i, 1, N ) rep( j, 1, M )
            write( E[i][j] ), putchar( j == M ? '\n' : ' ' );
    }
    return 0;
}