import networkx as nx
import matplotlib.pyplot as plt
import random
最小木問題とは
木 (tree)とは，閉路を含まない連結グラフを指す． ここで，閉路とは始点と終点が一致するパスであり，連結グラフとは，すべての点の間にパスが存在するグラフである． また，与えられた無向グラフのすべての点を繋ぐ木を，特に全域木 (spanning tree)とよぶ． 全域木の概念を用いると，最小木問題は以下のように定義できる．
最小木問題 （minimum spanning tree problem）
$n$ 個の点から構成される点集合 $V$， $m$ 本の枝から構成される枝集合 $E$， 無向グラフ $G=(V,E)$, 枝上の費用関数 $c: E \rightarrow \mathbf{R}$ が与えられたとき， 枝の費用の合計を最小にする全域木を求めよ．
枝の費用の合計を最小にする全域木を最小木 （minimum spanning tree）とよぶ． 本来ならば，最小費用全域木問題と訳すべきであるが，簡単のため最小木問題と よぶことが多いので，ここでも慣例にならうものとする．
最小木問題は，人間の自然な摂理（貪欲性と改善性）にしたがうことによって，簡単に解くことができる．
はじめに，貪欲性に基づくものを紹介しよう． まず，枝上で定義される費用を小さい順に並べ替える． 次に，木を空集合に初期設定して， 費用の小さい順に，枝を加えていく． この際，枝を加えることによって閉路ができてしまう場合には，木に加えない．（正確にはアルゴリズムの途中では連結とは限らないでので，木ではなく，森（forest）とよぶ．） すべての点が木に含まれたら終了する．このときの全域木が，最小木になっている．
この方法は，貪欲アルゴリズム（greedy algorithm），もしくは発案者の名前をとって Kruskal法（Kruskal method）とよばれ，最小木問題の最適解を与える．
次に，人間のもう1つの自然な摂理である改善性に基づくものを紹介しよう． まず，適当な全域木からはじめる．次に，木に含まれていない枝を追加したときにできる閉路から， 最も費用の大きい枝を除くことを試みる（これを改善操作とよぶ）． もし，改善操作によって費用が減少するなら，その操作を実行する． すべての木に含まれていない枝に対して，改善操作が費用の減少をもたらさなかったら終了する． このときの全域木が，最小木になっている．
閉路除去定式化
最初の定式化は，点の部分集合に対して閉路を含まないことを表現したものであり， 閉路除去定式化（circuit elimination formulation）とよばれる．
定式化に必要な記号を導入しておく． 点集合 $V$ の要素の数（位数）を $n$ とする． 枝集合 $E$ の要素の数（位数）を $m$ とする． 点の部分集合 $S$ に対して，両端点が $S$ に含まれる枝の集合を $E(S)$ と書く． また，枝 $e \in E$ が全域木に含まれるとき $1$，それ以外のとき $0$ を表す $0$-$1$ 変数 $x_e$ を導入する．
木は，閉路を含まないグラフであったので， 点集合の任意の部分集合 $S \subset V$ に対して，$S$ に両端点が含まれる枝の本数は， 点の数 $|S|$ から $1$ を減じた値以下である必要がある． また，全域木であるためには，枝数の合計はちょうど $n-1$ 本でなければならない．
上の議論から，以下の定式化を得る． $$ \begin{array}{l l l } minimize & \sum_{e \in E} c_e x_e & \\ s.t. & \sum_{e \in E} x_e = n-1 & \\ & \sum_{e \in E(S)} x_e \leq |S|-1 & \forall S \subset V \\ & x_e \in \{0,1\} & \forall e \in E \end{array} $$
カットセット定式化
点の部分集合 $S$ に対して，$\delta(S)$ を端点の1つが $S$ に含まれ，もう1つの端点が $S$ に含まれない枝の集合とする． 枝集合は，$S$ 内に両端点をもつ枝の集合 $E(S)$， $V \setminus S$ 内に両端点をもつ枝の集合 $E(V \setminus S)$ と $\delta(S)$ の和集合である． すなわち， $$ E= E(S) \cup E(V \setminus S) \cup \delta(S) $$ が成立する． 全域木の特性ベクトル $x$ に対しては， $\sum_{e \in E} x_e =n-1$ であり，さらに $S$ と $V \setminus S$ に対する部分閉路除去を表す制約$\sum_{e \in E(S)} x_e \leq |S|-1$ ならびに $\sum_{e \in E(V \setminus S)} x_e \leq |V \setminus S|-1$ であるので，制約$\sum_{e \in \delta(S)} x_e \geq 1$ が妥当不等式であることが分かる．これをカットセット制約とよぶ．
カットセット制約を用いた最小木問題の定式化は，以下のように書ける．
$$ \begin{array}{l l l } minimize & \sum_{e \in E} c_e x_e & \\ s.t. & \sum_{e \in E} x_e = n-1 & \\ & \sum_{e \in \delta(S)} x_e \geq 1 & \forall S \subset V \\ & x_e \in \{0,1\} & \forall e \in E \end{array} $$閉路除去制約は，カットセット制約より「真に」強い定式化である．
単品種流定式化
上の2つの定式化では，入力サイズの指数オーダーの制約式を必要とした． ここでは，多項式オーダーの制約式から構成される定式化を考える． カットセット制約は，グラフの任意のカットが $1$ 以上の容量をもつことを規定していたが， 最大フロー・最小カット定理から，フローを用いた定式化を自然に導くことができる．
いま，グラフ $G=(V,E)$ 内の特定の点 $0 (\in V)$ から，他のすべての点に $1$単位のフローを流すことを考える． 枝 $e =(i,j)$ を $i$ から $j$ へ流れるフロー量を $f_{ij}$ とする． 点 $0$ から出る（供給する）フロー量は $n-1$ 単位であり，それを他の各点で $1$単位ずつ消費するものとする． 点 $i$ から点$j$ もしくは点 $j$ から点$i$ へ，いずれかの方向にフローが流れているときに， 枝 $e=(i,j)$ が最小木に含まれものとすると， 以下の単品種流定式化（single commodity flow formulation）を得ることができる．
$$ \begin{array}{l l l } minimize & \sum_{e \in E} c_e x_e & \\ s.t. & \sum_{e \in E} x_e = n-1 & \\ & \sum_{j: (0,j) \in \delta(\{0\})} f_{0j} = n-1 & \\ & \sum_{(j,i) \in \delta (\{i\})} f_{ji} -\sum_{(i,j) \in \delta (\{i\})} f_{ij} = 1 & \forall i \in V \setminus \{0\} \\ & f_{ij} \leq (n-1) x_e & \forall e=(i,j) \in E \\ & f_{ji} \leq (n-1) x_e & \forall e=(i,j) \in E \\ & x_e \in \{0,1\} & \forall e \in E \\ & f_{ij} \geq 0, f_{ji} \geq 0 & \forall e=(i,j) \in E \end{array} $$この定式化は簡潔ではあるが弱い． 制約が付加された問題に対して対処しやすいという利点もある。
多品種流定式化
カットセット制約と同等の強さをもつ定式化を得るためには， 特定の点 $0 (\in V)$ から，他の点 $k (\in V \setminus \{0\})$ に流すフローを区別する必要がある． フローの出先は発地（origin），行き先は着地（destination）とよばれる． この場合は，発地はすべて特定の点$0$ であり，着地は他のすべての点である． 一般に，異なる発地と着地をもつフローは区別して扱う必要があり，これを品種（commodity）とよぶ． ここで考える定式化は，複数の品種をもつ問題を用いるので， 多品種流定式化（multi-commodity flow formulation）とよばれる．
点$0$ から点$k$ へ流れるフローが，枝$(i,j)$ 上を $i$ から $j$ の向きで通過する量を表す実数変数を $f_{ij}^k$ と したとき，多品種流定式化は，以下のようになる．
$$ \begin{array}{l l l } minimize & \sum_{e \in E} c_e x_e & \\ s.t. & \sum_{e \in E} x_e = n-1 & \\ & \sum_{(j,i) \in \delta (\{i\})} f_{ji}^k -\sum_{(i,j) \in \delta (\{i\})} f_{ij}^k = \left\{ \begin{array}{l } -1 \\ 0 \\ 1 \\ \end{array} \right. & \begin{array}{l } i=0, \forall k \in V \setminus \{0\} \\ \forall i \in V \setminus \{0,k\}, k \in V \setminus \{0\} \\ i=k, \forall k \in V \setminus \{0\} \\ \end{array} \\ & f_{ij}^k +f_{ji}^{k'} \leq x_e & \forall k,k' \in V \setminus \{0\}, e=(i,j) \in E \\ & x_e \in \{0,1\} & \forall e \in E \\ & f_{ij}^k \geq 0, f_{ji}^k \geq 0 & \forall k \in V \setminus \{0\}, e=(i,j) \in E \end{array} $$この定式化の線形最適化緩和問題は，（指数オーダーの制約をもつ）閉路除去定式化の線形最適化緩和問題と等しいことを示すことができる．
networkXの利用
以下の関数で最適化ができる．
- minimum_spanning_tree(G)は無向グラフGの最小重みの全域木（最小木）をグラフとして返す．
- minimum_spanning_edges(G)は無向グラフGの最小木を枝の集合として返す．
引数のalgorithmでアルゴリズムを選択できる．既定値は‘kruskal’でKruskal法（貪欲解法）である． 他には，‘prim’（Prim法）と‘boruvka’（Boruvka法）が選択できる．
以下の例題をminimum_spanning_tree関数を用いて解いてみる．
あなたは，あなたの母国と冷戦下にある某国に派遣されているスパイだ． いま，この国に派遣されているスパイは全部で $5$人いて，それぞれが 偽りの職業について諜報活動をしている．スパイ同士の連絡には秘密の 連絡法がそれぞれ決まっていて， 本国からの情報によれば，連絡にかかる費用は与えられている．
いま，あなたが得た新しい極秘情報を，他の $4$人のスパイに連絡せよという指令が伝えられた． ただし，昨今の不景気風はスパイ業界にも吹いているようで， なるべく費用を安くしなければならないというおまけつきである． どのように連絡をとれば，最小の費用で極秘情報を仲間に連絡できるだろうか?
G = nx.Graph()
G.add_edge("Arigator", "WhiteBear", weight=2)
G.add_edge("Arigator", "Bull", weight=1)
G.add_edge("Bull", "WhiteBear", weight=1)
G.add_edge("Bull", "Shark", weight=3)
G.add_edge("WhiteBear", "Condor", weight=3)
G.add_edge("WhiteBear", "Shark", weight=5)
G.add_edge("Shark", "Condor", weight=4)
print(nx.minimum_spanning_tree(G).edges())
m, n = 5, 5
lb, ub = 1, 20
random.seed(1)
G = nx.grid_2d_graph(m, n)
for (i, j) in G.edges():
G[i][j]["weight"] = random.randint(lb, ub)
pos = {(i, j): (i, j) for (i, j) in G.nodes()}
edges = list(nx.minimum_spanning_edges(G))
plt.figure()
nx.draw(G, pos=pos, node_size=100)
edge_labels = {}
for (i, j) in G.edges():
edge_labels[i, j] = f"{ G[i][j]['weight'] }"
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)
nx.draw(G, pos=pos, width=5, edgelist=edges, edge_color="orange")
plt.show()
weight = []
for (i, j, w) in edges:
weight.append((w["weight"], i, j))
weight.sort(reverse=True)
print("weight=", weight)
print("max distance=", weight[3 - 1][0])
G1 = nx.Graph()
for (w, i, j) in weight[3:]:
G1.add_edge(i, j)
nx.draw(G, pos=pos, node_size=100)
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)
nx.draw(G1, pos=pos, node_size=100, width=10, edge_color="orange")
有向最小木
有向木 (arborescence)とは，入次数が高々1の連結有向グラフである． また，与えられた無向グラフのすべての点を繋ぐ木を，特に全域有向木 (spanning arborescence)とよぶ． 全域木の概念を用いると，最小木問題は以下のように定義できる．
最小有向木問題 （minimum spanning arborescence problem）
$n$ 個の点から構成される点集合 $N$， $m$ 本の枝から構成される有向枝集合 $A$， 有向グラフ $D=(N,A)$, 枝上の費用関数 $c: A \rightarrow \mathbf{R}$ が与えられたとき， 枝の費用の合計を最小にする全域有向木 (arborescence) を求めよ．
枝の費用の合計を最小にする全域木を最小有向木 （minimum spanning arborescence）とよぶ． 有向最小木は，最小木のように貪欲解法では求めることはでいないが，Edmond法で多項式時間で求めることができる．
Edmondsはクラスであり，コンストラクタで有向グラフを与える． その後で，find_optimumメソッドで解を求める． 引数で最小化か最大化か，連結していない有向木(branching) を求めるか，有向木を求めるかが指定できる．
- attr : 最小化するための枝属性名を与える．既定値は 'weight'．
- kind : 最小化か最大化かを指定する． 'min'で最小化，'max'（既定値）で最大化する．
- style : 連結な有向木'arborescence'か，非連結を許す'branching'（既定値）かを指定する．
m, n = 3, 3
lb, ub = 1, 20
G = nx.grid_2d_graph(m, n)
D = G.to_directed()
for (i, j) in D.edges():
D[i][j]["weight"] = random.randint(lb, ub)
pos = {(i, j): (i, j) for (i, j) in G.nodes()}
plt.figure()
nx.draw(D, pos=pos, with_labels=True, node_size=1000, node_color="yellow")
edge_labels = {}
for (i, j) in G.edges():
edge_labels[i, j] = f"{D[i][j]['weight']} \n{D[j][i]['weight']}"
nx.draw_networkx_edge_labels(D, pos, edge_labels=edge_labels)
plt.show()
edmonds = nx.tree.Edmonds(D)
sol = edmonds.find_optimum(attr="weight", kind="min", style="arborescence")
plt.figure()
edge_labels = {}
for (i, j) in sol.edges():
edge_labels[i, j] = f"{D[i][j]['weight']}"
nx.draw_networkx_edge_labels(D, pos, edge_labels=edge_labels)
nx.draw(D, pos=pos, width=5, edgelist=sol.edges(), edge_color="orange")
