ymatsux@競技プログラミング

2022年11月26日土曜日

HTTF 2023 予選 (AHC 016) 3位解法の解説

2022年11月に開催された HACK TO THE FUTURE 2023 予選 (AtCoder Heuristic Contest 016) の3位解法の解説をします。

$ \epsilon = 0 $ の場合

この場合は厳密に解けるので、特別扱いします。頂点数毎の非同型なグラフの個数は OEIS で調べられます（リンク）。これによると、頂点数4では11個、頂点数5では34個、頂点数6では156個のグラフが存在します。よって、必要なグラフの個数 M に合わせて、頂点数 N を選んで、非同型なグラフを M 個構築すれば良いでしょう。

非同型なグラフを列挙したり、クエリに答える際に、グラフの同型判定が必要になりますが、同型なグラフの「標準系」を定義するのが簡便でしょう。グラフの文字列表現は問題文で定義されていますが、文字列表現が辞書順最小になるような頂点の番号付けを標準系とすれば良いです。

解法1: 次数列で分類する

まず、グラフの次数列（各頂点の次数を降順に並べたもの）に注目する解法を考えましょう。これは最終的な解法ではありませんが、部品として利用します。

基本的な考え方は、グラフの次数列が $ (d_1, d_2, \ldots , d_N) $ であるならば、ノイズを与えたあとのグラフの次数列の期待値は $ ( (1 - \epsilon) d_1 + \epsilon (N - 1 - d_1), \ldots, (1 - \epsilon) d_N + \epsilon (N - 1 - d_N) ) $ になるだろうということです。

注:この考え方は実は正しくありません。yunix さんのこちらのツイートを参照ください。ただ、ある程度の近似にはなるので、この記事では「ある程度の近似になっている」として進めます。

とりあえず、できるだけ多様な次数列を持つグラフを M 個構築したとして（方法は後述します）、クエリの答え方を考えます。

構築したグラフそれぞれについて、ノイズを与えたあとのグラフの次数列の期待値を求めておきます。そして、与えらえたグラフの次数列を計算し、ノイズを与えた後の次数列の期待値がそれに最も「近い」グラフを返せば良いです。「近い」の定義ですが、次数を降順にソートした後の次数列をベクトルと見てユークリッド距離を計算すれば良いでしょう。

この方法がうまくいくためには、できるだけ距離が離れている次数列を持つようなグラフを構築する必要があります。私は次のような方法でグラフを構築しました。

ランダムなグラフを大量に準備する。
まず、空グラフと完全グラフを選択するグラフの集合 S に追加する。
次の操作を S のサイズが M になるまで繰り返す : 準備したランダムなグラフのうち、S に追加されているグラフの次数列との距離の最小値が最大になるようなグラフを S に追加する。

ランダムなグラフの準備の仕方ですが、単純にランダムに辺を選択すると、なだらかな次数列を持つグラフばかりができてしまうので、偏りがある次数列を持つようなグラフをたくさん作ることが必要です。そのためには、様々な大きさのクリークや独立集合を組み合わせたグラフを作ると良いでしょう。行列表示だと以下のようなグラフです。

解法2: 辺の有無の差が最小になるようなグラフを選択する

解法1ではグラフの情報を次数列しか使っていません。次数列が同じグラフはたくさんあるので、これは多くの情報を使っていないことになります。

そこで、与えられたグラフに対して、辺の有無の差が最小になるようなグラフを選択することを考えます。

辺の有無の差を計算するには、ターゲットとなるグラフを固定した上で、頂点の対応付けを求める必要があります。これを簡単にするために、構築するグラフは次のような 3x3 のブロック行列に対応するようなグラフに限定することにします。

そうすると、頂点の対応付けは、各頂点がどのブロックに対応しているかだけを決めれば良いことになります。そのために、次のような山登りアルゴリズムを実装しました。

まず、頂点を次数順にソートして、ターゲットとなるグラフの次数順にソートしたブロックに対応させる
違うブロックに対応している任意の二頂点を選んで、入れ替えたときに一致する辺が増えるかどうかを判定し、増えるならば入れ替える。これを繰り返す。

このアルゴリズムを M 個のグラフに対して実行して、辺の有無の差が最小になるようなグラフを選んで出力すれば良いのですが、そうすると山登りに十分な時間を使うことができず、うまくいきません。

そこで、解法1を利用して、次数列だけに注目して候補を10個に絞ってから、上述の山登りアルゴリズムを実行して、辺の有無の差が最小になるグラフを選ぶという方法を実装しました。

