Advent Calendar 3日目の記事 組合せ最適化でボンバーパズルを解く
Advent Calendar 5日目の記事 組合せ最適化でステンドグラスを解く
これなに
エデンの園配置をしっていますか？
エデンの園配置とは、セル・オートマトンにおいて他のいかなる配置からも到達できない配置を指す。
//ja.wikipedia.org/wiki/%E3%82%A8%E3%83%87%E3%83%B3%E3%81%AE%E5%9C%92%E9%85%8D%E7%BD%AE)より
セル・オートマトンの一種である[ライフゲーム](https://ja.wikipedia.org/wiki/%E3%83%A9%E3%82%A4%E3%83%95%E3%82%B2%E3%83%BC%E3%83%A0)の下記のエデンの園配置が本当にそうか確かめてみましょう。
- 参考：[Achim's Game of Life Page](http://wwwhomes.uni-bielefeld.de/achim/orphan_9th.html)
組合せ最適化で使われるソルバーは、全ての組合せを調べてくれます。ソルバーで答えが出ないことがわかれば、エデンの園配置であることが証明できます。
## 考え方
ライフゲームで許される条件を数式で表現します。
- 現在＝生：
- 過去＝生、過去の周り8マス＝2
- 過去の周り8マス＝3
- 現在＝死：
- 過去の周り8マス<=1
- 過去＝死、過去の周り8マス＝2
- 過去の周り8マス>＝4
OR条件は、線形式で表現できないので、0-1変数を使って場合分けします。
- 参考：[パズルでみる組合せ最適化のテクニック - IF条件](https://qiita.com/SaitoTsutomu/items/f05f4ff05d4ce6099ff3#if%E5%88%B6%E7%B4%84)
## Pythonのプログラム
1つ前の状態を計算するプログラム(`solve`)は下記のようになります。
```py3: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つ前の状態を計算できているか、簡単な例で確かめてみましょう。
（下記は、エデンの園配置ではなく[グライダー](https://ja.wikipedia.org/wiki/%E3%82%B0%E3%83%A9%E3%82%A4%E3%83%80%E3%83%BC_(%E3%83%A9%E3%82%A4%E3%83%95%E3%82%B2%E3%83%BC%E3%83%A0))です）
```py3:python
solve("""\
.#.
#..
###""".splitlines())
>>>
Optimal
```
`Optimal`なので、解が存在します。確かに、図の状態から1ステップ進めると、入力の状態になることが確認できます。`solve`が1つ前の状態（の1つ）を計算できているのが確認できました。
## 確認その2
```py3:python
solve("""\
.##..#.##...
#..##..#.#.#
.#.#.##.#.#.
#....##..##.
.###...#....
..#.#.##.#..
.#.##...#.#.
#....#.#....""".splitlines())
>>>
Infeasible
```
今度は、`Infeasible`[^1]です。確かに、エデンの園配置のようです。
[^1]: 2018/10/25、筆者のPRでバグが修正されて「Infeasible」と表示されるようになりました。 [fix: issue-171, add MIP Infeasible status](https://github.com/coin-or/pulp/pull/172#event-1924585344)
以上