お絵描きロジック(ピクロス)をpulpで解く

## お絵かきロジックとは
以下のルールでマスを埋めていくパズルの一種です．
イラストロジック，ののぐらむ，ピクロスと呼ばれることもあるようです.

1. タテ，ヨコ各列，数字の数だけマスを連続して黒く塗る
1. 数字が複数ある場合は，右から(or上から)数字の分だけ黒く塗る
1. ただし，数字同士がつながらないように間に少なくとも一つ空白を入れる

例えば，この問題は
![スクリーンショット 2022-03-13 20.01.56.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2596013/9e92c58b-0643-46f3-d5ce-fd764b5a17fe.png)
↓
このようになります．
![スクリーンショット 2022-03-13 20.20.24.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2596013/cbdd417a-47db-bfee-0e1a-47ebfd24487c.png)
引用問題: [https://www.minicgi.net/logic/logic.html?num=44396](https://www.minicgi.net/logic/logic.html?num=44396)

これをpythonの線形最適化パッケージであるpulpを使って解いてみましょう！

## 定式化

今回は目的関数なしの定式化になります．
まずは，問題を入力します．
数字がない場合は代わりに```[0]```を入力します．
```
row = [
    [4, 4],
    [3, 2, 3],
    [2, 4, 2],
    [1, 1],
    [1, 4, 1],
    [2, 4, 2],
    [2, 4, 2],
    [2, 4, 2],
    [2, 2],
    [10],
]

column = [
    [4, 5],
    [3, 6],
    [2, 1],
    [1, 1, 4, 1],
    [2, 4, 1],
    [2, 4, 1],
    [1, 1, 4, 1],
    [2, 1],
    [3, 6],
    [4, 5],
]
```

#### 変数
今回，連続して黒く塗る部分を右から（上から）セット1，セット2，...　と呼ぶことにし，次のように変数を用意します．
```math
x_{ij} = \left\{
\begin{array}{ll}
1 & (\ i\ 行\ j\ 列をマスを黒く塗る) \\
0 & (それ以外)
\end{array}
\right.
```
```math
y_{ijk} = \left\{
\begin{array}{ll}
1 & (\ i\ 行目の第\ k\ セットは\ j\ 列目から始まる) \\
0 & (それ以外)
\end{array}
\right.
```
```math
z_{ijk} = \left\{
\begin{array}{ll}
1 & (\ j\ 列目の第\ k\ セットは\ i\ 行目から始まる) \\
0 & (それ以外)
\end{array}
\right.
```

```
from pulp import LpProblem, LpBinary, LpVariable

problem = LpProblem("Picross")

x = [
    [LpVariable(f"x({i})({j})", cat=LpBinary) for j in range(len(column))]
    for i in range(len(row))
]

y = [
    [
        [LpVariable(f"y({i})({j})({k})", cat=LpBinary) for k in range(len(row[i]))]
        for j in range(len(column))
    ]
    for i in range(len(row))
]
z = [
    [
        [LpVariable(f"z({i})({j})({k})", cat=LpBinary) for k in range(len(column[j]))]
        for j in range(len(column))
    ]
    for i in range(len(row))
]
```

#### 制約条件1

第 $i$ 行において，第 $k$ セットが現れるのは1回のみである
第 $j$ 列において，第 $k$ セットが現れるのは1回のみである

```math
\sum_{j=1}^n y_{ijk}  = 1 \qquad for\ \ ^\forall i,\ 
^\forall k 
```
```math
\sum_{i=1}^n z_{ijk}  = 1 \qquad for\ \ ^\forall j,\ 
^\forall k 
```
ただし， 数字がないときは
```math
\sum_{j=1}^n x_{ij}  = 0 \qquad for\ \ ^\forall i
```
のようにします

```
from pulp import lpSum

for i in range(len(row)):
    if row[i] == [0]:
        problem += lpSum(x[i][j] for j in range(len(column))) == 0
    else:
        for k in range(len(row[i])):
            problem += lpSum(y[i][j][k] for j in range(len(column))) == 1

for j in range(len(column)):
    if column[j] == [0]:
        problem += lpSum(x[i][j] for i in range(len(row))) == 0
    else:
        for k in range(len(column[j])):
            problem += lpSum(z[i][j][k] for i in range(len(row))) == 1
```

#### 制約条件2

第 $i$ 行において，第　$k$　セットは第 $k+1$ セットの左に存在し，その間には一つ以上の空白が存在する
また，末尾のセットが長さ 2 以上の時は，外にはみ出さないようにする

例えば，```row[0]=[4,4]```なので
```math
y_{1j1} \leq \sum_{m=j+4+1}^{n} y_{1m2} \\
\sum_{m=n-4+2}^{n} y_{1m2} = 0
```
```
for i in range(len(row)):
    for j in range(len(column)):
        for k in range(len(row[i]) - 1):
            problem += y[i][j][k] <= lpSum(
                y[i][m][k + 1] for m in range(j + row[i][k] + 1, len(row))
            )
        problem += (
            lpSum(y[i][j][-1] for j in range(len(row) - row[i][-1] + 1, len(row))) == 0
        )

for i in range(len(row)):
    for j in range(len(column)):
        for k in range(len(column[j]) - 1):
            problem += z[i][j][k] <= lpSum(
                z[m][j][k + 1] for m in range(i + column[j][k] + 1, len(column))
            )
        problem += (
            lpSum(
                z[i][j][-1] for i in range(len(column) - column[j][-1] + 1, len(column))
            )
            == 0
        )
```
#### 制約条件3

マス$(i,j)$が黒く塗られる時， $i$ 行 $j$ 列のいづれかのセットに含まれる

例えば，```row[0]=[4,4]```のとき，
```math
x_{1j} \leq \sum_{m=max(0, j-4+1)}^{j} y_{1m1} + \sum_{m=max(0, j-4+1)}^{j} y_{1m2}
```
ここで，右辺の第一項が1になる時， $x_{ij} = 1$ は第1セットに含まれ，第2項が1になる時， $x_{ij} = 1$は第2セットに含まれることを意味します．

```
for i in range(len(row)):
    for j in range(len(column)):
        y_set = list()
        for k in range(len(row[i])):
            y_set += [y[i][m][k] for m in range(max(0, j - row[i][k] + 1), j + 1)]
        problem += x[i][j] <= lpSum(y_set)

for i in range(len(row)):
    for j in range(len(column)):
        y_set = list()
        for k in range(len(column[j])):
            y_set += [z[m][j][k] for m in range(max(0, i - column[j][k] + 1), i + 1)]
        problem += x[i][j] <= lpSum(y_set)
```

#### 制約条件4

$x_{ij} = 0$ の時，どのセットにも属しない

例えば，```row[0]=[4,4]```のとき，
```math
x_{ij} \geq \sum_{m=j-4+1}^{j}y_{ijk}
```

```
for i in range(len(row)):
    for j in range(len(column)):
        for k in range(len(row[i])):
            problem += x[i][j] >= lpSum(
                y[i][m][k] for m in range(j - row[i][k] + 1, j + 1)
            )

for i in range(len(row)):
    for j in range(len(column)):
        for k in range(len(column[j])):
            problem += x[i][j] >= lpSum(
                z[m][j][k] for m in range(i - column[j][k] + 1, i + 1)
            )
```

#### 求解

最後に解を求めます

```
from pulp import PULP_CBC_CMD, LpStatus

solver = PULP_CBC_CMD()
result = problem.solve(solver=solver)
status = LpStatus[result]

if status == "Optimal":
    for i in range(len(row)):
        box = ""
        for j in range(len(column)):
            if x[i][j].value() == 1:
                box += "■ "
            else:
                box += "□ "
        print(box)
else:
    print("解なし")
```

## 結果

無事求めることができました
```
■ ■ ■ ■ □ □ ■ ■ ■ ■ 
■ ■ ■ □ ■ ■ □ ■ ■ ■ 
■ ■ □ ■ ■ ■ ■ □ ■ ■ 
■ □ □ □ □ □ □ □ □ ■ 
□ ■ □ ■ ■ ■ ■ □ ■ □ 
■ ■ □ ■ ■ ■ ■ □ ■ ■ 
■ ■ □ ■ ■ ■ ■ □ ■ ■ 
■ ■ □ ■ ■ ■ ■ □ ■ ■ 
■ ■ □ □ □ □ □ □ ■ ■ 
■ ■ ■ ■ ■ ■ ■ ■ ■ ■ 
```

少し大きな問題でも1s以内で解けています
```
Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.86   (Wallclock seconds):       0.88

□ ■ ■ ■ ■ ■ □ ■ ■ ■ ■ ■ ■ ■ ■ □ ■ ■ ■ □ 
■ ■ □ □ □ ■ ■ ■ □ □ □ □ □ □ □ ■ □ □ ■ ■ 
■ □ □ □ □ □ ■ □ □ ■ ■ □ □ ■ ■ □ ■ □ □ ■ 
■ □ □ □ □ □ □ □ □ ■ □ □ □ ■ □ □ □ □ □ ■ 
■ □ □ □ □ □ □ □ □ ■ ■ □ □ ■ ■ □ □ □ ■ ■ 
■ ■ □ □ □ □ □ □ □ ■ ■ □ □ ■ ■ □ □ □ □ ■ 
□ ■ ■ □ □ ■ □ ■ □ □ □ □ □ □ □ □ ■ □ ■ ■ 
□ ■ ■ □ □ ■ □ ■ □ □ ■ ■ ■ ■ □ □ ■ □ ■ ■ 
□ □ ■ □ □ □ □ □ □ □ ■ ■ ■ ■ □ □ □ □ □ ■ 
□ □ ■ □ □ □ □ □ □ □ □ ■ ■ □ □ □ □ □ □ ■ 
□ □ ■ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ ■ 
□ □ ■ ■ □ □ □ □ □ □ □ □ □ □ □ □ □ □ ■ ■ 
□ □ □ ■ □ □ □ □ □ □ □ □ □ □ □ □ □ □ ■ □ 
□ □ □ ■ ■ □ □ □ □ □ □ □ □ □ □ □ □ ■ ■ □ 
□ □ ■ ■ ■ ■ ■ □ □ □ □ □ □ □ □ □ ■ ■ □ □ 
□ ■ ■ ■ ■ ■ ■ ■ ■ □ □ □ □ ■ ■ ■ ■ □ □ □ 
□ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ □ □ □ 
□ ■ ■ ■ ■ ■ ■ ■ □ □ □ ■ ■ ■ ■ ■ ■ □ □ □ 
□ ■ ■ ■ ■ ■ ■ □ □ □ □ ■ ■ ■ ■ ■ ■ □ □ □ 
□ □ ■ ■ ■ ■ □ □ □ □ □ □ ■ ■ ■ ■ □ □ □ □ 
```
結果が目に見えるっていいですね！
これはなんでしょう...?

それでは！