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

# <i class='fa fa-question-circle' /> 因子の部屋を解く
**因子の部屋**とは、数独に似たパズルです。このパズルも[組合せ最適化](http://qiita.com/Tsutomu-KKE@github/items/bfbf4c185ed7004b5721)を使って解けます。
参考
- [因子の部屋 wikipedia](https://ja.wikipedia.org/wiki/%E5%9B%A0%E5%AD%90%E3%81%AE%E9%83%A8%E5%B1%8B)
- [因子の部屋 ニコリ](http://www.nikoli.co.jp/ja/puzzles/inshi_no_heya/)
- [数独を組合せ最適で解く](http://qiita.com/Tsutomu-KKE@github/items/4f919f453aae95b3834b)

# <i class='fa fa-question-circle' /> 問題例
![image](https://qiita-image-store.s3.amazonaws.com/0/13955/4d7c337b-b419-93c9-145b-97cbf8384624.png)

ニコリ様から承諾頂き、例題をお借りしました。左が問題で右が回答です。
問題では、太い線で囲まれた部屋に分かれており、左上にヒントの数字があります。
ヒントは、その部屋内の数字のかけ算になっています。
5×5のサイズならば、使える数字は1から5です。数独と同じく縦1行や横1列で、同じ数字は使えません。

# <i class='fa fa-question-circle' /> 定式化
組合せ最適化のモデルで大事なのは、なるべく1次式で表すことです。
単純に変数のかけ算をモデル化すると、非線形最適化となり解くのが困難になります。
今回は、対数を使って1次式にします。つまり、下記のようにとらえます。

$2 \times 3 = 6$  →  $\log(2) + \log(3) = \log(6)$

このようにすれば、1次式で表現できます。ただし、このままですと無理数を使うので計算誤差が発生します。そこで、制約式を等号ではなく、微少の範囲に入るように指定します。

定式化は以下のようになります。

<table>
<tr><td rowspan="2">$\mbox{variables}$</td><td>$x_{ijk} \in \{0, 1\} ~ \forall i, j, k$</td><td>マスi,jが数字k+1か (1)</td></tr>
<tr><td>$y_{ij} \in \{1 \cdots n\} ~ \forall i, j$</td><td>マスi,jの数字 (2)</td></tr>
<tr><td>$\mbox{subject to}$</td><td>$\sum_k{x_{ijk}} = 1 ~ \forall i, j$</td><td>数字は1つ (3)</td></tr>
<tr><td></td><td>$\sum_k{x_{ikj}} = 1 ~ \forall i, j$</td><td>縦に同じ数字はない (4)</td></tr>
<tr><td></td><td>$\sum_k{x_{kij}} = 1 ~ \forall i, j$</td><td>横に同じ数字はない (5)</td></tr>
<tr><td></td><td colspan="2">$y_{ij}をx_{ijk}で表す$ (6)</td></tr>
<tr><td></td><td colspan="2">マスの積がヒントに等しい (7)</td></tr>
</table>

# <i class='fa fa-question-circle' /> Pythonで解く
[pulp](http://qiita.com/Tsutomu-KKE@github/items/bfbf4c185ed7004b5721#%E3%82%BD%E3%83%95%E3%83%88%E3%81%AE%E3%82%A4%E3%83%B3%E3%82%B9%E3%83%88%E3%83%BC%E3%83%AB)とpandasを使います。

問題は、文字列に入っているとします。

```py3:python
ch = """
ABBCD
AEEFD
GGHFD
IJHKK
ILHMM
""".strip().split('\n')
rooms = [6, 15, 1, 12, 20, 8, 10, 6, 4, 4, 15, 1, 10]
```

準備をします。

```py3: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
```

変数の表を作って見てみましょう。

```py3: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}$ に対応します。

<table>
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>縦</th>
      <th>横</th>
      <th>Vars</th>
      <th>Var</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0</td>
      <td>0</td>
      <td>[v1, v2, v3, v4, v5]</td>
      <td>v6</td>
    </tr>
    <tr>
      <th>1</th>
      <td>0</td>
      <td>1</td>
      <td>[v7, v8, v9, v10, v11]</td>
      <td>v12</td>
    </tr>
    <tr>
      <th>2</th>
      <td>0</td>
      <td>2</td>
      <td>[v13, v14, v15, v16, v17]</td>
      <td>v18</td>
    </tr>
  </tbody>
</table>

定式化して解きます。

```py3: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() # 求解
```

結果の表示をします。

```py3: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.]]
```

以上