最適化問題の最短道問題のコードを読み解くメモ

#最適化問題を勉強中
数学は得意なので内容はすっと入ってくるんですが、そのコードを読み解くのって難しいですね。。。
最短道問題をハンズオンするために、[こちら](https://qiita.com/SaitoTsutomu/items/bfbf4c185ed7004b5721)のコードを読み解いています。
僕自身初心者なので理解するのに時間がかかるけど記録のためにメモメモ。

#コード
フルのコードは[ここから](https://qiita.com/SaitoTsutomu/items/bfbf4c185ed7004b5721#%E6%9C%80%E7%9F%AD%E8%B7%AF%E5%95%8F%E9%A1%8C)

##問題設定
ノード数8で0.26の確率で道(エッジと呼ぶ)が作られる問題設定。
LpProblemは引数無しだと最小化問題になる (引用は[ドキュメント](https://docs.pyq.jp/python/math_opt/pulp.html)から)
fast_gnp_random_graphメソッドのドキュメントは[ここから](https://networkx.org/documentation/networkx-1.10/reference/generated/networkx.generators.random_graphs.fast_gnp_random_graph.html)

```py
from pulp import *
import networkx as nx
g = nx.fast_gnp_random_graph(8, 0.26, 1).to_directed()
source, sink = 0, 2 # 始点, 終点
```
描画用のコードは[この記事](https://qiita.com/SaitoTsutomu/items/4d5715f6281be39f51c6)から引用

```py
import random, numpy as np, pandas as pd, networkx as nx, matplotlib.pyplot as plt
from itertools import chain, combinations
from pulp import *

def draw(g):
    """描画"""
    nx.draw_networkx_labels(g, pos=pos)
    nx.draw_networkx_nodes(g, node_color='w', pos=pos)
    nx.draw_networkx_edges(g, pos=pos)
    plt.show()
```
先ほどの問題設定を描画するとこんな感じ。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/666990/cbec2362-2e08-7896-5dab-4afe3c5da4d0.png)
fast_gnp_random_graphメソッドの第二引数を1（すべてのノードを道でつなぐ）にすると
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/666990/05ea252c-4c2f-e67a-f233-553c779bebaa.png)

##変数設定と目的関数のセット
この記事では0から1の間の変数とだけ定義しているから、整数とは限らなそう。
LpVaribaleのcat引数をIntegerにすれば整数問題にもできる様子。
参考はpulpの[ドキュメント](https://www.coin-or.org/PuLP/pulp.html)

```py
r = list(enumerate(g.edges()))
m = LpProblem() # 数理モデル
x = [LpVariable('x%d'%k, lowBound=0, upBound=1) for k, (i, j) in r] # 変数(路に入るかどうか)
m += lpSum(x) # 目的関数。xをぜえんぶ足す
```
ちなみに

```py
list(enumerate(g.edges())) #エッジ番号, (始点、終点)
>>>[(0, (0, 1)),
 (1, (1, 0)),
 (2, (1, 4)),
 (3, (1, 6)),
 (4, (1, 7)),
 (5, (2, 5)),
 (6, (3, 5)),
 (7, (3, 6)),
 (8, (4, 1)),
 (9, (5, 2)),
 (10, (5, 3)),
 (11, (6, 1)),
 (12, (6, 3)),
 (13, (7, 1))]
[LpVariable('x%d'%k, lowBound=0, upBound=1) for k, (i, j) in r]
>>>[x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13]
```

##制約条件足す
ここが難しかった。もっと勉強せねば・・・

```py
for nd in g.nodes():
    m += lpSum(x[k] for k, (i, j) in r if i == nd) \
      == lpSum(x[k] for k, (i, j) in r if j == nd) + {source:1, sink:-1}.get(nd, 0) # 制約
```
まず

```py
lpSum(x[k] for k, (i, j) in r if i == nd)) 
#ndというnodeから出るエッジ（道）の総和。1からでる道は4本でx1 + x2 + x3 +x4
>>>x0
x1 + x2 + x3 + x4
x5
x6 + x7
x8
x10 + x9
x11 + x12
x13
```
そして

```py
lpSum(x[k] for k, (i, j) in r if j == nd)
#ndというnodeへ行くエッジの総和
>>>x1
x0 + x11 + x13 + x8
x9
x10 + x12
x2
x5 + x6
x3 + x7
x4
```
経路において、通過するノードは出るノードから出る本数と入る本数が一致するはず。
始点の場合は出る本数のほうが一本多く、終点の場合は入る本数のほうが一本多いはず。
なので

```py
lpSum(x[k] for k, (i, j) in r if i == nd) \
      == lpSum(x[k] for k, (i, j) in r if j == nd) + {source:1, sink:-1}.get(nd, 0)
```
なお、辞書のgetメソッドはキーを探すためのメソッド。[この記事](https://note.nkmk.me/python-dict-get/)を参照。
終点始点の場合分けをしており、左辺（出る本数）と右辺（入る本数）に差をつけるための項。

##解く
valueはfrom pulp import * の*の中にあるモジュールのよう。
入った数と出た数の関係のみを参考にしているので、同じ場所は何回でも通ってよい＆行ってないノードの存在を許すという問題。

```py
m.solve()
print([(i, j) for k, (i, j) in r if value(x[k]) > 0.5])
```

#最後に
楽しい！