Advent Calendar 3日目の記事組合せ最適化でボンバーパズルを解く
Advent Calendar 5日目の記事組合せ最適化でステンドグラスを解く

これなに

エデンの園配置をしっていますか？

エデンの園配置とは、セル・オートマトンにおいて他のいかなる配置からも到達できない配置を指す。

エデンの園配置 - wikipediaより

セル・オートマトンの一種であるライフゲームの下記のエデンの園配置が本当にそうか確かめてみましょう。

参考：Achim's Game of Life Page

組合せ最適化で使われるソルバーは、全ての組合せを調べてくれます。ソルバーで答えが出ないことがわかれば、エデンの園配置であることが証明できます。

考え方

ライフゲームで許される条件を数式で表現します。

現在＝生：
- 過去＝生、過去の周り8マス＝2
- 過去の周り8マス＝3
現在＝死：
- 過去の周り8マス<=1
- 過去＝死、過去の周り8マス＝2
- 過去の周り8マス>＝4

OR条件は、線形式で表現できないので、0-1変数を使って場合分けします。

参考：パズルでみる組合せ最適化のテクニック - IF条件

Pythonのプログラム

1つ前の状態を計算するプログラム(solve)は下記のようになります。

python

import pandas as pd, matplotlib.pyplot as plt
from pulp import LpProblem, LpStatus, lpSum, value
from ortoolpy import addbinvar, addbinvars
def solve(data):
    ni,nj = len(data),len(data[0])
    a = pd.DataFrame([(i,j,data[i][j]!='.') for i in range(ni)
        for j in range(nj)], columns=list('行列値'))
    a['Var'] = addbinvars(len(a))
    m = LpProblem()
    for _,r in a.iterrows():
        v = lpSum(a.query(f'{r.行-1}<=行<={r.行+1}&{r.列-1}<=列<={r.列+1}').Var)-r.Var
        if r.値: # 3 <= v+x, 2v+x <= 7
            m += v + r.Var >= 3
            m += 2*v + r.Var <= 7
        else: # v+x <= 2 or >=4
            y = addbinvar()
            m += v + r.Var <= 2 + 7*y # y==0 → v+x <= 2
            m += v         >= 4*y     # y==1 →   v >= 4
    m.solve()
    print(LpStatus[m.status])
    if m.status==1:
        a['Val'] = a.Var.apply(value)
        plt.imshow((a.Val<0.5).values.reshape(ni,nj), cmap='gray', interpolation='none')
        plt.show()

確認その1

最初に、solveが1つ前の状態を計算できているか、簡単な例で確かめてみましょう。
（下記は、エデンの園配置ではなくグライダーです）

python

solve("""\
.#.
# ..
### """.splitlines())
>>>
Optimal

Optimalなので、解が存在します。確かに、図の状態から1ステップ進めると、入力の状態になることが確認できます。solveが1つ前の状態（の1つ）を計算できているのが確認できました。

確認その2

python

solve("""\
.##..#.##...
# ..##..#.#.#
.#.#.##.#.#.
# ....##..##.
.###...#....
..#.#.##.#..
.#.##...#.#.
# ....#.#....""".splitlines())
>>>
Infeasible

今度は、Infeasible¹です。確かに、エデンの園配置のようです。

2018/10/25、筆者のPRでバグが修正されて「Infeasible」と表示されるようになりました。 fix: issue-171, add MIP Infeasible status
以上 ↩