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

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

AOJ 2171 Strange Couple (AOJ-ICPC 600 点)

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

問題へのリンク

問題概要

 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();
    }
}