このように、精度が低いが高速な方法で候補を絞ってから、精度が高い方法でランク付けをする手法は、実際の検索システムや推薦システムでよく使われている手法で、一般的に有効な考え方と言えるでしょう。

次に、M 個のグラフを構築する方法を考えます。解法1と同じように、集合から一番「遠い」グラフを逐次追加していく貪欲法を行いますが、今度は距離の定義を「辺の有無の差の数」で定義することにします。

ここで上述山登りアルゴリズムを使うと時間がかかりすぎてうまくいきません。そこで、ブロックグラフ同士のみを比較することに注目して、ブロックの並べる順序を $ 3! \times 3! $ 通り試して、辺の有無の差の数が最小値を計算しました（両方のグラフのブロック順序を反転しても変わらないことに注目すると半分にできます）。

その他の工夫

頂点のマッチングが成功しやすいグラフの作り方

上述の山登りアルゴリズムが成功するためには、次数順に頂点を対応させる初期状態がそれなりに良い必要があります。そのため、ブロック毎の次数の差が小さいグラフは採用しないようにしました。具体的には、少なくとも次数の差が N / 10 以上あるグラフのみを採用することにしました。

空グラフ・完全グラフとの混同が起きやすい問題

公式解説でも言及されていますが、空グラフにノイズを与えたグラフに、一部密度が高い部分ができてしまって、他のグラフと混同しやすいという現象があります。完全グラフについても同様です。

この問題に対処するために、構築するグラフ間の距離を計算する際に、片方が空グラフまたは完全グラフになる場合には、距離を0.5倍扱いするというハックを入れました。非常にアドホックなハックですが、スコアには有効で、コンテスト終了2日前に7位から1位に浮上することができました。

3x3 ブロック行列以外のパターン

$ \epsilon $ が小さい場合は 4x4 ブロック行列が有効な場合があったので、採用しました。また、$ \epsilon $ が大きい場合は 2x2 ブロック行列のみを使った方がスコアが高くなったので、そのようにしました。$ \epsilon $ が大きい場合はブロックが3個あると頂点のマッチングに失敗する確率が上がりすぎてしまったからだと思われます。

当たりの乱数ガチャを固定

各 $ M, \epsilon $ に対して最もスコアの期待値が高くなる N は事前計算で求めておいて固定します。これに加えて、スコアの期待値が高くなる乱数も固定します。上述のアルゴリズムでも、距離が同じグラフが複数ある場合にどのグラフを選ぶかで最終的に選ばれる M 個のグラフは変わってきます。そこで、乱数の seed に対して決定的に動作するようにアルゴリズムを実装した上で、各 $ M, \epsilon $ に対して最もスコアの期待値が高くなる seed を事前計算しました。これは細かい改善ですが計算リソースを必要とするため、最後の一押しとして実行しました。

2022年11月1日火曜日

AHC015 プレイアウトを使ったモンテカルロ法による3位解法の解説

2022年10月に開催されたトヨタ自動車プログラミングコンテスト2022 (AtCoder Heuristic Contest 015) の3位解法の解説をします。

まず、シンプルなルールベースによる解法を解説します。これだけで本番100位以内相当のスコアを出すことが可能です。

解法1: 上下に分割する

まず、イチゴ味（タイプ1）のキャンディーを上に集めることを目指します。

このためには、次のキャンディーがイチゴ味の場合は、下に傾けると良いです。すべてのキャンディーが下側に集められ、イチゴ味のキャンディーがその上に配置されるからです。

逆に、次のキャンディーがイチゴ味以外の場合は、上に傾けると良いでしょう。

このルールを実装すると次のような挙動になり、イチゴ味のキャンディーを上に集められていることがわかります。

この解法によって 104M のスコアを出すことができました（コード例）。これは本番286位相当のスコアで、これだけで正の得点を獲得した778人中の上半分に入ることが可能です。

解法2: 上・左下・右下に分割する

よりスコアを高めるためには、イチゴ味のキャンディーだけでなく、スイカ味（タイプ2）とパンプキン味（タイプ3）のキャンディーもできるだけ同じ場所に集めることが必要です。そこで、イチゴ味のキャンディーを上、スイカ味のキャンディーを左下、パンプキン味のキャンディーを右下に、できる限り集めることを考えます。

解法1と同じく、次のキャンディーの種類によって、傾ける方向を決めます。

次のキャンディーがイチゴ味の場合は、解法1と同様に下に傾ければ良いです。

次のキャンディーがスイカ味の場合は注意が必要です。

