2021年5月28日金曜日

A Puzzle A Day のソルバーを書いてみた(6/20更新)

更新内容

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 未満とかなり小さく、私みたいに総当たりしてもそれほど難しくないので、プログラミング経験が浅くてもソルバーを実装にチャレンジするにはお手頃な問題ではないかと思う。


1 件のコメント:

  1. The moulding could free fall into a collection field or onto a switch conveyer, or may be be} eliminated by an automated robotic. In semi-automatic mode, the operator could intervene at this point within the cycle to take away the moulding manually. Once the moulding is clear from the mould software, the entire moulding cycle could be repeated. The second kind of general clamping arrangement is known as the Toggle Lock. Lastly, let’s look at Comfortable Underwear for Women at|have a glance at} how you can determine a plastic injection molder that's right in your project. There are many elements which are be} considered when determining the dimensions of the press.

    返信削除