maxflow_3.cc

#include <bits/stdc++.h>
using namespace std;
#define rep(i, n) for(int i = 0; i < (n); ++i)
#define REP(i, b, n) for(int i = (b); i < (n); ++i)
#define let(v, x) __typeof(x) v = (x)
#define foreach(i,v) for(let(i, (v).begin());i!=(v).end();i++)

/**
 * 割とどうしてくれてもいいけど, 理解 && バグ潰ししてから使うべきです.
 *
 */

/**
 * Dinic と同じように,
 * - 頂点数を取るコンストラクタ
 * - add_edge(src, dst, cap)
 * - run(s, t)
 * の 3 つを備えるクラス (例えば MyDinic) を作って,
 *
 * Benchmark<Dinic> bm
 * を
 * Benchmark<Dinic, MyDinic> bm
 * とかに変えれば, 比較出来る.
 */

struct DinicLC{//{{{
    typedef int Cap;
    static const Cap INF = 1<<29;
    struct LCNode{//{{{
        typedef LCNode* Ptr;
        Ptr ch[2], par, min_node;
        Cap val, min_val, lazy;
        int name;
        LCNode() : par(0), min_node(this), val(INF), min_val(INF), lazy(0) {
            ch[0] = ch[1] = 0;
        }
        inline void push(){//{{{
            if(lazy == 0) return;
            val += lazy; min_val += lazy;
            rep(i, 2) if(ch[i]) ch[i]->lazy += lazy;
            lazy = 0;
        }//}}}
        inline void update(){//{{{
            push();
            if(ch[0]){
                min_val = ch[0]->min_val + ch[0]->lazy;
                min_node = ch[0]->min_node;
            }
            if(!ch[0] || min_val > val){
                min_val = val; min_node = this;
            }
            if(ch[1] && min_val > ch[1]->min_val + ch[1]->lazy){
                min_val = ch[1]->min_val + ch[1]->lazy;
                min_node = ch[1]->min_node;
            }
        }//}}}
        inline void set_value(Cap v){ expose(); val = v; update(); }
        inline Cap get_value(){ expose(); return val; }
        inline Ptr get_min_node_on_path(){ expose(); return min_node; }
        inline void add_to_path(Cap c){ expose(); lazy += c; push(); }
        inline int state(){//{{{
            if(par && par->ch[0] == this) return -1;
            if(par && par->ch[1] == this) return +1;
            return 0;
        }//}}}
        inline void rotate(){//{{{
            Ptr p = par; p->push(); push();
            bool lr = state()+1;
            if(p->state()) p->par->ch[p->state()>0] = this;
            par = p->par;
            p->ch[lr] = ch[!lr];
            if(ch[!lr]) ch[!lr]->par = p;
            (p->par = this)->ch[!lr] = p;
            p->update(); update();
        }//}}}
        inline void splay(){//{{{
            while(state()){
                int s = state() * par->state();
                if(s) (s == 1 ? par : this)->rotate();
                rotate();
            }
            push();
        }//}}}
        inline void expose(){//{{{
            if(!par) return;
            for(Ptr p = this; p; p = p->par) p->splay();
            for(Ptr p = this; p->par; p = p->par) p->par->ch[1] = p;
            splay();
        }//}}}
        inline void cut(){//{{{
            expose();
            if(!ch[0]) return;
            ch[0]->par = 0;
            ch[0]->push();
            ch[0] = 0;
            update();
        }//}}}
        inline void link_to(Ptr p){//{{{
            expose();
            par = p;
        }//}}}
        inline Ptr find_root(){//{{{
            expose();
            Ptr p = this;
            while(p->ch[0]) p = p->ch[0];
            p->expose();
            return p;
        }//}}}
    };//}}}

    struct E{//{{{
        int dst;
        Cap cap;
        int rev;
        E(int dst, Cap cap, int rev) : dst(dst), cap(cap), rev(rev) {}
    };//}}}
    int n;
    vector<vector<E> > g;
    DinicLC(int n) : n(n), g(n) {}
    inline void add_edge(int src, int dst, Cap cap){//{{{
        if(src == dst) return;
        g[src].push_back(E(dst,cap,g[dst].size()));
        g[dst].push_back(E(src, 0 ,g[src].size() - 1));
    }//}}}
    inline void add_undirected_edge(int src, int dst, Cap cap){//{{{
        if(src == dst) return;
        g[src].push_back(E(dst,cap,g[dst].size()));
        g[dst].push_back(E(src,cap,g[src].size() - 1));
    }//}}}

    vector<int> level, active;