現在キャンディーが下側に集まっている場合（つまり、最後に下に傾けた後、一度も上に傾けていない場合）は、スイカ味のキャンディーをイチゴ味のキャンディーの上側に配置することを防ぐため、上に傾けましょう。この場合、スイカ味をパンプキン味よりも左側に配置されるようにすることはできませんが、仕方ありません。左右の関係よりも上下の関係を優先することにします。

一方、現在キャンディーが上側に集まっている場合（つまり、最後に上に傾けた後、一度も下に傾けていない場合）は、右側に傾けることで、スイカ味のキャンディーが左下に配置されるようにすることができます。

次のキャンディーがパンプキン味の場合も、スイカ味の場合と同様に、現在キャンディーが下側に集まっているか上側に集まっているかで場合分けしましょう。

このルールを実装すると次のような挙動になり、イチゴ味を上、スイカ味を左下、パンプキン味を右下に、おおまかにですが集められていることがわかります。

この解法によって 133M のスコアを出すことができました（コード例）。これは本番97位相当のスコアで、今回はこのヒューリスティックを発見できれば非常に短いコードで100以内に入ることが可能な回でした。何と盤面のシミュレーションすらしていないのです！

準備: 盤面のシミュレーションとスコア計算

より高度な解法に移る前に、盤面のシミュレーションとスコア計算が必要です。

盤面を傾ける操作は愚直に行ごとまたは列ごとに処理するしかないのですが、最終的にボトルネックになるため、不必要な遅い操作（メモリ確保など）が無いように気を付けましょう。

スコア計算では連結成分の計算が必要になりますが、深さ優先探索または素集合データ構造のどちらでも良いでしょう。

解法3: プレイアウトを使ったモンテカルロ法

では、3位解法となる、プレイアウトを使ったモンテカルロ法による解法を解説します。

ある時点で、次に盤面を傾ける向きを決めるときに、それぞれの向きに傾けたときの最終スコアの「期待値」を計算し、一番期待値が高い向きに傾けることを考えます。もちろん厳密に期待値を計算することはできないので、近似する方法を考えます。

まず、残りのキャンディの個数が N 個ある場合、有り得る残りの入力列は N! 通りあります。もちろんすべてを試すことはできないので、ランダムな入力列をいくつか生成してそれに対するスコアを計算し、平均値を取ることにします。

入力列を固定しただけではスコアを計算することができません。出力列（傾ける方向）も決める必要があります。今の目的は次の出力を決めることなので、次の出力だけを固定し、残りはランダムに選ぶという方針が考えられます。しかし、これではあまりにランダム性が高く、実際のスコアより大幅に悪いスコアになってしまうので、次の出力候補を比較するにはあまり役に立ちません。

そこで、次の出力だけを固定し、残りは解法2のアルゴリズムによって出力を決めることにしましょう。そうすると、高い精度で次の出力を固定したときのスコアの期待値を推定することができるようになります。これが可能なのは「解法2が非常に高速であり、かつ最適解にそれなりに近い」ためです。このように、高速なヒューリスティック（場合によってはランダム）によって出力を決定して最終状態のスコアを計算することを「プレイアウト」と呼びます。プレイアウトは、囲碁 AI などに用いられているモンテカルロ木探索において重要な概念です。興味のある方は調べてみてください。プレイアウトは、途中状態で盤面の良さを評価する評価関数を作るのが難しい場合に有効です。

この解法では、何個のランダム生成列を処理できるかによって性能が変わってきます。また、ゲームの序盤・中盤・終盤でそれぞれどの程度時間を使うかという時間管理も重要でしょう。私は毎回固定で150個のランダム生成列を処理するという方針を実装し、156M のスコアを出して3位を取ることができました。終了後にリファクタリングした後のコード例はこちらになります。

2021年6月6日日曜日

AtCoder Heuristic Contest 003: 10位解法の解説

2021年5月に開催された AtCoder Heuristic Contest 003 に参加し、10位を取ることができました。私の解法について解説します。

最終テストのスコアは2.917T、暫定テストのスコアは97.1Bでした。

問題

問題リンク

基本方針: 正規分布で近似して Maximum Likelihood Estimation

M = 2 のケースは複雑なので、まずは M = 1 のケースを考えることにします。また、D については D = 1000 と固定されていると仮定しておくことにします。

水平方向の辺について、行毎のパラメータが30個、辺毎のパラメータが30*29=870個あり、同様に垂直方向の辺について、列毎のパラメータが30個、辺毎のパラメータが870個あります。パスの長さはこれら1800個のパラメータの線形和として書けることに着目して、次のような確率モデルを考えます。

