AtCoder ABC044 C 高橋君とカード動的計画法

　前にも書いたように最初まったくわからず解法にあった問題を尋ねて遠まわりしてたがようやくできたので書く

問題
迷宮入り
動的計画法
解法1
解法2
応用篇:オレ解法

問題

数列から選んだ合計の平均がAになる組はいくらあるか。数列の規模50以下。数列は50以下の正の整数。

迷宮入り

　お、これは組み合わせだな？ということで、早速数列から平均値を引いたリストをつくって、とりあえず辞書につっこんで数字別に数えあげたものの、そこから先がよくわからない。はて。そこから0を挟んで対称となる組み合わせと0の組み合わせをかけていけばよさそうだが、さてどうやってそれを出せばいいのか。うぅーむ.... わからん。

f:id:inudaisho:20181011155434p:plain

　長考してもわからないのでとりあえず解法をみてみると解法1として bit演算とか書いてある。で、それの参考になるというので遡って解いたりしたあげくグラフの初歩にひっかかったのだが、それをとりあえずなんとかしてからまた戻ってきた。ためしにビット演算でやってみるとPython3では例4を手元で計算させるだけでもTLEする。これはやってもしかたない。そこで解法2,解法3をみると、どうもこちらはよく聞く動的計画法らしい。

動的計画法

　しばらく眺めてこれでどうして組み合わせがでてくるのかよくわからず、いろいろいじっているうちに、sum(k=1→n) nCk が2ⁿ -1 と等しいらしいということを見付けた。例4がそれだったのはたぶんそれもわかる人にわかるようにしてあるヒントなんだろう。まぁ自分はみただけではちっともわからなかったが。それである程度は積算していったら何かそれらしいものが出るらしいことがわかったが、ではこういうふうにする漸化式をどうやってひねり出すのか。

　動的計画法はリチャード・ベルマン(アメリカ)が1953年発表したものらしい。なんか詳しい本がありそうだったが夜だったので待てずにネットでなんとかした。幸いなことに Wikipedia がわりと詳しく書いてあった。ちなみに動的計画法 - Wikipediaではない。これはプログラム屋が動的計画法として知られているプログラム技法を紹介しているだけなので、専門の穴にハマっている(井戸の底から見る空は狭い)。大枠を知った後で読むものだ。最初に読むべきはこっちだ。

　このベルマン方程式の解説で役に立つのが、「動的計画法の考え方」だ。目標を数式化した目的関数。表現するための状態変数・制御変数。コントロールする政策関数。それらに別軸で評価をつける価値関数。この価値関数の最適化が動的計画法でやることらしい。それらが最善手で進むという前提で最適化の漸化式と初期値を立てるとあとは決まるということらしい。ベルマン方程式の下の方をみるとマルコフ過程とか書いてある。あ、そういうのの基礎か。

R. ベルマンと動的計画法あらがき - 日本オペレーションズ・リサーチ学会

　これは1985年(昭和60)のベルマン追悼記事だが、その中で動的計画法についてこんなふうにまとめてある。

特に重要なことは多くの場合DPの手法は古典的手法のような数学的訓練を必要とせぬという事実である。さらに古典的理論の大部分は数値計算に適するようには作られていない。しかしコンピュータが発達したので、コンピュータの能力にもとづいた新しい数学が創られつつある。まさにDPは今世紀が生んだ嫡子である。

DPの本質は政策の概念であり、それはシステムの状態の項から作った決定を語る法則である。

政策の思想は実験による学習の概念に通じている。

高水準の意思決定過程の理論はこれまでの動的計画法を含み、しかもこれから出現する高性能のコンピュータに相応ずる理論でなければならない。

　なるほど、この単純な漸化式によるミニ試行が、複雑な問題に答えを出してしまうことが、今のAI大流行のような事態を導いたわけだな。動的計画法の単純な升目がたぶん当時の人には今のニューラルネットワークのブラックボックスみたいに見えたんだろうな。競技プログラミングをやると強くなるかどうかは別にして、今現在進行中の技術の基礎を押さえれるわけだ。なるほど。

