因子の部屋を解く

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

問題例

ニコリ様から承諾頂き、例題をお借りしました。左が問題で右が回答です。
問題では、太い線で囲まれた部屋に分かれており、左上にヒントの数字があります。
ヒントは、その部屋内の数字のかけ算になっています。
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)
$\mbox{variables}$	$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	[v1, v2, v3, v4, v5]	v6
1	1	[v7, v8, v9, v10, v11]	v12
2	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.]]

以上