    vector<LCNode> lc;
    inline void unuse_edge(E &e){//{{{
        E &ee = g[e.dst][e.rev]; int u = ee.dst;
        Cap use = ee.cap - lc[u].get_value();
        e.cap += use; ee.cap -= use;
        active[u] = false;
        lc[u].cut(); lc[u].set_value(INF);
    }//}}}
    inline void init(){//{{{
        lc.resize(n);
        active.assign(n, false);
        rep(u, n) lc[u].name = u;
    }//}}}
    Cap dfs(const int &s, const int &t){//{{{
        vector<int> iter(n); rep(u, n) iter[u] = (int)(g[u].size()) - 1;
        Cap res = 0;

        while(true){
            const int &u = lc[t].find_root()->name;
            if(u == s){
                Cap f = lc[t].get_min_node_on_path()->get_value();
                res += f;
                lc[t].add_to_path(-f);
                while(true){
                    int v = lc[t].get_min_node_on_path()->name;
                    Cap f = lc[v].get_value();
                    if(f > 0) break;
                    unuse_edge(g[v][iter[v]]);
                }
            }else{
                for(int &i = iter[u]; i >= 0; --i){
                    E &e = g[u][i], &ee = g[e.dst][e.rev]; int v = e.dst;
                    if(level[u]-1 != level[v] || ee.cap <= 0) continue;
                    active[u] = true;
                    lc[u].set_value(ee.cap);
                    lc[u].link_to(&lc[v]);
                    break;
                }
                if(active[u]) continue;
                if(u == t) break;
                level[u] = -1;
                foreach(e, g[u]) if(active[e->dst] && iter[e->dst] == e->rev)
                    unuse_edge(g[e->dst][e->rev]);
            }
        }
        rep(u, n) if(active[u]) unuse_edge(g[u][iter[u]]);
        return res;
    }//}}}

    Cap run(int s, int t){//{{{
        init();
        vector<int> q(n);
        for(Cap flow = 0; ;){
            level.assign(n, -1);
            int ql = 0, qr = 0;
            level[q[qr++] = s] = 0;
            while(ql != qr && level[t] == -1){
                int u = q[ql++];
                foreach(e, g[u]) if(e->cap > 0 && level[e->dst] == -1)
                    level[ q[qr++] = e->dst ] = level[u] + 1;
            }
            if(level[t] == -1) return flow;
            flow += dfs(s, t);
        }
    }//}}}
};//}}}

struct Wave{//{{{
    typedef int Cap;
    static const Cap INF = 1<<29;
    struct E{//{{{
        int src, dst;
        Cap cap, _cap;
        int rev;
        E(int src, int dst, Cap cap, int rev) :
            src(src), dst(dst), cap(cap), _cap(0), rev(rev) {}
    };//}}}
    int n;
    vector<vector<E> > g;
    Wave(int n) : n(n), g(n) {}
    inline void add_edge(int src, int dst, Cap cap){//{{{
        if(src == dst) return;
        g[src].push_back(E(src, dst,cap,g[dst].size()));
        g[dst].push_back(E(dst, src, 0 ,g[src].size() - 1));
    }//}}}