　それからWikipediaのベルマンの項に逸話のようにして引いてある「次元の呪い」、これは結局機械でぶんまわすといっても要素が多いとつらくなるということらしい。

解法1

　なんだ結局場合の数は出せず試行で足しあわせて出すのか、とおもいながら動的計画法にとりかかる。ちなみにこういう問題は動的計画法の典型的な問題らしい。

　解法1のDPは縦N(数列j番目)横N(k枚)高さN*X (Xは数列とAの最大値)(そのときの合計)の場合の数を入れる立方体の箱をつくり、k枚のときKAとなる場合の数を合計するということで、自分が最初考えていて詰まった解法よりも原始的である。計算量がものすごいというのはわかるが、まだよくわからず、とりあえず解法にあるとおりにループを組んでみると、手元ではたしかに見事に解ける。なるほど。感心してからAtCoder の方に投げてみるとTLE。式をみなおして無駄な処理を減らしてもTLEだったので Python3 はこういうことするための言語じゃないなとおもいながら解法2に移った。

　ちなみにいろいろやってかなり飲みこんできてから、後で見直したのだが、

dp[0][0][0] = 1
for j in range(1,iN+1):
#   for k in range(0,iN+1):
    for k in range(0,j+1):
#       for s in range(0,iN*iMX + 1):
        for s in range(0,j*iMX + 1):
            if  s < aX[j-1]:

　こうやって立方体を斜めに二回切り取ることができる。そうすると解法1の愚直な for ループでも通る。

解法2

　解法2は最初自分が思いついたように差分をとって座標を0に持ってくるやりかただ。まぁよくある手だから当然だが、そこから二次元の表を作って全部計算するのは違う。縦N(J枚目)横2NX(t-NXのt)。0を中心にしたのはいいが、結局正と負の二面作戦になるので両面にNXが必要になるらしい。真ん中に0を持ってくるために t-NXという小細工をしたわけだ。

f:id:inudaisho:20181011151842p:plain

　最後の結果をdp[N][NX]でとるというのはつまりそこに結果0の場合の数が集まっているからだ。(真ん中の0の列がdpのNXの列)同様に最下段には-NX...-2,-1,1,2,...NXのときの場合の数が入っている(一枚も選ばなかった場合も含んでいるので答えを出すときは1引く)ので自分の思いついた解法にも応用できるが、このまま使うと当然遅くなるのでこのままの応用はしなかった。適当に実装するといきなり300msくらいで計算できた。実は実装するのは簡単で、なぜこれで答えが出るのかわからず、結局dpのリストに何が入っているのか出力させてみてようやくイメージがつかめたのだが、それはさておき。他の人の解答をみると、2桁msの速さで計算させている人たちがいる。その中身をみると collections の Counter をつかっている。これはつっこんだリストの要素を数えてハッシュの形で保持してくれるもので、便利な性質として、無い要素はエラーではなく0で返してくれる事と、辞書(ハッシュ、連想配列)で渡すとキーを合わせて足し算してくれる事の二つがある。Pythonはリストをぶんまわすのが遅いので動的計画法でも表の最後の行だけ使う今回のようなものであればそういう仕組みで足し算していってもいいわけだ。なるほど賢い。とはおもったが解法2をそのまま python のライブラリつかって高速化してもおもしろくないので、今回は最初に思いついたロジックを実装することにした。

応用篇:オレ解法

　最初に書いた通り、ゼロだけ分けて組み合わせを計算し、正と負の数列に分けて正負同じ結果になる組み合わせを探すというやりかたについてだが、これに動的計画法が適用できそうだ。解法だけ見ても結局自分の頭を使って応用しないと身につかない。ということで早速応用篇になる。

