kitsune- ICPC Minimum Cost Flow Library

2007-09-19 / /

Features

Team kitsune- による C++ フルスクラッチ実装の最小費用流ライブラリです。 本ライブラリには以下のような特徴があります:
高速である
アルゴリズムの工夫と実装の最適化により高速に動作します。 UVA 3562 にて 0.205 秒 (2位)。 速度チューニングのための無駄なスイッチ類があります。
データ構造が隣接リスト
隣接行列でないので、点数の多い疎なグラフや平行枝などにも対応しています。
スケーリングに対応
最小重みマッチングなど、条件を満たす場合にはスケーリングが可能です。 (ICPCでは使うことはないと思いますが)
長い
速度が足りなくて悔しい思いをするよりは長いのをがんばって打ち込むという方針にしたので、 長いです。

Algorithm

Primal-Dual 法。 詳細はとりあえず今は省きますので Reference 等を見てください。 以下おおまかな説明。
[定義] フロー f と ポテンシャル p は、reduced cost を枝長としたグラフ上で負閉路がない時に equibilium であるという。
[定理] (f, p)が equibilium にあるならば、フロー f は現在流れているフロー量に関して最小費用流である(証明略)。
最初に equibilium にある (f, p) から出発し、 equibilium を保つようにポテンシャルを更新しフローを流していく。 最終的に要求されたフロー量まで到達したら、上の定理より、その時のフロー f が最小費用流となる。
  1. 初期の equibilium な (f, p) を作る。
    本ライブラリでは、 p = 0、 負コスト枝には容量いっぱいのフローを流し、 それ以外の枝のフローは0、 とすることで equibilium としている。
  2. reduced cost を枝長としたグラフ上で Dijkstra して、ポテンシャル p を更新する。
  3. reduced cost が 0 になった枝にフローを流す。 普通のアルゴリズムでは、ここであるsourceからあるsinkまでの最短路に沿ってフローを流す。 本ライブラリでは複数のsource-sink間のパスに沿って複数のフローを流すことで、多くの場合高速化できる。
  4. 目指すフロー量 b に達するまで、2, 3 を繰り返す。

Complexity

計算量は、流すフロー総計をB、枝数e、点数nとして O(B * e log n)。 スケーリング時は O((log B) * n * e log n)。

Usage

Function Prototype

template<class V, class COST>
pair<bool, COST>
mincost(
int n,
const V *b0, V *b,
const vector< vector< int > > &adj,
edge<V,COST> *e, int m,
COST inf,

COST *d, COST *dr, COST *p,
int *prev, int *visit
);

Parameters

V - 容量/フローの型。整数型限定。
COST - コストの型。
n [IN] 点数。
b0 [IN] 流したいフロー。 b0[i] > 0 : source, b0[i] < 0 : sink。
b [OUT]流したいフローの残り。全部流せた場合には全0。
adj [IN] adj[i] : 点iに接続する枝IDのリスト。
e [IN] 枝。枝ID=jは e[j]
m [IN] 枝数。
inf [IN] 十分大きな数。
d , dr , p , prev , visit [BUF]作業用領域。各々 n 要素分確保すること。初期化は必要ない。
枝のデータ構造は以下の通り。点 i →点 j の、コスト cost2 、容量 cap の有向枝。 flow は初期化する必要はない。実行後、実際に流したフローが入る。
class edge{
public:
    int i, j;
    COST cost2;
    V flow, cap;
};

Return Value

(全部流せたか, 最小コスト) のペアを返す。

Notes

Switches

WL
1 = Dijkstra に priority_queue を使用する(デフォルト)。O(e log n)
0 = Dijkstra に配列を用いる。グラフが密な場合には速いことがある O(n^2)
SCALE
1 = スケーリング有り。但し、全枝容量無限大、強連結、コスト非負であること。
0 = スケーリング無し。
SINGLEFLOW
1 = 1回のダイクストラで複数フローを流す
0 = 1回のダイクストラでフローを1つだけ流す(普通のPrimal-Dual)
OPT
プログラムには2バージョンあり、それを切り替える。
1 = 最適化版。速度最優先、長い。
0 = short版。やや速度を落としてコードサイズ削減。

Usage Example

UVA 10594: Data Flow
// paste here the library

