AHC045 参加記
THIRD プログラミングコンテスト2025（AtCoder Heuristic Contest 045） お疲れ様でした。
システムテストの得点が 2,972,941,373,137 点で 2 位でした。
問題概要
個の点が二次元平面上にあり、各頂点のおおよその位置が長方形で与えられます。二点間の距離はユークリッド距離で定義されます。
最大 回のクエリを投げることができます。各クエリでは最大
個の頂点を指定することができ、ジャッジから指定された頂点の最小全域木で使われる辺の集合が返されます。
グループ数と各グループの頂点数が与えられます。頂点をグループに分け、各グループでできるだけ距離の総和が小さい木を構成してください。
詳しくは公式の問題文を参考にしてください。
解法
クエリを使って座標を推定する「推定パート」と、グループ分けのコストを最小化する「構築パート」で分けて説明します。
表記
- パラメータについて
のように書いたら、
は
の、
は
のパラメータを表すことにします。
- ablation study のように、アイディアごとの寄与を求めました。その要素を抜いた 300 ケースの相対スコアの悪化度です。目安として 0.1 % ぐらいずれることがあります。
推定パート
サンプリングを行いました。サンプル数は 64 ~ 8 です。
サンプリングを行うことで、二点間の距離の期待値などを簡単に推定することができます。
座標を推定する部分と、クエリの頂点集合を最適化する部分で分けて説明します。
座標の推定
今までのクエリ結果を利用して、サンプルごとに座標の更新を行います。
更新の流れを以下に記します。
詳細について解説します。
ランダムに動かす (寄与 3 %)
更新前の 座標を
、更新の度合いを
として
から一様ランダムに選んだ値で更新します。
座標も同様です。現在の座標
と長方形の共通部分から一様ランダムに選ぶことになります。
全ての制約を満たしたときに、勾配がなくなって更新が停止するよりも、ランダムウォークしたほうが良質なサンプリングになると考えてこのような処理を入れました。
今回の勾配法には局所解が存在するため、局所解からの kick という意味もありそうです。
不等式制約の生成
全域木 が「最小全域木である」ことは、「
で使われていない任意の辺
は、
を追加して生成される閉路の任意の辺より長い（短くない）」や「
で使われる任意の辺
は、
を削除してできる 2 つの連結成分を結ぶ任意の辺より短い（長くない）」ことと同値です。自明ではないと思いますが、そこそこ有名な定理です。
この性質を利用すると、木の高さが定数で抑えられると仮定して 個の 2 辺の長さに関する不等式が得られます。
余計な不等式はあらかじめ削除します。具体的には、すでに得られた制約と等しい制約は追加しないようにします。さらに、全てのサンプルが満たしている制約も追加しないようにします。（制約削除によるスコアへの寄与は確認できませんでした。）
に関する不等式制約をシャッフルする (寄与 0.05 %)
制約1つずつに対して勾配法を適用するため、制約を見る順序によって結果が変わります。順序が固定なのはよくないかもしれないと考え、制約をシャッフルしました。
勾配法による更新
最急降下法を使いました。目的関数の勾配を求め、その方向に だけ進みます。制約を満たしているときは、勾配がないので座標を更新しません。
座標を更新した後に が入力の長方形から出た場合、 clamp して座標を長方形の辺上に移動させました。
を基準とした不等式制約として、3 種類のケースがあります。
の場合
が
よりも
に近いという制約です。
制約を満たしていないとき、線分 の垂直二等分線からの距離をコストとしました。勾配の方向は
の微分と等しいです。3
数式で表すと以下のようになります。
の場合
右辺は定数なので、 とします。
を中心とする半径
の円の内側に
が存在するという制約です。
制約を満たしていないとき、 から円までの距離 (=
から
までの距離 -
) をコストとしました。
数式で表すと以下のようになります。4
の場合
を中心とする半径
の円の外側に
が存在するという制約です。
制約を満たしていないとき、 から円までの距離 (=
-
から
までの距離) をコストとしました。
非凸であることが気になりますが、ランダムに動かすところである程度の対処がされたと考えています。5
数式で表すと以下のようになります。
SIMD による高速化
AVX を使って勾配法を高速化しました。並列化の軸はサンプルです。同じ不等式制約に対して 8 つのサンプルを同時に更新します。
データの持ち方について触れておきます。8 つのサンプルごとに points[vertex ID][x or y][sample ID] のような固定長配列を保持しました。sample ID の軸を最後にすることで、データの読み込みや書き込みも並列化されます。SIMD を使わない場合でも、キャッシュ効率が改善されることが期待されます。
seed = 0 で実験したところ、SIMD を使わない場合の勾配法の回数は 回、使う場合は
回でした。6
SIMD を使わない場合、不等式制約を満たしていれば後の更新処理を行わないため、並列数倍ほどの高速化は望めないと思います。
クエリの最適化
以下の流れで、クエリで聞く頂点を選びます。
- サンプルにおける分散が大きい頂点を選択し
とする。
- 各頂点に対して
- 距離の期待値が最も小さい頂点を
とし、最初の頂点集合を
とする。
- 距離の期待値が小さい 30 個の頂点を候補とする。
- 後述する評価値を最大化するように、頂点を貪欲に 1 つずつ追加する。
詳細について解説します。
の選び方
サンプルを用いて各頂点の位置の分散を求めることができます。
まだ として 1 回も選ばれていない頂点のうち、分散が最も大きい頂点を
としました。
として選ばれる回数を 1 回に限定することで、選ぶ頂点が多様化されると考えています。（1 回に限定する寄与: 0.9 %）
ただし、 の場合だけ、3 回まで同じ頂点を選択できるようにしました。
頂点を追加する貪欲法
評価値を最大化する頂点集合を求めることを目標としました。評価値の計算コストが無視できないため、 に近い点 30 個に限定して、貪欲法で 1 つずつ頂点を追加しました。近い点に限定したのは、構築パートで欲しいのが近い点に関する情報だと考えたからです。
評価関数は、各サンプルの最小全域木と実際の最小全域木のコストの差の期待値を最大化するイメージで設計しました。
番目のサンプルにおける最小全域木を
、
番目のサンプルにおける木
のコストを
とします。
評価関数は以下のようにしました。
は定数とします。
だと差の期待値になり、
だとコストの差の期待値の割合のようなものになります。
貪欲法の高速化
の部分については、各辺について最小全域木に使われた回数を数えることで高速化できます。
貪欲法のボトルネックは、各サンプルの最小全域木 を求める部分です。頂点数を
として、愚直に Kruskal 法を適用すると
、Prim 法を適用すると
の時間計算量になります。
ここで、選ばれることが確定した 頂点の最小全域木がすでに求まっていることを利用します。
頂点の最小全域木として使われなかった辺は、頂点が追加されても使われることはありません。そのため、辺の候補は
頂点の最小全域木として使われた
辺と、追加した頂点と他の
点を結ぶ
辺のみになります。辺の数が
個に絞られたので、Kruskal 法で
で最小全域木を構成できます。
構築パート
クエリと勾配法で 1.5 秒ほど使った後、以下のような処理を行います。
- 全ての2点間の距離の期待値を求める。
- 使用可能な辺を限定する。
- 初期解を生成する。
- グループ数とグループサイズを合わせる焼きなまし法を行う。
- コストを減少させる焼きなまし法を行う。
- クエリ結果を利用して辺を置き換える。
詳細について解説します。
2点間の距離の期待値
サンプリングの直接的な目的は、点の位置を推定することではなく、2点間の距離の期待値を求めることです。不確かな点からは距離の期待値が小さい頂点が存在しないことを想像すると違いがわかると思います。
点の位置の期待値を使う場合と比較した寄与は 1.5 % でした。
使用可能な辺の限定
この後に焼きなましをする都合で、 辺全てを扱うのは計算コストが重いため、使える辺を限定しました。
端的に表現すると、「2 つの最小全域木」で使われる 辺に限定しました。1 つ目の最小全域木は普通の最小全域木で、2 つ目の最小全域木は、1 つ目の最小全域木で使われる辺を削除した後の最小全域木です。
初期解の生成
Kruskal 法を、連結成分の個数が 個になるまで行って初期解を生成しました。一般にはグループサイズの制約を満たしていないため、無効な解となります。
運良くグループサイズの制約が満たされた場合には、構築パートの最適解が得られることになるため、このような初期解生成を行いました。
グループ数とグループサイズを合わせる焼きなまし
有効な解が見つかるまで焼きなまし法を行います。大体 0.02 秒もかからずに有効な解が見つかるため、実行時間は短いです。
遷移は次の 2 つです。
- merge: 異なる 2 つのグループを結ぶ辺を追加する。
- split: 使われている 1 つの辺を削除する。
評価値は、
- 辺のコストの総和
- グループサイズのペナルティ
の和としました。グループサイズのペナルティというのは、サイズ のグループが
個過剰にある（
なら不足）とき
になります。
はペナルティの係数で 1000 としました。
評価値の差分計算は で、更新は愚直に行いました。各木について適当に根を指定して有向木とし、各頂点について、
- 親
- 子供
- 部分木の大きさ
- 根
の情報を持たせます。根の部分木の大きさは木の大きさとなります。木の大きさと部分木の大きさがわかるので merge も split も で評価できます。
コストを減少させる焼きなまし (寄与: 10.3 %)
グループサイズの制約を満たしたままコストの減少を試みる焼きなまし法を行います。
評価は辺のコストの総和としました。
遷移は 2 種類です。
遷移1: merge して split
merge した後に split します。merge したグループを split することもあれば、異なるグループを split することもあります。
merge では使われていない辺を適当に選びます。split では、使われている辺を 1 つずつ見ていき、削除することでサイズ制約が満たされ、評価値の差分が許容範囲内であれば更新を行いました。
遷移2: グループ内の辺の置換 (寄与: 0.06 %)
同一グループ成分内の 2 頂点を結ぶ辺を追加します。閉路が 1 つでき、閉路を構成する辺（追加した辺を除く）から最も長い辺を削除します。
部分木のサイズを見ながら親を辿っていくことで、余計な辺を見ることなく閉路を求めることができます。
各グループの最小全域木を求める (寄与: 0.15 %)
各グループで使われる辺が最小全域木になっているとは限らないため、最小全域木で辺を置換しました。
コストを減少させる焼きなましで best な解が更新される度にこの処理を行いました。
クエリ結果を利用して辺を置き換える (寄与: 0.06 %)
焼きなまし法が終わってから、クエリ結果を利用した解の改善を行います。
不等式制約における長い辺が使われていて、短い辺で置き換えられるとき、辺を置換しました。
サンプルが全ての不等式制約を満たしていればこの処理は不要なので、制約を満たさないサンプルがいくつかあったと考えています。
その他の工夫
が小さく 
が大きいときに一部の点の推定を諦める (寄与: 0.02 %)
が大きいとき、孤立点となる頂点が存在します。
が小さいとき、入力の長方形の面積が大きい頂点は推定が難しく、推定しても距離の期待値が近い頂点は存在しないと予想されます。そこで、そのような頂点は孤立点になると決め打って最初から推定を諦めました。
提出コード
コンテスト中の提出を微修正しました。
最後に
できれば優勝したかったという気持ちはあるものの、久しぶりに長期 AHC で上位に入れてよかったです。
過去の AHC から得た知見を活かせたのもよかったです。サンプリングの強さは AHC025 で実感していて、SIMD を使ったサンプリング並列化は AHC040 の terry_u16 さんの記事 で学びました。また、最小全域木の性質は去年の大学の授業に出てきたので憶えていました。
コンテストを主催してくださった THIRD 様、Writer の Shun_PI さん、コンテスト参加者の皆様、最後まで読んでくださった読者の皆様、ありがとうございました。
AHC033 参加記
トヨタ自動車プログラミングコンテスト2024#5（AtCoder Heuristic Contest 033 お疲れ様でした。
システムテストの得点が 1,971,936,410,473 点で 3 位でした。
問題概要
マスの空間で、5 台のクレーンを操作し、左の搬入口から右の搬出口にコンテナを運ぶ問題です。
詳しくは公式の問題文を参考にしてください。
最終提出の方針
ビームサーチをしました。
ビームサーチの 1 ステップは以下の操作です。
- 1 台のクレーンが、コンテナの場所まで移動してコンテナを掴み、コンテナを運んで離す。
40 ステップほどで全てのコンテナを搬出できます。
ビーム幅は 25,000 ぐらいです1。
遷移
盤面上にある各コンテナについて、その場所に最も早く到着できるクレーンがそのコンテナを運び出します2。
搬出口以外の場所にコンテナを置くことを「スタックする」と定義することにして、1 つのコンテナがスタックされるのは 1 回までとしました。すなわち、スタックされているコンテナを動かすときは必ず搬出します。
スタック位置は 2 列目と 3 列目のマスのみに限定しました。1 列目に置くことを許すと、その左の搬入口を小クレーンが使えなくなるからです3。
経路は BFS の結果から 1 通りだけ復元しました。ちなみに 2 通り以上の復元も試したのですが、多様性が失われて悪化したのでやめました。経路の過去改変については、計算コストや実装コストが重そうだったので試していません。
遷移の実装
x 座標、 y 座標、ターン数の 3 次元空間を幅優先探索 (BFS) します。
BFS はビットボードを使って実装しました。上下左右の遷移はビットシフトとビットマスクで計算できます。
コンテナを拾いに行くときは各クレーンから BFS する必要があり、これも同時に行いました。つまり、拾いに行く BFS は状態ごとに 1 回だけ行えば十分です。
コンテナを置きに行くときはコンテナごとに 1 回 BFS すれば十分で4、経路の復元は採用されることが分かってから行うようにして少し高速化できました。
到着条件
経路を決定する順序とマスを使用する順序が逆転しうるため、クレーンの到着条件に注意して実装しました。
- コンテナを搬入するとき
- その搬入口から最後にコンテナを搬入した時刻より後
- コンテナを搬出するとき
- その搬出口から最後にコンテナを搬出した時刻より後
- スタックされたコンテナを拾いに行くとき
- そのコンテナがスタックされた時刻より後
- コンテナをスタックするとき
- 小クレーンがコンテナを掴んで通った最後の時刻より後5
いずれの場合についても、クレーンが到着した直後に操作を行うため、到着した直後のターンに他のクレーンが使用していないことも必要条件となります。
QP の削除
コンテナを離して置くという行動を連続して行うのは無駄なので省略できるようにしました。ある種の過去改変です。
大クレーンが小クレーンを抜かすときなどに効果があると思います。
評価
評価関数はスコアの下界をベースにして実装しました。
スコアの下界
ある種の緩和問題を解けばよいです。
1 つのクレーンの 1 つの行動を労働力 1 と定義します。
全てのコンテナを搬出するには、最低でもどれだけの労働力が必要かを考えます。
まず、y 軸方向の移動について考えます。コンテナを運ぶときは、y = 0 から y = 4 まで運ぶ必要があり、コンテナを掴んで離す行動も含めると、PRRRRQ で 6 の労働力が必要です。コンテナを拾いに行くときに右端から左端に移動するため、最初や最後を除いて LLLL で 4 の労働力が必要です。
次に、x 軸方向の移動について考えます。コンテナを運ぶときは、搬入口と搬出口の x 座標の差の絶対値だけ移動する必要があります。
最後に、スタックのコストについて考えます。スタックする場合、PQ で 2 の労働力が余計にかかります。総スタック数の最小値は動的計画法で求まります。搬入口 から搬入したコンテナの数を
として
を求めればよいです。
今までの評価についてまとめると次のようになります。
- 評価するもの
- コンテナを掴んで離す行動
- y 軸方向の移動
- コンテナを運ぶときの x 軸方向の移動
- 評価しないもの
- コンテナを拾いに行くときの x 軸方向の移動
- スタックされたコンテナによって経路が伸びること
- クレーン同士の干渉
- etc
クレーンが 5 台あるため、必要な労働力を 5 で割ることでスコアの下界が求まります。
労働力の評価
今までに使った労働力と、コンテナの搬出が完了するまでに必要な労働力の和を評価します。
具体的には以下の項目を評価すればよいです。
- 搬入したコンテナの集合
- スタックしたコンテナを搬出するコスト
- クレーンの x 座標
労働力の和を 5 で割った値と、クレーンの最大行動長の max をとってスコアの下界を求めました。
スコアの下界が評価値のベースで、この時点で評価できていないもので評価値を修正します。
スタック位置の評価
スタックされたコンテナによって小クレーンの搬出経路が長くなることがあるため、スタック位置によってペナルティを課しました。
以下の項目を評価しました。
- スタック数
- 通れなくなった行数
- 隣り合うコンテナの相性
スタック数はその時点でスタックしているコンテナの数のことです。
通れなくなった行数は、スタックされているコンテナを含む行数のことです。1 列目から 4 列目まで直線的に移動できる行が少ないと搬出経路が長くなりやすく、さらにクレーン同士の干渉が増えると考えました。通れなくなった行数が 0 行から 1 行に増えても少ししか変わらない一方で、4 行から 5 行になると小クレーンにとって非連結になって影響が大きいことなどを考慮し、2 乗のオーダでペナルティを与えました。
隣り合うコンテナの相性は、同じ行に 2 つのコンテナがスタックされているとき、その 2 つのコンテナの搬出先が同じで左に番号が小さいコンテナが置いてあるときにペナルティを与えました。特にスタックされている行と搬出先の行が一致するときは、小クレーンで搬出するときに x 軸方向に動かす必要があります。この評価はスライドパズルの Linear Conflict から着想を得ました。
評価の比重としては、通れなくなった行数が最も大きいです。
大クレーンの評価
大クレーンの労働力の大きさを小クレーンよりも少し大きくしました。
多様性の確保
Zobrist Hash が一致した状態の中で最良のものだけを残すことで多様性を確保しました。
次の項目をハッシュ値の生成に使用しました。
- 搬入したコンテナの集合
- スタックしたコンテナ番号と位置
- 大クレーンの位置
- 小クレーンの位置の集合6
どれもパターン数が少なかったので、事前に生成した 64 bit の乱数を格納した配列にアクセスしてハッシュ値を得ました。
直近の経路までハッシュ化すれば状態の完全一致を確認できると思いますが、経路を考慮しないほうが多様性を確保できてスコアがよかったです。
提出コード
立ち回り
自分のコンテスト期間内の立ち回りについて書きます。暇な人向けです。
最初の 7 日間はタスク列の焼きなましを試していました。各タスクがどの場所にあるコンテナをどこに運ぶかを表し、ルールベースでタスクを処理して評価します。タスクのクレーンへの割り当てはルールベース側で行いました。ルールベースの部分が特に複雑で、実装には多くの時間がかかりました。
スタック位置の影響が後で出てくるという問題の性質を踏まえ、時系列順に処理する手法よりも焼きなまし法のほうが強いと考えていました。
7 日目の金曜日あたりに平均 69 点ぐらい出せるようになったので提出しようと思い、上位のスコアを概算したところ、上位は 65 点ぐらい出していることに気づきました。提出するからには 1 位をとりたいという思いがあり、提出を控えました。
タスク列の焼きなましでは 65 点あたりを出せるとは思えず、別の方針を考えました。スコアの下界の見積もりを雑に行うと 60 点ぐらいになることには序盤から気づいており、65 点との乖離はあまり大きくないことから、下界を評価関数にしてビームサーチできそうだと思いました。
残り期間が 3 日しかなく、方針を変えるか迷いました。頑張って実装した解法を捨てたくないという思いもありましたが、結局、直感的によさそうだったビームサーチを試すことにしました。
3 日で上位に入れる自信はあまりなかったですが、実装 1 日目でビームサーチを動かすことができ、平均 68 点台になって元の解法を上回り、ビームサーチ方針が強いことを確信しました。
最後の 2 日はビームサーチ方針を改善して平均 64.8 ぐらいまで改善できました。
コンテスト後に焼きなまし法を使ったコードを提出してみたところ、本番 17 位相当のスコアが出ていました。方針を変えたおかげで 3 位まで上げられたと思うと、方針を変えた決断は正しかったんだろうと思いました。
なぜ最初は強そうに思えた焼きなまし法よりもビームサーチが強かったのか自分なりに考えてみました。スタック位置の影響が後から出てくる点については、ビーム幅の太さと多様性の確保でカバーできていると思いました。盤面サイズが小さいので評価しづらくても誤魔化せるという感じでしょうか。焼きなまし法だとスコアの差分が離散的すぎてランダムウォークっぽくなるからうまくいかないというのもありそうです。
最後に
AHC031 はうまくいかず、AHC032 は Writer で出られず、そろそろ成功したいと思っていたところで 3 位を取れてほっとしています。非インタラクティブ問題でも優勝できるように頑張りたいです。
コンテストを主催してくださったトヨタ自動車様、Writer の chokudai さん、コンテスト参加者の皆様、最後まで読んでくださった読者の皆様、ありがとうございました。
- 実行時間によってビーム幅を調整しているため、ビーム幅は一定ではありません。↩
- 行動数が最小のクレーンを動かすように固定するよりも、コンテナに最初に到着できるクレーンを動かすほうがよいスコアになりました。↩
- 他の上位の方は 1 列目も使っていたので、自分の実装が下手なようでした。反省点です。↩
- 「十分」と書いていますが、コンテスト中にこのことに気づかず、コンテナの行き先ごとに BFS しました。↩
- 実はこの条件だけだと不十分で、コンテナを置いた後、大クレーンが同じ場所にコンテナをスタックし、そのコンテナを再び大クレーンが掴んで運ぶ場合にコンテナが衝突します。コンテスト後に気づきました。↩
- クレーンのターン数まで含めてハッシュ化しましたが、ターン数はハッシュ化しないほうがよかったかもしれません。↩
第一回マスターズ選手権 -決勝- 延長戦
https://atcoder.jp/contests/masters2024-final を復習したところ、A問題、B問題、C問題でそれぞれ486504点、561814点、722818点とれました。
はじめに
- この記事で紹介するのは延長戦の解法であり、コンテスト中の解法ではありません。
- 解法もビジュアライザも Rust で実装しましたが、私は Rust を使った経験がほとんどないため、読みにくかったり非効率な処理があったりするかもしれません。
問題概要
ドローンを操作して目的地を通過する問題です。
風が吹いているため1、ドローンの位置や速度が正確にはわかりませんが、「加速」の代わりに「計測」を行うことができます。
詳しくは公式の問題文を参考にしてください。
方針
粒子フィルタを用いてドローンの位置と速度を推定しながら、あらかじめ決めた順路を辿るようにします。
基本的には加速と計測を交互に行いますが、風がないケースについては加速のみを行います。
粒子フィルタ
目的地を目指したいわけですが、そのためにはドローンの位置と速度を知りたくなります。ドローンの位置と速度の期待値がわかるとよく、欲を言えば位置と速度（のペア）の分布がわかるとよいです。
過去の加速度、計測値。移動結果から推定することになります。解析的に分布を求めることは難しそうなので、ある種のサンプリングを行います。
位置と速度と重みをもった 個の粒子で分布を近似します。私の実装だと
になっています。
粒子の更新方法について説明します。
加速
全ての粒子の速度に加速度を加算します。
計測
粒子ごとに、計測結果が得られる尤度を求め、重みに乗算します。
移動結果およびサンプリング
個の次の粒子が生成されるまで以下の処理を繰り返します。
- 重みに比例した確率で粒子を復元抽出し、次の状態をシミュレーションする。すなわち、風の加速度をランダムに生成して速度を更新し、位置を更新する。
- 移動結果と矛盾がなければ次の粒子として採用する。このとき、新しい粒子の重みを
とする。
ただし、残り時間が短くなった場合には 個の粒子が生成される前に処理を打ち切ります。
稀に粒子が1つも生成されず、推定に失敗することがあります。目的地を通過したときは位置を目的地で更新するなど、適当にごまかしました。
粒子が更新される様子を可視化してみました。
青緑色の点は計測前の粒子の位置を表しています。
上方向に向かって計測すると、かなり小さい距離が測定され、粒子の位置が真のドローンの位置に寄っていることがわかると思います。
ちなみにコンテストで優勝したチーム Imusu さんが粒子フィルタを使ったらしく、私は今回のコンテストで初めて粒子フィルタを知りました。
経路の決定
最初にドローンの大まかな経路を決定します。
問題を簡単にするため、ドローンは直線的な移動しかできないものとします。
スタート地点と目的地のみの場合、壁によって直接移動できないケースが多く、しばしば非連結になります。そこで、壁の端点と別の壁の中で最も近い点との中点などを追加しました。
二点間の移動コストをユークリッド距離の平方根で定義します。このコストは、二点間の移動に必要なターン数と大体比例します。
さらに、狭いところはなるべく通りたくないため、二点を結ぶ線分と最も近い壁との距離を計算し、距離が短すぎる場合にはペナルティを課しました。
Dijkstra 法でスタート地点と各目的地から他の目的地までの最小移動コストとそのときの経路を求め、そのコストを用いて TSP の厳密解を求めました。経由点は Dijkstra 法では扱っていますが、TSP では扱わなくてよいため、TSP の頂点数は であることに注意してください。
加速度の決定
次に目指す場所を目標とし、目標の座標を とします。
ドローンの現在位置と速度の推定値を とします。
速度について、直線 と垂直な成分と平行な成分に分けて考えます。
垂直な成分
新しい速度の垂直な成分はできるだけ0になるように加速度を設定します。
平行な成分
理想的な速度に近づくように設定します。
理想的な速度というのは、目標に近づいたときに最大限減速してちょうど 目標で止まれる速度のことです。
理想的な速度は1ターンの平均加速度の大きさを 、目標までの距離を
としたとき、およそ
となります。実際には風などの影響を考慮して0.9倍した値を使用しています。
垂直な成分を決めた後に、平行な成分をいくらまで設定できるかは三平方の定理を使用して簡単に計算できます。
計測方向の決定
ドローンの現在位置の期待値を とします。
を壁に正射影した点および壁の端点のうち、
から直線的に移動できる点
について考えます。
から
への方向に向かって計測したときの期待値の標準偏差を
と
の距離で割った値が大きくなるように計測方向を決定しました。距離の分散は粒子フィルタから計算可能です。
ただし、このままだと同じ方向に連続して計測することがあったため、同じ点や同じ壁に正射影した点への計測を2回連続で行わないようにしました。
目標の変更
目的地間を移動するときに経由点を通る場合、現在位置から次の目標に直線的に移動できるときに経由点に到着したとみなし、次の目標に変更しました。
また、目標に直線的に移動できないときは目的地までの経路を生成し直しました。Dijkstra 法の結果を再利用すると簡単に計算できます。
実装
オープンコンテストでチームを組んでくださった cuthbert さんに共有していただいた幾何ライブラリの一部を移植して使用しました。
座標を複素数で管理していて、少し実装が楽になった気がします。
とはいえ、norm 関数の仕様が C++ と Rust で異なっていたり、数値誤差が問題になったり、実装にはそれなりに苦労しました。
提出コード
ビジュアライザ
yunix さんが作成したツールを使わせていただきました。
壁、目的地、ドローンの描画に加え、以下のコメントを元に描画する機能もつけました。
#e x y: 位置の推定#d distance: 計測結果#g q: 目的地の訪問イベント
作成したビジュアライザのコードを載せておきます。util.rs のビジュアライザ部分が本質的です。
lib.rs
use wasm_bindgen::prelude::*;
mod util;
const PROBLEM_ID: char = 'C';
#[wasm_bindgen]
pub fn gen(seed: i32) -> String {
util::gen(seed as u64, PROBLEM_ID).to_string()
}
#[wasm_bindgen(getter_with_clone)]
pub struct Ret {
pub score: i64,
pub err: String,
pub svg: String,
}
#[wasm_bindgen]
pub fn vis(_input: String, _output: String, turn: usize) -> Ret {
let input = util::parse_input(&_input);
let output = util::parse_output(&input, &_output).unwrap();
let (score, err, svg) = util::vis(&input, &output, turn);
Ret {
score,
err: err.to_string(),
svg: svg.to_string(),
}
}
#[wasm_bindgen]
pub fn get_max_turn(_input: String, _output: String) -> usize {
let input = util::parse_input(&_input);
let output = util::parse_output(&input, &_output).unwrap();
output.out.len() - 1
}
util.rs
// 公式ツールを一部改変したもの
#![allow(non_snake_case, unused_macros)]
use proconio::input;
use rand::prelude::*;
use std::ops::RangeBounds;
pub trait SetMinMax {
fn setmin(&mut self, v: Self) -> bool;
fn setmax(&mut self, v: Self) -> bool;
}
impl<T> SetMinMax for T
where
T: PartialOrd,
{
fn setmin(&mut self, v: T) -> bool {
*self > v && {
*self = v;
true
}
}
fn setmax(&mut self, v: T) -> bool {
*self < v && {
*self = v;
true
}
}
}
#[macro_export]
macro_rules! mat {
($($e:expr),*) => { Vec::from(vec![$($e),*]) };
($($e:expr,)*) => { Vec::from(vec![$($e),*]) };
($e:expr; $d:expr) => { Vec::from(vec![$e; $d]) };
($e:expr; $d:expr $(; $ds:expr)+) => { Vec::from(vec![mat![$e $(; $ds)*]; $d]) };
}
const MAX_T: usize = 5000;
#[derive(Clone, Debug)]
pub struct Input {
pub eps: f64,
pub delta: f64,
pub s: (i64, i64),
pub ps: Vec<(i64, i64)>,
pub walls: Vec<(i64, i64, i64, i64)>,
pub fs: Vec<(i64, i64)>,
pub alphas: Vec<f64>,
}
impl std::fmt::Display for Input {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
writeln!(f, "{} {} {:.2} {:.2}", self.ps.len(), self.walls.len(), self.eps, self.delta)?;
writeln!(f, "{} {}", self.s.0, self.s.1)?;
for i in 0..self.ps.len() {
writeln!(f, "{} {}", self.ps[i].0, self.ps[i].1)?;
}
for i in 0..self.walls.len() {
writeln!(
f,
"{} {} {} {}",
self.walls[i].0, self.walls[i].1, self.walls[i].2, self.walls[i].3
)?;
}
for i in 0..MAX_T {
writeln!(f, "{}", self.alphas[i])?;
}
for i in 0..MAX_T {
writeln!(f, "{} {}", self.fs[i].0, self.fs[i].1)?;
}
Ok(())
}
}
pub fn parse_input(f: &str) -> Input {
let f = proconio::source::once::OnceSource::from(f);
input! {
from f,
N: usize, M: usize, eps: f64, delta: f64,
s: (i64, i64),
ps: [(i64, i64); N],
walls: [(i64, i64, i64, i64); M],
alphas: [f64; MAX_T],
fs: [(i64, i64); MAX_T],
}
Input {
eps,
delta,
s,
ps,
walls,
fs,
alphas,
}
}
pub fn read<T: Copy + PartialOrd + std::fmt::Display + std::str::FromStr, R: RangeBounds<T>>(
token: Option<&str>,
range: R,
) -> Result<T, String> {
if let Some(v) = token {
if let Ok(v) = v.parse::<T>() {
if !range.contains(&v) {
Err(format!("Out of range: {}", v))
} else {
Ok(v)
}
} else {
Err(format!("Parse error: {}", v))
}
} else {
Err("Unexpected EOF".to_owned())
}
}
// ビジュアライズ用に変更
pub struct Output {
pub out: Vec<(char, i64, i64)>,
pub pos: Vec<(i64, i64)>,
pub vel: Vec<(i64, i64)>,
pub est: Vec<Vec<(i64, i64)>>,
pub dst: Vec<i64>,
pub gol: Vec<Vec<bool>>,
}
pub fn parse_output(_input: &Input, f: &str) -> Result<Output, String> {
let mut out = Vec::new();
let mut pos = Vec::new();
let mut vel = Vec::new();
let mut est = vec![Vec::new(); 5000];
let mut dst = vec![-1; 5000];
let mut gol = vec![vec![false; _input.ps.len()]; 5000];
let mut turn = 0;
for line in f.lines() {
if line.starts_with('#') {
// ビジュアライズ用に追加
let mut it = line.split_whitespace();
let a = it.next().unwrap();
if a.len() == 1 {
continue;
}
match a.chars().nth(1).unwrap() {
'p' => {
let x = read(it.next(), -100000..=100000)?;
let y = read(it.next(), -100000..=100000)?;
pos.push((x, y));
},
'v' => {
let x = read(it.next(), i64::MIN..=i64::MAX)?;
let y = read(it.next(), i64::MIN..=i64::MAX)?;
vel.push((x, y));
},
'e' => {
if turn == 0 {
continue;
}
let x = read(it.next(), -100000..=100000)?;
let y = read(it.next(), -100000..=100000)?;
est[turn - 1].push((x, y));
},
'd' => {
let d = read(it.next(), 0..300000)?;
if turn == 0 {
return Err(format!("Distance before measure"));
}
dst[turn - 1] = d;
},
'g' => {
let g = read(it.next(), 0.._input.ps.len())?;
for t in turn..5000 {
gol[t][g] = true;
}
}
_ => return Err(format!("Invalid comment: {}", a)),
}
continue;
}
let mut it = line.split_whitespace();
let a = read(it.next(), 'A'..'Z')?;
let x = read(it.next(), -100000..=100000)?;
let y = read(it.next(), -100000..=100000)?;
if a != 'A' && a != 'S' {
return Err(format!("Invalid action: {}", a));
} else if a == 'A' && x * x + y * y > 500 * 500 {
return Err(format!("Out of range: ({}, {})", x, y));
} else if a == 'S' && x * x + y * y > 10000000000 {
return Err(format!("Out of range: ({}, {})", x, y));
} else if a == 'S' && (x, y) == (0, 0) {
return Err(format!("Out of range: ({}, {})", x, y));
}
out.push((a, x, y));
turn += 1;
}
if out.len() > MAX_T {
return Err(format!("Too many actions: {}", out.len()));
}
Ok(Output { out, pos, vel, dst, gol, est })
}
pub fn gen(seed: u64, problem: char) -> Input {
let mut rng = rand_chacha::ChaCha20Rng::seed_from_u64(seed);
let N = 10;
let (M, eps, delta) = match problem {
'A' => (0, rng.gen_range(1..=100) as f64, rng.gen_range(1..=20) as f64 * 0.01),
'B' => (10, rng.gen_range(0..=1) as f64, 0.01),
'C' => (
rng.gen_range(1..=10i32) as usize,
rng.gen_range(1..=100) as f64,
rng.gen_range(1..=20) as f64 * 0.01,
),
_ => {
panic!("Unknown problem: {}", problem)
}
};
let s = (rng.gen_range(-99999..=99999), rng.gen_range(-99999..=99999));
let mut ps: Vec<(i64, i64)> = vec![];
while ps.len() < N {
let p = (rng.gen_range(-100000..=100000), rng.gen_range(-100000..=100000));
if ps
.iter()
.chain(&[s])
.any(|&q| (p.0 - q.0) * (p.0 - q.0) + (p.1 - q.1) * (p.1 - q.1) < 5000 * 5000)
{
continue;
}
ps.push(p);
}
let mut walls: Vec<(i64, i64, i64, i64)> = vec![];
while walls.len() < M {
let x1 = rng.gen_range(-90000..=90000);
let y1 = rng.gen_range(-90000..=90000);
let x2 = x1 + rng.gen_range(-100000..=100000);
let y2 = y1 + rng.gen_range(-100000..=100000);
if (x2 < -100000 || 100000 < x2) && (y2 < -100000 || 100000 < y2) || (x2, y2) == (0, 0) {
continue;
}
let x2 = x2.min(100000).max(-100000);
let y2 = y2.min(100000).max(-100000);
if walls.iter().all(|w| {
!P::crs_ss(
(P(x1 as f64, y1 as f64), P(x2 as f64, y2 as f64)),
(P(w.0 as f64, w.1 as f64), P(w.2 as f64, w.3 as f64)),
)
}) {
if !P::crs_sp((P(x1 as f64, y1 as f64), P(x2 as f64, y2 as f64)), P(s.0 as f64, s.1 as f64)) {
walls.push((x1, y1, x2, y2));
}
}
}
let alphas = (0..MAX_T)
.map(|_| loop {
let t = 1.0 + rng.sample::<f64, _>(rand_distr::StandardNormal) * delta;
if t > 0.0 {
break t;
}
})
.collect();
let fs = (0..MAX_T)
.map(|_| {
(
(rng.sample::<f64, _>(rand_distr::StandardNormal) * eps).round() as i64,
(rng.sample::<f64, _>(rand_distr::StandardNormal) * eps).round() as i64,
)
})
.collect();
Input {
eps,
delta,
s,
ps,
walls,
alphas,
fs,
}
}
pub fn compute_score(input: &Input, out: &Output) -> (i64, String) {
let (mut score, mut err, (_, _, visited)) = compute_score_details(input, &out.out);
if visited.iter().any(|&b| !b) {
err = "Unexpected EOF".to_owned();
}
if err.len() > 0 {
score = 0;
}
(score, err)
}
struct Sim {
visited: Vec<bool>,
score: i64,
crt_score: i64,
p: P,
v: P,
t: usize,
}
impl Sim {
fn new(input: &Input) -> Self {
let visited = vec![false; input.ps.len()];
let score = 0;
let crt_score = 0;
let p = P(input.s.0 as f64, input.s.1 as f64);
let v = P(0.0, 0.0);
Self {
visited,
score,
crt_score,
p,
v,
t: 0,
}
}
fn query(&mut self, input: &Input, mv: char, x: i64, y: i64) -> (i32, Vec<usize>, i64) {
let mut ret = -1;
match mv {
'A' => {
self.v = self.v + P(x as f64, y as f64);
}
'S' => {
let mut d = 1e9;
for wall in input.walls.iter().chain(
[
(-100000, -100000, -100000, 100000),
(-100000, 100000, 100000, 100000),
(100000, 100000, 100000, -100000),
(100000, -100000, -100000, -100000),
]
.iter(),
) {
let dir = P(x as f64, y as f64);
let w1 = P(wall.0 as f64, wall.1 as f64);
let w2 = P(wall.2 as f64, wall.3 as f64);
if let Some(p) = P::pi_ll((self.p, self.p + dir), (w1, w2)) {
if sig(dir.det(w1 - self.p)) * sig(dir.det(w2 - self.p)) <= 0 && (p - self.p).dot(dir) >= 0.0 {
d.setmin((p - self.p).abs2().sqrt());
}
}
}
d *= input.alphas[self.t];
ret = d.round() as i64;
}
_ => {
unreachable!()
}
}
self.v = self.v + P(input.fs[self.t].0 as f64, input.fs[self.t].1 as f64);
self.crt_score -= 2;
self.t += 1;
let q = self.p + self.v;
if q.0 < -100000.0
|| 100000.0 < q.0
|| q.1 < -100000.0
|| 100000.0 < q.1
|| input
.walls
.iter()
.any(|&(x1, y1, x2, y2)| P::crs_ss((P(x1 as f64, y1 as f64), P(x2 as f64, y2 as f64)), (self.p, q)))
{
self.crt_score -= 100;
self.v = P(0.0, 0.0);
return (1, vec![], ret);
} else {
let mut hit = vec![];
for i in 0..input.ps.len() {
if !self.visited[i] && P::dist2_sp((self.p, q), P(input.ps[i].0 as f64, input.ps[i].1 as f64)) <= 1000000.0 {
self.visited[i] = true;
self.crt_score += 1000;
hit.push(i);
}
}
self.p = q;
self.score.setmax(self.crt_score);
(0, hit, ret)
}
}
}
pub fn compute_score_details(input: &Input, out: &[(char, i64, i64)]) -> (i64, String, (P, P, Vec<bool>)) {
let mut sim = Sim::new(input);
for &(mv, x, y) in out {
sim.query(input, mv, x, y);
}
(sim.score, String::new(), (sim.p, sim.v, sim.visited))
}
use std::cmp::Ordering;
use std::ops::*;
#[derive(Clone, Copy, Default, Debug, PartialEq, PartialOrd)]
pub struct P(pub f64, pub f64);
impl Add for P {
type Output = P;
fn add(self, a: P) -> P {
P(self.0 + a.0, self.1 + a.1)
}
}
impl Sub for P {
type Output = P;
fn sub(self, a: P) -> P {
P(self.0 - a.0, self.1 - a.1)
}
}
impl Mul<f64> for P {
type Output = P;
fn mul(self, a: f64) -> P {
P(self.0 * a, self.1 * a)
}
}
impl P {
pub fn dot(self, a: P) -> f64 {
(self.0 * a.0) + (self.1 * a.1)
}
pub fn det(self, a: P) -> f64 {
(self.0 * a.1) - (self.1 * a.0)
}
pub fn abs2(self) -> f64 {
self.dot(self)
}
}
fn sig<T>(x: T) -> i32
where
T: Default + PartialOrd,
{
match x.partial_cmp(&T::default()) {
Some(Ordering::Greater) => 1,
Some(Ordering::Less) => -1,
_ => 0,
}
}
impl P {
pub fn dist2_sp((p1, p2): (P, P), q: P) -> f64 {
if (p2 - p1).dot(q - p1) <= 0.0 {
(q - p1).abs2()
} else if (p1 - p2).dot(q - p2) <= 0.0 {
(q - p2).abs2()
} else {
P::dist2_lp((p1, p2), q)
}
}
pub fn dist2_lp((p1, p2): (P, P), q: P) -> f64 {
let det = (p2 - p1).det(q - p1);
det * det / (p2 - p1).abs2()
}
pub fn crs_sp((p1, p2): (P, P), q: P) -> bool {
P::crs_lp((p1, p2), q) && (q - p1).dot(q - p2) <= 0.0
}
pub fn crs_lp((p1, p2): (P, P), q: P) -> bool {
(p2 - p1).det(q - p1) == 0.0
}
pub fn crs_ss((p1, p2): (P, P), (q1, q2): (P, P)) -> bool {
let sort = |a, b| {
if a < b {
(a, b)
} else {
(b, a)
}
};
let (lp0, up0) = sort(p1.0, p2.0);
let (lq0, uq0) = sort(q1.0, q2.0);
let (lp1, up1) = sort(p1.1, p2.1);
let (lq1, uq1) = sort(q1.1, q2.1);
if up0 < lq0 || uq0 < lp0 || up1 < lq1 || uq1 < lp1 {
return false;
}
return sig((p2 - p1).det(q1 - p1)) * sig((p2 - p1).det(q2 - p1)) <= 0
&& sig((q2 - q1).det(p1 - q1)) * sig((q2 - q1).det(p2 - q1)) <= 0;
}
pub fn pi_ll((p1, p2): (P, P), (q1, q2): (P, P)) -> Option<P> {
let d = (q2 - q1).det(p2 - p1);
if d == 0.0 {
return None;
}
let r = p1 * d + (p2 - p1) * (q2 - q1).det(q1 - p1);
Some(P(r.0 / d, r.1 / d))
}
}
// ------------------------------------------------------------------------------
// ビジュアライザ
use svg::node::element::{Circle, Line, Rectangle,};
const MARGIN: usize = 5;
const CANVAS_SIZE: usize = 500;
const STROKE_WIDTH: usize = 1;
const POINT_RADIUS: usize = 3;
const AREA_SIZE: i64 = 100_000;
// (score, err, svg)
pub fn vis(input: &Input, output: &Output, turn: usize) -> (i64, String, String) {
let (score, err) = compute_score(input, output);
let mut doc = svg::Document::new()
.set("id", "vis")
.set("viewBox", (-(MARGIN as i64), -(MARGIN as i64), to_canvas(AREA_SIZE) + 2 * MARGIN, to_canvas(AREA_SIZE) + 2 * MARGIN))
.set("width", to_canvas(AREA_SIZE) + MARGIN)
.set("height", to_canvas(AREA_SIZE) + MARGIN)
.set("style", "background-color:white");
// 目的地
for g in 0..input.ps.len() {
let (x, y) = input.ps[g];
let color = if output.gol[turn][g] { "lime" } else { "red" };
doc = doc.add(make_circle(to_canvas(x), to_canvas(y), POINT_RADIUS, color));
}
// 壁
let corners = [(-AREA_SIZE, -AREA_SIZE), (-AREA_SIZE, AREA_SIZE), (AREA_SIZE, AREA_SIZE), (AREA_SIZE, -AREA_SIZE)];
for i in 0..4 {
let (x1, y1) = corners[i];
let (x2, y2) = corners[(i + 1) % 4];
doc = doc.add(make_line(to_canvas(x1), to_canvas(y1), to_canvas(x2), to_canvas(y2), "black"));
}
for (x1, y1, x2, y2) in input.walls.iter() {
doc = doc.add(make_line(to_canvas(*x1), to_canvas(*y1), to_canvas(*x2), to_canvas(*y2), "black"));
}
// ドローンの軌跡
for t in 0..=turn {
let (x, y) = output.pos[t];
doc = doc.add(make_circle(to_canvas(x), to_canvas(y), 1, "grey"));
}
let (px, py) = output.pos[turn];
// 測定
if output.out[turn].0 == 'S' && output.dst[turn] >= 0 {
let dx = output.out[turn].1 as f64;
let dy = output.out[turn].2 as f64;
let c = output.dst[turn] as f64 / (dx * dx + dy * dy).sqrt();
let x = px + (c * dx) as i64;
let y = py + (c * dy) as i64;
doc = doc.add(make_line(to_canvas(px), to_canvas(py), to_canvas(x), to_canvas(y), "navy"));
}
// 推定
for (x, y) in output.est[turn].iter() {
doc = doc.add(make_circle(to_canvas(*x), to_canvas(*y), 1, "teal"));
}
// ドローン
doc = doc.add(make_circle(to_canvas(px), to_canvas(py), POINT_RADIUS, "blue"));
(score, err, doc.to_string())
}
pub fn to_canvas(x: i64) -> usize {
(x + AREA_SIZE) as usize / (2 * AREA_SIZE as usize / CANVAS_SIZE)
}
#[allow(dead_code)]
pub fn make_circle(x: usize, y: usize, r: usize, fill: &str) -> Circle {
Circle::new()
.set("cx", x)
.set("cy", y)
.set("r", r)
.set("fill", fill)
}
#[allow(dead_code)]
pub fn make_line(x1: usize, y1: usize, x2: usize, y2: usize, color: &str) -> Line {
Line::new()
.set("x1", x1)
.set("y1", y1)
.set("x2", x2)
.set("y2", y2)
.set("stroke", color)
.set("stroke-width", STROKE_WIDTH)
.set("stroke-linecap", "round")
}
#[allow(dead_code)]
pub fn make_rectangle(x: usize, y: usize, w: usize, h: usize, fill: &str) -> Rectangle {
Rectangle::new()
.set("x", x)
.set("y", y)
.set("width", w)
.set("height", h)
.set("fill", fill)
}
Cargo.toml
[package]
name = "rust"
version = "0.1.0"
edition = "2021"
# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
[lib]
crate-type = ["cdylib"]
[dependencies]
wasm-bindgen = "0.2.89"
rand = { version = "0.8.5", features = ["small_rng", "min_const_gen"] }
rand_chacha = "0.3.1"
rand_distr = "0.4.3"
itertools = "0.11.0"
proconio = { version = "0.3.6", features = ["derive"] }
clap = { version = "4.0.22", features = ["derive"] }
svg = "0.9.0"
delaunator = "1.0.1"
web-sys = { version = "0.3.44", features = ["console"] }
getrandom = { version = "0.2", features = ["js"] }
最後に
そういえば、まだ焼肉を奢ってもらっていないですね。 2
AHC030 参加記
THIRD プログラミングコンテスト2023（AtCoder Heuristic Contest 030） お疲れ様でした。
システムテストの得点が 2,050,795,126,975 で 3 位でした。
問題概要
クエリを通して、石油がある場所を特定する問題です。
油田の個数と形は入力で与えられますが、各油田がどこにあるかは与えられません。
詳しくは公式の問題文を参考にしてください。
方針
毎ターン次のことを行いました。
- 焼きなまし法を用いて、油田配置の中で対数尤度が大きいものを 32 個ほど得る。
- 得られた油田配置の中に実際の配置が含まれていると仮定し、各油田配置の事後確率を求める。
- 得られた事後分布を利用して占うマスを選ぶ。
他の上位者と大まかな流れは一緒だと思います。焼きなまし法の部分が私は得意だったようなので、この記事では焼きなまし法で油田配置を推定する部分を中心に説明します。
ベイズ推定
占う場所を選ぶときに油田配置の事後分布が欲しくなるので、油田配置の事後分布を求めることを目標とします。
占いの履歴を とし、各油田配置を
と表記します。
が得られたときに実際の油田配置が
である事後確率
を求めたいです。入力生成方法を読むと事前確率
が全て等しいことがわかります。ベイズの定理を踏まえると次の式が成り立ちます。
を求める代わりに
を求めればよいことがわかりました。
です。
は正規分布の一部を積分した値であり、正規分布の累積分布関数が
であることから計算できます。
油田の個数 がある程度大きいとき、全ての油田配置
について
を計算する時間がありません。そこで、
が比較的大きいものをいくつか集めた集合を
とし、
に実際の油田配置が含まれていると仮定して
で分布を近似することにしました。
と過去の占いの情報量が十分大きければ、高い精度で近似できると思います。
焼きなまし法
を求める手法はいくつか考えられ、大きく 2 つに分類できると思います。1 種類目は枝刈り全探索、MCTS、ビームサーチなどの全探索を改善したもので、2 種類目は山登り法や焼きなまし法などの局所探索です。
結論として、今回は焼きなまし法を使用しました。 が大きいときに探索空間が広すぎて全探索系の手法が弱そうで、また、局所性がある近傍を設定して焼きなますことができそうだったからです。
焼きなまし法の評価値として、対数尤度 を使用しました。
評価値が高い状態を 32 個ほど優先度付きキューで管理し、最後に保持していたものを としました。同じ状態やクエリで推測して間違えた配置は省くようにします。
近傍は 3 種類です。
近傍 1: 1 マス移動（選択率: 9 / 16）
ランダムに 1 つの油田を選び、前後左右のランダムな方向に 1 マスずらします。
近傍 2: ランダムな場所に移動（選択率: 3 / 16）
ランダムに 1 つの油田を選び、ランダムに選んだ場所に移動させます。
近傍 3: 2 つの油田の入れ替え（選択率: 4 / 16）
ランダムに 2 つの油田を選び、2 つの油田の位置を入れ替えます。
石油埋蔵量が変化するマスの数が最も少なくなるような入れ替えのみを近傍として設定しました。複数通りある場合はそのいずれかをランダムに選びます。入れ替え方は入力を読み込んだ直後に事前計算します。
状態の持ち方
状態は、各油田の左上の座標と、各占いに対する石油埋蔵量を保持するようにしました。 の二次元グリッドは使いません。
油田が 個あって、それぞれの置き方はたかだか
通りしかありません。すなわち、1 つの油田を選んで配置する場合の数は
です。そこで、占いの結果を得るたびに
通りそれぞれに対して占った場所の個数を計算するようにしました。この処理によって石油埋蔵量を
で更新できます。
また、占いの結果を得るたびに、実際の石油埋蔵量が だったときにその占いの観測値が得られる確率も計算しておきます。焼きなまし法の最中に誤差関数を呼ぶ必要がなくなり、尤度の計算を
で実現できます。
実装では石油埋蔵量が変わる占いのみを取り出して更新していますが、石油埋蔵量が変わる占いの割合が高そうなので高速化されているかはわかりません。
多点スタート
8 つの状態を並行して焼きなまし、少しずつ状態を減らすようにしました。
状態を減らすとき、基本的にはその時点での評価値が最も低いものが削除されます。ただし、評価値の悪化を許容した直後に削除されるようなケースもあるため、評価値の指数移動平均をとったり、遷移の受容確率を評価したりしました。改善の寄与は小さいです。
初期状態
焼きなまし法の 1 ターン目の初期状態はランダムに生成しますが、2 ターン目以降の初期状態として、1 つ前のターンで得られた油田配置の中から尤度が高いものを選ぶようにしました。
初期解の尤度が高いと焼きなまし法の収束が速くなり、探索性能が改善されるようでした。
8 つの状態から 32 個の状態を生成し、次のターンは 32 個の状態から評価値が高い 8 つの状態を選ぶため、焼きなまし法を遷移としたビームサーチと捉えることもできます。
温度の自動調節
基本的には、初期温度を 6.0、最終温度を 0.6 にして温度を指数的に減少させます。
ただし、温度を上記のように固定すると、占いの回数が増えるにつれて焼きなまし法の遷移採択率が低くなるようだったので、採択率が 0.04 を下回ったときに初期温度と最終温度を 1.05 倍しました。採択率が 0.04 以上になるように温度を自動調節しているつもりです。難しいケースだと初期温度が 10.0 程度まで上がるようでした。
近傍における油田の選び方
近傍の説明でランダムに油田を選ぶと書きましたが、油田を等しい確率で選ぶのではなく、過去の遷移における採択率が高かったものを優先して選びました。
改善の寄与は小さいです。
対数尤度を整数で表現する
浮動小数点演算よりも整数演算のほうが速いです。そこで、対数尤度に をかけたものを 64 ビット整数で管理しました。焼きなまし法の中で浮動小数点数は扱いません。
速度を測っていないので確実なことはわからないですが、あまり速くなっていない気がします。
占うマスの決め方
最尤油田配置とそうでない油田配置を分離することを考えました。
焼きなまし法で求めた擬似的な事後分布から、各マスの石油埋蔵量の期待値 を求め、最尤油田配置における
が
以上なら占うマスとして暫定的に採用します。
他の候補同士も分離したかったので、適当に評価関数を作って山登り法で占うマスの集合を改善しているつもりですが、評価関数の質が悪かったからか、思うようには改善されませんでした。
感想戦でエントロピー最小化（占いとしては相互情報量の最大化）をすればよいことを知って納得しました。情報理論が身についていないのが敗因です。
序盤の占い方
最初の数ターンはルールベースで占うマスを決めました。事後分布は使いません。
が小さいときには損だったかもしれません。
時間管理
必要なターン数を雑に見積もったうえで、序盤のターンほど時間が長くなるように調整しました。
非常に難しいケースでは各ターンが使う時間を等しくしました。
提出コード
最後に
天才的なアイディアが必要というよりは、理論とヒューリスティック最適化の基礎的な理解と応用力が問われる問題だと思っていて、かなり自分好みの問題でした。
優勝できなかったのは悔しいですが、情報理論の大切さや MCMC との関連性など、多くのことを学べてよかったです。
コンテストを主催してくださった THIRD 様、writer の wata さん、コンテスト参加者の皆様、最後まで読んでくださった読者の皆様、ありがとうございました。
差分更新ビームサーチライブラリの実装 (C++)
この記事について
- 差分更新型ビームサーチライブラリの実装例について説明します。
- 差分更新型のビームサーチについては 高速なビームサーチが欲しい!!! などで既に解説されているため、被る部分については詳しく説明しません。
- この記事に書いたソースコードをプログラミングコンテストで自由に使っていただいて構いませんが、当ブログは損害などに対する責任を負いかねますのでご了承ください。
前提知識
差分更新型のビームサーチとは
木上のビームサーチ、Euler Tour ビームサーチとも呼ばれています1。
ビームサーチは、幅優先探索に枝刈りを取り入れた手法です。探索木における深さ の頂点集合から深さ
の頂点集合を生成し、その中から評価値が高い上位
個を選択します。
のことをビーム幅と呼びます。
、赤色が探索中のノード、橙色が採用されたノード、灰色が不採用のノード）
愚直なビームサーチでは、頂点ごとに独立した状態を作成します。実装は比較的楽だと思いますが、状態や履歴をコピーする必要が出てくるため、実行速度が遅くなりやすいです2。
そこで、探索木を明示的に作成し、Euler Tour の順序で1つの状態を更新するようにしたものが差分更新型のビームサーチです。探索木の葉を訪れたときに新しい葉の候補を生成します。
状態遷移の履歴は探索木から復元できるため、状態のコピーだけでなく履歴のコピーも省略できます。
深さが の1つの葉を探索するのに最悪で
回の遷移を必要としますが、実際には多くの葉で近い先祖を共有するため、遷移回数は平均して
回よりもかなり少なくなります。
差分更新型ビームサーチは愚直なビームサーチよりも高速に動作する場合が多く3、実装がやや複雑なため、ライブラリを作成しました。
実装
ライブラリ全体を折りたたみに記載しました。コードを読む場合は下の方にある beam_search 関数を最初に読むと分かりやすいかもしれません。
ビームサーチライブラリ
#include <bits/stdc++.h>
#include <atcoder/all>
using namespace std;
namespace beam_search {
// メモリの再利用を行いつつ集合を管理するクラス
template<class T>
class ObjectPool {
public:
// 配列と同じようにアクセスできる
T& operator[](int i) {
return data_[i];
}
// 配列の長さを変更せずにメモリを確保する
void reserve(size_t capacity) {
data_.reserve(capacity);
}
// 要素を追加し、追加されたインデックスを返す
int push(const T& x) {
if (garbage_.empty()) {
data_.push_back(x);
return data_.size() - 1;
} else {
int i = garbage_.top();
garbage_.pop();
data_[i] = x;
return i;
}
}
// 要素を（見かけ上）削除する
void pop(int i) {
garbage_.push(i);
}
// 使用した最大のインデックス(+1)を得る
// この値より少し大きい値をreserveすることでメモリの再割り当てがなくなる
size_t size() {
return data_.size();
}
private:
vector<T> data_;
stack<int> garbage_;
};
// 連想配列
// Keyにハッシュ関数を適用しない
// open addressing with linear probing
// unordered_mapよりも速い
// nは格納する要素数よりも4~16倍ほど大きくする
template <class Key, class T>
struct HashMap {
public:
explicit HashMap(uint32_t n) {
n_ = n;
valid_.resize(n_, false);
data_.resize(n_);
}
// 戻り値
// - 存在するならtrue、存在しないならfalse
// - index
pair<bool,int> get_index(Key key) const {
Key i = key % n_;
while (valid_[i]) {
if (data_[i].first == key) {
return {true, i};
}
if (++i == n_) {
i = 0;
}
}
return {false, i};
}
// 指定したindexにkeyとvalueを格納する
void set(int i, Key key, T value) {
valid_[i] = true;
data_[i] = {key, value};
}
// 指定したindexのvalueを返す
T get(int i) const {
assert(valid_[i]);
return data_[i].second;
}
void clear() {
fill(valid_.begin(), valid_.end(), false);
}
private:
uint32_t n_;
vector<bool> valid_;
vector<pair<Key,T>> data_;
};
using Hash = uint32_t; // TODO
// 状態遷移を行うために必要な情報
// メモリ使用量をできるだけ小さくしてください
struct Action {
// TODO
Action() {
// TODO
}
};
using Cost = int; // TODO
// 状態のコストを評価するための構造体
// メモリ使用量をできるだけ小さくしてください
struct Evaluator {
// TODO
Evaluator() {
// TODO
}
// 低いほどよい
Cost evaluate() const {
// TODO
}
};
// 展開するノードの候補を表す構造体
struct Candidate {
Action action;
Evaluator evaluator;
Hash hash;
int parent;
Cost cost;
Candidate(Action action, Evaluator evaluator, Hash hash, int parent, Cost cost) :
action(action),
evaluator(evaluator),
hash(hash),
parent(parent),
cost(cost) {}
};
// ビームサーチの設定
struct Config {
int max_turn;
size_t beam_width;
size_t nodes_capacity;
uint32_t hash_map_capacity;
};
// 削除可能な優先度付きキュー
using MaxSegtree = atcoder::segtree<
pair<Cost,int>,
[](pair<Cost,int> a, pair<Cost,int> b){
if (a.first >= b.first) {
return a;
} else {
return b;
}
},
[]() { return make_pair(numeric_limits<Cost>::min(), -1); }
>;
// ノードの候補から実際に追加するものを選ぶクラス
// ビーム幅の個数だけ、評価がよいものを選ぶ
// ハッシュ値が一致したものについては、評価がよいほうのみを残す
class Selector {
public:
explicit Selector(const Config& config) :
hash_to_index_(config.hash_map_capacity)
{
beam_width = config.beam_width;
candidates_.reserve(beam_width);
full_ = false;
st_original_.resize(beam_width);
}
// 候補を追加する
// ターン数最小化型の問題で、candidateによって実行可能解が得られる場合にのみ finished = true とする
// ビーム幅分の候補をCandidateを追加したときにsegment treeを構築する
void push(Action action, const Evaluator& evaluator, Hash hash, int parent, bool finished) {
Cost cost = evaluator.evaluate();
if (finished) {
finished_candidates_.emplace_back(Candidate(action, evaluator, hash, parent, cost));
return;
}
if (full_ && cost >= st_.all_prod().first) {
// 保持しているどの候補よりもコストが小さくないとき
return;
}
auto [valid, i] = hash_to_index_.get_index(hash);
if (valid) {
int j = hash_to_index_.get(i);
if (hash == candidates_[j].hash) {
// ハッシュ値が等しいものが存在しているとき
if (cost < candidates_[j].cost) {
// 更新する場合
candidates_[j] = Candidate(action, evaluator, hash, parent, cost);
if (full_) {
st_.set(j, {cost, j});
}
}
return;
}
}
if (full_) {
// segment treeが構築されている場合
int j = st_.all_prod().second;
hash_to_index_.set(i, hash, j);
candidates_[j] = Candidate(action, evaluator, hash, parent, cost);
st_.set(j, {cost, j});
} else {
// segment treeが構築されていない場合
hash_to_index_.set(i, hash, candidates_.size());
candidates_.emplace_back(Candidate(action, evaluator, hash, parent, cost));
if (candidates_.size() == beam_width) {
// 保持している候補がビーム幅分になったとき
construct_segment_tree();
}
}
}
// 選んだ候補を返す
const vector<Candidate>& select() const {
return candidates_;
}
// 実行可能解が見つかったか
bool have_finished() const {
return !finished_candidates_.empty();
}
// 実行可能解に到達する「候補」を返す
vector<Candidate> get_finished_candidates() const {
return finished_candidates_;
}
void clear() {
candidates_.clear();
hash_to_index_.clear();
full_ = false;
}
private:
size_t beam_width;
vector<Candidate> candidates_;
HashMap<Hash,int> hash_to_index_;
bool full_;
vector<pair<Cost,int>> st_original_;
MaxSegtree st_;
vector<Candidate> finished_candidates_;
void construct_segment_tree() {
full_ = true;
for (size_t i = 0; i < beam_width; ++i) {
st_original_[i] = {candidates_[i].cost, i};
}
st_ = MaxSegtree(st_original_);
}
};
// 深さ優先探索に沿って更新する情報をまとめたクラス
class State {
public:
explicit State(/* const Input& input */) {
// TODO
}
// 次の状態候補を全てselectorに追加する
// 引数
// evaluator : 今の評価器
// hash : 今のハッシュ値
// parent : 今のノードID（次のノードにとって親となる）
void expand(const Evaluator& evaluator, Hash hash, int parent, Selector& selector) {
// TODO
}
// actionを実行して次の状態に遷移する
void move_forward(Action action) {
// TODO
}
// actionを実行する前の状態に遷移する
// 今の状態は、親からactionを実行して遷移した状態である
void move_backward(Action action) {
// TODO
}
private:
// TODO
};
// 探索木（二重連鎖木）のノード
struct Node {
Action action;
Evaluator evaluator;
Hash hash;
int parent, child, left, right;
// 根のコンストラクタ
Node(Action action, const Evaluator& evaluator, Hash hash) :
action(action),
evaluator(evaluator),
hash(hash),
parent(-1),
child(-1),
left(-1),
right(-1) {}
// 通常のコンストラクタ
Node(const Candidate& candidate, int right) :
action(candidate.action),
evaluator(candidate.evaluator),
hash(candidate.hash),
parent(candidate.parent),
child(-1),
left(-1),
right(right) {}
};
// 二重連鎖木に対する操作をまとめたクラス
class Tree {
public:
explicit Tree(const State& state, size_t nodes_capacity, const Node& root) :
state_(state)
{
nodes_.reserve(nodes_capacity);
root_ = nodes_.push(root);
}
// 状態を更新しながら深さ優先探索を行い、次のノードの候補を全てselectorに追加する
void dfs(Selector& selector) {
update_root();
int v = root_;
while (true) {
v = move_to_leaf(v);
state_.expand(nodes_[v].evaluator, nodes_[v].hash, v, selector);
v = move_to_ancestor(v);
if (v == root_) {
break;
}
v = move_to_right(v);
}
}
// 根からノードvまでのパスを取得する
vector<Action> get_path(int v) {
// cerr << nodes_.size() << endl;
vector<Action> path;
while (nodes_[v].parent != -1) {
path.push_back(nodes_[v].action);
v = nodes_[v].parent;
}
reverse(path.begin(), path.end());
return path;
}
// 新しいノードを追加する
int add_leaf(const Candidate& candidate) {
int parent = candidate.parent;
int sibling = nodes_[parent].child;
int v = nodes_.push(Node(candidate, sibling));
nodes_[parent].child = v;
if (sibling != -1) {
nodes_[sibling].left = v;
}
return v;
}
// ノードvに子がいなかった場合、vと不要な先祖を削除する
void remove_if_leaf(int v) {
if (nodes_[v].child == -1) {
remove_leaf(v);
}
}
// 最も評価がよいノードを返す
int get_best_leaf(const vector<int>& last_nodes) {
assert(!last_nodes.empty());
int ret = last_nodes[0];
for (int v : last_nodes) {
if (nodes_[v].evaluator.evaluate() < nodes_[ret].evaluator.evaluate()) {
ret = v;
}
}
return ret;
}
private:
State state_;
ObjectPool<Node> nodes_;
int root_;
// 根から一本道の部分は往復しないようにする
void update_root() {
int child = nodes_[root_].child;
while (child != -1 && nodes_[child].right == -1) {
root_ = child;
state_.move_forward(nodes_[child].action);
child = nodes_[child].child;
}
}
// ノードvの子孫で、最も左にある葉に移動する
int move_to_leaf(int v) {
int child = nodes_[v].child;
while (child != -1) {
v = child;
state_.move_forward(nodes_[child].action);
child = nodes_[child].child;
}
return v;
}
// ノードvの先祖で、右への分岐があるところまで移動する
int move_to_ancestor(int v) {
while (v != root_ && nodes_[v].right == -1) {
state_.move_backward(nodes_[v].action);
v = nodes_[v].parent;
}
return v;
}
// ノードvの右のノードに移動する
int move_to_right(int v) {
state_.move_backward(nodes_[v].action);
v = nodes_[v].right;
state_.move_forward(nodes_[v].action);
return v;
}
// 不要になった葉を再帰的に削除する
void remove_leaf(int v) {
while (true) {
int left = nodes_[v].left;
int right = nodes_[v].right;
if (left == -1) {
int parent = nodes_[v].parent;
if (parent == -1) {
cerr << "ERROR: root is removed" << endl;
exit(-1);
}
nodes_.pop(v);
nodes_[parent].child = right;
if (right != -1) {
nodes_[right].left = -1;
return;
}
v = parent;
} else {
nodes_.pop(v);
nodes_[left].right = right;
if (right != -1) {
nodes_[right].left = left;
}
return;
}
}
}
};
// ビームサーチを行う関数
vector<Action> beam_search(const Config& config, State state, Node root) {
Tree tree(state, config.nodes_capacity, root);
// 探索中のノード集合
vector<int> curr_nodes;
curr_nodes.reserve(config.beam_width);
// 本来は curr_nodes = {state.root_} とすべきだが, 省略しても問題ない
// 新しいノードの集合
vector<int> next_nodes;
next_nodes.reserve(config.beam_width);
// 新しいノード候補の集合
Selector selector(config);
for (int turn = 0; turn < config.max_turn; ++turn) {
// Euler Tour で selector に候補を追加する
tree.dfs(selector);
if (selector.have_finished()) {
// ターン数最小化型の問題で実行可能解が見つかったとき
Candidate candidate = selector.get_finished_candidates()[0];
vector<Action> ret = tree.get_path(candidate.parent);
ret.push_back(candidate.action);
return ret;
}
// 新しいノードを追加する
for (const Candidate& candidate : selector.select()) {
next_nodes.push_back(tree.add_leaf(candidate));
}
if (next_nodes.empty()) {
// 新しいノードがないとき
cerr << "ERROR: Failed to find any valid solution" << endl;
return {};
}
// 不要なノードを再帰的に削除する
for (int v : curr_nodes) {
tree.remove_if_leaf(v);
}
// ダブルバッファリングで配列を使い回す
swap(curr_nodes, next_nodes);
next_nodes.clear();
selector.clear();
}
// ターン数固定型の問題で全ターンが終了したとき
int best_leaf = tree.get_best_leaf(curr_nodes);
return tree.get_path(best_leaf);
}
} // namespace beam_search
それぞれの構造体やクラスについて見ていきます。
Object Pool
ビームサーチと直接の関係がないクラスです。
配列にオブジェクトを保存し、削除したオブジェクトの場所を再利用します。
また、std::vector などと同様に reserve でメモリを確保できるようにしました4。
最初は長さ4の空の配列とします。 3, 1, 4 を順に追加します。追加した場所である 0, 1, 2 を順に返します。 a[1] を削除します。 5 を追加します。a[1] と a[3] のどちらに追加してもよいのですが、最後に使用した場所を優先的に使用することにします。今回の例だと a[1] です。追加した場所である 1 を返します。 削除した場所のインデックスを保持することでこのような挙動を実装することができます。例を用いた説明
a[0]
a[1]
a[2]
a[3]
a[0]
a[1]
a[2]
a[3]
3
1
4
a[0]
a[1]
a[2]
a[3]
3
4
a[0]
a[1]
a[2]
a[3]
3
5
4
差分更新型のビームサーチでは、ノードの追加と削除を頻繁に繰り返します。Object Pool を使うことで次のようなメリットを享受できます。
- メモリの再利用により空間計算量を削減できる。
- メモリ上の連続した領域を使用するので、データがキャッシュに乗りやすくなる。
- 十分な大きさのメモリを最初に確保することで、メモリの再割り当てをなくすことができる。
HashMap
Selector でハッシュ値が重複した候補を除去するところで使います。
open addressing を使用し、インデックスが衝突したときは linear probing を行いました。（ハッシュ値に対する）ハッシュ関数がなく、挿入と一括削除しか行わないので、std::unordered_set よりも単純で高速に動作します。
Hash が整数型でないときは整数型に変換するか、HashMap を std::unordered_set に変更する必要があります。
ちなみに、元々私は std::unordered_set を使っていて、saharan さんのツイート を読んで参考にしたら速度が少し上がりました5。
Action
状態遷移に必要な情報をまとめます。ここでいう状態遷移は、親からの移動と親への移動の両方を指しています。
Evaluator
評価値を計算するための構造体です。
基本的な実装は次のようなものです。
struct Evaluator { Cost cost; Evaluator(Cost cost) : cost(cost) {} Cost evaluate() const { return cost; } };
Evaluator::evaluate でコストを返します。コストが低いほど採用されやすくなります。
補足: Action や Evaluator で何を保持すべきか
既に述べたように、状態遷移に必要な情報を Action で保持し、評価結果の比較に必要な情報を Evaluator で保持することを想定しています。
一方で、Action と Evaluator に、より多くの情報を持たせることもできます。
状態が変数 x をメンバとして保持し、状態を更新するときに毎度 x を更新するものとします。このとき、x を状態ではなく Action や Evaluator が保持するようにすれば、Euler Tour における x の更新を省略できます。状態の関数内で x を使用したいときには、Action や Evaluator のメンバにアクセスすればよいです。x の後退処理（Euler Tour における子から親への遷移）を実装する必要がなくなるというメリットもあります。
例えば、複数の評価項目があるときに、Evaluator で各評価項目の値を保持するということが考えられます。また、後退処理がなくなるため、浮動小数点数も扱いやすくなります6。
一方で、Action や Evaluator の使用メモリが小さいほどよいという側面もあるため、全ての変数を Action や Evaluator に保持すればいいわけではありません。更新が面倒で使用メモリが少ない変数だけを Action や Evaluator に保持するとよいと思います。
例えば、状態内で探索木の深さを管理する場合、深さを更新するときにインクリメントやデクリメントという非常に軽い処理しか行われないため、深さは状態に保持すればよいと思います。
ちなみに、Action と State を空にして Evaluator で全ての情報を保持するようにすると、愚直なビームサーチらしくなります。
Candidate
新しいノードの候補を表現します。
Selector
Candidate の中から新しいノードとして採用するものを選びます。採用された Candidate を元に新しいノードが追加されます。
新しいノードを生成してから不要なものを削除するという実装も考えられますが、木に対する操作は定数倍が重いため、Candidate 構造体を経由する実装になっています。
コードのコメントに書かれているように、ハッシュ値が一致した候補については評価が最もよいものに限定し、その中から評価がよい候補をビーム幅の個数だけ選びます。ハッシュ値の重複を検出するために自作の HashMap を使ったり、ソートをなくすために atcoder::segtree を削除可能な優先度付きキューとして使ったり、高速化を意識して実装しました。
ハッシュ値以外で多様性を確保したい場合、Selector を実装し直す必要があります。例えばスライドパズルなどで「現在のマスの座標が一致するものの中から上位 個を選ぶ」という場合には書き直す必要があります。
State
Euler Tour に沿って更新する情報をまとめたクラスです。問題ごとに各メソッドを実装する必要があります。
ビームサーチの最中に State がコピーされることはないため、空間計算量は大きくても構いません。一方で、各メソッドは速いほどよいです。
Tree
探索木を二重連鎖木で表現し、木に対する操作をまとめたクラスです。
- 親
- 1つ上の兄
- 1つ下の弟
- 最も上の子供
状態の更新順序は Euler Tour と一緒です。
一方で、二重連鎖木上では兄弟間を直接移動するようにします。
二重連鎖木のノードは配列を使用しないため、Euler Tour もノードの追加や削除も簡潔に実装できます。
不要なノードは全て削除します。不要なノードというのは、子ができなかった、あるいは子が全て削除されたノードのことです。
高速なビームサーチが欲しい!!! で紹介されているように、根から一本道の部分は反復しないようにします。
上図では探索したノードが全て描かれていますが、実際には灰色のノードは作成されず、さらに赤色のノードを子孫として持たない6つの橙色のノードは削除されていることに注意してください。
beam_search
ビームサーチを実行する関数です。ライブラリの外からこの関数を呼び出します。
使用例
TOYOTA Programming Contest 2023 Summer（AtCoder Heuristic Contest 021） の実装例を紹介します。
大まかな方針
番号が小さいボールから揃えます。
番号が最小のボール、あるいはその左上または右上のボールを、左上か右上に移動させます。
評価関数は次のように設定しました。小さいほうがよいです。
ハッシュ値は、揃えているボールの位置と、既に揃えたボールの位置の集合の2つから生成しました。
ライブラリの使い方を紹介することが目的なので、考察などは省略します。
ハッシュ関数 （2025/12/16 修正）
using Hash = uint32_t; constexpr Hash hash_mask = ((Hash(1) << 22) - Hash(1)) << 10; // 上位22ビット constexpr uint32_t xorshift32(uint32_t& x) { x ^= x << 13; x ^= x >> 17; x ^= x << 5; return x; } template <size_t N> constexpr array<uint32_t, N> make_position_hashes(uint32_t seed) { array<uint32_t, N> hashes{}; uint32_t x = seed; for (size_t i = 0; i < N; ++i) { hashes[i] = (xorshift32(x) & hash_mask); } return hashes; } constexpr array<uint32_t, n * n> position_hashes = make_position_hashes<n * n>(123456); // 下位10ビットをターゲット位置とする inline Hash update_target_position(Hash hash, int x, int y) { return (hash & hash_mask) | (x * n + y); } // 上位22ビットをソート済み位置のZobrist hashとする inline Hash update_sorted_position(Hash hash, int x, int y) { return hash ^ position_hashes[x * n + y]; }
下位10ビットで揃えているボールの位置を保持し、上位22ビットで既に揃えたボールの位置の集合の Zobrist hash を保持しました。
Action
struct Action { int xyxy; Action(int x1, int y1, int x2, int y2) { xyxy = x1 | (y1 << 8) | (x2 << 16) | (y2 << 24); } tuple<int,int,int,int> decode() const { return {xyxy & 255, (xyxy >> 8) & 255, (xyxy >> 16) & 255, xyxy >> 24}; } };
AHC021の場合、スワップ位置が分かれば盤面の更新が行えるため、スワップする2つの位置を保持しました。
4つの整数を1つの int にエンコードし、メモリ使用量を減らしています8。
Evaluator
using Cost = int; constexpr int target_coefficient = 600; struct Evaluator { int target_ball; int potential; Evaluator(int target_ball, int potential) : target_ball(target_ball), potential(potential) {} Cost evaluate() const { return potential - target_coefficient * target_ball; } };
2つの評価項目を保持しました。
State
更新部分だけ説明します。他の部分を読みたい場合は提出コードをご覧ください。
class State { public: void move_forward(Action action) { auto [x1, y1, x2, y2] = action.decode(); swap_balls(x1, y1, x2, y2); } void move_backward(Action action) { auto [x1, y1, x2, y2] = action.decode(); swap_balls(x1, y1, x2, y2); } private: vector<vector<int>> b_; array<pair<int,int>,m> positions_; void swap_balls(int x1, int y1, int x2, int y2) { int b1 = b_[x1][y1]; int b2 = b_[x2][y2]; b_[x1][y1] = b2; b_[x2][y2] = b1; positions_[b2] = {x1, y1}; positions_[b1] = {x2, y2}; } };
b_ は盤面を表しています。positions_[i] はボール i の位置を表します。 positions_[b_[x][y]] = {x, y} が成り立ちます。positions_ は、あるボールを揃えた後に次のボールの位置を得るときなどに使われます。
Evaluator が揃えたボールの数を保持するため、現在揃えているボールの番号を保持・更新する必要がないことに注意してください。
今回は2つのボールを入れ替えるだけなので move_forward と move_backward が等しくなっていますが、一般的には異なります。
提出コード
ビーム幅を3500に設定しました。
提出言語は C++ 20 (Clang 16.0.6) です。私のビームサーチライブラリ（の主に Selector）は GCC よりも Clang のほうが高速に動作するようでした。
他の使用例
ゲーム実況者Xの挑戦 の提出コードです。解説などは省略します。
ターンが飛ぶビームサーチの実装 (2025/03/16 追記)
1つの Action で2ターン以上後にも遷移できるビームサーチを実装しました。
ターンが等しいもの同士を比較して候補を選びます。
追加した候補とそれらの先祖を active に設定し、 active なノードだけ探索しています。
ターンを飛ばせない実装よりも一般化されていますが、速度に関して定数倍が重くなっているので、使い分けたほうがいいと思います。
（2025/07/29 追記）ノードの削除に関するバグを修正し、使用例を更新しました。各ノードに対して子が追加されうる最も遅いターンを管理し、ノードの削除が可能かどうかを判定するようにしました。
ビームサーチライブラリ
#include <bits/stdc++.h>
#include <atcoder/segtree>
using namespace std;
namespace beam_search {
// メモリの再利用を行いつつ集合を管理するクラス
template<class T>
class ObjectPool {
public:
// 配列と同じようにアクセスできる
T& operator[](int i) {
return data_[i];
}
// 配列の長さを変更せずにメモリを確保する
void reserve(size_t capacity) {
data_.reserve(capacity);
}
// 要素を追加し、追加されたインデックスを返す
int push(const T& x) {
if (garbage_.empty()) {
data_.push_back(x);
return data_.size() - 1;
} else {
int i = garbage_.top();
garbage_.pop();
data_[i] = x;
return i;
}
}
// 要素を（見かけ上）削除する
void pop(int i) {
garbage_.push(i);
}
// 使用した最大のインデックス(+1)を得る
// この値より少し大きい値をreserveすることでメモリの再割り当てがなくなる
size_t size() {
return data_.size();
}
private:
vector<T> data_;
stack<int> garbage_;
};
// 連想配列
// Keyにハッシュ関数を適用しない
// open addressing with linear probing
// unordered_mapよりも速い
// nは格納する要素数よりも4~16倍ほど大きくする
template <class Key, class T>
struct HashMap {
public:
explicit HashMap(uint32_t n) {
n_ = n;
valid_.resize(n_, false);
data_.resize(n_);
}
// 戻り値
// - 存在するならtrue、存在しないならfalse
// - index
pair<bool,int> get_index(Key key) const {
Key i = key % n_;
while (valid_[i]) {
if (data_[i].first == key) {
return {true, i};
}
if (++i == n_) {
i = 0;
}
}
return {false, i};
}
// 指定したindexにkeyとvalueを格納する
void set(int i, Key key, T value) {
valid_[i] = true;
data_[i] = {key, value};
}
// 指定したindexのvalueを返す
T get(int i) const {
assert(valid_[i]);
return data_[i].second;
}
void clear() {
fill(valid_.begin(), valid_.end(), false);
}
private:
uint32_t n_;
vector<bool> valid_;
vector<pair<Key,T>> data_;
};
using Hash = uint64_t; // TODO
// 状態遷移を行うために必要な情報
// メモリ使用量をできるだけ小さくしてください
struct Action {
// TODO
Action() {
// TODO
}
};
using Cost = int; // TODO
// 状態のコストを評価するための構造体
// メモリ使用量をできるだけ小さくしてください
struct Evaluator {
// TODO
Evaluator() {
// TODO
}
// 低いほどよい
Cost evaluate() const {
// TODO
}
};
// 展開するノードの候補を表す構造体
struct Candidate {
Action action;
Evaluator evaluator;
Hash hash;
int parent;
Cost cost;
Candidate(Action action, Evaluator evaluator, Hash hash, int parent, Cost cost) :
action(action),
evaluator(evaluator),
hash(hash),
parent(parent),
cost(cost) {}
};
// ビームサーチの設定
struct Config {
int max_turn;
size_t beam_width;
size_t nodes_capacity;
uint32_t hash_map_capacity;
};
// 削除可能な優先度付きキュー
using MaxSegtree = atcoder::segtree<
pair<Cost,int>,
[](pair<Cost,int> a, pair<Cost,int> b){
if (a.first >= b.first) {
return a;
} else {
return b;
}
},
[]() { return make_pair(numeric_limits<Cost>::min(), -1); }
>;
// ノードの候補から実際に追加するものを選ぶクラス
// ビーム幅の個数だけ、評価がよいものを選ぶ
// ハッシュ値が一致したものについては、評価がよいほうのみを残す
class Selector {
public:
explicit Selector(const Config& config) :
hash_to_index_(config.hash_map_capacity)
{
beam_width = config.beam_width;
candidates_.reserve(beam_width);
full_ = false;
st_original_.resize(beam_width);
}
// 候補を追加する
// ターン数最小化型の問題で、candidateによって実行可能解が得られる場合にのみ finished = true とする
// ビーム幅分の候補をCandidateを追加したときにsegment treeを構築する
bool push(Action action, const Evaluator& evaluator, Hash hash, int parent, bool finished) {
Cost cost = evaluator.evaluate();
if (finished) {
finished_candidates_.emplace_back(Candidate(action, evaluator, hash, parent, cost));
return true;
}
if (full_ && cost >= st_.all_prod().first) {
// 保持しているどの候補よりもコストが小さくないとき
return false;
}
auto [valid, i] = hash_to_index_.get_index(hash);
if (valid) {
int j = hash_to_index_.get(i);
if (hash == candidates_[j].hash) {
// ハッシュ値が等しいものが存在しているとき
if (cost < candidates_[j].cost) {
// 更新する場合
candidates_[j] = Candidate(action, evaluator, hash, parent, cost);
if (full_) {
st_.set(j, {cost, j});
}
return true;
}
return false;
}
}
if (full_) {
// segment treeが構築されている場合
int j = st_.all_prod().second;
hash_to_index_.set(i, hash, j);
candidates_[j] = Candidate(action, evaluator, hash, parent, cost);
st_.set(j, {cost, j});
} else {
// segment treeが構築されていない場合
hash_to_index_.set(i, hash, candidates_.size());
candidates_.emplace_back(Candidate(action, evaluator, hash, parent, cost));
if (candidates_.size() == beam_width) {
// 保持している候補がビーム幅分になったとき
construct_segment_tree();
}
}
return true;
}
// 選んだ候補を返す
const vector<Candidate>& select() const {
return candidates_;
}
// 実行可能解が見つかったか
bool have_finished() const {
return !finished_candidates_.empty();
}
// 実行可能解に到達する「候補」を返す
vector<Candidate> get_finished_candidates() const {
return finished_candidates_;
}
// 最も評価がよい候補を返す
Candidate calc_best_candidate() {
int best = 0;
for (size_t i = 0; i < candidates_.size(); ++i) {
if (candidates_[i].cost < candidates_[best].cost) {
best = i;
}
}
return candidates_[best];
}
void clear() {
candidates_.clear();
hash_to_index_.clear();
full_ = false;
}
private:
size_t beam_width;
vector<Candidate> candidates_;
HashMap<Hash,int> hash_to_index_;
bool full_;
vector<pair<Cost,int>> st_original_;
MaxSegtree st_;
vector<Candidate> finished_candidates_;
void construct_segment_tree() {
full_ = true;
for (size_t i = 0; i < beam_width; ++i) {
st_original_[i] = {candidates_[i].cost, i};
}
st_ = MaxSegtree(st_original_);
}
};
// ターン毎に候補を管理する
class MultiSelectors {
public:
explicit MultiSelectors(const Config& config) :
config_(config)
{
step_max_ = 1;
}
// 候補を追加する
// ターン数最小化型の問題で、candidateによって実行可能解が得られる場合にのみ finished = true とする
// stepは何ターン後に遷移するかを表す（普通のビームサーチなら1）
bool push(Action action, const Evaluator& evaluator, Hash hash, int parent, bool finished, size_t step) {
while (selectors_.size() < step) {
selectors_.emplace_back(Selector(config_));
}
if (selectors_[step - 1].push(action, evaluator, hash, parent, finished)) {
if (step > step_max_) {
step_max_ = step;
}
return true;
}
return false;
}
// expandの直前に呼ぶ
void reset_step_max() {
step_max_ = 1;
}
size_t get_step_max() const {
return step_max_;
}
// 次の候補を保持したSelectorを取り出す
Selector pop_selector() {
Selector ret = move(selectors_.front());
selectors_.pop_front();
return ret;
}
// Selectorを使い回す
void push_selector(Selector&& selector) {
selector.clear();
selectors_.push_back(move(selector));
}
private:
Config config_;
deque<Selector> selectors_;
size_t step_max_;
};
// 深さ優先探索に沿って更新する情報をまとめたクラス
class State {
public:
explicit State(/* const Input& input */) {
// TODO
}
// 次の状態候補を全てselectorに追加する
// 引数
// evaluator : 今の評価器
// hash : 今のハッシュ値
// parent : 今のノードID（次のノードにとって親となる）
void expand(const Evaluator& evaluator, Hash hash, int parent, MultiSelectors& multi_selectors) {
// TODO
}
// actionを実行して次の状態に遷移する
void move_forward(Action action) {
// TODO
}
// actionを実行する前の状態に遷移する
// 今の状態は、親からactionを実行して遷移した状態である
void move_backward(Action action) {
// TODO
}
private:
// TODO
};
// 探索木（二重連鎖木）のノード
struct Node {
Action action;
Evaluator evaluator;
Hash hash;
int parent, child, left, right;
bool active;
int remove_check_turn;
// 根のコンストラクタ
Node(Action action, const Evaluator& evaluator, Hash hash) :
action(action),
evaluator(evaluator),
hash(hash),
parent(-1),
child(-1),
left(-1),
right(-1),
active(true),
remove_check_turn(-1) {}
// 通常のコンストラクタ
Node(const Candidate& candidate, int right) :
action(candidate.action),
evaluator(candidate.evaluator),
hash(candidate.hash),
parent(candidate.parent),
child(-1),
left(-1),
right(right),
active(true),
remove_check_turn(-1) {}
};
// 二重連鎖木に対する操作をまとめたクラス
class Tree {
public:
explicit Tree(const State& state, size_t nodes_capacity, const Node& root) :
state_(state)
{
nodes_.reserve(nodes_capacity);
root_ = nodes_.push(root);
}
// 状態を更新しながら深さ優先探索を行い、次のノードの候補を全てselectorに追加する
void dfs(MultiSelectors& multi_selectors, int turn) {
remove_useless_nodes(turn);
update_root(turn);
int v = root_;
if (!nodes_[v].active) {
// activeなノードがないとき
return;
}
while (true) {
v = move_to_leaf(v);
multi_selectors.reset_step_max();
state_.expand(nodes_[v].evaluator, nodes_[v].hash, v, multi_selectors);
while (remove_nodes_.size() < multi_selectors.get_step_max()) {
remove_nodes_.emplace_back();
}
// 削除可能か確認するターンを設定する
remove_nodes_[multi_selectors.get_step_max() - 1].push_back(v);
nodes_[v].remove_check_turn = turn + multi_selectors.get_step_max();
v = move_to_ancestor(v);
if (v == root_) {
break;
}
}
}
// 根からノードvまでのパスを取得する
vector<Action> get_path(int v) {
// cerr << nodes_.size() << endl;
vector<Action> path;
while (nodes_[v].parent != -1) {
path.push_back(nodes_[v].action);
v = nodes_[v].parent;
}
reverse(path.begin(), path.end());
return path;
}
// 新しいノードを追加する
int add_leaf(const Candidate& candidate) {
int parent = candidate.parent;
int sibling = nodes_[parent].child;
int v = nodes_.push(Node(candidate, sibling));
nodes_[parent].child = v;
if (sibling != -1) {
nodes_[sibling].left = v;
}
// 祖先をactivateする
int u = parent;
while (!nodes_[u].active) {
nodes_[u].active = true;
if (u == root_) {
break;
}
u = nodes_[u].parent;
}
return v;
}
// 最も評価がよいノードを返す
int get_best_leaf(const vector<int>& last_nodes) {
assert(!last_nodes.empty());
int ret = last_nodes[0];
for (int v : last_nodes) {
if (nodes_[v].evaluator.evaluate() < nodes_[ret].evaluator.evaluate()) {
ret = v;
}
}
return ret;
}
private:
State state_;
ObjectPool<Node> nodes_;
int root_;
deque<vector<int>> remove_nodes_;
// 根から一本道の部分は往復しないようにする
void update_root(int turn) {
int child = nodes_[root_].child;
// 後で子供が追加されうるノードはスキップしないようにする
while (child != -1 && nodes_[child].right == -1 && nodes_[root_].remove_check_turn <= turn) {
root_ = child;
state_.move_forward(nodes_[child].action);
child = nodes_[child].child;
}
}
// ノードvの子孫で、最も左にある葉に移動する
int move_to_leaf(int v) {
int child = nodes_[v].child;
while (child != -1) {
// activeなノードが見つかるまで右に移動する
while (!nodes_[child].active) {
child = nodes_[child].right;
}
nodes_[v].active = false;
v = child;
state_.move_forward(nodes_[child].action);
child = nodes_[child].child;
}
nodes_[v].active = false;
return v;
}
// ノードvの先祖で、右への分岐があるところまで移動する
int move_to_ancestor(int v) {
while (v != root_) {
state_.move_backward(nodes_[v].action);
// activeなノードが見つかるまで右に移動する
int u = nodes_[v].right;
while (u != -1) {
if (nodes_[u].active) {
state_.move_forward(nodes_[u].action);
return u;
}
u = nodes_[u].right;
}
v = nodes_[v].parent;
}
return root_;
}
// 不要になったノードを全て削除する
void remove_useless_nodes(int turn) {
if (remove_nodes_.empty()) {
return;
}
for (int v : remove_nodes_.front()) {
if (nodes_[v].child == -1) {
remove_leaf(v, turn);
}
}
remove_nodes_.front().clear();
// 先頭の要素を末尾に移動
remove_nodes_.push_back(move(remove_nodes_.front()));
remove_nodes_.pop_front();
}
// 不要になった葉を再帰的に削除する
void remove_leaf(int v, int turn) {
while (true) {
if (nodes_[v].remove_check_turn > turn) {
// 子供が追加される可能性があるので削除しない
return;
}
int left = nodes_[v].left;
int right = nodes_[v].right;
if (left == -1) {
int parent = nodes_[v].parent;
if (parent == -1) {
cerr << "ERROR: root is removed" << endl;
exit(-1);
}
nodes_.pop(v);
nodes_[parent].child = right;
if (right != -1) {
nodes_[right].left = -1;
return;
}
v = parent;
} else {
nodes_.pop(v);
nodes_[left].right = right;
if (right != -1) {
nodes_[right].left = left;
}
return;
}
}
}
};
// ビームサーチを行う関数
vector<Action> beam_search(const Config& config, State state, Node root) {
Tree tree(state, config.nodes_capacity, root);
// 新しいノード候補の集合
MultiSelectors multi_selectors(config);
for (int turn = 0; turn < config.max_turn; ++turn) {
// Euler Tour で selector に候補を追加する
tree.dfs(multi_selectors, turn);
Selector selector = multi_selectors.pop_selector();
if (selector.have_finished()) {
// ターン数最小化型の問題で実行可能解が見つかったとき
Candidate candidate = selector.get_finished_candidates()[0];
vector<Action> ret = tree.get_path(candidate.parent);
ret.push_back(candidate.action);
return ret;
}
if (turn == config.max_turn - 1) {
// 最終ターン
Candidate candidate = selector.calc_best_candidate();
vector<Action> ret = tree.get_path(candidate.parent);
ret.push_back(candidate.action);
return ret;
}
// 新しいノードを追加する
for (const Candidate& candidate : selector.select()) {
tree.add_leaf(candidate);
}
// Selector を使い回す
multi_selectors.push_selector(move(selector));
}
unreachable();
}
} // namespace beam_search
スコアは改善していないですが、一応使用例を貼っておきます。
Euler Tour の辺を保持する実装 (2024/02/07 追記)
二重連鎖木や Object Pool を使うのではなく、探索木の有向辺を Euler Tour の順序で保持するほうが速いという話があり、実装してみました。
状態を差分計算するときは Euler Tour の配列に前から順にアクセスします。
探索木を更新するときは、辺の追加や削除をしながら、別の配列に Euler Tour をコピーしています。
ビームサーチライブラリ
#include <bits/stdc++.h>
using namespace std;
namespace beam_search {
// ビームサーチの設定
struct Config {
int max_turn;
size_t beam_width;
size_t tour_capacity;
uint32_t hash_map_capacity;
};
// 連想配列
// Keyにハッシュ関数を適用しない
// open addressing with linear probing
// unordered_mapよりも速い
// nは格納する要素数よりも16倍ほど大きくする
template <class Key, class T>
struct HashMap {
public:
explicit HashMap(uint32_t n) {
if (n % 2 == 0) {
++n;
}
n_ = n;
valid_.resize(n_, false);
data_.resize(n_);
}
// 戻り値
// - 存在するならtrue、存在しないならfalse
// - index
pair<bool,int> get_index(Key key) const {
Key i = key % n_;
while (valid_[i]) {
if (data_[i].first == key) {
return {true, i};
}
if (++i == n_) {
i = 0;
}
}
return {false, i};
}
// 指定したindexにkeyとvalueを格納する
void set(int i, Key key, T value) {
valid_[i] = true;
data_[i] = {key, value};
}
// 指定したindexのvalueを返す
T get(int i) const {
assert(valid_[i]);
return data_[i].second;
}
void clear() {
fill(valid_.begin(), valid_.end(), false);
}
private:
uint32_t n_;
vector<bool> valid_;
vector<pair<Key,T>> data_;
};
using Hash = uint32_t; // TODO
// 状態遷移を行うために必要な情報
// メモリ使用量をできるだけ小さくしてください
struct Action {
// TODO
Action() {
// TODO
}
bool operator==(const Action& other) const {
// TODO
}
};
using Cost = int;
// 状態のコストを評価するための構造体
// メモリ使用量をできるだけ小さくしてください
struct Evaluator {
// TODO
Evaluator() {
// TODO
}
// 低いほどよい
Cost evaluate() const {
// TODO
}
};
// 展開するノードの候補を表す構造体
struct Candidate {
Action action;
Evaluator evaluator;
Hash hash;
int parent;
Candidate(Action action, Evaluator evaluator, Hash hash, int parent) :
action(action),
evaluator(evaluator),
hash(hash),
parent(parent) {}
};
// ノードの候補から実際に追加するものを選ぶクラス
// ビーム幅の個数だけ、評価がよいものを選ぶ
// ハッシュ値が一致したものについては、評価がよいほうのみを残す
class Selector {
public:
explicit Selector(const Config& config) :
hash_to_index_(config.hash_map_capacity)
{
beam_width = config.beam_width;
candidates_.reserve(beam_width);
full_ = false;
costs_.resize(beam_width);
for (size_t i = 0; i < beam_width; ++i) {
costs_[i] = {0, i};
}
}
// 候補を追加する
// ターン数最小化型の問題で、candidateによって実行可能解が得られる場合にのみ finished = true とする
// ビーム幅分の候補をCandidateを追加したときにsegment treeを構築する
void push(const Candidate& candidate, bool finished) {
if (finished) {
finished_candidates_.emplace_back(candidate);
return;
}
Cost cost = candidate.evaluator.evaluate();
if (full_ && cost >= st_.all_prod().first) {
// 保持しているどの候補よりもコストが小さくないとき
return;
}
auto [valid, i] = hash_to_index_.get_index(candidate.hash);
if (valid) {
int j = hash_to_index_.get(i);
if (candidate.hash == candidates_[j].hash) {
// ハッシュ値が等しいものが存在しているとき
if (full_) {
// segment treeが構築されている場合
if (cost < st_.get(j).first) {
candidates_[j] = candidate;
st_.set(j, {cost, j});
}
} else {
// segment treeが構築されていない場合
if (cost < costs_[j].first) {
candidates_[j] = candidate;
costs_[j].first = cost;
}
}
return;
}
}
if (full_) {
// segment treeが構築されている場合
int j = st_.all_prod().second;
hash_to_index_.set(i, candidate.hash, j);
candidates_[j] = candidate;
st_.set(j, {cost, j});
} else {
// segment treeが構築されていない場合
int j = candidates_.size();
hash_to_index_.set(i, candidate.hash, j);
candidates_.emplace_back(candidate);
costs_[j].first = cost;
if (candidates_.size() == beam_width) {
// 保持している候補がビーム幅分になったときにsegment treeを構築する
full_ = true;
st_ = MaxSegtree(costs_);
}
}
}
// 選んだ候補を返す
const vector<Candidate>& select() const {
return candidates_;
}
// 実行可能解が見つかったか
bool have_finished() const {
return !finished_candidates_.empty();
}
// 実行可能解に到達するCandidateを返す
vector<Candidate> get_finished_candidates() const {
return finished_candidates_;
}
// 最もよいCandidateを返す
Candidate calculate_best_candidate() const {
if (full_) {
size_t best = 0;
for (size_t i = 0; i < beam_width; ++i) {
if (st_.get(i).first < st_.get(best).first) {
best = i;
}
}
return candidates_[best];
} else {
size_t best = 0;
for (size_t i = 0; i < candidates_.size(); ++i) {
if (costs_[i].first < costs_[best].first) {
best = i;
}
}
return candidates_[best];
}
}
void clear() {
candidates_.clear();
hash_to_index_.clear();
full_ = false;
}
private:
// 削除可能な優先度付きキュー
using MaxSegtree = atcoder::segtree<
pair<Cost,int>,
[](pair<Cost,int> a, pair<Cost,int> b){
if (a.first >= b.first) {
return a;
} else {
return b;
}
},
[]() { return make_pair(numeric_limits<Cost>::min(), -1); }
>;
size_t beam_width;
vector<Candidate> candidates_;
HashMap<Hash,int> hash_to_index_;
bool full_;
vector<pair<Cost,int>> costs_;
MaxSegtree st_;
vector<Candidate> finished_candidates_;
};
// 深さ優先探索に沿って更新する情報をまとめたクラス
class State {
public:
explicit State() {
// TODO
}
// EvaluatorとHashの初期値を返す
pair<Evaluator,Hash> make_initial_node() {
// TODO
}
// 次の状態候補を全てselectorに追加する
// 引数
// evaluator : 今の評価器
// hash : 今のハッシュ値
// parent : 今のノードID（次のノードにとって親となる）
void expand(const Evaluator& evaluator, Hash hash, int parent, Selector& selector) {
// TODO
}
// actionを実行して次の状態に遷移する
void move_forward(Action action) {
// TODO
}
// actionを実行する前の状態に遷移する
// 今の状態は、親からactionを実行して遷移した状態である
void move_backward(Action action) {
// TODO
}
private:
// TODO
};
// Euler Tourを管理するためのクラス
class Tree {
public:
explicit Tree(const State& state, const Config& config) :
state_(state)
{
curr_tour_.reserve(config.tour_capacity);
next_tour_.reserve(config.tour_capacity);
leaves_.reserve(config.beam_width);
buckets_.assign(config.beam_width, {});
}
// 状態を更新しながら深さ優先探索を行い、次のノードの候補を全てselectorに追加する
void dfs(Selector& selector) {
if (curr_tour_.empty()) {
// 最初のターン
auto [evaluator, hash] = state_.make_initial_node();
state_.expand(evaluator, hash, 0, selector);
return;
}
for (auto [leaf_index, action] : curr_tour_) {
if (leaf_index >= 0) {
// 葉
state_.move_forward(action);
auto& [evaluator, hash] = leaves_[leaf_index];
state_.expand(evaluator, hash, leaf_index, selector);
state_.move_backward(action);
} else if (leaf_index == -1) {
// 前進辺
state_.move_forward(action);
} else {
// 後退辺
state_.move_backward(action);
}
}
}
// 木を更新する
void update(const vector<Candidate>& candidates) {
leaves_.clear();
if (curr_tour_.empty()) {
// 最初のターン
for (const Candidate& candidate : candidates) {
curr_tour_.push_back({(int)leaves_.size(), candidate.action});
leaves_.push_back({candidate.evaluator, candidate.hash});
}
return;
}
for (const Candidate& candidate : candidates) {
buckets_[candidate.parent].push_back({candidate.action, candidate.evaluator, candidate.hash});
}
auto it = curr_tour_.begin();
// 一本道を反復しないようにする
while (it->first == -1 && it->second == curr_tour_.back().second) {
Action action = (it++)->second;
state_.move_forward(action);
direct_road_.push_back(action);
curr_tour_.pop_back();
}
// 葉の追加や不要な辺の削除をする
while (it != curr_tour_.end()) {
auto [leaf_index, action] = *(it++);
if (leaf_index >= 0) {
// 葉
if (buckets_[leaf_index].empty()) {
continue;
}
next_tour_.push_back({-1, action});
for (auto [new_action, evaluator, hash] : buckets_[leaf_index]) {
int new_leaf_index = leaves_.size();
next_tour_.push_back({new_leaf_index, new_action});
leaves_.push_back({evaluator, hash});
}
buckets_[leaf_index].clear();
next_tour_.push_back({-2, action});
} else if (leaf_index == -1) {
// 前進辺
next_tour_.push_back({-1, action});
} else {
// 後退辺
auto [old_leaf_index, old_action] = next_tour_.back();
if (old_leaf_index == -1) {
next_tour_.pop_back();
} else {
next_tour_.push_back({-2, action});
}
}
}
swap(curr_tour_, next_tour_);
next_tour_.clear();
}
// 根からのパスを取得する
vector<Action> calculate_path(int parent, int turn) const {
// cerr << curr_tour_.size() << endl;
vector<Action> ret = direct_road_;
ret.reserve(turn);
for (auto [leaf_index, action] : curr_tour_) {
if (leaf_index >= 0) {
if (leaf_index == parent) {
ret.push_back(action);
return ret;
}
} else if (leaf_index == -1) {
ret.push_back(action);
} else {
ret.pop_back();
}
}
unreachable();
}
private:
State state_;
vector<pair<int,Action>> curr_tour_;
vector<pair<int,Action>> next_tour_;
vector<pair<Evaluator,Hash>> leaves_;
vector<vector<tuple<Action,Evaluator,Hash>>> buckets_;
vector<Action> direct_road_;
};
// ビームサーチを行う関数
vector<Action> beam_search(const Config& config, const State& state) {
Tree tree(state, config);
// 新しいノード候補の集合
Selector selector(config);
for (int turn = 0; turn < config.max_turn; ++turn) {
// Euler Tourでselectorに候補を追加する
tree.dfs(selector);
if (selector.have_finished()) {
// ターン数最小化型の問題で実行可能解が見つかったとき
Candidate candidate = selector.get_finished_candidates()[0];
vector<Action> ret = tree.calculate_path(candidate.parent, turn + 1);
ret.push_back(candidate.action);
return ret;
}
assert(!selector.select().empty());
if (turn == config.max_turn - 1) {
// ターン数固定型の問題で全ターンが終了したとき
Candidate best_candidate = selector.calculate_best_candidate();
vector<Action> ret = tree.calculate_path(best_candidate.parent, turn + 1);
ret.push_back(best_candidate.action);
return ret;
}
// 木を更新する
tree.update(selector.select());
selector.clear();
}
unreachable();
}
} // namespace beam_search
元の二重連鎖木を使った実装よりも高速に動作するようでした。
最後に
この記事では私の実装のみを紹介しました。実装した人によって異なる部分があるので調べてみると面白いかもしれません。
最後まで読んでくださりありがとうございました。
- 若干意味合いが異なるかもしれません。私は同一視しています。↩
- 履歴は永続配列を使用することで高速化できます。↩
-
状態遷移の計算量が状態をコピーする計算量と同程度の場合には愚直なビームサーチでよいと思います。過去のAHCを見る限りだと、状態遷移が
の場合が多く、愚直なビームサーチが強い問題は少ないです。↩
-
内部的には
std::vector::reserveを呼び出しているだけです。↩ - open addressing と linear probing について参考にさせていただきました。私の実装では、遅延削除は行っていません。↩
-
浮動小数点演算を含む関数に対して、正確な逆関数を定義することは難しいことが多いです。例えば
として
から
を復元するときに、
を行うと思いますが、数値誤差を考慮すると
- 正確には ObjectPool におけるインデックスです。↩
- 16ビットにエンコードすることも可能ですが、個人的には限界まで減らすモチベーションがなかったです。↩
AHC029 参加記
RECRUIT 日本橋ハーフマラソン 2024冬（AtCoder Heuristic Contest 029） お疲れ様でした。
システムテストの得点率が 73.4 % で優勝しました！
問題概要
所持金の最大化を行うカードゲームです。
詳細は公式の問題文を参考にしてください。
方針
貪欲法をモンテカルロ法によって改善しました。
コンテスト上位者の中で私の貪欲法の性能は悪かったため、貪欲法についてはあまり参考にならないかもしれません。
貪欲法の準備
貪欲法を説明する前に、貪欲法の説明に必要なものについて述べます。
各ターンの流れ
各ターンの流れについて、問題文では
- 手札からカードを1枚選んで使う。
- 必要に応じてプロジェクが補充される。
- 補充するカードの候補が提示される。
- 補充するカードを1枚選んで補充する。
の順序で書かれています。
しかし、補充するカードを選んだ直後に手札のカードを使うため、これら2つの処理をまとめて考えることができます。
そこで、各ターンの流れを
- 補充するカードの候補が提示される。
- 補充するカードを1枚選んで補充し、手札からカードを1枚選んで使う。
- 必要に応じてプロジェクトが補充される。
のように捉えることにします。
1ターン目は、手札の1枚目のカードが欠けていて、補充するカードの候補として手札の1枚目のカードだけが提示されていることにすると実装しやすいです。また、問題文における最終ターンでは、補充するカードとして1枚目を選ぶのが最適なので、最後のカードを選ぶ手順は省略して考えることができます。
行動の定義
（補充するカード、使用するカードと対象のプロジェクト）という組を行動と定義します。問題文中の記号で表すと です。
貪欲法
行動価値関数を定義し、行動価値が最大になる行動を選択します。
直後に取りうる行動を全探索するため、ループ内部の時間計算量を とすると、全体の計算量は
になります。
行動価値関数は、後述する状態価値関数とは直接的な関係がないことに注意してください。（強化学習における行動価値関数とは定義が異なると思います。）
行動価値関数は、以下の項目の変化を評価します。
- 所持金
- プロジェクト
- 手札のカード
- 使用した投資カードの枚数
単位は所持金に合わせることにします。すなわち、所持金以外の項目の変化をお金に換算します。所持金ベースで考えることで、定量的な考察がしやすくなると思います。
所持金の変化
カードのコスト と、カードを使うことによって完了するプロジェクトの価値
を用いて計算することができます。
プロジェクトの残務量を減らす価値
プロジェクトが完了する場合 、所持金の変化のみを考えることにします。すなわち、新しいプロジェクトの評価は0とします。
プロジェクトが完了しない場合 、
をベースにして評価しました。減らした残務量の割合に応じてプロジェクトの価値を得られるイメージです。
実際にはお金を稼げないことや、累積和が よりも大きくなることを踏まえると、
を下方修正したほうがよさそうです。
また、カード生成方法を見ると、所持金と減らせる残務量が大体同じであることがわかります。そのため、所持金よりも大きい残務量のプロジェクトを完了するには時間がかかることが予想されます。
そこで、 が大きいほど、
を下方修正するようにしました。
プロジェクトを取りやめる価値
をベースとしました。
ただし、 が小さいほど、取りやめる価値を下げました。以下の2つの理由があります。
が小さければすぐに完了できるから。
- 途中まで進めたプロジェクトを業務転換カードで換えないようにするため。
カードの変化
カードの価値を計算すればよいです。
ゲーム終了時に強いカードを残すのは無駄なので、最後20ターンほどからカードの価値を下げるようにしています。
通常労働カードの価値
としました。プロジェクトの残務量を0未満にする（overkill と呼ぶことにします）ことがあるため、
よりも少しだけ小さくしています。
全力労働カードの価値
としました。overkill が通常労働よりも多発することや、効率の悪いプロジェクトも進めるため、
よりも小さくしています。
キャンセルカードと業務転換カードの価値
としました。
増資について
900ターンぐらいまで、増資カードは買えるときに必ず買うようにし、買ったらすぐに使用しました。評価値が無限大と考えることもできます。
実は増資カードが高いときには買わないべきで、増資ストックも検討するべきなのですが、コンテスト中にうまく実装できませんでした。今後暇なときに実装したいです。
増資カードを購入し、しばらく使わずに手札に残すことを増資ストックと呼ぶことにします。 増資すると多くの変数が2倍になりますが、所持金は2倍になりません。相対的に所持金が半分になると考えることができます。 所持金が少ないと効率が悪くなりやすいため、所持金が少ない期間は短いほうがよいです。 以上より、所持金が少ないときに増資カードを一気に使いたいというモチベーションが生まれます。 一方で、所持金を一時的でも極端に減らしていいのか、増資ストック中に使える手札のカード枚数を減らしていいのかなど、増資ストックのデメリットに関しては自分の中で結論が出てないです。 （2024/01/08 追記: 増資ストックではなく増資スタックでした。同じアイテムを複数まとめることをスタックというらしいです。）増資ストックとは
貪欲法の高速化
モンテカルロ法におけるボトルネックは貪欲法になります。そのため、貪欲法は速ければ速いほどよいです。
処理をまとめる
基本的に可能な行動全てについて行動価値を計算しますが、補充したカードの価値など、いくつかの行動で共通する要素があるので、それらはなるべく1回しか計算しないようにしました。
行動価値関数として実装しないため、実装が汚くなりました。
貪欲法における枝刈り
行動を全探索すると書きましたが、明らかに無駄な探索は枝刈りします。
枝刈りするのは以下のケースです。
- 下位互換のカードは補充しない。すなわち、補充するカード候補の中に、同じ種類でコスト
が小さくて効能がそれ以上のものがあるようなカードを補充することを禁止する。
- 手札に存在する効能が同じカードを同一視する。
枝刈りは貪欲法を高速化するだけでなく、モンテカルロ法でプレイアウトを行う行動の候補を選ぶときにも効果的に働きます。
状態評価関数
モンテカルロ法で状態を評価するときに使います。
状態評価値の比較は同じターン同士でしか行われないため、同じターンで評価基準が揃っていればよいものとします。
ゲームが終了している場合、所持金を評価しました。以降、ゲームが終了していない場合について説明します。
基本的には行動価値関数と同じ考え方で状態を評価します。
所持金の評価
所持金が より大きいときは増資カードを買える可能性があるので所持金に応じて高く評価しました。
カードの評価
行動価値関数とほぼ一緒です。
プロジェクトの評価
をベースとしました。
ただし、キャンセルカードや業務転換カードを保持しているときは の場合のペナルティを緩和しました。
増資カード使用回数の評価
としました。
その他
最後に評価値を 0.95 乗して極端な値の影響を下げました。スコアはあまり変わらなかったと思います。
モンテカルロ法
可能な行動について、貪欲法によるプレイアウトで評価値の期待値を求め、期待値が最大の行動を選択します。
処理の流れは以下の通りです。
- 行動価値が大きい7つの行動を候補として取得する。
- およそ2msec経過するまで以下の処理を繰り返す。
- プレイアウトに必要な情報を生成する。
- 各候補について9ターンの貪欲プレイアウトを行う。
- 一定時間が経過した場合、期待値が最小の候補を削除する。
- 期待値が最大の候補を選択する。
重要そうな部分について説明します。
シミュレーション回数
より正確な期待値を求めるにはプレイアウトの回数を増やすことが大切です。
そこで、プレイアウト回数を増やす工夫をいくつか行いました。
プレイアウトを行う行動の選択
可能な行動が 通りあり、理想的にはそれら全てについて十分な回数のプレイアウトを行いたいです。
しかし、実際には実行時間が非常に限られているため、可能な行動全てに対してプレイアウトを行うとプレイアウト回数が減ってしまいます。
そこで、行動価値が大きい行動についてのみプレイアウトを行うようにしました。行動価値関数は完璧ではありませんが、上位何個か選べば最良のものが含まれているだろうという感覚です。
また、プレイアウトを行う行動候補を選ぶときに、補充するカードと使用するカードの組が等しいものは1個までとしました。すなわち、候補の中にカードを適用するプロジェクトだけが異なるような行動がないようにしました。おそらく候補の多様性が重要なのだと思います。
データ生成
基本的には問題文と同様の方法でプロジェクトとカードについてのデータを生成します。
プロジェクトの生成には隠しパラメータが使用されていないため、一度生成したプロジェクトを別のターンでも使いまわしました。
一方で、カードの生成には隠しパラメータが使用されているため、今までに出現したカードの統計データをとり、毎ターンその分布にしたがってカードを生成しました。
統計データは各カードの種類について出現した回数を数えるだけですが、出現した回数の初期値は 21, 11, 11, 6, 4 としました。序盤はデータが少ないので無情報の期待値を参考にできるということです。スコアへの寄与は小さいと思います。
プレイアウト
狭義的には最終状態までシミュレーションすることをプレイアウトと言いますが、途中までシミュレーションすることもプレイアウトと呼ぶことにします。
最終状態までシミュレーションした場合、序盤や中盤では非常に長い時間がかかるうえ、ランダム性が大きくなって評価したい行動のよさが計りにくくなります。そこで、シミュレーションを9ターン後まで行い、9ターン後の状態価値を最大化することを目標としました。ただし、9ターン後までにゲームが終了する場合はゲームを終了させ、そのときの評価は所持金になります。
候補の削除
7つの行動候補を段階的に減らし、最終的に2つの候補についてプレイアウトを行うようにしました。一定時間が経過するごとに最も期待値が低い候補を1つ削除します。
最良でなさそうな候補を途中で削除することにより、よさそうな候補のシミュレーション回数を増やすことができます。
感想
全体
モンテカルロ法について、AHCの過去問で何度も練習していて、今回結果を出せて非常に嬉しく思います。
問題文を読み終えた時点で、数ターンのプレイアウトを行うモンテカルロ法が強いことを確信し、自分の得意ジャンルであることから優勝するチャンスだと思いました。
大まかな方針はすぐに決まりましたが、「数ターンのプレイアウト」の「数ターン」は20ターンぐらいだと予想していましたし、モンテカルロ法によって、スコアはせいぜい2倍ぐらいにしかならないだろうと思っていたので、その辺りの感覚はずれていました。（スコアは平均して32倍程度、順位表のスコアは12倍程度でした。）
優勝できたのはよかったですが、感想戦において貪欲法で大きく負けていたことがわかり、貪欲法については反省点が多そうです。貪欲法（と状態評価関数）が基本で、モンテカルロ法を行う場合でも、最終的なスコアは貪欲法に大きく依存すると思っていたのでかなり意外な結果でした。
相対評価システム
自分がコンテストを荒らしている感じがして楽しかったです。
今回の相対評価が気に入らない人がいるみたいなので個人的な意見を少し書きます。
順位スコアはテストケースにおける圧倒的1位の人にとって嬉しくないので基本的にはよくないと思っています。元のスコアが指数的な分布なら、対数をとって絶対スコアか相対スコアだとバランスがいいかもしれないです。指数的な分布でないとき（一様分布など）に対数をとると、僅差になって非ACのペナルティが大きくなるのでよくなさそうです。
バグ
自分の実装において、直すとスコアが悪化するバグが2つあったので紹介します。
状態評価関数におけるカードの評価
私の状態評価関数は、カードを使ってプロジェクトが補充された直後の状態に適用するようになっています。カードが1枚欠けた状態です。
バグで欠けた1枚のカード、すなわち最後に使用したカードも評価に入れていることに気づき、バグを直したところスコアが少し下がりました。
原因に心当たりはありません。単なるスコアのぶれかもしれないです。
カードの種類の統計の取り方
カードの種類の統計をとるとき、提示されたカードの1枚目は必ず通常労働カードで統計データに含めてはいけないのですが、バグで統計に含まれていました。バグを直したところ、スコアが少し下がりました。
が小さいときに大きく影響しそうですが、通常労働カードは基本的なカードなので、少し出やすくしたほうがいいのでしょうか。納得はしていませんが、致命的なバグではないような気がしています。
最後に
コンテストを主催してくださったリクルート様、学生賞金があって大変感謝しております。長期AHCとしては珍しいタイプの問題で、とても面白かったです。
AtCoderの皆様、参加者の皆様、最後まで読んでくださった読者の皆様、ありがとうございました。
ALGO ARTIS プログラミングコンテスト2023 冬（AtCoder Heuristic Contest 028） も参加する予定なので対戦よろしくお願いします。
AHC027 参加記
HACK TO THE FUTURE 2024 （AtCoder Heuristic Contest 027） お疲れ様でした。
システムテストの得点率が 98.44 % で準優勝しました！
前提知識
- 木上のビームサーチ
- バックトラック法
- low-link
問題概要
マスの二次元グリッドがあり、各マス
の汚れは毎ターン
だけ増加します。ロボット掃除機が訪れたマスは汚れが 0 になります。ロボットの掃除ルートで、平均汚れができるだけ小さいものを求めてください。
詳細は公式の問題文を参考にしてください。
制約
方針
ビームサーチで初期解を生成し、焼きなまし法を適用しました。
私が調べた限りでは、bowwowforeach さんや siman さんと大まかな方針は一緒で、他の上位の人は方針が違いそうでした。
ビームサーチ
深さ
深さ1を に設定しました。
ただし、深さが偶数になるように必要に応じてインクリメントしました。深さを偶数にすることで始点と終点が同じものを得ることができます。
多様性の確保
(現在のマス, 直前にいたマス) をハッシュ値の代わりに使用しました。ビーム幅は 弱ということになります。
次のノードを選ぶときにソートする必要がなくなるため、ビームサーチが高速になります。
評価関数
4つの項目を評価しました。
影響が大きそうなものから順に説明します。
評価項目1: 汚れの累積和
後述する汚れの総和を毎ターン加算したものです。
評価項目2: 汚れの総和
現在の各マスの汚れの総和です。初期状態における汚れは で初期化しました。
汚れの総和は、直前の汚れの総和と、そのターンに回収した汚れの量から計算することができます。回収した汚れの量を求めるには、各マスについて最後に訪れたターンをメモしておけばよいです。
の代わりに
を使うとスコアがよくなりました。
をそのまま使うと
が大きいマスばかり訪問するのに対し、修正したものを使うと比較的バランスよく訪問していました。
評価項目3: 移動平均
指数移動平均を計算し、絶対値が大きいほどよい評価をしました。
ここでいう指数移動平均とは、 ターン目の座標の変化量を
として
を表すものとします。平滑化係数は としました。
更新は
のように行います。
が大きい狭い領域に居続けるのがよくないので、移動平均の絶対値が大きいものを優先するとよいと考えられます。
例えば Seed = 25 は、 が大きい領域が非常に限られており、評価関数によっては中心やや下の領域から出なくなります。
評価項目4: 周期解としての生スコア
汚れの累積和では、掃除ルートが繰り返されることが考慮されていません。そこで、 ターン目まで何もしなかった場合の生スコアを評価しました。
各マスを訪問するたびに、訪問時刻を挿入するような感じで差分計算を行うことができます。
未訪問のマスは ターンの周期で訪問されるものとして計算しました。
差分更新
木上のビームサーチでは、評価の更新に必要な情報を深さ優先探索に沿って更新します。私の実装だと、評価値はノードに保存されるので深さ優先探索では更新する必要がありません。
評価の更新に必要な情報で動的なものは、現在のターンと各マスを最後に訪れた時刻の2つです。これらは で更新できます。
簡易的な C++ のコードを載せておきます。座標が一次元に平坦化されていることに注意してください。
struct State { // 経過したターン数 int turn; // timestamps[i] = マスiを訪問した時刻の配列 vector<vector<int>> timestamps; // 深さ優先探索における子孫への移動 void move_forward(int position) { timestamps[position].push_back(++turn) } // 深さ優先探索における親への移動 void move_backward(int position) { assert(timestamps[position].back() == turn) timestamps[position].pop_back(); --turn; } }
前回訪問した時刻をノードに保存することで、一次元配列のみで実装することも可能です。2
未訪問頂点の対応
未訪問頂点が減るように評価関数を工夫しているつもりですが、それでも未訪問の頂点が残る場合があります。
なるべく低コストで未訪問頂点を挿入したいので、2マス拡張（後述）を基本として、2マス拡張ができないときは1マス拡張（後述）を行いました。
元々
だった経路を
のように拡張するのが1マス拡張で
のように拡張するのが2マス拡張です。
ビームサーチの結果
コンテスト後にビームサーチのみのコードを提出したところ、44Gほどのスコアが出ていました。
焼きなまし法
3種類の近傍で焼きなまします。
近傍1: L字型の変更
a, b, c を a, d, c のように変更します。L字型の部分を反対側に曲げる感じです。
経路長が変わらないので高速に計算できます。
近傍2: 十字型の解消
掃除ルートが同じマスでクロスしていた場合に、区間の訪問順序を反転させてクロスを解消します。2-optをイメージするとよいかもしれません。
クロスを解消した後にL字型の変更が行われると1マス分得をすることができます。
例えば
の区間を反転させて、クロスしていたところにL字型の変更を施すと
になります。
反転区間の長さが約150以上の場合にはクロスの解消を行いませんでした。区間が長すぎると訪問時間の変更によるデメリットのほうが大きくなります。
近傍3: バックトラック法による破壊再構築
掃除ルートの連続する 4 ~ 8 マスを破壊し、バックトラック法で生成した経路で置き換えます。
バックトラック法の工夫1: 同じマスを通らない
よくあるバックトラック法と同様に、同じマスを複数回通らないようにしました。
ただし、関節点については2回まで同じマスを通ってもよいものとしました。関節点も1回しか通れないようにした場合、関節点の削除はできても挿入はできないので、操作が不可逆的になります。
関節点については、low-link で事前に列挙しました。
バックトラック法の工夫2: 経路長の制限
経路長を無制限にすると探索空間が膨大すぎて計算が終わりません。そこで、経路長は「元の経路長 + 2」までとしました。
バックトラック法の工夫3: 経路長による枝刈り
経路長の最大値を設定したため、経路長による枝刈りができます。
具体的には、現在のパスの長さと現在位置から終点までの距離を足したものが、経路長の最大値よりも大きいならば枝刈りできます。
全頂点から幅優先探索をするのは計算量的に少し重そうだったので、距離はマンハッタン距離で代用しました。
バックトラック法の工夫4: よさそうな頂点から探索する
次の頂点を選ぶときに、訪問したときに得られるスコアの改善量が大きい頂点を優先しました。
具体的にはスコアの改善量を として、
が大きいものから順に探索しました。
工夫できなかったこと
経路長が変わる場合、変更部分以降の時刻が全てずれ、更新の計算量が大きくなりました。
ただ更新するだけだともったいなかったので、経路全体を逆順にしたり、掃除ルートのスタート地点を変更したりしました。ときどきこの処理が行われることを前提とし、各近傍では掃除ルートの終わりと始まりをまたぐような遷移は考えませんでした。
各近傍の選択比率
近傍1, 2, 3 の選択比率は 6 : 3 : 1 としました。選択比率と実行時間は異なることに注意してください。
途中で掃除ルートを2倍にする
焼きなましの実行時間の8割が経過したときに、同じ掃除ルートを2回繰り返したものに変更しました。
掃除ルートが長くなって焼きなましにくくなるデメリットがある一方で、全頂点を2回訪れていると制約が緩くなって焼きなましやすくなったり、汚れやすさが極端に低い頂点の掃除頻度を半減できたりするメリットがありそうです。
試行回数
Seed = 0 () で約 2,500,000 回のループを回せていました。
日記
暇な人向けです。
12/1
- 問題を読む。
- 適当に初期解を作り、破壊再構築の焼きなましをしたい。
- 経路長が変わると差分計算が面倒。
- 遺伝的アルゴリズムとか使えないかな。無理か。
- 移動が自由なら Introduction to Heuristic Contest に似ているな。新ジャッジコン で使った貪欲とか活かせるといいな。
- 掃除ルートの周期性を考えなくても、経路が十分長ければ接続部分は無視できそう。
12/2
- BFSに沿ったDPをベースとした貪欲とかがいいのかな。
- 掃除できた汚れの量を評価関数にしてビームサーチを打つほうがいいか。
- 今のマスをハッシュ値として多様性を確保できる。
- Sliding Tree Puzzle と一緒だ。
- 差分更新ビームサーチライブラリの出番が来た！
- 今のマスをハッシュ値として多様性を確保できる。
- ビームサーチを実装する。
- 初提出で38G。悪くない。
- 最初と最後の接続部分を何とかしたい
12/3
- 汚れの累積値という生スコアを評価していないことに気づいて評価する。
- その他にも最適化をして41G。
- 局所探索で改善できそうなことを確認する。
12/4
- 未訪問箇所の挿入を改善する。
- 44G。
12/5
- 平日で唯一の休み。
- 焼きなましの遷移をいくつか考える。
- L字型変更と十字型解消を実装して45G。
- バックトラック法といくつか最適化をして49Gで暫定2位。よし！
12/6, 12/7
- ビームサーチで子孫の数を制限して悪化したのでやめる。
- 過去改変機能をビームサーチに組み込んでも悪化するのでやめる。
- 過去改変した情報を保持しておけば木上のビームサーチでも過去改変できる。
- アイディアとしては面白いと思うが、うまくいかないので諦める。
12/8
- プレテストの
の分布を雑に確認したところ、
- システムテストで下がるパターンだ。
- 移動平均をビームサーチの評価関数に組み込んでスコアが上がる。
12/9
- 朝から微熱が出ていて嫌な気分。
- バックトラック法で関節点を2回通れるようにしたらスコアが上がった。
- 午後に40度の熱が出て寝込む。発熱してもコンテストに参加するつもりだったが、40度までいくと自分には無理だった。
- 調子いいときに発熱するのはやめてほしい。
- まだパラメータ調整してないので明日はコンテストに取り組みたい。
12/10
- 38度弱だったので問題なくコンテストに参加できた。
- リファクタリングとパラメータ調整をやって暫定トップに立つ。
- おそらく潜伏していた cuthbert さんに抜かれ、2位でフィニッシュ。
- 最後やりきれない感じだったが悪くはなさそう。
- 解説放送を見ずに寝る。
12/11
- システムテストが終わって、ギリギリ bowwowforeach さんに勝っていて嬉しかった。
- エラーが起きなかったのもよかった。
- 解説放送を見て、
に比例した割合で訪問するなどの話を聞いて勉強になった。
- 検査によりインフルエンザに罹患していたことがわかる。
感想
木上のビームサーチも、バックトラック法を遷移とする焼きなましも、過去のAHCで何度か実装していたので、非常にスムーズに実装することができました。過去問の復習を活かせたという実感があって非常に嬉しいです。
一方で、優勝するにはもう少し adhoc な考察や実装をする必要があった気がしていて反省もしています。
最後に
コンテストを主催してくださったフューチャー様、writer の wata さん、コンテスト参加者の皆様、そして最後まで読んでくださった読者の皆様、ありがとうございました。
