- 本記事はProjectEulerの「100番以下の問題の説明は記載可能」という規定に基づいて回答のヒントが書かれていますので、自分である程度考えてみてから読まれることをお勧めします。
問題 96. 数独
問題の要約：ファイルから数独の問題50問を読み込んで、その答えの左上の3桁の数字の合計を求めよ
PuLPを用いた数理モデルの解法
数独と解くプログラムは世の中にたくさんあると思われますが、その中でも印象的だった「PuLPを用いた数理モデルの解法」を紹介します。詳細は「パズルでみる組合せ最適化のテクニック」等、Tsutomu Saito氏の投稿に詳しい説明がありますので参考にして見てください。
アルゴリズムの簡単な説明
数理モデルの考え方はちょっと独特なので慣れが必要です。ここでは簡単な例として$2 \times 2$の数独を考えてみます。入る数字は１と２のみで各行・列に１と２があるように配置するという超簡単な例です。
| 1 | - |
|---|---|
| - | - |
これに対して以下の8個の変数を作ります。例えば変数V_0_0_1は0行0列の数が1になる確率とみなし、値が1ならば数は1、値が0ならば1以外と考えます。
| 変数 | 1になるのは |
|---|---|
| V_0_0_1 | 0行0列の数が1 |
| V_0_0_2 | 0行0列の数が2 |
| V_0_1_1 | 0行1列の数が1 |
| V_0_1_2 | 0行1列の数が2 |
| V_1_0_1 | 1行0列の数が1 |
| V_1_0_2 | 1行0列の数が2 |
| V_1_1_1 | 1行1列の数が1 |
| V_1_1_2 | 1行1列の数が2 |
これらの変数に対して条件を連立方程式にしていきます。
#--- 行のルール
V_0_0_1+V_0_1_1 == 1 # 0行目には1は1つしかない
V_0_0_2+V_0_1_2 == 1 # 0行目には2は1つしかない
V_1_0_1+V_1_1_1 == 1 # 1行目には1は1つしかない
V_1_0_2+V_1_1_2 == 1 # 1行目には2は1つしかない
#--- 列のルール
V_0_0_1+V_0_0_1 == 1 # 0列目には1は1つしかない
V_1_0_2+V_1_0_2 == 1 # 0列目には2は1つしかない
V_0_1_1+V_0_1_1 == 1 # 1列目には1は1つしかない
V_1_1_2+V_1_1_2 == 1 # 1列目には2は1つしかない
#--- 一つのマスには一つの数
V_0_0_1+V_0_0_2 == 1 # 0行0列には1か2かどちらかが入る
V_0_1_1+V_0_1_2 == 1 # 0行1列には1か2かどちらかが入る
V_1_0_1+V_1_0_2 == 1 # 1行0列には1か2かどちらかが入る
V_1_0_1+V_1_1_2 == 1 # 1行1列には1か2かどちらかが入る
#---
V_0_0_1 == 1 # 0行0列が1になる確率は１
これらの連立方程式をsolveに渡すと答えを返してくれるという仕組みです。
数独はこれを$9 \times 9$にして$3 \times 3$のマスのルールを足したものなので、以下のようなコードになります。数独を解くプログラムがこれだけ簡潔に書けるのは凄いと思います。
- 変数はaddvinvarsを使って自動的に生成している
- 連立方程式の足し算はlpSum関数を使う。
- 結果の値はvalue関数をnumpyのarrayにvectorize関数を適用して取り出します
Google ColabでPuLPを使うためにはインストールが必要なので下の方にコードを添付しました。ファイルの読み込みも添付してあります。
from ortoolpy import addbinvars
from pulp import LpProblem, LpVariable, LpBinary, lpSum, lpDot, value
import numpy as np
import itertools
def SuDoku(sd):
m = LpProblem()
gd = np.array(addbinvars(9,9,9)) # Variable array
for i,j in itertools.product(range(9),range(9)):
m += lpSum(gd[:,i,j]) == 1 # all numbers in a column
m += lpSum(gd[i,:,j]) == 1 # all numbers in a row
m += lpSum(gd[i,j,:]) == 1 # 1 number in the grid
k, l = i//3*3, i%3*3
m += lpSum(gd[k:k+3,l:l+3,j]) == 1 # all numbers in 3x3
c = sd[i*9+j]
if int(c)>0:
m += gd[i,j,int(c)-1] == 1 # if a number already assigned
m.solve() # Call solver
ret = np.vectorize(value)(gd).astype(int).dot(range(1,10))
return "".join(map(str, list(itertools.chain.from_iterable(ret)))) # flatten the result
print(f"Answer: {sum([int(SuDoku(sd)[:3]) for sd in sdata])}")
PuLpとortoolpyのインストール
import pip
!pip install ortoolpy
!pip install pulp
数独問題ファイルの読み込み
# show upload dialog
from google.colab import files
uploaded = files.upload()
#---- read numbers from file to name[]
f = open("p096_sudoku.txt")
Flines = f.readlines()
f.close()
sdata = ["".join([Flines[g*10+i].strip() for i in range(1,10)]) for g in range(50)]
print(sdata[0])
（開発環境：Google Colab）