読者です 読者をやめる 読者になる 読者になる

tatanaideyoの備忘録Ⅱ

競技プログラミンはtshitaでやっています

双方向ダイクストラ

[概要]

双方向ダイクストラを初めて実装したので備忘録としてまとめておきます.
時間がないのでざっくりとしか書いてません.
双方向ダイクストラ法を知らなかったのでMisawaさんに教えていただきました.有り難うございます.

無向グラフ上で始点sと終点tが与えられたときのsからtへの最短路を求める問題を考えます.この問題は2頂点対最短経路問題と呼ばれています.グラフの上で最短経路を求める問題は一般的に最短経路問題と呼ばれておりいろいろな種類があります.
最短経路問題 - Wikipedia

特にsから各頂点への最短経路を求める問題は単一始点最短経路問題と呼ばれており,その問題を解く有名なアルゴリズムダイクストラ法(ダイクストラ法 - Wikipedia)があります.
ダイクストラ法の詳しい説明はしませんがイメージとしてはsを中心に近い頂点から最短路を確定していきます.
f:id:tatanaideyo:20151031190922p:plain

単一始点最短経路問題ではsからすべての頂点への最短経路を求める問題ですが,途中で頂点tへの最短経路が確定したらアルゴリズムを終了するとして2頂点対最短経路問題としても考えることが出来ます.
sから広がっていくダイクストラ法ではsを中心に近い方向へ優先的に広がっていくのでtに到達するまでに余分な頂点を見ている気がします.
そこで,sからのみではなくtからも同じくダイクストラ法をして,お互いからの最短路が確定した頂点が見つかったときに終了すると余分な頂点が少なく済みそうな気がします.
この双方向から始めるダイクストラ法を双方向ダイクストラ法(Bidirectional search - Wikipedia, the free encyclopedia)と呼びます.注意が必要なのですがダイクストラ法と双方向ダイクストラ法の最悪計算量は同じです.下の図はイメージ図です(下の参考のp.16の図が分かりやすかったのでそのまま同じです).
f:id:tatanaideyo:20151031192716p:plain

[双方向ダイクストラ法のざっくり説明]

ダイクストラ法を知っている前提で説明します.
sとtからダイクストラ法を行います.頂点vにおいて初めてsからの最短路距離とtからの最短距離が確定されたとします.このとき,2つのダイクストラ法を終了します.
そして,各頂点vに対して(sからvへの最短距離)+(tからvへの最短距離)が最小となるものが最短距離となります.
最後の部分が必要な理由があまり分かっていないので分かったら書きます.
(下のソースコードの方が伝わるかもしれません)

[実験]

実験内容

最悪計算量は同じなので実際のインスタンで時間計測をすることによって比較することにします.用いるデータはTSP(巡回セールスマン問題 - Wikipedia)のベンチマークとして知られているアメリカの都市のデータです.
  データ・セット:USA Traveling Salesman Problem

データは二次元平面上の点の座標として与えられています.都市の数は115475です.各々の都市を頂点とした完全グラフで頂点間の距離はユークリッド距離として無向グラフを作成します(実際には頂点数が多いので陽にグラフを持たない).
グラフの描画はOpenFrameworks(openFrameworksJp)を使っています.
f:id:tatanaideyo:20151101013508p:plain

実験結果

一様ランダムにsとtを選びs-t最短路をダイクストラ法と双方向ダイクストラ法で求めます.これを10セット行い時間を計測しました.lenは求めた距離で2つのアルゴリズムが同じ値を出力するかのチェックで表示させています.
それぞれのデータ・セットに対して最短経路が確定した頂点を色付けしています.sとtから最短経路が確定している頂点をそれぞれ赤色と青色とします.sとtは少し大きく表示しています.

(1) Find 107983-37236 shortest path
               Dijkstra : len = 37341.5, time = 33秒
 Bidirectional Dijkstra : len = 37341.5, time = 43秒

f:id:tatanaideyo:20151101024812p:plain:w270 f:id:tatanaideyo:20151101024824p:plain:w270



