けんちょんの競プロ精進記録

競プロの精進記録や小ネタを書いていきます

TopCoder SRM 694 D1H SRMDiv0Easy

二段階単体法のデバッグに苦労しました。

問題概要

初期状態が全要素が 0 であるような N 次元ベクトルに対し、
以下のような Q 個の区間クエリを実施して得られる結果が
全要素が等しい N 次元ベクトルとなるようにしたい。

区間クエリ i:
区間 [ L[i], R[i] ] に X[i] 以上 Y[i] 以下のいずれかの整数 x[i] を一律に加算する

区間クエリにおいて加算する x[i] の値を調整して、
得られる N 次元ベクトルの要素の値を最大化せよ。
そのようなものが存在しない場合は -1 をリターンせよ。

(制約)
・1 <= Q, N <= 200

解法 1 (単体法)

(i, j) 成分が L[j] <= i <= R[j] ならば 1, それ以外ならば 0 となるような行列を A,
等しい値を z として

max  z
s.t.
    X[i] <= x[i] <= Y[i]
    x[i] は整数
    Ax = [z, z, …, z]'

これは整数計画問題であるが、
A が完全単模な区間行列なので整数ベクトルな最適解が存在する。
(区間行列は完全単模であることが知られている)

あとは、これを単体法ライブラリの形 (自分のは等式標準形) に直して
ライブラリに投げる。実行時間は topcoder サーバ上で 800ms 程度。

解法 2 (フロー)

解法 1 で登場したような "区間行列" は、
差分をとることで、あるグラフの接続行列になることが知られている。
このことからフロー解が存在することが示唆される。

Ax = [z, z, …, z]'

について、各行について差分をとると、

B[i][j] =
 1 (j = L[i]) 
 -1 (j = R[i]+1)
 0 (それ以外)

となるような行列 B (A が行数 N であるのに対し、B の行数は N+1) に対して、

Bx = [z, 0, 0, …, 0, -z]'

となる。
これは、

・ノード  : (0, 1, 2, …, N)
・エッジ  : L[i] -> R[i]+1 (ここを流れるフローが x[i])
・始点と終点: 0, N
・フロー値: z

となるネットワークのフロー制約に他ならない。
容量制約は、

X[i] <= x[i] <= Y[i]

で与えられ、この容量制約の下でフロー値 z を最大化する最大流問題である。
下限容量制約つきの最大フロー問題は、
最小費用流や、二分探索 + 最大流によって解くことができる。

コード (解法1)

#include <iostream>
#include <sstream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <cstring>
#include <string>
#include <vector>
#include <stack>
#include <queue>
#include <deque>
#include <map>
#include <set>
#include <bitset>
#include <numeric>
#include <utility>
#include <iomanip>
#include <algorithm>
#include <functional>
using namespace std;
#define COUT(x) cout << #x << " = " << (x) << " (L" << __LINE__ << ")" << endl
template<class T1, class T2> ostream& operator << (ostream &s, pair<T1,T2> P)
{ return s << '<' << P.first << ", " << P.second << '>'; }
template<class T> ostream& operator << (ostream &s, vector<T> P)
{ for (int i = 0; i < P.size(); ++i) { if (i > 0) { s << " "; } s << P[i]; } return s; }

// min cx s.t. Ax = b, x >= 0
template<class T> struct Matrix {
    vector<vector<T> > val;
    Matrix(int n = 1, int m = 1) {val.clear(); val.resize(n, vector<T>(m));}
    Matrix(int n, int m, T x) {val.clear(); val.resize(n, vector<T>(m, x));}
    void init(int n, int m, T x = 0) {val.clear(); val.resize(n, vector<T>(m, x));}
    void resize(int n) {val.resize(n);}
    void resize(int n, int m, T x = 0) {val.resize(n); for (int i = 0; i < n; ++i) val[i].resize(m, x);}
    int size() {return val.size();}
    inline vector<T>& operator [] (int i) {return val[i];}
    friend ostream& operator << (ostream& s, Matrix<T> M) {s << endl;
        for (int i = 0; i < M.val.size(); ++i) s << M[i] << endl; return s;}
};

