更新内容
面白そうなパズルとの出会い
Twitter で A Puzzle A Day というカレンダーパズルを見かけた。
明日の答えです（閲覧注意） pic.twitter.com/jzU5CbZWDU
— Takeshi HASEGAWA (@hasegaw) May 26, 2021
このパズルは1月〜12月、1日〜31日までのセルがあり、任意の日付（たとえば5月28日だとか6月1日だとか）をパズルを通して表現できる、というもののようだ。テトリスにでてきそうなブロックを組み合わせて、目的の日付にかかわる部分以外をキレイで埋められれば完成。毎日パズルが楽しめるし、1年後に同じ日付の答えを得ようとしても改めて楽しめそうで、とても面白そうだ。
久々に「このパズルは私も欲しい！」と思った。
国内の販売店が普通に売ってるようなものだったらうっかりオーダーしていたかもしれない。どうやらノルウェーのほうでレーザープリンタを用いて作ってるような感じで、持っている人は少ないのではないか。 Twitter タイムラインをみると紙やら木材のレーザーカットやら思い思いの方法で手作りをして楽しんでいる人もいるようで、私も何か作ろうと思った。
その結果として、ソルバーができた。
https://github.com/hasegaw/puzzle-a-day-solver
ソルバー？
ソルバー（Solver）というのは答えを求めるツールだ。たとえば今回のソルバーの場合、5月28日であればこうすると良い感じにブロックがはまりますよ、というのを教えてくれる。このパズルは365通りの目標があるほか、それぞれに複数の解法があるので、趣味などとしてソルバーを書いてみても楽しいだろうと思って、作ってみた。
私はコンピュータサイエンスを学んだエンジニアでも競プロ勢でもアルゴリズムに強いわけでもないので、ここで紹介するアルゴリズムについてはツッコミどころもあると思う。この記事の内容は、とりあえず趣味で簡単なプログラムを書いているだけでも、このパズルのソルバーはこれぐらいの考えで作れる内容だよ、という紹介ということでご理解いただきたい。
ルールの確認
このパズルのボードは下記の構造になっている。
左が、パズルの舞台となる盤面で、ここにブロックを詰めていくことになる。黒色のセルはブロックが置けない場所で、オレンジ色のセルは、一例として5月28日を表現する場合においてブロックで隠れてはいけないセルだ。
右側に示した 8つのブロックを盤面に並べていく。各ひとつ使えるのだが、各ブロックは回転させても良いし、反転させても良いようだ。なので、ひとつのブロックでも盤面上では最大で8つの姿をとりうる。
盤面の黒色セルおよびオレンジのセルに重ならないように、これらの8つのブロックを配置できたらパズルが解けたことになる。
基本的な考え方
ボードサイズは7x7のマトリックスで表現できる。ここで、ブロックを置いてはいけないセルを 1 、それ以外を 0 としたマトリックスを作る。ルール上絶対にブロックが乗ってはいけない部分のみ 1 としたマトリックスの例、それに加えて目的の日付セルも 1 としたマトリックスの例を示す。
ここで黒色セルと日付セルをわけて考えてもよいのだが、候補を見つけていくには、黒色セルおよび日付セルの両方がフラグされている後者のマトリックスが基点になる。
これに対して、特定のブロックを特定の回転・反転・位置の条件で配置する場合に、どのセルにあたるかを示す行列（以後、候補とよぶ）を作成する。ここではふたつの例を示す。
候補がボードにはまるかどうかは、現状のボードの状況を示すマトリックスと候補のマトリックスを加算してみればわかる。整数値の加算結果として 1 を超えるセルが発生する場合は、ブロックと何かがぶつかってしまい、そこに当てはめられない事を意味する。（ここで、後の最適化で触れる論理積を使ってもいい）
当てはめられる候補を見つけたら、再帰的に他のブロックも当てはまる所を探していって、全てのブロックを使い切ったら解法が見つかったことになる。
実際のところ、これを適当に置いていってもブロックが素直に収まってくれる可能性はかなり低い。黒色セルとオレンジの日付セルを除いて全てのセルをブロックで埋めないと答えにたどり着かないことを利用して、マトリックスのふち側のセルから優先して埋める。これにより、隙間の発生により絶対に解法にたどり着けないことが自明な候補を探索対象から外すことができる。
私の実装では左上から右へ、行が埋まったら次の行へ、という順番でセルを埋めるようにしている。次に埋めるセルを見つけるには、ボードの最新盤面におけるゼロのセルを見つければよい。候補がそのセルを埋めてくれるかは、その候補において目的のセルが 1 になっているかを見ればよい。候補が目的セルを埋めなくても解法に繋がる可能性はあるが、後に別のセルを埋める際に改めて評価されるので、ここでは候補から除外してよい。とにかくフチから埋めて、ボード上に隙間ができない候補から試したほうがよいからだ。
これを繰り返していくと、ボード上のセルでゼロがなくなり（すべて 1 になる）、未使用ブロックもなくなっているはずだ。この方法で総当たりしていけば、解法にたどり付ける。
実装については Pyhthon + NumPy を利用した。 NumPy で上記において必要となるマトリックス演算のほとんどが実現できるほか、この後に触れるベクタライゼーションを意識した記述にすることで、さらなる高速化が期待できる。
ベクタライゼーション (Vectorization)
上記のアルゴリズムは、どのようなプログラムで書いてもよい。私は Python を使用したが、別に Ruby を使ってもいいし、 JavaScript/TypeScript でも、 Perl でも、 Golang でも C でも Pascal でも、 C# でも、 Java でも、 BASIC でもなんでもよい。
ただ、このような処理を愚直にループ処理で書いても、コンピュータはこれを速く処理ができない。このため、コンピュータが高速に処理できるような、まとまった行列演算にできるだけ落とし込んでいく。たとえば盤面に対して次の候補を探すためにマトリックスを加算する場合に、まとめて加算をしたほうが速い。
もちろん加算すべきは ボードの現状 × 候補の数 なので計算すべき量は一緒なのだが、 Python コードとして繰り返しを書くよりも、 NumPy などの C 言語ベースで書かれたライブラリ内でループをしたほうが速い。場合によっては、バックエンドのライブラリ内において、 CPU の SIMD 命令で計算することにより、より実行時間を短縮してくれる可能性がある。
具体的には、 10x10 のボードに対して 10x10 の候補が 10,000 あった場合、 Python レベルで10,000回加算するよりも、 10000x10x10 の候補に対して 10x10 のマトリックスを加算するほうが速く実行できる。
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import cProfile
import numpy as np
mat1 = np.zeros((10000, 10, 10), dtype=np.uint8)
mat2 = np.zeros((10, 10), dtype=np.uint8)
mat1[:, :, :] = 3
mat2[::, :] = 3
def func1():
mat1 + mat2
def func2():
for m in mat1:
m + mat2
if __name__ == '__main__':
cProfile.run('func1()')
cProfile.run('func2()')% python bench.py
4 function calls in 0.000 seconds
Ordered by: standard name
ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 0.000 0.000 :1()
1 0.000 0.000 0.000 0.000 bench.py:12(func1)
1 0.000 0.000 0.000 0.000 {built-in method builtins.exec}
1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects}
4 function calls in 0.005 seconds
Ordered by: standard name
ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 0.005 0.005 :1()
1 0.005 0.005 0.005 0.005 bench.py:15(func2)
1 0.000 0.000 0.005 0.005 {built-in method builtins.exec}
1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects}
また、大抵の場合、 for ループなどを使用し最小値を探すよりも、 NumPy.argmin() などのライブラリ関数を活用したほうが速い。
このような観点でベクタライゼーションをある程度進めたところ、当初手元のコンピュータ(MacBook Pro 2019, 2.5GHz)では 600秒(10分)かかっていた探索が60秒以内で完了するようになった 。
うっかり vectorize して、昨晩から20倍ぐらい速くなった pic.twitter.com/Z3IJzyMgaJ
— Takeshi HASEGAWA (@hasegaw) May 26, 2021
バイナリ化
ベクタライゼーションにより特定の日付向けの解法が60秒以内に得られるようになったが、もう少し速くしてみようかなと思い、内部表現をバイナリ化してみることにした。
ここまでの作業は (7, 7) の大きさのマトリックスで、各要素は 8bit の整数値として扱ってきた。しかし 7x7 というサイズで、上記アルゴリズムでは基本的に 0 と 1 のみ（衝突検出のみ 2 が登場）で成り立っているので、2進数で表現できる。これにより 7x7=49バイトで表現していた盤面や候補が8バイトで表現できるようになり、メモリ上の表現が7分の1になること、uint64（64ビット符号なし整数）表現でのビット演算は現在のプロセッサだと1命令で実行できる。大幅な高速化が望めるだろう。
盤面の各ビットを以下のように割り当てることにした。最上位ビットから利用して最下位ビット側を未使用にしているのは若干気持ち悪く見えるかもしれないが、当初はブロックをシフトした際の範囲外判定のために右側のビットを残しておこうと思ったのでこのような割り当てになっている（結局使わなかったが）。気に入らなければ後述するルックアップテーブルの値を調整すればよい。
既存のソースコードの (7, 7) のマトリックスの内部表現からこのバイナリ64ビットの内部表現に変換するため、以下のような関数を定義した。本格的なループ処理の前後で使われるだけなので、かなり適当な書き方だが、ここでの実行時間は大したことないので問題ない。
def generate_mat2bin_lut():
lut = np.zeros((7, 7), dtype=np.uint64)
for y in range(7):
for x in range(7):
b = 8 * y + x
lut[y, x] = base_bit >> (63- 8 * y - x)
return lut
def mat2bin(mat):
assert mat.shape[0] == 7 and mat.shape[1] == 7
lut = mat2bin_lut.copy()
lut[mat == 0] = 0
return np.sum(lut) # as np.unit64
def bin2mat(b):
mat = np.zeros((7, 7), dtype=np.uint8)
for y in range(7):
for x in range(7):
if (b & mat2bin_lut[y, x]):
mat[y, x] = 1mat2bin_lut = generate_mat2bin_lut()
二進数について理解している人なら（もっといいやりかたがあるよという指摘はいくらでもあるだろうが）何をしているか判ると思うので、詳細は省略する。
これで mat2bin() によりマトリックス表現からバイナリ表現に変換し、 bin2mat() によりバイナリ表現からマトリックス表現に変換できるので、既存コードにある初期値や出力用の関数はほぼそのまま使い回せる。
バイナリ表現化したことで盤面のほか候補についても64bitで表現できるようになったし、いっそのこと、あらかじめ8つのブロックの方向・反転・位置スライドの全パターンをメモリ上に展開しておくことにした。たとえばひとつのブロックの候補は以下のような196のuint64値で表現できる。
この値を bin2mat でデコードしてみるとブロックが見えてくる。
% python Python 3.8.1 (default, Mar 12 2020, 17:27:26) [Clang 11.0.0 (clang-1100.0.33.17)] on darwin Type "help", "copyright", "credits" or "license" for more information. >>> import puzzle_a_day_solver >>>>>> puzzle_a_day_solver.bin2mat([0x000000000000030e]) array([[0, 1, 1, 1, 0, 0, 0], [1, 1, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0]], dtype=uint8) >>> puzzle_a_day_solver.bin2mat([0x000000000000061c]) array([[0, 0, 1, 1, 1, 0, 0], [0, 1, 1, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0]], dtype=uint8)
8ブロック合計では、重複を除くと合計1196パターンあり、正味9568バイトのデータ量に収まるので、全てのパターンを予め展開しておいても、Core 2 DuoのプロセッサでももしかするとL1キャッシュに収まってしまうかもしれないデータ量だといえる（プロファイリングしていないので実際にどうなっているかは判らない）。
ブロックの衝突チェックは、従来方法では足し算をしていたが、今回は論理積（AND）を使用する。盤面 AND 候補を計算した結果、重複ビットがあるとゼロ以外の値になるので、論理積の結果がゼロ以外の候補は除外する。ここでも NumPy による vectorization を行う。
blocks_f2_masked = blocks_f1_uint64 & mat_uint64
blocks_f2_bool =blocks_f2_masked == 0
blocks_f2_uint64 = blocks_f1_uint64[blocks_f2_bool]特定セルを埋めてくれる候補かどうかを判断するためには、特定セルを示すビットが 1 である候補に絞り込めばよい。 mat2bin/bin2mat 向けに作った、各セルに対するビットを示すルックアップテーブルから 64bit のビットマスクを得て、先と同様にフィルタリングする。
bit_pos = mat2bin_lut[pos[0], pos[1]]
blocks_uint64_masked = blocks_uint64 & bit_pos
blocks_f1_bool =blocks_uint64_masked != 0
blocks_f1_uint64 = blocks_uint64[blocks_f1_bool]基本的なアルゴリズムは変えずに内部表現だけ変更しているので、得られる結果は一緒だが、49秒かかっていた処理が3.5秒まで短縮できた。
単純な計算部分だけで考えればもっと速くなる（1秒は切れるだろう）と思っていたので、思ったより遅かった。
うっかり今晩も最適化して、昨晩から14倍速くなった。私にはこれぐらいが限界かも。 pic.twitter.com/moY3RjJs4Q
— Takeshi HASEGAWA (@hasegaw) May 27, 2021
C言語で再実装してみた
翌週、実はC言語で再実装してみた。
A Puzzle of the Day のソルバーが 10分 → 60秒 → 4秒 → 1.3秒 → 0.04秒 イマココ pic.twitter.com/LhYqLqnxeR
— Takeshi HASEGAWA (@hasegaw) June 5, 2021
もうこれ以上は最適化なんてしないって思っていた。でも1秒を切れなかったのがやはり心残りだったんだと思う。
基本的なアルゴリズムは Python バージョンと一緒だが、本質的ではないところで 楽をしたり、異なる実装をしたりている。
まず Python バージョンでは盤面の制約やブロックの形状を配列として表現しておき、実行時に 64bit のバイナリ表現を得ていた。 C 言語でこれを書き直そうとするとちょっと面倒なので、 Python バージョンが得られたバイナリ表現をソースコード中にLUTとして埋め込みしておく。
ブロックのLUTと、ブロックのメタデータ（LUT内の開始位置と数量）を独立したデータとして持っている。各ブロックの大きさや形状の違いからブロック数量が異なる。このためLUT自体はひとつ一次元配列で持っているので、CPU L1キャッシュを無駄にしない。
現状のコードは日付の指定や結果の出力機能を持っていない。 Python バージョンと同じデータ構造なので、 Python モジュールとしてコンパイルしてしまい、 Python バージョンを宿主として足りない部分は既存コードを使い回そうかなと思っている。
C言語で書くと有り難いのが次に攻略すべきセルを選択する部分で、ここはどうしても Python 上で速くする方法を思いついていなかった。いくつかのパターンを試した感じでは NumPy に任せておいたほうが速かったので、そうなっていた。
# バイナリ化前# バイナリ化後 def next_cell(mat): mat = mat2bin_lut & mat pos = np.unravel_index(np.argmin(mat), mat.shape) minval = mat[pos[0], pos[1]] if minval != 0: return None return pos
この部分はC言語であれば素直に書ける。
/* 盤面 board に対して次のセルを選んで、その bit の値を返す */
inline uint64_t next_cell_bit(uint64_t board) {
for (int y = 0 ; y < 7; y++)
for (int x = 0 ;x < 7; x++) {
uint64_t done = board & mat2bin_lut[y][x];
if (! done)
return mat2bin_lut[y][x];
}
return 0;
}ある時点のC言語バージョンでは gcc の最適化なしで 40ms (0.04秒）、 -O3 で 20ms まで実行時間を短縮できる事を確認できた。M1 Macでは10msで完了したという報告ももらった（速っ！）。現時点のものは、もう少し最適化が進んでいる。
キャッシュミス分析
このソルバーは約10KBのLUTを使って、最大8回の再帰呼び出しするだけなので、この規模であれば、ここ10年のプロセッサーではCPUのL1 Cacheに載りきってしまうことが想像できる。以下は Ryzen 9 3900X の Linux box において cachegrind でこのソルバーをプロファイルしてみた結果である。
$ valgrind --tool=cachegrind ./solver ==114519== Cachegrind, a cache and branch-prediction profiler ==114519== Copyright (C) 2002-2017, and GNU GPL'd, by Nicholas Nethercote et al. ==114519== Using Valgrind-3.15.0 and LibVEX; rerun with -h for copyright info ==114519== Command: ./solver ==114519== --114519-- warning: L3 cache found, using its data for the LL simulation. 50 solutions found ==114519== ==114519== I refs: 100,596,976 .... 命令の参照数 ==114519== I1 misses: 1,140 .... L1 キャッシュミス ==114519== LLi misses: 1,116 .... L3 キャッシュミス ==114519== I1 miss rate: 0.00% ==114519== LLi miss rate: 0.00% ==114519== ==114519== D refs: 19,428,245 (18,233,828 rd + 1,194,417 wr) .... データ読み書き数 ==114519== D1 misses: 3,401 ( 2,718 rd + 683 wr) .... L1 キャッシュミス ==114519== LLd misses: 2,817 ( 2,213 rd + 604 wr) .... L3 キャッシュミス ==114519== D1 miss rate: 0.0% ( 0.0% + 0.1% ) ==114519== LLd miss rate: 0.0% ( 0.0% + 0.1% ) ==114519== ==114519== LL refs: 4,541 ( 3,858 rd + 683 wr) ==114519== LL misses: 3,933 ( 3,329 rd + 604 wr) ==114519== LL miss rate: 0.0% ( 0.0% + 0.1% )
I1 miss rate および D1 miss rate が 0.0% となっていることで、プログラム自体もデータも L1 キャッシュに収まっていることがわかる。念のためソルバーの実行回数を増やしたものを同じくプロファイラにかけてみたりもしたが、大きくキャッシュミスが増えたり、キャッシュミス率が上昇したりする状況を確認できなかった。
LLd miss rate の 2213 rd は恐らく 64bit のキャッシュラインへの読み込みが2213回発生したということだと思うので、 2,213*8=17,704 [bytes] で、コード中で使っている LUT サイズが 10KB 前後、あとはヒープや自分が書いたものではない関数内のメモリ操作をあわせると、こんなものなんだなと思っている。
知識として CPU のキャッシュの事は知っていたけども、こういったデータを具体的に見る機会は今までなかったので、パズルついでにここまで楽しめるのはお得感があった。
これで終わりだと思うけど
- mat2bin_lut の unroll
試した範囲では認識できるほどの差がないことがわかっている - mat2bin_lut を引くのをやめて毎回計算で求める (のとどちらが速い？)
- AVX512
おわりに
とりあえず数日間でソルバーを書いて遊んでみた内容をまとめてみた。
A Puzzle A Day のパズルはアナログでありながら、365日飽きもせずに遊べそうで、とても面白いと思う。同時に盤面は 7x7 未満とかなり小さく、私みたいに総当たりしてもそれほど難しくないので、プログラミング経験が浅くてもソルバーを実装にチャレンジするにはお手頃な問題ではないかと思う。
今日はパターンが多く難易度低めのようなので、持ってる人はパズル楽しんで欲しい pic.twitter.com/ZQ8BtgLcf2
— Takeshi HASEGAWA (@hasegaw) May 27, 2021
