【R】数独で遊ぼう

## はじめに

あけましておめでとうございます。お正月に乗り遅れました。

正月と言えば数独です。
手で解くのは苦手なので、コンピュータ君に解いてもらいましょう。

## 数独のルール

- 任意の行、列内に同じ数字が入ってはいけない
- 太線で囲まれた$3\times3$のブロック内に同じ数字が入ってはいけない

シンプルですね。

## どうやって解くの？
整数計画法として定式化し、Rのパッケージ{lpSolveAPI}を使用して解を求めます。

参考：
[はじめよう整数計画法](https://web.tuat.ac.jp/~miya/fujie_OR2014spring.pdf)

### 定義
まずは、数独のルールを定式化するために変数$x_{ijk}$を定義します。[^1]
[^1]: $x_{ij}=k$というようにしても良いのですが、制約数が多くなるため計算量が大きくなります。
$x_{ijk}$は$9\times9\times9$の$3$次元配列をイメージすると分かりやすいです。

```math
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$となります。
![x_ijk_example.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/182504/9650743c-3473-34b1-b635-366c95def35e.jpeg)

###定式化
目的関数は指定する必要は無く、任意です。
目的関数、および制約は次の通りです。

```math
\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$のブロック内に同じ数字が入ってはいけない」に対応しています。

## スクリプト
[パッケージ{lpSolveAPI}の使い方](https://cran.r-project.org/web/packages/lpSolveAPI/lpSolveAPI.pdf)

```r:lpSolveAPI_sudoku.R
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"は不適？

```r:result
> 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