今回の内容
前回の記事で，ペットのキャンディー配置を最適化するために問題を数式として書き下す（定式化する）ところまで行いました．
今回は，実際にコードを書いて問題を解くところまで行います．
使用言語はPython2.x系で，主に使用する外部モジュールはPuLPです．
PythonやPuLPのインストール方法などは以下のページを参考にして下さい（日本語がいい場合は検索するとすぐ出てきます）．
* https://www.python.org/
* https://github.com/coin-or/pulp
PuLPとは
GitHubのページにもありますが，Pythonで書かれた線形計画モデラです（といいつつ整数変数も扱えるため、整数計画、混合整数計画もいけます）．
つまり，前回の記事で書いたような定式化を，実際にソルバー（最適化問題を解いてくれるプログラム，アプリケーション）に解かせるために分かる形で記述してやる（プログラミング）ためのものです．
ということは，ソルバーは別途用意しなければならないはずですが，実はPuLPはインストールした時点で使えるソルバーが用意されているため，今回のような小規模な問題を解くのにはそのデフォルトのソルバーで問題ありません．
もちろん，自前で用意したソルバーを指定して解くこともできます（ソルバーの指定方法などはマニュアルを参照して下さい）．学術機関に所属されている方は商用最速のGurobiというソルバーが無料で使えるし圧倒的に速いと思われます．モデリングもgurobipyという専用のモジュールが用意されており，モデリング自体もPuLPを経由してGurobiを呼び出すより高速です．
前回の定式化おさらい
前回の定式化は以下の様なものでした．
\begin{align}
\mathbb{maximize} & & \sum_{i \in I} \sum_{j \in J} \sum_{p \in P} v_{p} x_{ijp} & \\
\mathbb{subject \ to} & & x_{ijp} \leq a_{ijp} & & \forall i \in I, j \in J, p \in P \\
& & \sum_{k \in I} \sum_{l \in J} \sum_{p \in P} c_{klp}^{ij} x_{klp} \leq 1 & & \forall i \in I, j \in J \\
& & b_{p} \leq \sum_{i \in I} \sum_{j \in J} x_{ijp} \leq n_{p} & & \forall p \in P \\
& & \sum_{i \in I} \sum_{j \in J} \sum_{p \in P_{s}} x_{ijp} \leq m_{s} & & \forall s \in S \\
& & x_{ijp} \in \{ 0, 1 \} & & \forall i \in I, j \in J, p \in P \\
\end{align}
プログラム
では，実際にコードを書いて解いていきたいと思います．
せっかくオブジェクト指向言語を使っていますが，流れがわかりやすいと思うので今回はクラス化したりしないで書いています．実際には，キャンディーボックスなどはクラス化した方がかなり扱いやすいと思います．
モジュール準備
まず，モジュールなどを読み込んでおきます．
futureは（主に私が）のちのちPython3系に移行するときのために備えているものです（つまりなくてもよいです）．
repoze.lruについては，後で出てくる関数の高速化のために使っています（なくても動きますが，あると速いのでインストールしておくことをおすすめします）．
# -*- coding: utf-8 -*-
from __future__ import unicode_literals, print_function
from repoze.lru import lru_cache
from pulp import *
データの準備
キャンディーやキャンディーボックスなどのデータを用意します．
キャンディーボックス
# ワンダのキャンディーボックス。0がキャンディーを置ける場所、1が置けない場所
candy_box = [[1, 0, 0, 1, 1, 1, 0, 1],
[0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 1],
[0, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 0],
[1, 0, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 0, 0],
[1, 0, 1, 1, 1, 0, 0, 1]]
candy_boxはそのまま，今回解きたいペットのキャンディーボックスのデータです．
他のペットの場合は以下のようになりますので，解きたいペットのものに入れ替えれば解けます．
trim = [[1, 0, 0, 0, 0, 1, 1, 1],
[0, 0, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 1],
[0, 0, 0, 0, 0, 1, 0, 0],
[1, 0, 0, 0, 1, 1, 1, 0],
[1, 1, 0, 0, 0, 1, 0, 0],
[1, 0, 0, 1, 0, 0, 0, 0]]
sally = [[1, 1, 0, 1, 1, 0, 1, 1],
[0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 1, 0, 0],
[0, 0, 0, 1, 1, 0, 0, 0],
[1, 1, 0, 0, 0, 0, 1, 1]]
marron = [[0, 0, 0, 1, 1, 0, 0, 1],
[0, 1, 0, 0, 0, 0, 1, 1],
[0, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 1],
[1, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 0],
[1, 1, 0, 0, 0, 0, 1, 0],
[1, 0, 0, 1, 1, 0, 0, 0]]
melon = [[1, 1, 0, 1, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 1, 1, 0],
[0, 0, 0, 0, 0, 0, 0, 0],
[1, 0, 0, 0, 0, 0, 0, 1],
[1, 0, 0, 0, 0, 0, 0, 1],
[0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 1, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 1, 0, 1, 1]]
rappy = [[0, 1, 0, 0, 0, 0, 1, 0],
[0, 1, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 1, 0],
[1, 0, 1, 0, 0, 1, 0, 1],
[0, 1, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 1, 1, 0, 0, 0],
[0, 0, 0, 1, 1, 0, 0, 0]]
viola = [[0, 1, 1, 0, 1, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 1, 0],
[0, 1, 0, 0, 0, 0, 0, 1],
[0, 0, 0, 0, 0, 0, 0, 1],
[1, 0, 0, 0, 0, 0, 0, 0],
[1, 0, 0, 0, 0, 0, 1, 0],
[0, 1, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 1, 0, 1, 1, 0]]
キャンディー
# キャンディーの一覧。辞書の辞書形式。名前、height、width、profitは必須、min_set、max_setは任意
candies = {"blast_parfait": {"height": 2, "width": 2, "profit": 100, "min_set": 0},
"megaton_parfait": {"height": 2, "width": 2, "profit": 100, "min_set": 0},
"hissatsu_parfait": {"height": 2, "width": 2, "profit": 100, "min_set": 0},
"pittari_parfait": {"height": 2, "width": 2, "profit": 100, "min_set": 0},
"gamushara_parfait": {"height": 2, "width": 2, "profit": 100, "min_set": 0},
"harikiri_parfait": {"height": 2, "width": 2, "profit": 100, "min_set": 0},
"girigiri_roll": {"height": 2, "width": 2, "profit": 50, "min_set": 0},
"funbari_roll": {"height": 2, "width": 2, "profit": 50, "min_set": 0},
"ouen_roll": {"height": 2, "width": 2, "profit": 50, "min_set": 0},
"muteki_roll": {"height": 2, "width": 2, "profit": 50, "min_set": 0},
"shikkari_roll": {"height": 2, "width": 2, "profit": 50, "min_set": 0},
"henkan_roll": {"height": 2, "width": 2, "profit": 50, "min_set": 0},
"migawari_roll": {"height": 2, "width": 2, "profit": 50, "min_set": 0},
"tsuyogari_roll": {"height": 2, "width": 2, "profit": 50, "min_set": 0},
"stamina_cookie": {"height": 2, "width": 2, "profit": 30, "min_set": 0},
"spirita_cookie": {"height": 2, "width": 2, "profit": 30, "min_set": 0},
"power_cookie": {"height": 2, "width": 2, "profit": 30, "min_set": 0},
"shoot_cookie": {"height": 2, "width": 2, "profit": 30, "min_set": 0},
"technique_cookie": {"height": 2, "width": 2, "profit": 30, "min_set": 0},
"arm_cookie": {"height": 2, "width": 2, "profit": 30, "min_set": 0},
"1stage_pancake": {"height": 2, "width": 2, "profit": 20, "min_set": 0},
"2stage_pancake": {"height": 2, "width": 2, "profit": 25, "min_set": 0},
"3stage_pancake": {"height": 2, "width": 2, "profit": 30, "min_set": 0},
"honoono_ramune": {"height": 1, "width": 1, "profit": 1, "min_set": 0},
"koorino_ramune": {"height": 1, "width": 1, "profit": 1, "min_set": 0},
"kaminarino_ramune": {"height": 1, "width": 1, "profit": 1, "min_set": 0},
"kazeno_ramune": {"height": 1, "width": 1, "profit": 1, "min_set": 0},
"hikarino_ramune": {"height": 1, "width": 1, "profit": 1, "min_set": 0},
"yamino_ramune": {"height": 1, "width": 1, "profit": 1, "min_set": 0},
"richna_crepe": {"height": 2, "width": 2, "profit": 20, "min_set": 0},
"otonano_crepe": {"height": 2, "width": 2, "profit": 20, "min_set": 0},
"lucky_crepe": {"height": 2, "width": 2, "profit": 20, "min_set": 0},
"spirita_candy": {"height": 1, "width": 4, "profit": 15, "min_set": 0},
"brow_resist_candy": {"height": 1, "width": 4, "profit": 15, "min_set": 0},
"shot_resist_candy": {"height": 1, "width": 4, "profit": 15, "min_set": 0},
"mind_resist_candy": {"height": 1, "width": 4, "profit": 15, "min_set": 0},
"ability_candy": {"height": 1, "width": 4, "profit": 15, "min_set": 0},
"body_sand": {"height": 1, "width": 2, "profit": 10, "min_set": 0},
"react_sand": {"height": 1, "width": 2, "profit": 10, "min_set": 0},
"mind_sand": {"height": 1, "width": 2, "profit": 10, "min_set": 0},
"stamina_gmi": {"height": 1, "width": 1, "profit": 1, "min_set": 0, "max_set": 20}
}
# キャンディーのタイプごとの配置数上限
max_set_by_types = {'parfait': 2, 'roll': 2, 'cookie': 3, 'pancake': 2,
'ramune': 1, 'crepe': 1, 'candy': 4, 'sand': 5, 'gmi': 20}
解く上で自分が考慮したいキャンディーの一覧を書いておきます．必要に応じて少なくすると良いでしょう．
それぞれのキャンディーの名前はなんでも良いのですが，種類名「＿gmi」や「＿ramune」などの右の「＿」以降の部分は，max_set_by_typesに書いてあるキャンディーの種類（9種類）の何れか正しい種類を付けておいて下さい．
楽をするための準備
後で楽をするために，上記データから加工してできるデータや，関数を定義しておきます．
# キャンディーボックスで配置可能（0）な座標のリスト
settable_pos = [(i, j) for i in xrange(len(candy_box)) for j in xrange(len(candy_box[i])) if candy_box[i][j] < .5]
@lru_cache(64 * 64)
def is_settable(i=0, j=0, height=1, width=1):
"""i行目、j列目のマスを左上の点として、指定の大きさのキャンディーをセットできるかどうかを返す
:param i: キャンディーをセットする最も左上のマスの行番号
:param j: キャンディーをセットする最も左上のマスの列番号
:param height: セットしたいキャンディーの縦の長さ
:param width: セットしたいキャンディーの横の長さ
:return: セット可能（True）かセット不可能（False）か
"""
global candy_box
try:
return 0 == sum(candy_box[i + k][j + l] for k in xrange(height) for l in xrange(width))
except Exception:
return False
- settable_posは，キャンディーボックスにおいてキャンディーを配置可能な座標$(i, j)$のリストです．
- 例えばワンダの場合，[(0, 1), (0, 2), (0, 6), ...]となっています．
- is_settableは，キャンディーボックスの座標$(i, j)$を左上の点として，縦の長さheight，横の長さwidthの大きさのキャンディーを配置可能かどうかを判定する関数です．
- 配置可能ならTrue，そうでなければFalseを返します．
モデルの作成
ではいよいよ数式を表現していきます．まずは，モデル自体を表すLpProblemクラスのオブジェクトmodelを生成しておきます．
# 最適化問題のモデル作成。価値の合計の最大化が目的
model = LpProblem('SetCandyProblem', LpMaximize)
- 'SetCandyProblem'というのはモデルの名前で，任意の文字列です．
- LpMaximizeというのは，今回の目的関数が最大化であることを示しています．最小化ならLpMinimizeです． これ以降，このmodelオブジェクトに対して制約を追加したり目的関数を追加していきます．
変数の作成
変数を表すクラスLpVariableのオブジェクトx[i, j, p]を生成しておきます．
# 変数xの作成
x = {}
for i, j in settable_pos:
for name, candy in candies.iteritems():
if is_settable(i, j, candy['height'], candy['width']):
# 0-1変数として定義
x[i, j, name] = LpVariable('x_({},{},{})'.format(i, j, name), cat=LpBinary)
else:
# そもそも配置不可能な場合は変数として定義しない（＝解きやすくなる）
# さらに、こうしておくと1番目の制約を書かなくてよくなる
x[i, j, name] = 0
ここで，LpVariable()に渡している引数の内前者はその変数の名前，後者は種類を表しており，LpBinaryとはその名の通り0-1変数を意味します．
名前は指定しない場合は適当な名前が付けられるため，省略可能です．
デフォルトの種類は連続変数（LpContinuous）のため，連続変数の場合は省略可能です．
その他，連続変数や整数変数（LpInteger）の場合には，その変数の上限（upBound，デフォルトはNoneで無限大），下限（lowBound，デフォルトはNoneで負の無限大）を指定することもできます．
2番目の制約
キャンディーボックスのある座標$(i,j)$に置けるキャンディーは1つのみである制約を追加します．
lin_expressions = {(i, j): LpAffineExpression() for i, j in settable_pos}
for i, j in settable_pos:
for name, candy in candies.iteritems():
if not is_settable(i, j, candy['height'], candy['width']):
continue
for k in xrange(candy['height']):
for l in xrange(candy['width']):
lin_expressions[i + k, j + l] += x[i, j, name]
for pos, lhs in lin_expressions.iteritems():
model += lhs <= 1, 'OnlySetConstraint:{}'.format(pos)
ここでは，制約は$(i,j)$の総ての組合せに対して存在するため，まず各組合せに対して線形式を表す（線形表現）LpAffineExpressionクラスのオブジェクトlin_expressions[i, j]を生成しています．これは，2番目の制約の左辺を表すためのオブジェクトです．
その後，全ての変数に対して着目し，その変数が1（つまりキャンディーを配置する）になった場合に専有する座標$(i+k, j+l)$の線形表現にその変数を追加指定ます．
最後に，各制約の左辺の作成が終わった段階で，それが1以下である（lhs <= 1）という制約をモデルに追加しています．「'OnlySetConstraint:{}'.format(pos)」というのは，その制約の名前を表していますが，なくても適当な名前が付けられます．
3番目の制約（左側）
キャンディーボックス全体でキャンディー$p$の配置個数は最低個数以上であるという制約を追加します．
for name, candy in candies.iteritems():
if candy.get('min_set', 0) <= 0:
continue
model += lpSum(x[i, j, name] for i, j in settable_pos) >= candy['min_set'], 'MinSetConstraint:{}'.format(name)
これはわかりやすいですね．
キャンディーごとにfor文を回し，対象キャンディーに関する変数全てを足したもの（つまりキャンディーを配置する個数）が下限（candy['min_set']）以上であることを表しています．
3番目の制約（右側）
同様に，キャンディーボックス全体でキャンディー$p$の配置個数は上限以下であるという制約を追加します．
for name, candy in candies.iteritems():
if candy.get('max_set') is None:
if name.endswith('_parfait'):
# パフェだけは上限1で決まっている
model += lpSum(x[i, j, name] for i, j in settable_pos) <= 1, 'MaxSetConstraint:{}'.format(name)
continue
model += lpSum(x[i, j, name] for i, j in settable_pos) <= candy['max_set'], 'MaxSetConstraint:{}'.format(name)
基本的には下限の場合と同様ですが，種類がパフェの場合だけは上限が1つにつき1個と仕様上決まっています．
4番目の制約
キャンディーボックス全体でキャンディーの種類ごとの配置個数が上限を超えないことを表す制約を追加します．
lin_expressions = {type_name: LpAffineExpression() for type_name in max_set_by_types.keys()}
for i, j in settable_pos:
for name, candy in candies.iteritems():
try:
lin_expressions[name[name.rfind('_')+1:]] += x[i, j, name]
except Exception:
print('Invalid Candy Name! Unknown candy type.: {}'.format(name))
for name, lhs in lin_expressions.iteritems():
model += lhs <= max_set_by_types[name], 'MaxSetByTypeConstraints:{}'.format(name)
2番目の制約ど同様に，今度は種類ごとに線形表現オブジェクトを作成しておきます．
その後，変数ごとに着目し，その変数のキャンディーの種類に合わせた線形表現オブジェクトにその変数を追加しています．
そして最後にモデルに制約として追加しています．
目的関数
目的関数を追加します．
# 目的関数
model += lpSum(candy['profit'] * x[i, j, name] for i, j in settable_pos for name, candy in candies.iteritems())
目的関数は，各キャンディーの価値candy['profit']とそれを置いた場合1，それ以外のとき0となる変数xの積の総和です．
そして以上でモデルの記述が終わりました．後は解いて結果を取ってくるだけです．
求解と結果の取得
# 求解とそのステータス（最適解がでたのかどうか）を取得
status = model.solve()
print('Status: {}'.format(LpStatus[status]))
print('CalculationTime: {:.3f}s'.format(model.solutionTime))
# 最適解が出ていなければ終了
if status != LpStatusOptimal:
exit(0)
# 最適解が出ていれば結果を出力していく
# 目的関数値出力
print('Optimal value: {}'.format(model.objective.value()))
# キャンディーとIDの対応出力
sorted_candy_names = sorted(candies.keys(), key=lambda name: name[name.rfind('_'):]+name[:name.rfind('_')])
print('xxx: Block (can not set)\nooo: Free space\n---: "Top-left" candy\'s space')
print('\n'.join('{:0>3}: {}'.format(i, name) for i, name in enumerate(sorted_candy_names)))
# キャンディーボックスの配置出力
solution_box = [['xxx' if e > 0.5 else 'ooo' for e in l] for l in candy_box]
for i, j in settable_pos:
for ID, name in enumerate(sorted_candy_names):
if value(x[i, j, name]) < 0.5:
continue
for k in xrange(candies[name]['height']):
for l in xrange(candies[name]['width']):
solution_box[i + k][j + l] = '---'
solution_box[i][j] = '{:0>3}'.format(ID)
print('\n-------------- Solution --------------')
print('\n'.join(', '.join(l) for l in solution_box))
print('-------------- End --------------')
以上のプログラムを実行すると，私の環境では以下の様な結果を得ました．
Status: Optimal
CalculationTime: 0.069s
Optimal value: 452.0
xxx: Block (can not set)
ooo: Free space
---: "Top-left" candy's space
000: ability_candy
001: brow_resist_candy
002: mind_resist_candy
003: shot_resist_candy
004: spirita_candy
005: arm_cookie
006: power_cookie
007: shoot_cookie
008: spirita_cookie
009: stamina_cookie
010: technique_cookie
011: lucky_crepe
012: otonano_crepe
013: richna_crepe
014: stamina_gmi
015: 1stage_pancake
016: 2stage_pancake
017: 3stage_pancake
018: blast_parfait
019: gamushara_parfait
020: harikiri_parfait
021: hissatsu_parfait
022: megaton_parfait
023: pittari_parfait
024: hikarino_ramune
025: honoono_ramune
026: kaminarino_ramune
027: kazeno_ramune
028: koorino_ramune
029: yamino_ramune
030: funbari_roll
031: girigiri_roll
032: henkan_roll
033: migawari_roll
034: muteki_roll
035: ouen_roll
036: shikkari_roll
037: tsuyogari_roll
038: body_sand
039: mind_sand
040: react_sand
-------------- Solution --------------
xxx, 005, ---, xxx, xxx, xxx, 029, xxx
014, ---, ---, 030, ---, 020, ---, 014
014, xxx, 014, ---, ---, ---, ---, xxx
014, xxx, 009, ---, 017, ---, 038, ---
038, ---, ---, ---, ---, ---, xxx, 014
xxx, 032, ---, 017, ---, 014, xxx, 014
014, ---, ---, ---, ---, 021, ---, 014
xxx, 014, xxx, xxx, xxx, ---, ---, xxx
-------------- End --------------
最適解が0.069秒で求まり，目的関数値は452であったということが読み取れます．
注意）あくまでサンプル的なデータで解いた結果のため，実際に使う際には自分でキャンディーの価値や最低配置個数などを調整する必要があります．
余力のある人向け
同じ目的関数値で他の解があるのではないか？と鋭い方は思ったかもしれません．その通りです．
余力のある人は，以下の様な内容に挑戦してみてはいかがでしょうか．
- 同じ目的関数値となるような解を全て出すようにする
- 同じ目的関数値であれば，できるだけ多くの種類のキャンディーを配置する
おわりに
さて，今回は前回の記事で示した定式化をもとに，実際にコードを書いて問題を解くところまで行いました．疲れてかなり端折ってしまったため，わかりにくい箇所もあるかと思います．
コメントなどでリクエスト頂ければ，時間のあるときにその箇所について掘り下げて追記しようと思いますので，質問などありましたらコメントをよろしくお願いいたします．
また，今回は説明しやすくするためにクラス化などは行っていませんが，そもそも最初に書いたものはキャンディーボックスなどはクラス化しているため，要望があればそちらも公開しようと思います．
お気づきの点などございましたらコメント等でご指摘をお願いいたします．