更新内容
6/20 ... C言語版について追記しました。
面白そうなパズルとの出会い
Twitter で A Puzzle A Day というカレンダーパズルを見かけた。
このパズルは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秒以内で完了するようになった 。
バイナリ化
ベクタライゼーションにより特定の日付向けの解法が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] = 1
mat2bin_lut = generate_mat2bin_lut()
二進数について理解している人なら(もっといいやりかたがあるよという指摘はいくらでもあるだろうが)何をしているか判ると思うので、詳細は省略する。
これで mat2bin() によりマトリックス表現からバイナリ表現に変換し、 bin2mat() によりバイナリ表現からマトリックス表現に変換できるので、既存コードにある初期値や出力用の関数はほぼそのまま使い回せる。
バイナリ表現化したことで盤面のほか候補についても64bitで表現できるようになったし、いっそのこと、あらかじめ8つのブロックの方向・反転・位置スライドの全パターンをメモリ上に展開しておくことにした。たとえばひとつのブロックの候補は以下のような196のuint64値で表現できる。
000000000000030e
000000000000061c
000000000000070c
0000000000000c07
0000000000000c38
0000000000000e03
0000000000000e18
000000000000180e
0000000000001870
0000000000001c06
0000000000001c30
000000000000301c
000000000000380c
0000000000003860
0000000000006038
0000000000007018
0000000000030e00
0000000000061c00
0000000000070c00
00000000000c0700
00000000000c3800
00000000000e0300
00000000000e1800
0000000000180e00
0000000000187000
00000000001c0600
00000000001c3000
0000000000301c00
0000000000380c00
0000000000386000
0000000000603800
0000000000701800
0000000001010302
0000000001030202
0000000002020301
0000000002020604
0000000002030101
0000000002060404
00000000030e0000
0000000004040602
0000000004040c08
0000000004060202
00000000040c0808
00000000061c0000
00000000070c0000
0000000008080c04
0000000008081810
00000000080c0404
0000000008181010
000000000c070000
000000000c380000
000000000e030000
000000000e180000
0000000010101808
0000000010103020
0000000010180808
0000000010302020
00000000180e0000
0000000018700000
000000001c060000
000000001c300000
0000000020203010
0000000020206040
0000000020301010
0000000020604040
00000000301c0000
00000000380c0000
0000000038600000
0000000040406020
0000000040602020
0000000060380000
0000000070180000
0000000101030200
0000000103020200
0000000202030100
0000000202060400
0000000203010100
0000000206040400
000000030e000000
0000000404060200
00000004040c0800
0000000406020200
000000040c080800
000000061c000000
000000070c000000
00000008080c0400
0000000808181000
000000080c040400
0000000818101000
0000000c07000000
0000000c38000000
0000000e03000000
0000000e18000000
0000001010180800
0000001010302000
0000001018080800
0000001030202000
000000180e000000
0000001870000000
0000001c06000000
0000001c30000000
0000002020301000
0000002020604000
0000002030101000
0000002060404000
000000301c000000
000000380c000000
0000003860000000
0000004040602000
0000004060202000
0000006038000000
0000007018000000
0000010103020000
0000010302020000
0000020203010000
0000020206040000
0000020301010000
0000020604040000
0000030e00000000
0000040406020000
000004040c080000
0000040602020000
0000040c08080000
0000061c00000000
0000070c00000000
000008080c040000
0000080818100000
0000080c04040000
0000081810100000
00000c0700000000
00000c3800000000
00000e0300000000
00000e1800000000
0000101018080000
0000101030200000
0000101808080000
0000103020200000
0000180e00000000
0000187000000000
00001c0600000000
00001c3000000000
0000202030100000
0000202060400000
0000203010100000
0000206040400000
0000301c00000000
0000380c00000000
0000386000000000
0000404060200000
0000406020200000
0000603800000000
0000701800000000
0001010302000000
0001030202000000
0002020301000000
0002020604000000
0002030101000000
0002060404000000
00030e0000000000
0004040602000000
0004040c08000000
0004060202000000
00040c0808000000
00061c0000000000
00070c0000000000
0008080c04000000
0008081810000000
00080c0404000000
0008181010000000
000c070000000000
000c380000000000
000e030000000000
000e180000000000
0010101808000000
0010103020000000
0010180808000000
0010302020000000
00180e0000000000
0018700000000000
001c060000000000
001c300000000000
0020203010000000
0020206040000000
0020301010000000
0020604040000000
00301c0000000000
00380c0000000000
0038600000000000
0040406020000000
0040602020000000
0060380000000000
0070180000000000
この値を 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秒は切れるだろう)と思っていたので、思ったより遅かった。
C言語で再実装してみた
翌週、実はC言語で再実装してみた。
もうこれ以上は最適化なんてしないって思っていた。でも1秒を切れなかったのがやはり心残りだったんだと思う。
基本的なアルゴリズムは Python バージョンと一緒だが、本質的ではないところで 楽をしたり、異なる実装をしたりている。
まず Python バージョンでは盤面の制約やブロックの形状を配列として表現しておき、実行時に 64bit のバイナリ表現を得ていた。 C 言語でこれを書き直そうとするとちょっと面倒なので、 Python バージョンが得られたバイナリ表現をソースコード中にLUTとして埋め込みしておく。
ブロックのLUTと、ブロックのメタデータ(LUT内の開始位置と数量)を独立したデータとして持っている。各ブロックの大きさや形状の違いからブロック数量が異なる。このためLUT自体はひとつ一次元配列で持っているので、CPU L1キャッシュを無駄にしない。
現状のコードは日付の指定や結果の出力機能を持っていない。 Python バージョンと同じデータ構造なので、 Python モジュールとしてコンパイルしてしまい、 Python バージョンを宿主として足りない部分は既存コードを使い回そうかなと思っている。
C言語で書くと有り難いのが次に攻略すべきセルを選択する部分で、ここはどうしても Python 上で速くする方法を思いついていなかった。いくつかのパターンを試した感じでは NumPy に任せておいたほうが速かったので、そうなっていた。
# バイナリ化前
def next_cell(context):
mat = context['mat']
pos = np.unravel_index(np.argmin(mat), mat.shape)
minval = mat[pos[0], pos[1]]
if minval != 0:
return None
return pos
# バイナリ化後
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 未満とかなり小さく、私みたいに総当たりしてもそれほど難しくないので、プログラミング経験が浅くてもソルバーを実装にチャレンジするにはお手頃な問題ではないかと思う。