    inline E &rev(const E& e){ return g[e.dst][e.rev]; }
    vector<int> level, iter;
    vector<int> blocked;
    vector<int> unbalance[2]; // { unblocked, blocked }
    vector<int> inS;
    vector<Cap> ex;
    inline void push(E &e){//{{{
        Cap f = min(e.cap, ex[e.src]);
        if(!inS[e.dst]) unbalance[blocked[e.dst]].push_back(e.dst);
        inS[e.dst] = true;
        e.cap -= f; rev(e).cap += f;
        ex[e.src] -= f; ex[e.dst] += f;
    }//}}}
    // dir = +1 ? s to t : t to s
    inline void discharge(const int& u, const int dir){//{{{
        for(int &i = iter[u]; i < g[u].size(); ++i){
            E &e = g[u][i];
            if(e.cap == 0 || level[e.src] + dir != level[e.dst]) continue;
            if(dir == +1 && blocked[e.dst]) continue;
            push(e);
            if(ex[u] == 0) return;
        }
        blocked[u] = true;
        if(!inS[u]) unbalance[1].push_back(u);
        inS[u] = true;
        iter[u] = 0;
    }//}}}
    Cap run(int s, int t){//{{{
        vector<int> q(n);
        vector<int> tmp; tmp.reserve(n);
        rep(i, 2) unbalance[i].reserve(n);
        rep(i, 2) unbalance[i].clear();
        inS.assign(n, false); inS[s] = inS[t] = true;
        for(Cap flow = 0; ;){
            level.assign(n, -1);
            int ql = 0, qr = 0;
            level[q[qr++] = s] = 0;
            while(ql != qr && level[t] == -1){
                const int &u = q[ql++];
                foreach(e, g[u]) if(level[e->dst] == -1 && e->cap + e->_cap > 0)
                    level[ q[qr++] = e->dst ] = level[u] + 1;
            }
            if(level[t] == -1) return flow;

            rep(i, qr){
                if(level[q[i]] == level[t] && q[i] != t){
                    level[q[i]] = -1; continue;
                }
                foreach(e, g[q[i]]){
                    if(level[e->src] + 1 == level[e->dst]){
                        e->cap += e->_cap; e->_cap = 0;
                    }else{
                        e->_cap += e->cap; e->cap = 0;
                    }
                }
            }

            iter.assign(n, 0);
            blocked.assign(n, false);
            ex.assign(n, 0); ex[s] = INF; discharge(s, +1); ex[s] = 0;

            while(!unbalance[0].empty() || !unbalance[1].empty()){
                rep(b, 2) while(!unbalance[b].empty()){
                    tmp.clear();
                    int l = level[unbalance[b].back()];
                    while(!unbalance[b].empty() && level[unbalance[b].back()] == l){
                        int v = unbalance[b].back(); unbalance[b].pop_back();
                        inS[v] = false; tmp.push_back(v);
                    }
                    foreach(v, tmp) discharge(*v, b ? -1 : +1);
                }
            }
            flow += ex[t];
        }
    }//}}}
};//}}}

struct Dinic{//{{{
    typedef int Cap;
    static const Cap INF = 1<<29;
    struct E{//{{{
        int dst;
        Cap cap;
        int rev;
        E(int dst, Cap cap, int rev) : dst(dst), cap(cap), rev(rev) {}
    };//}}}
    int n;
    vector<vector<E> > g;
    Dinic(int n) : n(n), g(n) {}
    inline void add_edge(int src, int dst, Cap cap){//{{{
        if(src == dst) return;
        g[src].push_back(E(dst,cap,g[dst].size()));
        g[dst].push_back(E(src, 0 ,g[src].size() - 1));
    }//}}}
    inline void add_undirected_edge(int src, int dst, Cap cap){//{{{
        if(src == dst) return;
        g[src].push_back(E(dst,cap,g[dst].size()));
        g[dst].push_back(E(src,cap,g[src].size() - 1));
    }//}}}

    vector<int> level, iter;
    Cap dfs(const int &s, const int &u, Cap crr){//{{{
        if(s == u || crr == 0) return crr;
        Cap sum = 0;
        for(int &i = iter[u]; i < g[u].size(); ++i){
            E &e = g[u][i], &ee = g[e.dst][e.rev];
            const int &v = e.dst; // s -- v - u -- t
            if(level[v] >= level[u] || ee.cap <= 0) continue;
            Cap f = dfs(s, v, min(crr - sum, ee.cap));
            if(f <= 0) continue;
            ee.cap -= f; e.cap += f;
            sum += f;
            if(sum == crr) break;
        }
        return sum;
    }//}}}
    Cap run(int s, int t){//{{{
        vector<int> q(n);
        for(Cap flow = 0; ;){
            level.assign(n, -1);
            int ql = 0, qr = 0;
            level[q[qr++] = s] = 0;
            while(ql != qr && level[t] == -1){
                int u = q[ql++];
                foreach(e, g[u]) if(e->cap > 0 && level[e->dst] == -1)
                    level[ q[qr++] = e->dst ] = level[u] + 1;
            }
            if(level[t] == -1) return flow;
            iter.assign(n, 0);
            flow += dfs(s, t, INF);
        }
    }//}}}
};//}}}

unsigned int xor128(){//{{{
    static unsigned int x=123456789,y=362436069,z=521288629,w=88675123;
    //static unsigned int x=129042, y=38923212, z=3829312, w=3893189;
    unsigned int t=x^(x<<11);
    x=y;y=z;z=w;
    return w=(w^(w>>19))^(t^(t>>8));
}//}}}
int rnd(int m){//{{{
    return xor128() % m;
}//}}}
struct Input{//{{{
    typedef int Cap;
    static const Cap C = 1<<10; // max capacity