typedef long long ll;
ll  d[1000], dr[1000], p[1000];
int prev[1000], visited[1000], b0[1000], b[1000];
edge<int, ll> e[10000];
main(){
    while( 1 ){
        int N, M, K, D;
        scanf( "%d%d", &N, &M );
        if( feof( stdin ) ) break;
        vector< vector< int > > adj( N );
        for( int i = 0; i < M; i ++ ){
            int a, b, j;
            scanf( "%d%d%d", &a, &b, &j ); a --, b --;
            { edge<int,ll> &est = e[i*2+0]; est.i = a, est.j = b, est.cap = 1, est.cost2 = j; }
            { edge<int,ll> &est = e[i*2+1]; est.i = b, est.j = a, est.cap = 1, est.cost2 = j; }
        }

        scanf( "%d%d", &D, &K );
        for( int i = 0; i < 2 * M; i ++ )
            e[i].cap = K, adj[e[i].i].push_back( i ), adj[e[i].j].push_back( i );

        for( int i = 0; i < N; i ++ ) b0[i] = 0;
        b0[0]   = +D; b0[N-1] = -D; // set the demand flow D from 0 to N-1

        pair<bool, ll> r = mincost<int,ll>( N, b0, b, adj, e, 2 * M, 1LL << 62, d, dr, p, prev, visited );
        if( r.first ) // success
            cout << r.second << endl;
        else // failed
            cout << "Impossible." << endl;
    }
}

Validation and Execution Time

実行時間の単位は秒。
SCALE WL SINGLEFLOW UVA 3562 UVA 10594 UVA 10746 PKU 2195
short opt short opt short opt short opt
0000.2710.2111.4901.264<0.012<0.0120.0150.015
0010.8180.6191.5181.199<0.012<0.0120.1090.078
0100.4160.3522.1071.750<0.012<0.0120.0150.031
0111.2771.0942.0681.744<0.012<0.0120.2340.375
1000.5230.371N/A N/A <0.012<0.0120.0150.015
1011.0290.738N/A N/A <0.012<0.0120.1250.171
1100.8420.666N/A N/A <0.012<0.0120.0310.015
1111.7171.406N/A N/A <0.012<0.0120.2650.421

References

Licence

NYSL Version 0.9982

A. 本ソフトウェアは Everyone'sWare です。このソフトを手にした一人一人が、
   ご自分の作ったものを扱うのと同じように、自由に利用することが出来ます。

  A-1. フリーウェアです。作者からは使用料等を要求しません。
  A-2. 有料無料や媒体の如何を問わず、自由に転載・再配布できます。
  A-3. いかなる種類の 改変・他プログラムでの利用 を行っても構いません。
  A-4. 変更したものや部分的に使用したものは、あなたのものになります。
       公開する場合は、あなたの名前の下で行って下さい。

B. このソフトを利用することによって生じた損害等について、作者は
   責任を負わないものとします。各自の責任においてご利用下さい。

C. 著作者人格権は 林崎弘成 に帰属します。著作権は放棄します。

D. 以上の3項は、ソース・実行バイナリの双方に適用されます。

Source Code

Download

Download the source code (graph_mincost.cpp)

short版

#include<stdio.h>
#include<iostream>
#include<stdlib.h>
#include<string.h>
#include<string>
#include<vector>
#include<set>
#include<map>
#include<queue>
#include<algorithm>
#include<utility>
#include<complex>
#include<float.h>
#include<limits.h>
#include<list>
using namespace std;

// (SCALE,WL,SINGLEFLOW) = 0, 0/1, 0 が普通最速, 0, 1, 1 が割と安全
#define SCALE      0
#define WL         1
#define SINGLEFLOW 1

#define OPT        1

/****************************************************************

    short: short版


****************************************************************/

#if OPT == 0

template<class V, class COST>
class edge{
public:
    int i, j;
    COST cost2;
    V flow, cap;
    int dest( int f){ return i == f ?
        j          : i;         }
    V residue(int f){ return i == f ?
#if SCALE
        INT_MAX    : flow;      }
#else
        cap - flow : flow;      }
#endif
    COST cost(int f){ return i == f ?
        cost2      : -cost2;    }
    void flow_add( int f, V r ){ i == f ?
        flow += r  : flow -= r; }
};