　正と負に分けるとその上限が正なり負なりのそれぞれの合計の小さい方だということがわかるので正負に分けると必要な計算量がものすごく減る。1以上の数字を足し合わせて行くだけだから、リミットを越えた数字は無視してよい。ということで、縦を正負いずれかの数列、横をさきほどの合計の大きい方で取る。これ箱の大きさだけは大きい方を採用したのは小さいのだと1になることがあってエラーが出る。

　漸化式は求めたい合計よりも数列の数字が小さければ前の数列のときの場合の数を取り、大きければ、前の数列のとき(不採用のとき)に加えて前回が今求めたい合計より数列の数字分引いた合計のとき(採用のとき)の場合の数を足す。それを順次やっていくのだが重要なのは合計0の箱に1をいれておくことで、それが最初にその数字が出てきた場合の種になる。

        for j in range(1,iLenPlus):
            for k in range(iUlim + 1):
                if k < aPlus[j] :
                    dp1[j][k] = dp1[j-1][k]
                else:
                    dp1[j][k] = dp1[j-1][k] + dp1[j-1][k-aPlus[j]]

　つくったdpの表の中身を出してみるとこうなる。(入力例3)

#--aPlus
[1, 3, 2, 1, 4]
#-dp1-
# 0,1,2, 3, 4, 5
[1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] # 1 -> 0と1に1
[1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0] # 3 -> 3と4に1 (0+3 , 1+3
[1, 1, 1, 2, 1, 1, 0, 0, 0, 0, 0, 0] # 2 ->という具合
[1, 2, 2, 3, 3, 2, 0, 0, 0, 0, 0, 0] # 1
[1, 2, 2, 3, 4, 4, 0, 0, 0, 0, 0, 0] # 4
#-aMinus-
[2, 3]
#-dp2-
[1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0]

　このそれぞれの最後の行の数字と0の場合の数(0の場合のときを考えて1足す)をかけて0の時の場合の数を足したら答が出る。ABCの方のPython3でこのやりかたで最速(25ｍs)を取ったのでちょっと自慢したくてこんなの書いたのだが、ARCの方も同じ問題出てるのに気付いてそっちを確認してみると普通にもっと速く書いてる人がいた。うぅーむ。

　とくにこのへん。

提出 #1780666 - AtCoder Regular Contest 060

提出 #1389625 - AtCoder Regular Contest 060

    dp = {0:1}
    for n in a:
        for k, v in dp.copy().items():
            dp[k+n] = dp.get(k+n, 0) + v

　すげー Counter とかいらんかったんや！

python3 の dict のgetはキーがないときエラーではなく、第2引数を返してくれる
copy()を挟むことによって辞書サイズ変更エラーがなくなる

　このテクをパク.....応用することによってARC/ABCを併せての最速(19ms)出しました(厚顔無恥

　ちなみにこれを応用するとどんなループになるかはこちら参照。上のリストと比較するとわかりやすいとおもいます。上のリストが上限切ってそれ以上計算していない以外は一緒ですね。

iX: 1
 k,v: 0 1 |k+iX: 1
iX: 3
 k,v: 0 1 |k+iX: 3
 k,v: 1 1 |k+iX: 4
iX: 2
 k,v: 0 1 |k+iX: 2
 k,v: 1 1 |k+iX: 3
 k,v: 3 1 |k+iX: 5
 k,v: 4 1 |k+iX: 6
iX: 1
 k,v: 0 1 |k+iX: 1
 k,v: 1 1 |k+iX: 2
 k,v: 3 2 |k+iX: 4
 k,v: 4 1 |k+iX: 5
 k,v: 2 1 |k+iX: 3
 k,v: 5 1 |k+iX: 6
 k,v: 6 1 |k+iX: 7
iX: 4
 k,v: 0 1 |k+iX: 4
 k,v: 1 2 |k+iX: 5
 k,v: 3 3 |k+iX: 7
 k,v: 4 3 |k+iX: 8
 k,v: 2 2 |k+iX: 6
 k,v: 5 2 |k+iX: 9
 k,v: 6 2 |k+iX: 10
 k,v: 7 1 |k+iX: 11