    struct E{//{{{
        int src, dst;
        Cap cap;
        E(){}
        E(int src, int dst, Cap cap) : src(src), dst(dst), cap(cap) {}
    };//}}}
    vector<E> edges;
    int n;
    int s, t;
    Input(){}
    Input(int n, int s, int t) : n(n), s(s), t(t) { }
    inline void add_edge(int src, int dst, Cap cap){ add_edge(E(src, dst, cap)); }
    inline void add_edge(E e){ edges.push_back(e); }
    inline void reset(){ n = s = t = 0; edges.clear(); }

    friend istream& operator>>(istream &is, Input& in){//{{{
        in.reset();
        int m;
        is >> in.n >> m;
        is >> in.s >> in.t;
        rep(i, m){
            int src, dst; Cap cap;
            is >> src >> dst >> cap;
            in.add_edge(src, dst, cap);
        }
        return is;
    }//}}}
    friend ostream& operator<<(ostream &os, const Input& in){//{{{
        os << in.n << " " << in.edges.size() << endl;
        os << in.s << " " << in.t << endl;
        foreach(it, in.edges)
            os << it->src << " " << it->dst << " " << it->cap << endl;
        return os;
    }//}}}

    static Input random(int n, int m){//{{{
        Input res(n, 0, n-1);
        rep(i, m){
            E e;
            e.src = rnd(n); e.dst = rnd(n); e.cap = rnd(C);
            while(e.src == e.dst) e.dst = rnd(n);
            res.add_edge(e);
        }
        return res;
    }//}}}

    static Input random_dense(int n){//{{{
        Input res(n, 0, n-1);
        rep(i, n) rep(j, n) if(i != j) if(rnd(2)) res.add_edge(i, j, rnd(C));
        return res;
    }//}}}
    static Input random_complete(int n){//{{{
        Input res(n, 0, n-1);
        rep(i, n) rep(j, n) if(i != j) res.add_edge(i, j, rnd(C));
        return res;
    }//}}}

    // http://deepblue.lib.umich.edu/handle/2027.42/29530
    // k+1 頂点, 2k-1 辺
    static Input worst_dinic1(int k){//{{{
        Input res(k+1, 0, k);
        rep(i, k) res.add_edge(i, k, 1);
        rep(i, k-1) res.add_edge(i, i+1, k+1);
        return res;
    }//}}}

    // http://community.topcoder.com/tc?module=Static&d1=tutorials&d2=maxFlowRevisited
    // 6k-2 頂点, 3k^2 + 4k - 4 辺
    static Input worst_dinic2(int k){//{{{
        const Cap INF = 1<<25;
        const int q = k-1, n = 6*k - 2;
        const int src = 0, dst = 1;
        int a = 2;
        const int sb = a, se = a+=k;
        const int tb = a, te = a+=k;
        const int ub = a, ue = a+=2*q;
        const int vb = a, ve = a+=2*q;
        Input res(n, src, dst);
        REP(s, sb, se) res.add_edge(src, s, k);
        REP(t, tb, te) res.add_edge(t, dst, k);
        REP(s, sb, se) REP(t, tb, te) res.add_edge(s, t, 1);
        res.add_edge(src, ub, INF); res.add_edge(vb, dst, INF);
        REP(u, ub, ue-1) res.add_edge(u, u+1, INF);
        REP(v, vb, ve-1) res.add_edge(v+1, v, INF);
        for(int u = ub+1; u < ue; u += 4) REP(t, tb, te)
            res.add_edge(u, t, k);
        for(int u = ub+3; u < ue; u += 4) REP(s, sb, se)
            res.add_edge(u, s, k);
        for(int v = vb+1; v < ve; v += 4) REP(s, sb, se)
            res.add_edge(s, v, k);
        for(int v = vb+3; v < ve; v += 4) REP(t, tb, te)
            res.add_edge(t, v, k);
        return res;
    }//}}}
};//}}}

struct Timer{//{{{
    clock_t t;
    Timer() {reset();}
    void reset(){t = clock();}
    double getTime(){return ((double)clock()-(double)t)/(double)CLOCKS_PER_SEC;}
    void print(){ printf("runtime: %.3f seconds\n", getTime()); }
};//}}}

namespace U{//{{{
    template<typename T> string to_s(T t){ //{{{
        stringstream ss;
        ss << t;
        return ss.str();
    } //}}}

#include<cxxabi.h>
    template<typename T> string classname(){//{{{
        char *p;
        int status = 0;
        p = abi::__cxa_demangle(typeid(T).name(), 0, 0, &status);
        string res = p;
        free(p);
        return res;
    }//}}}
};//}}}

