二段階単体法のデバッグに苦労しました。
問題概要
初期状態が全要素が 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; } };