パズルでみる組合せ最適化のテクニック

追記：参考リンク

- [数理最適化によるパズルの解法 - Qiita](https://qiita.com/SaitoTsutomu/items/0c0db8d22979fc9de8f4)
- [Jupyter Notebook- GitHub](https://github.com/SaitoTsutomu/opt4puzzle)
- [パズル Advent Calendar 2017](https://qiita.com/advent-calendar/2017/puzzle)
- [PyCon JP 2014](https://pycon.jp/2014/schedule/presentation/23/)
  - [動画](https://www.youtube.com/watch?v=DCrDBmEWAzQ)
- [gist](http://nbviewer.jupyter.org/gist/Tsutomu-KKE/1b526ed7ef2ab111c829)
- [Docker](https://hub.docker.com/r/tsutomu7/puzzle/)

# はじめに
[組合せ最適化](bfbf4c185ed7004b5721)(いわゆる混合整数最適化)では、さまざまなルールをシンプルに記述することができます。
ここでは、パズルを例題に、いくつかのテクニックを簡単にご紹介します。
これらのテクニックのいくつかは、Pythonによる部分もあります。数理最適化のモデル化においてPythonは相性がよいといえるでしょう。
## 環境構築
- サクッと試したい場合：Binderというサービスを使うことにより、ブラウザだけで試せます。詳しくは、[無料JupyterサービスのBinderの紹介](821ebd608c412e57382d)をご覧ください。
- きちんと試したい場合：[Anaconda](https://www.continuum.io/downloads)を導入後、下記を実行してください。

```
pip install pulp ortoolpy unionfind more_itertools
```

## 準備
モデラーとしては、PuLPを用います。([PuLPについて](070ca9cb37c6b2b492f0))
後述のコード例では、下記を省略します。

```py3:python
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt
from itertools import groupby
from pulp import *
from unionfind import unionfind
from ortoolpy import addvar, addvars, addbinvar, addbinvars
m = LpProblem()
```

メソッド|説明
:--|:--
groupby|キーが同じものをグループ化する
unionfind|[素集合データ構造](https://ja.wikipedia.org/wiki/%E7%B4%A0%E9%9B%86%E5%90%88%E3%83%87%E3%83%BC%E3%82%BF%E6%A7%8B%E9%80%A0)用のクラス
addvar|変数を1つ作成
addvars|変数のリストを作成
addbinvar|0-1変数を1つ作成
addbinvars|0-1変数のリストを作成
LpProblem|PuLPの数理最適化モデル

## 対象パズル
[SaitoTsutomu/OptForPuzzle](https://github.com/SaitoTsutomu/OptForPuzzle)に用意しました。

# テクニック

## 変数をnp.arrayで作成
効用：演算の高速化、多彩なアクセス
対象パズル：「数独」など
例：

```py3:python
v = np.array(addbinvars(9, 9, 9)) # (1)
m += lpSum(v[i,:,k]) == 1 # (2)
m += lpSum(v[i:i+3,j:j+3,k]) == 1 # (3)
w = addvars(9, 9) # (4)
m += lpDot(range(9), v[i,j])+1 == r[i][j] # (4)
```

- (1)では、9x9x9の多次元配列の0-1変数を作成しています。それぞれの次元は、行、列、数を表しています。
- (2)では、$i$行目に数字$k+1$が1つだけを意味しています。
- (3)では、左上が$(i,j)$となる3x3の領域で、数字$k+1$が1つだけを意味しています。
- また、この$v$を使うと、(4)のように、マス$(i,j)$の数字を変数$w_{ij}$のように表現できます。

## 結果をnp.vectorizeで取得
効用：演算の高速化、多彩なアクセス
対象パズル：「数独」など(可視化は「黒どこ」など)
例：

```py3:python
r = np.vectorize(value)(v).astype(int) # (1)
plt.imshow(1-r, cmap='gray', interpolation='none') # (2)
```

変数をnp.arrayで作成し最適化を解いた結果は、上記の(1)のように結果$r$を取得すると便利です。こうすると、高速に簡単に結果を得ることができ、さらにnumpyの多彩な機能を使って処理を続けられます。
2次元の白黒の結果を図として確認する場合は、(2)のようにmatplotlibを使うと簡単に可視化できます。

なお、変数をpandasのDataFrameのSeriesで管理しているときは、applyで同様にできます。

## すべて同じ
効用：効率的なモデル化
対象パズル：「ペイントエリア」など
例：

```py3:python
for vi, vj in zip(v, v[1:]):
    m += vi == vj
```

変数の1次元配列$v$のすべての要素が同じ値を取らないといけない場合、上記のように記述できます。
(変数そのものを置き換える方法もあります)

## 周りの数
効用：効率的なモデル化
対象パズル：「クリーク」「シャカシャカ」「数コロ」「のりのり」「ペイントエリア」「ボンバーパズル」

マスの変数$v[i,j]$の上下左右の変数の和を使いたいとき、下記のように$w$を使うことによりができます。
例：

```py3:python
u = np.zeros((nh+2, nw+2), dtype=object)
v = u[1:-1,1:-1] = np.array(addbinvars(nh, nw))
w = u[:-2,1:-1]+u[2:,1:-1]+u[1:-1,:-2]+u[1:-1,2:]
```

これは、vの周りに1周多い配列$u$を用意し、うまくスライスを使った例です。


## IF制約
効用：条件によって成り立つ場合を表現
対象パズル：「カナオレ」など
変数$x,y$があったとき、「$x == 1$」ならば「$y \le a$」としたいときは、十分に大きな数$M$を使用して、下記のように書くことができます。
例：

```py3:python
m += y <= a + M(1-x)
```

また、「$x==1$」ならば「$y=a$」としたいときは、下記のように書くことができます。
例：

```py3:python
m += y <= a + M(1-x)
m += y >= a - M(1-x)
```

## 連結制約
効用：簡単なモデル化
対象パズル：「黒どこ」など

「黒どこ」のように「全ての白マスがつながっていること」という制約があります。このような制約を数式で表そうとすると、結構大変だったりします。そこで、一旦無視して、接続していなかったら、その解を禁止して解きなおすやり方があります。

例：連結性の確認

```py3:python
from unionfind import unionfind
r = np.array(
    [[0,0,0],
     [1,1,1],
     [0,0,0]])
print(unionfind.isconnected(1-r))
r[1][1] = 0
print(unionfind.isconnected(1-r))
>>>
False
True
```

結果の$r$は、0:白、1:黒を表しているとします。unionfind.isconnected(1-r) で白が連結しているかを確認できます。
これを踏まえると、下記のようにして連結になるまで、モデルを変更していくことができます。

```py3:python
while True:
    m.solve()
    r = np.vectorize(value)(v).astype(int)
    if unionfind.isconnected(1-r):
        break
    m += lpSum(v[r==1]) <= r.sum() - 1
```

上記のコードはシンプルですが、パズルによっては繰り返しが長くなることがあります。その場合は、「結果の全黒マスを禁止」するのではなく、黒マスで分断された白のエリアごとに「その白を囲っている黒マスを禁止」することにより、効率よく解くことができます。

## エリア(部屋や国)ごとの処理
効用：シンプルな入力とプログラム
対象パズル：「因子の部屋」「カントリーロード」「さとがえり」「スターバトル」「タイルペイント」「チョコナ」「のりのり」「へやわけ」「ペイントエリア」
例：

```py3:python
data = """\
AABBB
AABCC
ADDDC
DDECC
EEEEC""".split()
v = np.array(addbinvars(nn, nn))
for _, d in groupby(sorted(zip(''.join(data), v.flat)), lambda x:x[0]):
    m += lpSum(c[1] for c in d) == 1
```

2次元のマスの表が上記 data のように、部屋に分かれています。このとき、「各部屋の中で黒マスは1つ」という条件は、groupbyを使うと上記のようにシンプルに書けます。

## 候補から選択
効用：簡単なモデル化
対象パズル：「カナオレ]「さとがえり]「四角に切れ]「ぬりかべ]「ののぐらむ]「フィルマット]「ブロックパズル]「ホタルビーム]

2つの0-1変数$x,y$があったとき、多くてもどちらか1方だけ成り立つ(=1)ようにしたければ、下記のようにできます。

例：

```py3:python
m += x + y <= 1
```

パズルのルールを数式で表すことが困難な場合、ルールを満たす組合せを列挙し、その中から選ばせることで簡単にモデル化できることがあります。このような定式化は[典型問題](0f6c1a4415d196e64314)の[集合分割問題](22ec0e42999141a0ba1e)にあたります。

「カナオレ」の候補の作成関数を見てみましょう。
例：

```py3:python
def makecand(lst, y, x, d, l, p, u):
    yy, xx = y+[0, -1, 0, 1][d], x+[-1, 0, 1, 0][d] # (1)
    if 0 <= yy < nh and 0 <= xx < nw and (yy,xx) not in u:
        if p == l - 1:
            lst.append(u + [(yy, xx)])
            return
        for k in range(4):
            makecand(lst, yy, xx, k, l, p + 1, u + [(yy,xx)])
```

引数のdが方向になっています。(1)で 1マス追加し、pが長さが達したらlstに追加して終了。そうでなければ、4方向に探索を繰り返します。
「カナオレ」をこのような方法以外でモデル化するのは、困難でしょう。
AMPLなどのモデリング言語では、柔軟な記述ができません。このようにモデル化でPythonを使うのは、非常に有用といえます。

以上