マラソンマッチにおける典型データ構造とアプローチ
TopCoderのマラソンマッチは大きく以下の二つの形式に分類することができます。
- 最適化
- 機会学習
何れの形式のマッチもそれぞれの魅力的がありますが、今回は特に最適化に主眼を置いたマッチを題材にお話しをします。最適化を主眼においたマッチがどのようなものであるか、2013年の12月に開催されたMM 82を用いて簡単に紹介します。ColorLinkerという問題です。
プログラムには正方形のグリッドが与えられます。グリッドのいくつかのセルは予め彩色がなされています。
我々に与えられた課題は、同じ色のセルが全て繋がるように空のセルを彩色していくことです。セルを彩色するにはコストがかかります。また、同一のセルを複数色に彩色するとコストは指数的に増加します。合計のコストが低ければ低いほどよい彩色であると評価されます。同じ色のセルが連結でない場合は無効となります。
この問題に対する典型的なアプローチの一つは最短経路探索、特に複数始点から同時にダイクストラ法を適用し焼きなましを行う、というものであったようです。
さて、TopCoderの最適化問題ではこのようにシンプルな二次元グリッドの問題がほどよい頻度で出題されます。そのため、二次元グリッドを効率よく扱うことはTopCoderのマラソンマッチに参加する上で必須のテクニックであると言えます。このエントリでは、TopCoderマラソンマッチにおける典型データ構造とアプローチを主題に以下のトピックに精選し、紹介したいと思います。
- 盤面の持ちかたと最短経路探索
- キュー
- 経路復元
- コスト付き最短経路探索
盤面の持ちかたと最短経路探索
まず初めに盤面の持ちかたについて考えてみましょう。実際のマラソンマッチの問題は過分に複雑なので、これを模した問題を用意することにします。
迷路が与えられます。迷路は通路と壁で構成されています。その後3秒の間、二点がランダムに、繰り返し与えられます。この二点間の最短経路を出来るだけ速く探索するメソッドを書くのを目標とすることとします。説明やサンプルプログラムを簡単にするために、このプログラムは最短距離の長さも経路も返しません。もはや問題とはなっていませんが、ご容赦ください。
さて、目的のクラス名をsolver_tとしましょう。このクラスは以下のような二つのメソッドを持つものとします。
void solver_t::init(const std::vector<std::string>& a); void solver_t::bfs(int sx, int sy, int ex, int ey);
initメソッドにはいわゆる「迷路」が渡されます。迷路がプログラムに渡されるのは初めの1回だけとします。データは、std::stringのstd::vectorとします。パウンド(#)が壁、ピリオド(.)が通路を表します。迷路の縦と横の長さは等しくNとします。Nの範囲は毎回異なり、15〜55としましょう。迷路生成の都合から、Nは奇数とします。範囲外には外壁があるものとします。
実際の迷路のデータを覗いてみましょう。競技プログラミングではよく見る感じのデータですので、慣れているかたも多いのではないでしょうか。
.....#.#...#.....#. .###.#.#.###.#.###. .#.......#...#.#... ##.###.###.#######. ...#...#...#.#..... .###.#.#.###.#.###. .#...#.#.....#.#... ##.###.#.#.###.#.## ...#...#.#.#.#.#.#. .#####.###.#.#.###. ...#.......#.....#. .#.#####.###.#####. .#.#.#.#........... ##.#.#.###.######## .#.#...#...#....... .#####.#.#.###.#### ...#.....#.....#... .###.#####.###.###. .......#...#.......
bfsメソッドには開始点と終了点が渡されます。そこから以下のような最短経路を探索したい、ということです。
この問題を計算機で解く際はグラフの問題として捉えると便利なことがあります。その様子を少し覗いておきましょう。通路を頂点と置き換え、隣接した頂点間に辺を張ります。
さて、bfsメソッドという名前が諮詢しているように、我々の仕事は速い幅優先探索を書くことです...
では早速実装を開始します。迷路のデータはSTLを使って渡されることですし、そのままデータとして使うことにします。
クラスの定義はこのような感じになるでしょうか。
class solver_t { pubilc: void init(const std::vector<std::string>& a); void bfs(int sx, int sy, int ex, int ey); private: std::vector<std::string> a_; int N_; std::vector<std::vector<bool>> visit_; };
solver_t::init()では盤面をそのまま保持します。また、盤面の大きさに伴ったvisit_変数を確保しておきましょう。
void solver_t::init(const std::vector<std::string>& a) { a_ = a; N_ = a_.size(); visit_.assign(N_, std::vector<bool>(N_)); };
最後にbfs()メンバ関数を記述します。とてもオーソドックスな幅優先探索の実装です。
void solver_t::bfs(int sx, int sy, int ex, int ey) { for (int i = 0; i < N_; i ++) for (int j = 0; j < N_; j ++) visit_[i][j] = false; visit_[sy][sx] = true; std::queue<int> queue; queue.push(XY(sx, sy)); while (! queue.empty()) { int p = queue.front(); queue.pop(); int x = X(p); int y = Y(p); for (int k = 0; k < 4; k ++) { int xx = x + dx[k], yy = y + dy[k]; if (0 <= xx && xx < N_ && 0 <= yy && yy < N_) { if (xx == ex && yy == ey) return; if (! visit_[yy][xx]) if (a_[yy][xx] == '.') { visit_[yy][xx] = true; queue.push(XY(xx, yy)); } } } } }
座標値はintにパックするようにしました。パックには以下のマクロを用いています。与えられる迷路のサイズが最大55であるということで、55以上の一番小さなの2の累乗の数(2nの数。2、4、8、16、...等)、64を選びました。これにより、ビットシフトやビットマスクを用いて座標値のパックやアンパックをすることができます。
#define X(x) ((x) & 63) // Unpacking #define Y(x) ((x) >> 6) #define XY(x, y) (((y) << 6) + (x)) // Packing
下位16ビットの使い方は以下のようになっています。
また、dx、dyには上左下右の位置に対応する差分値が格納されています。SRMでもよく見る典型的な手法です。
const int dy[] = {-1, 0, 1, 0}; const int dx[] = { 0, -1, 0, 1};
ここで、この実装のベンチマークを取ってみましょう。initメソッド経由で55 × 55の盤面を与え、3秒の間にどれだけbfsメソッドを呼べるかを計測してみました。私の端末では111,410回呼び出すことができました(OS X 10.9.5 + clang + -O2フラグを使用)。
さて、TopCoderのマラソンマッチではアルゴリズムの善し悪しが勝負を決めます。これは決定的で揺るぎない事実です。そのため、闇雲に高速化を行うことは一般的に悪手とされています。しかし、思いついたアルゴリズムがトップ争いを征するポテンシャルを備えていたとしても、アルゴリズムを計算する速度が遅いと勝つことができません。これは勿体ない。
例えばここで説明をしている単純な幅優先探索等はその動作が速ければ速いほど有利になるかもしれません。少しばかり時間をかけ、チューンアップすることを考えてみましょう。
まずSTLのパフォーマンスを疑ってみましょう。盤面をCの二次元配列で持つように直してみます。メンバ変数の宣言と初期化部分の実装を変更するだけですので、変更箇所は少ないです。
class solver_t { private: char a_[64][64]; int N_; bool visit_[64][64]; };
void solver_t::init(const std::vector<std::string>& a) { N_ = a.size(); for (int i = 0; i < N_; i ++) for (int j = 0; j < N_; j ++) a_[i][j] = a[i][j]; }
二次元配列に置き換える際に考慮すべきこととして、各要素の要素数が上げられます。座標値をパックした際と同様に、まず2の累乗の数を検討してみるのがよいでしょう。ここでも55以上の一番小さな2の累乗の数、つまり64を選択しました。
たったこれだけの変更ですがbfsを呼び出す回数は164,900回に増えました。約1.5倍の速度改善です。STLのコンテナは大変便利ですが、処理速度を重視する必要がある場合、Cの原始的な配列に置き換えることを検討するのは良手の一つです。
ではこの実装をさらに速度改善する方法を考えてみます。一次元配列に直してみましょう。メンバ変数は以下のようになります。盤面を表すa_メンバ変数だけではなく、visit_メンバ変数にも一次元配列を採用していることに注意してください。
class solver_t { private: int N_; int LASTROW_; char a_[4096]; // 64 * 64 = 4096 bool visit_[4096]; };
LASTROW_は大変便利な変数です。盤面の範囲をチェックする際に用います。これに関しては後ほど説明をします。
盤面の初期化は以下のようになります。
void solver_t::init(const std::vector<std::string>& a) { N_ = a.size(); LASTROW_ = (N_ - 1) << 6; for (int i = 0, p = 0; i < N_; i ++, p += 64) for (int j = 0, pp = p; j < N_; j ++, pp ++) a_[pp] = a[i][j]; };
また、bfsの実装は以下のようになります。
void solver_t::bfs(int sx, int sy, int ex, int ey) { memset(visit_, 0, sizeof(visit_)); int s = XY(sx, sy), e = XY(ex, ey); std::queue<int> queue; queue.push(s); visit_[s] = true; while (! queue.empty()) { int p = queue.front(), pp; queue.pop(); int x = X(p); if (p >= 64 && ! visit_[pp = p - 64] && a_[pp] == '.') { if (pp == e) return; visit_[pp] = true; queue.push(pp); } if (x > 0 && ! visit_[pp = p - 1] && a_[pp] == '.') { if (pp == e) return; visit_[pp] = true; queue.push(pp); } if (p < LASTROW_ && ! visit_[pp = p + 64] && a_[pp] == '.') { if (pp == e) return; visit_[pp] = true; queue.push(pp); } if (x < N_ - 1 && ! visit_[pp = p + 1] && a_[pp] == '.') { if (pp == e) return; visit_[pp] = true; queue.push(pp); } } }
何かすごいことになってしまった感もありますが、これまでと同様に3秒間の呼び出し回数を計測してみると282,420回となりました。Cの二次元配列を用いていたときから2倍弱の速度改善で、STLに比べると3倍弱の速度改善、といったところでしょうか。
何故このような速度改善が達成できるかは後ほど考えてみるとして、まず実装を見て行きましょう。whileによる反復部分です。pはキューから取得した現在の位置です。この値はx座標値とy座標値をintにパックしたものでした。
ここから、以下の処理を上左下右のそれぞれの方向に対して適用しています。
- 次の位置に進めるかを判定し、
- ppに次の位置を保持しつつすでに訪れているか否かを判定し、
- 次の位置が壁でなければ、
- 目的の位置であれば、探索を終了する
- 目的の位置でなければ、次の位置をキューに積む。この際、次の頂点に訪れることを記録する
「キューに積む際に次の頂点に訪れることを記録する」のはとても大事なので特に注意しましょう。これにより、同じ頂点が再びキューに積まれることを抑制することができます。
さて、位置pから上の方向に訪れることができるかどうかは、位置pが2行目以降にあればよいことになります。これは、pと64を比較すれば判定できます。
では位置pから下の方向に訪れることができるかどうかはどのように判定すればよいでしょうか? 予め計算しておいたLASTROW_をここに用います。LASTROW_はこのように計算していますので、
LASTROW_ = (N_ - 1) << 6;
以下の図に示す位置を表すことになります。そのため、現在の位置pがLASTROW_より小さければ、下方向に探索を進めても盤面からはみ出ないことが保証されるのです。
余談ですが、この変数の使い方はTopCoderマラソンマッチでは必ず上位に食い込んでいるwleiteさんの実装を眺めていて見つけました。とてもスマートな方法ですよね。
水平方向に関しては現在のx座標値をデコードしたほうが判定がしやすいのでそのようにします。これをまとめると以下となります。
さて、何故このような速度改善が達成できたのでしょう? メモリ構造が異なるのでしょうか? いえ、Cの二次元配列と一次元配列は何れも連続して確保されるはずです。例えば、char a[5][10]として配列を確保した場合とchar a[50]として配列を確保した場合、アドレッシングの差はあれど、メモリの確保のされかたは等しくなるはずです(大きな図で見たい場合は画像をクリックしてください)。
そのため、アドレッシングをする際にかけ算かビットシフトが入ることが性能低下の原因であると考えるのが自然です。
ここで少しだけアセンブラを眺めてみましょう。64 × 64の二次元グリッドにランダムに数値を設定し、ランダムにピックアップした頂点の値のxorを取る関数を書いてみます。繰り返し数は10000回です。
int rand64(); // ランダムに[0, 64)を返す関数が外部定義してあるものとする int rand4096(); // ランダムに[0, 4096)を返す関数が外部定義してあるものとする static int a[64][64]; int test_a() { int x = 0; for (int i = 0; i < 10000; i ++) x ^= a[rand64()][rand64()]; return x; } int test_b() { int* p = reinterpret_cast<int*>(&a[0][0]); int x = 0; for (int i = 0; i < 10000; i ++) x ^= p[rand4096()]; return x; }
二次元配列を使ったtest_a関数のアセンブラ(の主要な部分)は以下のようになります。
xorl %r15d, %r15d # %r15dレジスタを0にする movl $10000, %ebx # %ebxレジスタを10000にする leaq __ZL1a(%rip), %r12 # 二次元配列aの先頭のアドレスを%r12レジスタに格納する LBB1_1: callq __Z2rand64v # rand64関数をコールする movl %eax, %r14d # 結果を%r14dレジスタに格納し、 andl $63, %r14d # 63とandを取る(これがxとなる) callq __Z2rand64v # rand64関数をコールする andl $63, %eax # 63とandを取る(これがyとなる) shlq $8, %rax # yを8ビットシフトする # (intが4バイト長であるため、6 + 2ビットシフトしている) addq %r12, %rax # aの先頭のアドレスを足す xorl (%rax,%r14,4), %r15d # さらに、xの4倍分を足したアドレスの値とxorを取る decl %ebx # %ebxレジスタから1を引く jne LBB1_1 # %ebxが0でなければ反復する
これに比べ、一次元配列を使ったtest_b関数のアセンブラは以下のようになります。反復の内側で実行されるニモーニックの数がとても少ないことが分かります。
xorl %r14d, %r14d # %r14dレジスタを0にする movl $10000, %ebx # %ebxレジスタを10000にする leaq __ZL1b(%rip), %r15 # 二次元配列aの先頭のアドレスを%r15レジスタに格納する LBB2_1: callq __Z2rand4096v # rand4096関数をコールする andl $4095, %eax # 4095とandを取る(これをpとする) xorl (%r15,%rax,4), %r14d # b + p × 4のアドレスの値とxorを取る decl %ebx # %ebxレジスタから1を引く jne LBB2_1 # %ebxが0でなければ反復する
たとえ二次元配列を2の累乗の数で確保していたとしても、アドレッシングの計算の際にビットシフトと加算、また多少のレジスタのコピー操作が入ってしまうのです。たったこれだけなのですが、速度的には相当分が悪いのです。
キュー
さて、さらなる速度改善を考えてみましょう。具体的にはstd::queueもより原始的なデータで置き換えられないでしょうか?
キューの先頭と末尾の位置を保持するインデックスを2つ用意すると、一次元配列でキューを作成することができます。例えばアルゴリズムCに掲載されているキューの実装は以下のようにとてもシンプルです。
#define max 100 static int queue[max+1],head,tail; put(int v) { queue[tail++] = v; if (tail > max) tail = 0; } int get() { int t = queue[head++]; if (head > max) head = 0; return t; } queueinit() { head = 0; tail = 0; } int queueempty() { return head == tail; }
このデータ構造はstd::queueと比べ、ずっとずっと原始的です。例えば、キューに溜めるべき要素の数がmax個を超えた場合、プログラムは破綻します。そのため、キューに溜まる最大の要素の数を予め考察しておくことが大切なのですが、なかなか難しく、テストを繰り返しながらパラメータを調整する手法に頼りがちです。
問題によっては発想を少し変えることによってキューの実装がより簡単で安全なものになります。今回の迷路の最短経路探索にて使用している幅優先探索はそのような問題の一つです。
幅優先探索が訪れる最大の頂点数がどのくらいか考えてみましょう。すでに訪れている頂点にはもう訪れることがないことを考えると高々N2です(Nは頂点数ではなく、迷路の幅であることを思い出してください)。そのため、要素数がN2の一次元配列を確保しておけばどのような状況であってもキューがはみ出さないことが分かります。ここではアライメントを考慮して、642分の要素を確保することにします。さらに、putやget関数を以下のように簡略化することができます。
#define max 4096 static int queue[max],head,tail; put(int v) { queue[tail++] = v; } int get() { return queue[head++]; }
if文と代入を取り除いているだけですが、速度に大きく影響することにも注意してください(自分の端末だと、10%ほど速度に影響を与えました)。
さて、では実装を見てみましょう。メンバ変数としてqueue_を用意します。
class solver_t { private: int N_; int LASTROW_; char a_[4096]; bool visit_[4096]; int queue_[4096]; };
bfsの実装は以下のようになります。
void solver_t::bfs(int sx, int sy, int ex, int ey) { memset(visit_, 0, sizeof(visit_)); int s = XY(sx, sy), e = XY(ex, ey); visit_[s] = true; int q1 = 0, q2 = 0; queue_[q1 ++] = s; while (q1 != q2) { int p = queue_[q2 ++], pp; int x = X(p); if (p >= 64 && ! visit_[pp = p - 64] && a_[pp] == '.') { if (pp == e) return; visit_[queue_[q1 ++] = pp] = true; } if (x > 0 && ! visit_[pp = p - 1] && a_[pp] == '.') { if (pp == e) return; visit_[queue_[q1 ++] = pp] = true; } if (p < LASTROW_ && ! visit_[pp = p + 64] && a_[pp] == '.') { if (pp == e) return; visit_[queue_[q1 ++] = pp] = true; } if (x < N_ - 1 && ! visit_[pp = p + 1] && a_[pp] == '.') { if (pp == e) return; visit_[queue_[q1 ++] = pp] = true; } } };
これまでと同様に3秒間の呼び出し回数を計測してみましょう。結果、373,280回となりました。初めに比べると3.5倍以上の速度改善を達成したことになります。
経路復元
さて、次は機能拡張をしてみましょう。bfsメソッドが幅優先探索によって発見した最短経路を返すように変更します。つまり、経路復元を行います(経路復元については@machyさんが去年のアドベントカレンダーで記載したスライドが大変秀逸ですのでご一読することをお勧めいたします)。
経路復元は以下のように行います。
- 「その頂点にはどの頂点から訪れたか?」という情報を別途記録しておく
- 終点を発見したら、始点までさかのぼる
「その頂点にはどの頂点から訪れたか?」という情報は相対位置情報(上から来た、下から来た等)でも構いませんが、ここでは絶対位置情報を記録することにします。
探索が終わったときの経路情報は以下のようなものとなります。見方を変えると、これは始点を根とした木であるとも言えます(薄く描画されている辺は、終点が発見されたことにより探索が打ち切られ、調べられなかった辺を表しています)。
さて、メンバ変数を追加しましょう。木を意識して、parent_という変数名を選択しました。
class solver_t { private: int N_; int LASTROW_; char a_[4096]; bool visit_[4096]; int queue_[4096]; int parent_[4096]; };
bfsの実装は以下のようになります。0はパックされた座標として有効ですので、-1で根を表すことにしました。また、探索を行う毎にparent_の内容を全て初期化する必要がないことに注意しましょう。探索と共に経路情報が木のように育っていることを意識すると分かりやすいかもしれません。
std::vector<int> solver_t::bfs(int sx, int sy, int ex, int ey) { memset(visit_, 0, sizeof(visit_)); visit_[XY(sx, sy)] = true; parent_[XY(sx, sy)] = -1; int e = XY(ex, ey); int q1 = 0, q2 = 0; queue_[q1 ++] = XY(sx, sy); while (q1 != q2) { int p = queue_[q2 ++], pp; int x = X(p); if (p >= 64 && ! visit_[pp = p - 64] && a_[pp] == '.') { parent_[pp] = p; if (pp == e) break; visit_[queue_[q1 ++] = pp] = true; } if (x > 0 && ! visit_[pp = p - 1] && a_[pp] == '.') { parent_[pp] = p; if (pp == e) break; visit_[queue_[q1 ++] = pp] = true; } if (p < LASTROW_ && ! visit_[pp = p + 64] && a_[pp] == '.') { parent_[pp] = p; if (pp == e) break; visit_[queue_[q1 ++] = pp] = true; } if (x < N_ - 1 && ! visit_[pp = p + 1] && a_[pp] == '.') { parent_[pp] = p; if (pp == e) break; visit_[queue_[q1 ++] = pp] = true; } } std::vector<int> path; for (int p = e; p != -1; p = parent_[p]) path.push_back(p); std::reverse(std::begin(path), std::end(path)); return path; };
さて、今回は絶対位置情報を使いましたが、何れのデータ構造が得であるかは問題と場合に依ります。相対位置情報で記録したほうが省メモリかもしれません。ですが、経路復元をする際に座標値の再計算が必要となるかもしれません。また、絶対位置情報では再計算を行わないかもしれませんが、相対的な方向を解として保持しなけらばならない場合、やはり再計算が必要となるかもしれません。
コスト付き最短経路探索
最後にコスト付き最短経路探索を考えてみましょう。頂点のコストを加味した最短経路探索はこれまで扱ってきたコストを考えない最短経路探索よりもずっと複雑ですが、より現実的な問題設定であると言えると思います。
ここでもこれまでに扱ってきた迷路のデータを使用しましょう。これまで壁は通れないものとして扱ってきましたが、通り抜けられるようにしてしまいましょう。そして、以下のようにコストを設定してみましょう。
- 通路を一つ進む度にコストが1かかる
- 壁を一つ乗り越える度にコストが10かかる
この条件下では壁を一つ乗り越えることで、10個以上の通路をスキップしてしまったほうがよい状況を考えなければいけなくなります。例えば、以下の例では始点のすぐ下で一回壁を乗り越えたほうが得であるようです(参考として、壁を乗り越えない場合の最短経路を薄い色の経路で表しています)。
随分と実際のマラソンマッチの様相を呈してきました。事実、この問題設定は冒頭で紹介したMM82のものにとても類似しています。
先ほども触れましたが、この問題設定はなかなか複雑で難しいものです。前半でやったのと同様に、迷路をグラフとして表現してみましょう。以下のような疎ではありますが、先ほどまでと比べるととても密度の高いグラフとなります(白い頂点は通路を表し、黒い頂点は壁を表しています)。
さて、この最短経路探索はどのようなアルゴリズムを用いるのが適切でしょうか? アルゴリズムに詳しいかたは迷いなくダイクストラ法を選択するでしょう。
二次元配列を使った実装から見て行くことにしましょう。優先度付きキューを伴ったダイクストラ法の実装は計算量を抑えることが知られています。std::priority_queueを使いましょう。コストが低い頂点から優先して処理を行うため、std::greaterを併用していることに気を付けてください。
コストの上限は64 × 64 × 10として40960ですので、16ビット用意すれば大丈夫そうです。また、コストは初期化時に12ビット左にシフトしてしまうことにしました。位置情報と共にパックする際の実装が簡潔になるからです。そのため、パック、アンパックをするためのマクロも改変します。
#define X(x) ( (x) & 63) #define Y(x) (((x) >> 6) & 63) #define Z(x) ( (x) & -4096) #define P(x) ( (x) & 4095) #define XY(x, y) (((y) << 6) + (x)) #define PZ(p, z) ((z) + (p)) #define XYZ(x, y, z) PZ(XY((x), (y)), (z))
ビットの使い方は以下のようになりました。
では実装を見て行きましょう。a_にはもはや文字を収めません。12ビット分左にシフトしたコストを記録します。そのため、char型の配列ではなく、int型の配列として確保します。
class solver_t { private: int a_[64][64]; int N_; int cost_[64][64]; };
迷路のデータを受け取り次第、12ビット分左にシフトします。
void solver_t::init(const std::vector<std::string>& a) { N_ = a.size(); for (int i = 0; i < N_; i ++) for (int j = 0; j < N_; j ++) a_[i][j] = (a[i][j] == '.' ? 1 : 10) << 12; };
実装は以下のようになります。
int solver_t::dijkstra(int sx, int sy, int ex, int ey) { for (int i = 0; i < N_; i ++) for (int j = 0; j < N_; j ++) cost_[i][j] = INF; std::priority_queue<long long, std::vector<long long>, std::greater<long long>> pq; pq.push(XYZ(sx, sy, cost_[sy][sx] = a_[sy][sx])); while (! pq.empty()) { int x = X(pq.top()); int y = Y(pq.top()); int z = Z(pq.top()); pq.pop(); if (cost_[y][x] < z) continue; for (int k = 0; k < 4; k ++) { int xx = x + dx[k], yy = y + dy[k]; if (0 <= xx && xx < N_ && 0 <= yy && yy < N_) { int zz = z + a_[yy][xx]; if (xx == ex && yy == ey) return zz; if (zz < cost_[yy][xx]) { cost_[yy][xx] = zz; pq.push(XYZ(xx, yy, zz)); } } } } return -1; };
3秒間の計測を行います。39,490回でした(参考までに、STLを使った場合の実装は30,900回でした)。
さて、こちらを一次元配列を使って書き直してみましょう。std::priority_queueはそのまま使用します。
class solver_t { private: int a_[4096]; int N_; int LASTROW_; int cost_[4096]; };
int solver_t::dijkstra(int sx, int sy, int ex, int ey) { for (int i = 0, p = 0; i < N_; i ++, p += 64) for (int j = 0, pp = p; j < N_; j ++, pp ++) cost_[pp] = INF; int s = XY(sx, sy), e = XY(ex, ey); std::priority_queue<int, std::vector<int>, std::greater<int>> pq; pq.push(PZ(s, cost_[s] = a_[s])); while (! pq.empty()) { int x = X(pq.top()); int p = P(pq.top()), pp; int z = Z(pq.top()), zz; pq.pop(); if (cost_[p] < z) continue; if (p >= 64) { pp = p - 64; if ((zz = z + a_[pp]) < cost_[pp]) { if (pp == e) { return zz; } else { pq.push(PZ(pp, cost_[pp] = zz)); } } } if (x > 0) { pp = p - 1; if ((zz = z + a_[pp]) < cost_[pp]) { if (pp == e) { return zz; } else { pq.push(PZ(pp, cost_[pp] = zz)); } } } if (p < LASTROW_) { pp = p + 64; if ((zz = z + a_[pp]) < cost_[pp]) { if (pp == e) { return zz; } else { pq.push(PZ(pp, cost_[pp] = zz)); } } } if (x < N_ - 1) { pp = p + 1; if ((zz = z + a_[pp]) < cost_[pp]) { if (pp == e) { return zz; } else { pq.push(PZ(pp, cost_[pp] = zz)); } } } } return -1; };
一次元配列を用いた場合、dijkstraメソッドは45,840回呼び出すことができました。STLを用いたものに比べ、1.5倍高速です。
さて、ここまで来ると優先順位付きキューを一次元配列を利用できないのかな? と考えてしまいますね。しかし、優先順位付きキューは単純なキューのように機械的に一次元配列等に変換することはできません。残念。
ここで、優先順位付きキューの代わりに単純なキューを使うと何が起こるか考えてみましょう。この問題では負のコストがないため、優先順位付きキューを使えば最もコストが少ない頂点から訪れることが保証されていました。しかし、単純なキューを使った場合、例えば終点に初めて訪れた際のコストが最小である保証はなされません。
ならば、キューが空になるまで繰り返し、それぞれの頂点に訪れたコストの中で最小のものを記録するようにしてみましょう。言い換えれば、一つ前の頂点からの緩和を繰り返す、ということになります。何かどこかで見たことがありますね...? そうです、これはキューを伴ったベルマンフォードアルゴリズムです。
int bellman_ford(int sx, int sy, int ex, int ey) { for (int i = 0, p = 0; i < N_; i ++, p += 64) for (int j = 0, pp = p; j < N_; j ++, pp ++) cost_[pp] = INF; int s = XY(sx, sy), e = XY(ex, ey); int q1 = 0, q2 = 0; queue_[q1 ++] = PZ(s, cost_[s] = a_[s]); int zm = INF; while (q1 != q2) { int q = queue_[q2 ++]; int x = X(q); int p = P(q), pp; int z = Z(q), zz; if (cost_[p] < z || z > zm) continue; if (p >= 64) { pp = p - 64; if ((zz = z + a_[pp]) < cost_[pp]) { if (pp == e) { zm = std::min(cost_[pp] = zz, zm); } else { queue_[q1 ++] = PZ(pp, cost_[pp] = zz); } } } if (x > 0) { pp = p - 1; if ((zz = z + a_[pp]) < cost_[pp]) { if (pp == e) { zm = std::min(cost_[pp] = zz, zm); } else { queue_[q1 ++] = PZ(pp, cost_[pp] = zz); } } } if (p < LASTROW_) { pp = p + 64; if ((zz = z + a_[pp]) < cost_[pp]) { if (pp == e) { zm = std::min(cost_[pp] = zz, zm); } else { queue_[q1 ++] = PZ(pp, cost_[pp] = zz); } } } if (x < N_ - 1) { pp = p + 1; if ((zz = z + a_[pp]) < cost_[pp]) { if (pp == e) { zm = std::min(cost_[pp] = zz, zm); } else { queue_[q1 ++] = PZ(pp, cost_[pp] = zz); } } } } return zm; };
実装部分はまあこれでよいとして、とても注意すべきことがあります。もはやキューの最大要素数はN2ではありません。緩和が終わるまでそれぞれの頂点に複数回訪れることを考えると不思議ではありません。
この問題でデータのアライメントを64と設定した場合、4,096の16倍、65,536個分を用意する必要がありました。理論的に計算できればカッコいいのですが、私はこれを実験的に決定しました。実際には50,000個ほどでも十分だと思います。
int queue_[65536];
優先順位付きキューを伴ったダイクストラ法では45,840回の探索でしたが、一次元配列によるキューを伴ったベルマンフォード法では59,600回の探索を行うことが出来ました。ほぼ1.3倍の速度改善です。STLを使ったものに比べると約2倍のパフォーマンス改善を行ったことになります。
この実装のポイントとして、以下の3種類の枝刈りを行っていることに触れておきましょう。
- 注目する頂点により低いコストで到達する経路がある場合
- 注目する頂点にたどり着いたコストが、現在までに見つかっている終点までのコストよりも大きい場合
- 隣接する次の候補の頂点により低いコストで到達する経路がある場合
2番目の枝刈りの条件は忘れやすいと思いますが(「z > zm」の条件)、特に重要です。これを忘れると、ダイクストラ法のパフォーマンスに遥かに及びません。
「ベルマンフォード法の計算量はダイクストラ法と比べ随分大きいんじゃないのか? 遅くなるんじゃないのか?」と考えたかたも多いのではないのでしょうか? 確かにその通りで、ベルマンフォード法を適用できるか否かは問題の性質に大きく依存します。例えば大部分が同じ設定であったとしても、盤面の大きさが999 × 999である場合は以下のことが起こります。
- キューの要素数を相当大きく確保しておく必要がある(例えば、67,108,864要素(= 1024 × 1024 × 64))
- 優先順位付きキューを伴ったダイクストラ法と比べ、半分以下の速度となってしまう。
最後に、この条件下での経路復元について簡単に触れておきます。キューを伴ったベルマンフォード法の場合も幅優先探索と変わらず、「その頂点にはどの頂点から訪れたか?」という情報を記録すれば十分に経路復元をすることが可能です。
こちらの問題設定では先ほどと比べ、多少複雑なグラフを扱うことになります(以下の図は再掲)。
しかし経路の情報は、やはり始点を根とした木の様相を成します。以下の図をよく見ると、薄い辺があります。これらの辺は緩和の課程で一度は接続されたものの他の最短な経路によって切断された辺です。ベルマンフォード法では同じ頂点に何度も訪れるため、親の頂点情報も何度も入れ替わるのです。
まとめ(または、回顧録)
今回はマラソンマッチにおける典型データ構造とアプローチという題材で多少のトピックを解説しました。それぞれのトピックを見直すと、以下の共通点があることが挙げられます。
- 原始的なデータ構造およびロジックは速い
- 二次元配列は意外に遅い
- 一次元配列は十分に単純であるため、速い
- メモリを動的に確保するのは遅い
- 問題に適したミニマムかつ効率的な設計を行う
「ひゃー、マラソンマッチではこんなに雑多な実装を行わなきゃいけないのか。つまらないじゃないか」と思われた方もいらっしゃるかもしれません。再掲となりますが、マラソンマッチで真に大切なのは問題にフィットしたロジックであり、アルゴリズムです。このエントリに記載されていることはある意味「おまけ」なのですが、トップ陣はそのおまけの部分も異様に強い。まるで魔術師を見ているようです。しかし、一つ一つのテクニックを丁寧に観察しひも解くと、実装の大部分に明確な理屈があり、ここで記載したような小さな工夫やチューニングが巧みに全体を構成しているのが見えてきます。このエントリがそんな「面倒臭い部分」に向き合うときの少しばかりの助けになれば大変嬉しく思います。
さて、この記事はCompetitive Programming Advent Calendar 2014 - PARTAKEのマラソン部門、12月6日分のエントリとして記載しました。これまでも、またこれからも続々とマラソンに関する記事が執筆されます。大変貴重なイベントだと思います。
アドベントカレンダーが始まってから実力者のエントリが続き、「問題をどう解くか」という核心の部分を熱く熱く執筆くださっています。素晴らしい記事を読むにつれ、自分の選んだ題材はあまりにも稚拙だよなあ? と迷いながらも、特に日本のマラソンマッチコンペティターの底上げに少しでも貢献できればと思い、なんとか書き上げました。
さて、そもそも今年の競技プログラミングアドベントカレンダーにマラソン部門が出来たことに関して少し触れておきたいと思います。私が日本に一時帰国した際にhotpepsiさんが開催してくださった飲み会の席でのことでした。突然、simezi_tanが熱く、そして激しく、「マラソンは何をしていいか分からない。競技プログラミングのアドベントカレンダーにマラソン部門が欲しい!」と声を上げたのです。tomerunさんやmachyさん達の援護が加わったこと、また、_tanzaku_さんが柔軟な企画、そして対応をしてくださったこともあり、気がついたときには見たこともない素晴らしいイベントとして運用されることになっていました。日本の競技プログラミングコミュニティの柔軟さには舌を巻くばかりです。
そんなsimezi_tanに伝えたいことがあります。この場を借りて記したいと思います。
「おいしめたん。そもそもの言い出しっぺなんだから参加しなはれ(何が分からないか書いてみるとか)」
明日、12月7日はyowaさんが「マラソンマッチで1位になる方法(最終順位とは言ってない)」という題材で執筆なさる予定です。お楽しみに。