数独をExact cover問題として解く（DLX package 使用）

[ペントミノをExact cover問題として解く（その１）](https://qiita.com/masa0599/items/7d5aa1318149f4ad9edf) ではペントミノをExact cover問題として解く方法を紹介しました、今回は数独も同様に解いてみたいと思います。Exact cover問題に馴染みのない方は前回の投稿を参照にしていただければと思います。


### 数独 N=4をExact cover問題に変換する

仕組みを理解するためにより簡単な例として数独N=4を考えます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1922159/0734c64e-7fdc-48c6-8c8d-53f08cf23f08.png)

この図のように各行、列、ブロックに1-4までの数字がすべてあるようにする問題です。これをExact cover問題として解くためには[Exact_cover#Sudoku (Wikipedia)](https://en.wikipedia.org/wiki/Exact_cover#Sudoku)に示されているように以下の4つの制約条件を考えます
 
1. 各セルには数が一つ入る (rxcy)
1. 各行には1-4の数字が一個づつ入る (rxny)
1. 各列には1-4の数字が一個づつ入る (cxny)
1. 各ブロックには1-4の数字が一個づつ入る (bxny)

このx,yには1-4が入れたものを制約条件としてDLXの列データに入力します。

### DLXのmoduleをインストールする

```python 
!pip install dlx
```
### DLXの列(column)のデータを作る

４つの制約条件のx,yに1-4が入るので各々16個で計64個のcolumnデータになります。
これをDLXに入力しますが、このデータの形式や、なぜcolidxを作るかは最初のリンクを参照してください。

```python
from itertools import product
from dlx import DLX
N, bs = 4, 2
constraints = [s[0]+str(x)+s[1]+str(y) for s, x, y in product(["rc","rn","cn","bn"], range(1,N+1),range(1,N+1))]
for i in range(0,64,16):  
  print(f"constraints[{i}-{i+15}] : {constraints[i:i+16]}")
columns = [(c, 0) for c in constraints]
print(f"columns = {columns}")
colidx = {cname[0]: i for i, cname in enumerate(columns)}    # column index
print(f"colidx = {colidx}")
dlx = DLX(columns)
# constraints[0-15] : ['r1c1', 'r1c2', 'r1c3', 'r1c4', 'r2c1', 'r2c2', 'r2c3', 'r2c4', 'r3c1', 'r3c2', 'r3c3', 'r3c4', 'r4c1', 'r4c2', 'r4c3', 'r4c4']
# constraints[16-31] : ['r1n1', 'r1n2', 'r1n3', 'r1n4', 'r2n1', 'r2n2', 'r2n3', 'r2n4', 'r3n1', 'r3n2', 'r3n3', 'r3n4', 'r4n1', 'r4n2', 'r4n3', 'r4n4']
# constraints[32-47] : ['c1n1', 'c1n2', 'c1n3', 'c1n4', 'c2n1', 'c2n2', 'c2n3', 'c2n4', 'c3n1', 'c3n2', 'c3n3', 'c3n4', 'c4n1', 'c4n2', 'c4n3', 'c4n4']
# constraints[48-63] : ['b1n1', 'b1n2', 'b1n3', 'b1n4', 'b2n1', 'b2n2', 'b2n3', 'b2n4', 'b3n1', 'b3n2', 'b3n3', 'b3n4', 'b4n1', 'b4n2', 'b4n3', 'b4n4']
# columns = [('r1c1', 0), ('r1c2', 0), ('r1c3', 0), ('r1c4', 0), ('r2c1', 0), ('r2c2', 0), ('r2c3', 0), ('r2c4', 0), ('r3c1', 0), ('r3c2', 0), ('r3c3', 0), ('r3c4', 0), ('r4c1', 0), ('r4c2', 0), ('r4c3', 0), ('r4c4', 0), ('r1n1', 0), ('r1n2', 0), ('r1n3', 0), ('r1n4', 0), ('r2n1', 0), ('r2n2', 0), ('r2n3', 0), ('r2n4', 0), ('r3n1', 0), ('r3n2', 0), ('r3n3', 0), ('r3n4', 0), ('r4n1', 0), ('r4n2', 0), ('r4n3', 0), ('r4n4', 0), ('c1n1', 0), ('c1n2', 0), ('c1n3', 0), ('c1n4', 0), ('c2n1', 0), ('c2n2', 0), ('c2n3', 0), ('c2n4', 0), ('c3n1', 0), ('c3n2', 0), ('c3n3', 0), ('c3n4', 0), ('c4n1', 0), ('c4n2', 0), ('c4n3', 0), ('c4n4', 0), ('b1n1', 0), ('b1n2', 0), ('b1n3', 0), ('b1n4', 0), ('b2n1', 0), ('b2n2', 0), ('b2n3', 0), ('b2n4', 0), ('b3n1', 0), ('b3n2', 0), ('b3n3', 0), ('b3n4', 0), ('b4n1', 0), ('b4n2', 0), ('b4n3', 0), ('b4n4', 0)]
# colidx = {'r1c1': 0, 'r1c2': 1, 'r1c3': 2, 'r1c4': 3, 'r2c1': 4, 'r2c2': 5, 'r2c3': 6, 'r2c4': 7, 'r3c1': 8, 'r3c2': 9, 'r3c3': 10, 'r3c4': 11, 'r4c1': 12, 'r4c2': 13, 'r4c3': 14, 'r4c4': 15, 'r1n1': 16, 'r1n2': 17, 'r1n3': 18, 'r1n4': 19, 'r2n1': 20, 'r2n2': 21, 'r2n3': 22, 'r2n4': 23, 'r3n1': 24, 'r3n2': 25, 'r3n3': 26, 'r3n4': 27, 'r4n1': 28, 'r4n2': 29, 'r4n3': 30, 'r4n4': 31, 'c1n1': 32, 'c1n2': 33, 'c1n3': 34, 'c1n4': 35, 'c2n1': 36, 'c2n2': 37, 'c2n3': 38, 'c2n4': 39, 'c3n1': 40, 'c3n2': 41, 'c3n3': 42, 'c3n4': 43, 'c4n1': 44, 'c4n2': 45, 'c4n3': 46, 'c4n4': 47, 'b1n1': 48, 'b1n2': 49, 'b1n3': 50, 'b1n4': 51, 'b2n1': 52, 'b2n2': 53, 'b2n3': 54, 'b2n4': 55, 'b3n1': 56, 'b3n2': 57, 'b3n3': 58, 'b3n4': 59, 'b4n1': 60, 'b4n2': 61, 'b4n3': 62, 'b4n4': 63}
```

### 行(row)のデータを作る(appendRow)

各行のデータをappendRowを使って入れていきます。行、列、数がすべて1-4の直積でrowデータを作るので総数は$4^3=64$です。


指定するのは列名ではなくて列indexなのでcolidxを使って変換します。

```python 
for i, (r,c,n) in enumerate(product(range(1,N+1), repeat=3)):
  b = ((r-1)//bs)*bs+((c-1)//bs)+1
  rc, rn, cn, bn = f"r{r}c{c}", f"r{r}n{n}",f"c{c}n{n}", f"b{b}n{n}"
  row = [colidx[rc], colidx[rn], colidx[cn], colidx[bn]]
  print(i, ":", f"r{r}c{c}n{n}", "---", rc, colidx[rc], rn,colidx[rn], cn, colidx[cn], bn, colidx[bn],row)
  dlx.appendRow(row)   # add row
# 0 : r1c1n1 --- r1c1 0 r1n1 16 c1n1 32 b1n1 48 [0, 16, 32, 48]
# 1 : r1c1n2 --- r1c1 0 r1n2 17 c1n2 33 b1n2 49 [0, 17, 33, 49]
# 2 : r1c1n3 --- r1c1 0 r1n3 18 c1n3 34 b1n3 50 [0, 18, 34, 50]
#      :       :       
# 63 : r4c4n4 --- r4c4 15 r4n4 31 c4n4 47 b4n4 63 [15, 31, 47, 63]
```
### 答えを得る(solve+getRowList)

とりあえずsolveで解答をすべて出してみると$288$個あることが分かりました。これは最初の数字が全く入っていない場合なので次に数字が入っている場合を考えます。

```python 
ans = []
for k, solution in enumerate(dlx.solve()):
  ans += [[sorted(dlx.getRowList(i))for i in solution]]
print(k+1)
# 288
```

### 数字の初期値がある場合

以下のように正の数値が入っている箇所は既に決まっているとすると、そこには他の数は入らないので、そういうRowをスキップするようにします。
例えばr1c1はn1なので(r1c1n2,r1c1n3,r1c1n4)はスキップします。

```
[1, 0, 0, 0]
[0, 0, 0, 2]
[0, 0, 4, 0]
[0, 3, 0, 0]
```
これは単純に最初に入っている数字を変数$fixed$で表し、スキップするコード2行をappendRowに渡すrowを作るコードのfor分のあとに加えるだけです。

```python 
fixed = "1000000200400300"
for i, (r,c,n) in enumerate(product(range(1,N+1), repeat=3)):
  x = int(fixed[(r-1)*N+c-1])     #--- if data[r,c] = x
  if x > 0 and n != x: continue   #---  remove data if n != x
    :      :
```

### 結果を表示する

これでdlx.solve()を呼んだ結果を表示します。無事答えが表示されました。

```python 
def getrcn(a):
  for s in a:
    if s[0] == 'r': r = int(s[1])
    if s[2] == 'c': c = int(s[3])
    if s[2] == 'n': n = int(s[3])
  return r, c, n

ansmtx = [[0 for _ in range(N)] for _ in range(N)]
for k, solution in enumerate(dlx.solve(),1):
  ans = [sorted(dlx.getRowList(i))for i in solution]
  for a in ans:
    r, c, n = getrcn(a)
    ansmtx[r-1][c-1] = n 
  print(f"----- Answer #{k}-----")
  for r in ansmtx:
    print(r)
# ----- Answer #1-----
[1, 2, 3, 4]
[3, 4, 1, 2]
[2, 1, 4, 3]
[4, 3, 2, 1]
```

### 9x9の数独を解く

基本的にコードは同じでN=9, bs=3とすれば良いだけです。[Project Euler Problem 96: Su Doku](https://projecteuler.net/problem=96)の例題を解くコードの全体を改めて示します。

あの数独かこれだけのコードで簡単に解けるのですから、かなりパワフルですよね。

```python 
from itertools import product
from dlx import DLX
def getrcn(a):
  for s in a:
    if s[0] == 'r': r = int(s[1])
    if s[2] == 'c': c = int(s[3])
    if s[2] == 'n': n = int(s[3])
  return r, c, n

def sudokuSolver(N, bs, data):
  #--------- Columns -----------
  constraints = [s[0]+str(x)+s[1]+str(y) for s, x, y in product(["rc","rn","cn","bn"], range(1,N+1),range(1,N+1))]
  columns = [(c, 0) for c in constraints]
  colidx = {cname[0]: i for i, cname in enumerate(columns)}    # column index
  dlx = DLX(columns)
  #--------- Row data ---------
  for i, (r,c,n) in enumerate(product(range(1,N+1), repeat=3)):
    x = int(data[(r-1)*N+c-1])       #--- if data[r,c] = x
    if x > 0 and  n != x: continue   #---  remove data if n != x
    b = ((r-1)//bs)*bs+((c-1)//bs)+1
    row = [colidx[f"r{r}c{c}"], colidx[f"r{r}n{n}"], colidx[f"c{c}n{n}"], colidx[ f"b{b}n{n}"]]
    dlx.appendRow(row)   # add row
  #------- Solve -------------
  ansmtx = [[0 for _ in range(N)] for _ in range(N)]
  for k, solution in enumerate(dlx.solve(),1):
    if k > 1: print("???? Not a single answer")
    for sol in solution:
      r, c, n = getrcn(dlx.getRowList(sol) )
      ansmtx[r-1][c-1] = n 
  return ansmtx
  
sdata = "003020600900305001001806400008102900700000008006708200002609500800203009005010300"
ansmtx = sudokuSolver(9, 3, sdata)
print(f"----- Answer -----")
for r in ansmtx:
  print(r)
# ----- Answer -----
# [4, 8, 3, 9, 2, 1, 6, 5, 7]
# [9, 6, 7, 3, 4, 5, 8, 2, 1]
# [2, 5, 1, 8, 7, 6, 4, 9, 3]
# [5, 4, 8, 1, 3, 2, 9, 7, 6]
# [7, 2, 9, 5, 6, 4, 1, 3, 8]
# [1, 3, 6, 7, 9, 8, 2, 4, 5]
# [3, 7, 2, 6, 8, 9, 5, 1, 4]
# [8, 1, 4, 2, 5, 3, 7, 6, 9]
# [6, 9, 5, 4, 1, 7, 3, 8, 2]
```

（開発環境：Google Colab）