const double INF = 1LL<<29;
const double EPS = 1e-10;

template<class T> int PreProcess(Matrix<T> &A, vector<T> &b) {
    int rank = 0;
    Matrix<T> M(A.size(), A[0].size()+1);
    for (int i = 0; i < A.size(); ++i) {
        for (int j = 0; j < A[0].size(); ++j) M[i][j] = A[i][j];
        M[i][A[0].size()] = b[i];
    }

    for (int i = 0; i < A[0].size(); ++i) {
        int pivot = rank;
        T max = 0;
        for (int j = rank; j < M.size(); ++j) {
            if (max < abs(M[j][i])) {
                pivot = j;
                max = abs(M[j][i]);
            }
        }
        if (max > EPS) {
            swap(M[pivot], M[rank]);
            T div = M[rank][i];
            for (int j = 0; j < M[0].size(); ++j) {
                M[rank][j] = M[rank][j] / div;
            }
            for (int j = 0; j < M.size(); ++j) {
                if (j != rank && abs(M[j][i]) > EPS) {
                    T fac = M[j][i];
                    for (int k = 0; k < M[0].size(); ++k) {
                        M[j][k] = M[j][k] - M[rank][k] * fac;
                    }
                }
            }
            ++rank;
        }
    }
    A.resize(rank);
    b.resize(rank);
    for (int i = 0; i < rank; ++i) {
        for (int j = 0; j < A[0].size(); ++j) A[i][j] = M[i][j];
        b[i] = M[i].back();
    }
    return rank;
}

enum STATE { OPTIMIZED, INFEASIBLE, UNBOUNDED };
template<class T> T Simplex(Matrix<T> A, vector<T> b, vector<T> c, vector<T> &res, STATE &state) {
    res = vector<T>();
    PreProcess(A, b);                                               // let A be row-fullrank

    // build tableau
    int n = A.size(), m = A[0].size();
    for (int i = 0; i < n; ++i) if (b[i] < 0) {b[i] *= -1; for (int j = 0; j < m; ++j) A[i][j] *= -1; }
    vector<int> base(n), non(m);
    for (int i = 0; i < n; ++i) base[i] = m+i;
    for (int j = 0; j < m; ++j) non[j] = j;

    A.resize(n+2, n+m+1);
    for (int i = 0; i < n; ++i) A[i][m+i] = 1, A[i][n+m] = b[i];
    for (int i = 0; i < n; ++i) { A[n][n+m] += A[i][n+m]; for (int j = 0; j < m; ++j) A[n][j] += A[i][j]; }
    for (int j = 0; j < m; ++j) A[n+1][j] = -c[j];

    // start simplex
    for (int phase = 0; phase < 2; ++phase) {
        while (true) {
            int nn = -1, nb = -1;
            for (int i = 0; i < non.size(); ++i) {
                if (non[i] >= m) continue; // We cannot let slack value move to the base
                if (A[n][non[i]] > EPS) {
                    if (nn == -1) nn = i;
                    else if (non[i] < non[nn]) nn = i; // Bland's smallest subscript rule
                    //if (chmax(FMax, A[n][non[i]])) nn = i; // max coefficient rule
                }
            }
            if (nn == -1) {
                if (phase == 1) break;                              // All done!
                if (A[n][A[0].size()-1] > EPS) {                    // No feasible solution!
                    state = INFEASIBLE;
                    return -1;
                }

                // detail postprocess of phase 0
                bool ok = true;
                for (int i = 0; i < base.size(); ++i) {
                    if (base[i] >= m) {
                        ok = false; nb = i; break;                  // If base doesn't contain slack, go to phase 1
                    }
                }
                if (ok) break;
                for (int i = 0; i < non.size(); ++i) {
                    if (non[i] < m && abs(A[nb][non[i]]) > EPS) {
                        nn = i; break; // slack has base, continue phase 0
                    }
                }
                if (nn == -1) {
                    break;                                          // It can't be happened (A is row-fullrank)
                }
            }
            int col = non[nn];

            if (nb == -1) {
                T min_ratio = INF;
                for (int i = 0; i < base.size(); ++i) if (A[i][col] > EPS) {
                    T tmp = A[i][A[0].size()-1] / A[i][col];
                    if (min_ratio > tmp + EPS) { min_ratio = tmp, nb = i; }
                    else if (min_ratio >= tmp - EPS) {
                        if (nb == -1) nb = i;
                        else if (base[i] < base[nb]) nb = i;        // Bland's smallest subscript rule
                    }
                }
                if (nb == -1) {                                     // It cannot happen at the 1st stage
                    state = UNBOUNDED;
                    return -1;
                }
            }
            int row = nb;

            T piv = A[row][col];
            for (int j = 0; j < A[0].size(); ++j) A[row][j] = A[row][j] / piv;
            for (int i = 0; i < A.size(); ++i) {
                if (i == row) continue;
                T pn = A[i][col];
                for (int j = 0; j < A[0].size(); ++j) A[i][j] = A[i][j] - A[row][j] * pn;
            }
            swap(base[nb], non[nn]);
        }

        if (phase == 0) {
            swap(A[n], A[n+1]);
            A.val.pop_back();
            for (int i = 0; i < n; ++i)
                for (int j = 0; j < A.size(); ++j)
                    A[j].erase(A[j].begin() + m);

            /*
            for (int i = 0; i < base.size(); ++i) {                 // base slack can be removed
                if (base[i] >= m) {
                    A.val.erase(A.val.begin() + i);
                    int eraseID = base[i];
                    base.erase(base.begin() + i);
                    --i;
                    --n;
                }
            } */


            for (int i = 0; i < non.size(); ++i)
                if (non[i] >= m)
                    non.erase(non.begin() + i--);
        }
    }
    res = vector<T>(m, 0);
    for (int i = 0; i < base.size(); ++i) res[base[i]] = A[i].back();
    state = OPTIMIZED;
    return A[n].back();
}

