数独を通して組合せ最適化を学ぼう

Advent Calendar 2日目の記事 [組合せ最適化でカックロを解く](https://qiita.com/SaitoTsutomu/items/76e13aea15be17b9fd06)

## 数独を通して組合せ最適化を学ぼう

### i. 目的

「様々な課題を組合せ最適化を使って解決できたら」と思った人に向けて、**数独**[^1]を題材にして説明していきます。
一般に、パズルは解く過程を楽しむものです。
組合せ最適化を使って解くパズルは、モデル化の過程を楽しむものです。

[^1]: 数独は[株式会社二コリ](http://nikoli.co.jp/)の[登録商標](http://www.nikoli.co.jp/ja/copyright/)です。

### ii. 組合せ最適化とは

数理モデルを使った最適化で、下記のような、いろいろな問題があります。

- 参考：[典型問題と実行方法](https://qiita.com/SaitoTsutomu/items/0f6c1a4415d196e64314)

もう少し詳しい説明は、下記をご覧ください。

- 参考：[組合せ最適化を使おう](https://qiita.com/SaitoTsutomu/items/bfbf4c185ed7004b5721)

組合せ最適化は、難しい問題が多いですが、ソルバーを使うと数理モデルを作成するだけで、ソルバーが答えを出してくれます。

最近では、ソルバーの性能が上がってきて、過去、解けなかったものも解けるようになってきています。

- 参考：[組合せ最適化ソルバーの威力](https://qiita.com/SaitoTsutomu/items/82831e01adc3f84c36f5)

### iii. 手順について

課題を解くための手順です。

1. 問題をまとめる
    - ヒアリングしたり、頭を整理したりして、必要な情報をまとめていきます。
- 数理モデルを考える
    - 変数、目的関数、制約条件を考えます。たくさんのアプローチ方法を考えることが重要です。この能力は、経験によって鍛えられます。
- 実装して解く
    - ソルバーで実行できるモデルを記述します。お勧めするのは Pythonで記述する方法です。解くのはソルバーを呼ぶだけです。エラーを取り除くデバッグも含まれます。時間がかかるなどして解けない場合は、2.に戻って別のアプローチを探しましょう。
- 結果を確認する
    - 結果を見てモデルの通りに実装されているか確認しましょう。テスト結果がおかしければ、3.に戻って直してください。モデル通りにできていたら、次にその課題の専門家に結果を見てもらいましょう。ここでは、図示が有効です。専門家に おかしいと言われたら、1.に戻って見直します。

### 1. 問題をまとめる

今回は、下図のような数独を解きます。
![](https://raw.githubusercontent.com/SaitoTsutomu/opt4puzzle/master/pic/sudoku.png)

- 9x9の全マスに、1～9の数字を必ず埋めます。(1)
- どの1行(2)、どの1列(3)、どの3x3(4)において、同じ数字は1回だけ使えます。
- 数字が埋まっている場所は、その数字を使います。(5)

### 2. 数理モデルを考える

下図のような、9x9x9の729個の箱を考えましょう。この3軸は、行、列、数字に対応します。

<img src="https://qiita-image-store.s3.amazonaws.com/0/13955/8068e5a9-8170-ca97-befc-e2c8e349a810.png" width="300px">

1つの箱には、選ばれているか／選ばれていないか のどちらかの状態を持ちます。
これを1と0の数字で表しましょう。行i 列j 数字k の箱が1の場合、i行j列のマスの数字がkであることを意味するものとします。この729個の箱が**0-1変数**になります。
制約条件は、下記の通りです。

- 1つのマスに同じ数字は1つ。1x1x9の9箱の合計が1。 …(1)
- 1行のマスに同じ数字は1つ。1x9x1の9箱の合計が1。 …(2)
- 1列のマスに同じ数字は1つ。9x1x1の9箱の合計が1。 …(3)
- 3x3のマスに同じ数字は1つ。3x3x1の9箱の合計が1。 …(4)
- 数字指定の場合、その数字。1x1x1の1箱の合計が1。 …(5)

今回、目的関数はありません（便宜上、式 0 になります）。

### 3. 実装して解く

Pythonで数理モデルを作成する方法は、いろいろありますが、よく使われる PuLPとpandasを利用します。
PuLPについては、下記を参照ください。

- 参考：[最適化におけるPython - PuLPの使い方](https://qiita.com/SaitoTsutomu/items/070ca9cb37c6b2b492f0#pulp%E3%81%AB%E3%81%A4%E3%81%84%E3%81%A6)
- 参考：[数理最適化モデラー(PuLP)チートシート](https://qiita.com/SaitoTsutomu/items/c0bbf6cf8873ccd7edf3)

#### インストール

Pythonはバージョン3.6をインストールしてください。
追加ライブラリのために、下記を実行してください(処理系によっては、pipの代わりにpip3になります)。

> pip install pandas pulp ortoolpy

#### 変数表の作成

pandasで下記のような表を作成します。1レコード1変数(Var)に対応させ729レコードになります（1レコードは表の1行に対応します）。
<table>
<tr><th></th><th>行</th><th>列</th><th>数</th><th>_3x3</th><th>固</th><th>Var</th></tr>
<tr><th>0</th><td>0</td><td>0</td><td>1</td><td>0</td><td>False</td><td>v000001</td></tr>
<tr><th>0</th><td>0</td><td>0</td><td>2</td><td>0</td><td>False</td><td>v000002</td></tr>
<tr><th>...</th><td></td><td></td><td></td><td></td><td></td><td></td></tr>
</table>

この表を使うと制約条件は、以下のように表せます（カッコ内の数字は、`2. 数理モデルを考える`と対応します）。

- 1つのマスに同じ数字は1つ。「行と列」が同じレコード集合の変数の和は1。 …(1)
- 1行のマスに同じ数字は1つ。「行と数」が同じレコード集合の変数の和は1。 …(2)
- 1列のマスに同じ数字は1つ。「列と数」が同じレコード集合の変数の和は1。 …(3)
- 3x3のマスに同じ数字は1つ。「3x3の種別と数」が同じレコード集合の変数の和は1。 …(4)
- 数字指定の場合、その数字。「固」がTrueなら変数は1。 …(5)

キーが同じ集合は pandasの groupで、できます。

#### プログラム

```py3:python
import re, pandas as pd
from itertools import product
from pulp import LpProblem, lpSum, value
from ortoolpy import addbinvars

s = ('. . 6 |. . . |. . 1 '
     '. 7 . |. 6 . |. 5 . '
     '8 . . |1 . 3 |2 . . '
     '------+------+------'
     '. . 5 |. 4 . |8 . . '
     '. 4 . |7 . 2 |. 9 . '
     '. . 8 |. 1 . |7 . . '
     '------+------+------'
     '. . 1 |2 . 5 |. . 3 '
     '. 6 . |. 7 . |. 8 . '
     '2 . . |. . . |4 . . ')
data = re.sub(r'[^\d.]','',s) # 数字とドット以外を削除
r = range(9)
a = pd.DataFrame([(i,j,(i//3)*3+j//3,k+1,c==str(k+1))
    for (i,j),c in zip(product(r,r),data) for k in r],
    columns=['行','列','_3x3','数','固'])
a['Var'] = addbinvars(len(a))
m = LpProblem()
for cl in [('行','列'),('行','数'),('列','数'),('_3x3','数')]:
    for _,v in a.groupby(cl):
        m += lpSum(v.Var) == 1
for _,r in a[a.固==True].iterrows():
    m += r.Var == 1
m.solve() # ソルバーで求解
```

### 4. 結果を確認する

結果を列として表に追加しましょう。
結果が1（＝選ばれた）のレコードを抜き出して、9x9に整形して表示します。

```py3:python
a['Val'] = a.Var.apply(value)
print(a[a.Val>0.5].数.values.reshape(9,9))
>>>
[[5 3 6 8 2 7 9 4 1]
 [1 7 2 9 6 4 3 5 8]
 [8 9 4 1 5 3 2 6 7]
 [7 1 5 3 4 9 8 2 6]
 [6 4 3 7 8 2 1 9 5]
 [9 2 8 5 1 6 7 3 4]
 [4 8 1 2 9 5 6 7 3]
 [3 6 9 4 7 1 5 8 2]
 [2 5 7 6 3 8 4 1 9]]
```

### A. Q&A

#### Q: `4. 結果を確認する`で、`a.Val>0.5`は、`a.Val==1`ではないのですか？

- A: ソルバーでは、整数変数も実数として計算しています。そのため、ごくたまに、0.99999999などの出力になることもあります。従って、等号ではなく不等号で判断しています。

#### Q: `lpSum`とは何ですか？`sum`ではダメですか？

- A: 機能的には同じですが、`lpSum`の方が効率的なので、`lpSum`を使いましょう。

  - 参考：[数理モデルにおける変数の和](https://qiita.com/SaitoTsutomu/items/6701841122acc3130a29)

### B. 参考

- [数理最適化によるパズルの解法](https://qiita.com/SaitoTsutomu/items/0c0db8d22979fc9de8f4)
- [初歩からの数理モデル](https://qiita.com/SaitoTsutomu/items/98167dc25c7f133899f7)
- [魔方陣を通して組合せ最適化を学ぶ](https://qiita.com/SaitoTsutomu/items/4ce4559d1d0f41613d2a)