まず、1800個のパラメータはそれぞれ一様分布に従って生成されますが、これを扱いやすくするために、平均・分散が同じ正規分布で生成されることにするという近似をします。また、各クエリの結果についても、真のパスの長さを中心とする正規分布で生成されることにするという近似をします。すると negative log-likelihood が

各パラメータについて [パラメータの値とパラメータの期待値の差の二乗 / パラメータの分散]
各クエリについて、[クエリ結果と真のパスの長さの差の二乗 / クエリ結果の分散]

の和として書けます。これを最小化するようにパラメータを決めます。真のパスの長さはわからないので、線形モデルによる予測値で置き換えることにします。すると、これは Ridge 正則化有りの線形回帰と等価であることがわかります。これは機械学習で一般的なモデルですから、解くことができそうです。

Ridge 正則化有りの線形回帰の解き方

Ridge 正則化有りの線形回帰の解は closed form がある（英語版 Wikipedia: Ridge regression）のですが、これには逆行列の計算が必要なため、計算量がパラメータ数について三乗になってしまってこの問題では現実的ではありません。

そこで、最急降下法で解くことにします。色々なケースについてうまく行く学習率の設定は容易でないため、直線探索を使って更新量を決めることにします。私の場合は、ベクトルの更新量のノルムの最小値をあらかじめ決めておいて、それを二倍にすることを繰り返し、コストの改善量が最大になる値を採用するという方法でやりました。

コンテスト終了後の他の方の議論によると、共役勾配法を使うのが良い方法だったようです。

ここまでを実装することで95.4Bを獲得することができました。

注意: 途中経過の得点については、自分の実装ノートを基に復元しているので、不正確な場合があります。

M = 2 の扱い: 線形関数による近似

ここまで M = 1 を仮定してきが、M = 2 の場合を扱うことを考えます。

M = 2 では区間の境界という扱いにくいパラメータがあるため、線形関数で近似することにします。つまり、例えば水平方向の辺については、$ r $ 行目の左側の区間の辺の長さの「基本値」は $ H_{r, 0} $、右側の区間の辺の長さの基本値は $ H_{r, 1} $ となるのですが、これを左側から $ k $ 番目の基本値は $ H_{r, 0} + (k / 28) H_{r, 1} $ であるというように近似します。M = 1 の場合は \( H_{r, 0} = H_{r, 1} )\ の場合であるとみなすことができるので、とりあえずすべてのケースを M = 2 として扱うことにします。

ここまでを実装することで95.8Bを獲得することができました。

M = 1 と M = 2 の分類

先ほどはすべてのケースを M = 2 として扱いましたが、M = 1 の場合は自由度が減るので、M = 1 と判断できればより正確に推定することが可能なはずです。

そこで、300ターンが経過した時点で $ (H_{r, 0} - H_{r, 1})^2 $ と $ (V_{c, 0} - V_{c, 1})^2 $ の平均値を調べ、これが一定の値未満ならば M = 1、以上ならば M = 2 と分類することにしました。300ターンというターン数や分類の閾値は実験の結果から決めました。

これを実装することで M = 1 の場合のスコアが上がり、95.9Bを獲得することができました。

コンテスト期間終盤には、$ (H_{r, 0} - H_{r, 1})^2 $ と $ (V_{c, 0} - V_{c, 1})^2 $ の平均値がとても小さいまたは大きい場合には、100ターン経過時点でも十分正確に判断できることがわかったので、そのようなロジックも追加しました。

強化学習の考え方:Exploration に対する報酬

ここまで開発してきたプログラムによる辺の長さの予測値と真の値の差を調べてみると、ある行または列がほとんど選ばれなかった場合に差が大きくなることがわかりました。特に、ある行の基本値が2000などとても小さい場合であったとしても、初期のフェーズにほとんど使われなかった場合、基本値の期待値の5000とみなされてしまってその行がパスに使われず、基本値が2000であることを知るチャンスが最後まで失われてしまっているケースがあることに気付きました。

この問題は強化学習の問題であると捉えることができ、Exploration vs Exploitation のバランスを取ることが重要なのです。情報を集めることが重要な序盤においては、その時点のパラメータ推定に基づく最短路よりも、まだ使っていない行または列を使う方が、最終的にはパラメータ推定が上がってスコアも上がる場合も多いのです。このアイデアを実現するために、各行・列において選んだ辺の個数が $ k \leq 15 $ であるとき、その行・列の辺を使うことに対して $ 3000 (15 - k)^2 / 15^2 $ の報酬を与える、つまりグラフにおける辺の長さを報酬分短くする、というロジックを実装しました。