template<typename ... Solvers> struct Benchmark{//{{{
    typedef int Cap;
    map<string, double> time;
    template<typename T> Cap _run(const Input& in){//{{{
        T solver(in.n);
        foreach(e, in.edges) solver.add_edge(e->src, e->dst, e->cap);
        Timer timer;
        Cap res = solver.run(in.s, in.t);
        time[U::classname<T>()] += timer.getTime();
        return res;
    }//}}}
    template<typename T0, typename T1, typename ... T> Cap _run(const Input& in){//{{{
        double a = _run<T0>(in);
        double b = _run<T1, T...>(in);
        assert(a == b);
        return a;
    }//}}}
    void run(const Input& in){ _run<Solvers...>(in); }
    void reset(){ time.clear(); }
    void print(){//{{{
        cout.setf(ios::fixed);
        cout.precision(3);
        foreach(it, time){
            cout << left << setw(30) << it->first << ": ";
            cout << left << setw(6) << it->second << " sec" << endl;
        }
    }//}}}
};//}}}

/*
|* ケース                     |* 頂点数 |* 辺数             |
| ランダムを80ケース(合計)    | 5000    | 1000000           |
| ランダムを40ケース(合計)    | 50000   | 1000000           |
| ランダムを13ケース(合計)    | 500000  | 1000000           |
| 密なランダムを5ケース(合計) | 5000    | 約 [tex:n(n-1)/2] |
| 完全グラフを15ケース(合計)  | 2000    | [tex:n(n-1)]      |
| ワーストケース1つ目         | 7000    | 13997             |
| ワーストケース2つ目         | 628     | 33491             |
*/


void test(int c){
    Benchmark<Dinic, DinicLC, Wave> bm;
    //Benchmark<Dinic, DinicLC> bm;
    Input in;
    int T;

    if(c == 0 || c == -1){ // ランダムを80ケース(合計) n=5000, m=1000000//{{{
        bm.reset();
        T = 80;
        rep(i, T){
            Input in = Input::random(5000, 1000000);
            cout << "case " << i << " (n = " << in.n << ", m = " << in.edges.size() << ")" << endl;
            bm.run(in);
        }
        bm.print();
    }//}}}
    if(c == 1 || c == -1){ // ランダムを40ケース(合計) n=50000, m=1000000//{{{
        bm.reset();
        T = 40;
        rep(i, T){
            Input in = Input::random(50000, 1000000);
            cout << "case " << i << " (n = " << in.n << ", m = " << in.edges.size() << ")" << endl;
            bm.run(in);
        }
        bm.print();
    }//}}}
    if(c == 2 || c == -1){ // ランダムを13ケース(合計) n=500000, m=1000000//{{{
        bm.reset();
        T = 13;
        rep(i, T){
            Input in = Input::random(500000, 1000000);
            cout << "case " << i << " (n = " << in.n << ", m = " << in.edges.size() << ")" << endl;
            bm.run(in);
        }
        bm.print();
    }//}}}
    if(c == 3 || c == -1){ // 密なランダムを5ケース(合計) n=5000 //{{{
        bm.reset();
        T = 5;
        rep(i, T){
            Input in = Input::random_dense(5000);
            cout << "case " << i << " (n = " << in.n << ", m = " << in.edges.size() << ")" << endl;
            bm.run(in);
        }
        bm.print();
    }//}}}
    if(c == 4 || c == -1){ // 完全グラフを15ケース(合計) n=2000 //{{{
        bm.reset();
        T = 15;
        rep(i, T){
            Input in = Input::random_complete(2000);
            cout << "case " << i << " (n = " << in.n << ", m = " << in.edges.size() << ")" << endl;
            bm.run(in);
        }
        bm.print();
    }//}}}
    if(c == 5 || c == -1){ // ワーストケース1つ目 n=7000 //{{{
        bm.reset();
        T = 1;
        rep(i, T){
            Input in = Input::worst_dinic1(6999);
            cout << "case " << i << " (n = " << in.n << ", m = " << in.edges.size() << ")" << endl;
            bm.run(in);
        }
        bm.print();
    }//}}}
    if(c == 6 || c == -1){ // ワーストケース2つ目 k=105 //{{{
        bm.reset();
        T = 1;
        rep(i, T){
            Input in = Input::worst_dinic2(105);
            cout << "case " << i << " (n = " << in.n << ", m = " << in.edges.size() << ")" << endl;
            bm.run(in);
        }
        bm.print();
    }//}}}
}

int main(int argc, char *argv[]){
    if(argc > 1){
        test(atoi(argv[1]));
    }else{
        test(0);
    }
    return 0;
}

// vim:set foldmethod=marker commentstring=//%s: