2019-03-20

AOJ 2171 Strange Couple (JAG 夏合宿 2009 day3-G) (600 点)

AOJ JAG AOJ-ICPC600点期待値ランダムウォーク期待値DP DP 連立一次方程式行列方程式最適化テク:解に確実に含まれる要素を列挙するグラフ最短路問題 Dijkstra法前処理自己ループを除く AOJ-ICPC JAG夏合宿そのまま覚えたいシンプル設定の中堅以上の典型問題回数の期待値 Gauss-Jordanの掃き出し法浮動小数点型を扱う問題

ランダムウォークな問題！
実数係数の連立方程式の練習に

問題へのリンク

問題概要

$N$ 頂点の重み付き無向グラフが与えられる。ただし自己ループを含み得る (多重辺はない)。二頂点 $s, t$ が指定されていて、 $s$ から $t$ へとランダムウォークによって辿り着きたい。

ただし $N$ 頂点にはそれぞれ「標識あり」「標識なし」の属性があり、

標識あり：今いる頂点から出ている辺のうち、 $t$ への最短路としてありうるものすべてを考えた時、使う可能性のある辺のみからランダムに一つ辺を選んで進む
標識なし：今いる頂点から出ている辺をランダムに選ぶ

となっている。 $s$ から $t$ へとたどり着くまでに要するコストの期待値を求めよ。ただしたどり着くことが不可能な場合は "impossible" と出力せよ。

制約

$1 \le N \le 100$

考えたこと

少し設定が複雑とは言え、ランダムウォークなので「蟻本上級編」に載っている通り、連立一次方程式で解くことができる。

「標識あり」の場合について、最短路に使われる可能性のある辺をすべて列挙する処理が必要になる。それは典型で、Dijkstra 法の復元をする気持ちで実装することができる。

頂点 v から出ている辺のうち、最短路として使われ得るものを挙げる

というのを効率良く実装するためには、Dijkstra 法を t 側から回しておくと便利である。さてそのようにすると結局、標識の有無を問わず

各頂点 $v$ について、次に行く可能性のある頂点の集合 $w_1, w_2, \dots, w_K$ が列挙されてそれらに等確率で行く

という状態になる。期待値 DP の式を立てると

$E_v$ = 頂点 $v$ から $t$ へと至るまでのコストの期待値

として

$E_t = 0$
$E_v = \frac{\sum_{w} (E_{w} + {\rm cost}(v, w))}{K}$ 　( $v \neq t$ のとき)

となる。「標識なし」の場合についてはさらに「自己ループ」がある場合はそれを除く必要がる。すなわち左辺にも右辺にも $E_v$ が登場しているので、それを式変形してあげる必要がある。式変形すると

$E_{v} = \frac{\sum_{w \neq v} (E_{w} + {\rm cost}(v, w)) + {\rm cost}(v, v)}{K-1}$

という風になる。あとはこれらを連立一次方程式として立てればよい。まとめると以下のようになる。

v = t の場合:

$E_{v} = 0$

それ以外で自己ループなしの場合:

$KE_{v} - \sum_{w} E_{w} = \sum_{w} {\rm cost}(v, w)$

それ以外で自己ループありの場合:

$(K - 1)E_{v} - \sum_{w \neq v} E_{w} = \sum_{w} {\rm cost}(v, w)$

計算量は $O(N^{3})$ 。

#include <iostream>
#include <vector>
#include <queue>
#include <cmath>
#include <iomanip>
#include <algorithm>
using namespace std;
#define COUT(x) cout << #x << " = " << (x) << " (L" << __LINE__ << ")" << endl


using D = double;
const D EPS = 1e-10;

template<class T> struct Matrix {
    vector<vector<T> > val;
    Matrix(int n, int m, T x = 0) : val(n, vector<T>(m, x)) {}
    void init(int n, int m, T x = 0) {val.assign(n, vector<T>(m, x));}
    size_t size() const {return val.size();}
    inline vector<T>& operator [] (int i) {return val[i];}
};

template<class T> ostream& operator << (ostream& s, Matrix<T> A) {
    s << endl; 
    for (int i = 0; i < A.size(); ++i) {
        for (int j = 0; j < A[i].size(); ++j) {
            s << A[i][j] << ", ";
        }
        s << endl;
    }
    return s;
}

template<class T> Matrix<T> operator * (Matrix<T> A, Matrix<T> B) {
    Matrix<T> R(A.size(), B[0].size());
    for (int i = 0; i < A.size(); ++i)
        for (int j = 0; j < B[0].size(); ++j)
            for (int k = 0; k < B.size(); ++k)
                R[i][j] += A[i][k] * B[k][j];
    return R;
}

template<class T> Matrix<T> pow(Matrix<T> A, long long n) {
    Matrix<T> R(A.size(), A.size());
    for (int i = 0; i < A.size(); ++i) R[i][i] = 1;
    while (n > 0) {
        if (n & 1) R = R * A;
        A = A * A;
        n >>= 1;
    }
    return R;
}

template<class T> int GaussJordan(Matrix<T> &A, bool is_extended = false) {
    int m = A.size(), n = A[0].size();
    int rank = 0;
    for (int col = 0; col < n; ++col) {
        if (is_extended && col == n-1) break;
        int pivot = -1;
        T ma = EPS;
        for (int row = rank; row < m; ++row) {
            if (abs(A[row][col]) > ma) {
                ma = abs(A[row][col]);
                pivot = row;
            }
        }
        if (pivot == -1) continue;
        swap(A[pivot], A[rank]);
        auto fac = A[rank][col];
        for (int col2 = 0; col2 < n; ++col2) A[rank][col2] /= fac;
        for (int row = 0; row < m; ++row) {
            if (row != rank && abs(A[row][col]) > EPS) {
                auto fac = A[row][col];
                for (int col2 = 0; col2 < n; ++col2) {
                    A[row][col2] -= A[rank][col2] * fac;
                }
            }
        }
        ++rank;
    }
    return rank;
}

template<class T> vector<T> linear_equation(Matrix<T> A, vector<T> b) {
    // extended
    int m = A.size(), n = A[0].size();
    Matrix<T> M(m, n + 1);
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < n; ++j) M[i][j] = A[i][j];
        M[i][n] = b[i];
    }
    int rank = GaussJordan(M, true);
    
    // check if it has no solution
    vector<T> res;
    for (int row = rank; row < m; ++row) if (abs(M[row][n]) > EPS) return res;

    // answer
    res.assign(n, 0);
    for (int i = 0; i < rank; ++i) res[i] = M[i][n];
    return res;
}



int N, s, t;
vector<int> q;
vector<vector<int> > a;

void solve() {
    // Dijkstra
    const int INF = 1<<29;
    vector<int> dist(N, INF);
    vector<bool> seen(N, 0);
    dist[t] = 0;
    for (int iter = 0; iter < N; ++iter) {
        int curd = INF;
        int v = -1;
        for (int i = 0; i < N; ++i) {
            if (seen[i]) continue;
            if (curd > dist[i]) {
                curd = dist[i];
                v = i;
            }
        }
        if (v == -1) break;
        for (int w = 0; w < N; ++w) {
            if (w == v) continue;
            if (a[v][w] == 0) continue;
            dist[w] = min(dist[w], curd + a[v][w]);
        }
        seen[v] = true;
    }
    if (dist[s] >= INF) {
        cout << "impossible" << endl;
        return;
    }

    // 連立一次方程式
    Matrix<D> A(N, N, 0); vector<D> b(N, 0);
    for (int v = 0; v < N; ++v) {
        if (v == t) {
            A[v][v] = 1;
            b[v] = 0;
        }
        else {
            vector<int> neigbor;
            for (int w = 0; w < N; ++w) {
                if (a[v][w] == 0) continue;
                if (q[v] == 1 && dist[w] + a[v][w] != dist[v]) continue;
                neigbor.push_back(w);
            }
            int K = neigbor.size();
            for (auto w : neigbor) {
                A[v][w] -= 1;
                b[v] += a[v][w];
            }
            A[v][v] += K;
        }
    }
    auto res = linear_equation(A, b); 
    if (res.empty()) cout << "impossible" << endl;
    else cout << fixed << setprecision(15) << res[s] << endl;
}

int main() {
    while (cin >> N >> s >> t, N) {
        --s, --t;
        q.resize(N);
        for (int i = 0; i < N; ++i) cin >> q[i];
        a.assign(N, vector<int>(N));
        for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) cin >> a[i][j];
        solve();
    }
}

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

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

AOJ 2171 Strange Couple (JAG 夏合宿 2009 day3-G) (600 点)

問題概要

制約

考えたこと