class SRMDiv0Easy {
public:
    int get(int N, vector <int> L, vector <int> R, vector <int> X, vector <int> Y) {
        int Q = L.size();
        Matrix<double> A(Q*2+N, Q*3+1, 0.0);
        vector<double> b(Q*2+N, 0.0);
        vector<double> c(Q*3+1, 0.0);

        for (int i = 0; i < Q; ++i) {
            for (int j = L[i]; j <= R[i]; ++j)
                A[j][i] = -1.0;
            A[N+i][i] = 1.0;
            A[N+i][Q+i] = -1.0;
            A[N+Q+i][i] = 1.0;
            A[N+Q+i][Q*2+i] = 1.0;
            b[N+i] = X[i];
            b[N+Q+i] = Y[i];
        }
        for (int j = 0; j < N; ++j) {
            A[j][Q*3] = 1.0;
            b[j] = 0.0;
        }
        c[Q*3] = -1.0;

        vector<double> x;
        STATE state;
        double res = Simplex(A, b, c, x, state);
        if (state != OPTIMIZED) return -1;

        return (int)(-res + EPS);
    }
};

コード (解法2)

#include <iostream>
#include <sstream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <cstring>
#include <string>
#include <vector>
#include <stack>
#include <queue>
#include <deque>
#include <map>
#include <set>
#include <bitset>
#include <numeric>
#include <utility>
#include <iomanip>
#include <algorithm>
#include <functional>
using namespace std;
#define COUT(x) cout << #x << " = " << (x) << " (L" << __LINE__ << ")" << endl
template<class T> ostream& operator << (ostream &s, vector<T> P)
{ for (int i = 0; i < P.size(); ++i) { if (i > 0) { s << " "; } s << P[i]; } return s; }
template<class T> ostream& operator << (ostream &s, vector<vector<T> > P)
{ for (int i = 0; i < P.size(); ++i) { s << endl << P[i]; } return s << endl; }


const int MAX_V = 2100;
const int INF = 1LL<<29;

