import networkx as nx
import matplotlib.pyplot as plt
import random
import math
from itertools import islice
from graphillion import GraphSet
n = 5
G = nx.grid_2d_graph(n, n)
pos = {(i, j): (i, j) for (i, j) in G.nodes()}
for (i, j) in G.edges():
G[i][j]["weight"] = random.randint(1, 10)
def path_length(G, path):
_sum = 0
i = path[0]
for count in range(1, len(path)):
j = path[count]
_sum += G[i][j]["weight"]
i = j
return _sum
k = 10
for path in islice(
nx.shortest_simple_paths(G, source=(0, 0), target=(n - 1, n - 1), weight="weight"),
k,
):
print(path, path_length(G, path))
edge_list = []
i = path[0]
for count in range(1, len(path)):
j = path[count]
edge_list.append((i, j))
i = j
plt.figure()
nx.draw(G, pos=pos, with_labels=False, node_size=100)
nx.draw(
G,
pos=pos,
with_labels=False,
node_size=100,
edgelist=edge_list,
edge_color="red",
width=10,
alpha=0.3,
)
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)
plt.show()
無向パス（閉路，森など）の列挙
Graphillionを使うと，無向部分グラフの列挙ができる．パスだけでなく，閉路，木，森，クリーク（完全部分グラフ）なども列挙もできる．
Graphillionについての詳細は，以下のサイトを参照されたい．
https://github.com/takemaru/graphillion
例として，上と同様の格子グラフにランダムな枝長を与えたグラフに対して，「すべて」のパスを数え上げる．パスの総数は $8512$ であることが分かる．
random.seed(1)
n = 5
G = nx.grid_2d_graph(n, n)
pos = {(i, j): (i, j) for (i, j) in G.nodes()}
weight = {}
for (i, j) in G.edges():
weight[i, j] = random.randint(1, 10)
GraphSet.set_universe(G.edges())
paths = GraphSet.paths(terminal1=(0, 0), terminal2=(n - 1, n - 1))
len(paths) # パスの総数
count = 0
for p in paths.min_iter(weight):
count += 1
if count >= 10:
break
print(p)
plt.figure()
nx.draw(G, pos=pos, with_labels=False, node_size=100)
nx.draw(
G,
pos=pos,
with_labels=False,
node_size=100,
edgelist=p,
edge_color="red",
width=10,
alpha=0.3,
)
nx.draw_networkx_edge_labels(G, pos, edge_labels=weight)
plt.show()
count = 0
for p in paths.max_iter(weight):
count += 1
if count >= 2:
break
print(p)
plt.figure()
nx.draw(G, pos=pos, with_labels=False, node_size=100)
nx.draw(
G,
pos=pos,
with_labels=False,
node_size=100,
edgelist=p,
edge_color="red",
width=10,
alpha=0.3,
)
nx.draw_networkx_edge_labels(G, pos, edge_labels=weight)
plt.show()
def distance(x1, y1, x2, y2):
"""distance: euclidean distance between (x1,y1) and (x2,y2)"""
return math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
n = 10
x = dict([(i, 100 * random.random()) for i in range(n)])
y = dict([(i, 100 * random.random()) for i in range(n)])
G = nx.Graph()
c = {}
pos = {}
for i in range(n):
pos[i] = x[i], y[i]
for j in range(n):
if j > i:
c[i, j] = distance(x[i], y[i], x[j], y[j])
G.add_edge(i, j, weight=c[i, j])
edges = list(G.edges())
GraphSet.set_universe(edges)
cycles = GraphSet.cycles()
gs = cycles.graph_size(4)
gs2 = gs.including(0) # include depot
cycle = gs2.choice()
nx.draw(G, pos=pos, with_labels=True, node_color="y")
nx.draw_networkx_edges(G, pos=pos, edgelist=cycle, edge_color="orange", width=5);
cycles = GraphSet.cycles(is_hamilton=True)
print(len(cycles))
for p in cycles.min_iter(c):
hamilton_cycle = p
break
nx.draw(G, pos=pos, with_labels=True, node_color="y")
nx.draw_networkx_edges(
G, pos=pos, edgelist=hamilton_cycle, edge_color="orange", width=5
);
多目的最短路問題
費用と時間の2つの重みをもつグラフに対して，意思決定者にとって便利な複数のパスを示す問題を考える． これは，多目的最適化問題になる．
まず，多目的最適化の基礎と用語について述べる．
以下に定義される $m$個の目的をもつ多目的最適化問題を対象とする．
解の集合 $X$ ならびに $X$ から $m$次元の実数ベクトル全体への写像 $f: X \rightarrow \mathbf{R}^m$ が与えられている． ベクトル $f$ を目的関数ベクトルとよび，その第 $i$ 要素を $f_i$ と書く． ここでは，ベクトル $f$ を「何らかの意味」で最小にする解（の集合）を求めることを目的とする．
2つの目的関数ベクトル$f,g \in \mathbf{R}^m$ に対して， $f$ と $g$ が同じでなく，かつベクトルのすべての要素に対して $f$ の要素が $g$ の要素以下 であるとき，ベクトル $f$ がベクトル $g$ に優越しているとよび，$f \prec g$ と記す．
すなわち，順序 $\prec$ を以下のように定義する． $$ f \prec g \Leftrightarrow f \neq g, f_i \leq g_i \ \ \ \forall i $$ たとえば，2つのベクトル $f=(2,5,4)$ と $g=(2,6,8)$ に対しては，第1要素は同じであるが，第2,3要素に対しては $g$ の 方が大きいので $f \prec g$ である．
2つの解 $x,y$ に対して，$f(x) \prec f(y)$ のとき，解$x$ は解$y$ に優越していると呼ぶ． 以下の条件を満たすとき，$x$ は非劣解(nondominated solution) もしくはPareto最適解(Pareto optimal solution)と呼ばれる． $$ f(y) \prec f(x) を満たす解 y \in X は存在しない $$
多目的最適化問題の目的は，すべての非劣解(Pareto最適解)の集合を求めることである． 非劣解の集合から構成される境界は， 金融工学における株の構成比を決める問題（ポートフォリオ理論） では有効フロンティア(efficient frontier)と呼ばれる． ポートフォリオ理論のように目的関数が凸関数である場合には， 有効フロンティアは凸関数になるが，一般には非劣解を繋いだものは凸になるとは限らない．
非劣解の総数は非常に大きくなる可能性がある． そのため，実際にはすべての非劣解を列挙するのではなく， 意思決定者の好みにあった少数の非劣解を選択して示すことが重要になる．
最も単純なスカラー化は複数の目的関数を適当な比率を乗じて 足し合わせることである．
$m$次元の目的関数ベクトルは，$m$次元ベクトル $\alpha$ を用いてスカラー化できる． 通常，パラメータ $\alpha$ は $$ \sum_{i=1}^m \alpha_i =1 $$ を満たすように正規化しておく．
この $\alpha$ を用いて重み付きの和をとることにより， 以下のような単一の（スカラー化された）目的関数 $f_{\alpha}$ に変換できる． $$ f_{\alpha}(x)= \sum_{i=1}^m \alpha_i f_i(x) $$
これを2目的の最短路問題に適用してみよう．格子グラフの枝に２つの重み（costとtime）を定義して，スカラー化を用いて，有効フロンティアを描画する．
m, n = 100, 100
lb, ub = 1, 100
G = nx.grid_2d_graph(m, n)
for (i, j) in G.edges():
G[i][j]["cost"] = random.randint(lb, ub)
G[i][j]["time"] = 100 / G[i][j]["cost"]
x, y = [], []
for k in range(100):
alpha = 0.01 * k
for (i, j) in G.edges():
G[i][j]["weight"] = alpha * G[i][j]["cost"] + (1 - alpha) * G[i][j]["time"]
pred, distance = nx.dijkstra_predecessor_and_distance(G, source=(0, 0))
# print("minimum cost=", distance[m-1,n-1])
j = (m - 1, n - 1)
cost = time = 0
while i != (0, 0):
i = pred[j][0]
cost += G[i][j]["cost"]
time += G[i][j]["time"]
j = i
# print(cost,time)
x.append(cost)
y.append(time)
plt.plot(x, y, color="green", marker="o")
import collections
def all_simple_paths_graph(G, source, targets, cutoff):
visited = collections.OrderedDict.fromkeys([source])
# visited に現在のパスの情報を追加保存する
# 必要がある．辞書の値に，積んである荷物量や発時刻などを保管しておく．
stack = [iter(G[source])]
while stack:
children = stack[-1] # 後続点の集合をイテレータとして得る
# print("children=", list(children) )
child = next(children, None) # 後続点から1つの点を選択する．なければNone
if child is None: # もう調べる後続点がなければ、スタックを減らす
stack.pop()
visited.popitem()
elif len(visited) < cutoff:
if child in visited: # すでに訪問いている点ならコンティニュー
continue
# ここでvisited[-1] からchildへ移動可能かどうかのcheckを入れる。
print(" from", list(visited)[-1], " to", child)
# 時間枠のcheckし、実行不可能や閾値を超えた待ち時間が生じるならなcontinueなど
if child in targets:
# visitedに現在の点を追加したものを出力
yield list(visited) + [child]
visited[child] = None
if targets - set(visited.keys()):
# expand stack until find all targets
stack.append(iter(G[child]))
# 次の後続点集合をスタックに追加（深さ優先探索）
else:
visited.popitem() # maybe other ways to child
else: # len(visited) == cutoff:
for target in (targets & (set(children) | {child})) - set(visited.keys()):
yield list(visited) + [target]
# visitedに現在の点を追加したものを出力
stack.pop() # カットオフ値に達したら、スタックを減らす
visited.popitem()
G = nx.complete_graph(4)
for path in all_simple_paths_graph(G, source=0, targets={3}, cutoff=len(G) - 1):
print(path)