はじめに
あけましておめでとうございます。お正月に乗り遅れました。
正月と言えば数独です。
手で解くのは苦手なので、コンピュータ君に解いてもらいましょう。
数独のルール
- 任意の行、列内に同じ数字が入ってはいけない
- 太線で囲まれた$3\times3$のブロック内に同じ数字が入ってはいけない
シンプルですね。
どうやって解くの？
整数計画法として定式化し、Rのパッケージ{lpSolveAPI}を使用して解を求めます。
参考：
はじめよう整数計画法
定義
まずは、数独のルールを定式化するために変数$x_{ijk}$を定義します。1
x_{ijk} := \begin{cases}
1 & マス(i, j)にkが入っている \\
0 & {\rm otherwise}
\end{cases}
定義の通り、$x_{ijk}$は$0$か$1$のどちらをとります。
###例
例として、$3\times3$のマスを考えます。
下図左の$3\times3$のマスを$x_{ijk}$で表現します。
下図右は$x_{ijk}$を表しており、各$k$ごとに$x_{ijk}$を$3\times3$で表しています。
$i = 2,\ j = 3$の場合を考えると、マス$(2, 3)$には4が入っているので、定義より、$x_{234}=1$となります。
###定式化
目的関数は指定する必要は無く、任意です。
目的関数、および制約は次の通りです。
\begin{align}
{\rm min} \quad &{\rm arbitrary}\\
{\rm s.t.} \quad &\sum_{k=1}^{n} x_{ijk} = 1 &(i, j = 1,\cdots,n)\\
&\sum_{j=1}^{n} x_{ijk} = 1 &(i, k = 1,\cdots,n)\\
&\sum_{i=1}^{n} x_{ijk} = 1 &(j, k = 1,\cdots,n)\\
&\sum_{i=3(r-1)+1}^{3r}\sum_{j=3(c-1)+1}^{3c} x_{ijk}=1 &(r, c=1,2,3;, k=1,\cdots,n )\\
&x_{ijk} = 1 &(マス(i, j)にkが初期配置されている) \\
& x_{ijk} \in \{0, 1\} &(i,j,k = 1,\cdots,n)
\end{align}
制約の4行目は、「太線で囲まれた$3\times3$のブロック内に同じ数字が入ってはいけない」に対応しています。
スクリプト
library(lpSolveAPI)
dim.x <- dim.y <- dim.z <- 9
N.var <- dim.x * dim.y * dim.z
# array for 0-1variable
x.array <- array(rep(0, N.var),
dim = c(dim.x, dim.y, dim.z))
# make lpSolve linear program model object
lprec <- make.lp(ncol = N.var)
set.type(lprec, 1:N.var, "integer")
set.objfn(lprec, rep(1, N.var))
# 初期配置
# マス(i, j)に数kが入っている
board.init <- matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 3, 0, 1, 0,
6, 9, 0, 0, 2, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 4, 3,
0, 2, 0, 7, 6, 0, 0, 0, 0,
5, 0, 0, 9, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 7, 0, 0,
0, 0, 0, 0, 0, 0, 6, 0, 0,
0, 0, 4, 0, 0, 8, 0, 0, 0), nrow = 9, byrow = T)
# constraint
# sigma_k x_ijk = 1 (i, j = 1,2,...,9)
for (i in 1:9) {
for (j in 1:9) {
x.array[i, j, ] <- 1 # 制約に使う変数に1を立てる
add.constraint(lprec, rep(1, 9), "=", 1, which(x.array == 1))
x.array[i, j, ] <- 0 # 元に戻す
}
}
# sigma_j x_ijk = 1 (i, k = 1,2,...,9)
for (i in 1:9) {
for (k in 1:9) {
x.array[i, , k] <- 1
add.constraint(lprec, rep(1, 9), "=", 1, which(x.array == 1))
x.array[i, , k] <- 0
}
}
# sigma_i x_ijk = 1 (j, k = 1,2,...,9)
for (j in 1:9) {
for (k in 1:9) {
x.array[, j, k] <- 1
add.constraint(lprec, rep(1, 9), "=", 1, which(x.array == 1))
x.array[, j, k] <- 0
}
}
# sigma_r sigma_c x_ijk = 1 (r, c = 1,2,3; k = 1,2,...,9)
for (r in 1:3) {
for (c in 1:3) {
for (k in 1:9) {
idx.i <- (3*(r-1)+1):(3*r)
idx.j <- (3*(c-1)+1):(3*c)
x.array[idx.i, idx.j, k] <- 1
add.constraint(lprec, rep(1, 9), "=", 1, which(x.array == 1))
x.array[idx.i, idx.j, k] <- 0
}
}
}
# initial location
for (i in 1:9) {
for (j in 1:9) {
if (board.init[i, j] != 0) {
x.array[i, j, board.init[i, j]] <- 1
add.constraint(lprec, 1, "=", 1, which(x.array == 1))
x.array[i, j, board.init[i, j]] <- 0
}
}
}
# solve
solve(lprec)
var.result <- get.variables(lprec)
# result plot
board.3d <- array(var.result, dim = c(dim.x, dim.y, dim.z))
board.2d <- matrix(rep(0,dim.x*dim.y), nrow = dim.y)
for (i in 1:9) {
for (j in 1:9) {
for (k in 1:9) {
if (board.3d[i, j, k] == 1) board.2d[i, j] <- k
}
}
}
print(board.2d)
set.type(lprec, 1:N.var, "integer")の"integer"を"binary"にすると一晩経っても計算が終わらなかった...
0-1整数計画のつもりで"binary"は不適？
> print(board.2d)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 4 1 5 8 7 6 3 2 9
[2,] 2 8 7 4 9 3 5 1 6
[3,] 6 9 3 1 2 5 4 7 8
[4,] 7 6 1 5 8 2 9 4 3
[5,] 3 2 9 7 6 4 8 5 1
[6,] 5 4 8 9 3 1 2 6 7
[7,] 1 3 6 2 4 9 7 8 5
[8,] 8 5 2 3 1 7 6 9 4
[9,] 9 7 4 6 5 8 1 3 2
さいごに
数独、楽しいですね。
気が向いたら、スクリプトをすっきりさせたり、マスの大きさを$9\times9$から大きくしたり、$3$次元にしてみたり、いろんな条件を調べてみたりいろいろ遊んでみようと思います（ほんとはここで遊びたかったけど時間切れ）。
$9\times9$だけでなくいろんなタイプの数独があるそう。
参考文献
https://web.tuat.ac.jp/~miya/fujie_OR2014spring.pdf
https://cran.r-project.org/web/packages/lpSolveAPI/lpSolveAPI.pdf
http://www.sudokugame.org/
http://lpsolve.sourceforge.net/5.5/
https://ja.wikipedia.org/wiki/%E6%95%B0%E7%8B%AC
-
$x_{ij}=k$というようにしても良いのですが、制約数が多くなるため計算量が大きくなります。
$x_{ijk}$は$9\times9\times9$の$3$次元配列をイメージすると分かりやすいです。 ↩