template<class V, class COST>
pair<bool, COST>
mincost(
int n,
const V *b0, V *b,
const vector< vector< int > > &adj,
edge<V,COST> *e, int m,
COST inf,

// バッファ領域; 各々 n 要素分確保すること
COST *d, COST *dr, COST *p,
int *prev, int *visit
)
{
    int n1 = n - 1;

// initial equibilium flow and potential
    for( int i = n1; i >= 0; -- i ){
        b[i] = b0[i];
        p[i] = 0;
    }
    for( int i = m - 1; i >= 0; -- i ){
        if( e[i].cost2 < 0 ){
            V r = e[i].cap;
            e[i].flow = r;
            b[e[i].i] -= r;
            b[e[i].j] += r;
        }
        else
            e[i].flow = 0;
    }

    V delta = 1;

#if SCALE
    V bmax = 0;
    for( int i = n1; i >= 0; -- i ){
        if( b[i] < 0 )
            bmax = max( bmax, - b[i] );
        else
            bmax = max( bmax, + b[i] );
    }
    while( delta <= (bmax>>1) ) delta <<= 1;
#endif

    for( ; delta >= 1; delta >>= 1 ){
    while( 1 ){
#if WL
        priority_queue< pair<COST, int> > wl;
#endif

        for( int i = n1; i >= 0; -- i ){
            d[i] = inf;
            if( b[i] >= delta ){
#if WL
                wl.push( pair<COST, int>( -0, i) );
#endif
                d[i] = 0;
            }
        }
#if WL
        if( wl.empty() )
            break;
#endif

        memset( prev,  0xff, sizeof(int)*n );
        memset( visit, 0x00, sizeof(int)*n );

        int t = -1;

#if WL
        while( !wl.empty() ){
            int i   =   wl.top().second;
            COST di = - wl.top().first;
            wl.pop();
            if( visit[i] )
                continue;
#else
        while( 1 ){
            int i = min_element( d, d + n ) - d;
            COST di = d[i];
            if( di == inf )
                break;
#endif
            COST dbase = di + p[i];
            dr[i]      = di;
            d[i]       = inf;
            visit[i] = 1;
#if SINGLEFLOW
            if( t < 0 && b[i] <= -delta )
                t = i;
#endif

            for( vector<int>::const_iterator
            J=adj[i].begin(), JE=adj[i].end();
            J != JE; ++ J ){
                edge<V,COST> &est = e[*J];
                int j = est.dest( i );
                if( visit[j] || est.residue(i)<=0 )
                    continue;
                COST dj = dbase+est.cost(i)-p[j];
                if( d[j] > dj ){
                    d[j] = dj;
                    prev[j] = *J;
#if WL
                    wl.push(pair<COST,int>(-dj, j));
#endif
                }
            }
        } // dij end

        for( int i = n1; i >= 0; -- i )
            if( visit[i] ) p[i] += dr[i];

#if SINGLEFLOW
        if( t >= 0 ){
#else
        for( int t2 = n1; t2 >= 0; -- t2 ){
            if( !visit[t2] || b[t2] > -delta )
                continue;
            t = t2;
#endif
            int i, j;

//--------
// r を決定
#if SCALE && SINGLEFLOW
            V r = delta;
#else
    #if SCALE
            V r = delta;
    #else
            V r = -b[t];
    #endif
            for( i = t; (j = prev[i]) >= 0; ){
                edge<V,COST> &est = e[j];
                i = est.dest( i );
                r = min( r, est.residue( i ) );
            }
            r = min( r, b[i] );
#endif
            if( r < delta )
                continue;

//--------
// r だけ流す
            for( i = t; (j = prev[i]) >= 0; ){
                edge<V,COST> &est = e[j];
                i = est.dest( i );
                est.flow_add( i, r );
            }
            b[t]  += r;
            b[i]  -= r;
        }

        if ( t == -1 ) // no feasible flow
            break;

    }}

// コスト総計を計算
    COST cost_sum = 0;
    bool is_b = true;
    for( int i = m - 1; i >= 0; -- i )
        cost_sum += e[i].cost2 * e[i].flow;
    for( int i = n - 1; i >= 0; -- i ){
        if( b[i] != 0 ){
            is_b = 0;
            break;
        }
    }

    return pair<bool, COST>(is_b, cost_sum);
}
#endif

最適化版

/****************************************************************

    opt: 最適化版


****************************************************************/

#if OPT == 1

template<class V, class COST>
class edge{
public:
    int i, j;
    COST cost2;
    V flow;
    V cap;
};

template<class V, class COST>
pair<bool, COST>
mincost(
int n,
const V *b0, V *b,
const vector< vector< int > > &adj,
edge<V,COST> *e, int m,
COST inf,

// バッファ領域; 各々 n 要素分確保すること
COST *d, COST *dr, COST *p,
int *prev, int *visited
)
{
    int n1 = n - 1;

// initial equibilium flow and potential
    for( int i = n1; i >= 0; -- i ){
        b[i] = b0[i];
        p[i] = 0;
    }
    for( int i = m - 1; i >= 0; -- i ){
        if( e[i].cost2 < 0 ){
            V r = e[i].cap;
            e[i].flow = r;
            b[e[i].i] -= r;
            b[e[i].j] += r;
        }
        else
            e[i].flow = 0;
    }

    V delta = 1;

#if SCALE
    V bmax = 0;
    for( int i = n1; i >= 0; -- i ){
        if( b[i] < 0 )
            bmax = max( bmax, - b[i] );
        else
            bmax = max( bmax, + b[i] );
    }
    while( delta <= (bmax >> 1) ) delta <<= 1;
#endif

    for( ; delta >= 1; delta >>= 1 ){ // until the flow becomes feasible

    int flag_t = 0, flag_s = 0;
    for( int i = n1; i >= 0; -- i ){
        if( b[i] >= delta )
            flag_s ++;
        else if( b[i] <= -delta )
            flag_t ++;
    }

    while( flag_t > 0 && flag_s > 0 ){
#if WL
        priority_queue< pair<COST, int> > wl;
#endif

        for( int i = n1; i >= 0; -- i ){
            d[i] = inf;
            if( b[i] >= delta ){
#if WL
                wl.push( pair<COST, int>( -0, i) );
#endif
                d[i] = 0;
            }
//            prev[i] = -1;
//            visited[i] = 0;
        }
#if WL
        if( wl.empty() )
            break;
#endif

        memset( prev, 0xff, sizeof(int)*n );
        memset( visited, 0x00, sizeof(int)*n );

        int t = -1;

#if WL
        while( !wl.empty() ){
            int i   =   wl.top().second;
            COST di = - wl.top().first;
            wl.pop();

            if( visited[i] )
                continue;
            COST dbase = di + p[i];
            dr[i]      = di;
#else
        while( 1 ){
            COST *jp = min_element( d, d + n );
            if( *jp == inf )
                break;
            int      i = jp - d;
            COST dbase = *jp + p[i];
            dr[i]      = *jp;
            *jp        = inf;
#endif

            visited[i] = 1;
#if SINGLEFLOW
            if( t < 0 && b[i] <= -delta ) // sink
                t = i;
#endif

            for( vector<int>::const_iterator J = adj[i].begin(), JE = adj[i].end(); J != JE; ++ J ){
                edge<V,COST> &est = e[*J];
                if( est.j == i ){
                    if( est.flow > 0 && !visited[est.i] ){
                        COST dj = dbase - est.cost2 - p[est.i];
                        if( d[est.i] > dj ){
                            d[est.i] = dj;
                            prev[est.i] = *J;
#if WL
                            wl.push( pair<COST, int>( -dj, est.i ) );
#endif
                        }
                    }
                }
                else{
#if SCALE
                    if( !visited[est.j] ){
#else
                    if( est.flow < est.cap && !visited[est.j] ){
#endif
                        COST dj = dbase + est.cost2 - p[est.j];
                        if( d[est.j] > dj ){
                            d[est.j] = dj;
                            prev[est.j] = *J;
#if WL
                            wl.push( pair<COST, int>( -dj, est.j ) );
#endif
                        }
                    }
                }
            }
        } // dij end

        for( int i = n1; i >= 0; -- i )
            if( visited[i] )
                p[i] += dr[i];

#if SINGLEFLOW
        {
            if( t >= 0 ){
                int t2 = t;
#else
        for( int t2 = n1; t2 >= 0 && flag_t > 0 && flag_s > 0; -- t2 ){
            if( visited[t2] && b[t2] <= -delta ){
                t = t2;
#endif
                int i, j;

//--------
// r を決定

#if SCALE && SINGLEFLOW
                V r = delta;
#else

#if SCALE
                V r = delta;
#else
                V r = -b[t2];
#endif
                for( i = t2; r > 0 && (j = prev[i]) >= 0; ){
                    edge<V,COST> &est = e[j];
                    if( est.i == i ){
                        r = min( r, est.flow );
                        i = est.j;
                    }
                    else{
#if !SCALE
                        r = min( r, est.cap - est.flow );
#endif
                        i = est.i;
                    }
                }
                r = min( r, b[i] );

#endif

//--------
// r だけ流す

                if( r >= delta ){
                    for( i = t2; (j = prev[i]) >= 0; ){
                        edge<V,COST> &est = e[j];
                        if( est.i == i ){
                            est.flow -= r;
                            i = est.j;
                        }
                        else{
                            est.flow += r;
                            i = est.i;
                        }
                    }
                    b[t2] += r;
                    b[i] -= r;
                    if( b[i] < delta )
                        flag_s --;
                    if( b[t2] > -delta )
                        flag_t --;
                }


            }
        }

        if ( t == -1 ) // no feasible flow
            break;

    }
    }

    COST cost_sum = 0;
    bool is_b = true;
    for( int i = m - 1; i >= 0; -- i )
        cost_sum += e[i].cost2 * e[i].flow;
    for( int i = n - 1; i >= 0; -- i ){
        if( b[i] != 0 ){
            is_b = 0;
            break;
        }
    }

    return pair<bool, COST>(is_b, cost_sum);
}
#endif