struct Edge {
    int rev, from, to, cap, icap;
    Edge(int r, int f, int t, int c) : rev(r), from(f), to(t), cap(c), icap(c) {}
    friend ostream& operator << (ostream& s, const Edge& E) {
        if (E.cap > 0) return s << E.from << "->" << E.to << '(' << E.cap << ')';
        else return s;
    }
};

struct Graph {
    int V;
    vector<Edge> list[MAX_V];
    
    Graph(int n = 0) : V(n) {for (int i = 0; i < MAX_V; ++i) list[i].clear();}
    void init(int n = 0) {V = n; for (int i = 0; i < MAX_V; ++i) list[i].clear();}
    void resize(int n = 0) {V = n;}
    void reset() {for (int i = 0; i < V; ++i) for (int j = 0; j < list[i].size(); ++j) list[i][j].cap = list[i][j].icap;}
    inline vector<Edge>& operator [] (int i) {return list[i];}
    
    Edge &redge(Edge e) {
        if (e.from != e.to) return list[e.to][e.rev];
        else return list[e.to][e.rev+1];
    }
    
    void addedge(int from, int to, int cap) {
        list[from].push_back(Edge(list[to].size(), from, to, cap));
        list[to].push_back(Edge(list[from].size()-1, to, from, 0));
    }
    
    void addbiedge(int from, int to, int cap) {
        list[from].push_back(Edge(list[to].size(), from, to, cap));
        list[to].push_back(Edge(list[from].size()-1, to, from, cap));
    }
    
    friend ostream& operator << (ostream& s, const Graph& G) {
        s << endl; for (int i = 0; i < G.V; ++i) {s << i << " : " << G.list[i] << endl;}return s;
    }
};

Graph G;


int level[MAX_V];
int iter[MAX_V];

void dibfs(Graph &G, int s) {
    memset(level, -1, sizeof(level));
    level[s] = 0;
    queue<int> que;
    que.push(s);
    while (!que.empty()) {
        int v = que.front();
        que.pop();
        for (int i = 0; i < G[v].size(); ++i) {
            Edge &e = G[v][i];
            if (level[e.to] < 0 && e.cap > 0) {
                level[e.to] = level[v] + 1;
                que.push(e.to);
            }
        }
    }
}

int didfs(Graph &G, int v, int t, int f) {
    if (v == t) return f;
    for(int &i = iter[v]; i < G[v].size(); ++i) {
        Edge &e = G[v][i], &re = G.redge(e);
        if (level[v] < level[e.to] && e.cap > 0) {
            int d = didfs(G, e.to, t, min(f, e.cap));
            if (d > 0) {
                e.cap -= d;
                re.cap += d;
                return d;
            }
        }
    }
    return 0;
}

int Dinic(Graph &G, int s, int t) {
    int res = 0;
    while (true) {
        dibfs(G, s);
        if (level[t] < 0) return res;
        memset(iter, 0, sizeof(iter));
        int flow;
        while ((flow = didfs(G, s, t, INF)) > 0) {
            res += flow;
        }
    }
}

int demand[1000];

class SRMDiv0Easy {
public:
    int get(int N, vector <int> L, vector <int> R, vector <int> X, vector <int> Y) {
        for (int i = 0; i < 1000; ++i) demand[i] = 0;
        for (int i = 0; i < L.size(); ++i) {
            demand[L[i]] += X[i];
            demand[R[i]+1] -= X[i];
        }
        
        int low = -1, high = INF;
        while (high - low > 1) {
            int mid = (low + high) / 2;
            demand[0] -= mid; demand[N] += mid;
            G.init(N+3);
            int s = N+1, t = N+2;
            for (int i = 0; i < L.size(); ++i) G.addedge(L[i], R[i]+1, Y[i]-X[i]);
            int flow = 0;
            for (int i = 0; i <= N; ++i) {
                if (demand[i] > 0) G.addedge(i, t, demand[i]), flow += demand[i];
                else if (demand[i] < 0) G.addedge(s, i, -demand[i]);
            }
            G.addedge(N, 0, INF);
            int res = Dinic(G, s, t);
            
            if (res < flow) high = mid;
            else low = mid;
            
            demand[0] += mid; demand[N] -= mid;
        }
        
        return low;
    }
};