双方向ダイクストラ

・（追記）2015年11月3日：Mi_Sawaさんとtmaeharaさんのソースコードへのリンクとtmaeharaさんのプログラムとの比較

[概要]

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

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

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

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

参考

大規模グラフアルゴリズムの最先端 from Takuya Akiba

www.slideshare.net

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

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

tmaeharaさんのツイートで終了条件を工夫すると速くなるそうです．

これまでに見たことのある (sからの距離) --(枝)-- (tまでの距離) の上界 μ を覚えておいて，キューの先頭の和がこれを越えたところで終了条件にしていい．

— ᴛᴀᴋᴀɴᴏʀɪ ᴍᴀᴇʜᴀʀᴀ (@tmaehara) 2015, 11月 2

[実験]

実験内容

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

データは二次元平面上の点の座標として与えられています．都市の数は115475です．各々の都市を頂点とした完全グラフで頂点間の距離はユークリッド距離として無向グラフを作成します（実際には頂点数が多いので陽にグラフを持たない）．

＊注意　変更（2015/11/3）
二次元平面上の点の集合からグラフを作るときに完全グラフとしましたが，完全グラフだとsとtの間に辺があり三角不等式から他の頂点を経由するよりもsとtの間の辺のみを使ったほうが最短路となります．
完全グラフに変換した下の実験結果はいろいろおかしなことになっています．点の広がり具合は正しいのですがおかしいです．点の集合からグラフを作成する方法を考える必要がありそうです．

グラフの描画はOpenFrameworks（openFrameworksJp）を使っています．

実験結果

一様ランダムにsとtを選びs-t最短路を次の3つのアルゴリズムで求めます．
　(a) ダイクストラ法
　(b) 双方向ダイクストラ法
　(c) tmaeharaさんの上界を考慮した双方向ダイクストラ法

10セットに対して各々4回時間を計測しました．それぞれのtimeは4回の平均時間です．解が正しいかどうかは各々のアルゴリズムが出力する経路長が等しいことで確認をしました．
(a)と(b)のアルゴリズムのそれぞれのデータ・セットに対して最短経路が確定した頂点を色付けしています．sとtから最短経路が確定している頂点をそれぞれ赤色と青色とします．sとtは少し大きく表示しています．

(1) Find 107983-37236 shortest path = 37341.5
               Dijkstra : time = 35.75秒
 Bidirectional Dijkstra : time = 50.75秒
  Maeharasan no Program : time = 0秒

　

(2) Find 30352-64115 shortest path = 14454.2
               Dijkstra : time = 47.5秒
 Bidirectional Dijkstra : time = 27.75秒
  Maeharasan no Program : time = 0秒

　

(3) Find 100118-98660 shortest path = 2203.29
               Dijkstra : time = 1秒
 Bidirectional Dijkstra : time = 0.25秒
  Maeharasan no Program : time = 0秒

　

(4) Find 14236-84642 shortest path = 25333.2
               Dijkstra : time = 49.25秒
 Bidirectional Dijkstra : time = 59秒
  Maeharasan no Program : time = 0秒

　

(5) Find 88274-17971 shortest path = 7347.34
               Dijkstra : time = 22.75秒
 Bidirectional Dijkstra : time = 5.25秒
  Maeharasan no Program : time = 0秒

　

(6) Find 15312-9902 shortest path = 45132
               Dijkstra : time = 58.25秒
 Bidirectional Dijkstra : time = 64.5秒
  Maeharasan no Program : time = 0秒

　

(7) Find 101765-80234 shortest path = 13710.6
               Dijkstra : time = 10.25秒
 Bidirectional Dijkstra : time = 5.25秒
  Maeharasan no Program : time = 0秒

　

(8) Find 66463-50876 shortest path = 12381.7
               Dijkstra : time = 11.25秒
 Bidirectional Dijkstra : time = 9.5秒
  Maeharasan no Program : time = 0秒

　

(9) Find 35090-62401 shortest path = 19442.1
               Dijkstra : time = 51.5秒
 Bidirectional Dijkstra : time = 31秒
  Maeharasan no Program : time = 0秒

　

(10) Find 43497-89411 shortest path = 15385.5
               Dijkstra : time = 21.5秒
 Bidirectional Dijkstra : time = 27.5秒
  Maeharasan no Program : time = 0秒

　

まとめ

右側の都市たちは密なので広がるスピードが遅く，左側の都市たちは疎なので広がるスピードが速く感じます．

アメリカの都市のデータは完全グラフでtmaeharaさんのアルゴリズムでは始点と終点から隣接頂点を見るステップで終了してるので0秒となっています．なので，計算時間での正しい比較ができていないので対策を考えたいと思います．

時間があれば上の参考資料に書かれている他の実装もやってみたいです．

[ソースコード]

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

Mi_Sawaさんのソースコード： https://ideone.com/UayEyc
tmaeharaさんのソースコード: https://ideone.com/J9awus

#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;
}