このようにすると、まだ使っていない行・列の基本値は報酬によって2000になります。すると、まだ使っていない行・列を使うために、大幅に遠回りになるパスが選ばれるという副作用によって序盤のスコアが大幅に下がることがわかりました。Exploration が大事とはいえ、パスに含まれる辺の個数が増えないような範囲で Exploration するのが良いと考え、序盤はすべての辺の長さに一定値を足すことで辺の個数が最小でないようなパスは選ばれにくいようにしました。具体的には、報酬を含めたグラフの辺の長さの平均値が5000になるように、すべての辺の長さに一定値を足すようにしました。

ここまでを実装することで96.9Bを獲得することができました。

M = 2 の場合の区間の境界を推定する

これまでのところ、M = 2 の場合に区間の境界を推定することなく、線形関数で近似してきました。問題設定に忠実に、区間の境界を推定し、それぞれの区間では基本値は定数であるというようにモデル化することを考えます。

以下、水平方向の行について考察します。垂直方向の列についても同様です。

境界の初期値は中央であるとして、そこから右または左に1ずつ動かしていく方針を考えます。境界の左のインデックスが $ x - 1 $、右のインデックスが $ x $ であるとします。また、$ H_{r, 0} < H_{r, 1} $、つまり右側の区間の基本値が高いと仮定します。

上の図に注目してください。真の境界が現在の境界より右にある場合は、境界の右側、つまりインデックス$ x $において真の基本値が現在の基本値の推定値より小さくなるので、その分のギャップが$ \delta_{r, p} $に反映されるとすると$ \delta_{r, p} < 0 $となりそうです。この考察に基づくと、$ \delta_{r, p} < 0 $である場合は、境界を右に動かして良さそうです。逆に、$ \delta_{r, p - 1} > 0 $ である場合は、境界を左に動かしても良さそうです。この直観に基づいて境界を1ずつ動かすことにします。両方の条件が成り立つ場合はどちらに動かすかを乱数で決めることにします。

このヒューリスティックでは真の境界まで動かすことができない場合も多いのですが、$ H_{r, 0} $と$ H_{r, 1} $の差が大きい場合はそれなりにうまく行くようです。そして、そのような場合こそが境界の推定が重要な場合なのです。

ここまでを実装することで97.1Bを獲得することができました。

より高得点な方々の記事

AHC003の2.926T解法+経緯 by 1位の eivour さん

Twitter 投稿 by 2位の yosss さん

2021年4月4日日曜日

全方位木DP (Re-Rooting) の抽象的定式化と実装：頂点情報と辺情報を扱う

この記事では、全方位木DP（英語では re-rooting と呼ばれることが多いようです）と呼ばれるアルゴリズムを、できるだけ多くの問題を扱えるように抽象的に定式化しつつ、ライブラリとして使える実装を示すことを目的とします。全方位木DPを知らない方はまず keymoon さんによる【全方位木DP】明日使える便利な木構造のアルゴリズムをご覧下さい。

全方位木DPによって解ける問題とは、どのように定式化できる問題でしょうか？

以下のように書ける問題は、全方位木DP によって解くことができます。無向木 $ T = (V, E) $ が与えられるとする。また、$ T $ を頂点 $ i $ を根とする根付き木として見たものを $ T_i $ とする。ここで、根付き木の部分木 $ S $ について次のように再帰的に定義される関数 $ F(S) $ を考えます。

$$ F(S) = \textrm{add_vertex} (G(S_1) \bullet \cdots \bullet G(S_m)), v) $$

$$ G(S_i) = \textrm{add_edge}(F(S_i), e_i) $$

但し

$ v $ は $ S $ の根
$ S_i $ は $ v $ の $ i $ 番目の子 $ v_i $ を根とする部分木
$ e_i $ は $ v $ と $ v_i $ を接続する辺

とします。

このとき $ F(T_i) $ をすべての頂点 $ i $ について求める問題が、全方位木DPによって解ける問題となります。

関数 $ F $ は次の三つの関数と一つの値によって定まることがわかります。辺情報を処理する関数 $ \textrm{add_edge} $、子部分木の情報をマージするモノイド上の演算 $ \bullet $ とその単位元 $ \textrm{id} $（モノイドについては上記の keymoon さんの記事を参照してください）、頂点情報を処理する関数 $ \textrm{add_vertex} $ です。従って、全方位木DPのアルゴリズムは、これら三つの関数及びモノイドの単位元を受け取るような関数として、ライブラリとして再利用できる形で実装することができます。実装は記事の最後に示します。