(2) Find 30352-64115 shortest path
               Dijkstra : len = 14454.2, time = 45秒
 Bidirectional Dijkstra : len = 14454.2, time = 25秒

f:id:tatanaideyo:20151101024843p:plain:w270 f:id:tatanaideyo:20151101024858p:plain:w270



(3) Find 100118-98660 shortest path
               Dijkstra : len = 2203.29, time = 1秒
 Bidirectional Dijkstra : len = 2203.29, time = 0秒

f:id:tatanaideyo:20151101024912p:plain:w270 f:id:tatanaideyo:20151101024921p:plain:w270



(4) Find 14236-84642 shortest path
               Dijkstra : len = 25333.2, time = 48秒
 Bidirectional Dijkstra : len = 25333.2, time = 51秒

f:id:tatanaideyo:20151101024935p:plain:w270 f:id:tatanaideyo:20151101024944p:plain:w270



(5) Find 88274-17971 shortest path
               Dijkstra : len = 7347.34, time = 20秒
 Bidirectional Dijkstra : len = 7347.34, time = 4秒

f:id:tatanaideyo:20151101024957p:plain:w270 f:id:tatanaideyo:20151101025005p:plain:w270



(6) Find 15312-9902 shortest path
               Dijkstra : len = 45132, time = 53秒
 Bidirectional Dijkstra : len = 45132, time = 55秒

f:id:tatanaideyo:20151101025023p:plain:w270 f:id:tatanaideyo:20151101025030p:plain:w270



(7) Find 101765-80234 shortest path
               Dijkstra : len = 13710.6, time = 8秒
 Bidirectional Dijkstra : len = 13710.6, time = 5秒

f:id:tatanaideyo:20151101025045p:plain:w270 f:id:tatanaideyo:20151101025054p:plain:w270



(8) Find 66463-50876 shortest path
               Dijkstra : len = 12381.7, time = 10秒
 Bidirectional Dijkstra : len = 12381.7, time = 9秒

f:id:tatanaideyo:20151101025103p:plain:w270 f:id:tatanaideyo:20151101025112p:plain:w270



(9) Find 35090-62401 shortest path
               Dijkstra : len = 19442.1, time = 46秒
 Bidirectional Dijkstra : len = 19442.1, time = 27秒

f:id:tatanaideyo:20151101025123p:plain:w270 f:id:tatanaideyo:20151101025131p:plain:w270



(10) Find 43497-89411 shortest path
               Dijkstra : len = 15385.5, time = 21秒
 Bidirectional Dijkstra : len = 15385.5, time = 26秒

f:id:tatanaideyo:20151101025143p:plain:w270 f:id:tatanaideyo:20151101025153p:plain:w270

まとめ

右側の都市たちは密なので広がるスピードが遅く,左側の都市たちは疎なので広がるスピードが速く感じます.なので,sとtの頂点がどこにあるのかでどちらのアルゴリズムが良さそうなのか言えそうですが,今回のデータは平面的グラフでかつ三角不等式が成り立つので一般グラフの場合にこの感じが成り立つのかは分かりません.
時間があれば上の参考資料に書かれている他の実装もやってみたいです.

[ソースコード]

双方向ダイクストラ法の実装では2つのプライオリティキューを使っています.Misawaさんの実装でpriority_queueをvectorに入れて管理していてとてもスッキリしていたので真似ました.

#include <bits/stdc++.h>

using namespace std;

typedef pair<double, double>  P;

vector<P> p; // 二次元平面上の点の集合

// 二点間の距離を返す
inline double Weight(const pair<double, double> &a, const pair<double, double> &b)
{
    return sqrt((a.first - b.first) * (a.first - b.first) +
                (a.second - b.second) * (a.second - b.second));
}

