因子の部屋を解く
因子の部屋とは、数独に似たパズルです。このパズルも組合せ最適化を使って解けます。
参考
- 因子の部屋 wikipedia
- 因子の部屋 ニコリ
- 数独を組合せ最適で解く
問題例
ニコリ様から承諾頂き、例題をお借りしました。左が問題で右が回答です。
問題では、太い線で囲まれた部屋に分かれており、左上にヒントの数字があります。
ヒントは、その部屋内の数字のかけ算になっています。
5×5のサイズならば、使える数字は1から5です。数独と同じく縦1行や横1列で、同じ数字は使えません。
定式化
組合せ最適化のモデルで大事なのは、なるべく1次式で表すことです。
単純に変数のかけ算をモデル化すると、非線形最適化となり解くのが困難になります。
今回は、対数を使って1次式にします。つまり、下記のようにとらえます。
$2 \times 3 = 6$ → $\log(2) + \log(3) = \log(6)$
このようにすれば、1次式で表現できます。ただし、このままですと無理数を使うので計算誤差が発生します。そこで、制約式を等号ではなく、微少の範囲に入るように指定します。
定式化は以下のようになります。
| $\mbox{variables}$ | $x_{ijk} \in \{0, 1\} ~ \forall i, j, k$ | マスi,jが数字k+1か (1) |
| $y_{ij} \in \{1 \cdots n\} ~ \forall i, j$ | マスi,jの数字 (2) | |
| $\mbox{subject to}$ | $\sum_k{x_{ijk}} = 1 ~ \forall i, j$ | 数字は1つ (3) |
| $\sum_k{x_{ikj}} = 1 ~ \forall i, j$ | 縦に同じ数字はない (4) | |
| $\sum_k{x_{kij}} = 1 ~ \forall i, j$ | 横に同じ数字はない (5) | |
| $y_{ij}をx_{ijk}で表す$ (6) | ||
| マスの積がヒントに等しい (7) | ||
Pythonで解く
pulpとpandasを使います。
問題は、文字列に入っているとします。
python
ch = """
ABBCD
AEEFD
GGHFD
IJHKK
ILHMM
""".strip().split('\n')
rooms = [6, 15, 1, 12, 20, 8, 10, 6, 4, 4, 15, 1, 10]
準備をします。
python
import pandas as pd
from collections import defaultdict
from pulp import *
from math import log
def addvar(lowBound=0, count=[0], *args, **kwargs):
count[0] += 1
return LpVariable('v%d'%count[0], lowBound=lowBound, *args, **kwargs)
nn, nb = len(ch), len(rooms) # 数字の個数、部屋の数
rn, rb = range(nn), range(nb)
lognn = [log(k + 1) for k in rn] # 1..nnのlog
logrm = [log(h) for h in rooms] # ヒントのlog
変数の表を作って見てみましょう。
python
a = pd.DataFrame([(i, j, [addvar(cat=LpBinary) for k in rn], addvar())
for i in rn for j in rn],
columns=['縦', '横', 'Vars', 'Var']) # (1), (2)
print(a[:3])
- 縦 $i$ 横 $j$ のVarsは [$x_{ijk} \forall k$] に、Varは $y_{ij}$ に対応します。
| 縦 | 横 | Vars | Var | |
|---|---|---|---|---|
| 0 | 0 | 0 | [v1, v2, v3, v4, v5] | v6 |
| 1 | 0 | 1 | [v7, v8, v9, v10, v11] | v12 |
| 2 | 0 | 2 | [v13, v14, v15, v16, v17] | v18 |
定式化して解きます。
python
m = LpProblem() # 数理モデル
dic = defaultdict(list) # 部屋ごとの変数(x)のリスト
for _, r in a.iterrows():
m += lpSum(r.Vars) == 1 # (3)
m += lpDot(rn, r.Vars) + 1 == r.Var # (6)
dic[ch[r.縦][r.横]].append(r.Vars)
for i in rn:
for k in rn:
m += lpSum(v[k] for v in a.query('縦==%d'%i).Vars) == 1 # (4)
m += lpSum(v[k] for v in a.query('横==%d'%i).Vars) == 1 # (5)
for h in rb:
c = lpSum(lpDot(lognn, v) for v in dic[chr(h + 65)]) # 数字の積のlog
m += c >= logrm[h] - 0.001 # (7)
m += c <= logrm[h] + 0.001 # (7)
m.solve() # 求解
結果の表示をします。
python
a['Val'] = a.Var.apply(value)
print(a.Val.reshape(5, -1))
>>>
[[ 2. 3. 5. 1. 4.]
[ 3. 5. 4. 2. 1.]
[ 5. 2. 1. 4. 3.]
[ 1. 4. 2. 3. 5.]
[ 4. 1. 3. 5. 2.]]
以上