具体的な問題にあてはめることによって、様々な問題をこの定式化によって扱うことができることを見ていきましょう。

問題1: AOJ 木の直径 (問題リンク)

重み付き木における木の直径を求めてください。

この問題は全方位木DPを使わなくても解くことができますが、練習のため全方位木DPで解いてみましょう。

関数 $ F(S) $ を $ S $ の根から最も遠い頂点までの距離としましょう。$ \max(\{ F(T_i) | i \in V \}) $ が木の直径となります。

$ \textrm{add_edge}(x, e) = x + [e\textrm{'s length}] $
$ x \bullet y = \max(x, y) $
$ \textrm{id} = 0 $
$ \textrm{add_vertex}(x, v) = x $

とすると、関数 $ F(S) $ を定義することができます。

この問題では頂点情報は必要ないため $ \textrm{add_vertex} $ は自明な関数となっています。

ライブラリを使うと以下のように全方位木DPを実行することができます。

ReRooting<int64, int64>(
      graph,
      /* add_edge */ [](int64 x, const Edge e) -> int64 { return x + e.w; },
      /* monoid */ [](int64 x1, int64 x2) -> int64 { return max(x1, x2); },
      /* monoid id */ 0,
      /* add_vertex */ [](int64 x, int v) -> int64 { return x; });

問題2: AtCoder EDPC V Subtree（問題リンク）

木が与えられます。各頂点を白または黒に塗りますが、黒い頂点が連結になるようにします。このとき、各頂点 $ v $ が黒になるような塗り方は何通りあるかを $ \mod M $ で求めてください。

関数 $ F(S) $ を「部分木 $ S $ の根を黒く塗り、かつ黒い頂点が連結になるような色の塗り方の数」とします。各頂点 $ v $ に対して $ F(T_v) $ が求める値となります。

$ \textrm{add_edge}(x, e) = x + 1 $
$ x \bullet y = x \cdot y $
$ \textrm{id} = 1 $
$ \textrm{add_vertex}(x, v) = x $

とすると、関数 $ F(S) $ を定義することができます（詳細な解説は kyopro_friends さんによる解説を参照して下さい）。

ライブラリを使うと以下のように全方位木DPを実行することができます。

ReRooting<ModInt, ModInt>(
     graph,
      /* add_edge */ [](ModInt x, Edge e) -> { return x + 1; },
      /* monoid */ [](ModInt x, ModInt y) -> { return x * y; },
      /* monoid id */ 1,
      /* add_vertext */ [](ModInt x, int v) -> { return x; });

TODO: 頂点情報を使う問題を発見したら例として足す

ライブラリの実装例

struct Edge {
  const int from, to;
  Edge(int from, int to) : from(from), to(to) {}
};

template<class T, class U> vector<T> ReRooting(
    const vector<vector<Edge>>& graph, function<U(T, Edge)> edge_fn,
    function<U(U, U)> merge_fn, U merge_id, function<T(U, int)> vertex_fn) {
  int n = graph.size();

  // BFS
  vector<int> vs({0}), ps(n, -2); ps[0] = -1; queue<int> q; q.push(0);
  while (!q.empty()) {
    int v = q.front(); q.pop();
    for (const auto& e : graph[v]) {
      if (e.to == ps[e.to] || ps[e.to] >= -1) continue;
      vs.push_back(e.to); ps[e.to] = v; q.push(e.to);
    }
  }

  // DP phase 1.
  vector<T> dp_1(n);
  for (int i = n - 1; i >= 0; i--) {
    int v = vs[i]; int p = ps[v];
    auto a = merge_id;
    for (const auto& e : graph[v]) {
      if (e.to == ps[v]) continue; a = merge_fn(a, edge_fn(dp_1[e.to], e));
    }
    dp_1[v] = vertex_fn(a, v);
  }

  // DP phase 2.
  vector<T> dp_2(n), dp_3(n); dp_2[0] = merge_id;
  for (int i = 0; i < n; i++) {
    int v = vs[i]; int p = ps[v]; int m = graph[v].size();
    vector<T> ls(m + 1), rs(m + 1);
    ls[0] = merge_id;
    for (int j = 0; j < m; j++) {
      const auto& e = graph[v][j];
      if (e.to == p) ls[j + 1] = ls[j];
      else ls[j + 1] = merge_fn(ls[j], edge_fn(dp_1[e.to], e));
    }
    rs[m] = merge_id;
    for (int j = m - 1; j >= 0; j--) {
      const auto& e = graph[v][j];
      if (e.to == p) rs[j] = rs[j + 1];
      else rs[j] = merge_fn(rs[j + 1], edge_fn(dp_1[e.to], e));
    }
    dp_3[v] = vertex_fn(merge_fn(ls[m], dp_2[v]), v);
    for (int j = 0; j < m; j++) {
      const auto& e = graph[v][j]; if (e.to == p) continue;
      for (const auto& re : graph[e.to]) {
        if (re.to != v) continue;
        T merged = merge_fn(merge_fn(ls[j], rs[j + 1]), dp_2[v]);
        dp_2[e.to] = edge_fn(vertex_fn(merged, v), re);
      }
    }
  }
  return dp_3;
}

2020年1月27日月曜日

AtCoder DISCO presents ディスカバリーチャンネルコードコンテスト2020 本戦問題B - Hawker on Graph

問題リンク

公式解説と同じ解法ではあるのですが、より素朴な発想で同じ解法に至ったので、記録しておきます。

まず、この問題の難しいところは「所持金が負となることはない」というところで、これが無ければ有名問題で、以下のようにして解けます。

$ T $ ターンかけて $ i $ から $ j $ に移動するときの利得を $ A_T [i][j] $ とします。ここで、$ A_{T_1} [i][j] $ と $ A_{T_2} [i][j] $ が各 $ i, j $ についてわかっているとします。すると、$ T_1 + T_2 $ ターンかけて $ i $ から $ j $ に移動するとき、$ T_1 $ ターン時点で $ k $ を経由する場合を各 $ k $ について考えれば良いので、$ A_{T_1 + T_2}[i][j] = \max_k \{ A_{T_1}[i][k] + A_{T_2}[k][j] \} $ と書けます。これによって、繰り返し二乗法を使って $ A_K $ を求めることができます。計算量は $ O(N^3 \log K) $ となり、$ N \leq 200 $ なので十分高速です。

では、「所持金が負となることはない」という条件を考えていきます。初期の所持金が $ W $ のとき、$ T $ ターンかけて $ i $ から $ j $ に移動することを考えます。ひとまず「所持金が負となることはない」という条件を無視して、辺の重みの累積和の推移をグラフにします。

すると、遷移後の所持金は次のように書けます。

途中で所持金が $ 0 $ にならない場合は $ W + A_T [i][j] $ になります
途中で所持金が $ 0 $ になる場合は、累積和が最小値になる頂点が最後の所持金が $ 0 $ になる頂点になることに注目します。そうすると、累積和の最小値を $ A_T [i][j] $ から引いた値を $ B_T [i][j] $ として、最終的な所持金は $ B_T [i][j] $ になります

したがって、遷移後の所持金の最大値は $ \max( W + A_T [i][j], B_T [i][j] ) $ となります。

あとは、$ A_K $ と同様に $ B_K $ が繰り返し二乗法で計算できれば OK です。$ B_{T_1} [i][j] $ と $ B_{T_2} [i][j] $ が各 $ i, j $ についてわかっているとします。$ T_1 $ ターン時点で $ k $ を経由して$ T_1 + T_2 $ ターンかけて $ i $ から $ j $ に移動するパスを考えます。

このとき $ B_{T_1 + T_2}[i][j] $ の値は

累積和の最小値が $ i - k $ 間にある場合、$ B_{T_1}[i][k] + A_{T_2}[k][j] $
累積和の最小値が $ k - j $ 間にある場合、$ B_{T_2}[k][j] $

になるので、$ B_{T_1 + T_2}[i][j] = \max_k(B_{T_1}[i][k] + A_{T_2}[k][j], B_{T_2}[k][j]) $ と書けます。これによって、$ A_K, B_K $ に対して同時に繰り返し二乗法を適用できることがわかりました。

#include <bits/stdc++.h>
 
using namespace std;
using int64 = long long;
using vvi64 = vector<vector<int64>>;
 
constexpr int DEBUG = 0;
 
template<typename T>
vector<vector<T>> Make2DVector(int d1, int d2, T default_value) {
  return vector<vector<T>>(d1, vector<T>(d2, default_value));
}
 
template<class T> inline bool UpdateMax(T& a, T b) {
  if (a < b) { a = b; return 1; } return 0;
}
 
int main() {
  ios::sync_with_stdio(false);
  cin.tie(0);
 
  int n, m;
  cin >> n >> m;
  int64 initial_score;
  cin >> initial_score;
  int s;
  cin >> s;
  s--;
  int64 k;
  cin >> k;
 
  constexpr int64 INVALID = -1E18;
  auto a = Make2DVector<int64>(n, n, INVALID);
  auto b = Make2DVector<int64>(n, n, INVALID);
  for (int i = 0; i < m; i++) {
    int from;
    int to;
    int64 w;
    cin >> from >> to >> w;
    from--;
    to--;
    a[from][to] = w;
    b[from][to] = max(w, 0LL);
  }
 
  // Computes A_{T_1 + T_2} and B_{T_1 + T_2} from
  // A_{T_1}, B_{T_1}, A_{T_2} and B{T_2}.
  auto merge_fn = [&](
      const vvi64& a1, const vvi64& b1,
      const vvi64& a2, const vvi64& b2) -> tuple<vvi64, vvi64> {
    auto a = Make2DVector<int64>(n, n, INVALID);
    auto b = Make2DVector<int64>(n, n, INVALID);
    
    for (int i = 0; i < n; i++) {
      for (int j = 0; j < n; j++) {
        for (int k = 0; k < n; k++) {
          if (a1[i][k] == INVALID) continue;
          if (a2[k][j] == INVALID) continue;
 
          UpdateMax(a[i][j], a1[i][k] + a2[k][j]);
          UpdateMax(b[i][j], b1[i][k] + a2[k][j]);
          UpdateMax(b[i][j], b2[k][j]);
        }
      }
    }
    return make_tuple(a, b);
  };
 
  // Computes A_p and B_p from A and B by doubling.
  function<tuple<vvi64, vvi64>(const vvi64&, const vvi64&, int64)>
      power_fn =
          [&](const vvi64& a1, const vvi64& b1, int64 p)
              -> tuple<vvi64, vvi64> {
    if (p == 1) return make_tuple(a1, b1);
    vector<vector<int64>> a2, b2;
    tie(a2, b2) = power_fn(a1, b1, p / 2);
    tie(a2, b2) = merge_fn(a2, b2, a2, b2);
    if (p % 2 == 1) {
      tie(a2, b2) = merge_fn(a2, b2, a1, b1);
    }
    return make_tuple(a2, b2);
  };
 
  vector<vector<int64>> ap, bp;
  tie(ap, bp) = power_fn(a, b, k);
 
  int64 ans = -1;
  for (int t = 0; t < n; t++) {
    if (ap[s][t] == INVALID) continue;
    UpdateMax(ans, initial_score + ap[s][t]);
    UpdateMax(ans, bp[s][t]);
  }
 
  cout << ans << endl;
}

登録: 投稿 (Atom)

ymatsux@競技プログラミング

2022年11月26日土曜日

HTTF 2023 予選 (AHC 016) 3位解法の解説

\( \epsilon = 0 \) の場合

解法1: 次数列で分類する

解法2: 辺の有無の差が最小になるようなグラフを選択する

その他の工夫

2022年11月1日火曜日

AHC015 プレイアウトを使ったモンテカルロ法による3位解法の解説

2021年6月6日日曜日

AtCoder Heuristic Contest 003: 10位解法の解説

2021年4月4日日曜日

全方位木DP (Re-Rooting) の抽象的定式化と実装：頂点情報と辺情報を扱う

2020年1月27日月曜日

AtCoder DISCO presents ディスカバリーチャンネルコードコンテスト2020 本戦問題B - Hawker on Graph

HTTF 2023 予選 (AHC 016) 3位解法の解説

不正行為を報告

2022年11月26日土曜日

HTTF 2023 予選 (AHC 016) 3位解法の解説

\( \epsilon = 0 \) の場合

解法1: 次数列で分類する

解法2: 辺の有無の差が最小になるようなグラフを選択する

その他の工夫

2022年11月1日火曜日

AHC015 プレイアウトを使ったモンテカルロ法による3位解法の解説

2021年6月6日日曜日

AtCoder Heuristic Contest 003: 10位解法の解説

2021年4月4日日曜日

全方位木DP (Re-Rooting) の抽象的定式化と実装：頂点情報と辺情報を扱う

2020年1月27日月曜日

AtCoder DISCO presents ディスカバリーチャンネル コードコンテスト2020 本戦 問題B - Hawker on Graph

HTTF 2023 予選 (AHC 016) 3位解法の解説

AtCoder DISCO presents ディスカバリーチャンネルコードコンテスト2020 本戦問題B - Hawker on Graph