// s-t 双方向ダイクストラ法
double BidirectionalDijkstra(int s, int t)
{
    const int n = p.size(); // 頂点数
    const double MAX_W = DBL_MAX / 2;

    // d[0][v] := sからvへの最短距離, d[1][v] := tからvへの最短距離
    vector<vector<double>> d(2, vector<double>(n, MAX_W));
    // visit[0][v] := sからvへ更新終了, visit[1][v] := tからvへ更新終了
    vector<vector<bool>> visit(2, vector<bool>(n, false));
    // pq[0] := sからの更新情報を管理,pq[1] := tからの更新情報を管理
    vector<priority_queue<P, vector<P>, greater<P>>> pq(2);

    pq[0].emplace(d[0][s] = 0, s);
    pq[1].emplace(d[1][t] = 0, t);

    bool loop_continue = true;
    while (loop_continue && !pq[0].empty() && !pq[1].empty()) {
        for (int i = 0; i < 2; ++i) {
            if (pq[i].empty())
                continue;

            int u;
            double w;

            tie(w, u) = pq[i].top();  pq[i].pop();
            if (w <= d[i][u]) {
                visit[i][u] = true;

                // 終了条件: 頂点uにおいてsとtからの最短経路が確定
                if (w != 0 && visit[0][u] && visit[1][u]) {
                    loop_continue = false;
                    break;
                }

                for (int v = 0; v < n; ++v) {
                    if (u == v)
                        continue;
                    double cost = Weight(p[u], p[v]) + w;
                    if (cost < d[i][v]) {
                        d[i][v] = cost;
                        pq[i].emplace(d[i][v] = cost, v);
                    }
                }
            }
        }
    }

    // 最短距離はすべての頂点vに対してs-v-tの最短距離の最小値
    double ans = MAX_W;
    for (int v = 0; v < n; ++v)
        ans = min(ans, d[0][v] + d[1][v]);

    return ans;
}

// s-t ダイクストラ法
double Dijkstra(int s, int t)
{
    const int n = p.size(); // 頂点数
    const double MAX_W = DBL_MAX / 2;

    // d[0][v] := sからvへの最短距離, d[1][v] := tからvへの最短距離
    vector<double> d(n, MAX_W);
    // sからの更新情報を管理
    priority_queue<P, vector<P>, greater<P> > que;

    d[s] = 0.0;
    que.push(P(0.0, s));

    while (!que.empty()) {
        P now = que.top();  que.pop();
        int u = now.second;

        if (now.first <= d[u]) {
            // 終了条件
            if (u == t)
                return d[u];

            for (int v = 0; v < n; ++v) {
                if (u == v)
                    continue;

                double cost = Weight(p[u], p[v]) + now.first;
                if (MAX_W <= d[v] || cost < d[v]) {
                    d[v] = cost;
                    que.push(P(cost, v));
                }
            }
        }
    }

    return d[t];
}


int main()
{
    int n, id;
    double x, y;

    // Input
    cin >> n;
    cout << "Vertex size = " << n << "\n\n";
    p.resize(n);

    for (int i = 0; i < n; ++i)
        cin >> id >> p[i].first >> p[i].second;

    // Output
    for (int cnt = 0; cnt < 10; ++cnt) {
        int s = rand() % n, t = rand() % n;

        cout << "(" << cnt + 1 << ") Find " << s << "-" << t << " shortest path\n";

        // Dijkstra
        auto start = std::chrono::high_resolution_clock::now();
        cout << "               Dijkstra : len = " << Dijkstra(s, t);
        auto end = std::chrono::high_resolution_clock::now();
        auto take_time = std::chrono::duration_cast<std::chrono::seconds>(end - start);
        cout << ", time = " << take_time.count() << "秒\n";

        // Bidirectional Dijkstra
        start = std::chrono::high_resolution_clock::now();
        cout << " Bidirectional Dijkstra : len = " << BidirectionalDijkstra(s, t);
        end = std::chrono::high_resolution_clock::now();
        take_time = std::chrono::duration_cast<std::chrono::seconds>(end - start);
        cout << ", time = " << take_time.count() << "秒\n\n";
    }

    return 0;
}