Edited at

「因子の部屋」を組合せ最適で解く

More than 3 years have passed since last update.


因子の部屋を解く

因子の部屋とは、数独に似たパズルです。このパズルも組合せ最適化を使って解けます。

参考

- 因子の部屋 wikipedia

- 因子の部屋 ニコリ

- 数独を組合せ最適で解く


問題例

image

ニコリ様から承諾頂き、例題をお借りしました。左が問題で右が回答です。

問題では、太い線で囲まれた部屋に分かれており、左上にヒントの数字があります。

ヒントは、その部屋内の数字のかけ算になっています。

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.]]

以上