以前、マインスイーパーの初手について、記事（マインスイーパーの初手 - pzdc）を書いた。
この問題について、プログラムを書いて調べてみる。今回は、カウンタの実装までについて、書く。
計算したい式
確率を計算する必要があるため、その確率の分子と分母を知る必要がある。そのためには、ある特定の盤面状態における解の数を知る必要がある。解の数さえわかれば、あとは木を探索するだけで良い。
もっと具体的に書いておく。
マインスイーパーの盤面状態Bに対して、ある着手cを行った時に、n（n=0～8）が出現する確率を、p(n|B, c)とする。また、nが出現したことで得られる新しい盤面状態をBc, nとする。そして、一般に盤面状態Bの（最適戦略を取った場合の）成功率をs(B)とする。すると、
s(B) = maxc {Σn p(n|B, c) s(Bc, n)}
「失敗する確率」については考えなくて良い。失敗するケースの「成功率」への確率寄与はゼロであるためである。
nが出現する確率p(n|B, c)は、次のように、解の数の比率で求まる。
p(n|B, c) = #(Bn, c) / #(B)
ここで、「#(B)」は、Bの解の数を表す。
この解の数を計算する部分が最初の目標地点である。今回の記事では、カウント処理を実装して、現実的な時間で解の数を計算できるようにする。
解の数の決定処理の実装
API
今回は、ソルバーを使わずに、Python（numpyを使用する）で処理を書く。まず、解の数の決定処理を、次のように整理する。
問題設定:
N個のブール変数(x0, x1, ... x_{N-1})がある
制約が、M個与えられる
それぞれの制約は、変数集合の中から任意個数の変数を指定して、
その中に含まれるTrueの個数を与える
問題の表現:
変数選択行列X: MxNの行列. np.ndarray, np.int32
値指定行列y: 長さMのベクトル. np.ndarray, np.int32
解の数を数えるCounterクラスを定義しよう。APIは、次のようにする。
class Counter: def count(self, X, y) -> int: ...
ベースライン実装
まずは、難しいことは気にせず、実装してみる。ここで作る実装をベースラインにして、少しずつ高速化などの調整をしていく。
class BaselineCounter: __slots__ = () def count(self, X, y): if len(X[0]) == 0: if all(y == 0): return 1 else: return 0 if any(y < 0): return 0 # 全ての変数の係数が0なのに値がノンゼロ、という制約があれば、解なし coefficient_zero_mask = (np.sum(X, axis=1) == 0) if any((y != 0) & coefficient_zero_mask): return 0 # 最初の変数x0がFalseであると仮定した場合の、残りの問題の表現 X_0F = X[:, 1:].copy() y_0F = y.copy() # 部分問題の解の数 num_solution_0F = self.count(X_0F, y_0F) # 最初の変数x0がTrueであると仮定した場合の、残りの問題の表現 X_0T = X[:, 1:].copy() y_0T = y.copy() - X[:, 0] # 部分問題の解の数 num_solution_0T = self.count(X_0T, y_0T) return num_solution_0F + num_solution_0T
愚直に一変数ずつ見て、FalseのケースとTrueのケースの両方について再帰呼び出しで求めた解の数の和を返している。
動作確認
ベースライン実装が期待通りに動作することを確認する。そのために、（例えば）次のようなテストを実行して、計算結果が期待通りになっていることを確認した（テストにはpytestを使用した）。
class TestCounter: """カウンターに対するテスト.""" def test_baseline_counter(self): # arrange X = np.array([ [1, 1, 1, 1, 0, 0, 0, 0], [0, 0, 0, 1, 1, 1, 1, 0], [1, 1, 1, 1, 1, 1, 1, 1], ]) y = np.array([2, 2, 5]) # act n = BaselineCounter().count(X, y) # assert assert n == 9
処理時間計測：タイマー
処理時間の計測に、次のタイマーを使う。
class Timer: """処理時間計測.""" __slots__ = ('reporter', 'start_time', 'spent_time') def __init__(self, reporter=None): self.reporter = print if reporter is None else reporter def __enter__(self): self.start_time = time.perf_counter() def __exit__(self, *a, **k): self.spent_time = time.perf_counter() - self.start_time self.reporter(self.spent_time)
処理時間計測：ベンチマーク作成
カウンタークラスを改善していくうえで、処理時間のベンチマークにする問題を用意しておきたい、と考えた。
下に載せるクラスProblemGeneratorは、メソッドgenerateの引数として、n_problems（生成する問題(X, y)の総数）、M_range（制約の数の生成範囲）、N_range（変数の数の生成範囲）、rng（numpyの乱数ジェネレータ）を受け取り、一連の(X, y)を生成するジェネレータを返す。
class ProblemGenerator: """処理時間ベンチマーク用の問題生成器. """ def generate(self, *a, **k): self._check_args(*a, **k) yield from self._generate_simple(*a, **k) def _check_args(self, n_problems, M_range, N_range, rng): assert n_problems > 0 assert len(M_range) == 2 assert M_range[0] > 0 assert M_range[1] > M_range[0] assert len(N_range) == 2 assert N_range[0] > 0 assert N_range[1] > N_range[0] def _generate_simple(self, n_problems, M_range, N_range, rng): for i_problem in range(n_problems): M = rng.integers(M_range[0], M_range[1]) N = rng.integers(N_range[0], N_range[1]) X = rng.integers(0, 2, (M, N), dtype=np.int32) y = rng.integers(0, N + 1, M, dtype=np.int32) yield X, y
非常に単純な実装だ。とりあえず使用はできそうだが、実装している最中に、「これで良いのか？」と疑問に感じる箇所が、2か所あった。Xは、0か1の値をとる行列であり、上の実装では行列の各成分について0と1を等確率で選んでいる。これは、マインスイーパーの問題設定とは異なっており、今回の用途には適さない可能性がある。yについては、0以上N以下の整数を等確率で選んでいるが、これも少々おかしい。N個の変数をランダムに選択した場合の、ターゲットの発見数の分布は、一様分布ではないだろう。
やはり、気になるので、もう少しマインスイーパーを意識した生成方法に変更する。
def _generate_like_minesweeper(self, n_problems, M_range, N_range, rng): """マインスイーパーの問題設定に近い条件で生成する. """ for i_problem in range(n_problems): M = rng.integers(M_range[0], M_range[1]) N = rng.integers(N_range[0], N_range[1]) if (N > 8) and (M > 1): X = np.zeros((M, N), dtype=np.int32) y = np.zeros(M, dtype=np.int32) for m in range(M - 1): target_count = rng.integers(1, 9) choice = rng.choice(np.arange(N), target_count) X[m, choice] = 1 y[m] = rng.integers(0, target_count + 1) X[M-1, :] = 1 y[M-1] = rng.integers(np.max(y), N + 1) else: X = rng.integers(0, 2, (M, N), dtype=np.int32) y = rng.integers(0, N + 1, M, dtype=np.int32) yield X, y
これも、それほど精密ではないが、先ほどのものよりは現実的な問題に近くなったのではないかと思う。
評価結果
次のようにして評価する。
rng = np.random.default_rng(42) for X, y in generator.generate(10, (1, 6), (5, 21), rng): ...
評価結果は次の通り。出力カラムは順番に、盤面サイズ、解の数、使用したカウンタクラス名、処理時間。
1x17 (#=10560) BaselineCounter: 1.85565680 [sec] 3x10 (#=0) BaselineCounter: 0.00072850 [sec] 3x19 (#=9009) BaselineCounter: 1.10790760 [sec] 1x20 (#=0) BaselineCounter: 5.47811750 [sec] 2x6 (#=0) BaselineCounter: 0.00078840 [sec] 5x15 (#=18) BaselineCounter: 0.15075170 [sec] 4x12 (#=0) BaselineCounter: 0.00213060 [sec] 3x13 (#=0) BaselineCounter: 0.02239800 [sec] 2x9 (#=1) BaselineCounter: 0.00606240 [sec] 5x9 (#=0) BaselineCounter: 0.00004180 [sec]
今回は詳細な処理時間の集計調査はしないで、ざっと目視で処理時間を見るだけにする。BaselineCounterは改善の余地が多すぎるくらいで、細かい統計はあまり意味がないと考えているためである。
上の評価では制約の数Mを(1, 6)、変数の数Nを(5, 21)の範囲で生成して評価した。実は、BaselineCounterは処理時間が長いため、問題のサイズを小さく設定している。本当は、もう少し大きい問題が解きたい。例えば、Mを(2, 11)、Nを(20, 51)の範囲で10問生成するようにすると、
待機...(200秒以上)
しばらくPCのファンが唸っているが、結果が出力される気配がない。BaselineCounterは明らかに処理に時間がかかりすぎている。ここでPCのファンを唸らせている最初の一問は、M=2, N=43の問題で、解の数は 15471286560 である。
X=[[0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0] [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]] y=[ 0 23]
後でこの問題が0.0005秒程度で解けるようになるのを見る。
無駄な計算のカット
まず、解の数が明らかにゼロの場合は、すぐに計算を終了することができる。そこで、次の処理を追加してみる。
if any(np.sum(X, axis=1) < y): return 0
もう一度、最初のベンチマークを回してみる。最初のカウンタ実装(BaselineCounter)の結果と今回の実装(CutCounter)の結果を交互に並べる。
1x17 (#=10560) BaselineCounter: 1.95172680 [sec] 1x17 (#=10560) CutCounter: 0.56069510 [sec] 3x10 (#=0) BaselineCounter: 0.00096370 [sec] 3x10 (#=0) CutCounter: 0.00096790 [sec] 3x19 (#=9009) BaselineCounter: 1.10942420 [sec] 3x19 (#=9009) CutCounter: 1.25317910 [sec] 1x20 (#=0) BaselineCounter: 5.86163090 [sec] 1x20 (#=0) CutCounter: 0.00002380 [sec] 2x6 (#=0) BaselineCounter: 0.00073960 [sec] 2x6 (#=0) CutCounter: 0.00001810 [sec] 5x15 (#=18) BaselineCounter: 0.14047630 [sec] 5x15 (#=18) CutCounter: 0.00419070 [sec] 4x12 (#=0) BaselineCounter: 0.00208820 [sec] 4x12 (#=0) CutCounter: 0.00051980 [sec] 3x13 (#=0) BaselineCounter: 0.02039230 [sec] 3x13 (#=0) CutCounter: 0.00019950 [sec] 2x9 (#=1) BaselineCounter: 0.00873410 [sec] 2x9 (#=1) CutCounter: 0.00103560 [sec] 5x9 (#=0) BaselineCounter: 0.00004610 [sec] 5x9 (#=0) CutCounter: 0.00002580 [sec]
全体的に改善している。
Xの事前ソート
Xは、列を並び替えても問題は本質的に変化しない。ハタンは木の上の方で判定できた方が嬉しいので、制約がなるべく最初の方に出現するように並び替えることにする。
これまでのcount(X, y)の実装を_countにリネームして、count(X, y)が呼ばれたときはXをソートしてから_count(X, y)を呼ぶようにする。
def count(self, X, y): M = len(X) X_sort = X.copy() # M-1行目から0行目までを順番にキーとして、降順ソート for m in range(M-2, -1, -1): order = np.argsort( -X_sort[m, :], axis=0, kind='stable', ) X_sort = X_sort[:, order] return self._count(X_sort, y)
評価する。
1x17 (#=10560) CutCounter: 0.56773230 [sec] 1x17 (#=10560) SortCounter: 0.60588010 [sec] 3x10 (#=0) CutCounter: 0.00095840 [sec] 3x10 (#=0) SortCounter: 0.00026760 [sec] 3x19 (#=9009) CutCounter: 1.28898780 [sec] 3x19 (#=9009) SortCounter: 0.55118530 [sec] 1x20 (#=0) CutCounter: 0.00001910 [sec] 1x20 (#=0) SortCounter: 0.00002980 [sec] 2x6 (#=0) CutCounter: 0.00001780 [sec] 2x6 (#=0) SortCounter: 0.00004680 [sec] 5x15 (#=18) CutCounter: 0.00425780 [sec] 5x15 (#=18) SortCounter: 0.00173800 [sec] 4x12 (#=0) CutCounter: 0.00048140 [sec] 4x12 (#=0) SortCounter: 0.00020580 [sec] 3x13 (#=0) CutCounter: 0.00014470 [sec] 3x13 (#=0) SortCounter: 0.00009480 [sec] 2x9 (#=1) CutCounter: 0.00099620 [sec] 2x9 (#=1) SortCounter: 0.00023290 [sec] 5x9 (#=0) CutCounter: 0.00001680 [sec] 5x9 (#=0) SortCounter: 0.00005490 [sec]
思ったよりは効果がなかった。しかし、比較的長い時間がかかっていた3x19の問題で処理時間の改善が見られる。
不要になった制約行を削除
行列Xから順番に変数を消していく際に、制約は不要になっていく。不要になった制約行を消すことで問題のサイズが小さくなるので、多少高速化することを期待した。
# 不要になった行を削除 shrink_mask = ~(np.sum(X, axis=1) == 0) if all(~shrink_mask): # 変数が残っているのに情報が一切ないなら、残りの変数はすべて自由 return np.power(2, len(X[0])) X_shrink = X[shrink_mask] y_shrink = y[shrink_mask]
結果の処理時間を、ひとつ前のCutCounterと比べる。
1x17 (#=10560) CutCounter: 0.55367210 [sec] 1x17 (#=10560) ShrinkCounter: 0.81817200 [sec] 3x10 (#=0) CutCounter: 0.00096930 [sec] 3x10 (#=0) ShrinkCounter: 0.00034090 [sec] 3x19 (#=9009) CutCounter: 1.15895720 [sec] 3x19 (#=9009) ShrinkCounter: 0.74828370 [sec] 1x20 (#=0) CutCounter: 0.00001740 [sec] 1x20 (#=0) ShrinkCounter: 0.00002800 [sec] 2x6 (#=0) CutCounter: 0.00001710 [sec] 2x6 (#=0) ShrinkCounter: 0.00004280 [sec] 5x15 (#=18) CutCounter: 0.00428470 [sec] 5x15 (#=18) ShrinkCounter: 0.00222600 [sec] 4x12 (#=0) CutCounter: 0.00046730 [sec] 4x12 (#=0) ShrinkCounter: 0.00015180 [sec] 3x13 (#=0) CutCounter: 0.00012890 [sec] 3x13 (#=0) ShrinkCounter: 0.00013460 [sec] 2x9 (#=1) CutCounter: 0.00096070 [sec] 2x9 (#=1) ShrinkCounter: 0.00031570 [sec] 5x9 (#=0) CutCounter: 0.00001550 [sec] 5x9 (#=0) ShrinkCounter: 0.00007650 [sec]
今回も、大きな改善は見られない。
最後の制約行の高速化
Xが1行になれば、もはや条件分岐で1列ずつ削る計算をする必要はない。制約の与えられている列の変数については、コンビネーション計算でターゲットの全配置が決まる。制約の与えられていない列については、ターゲットを入れても入れなくても良いため、2のべき乗通りの配置がありうる。これらをかけ合わせればよい。したがって、次のように計算できる。
# 最後の行だけは工夫 if len(X) == 1: X_one_count = np.sum(X) X_zero_count = len(X[0]) - X_one_count y_value = y[0] # X_one_count個数の中からy_value個数を選ぶ場合の数 select_count = math.comb(X_one_count, y_value) rest_count = np.power(2, X_zero_count) return select_count * rest_count
評価する。
1x17 (#=10560) ShrinkCounter: 0.86067220 [sec] 1x17 (#=10560) LastSpeedupCounter: 0.00007200 [sec] 3x10 (#=0) ShrinkCounter: 0.00033490 [sec] 3x10 (#=0) LastSpeedupCounter: 0.00032730 [sec] 3x19 (#=9009) ShrinkCounter: 0.78144230 [sec] 3x19 (#=9009) LastSpeedupCounter: 0.00052140 [sec] 1x20 (#=0) ShrinkCounter: 0.00002110 [sec] 1x20 (#=0) LastSpeedupCounter: 0.00002820 [sec] 2x6 (#=0) ShrinkCounter: 0.00003390 [sec] 2x6 (#=0) LastSpeedupCounter: 0.00004110 [sec] 5x15 (#=18) ShrinkCounter: 0.00257470 [sec] 5x15 (#=18) LastSpeedupCounter: 0.00073990 [sec] 4x12 (#=0) ShrinkCounter: 0.00013630 [sec] 4x12 (#=0) LastSpeedupCounter: 0.00014680 [sec] 3x13 (#=0) ShrinkCounter: 0.00011880 [sec] 3x13 (#=0) LastSpeedupCounter: 0.00012110 [sec] 2x9 (#=1) ShrinkCounter: 0.00031030 [sec] 2x9 (#=1) LastSpeedupCounter: 0.00024280 [sec] 5x9 (#=0) ShrinkCounter: 0.00004680 [sec] 5x9 (#=0) LastSpeedupCounter: 0.00005570 [sec]
あまり改善していないように見える。
しかし、もっと大きな問題サイズ（Mを(2, 11)、Nを(20, 51)）に対して評価してみると、「計算できる」ことがわかる。
2x43 (#=15471286560) LastSpeedupCounter: 0.00042050 [sec] 10x42 (#=0) LastSpeedupCounter: 0.00034120 [sec] 5x30 (#=0) LastSpeedupCounter: 0.00136400 [sec] 5x47 (#=1113024) LastSpeedupCounter: 0.00184530 [sec] 7x34 (#=0) LastSpeedupCounter: 0.00009870 [sec] 8x34 (#=0) LastSpeedupCounter: 0.00295910 [sec] 6x46 (#=0) LastSpeedupCounter: 0.00033120 [sec] 3x29 (#=74256) LastSpeedupCounter: 0.00393230 [sec] 5x49 (#=630591) LastSpeedupCounter: 0.00538700 [sec] 4x37 (#=65780) LastSpeedupCounter: 0.00080640 [sec]
このサイズで計算できるカウンター実装はこれが初めてなので、明らかに改善している。例えば、最初の2x43の問題（悪魔的に長い処理時間を要していた問題）は、0.0005秒程度で計算できている。
これで、ひとまずは解の数を実用的な時間で計算できるようになった。
いったんここまで
記事が長くなるので、いったんここで一区切りしておく。
次回の記事では、この実装を使ったマインスイーパーの成功確率の計算について、書く予定である。
