PuLP で数独を解く

# なぜ数独なのか

前回紹介した Python の線形計画法ライブラリ PuLP
http://qiita.com/mzmttks/items/82ea3a51e4dbea8fbc17
で、数独が解けるので解いてみる。

PyConJP のトーク 数理最適化によるパズルの解法
https://pycon.jp/2014/schedule/presentation/23/
や、ドキュメント
https://pythonhosted.org/PuLP/CaseStudies/a_sudoku_problem.html
に書かれているので自分でもやってみた。
という話。


# 数独って何？
1-9 までの数字を制約に基づいて埋めるパズルゲーム。

これらが制約
1. 横一列に同じ数字はこない
2. 縦一列に同じ数字はこない
3. 3x3 のブロック内に同じ数字はこない

ドキュメントの例
[数独](https://pythonhosted.org/PuLP/_images/wikisudokuproblem.jpg)

これを解くと、
[数独解](https://pythonhosted.org/PuLP/_images/wikisudokusolution.jpg)
こうなる。で、これが出れば OK

# 数独を線形計画問題に変換する
いわゆるモデル化。これがこの手の問題のキモで、結構感動した。

まずは PuLP で問題オブジェクトを作る

```py
import pulp
prob = pulp.LpProblem('Sudoku', pulp.LpMinimize)
prob += 0 # No objective function
```
目的関数は特にないので、定数の 0　とか設定しておく。

## 変数の設計
キーアイデアは、 **0 か 1 の値のみとるバイナリ変数を大量に使う** ということ。

具体的には、横9マス, 縦9マス, 入る値が 1-9 の 9 種類なので、
9 x 9 x 9 = 729 個の変数を用意して、
横 x, 縦 y の位置のマスの数字が v のとき変数 `Cell_v_x_y = 1`
そうでないとき `Cell_v_x_y = 0` とする。

つまり、1つのマスごとに 9個の変数を用意しておいて、
それらのうち 1個だけが 1, それ以外の 8個が 0 になる。


変数たちを生成するコードは下記のとおり

```py
numbers = range(1, 10)
xs = range(1, 10)
ys = range(1, 10)
choices = pulp.LpVariable.dicts("Cell",(xs, ys, numbers),0,1,pulp.LpInteger)
```

これで、Cell_1_1_1, ..., Cell_9_9_9 の 729 個の変数たちが生成されて、
`choices[1][1][1]` のような感じでアクセスできるようになった。
変数は整数で、0か1の値をとる。

## ゲームのルールから制約条件の設計
### 1 つのマスに入る値は 1 つだけ
たとえば、座標 (2, 3) のマスに入るのは 1-9 のうち 1つだけなので、
このような式になる。

```tex
Cell_1_2_3 + Cell_2_2_3 + Cell_3_2_3 +
Cell_4_2_3 + Cell_5_2_3 + Cell_6_2_3 +
Cell_7_2_3 + Cell_8_2_3 + Cell_9_2_3 = 1
```

コードで書くと

```py
ys = range(1, 10)
xs = range(1, 10)
numbers = range(1, 10)
for y in ys: # 各列と
    for x in xs: # 各行について、
        prob += pulp.lpSum([choices[v][x][y] for v in numbers]) == 1 
        # すべての数字を足すと 1 になる。 
```

ここで、`lpSum`は、配列の変数をすべて足した式を表す。
総和を表すときにシンプルになるので多用される。

### 既に数字が埋められているときは必ずその数字のみがマスに入る。
`board[y][x]` に、既に埋められているマスの数字が入っていて、
空だったら 0 が入っているとする。

```py
for x in range(width): # 各列と
    for y in range(height): # 各列を調べて、
        if board[y][x] > 0: # もし数字が 0 より大きかったら、
            prob += choices[board[y][x]][x+1][y+1] == 1
            # その場所の数字が 1 になるような制約を追加する。
```
すでに、同じマスで 1 になるのは 1 つの値だけという制約を入れたので、
値が 1 のときのことだけ考えればよい。


### 同じ列・行・ボックスには同じ数字は入らない

#### 行の制約
例：1列目(y)で、数字(v)が 1 になるのは、すべての行のうちどれか1つだけ。
つまり、それらを全部足すと必ず 1 になる。
`Cell_v_x_y` の変数で v=1, y=1 のとき、下記のような制約になる

```latex
Cell_1_1_1 + Cell_1_2_1 + Cell_1_3_1 + 
Cell_1_4_1 + Cell_1_5_1 + Cell_1_6_1 + 
Cell_1_7_1 + Cell_1_8_1 + Cell_1_9_1  = 1
```

```py
for v in numbers:
    for y in ys:
        prob += pulp.lpSum([choices[v][x][y] for x in xs]) == 1
```

列とボックスは省略。

# コード
```py
import pulp
import numpy

# make problem
board = """
530070000
600195000
098000060
800060003
400803001
700020006
060000280
000419005
000080079"""

# board = map(lambda line: map(int, line.rstrip()), open("sudoku.txt").readlines())
board = [map(int, list(line)) for line in board.strip().split("\n")]
print board
width = len(board[0])
height = len(board)


## initialize
# Create a problem
prob = pulp.LpProblem('Sudoku', pulp.LpMinimize) # or minimize

## objective
# No objective function
prob += 0

# make 
numbers = range(1, 10)
xs = range(1, 10)
ys = range(1, 10)
choices = pulp.LpVariable.dicts("Choice",(xs, ys, numbers),0,1,pulp.LpInteger)


## constraints
# Add given-value constraints
for x in range(width):
    for y in range(height):
        if board[y][x] > 0:
            prob += choices[board[y][x]][x+1][y+1] == 1

# one cell must have one value
for y in ys:
    for x in xs:
        prob += pulp.lpSum([choices[v][x][y] for v in numbers]) == 1, ""

# horizontal line must have different values
for v in numbers:
    for y in ys:
        prob += pulp.lpSum([choices[v][x][y] for x in xs]) == 1

    for x in xs:
        prob += pulp.lpSum([choices[v][x][y] for y in ys]) == 1

    for x in [1, 4, 7]:
        vs = []
        for y in [1, 4, 7]:
            print [[choices[v][x+i][y+j] for i in range(3) for j in range(3)]]
            prob += pulp.lpSum([[choices[v][x+i][y+j] for i in range(3) for j in range(3)]]) == 1


# print prob

s = prob.solve()
# solve it

print pulp.LpStatus[s]

# print result
for y in ys:
    for x in xs:
        for v in numbers:
            if choices[v][x][y].value() == 1:
                print v,
    print
```

実行結果は こんな感じになって、解けましたとさ。

```py
5 3 4 6 7 8 9 1 2
6 7 2 1 9 5 3 4 8
1 9 8 3 4 2 5 6 7
8 5 9 7 6 1 4 2 3
4 2 6 8 5 3 7 9 1
7 1 3 9 2 4 8 5 6
9 6 1 5 3 7 2 8 4
2 8 7 4 1 9 6 3 5
3 4 5 2 8 6 1 7 9
```