第k最短路とGraphillionによる無向グラフの列挙と多目的最短路

準備

import networkx as nx
import matplotlib.pyplot as plt
import random
import math
from itertools import islice
from graphillion import GraphSet

第k最短路

networkXにYenの第 $k$ 最短路のアルゴリズムが含まれている.これは,最短路を短い(費用の小さい)順に列挙してくれる.

例として,格子グラフにランダムな枝長を与えたグラフに対して,第10番目までの最短路を列挙し,10番目に短いパスを図示する.

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))
[(0, 0), (1, 0), (1, 1), (1, 2), (2, 2), (3, 2), (3, 3), (3, 4), (4, 4)] 28
[(0, 0), (0, 1), (1, 1), (1, 2), (2, 2), (3, 2), (3, 3), (3, 4), (4, 4)] 29
[(0, 0), (1, 0), (1, 1), (2, 1), (2, 2), (3, 2), (3, 3), (3, 4), (4, 4)] 29
[(0, 0), (0, 1), (1, 1), (2, 1), (2, 2), (3, 2), (3, 3), (3, 4), (4, 4)] 30
[(0, 0), (1, 0), (2, 0), (2, 1), (2, 2), (3, 2), (3, 3), (3, 4), (4, 4)] 32
[(0, 0), (1, 0), (2, 0), (2, 1), (1, 1), (1, 2), (2, 2), (3, 2), (3, 3), (3, 4), (4, 4)] 33
[(0, 0), (1, 0), (1, 1), (1, 2), (2, 2), (3, 2), (4, 2), (4, 3), (4, 4)] 34
[(0, 0), (0, 1), (1, 1), (1, 2), (2, 2), (3, 2), (4, 2), (4, 3), (4, 4)] 35
[(0, 0), (1, 0), (1, 1), (2, 1), (2, 2), (3, 2), (4, 2), (4, 3), (4, 4)] 35
[(0, 0), (1, 0), (1, 1), (1, 2), (2, 2), (3, 2), (3, 3), (4, 3), (4, 4)] 36
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)  # パスの総数
8512

(最短)パスの列挙

得られたpathsオブジェクトのmin_iterメソッドでパス長の短い順に列挙できる.引数として,枝上に定義された重みを返す辞書を与え,10番目のパスまでを列挙し,10番目のパスを図示する.

count = 0
for p in paths.min_iter(weight):
    count += 1
    if count >= 10:
        break
    print(p)
[((0, 0), (1, 0)), ((1, 0), (2, 0)), ((2, 0), (2, 1)), ((2, 1), (2, 2)), ((2, 2), (2, 3)), ((2, 3), (3, 3)), ((3, 3), (4, 3)), ((4, 3), (4, 4))]
[((0, 0), (1, 0)), ((1, 0), (1, 1)), ((1, 1), (1, 2)), ((1, 2), (2, 2)), ((2, 2), (2, 3)), ((2, 3), (3, 3)), ((3, 3), (4, 3)), ((4, 3), (4, 4))]
[((0, 0), (1, 0)), ((1, 0), (1, 1)), ((1, 1), (2, 1)), ((2, 1), (2, 2)), ((2, 2), (2, 3)), ((2, 3), (3, 3)), ((3, 3), (4, 3)), ((4, 3), (4, 4))]
[((0, 0), (1, 0)), ((1, 0), (2, 0)), ((2, 0), (2, 1)), ((2, 1), (2, 2)), ((2, 2), (3, 2)), ((3, 2), (4, 2)), ((4, 2), (4, 3)), ((4, 3), (4, 4))]
[((0, 0), (1, 0)), ((1, 0), (1, 1)), ((1, 1), (1, 2)), ((1, 2), (2, 2)), ((2, 2), (3, 2)), ((3, 2), (4, 2)), ((4, 2), (4, 3)), ((4, 3), (4, 4))]
[((0, 0), (0, 1)), ((0, 1), (1, 1)), ((1, 1), (1, 2)), ((1, 2), (2, 2)), ((2, 2), (2, 3)), ((2, 3), (3, 3)), ((3, 3), (4, 3)), ((4, 3), (4, 4))]
[((0, 0), (1, 0)), ((1, 0), (2, 0)), ((2, 0), (2, 1)), ((2, 1), (2, 2)), ((2, 2), (2, 3)), ((2, 3), (3, 3)), ((3, 3), (3, 4)), ((3, 4), (4, 4))]
[((0, 0), (1, 0)), ((1, 0), (1, 1)), ((1, 1), (1, 2)), ((1, 2), (2, 2)), ((2, 2), (2, 3)), ((2, 3), (3, 3)), ((3, 3), (3, 4)), ((3, 4), (4, 4))]
[((0, 0), (1, 0)), ((0, 1), (0, 2)), ((0, 1), (1, 1)), ((0, 2), (1, 2)), ((1, 0), (1, 1)), ((1, 2), (2, 2)), ((2, 2), (2, 3)), ((2, 3), (3, 3)), ((3, 3), (4, 3)), ((4, 3), (4, 4))]
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()

最長パスの列挙(最長路問題)

max_iterメソッドを使うと,反復を短い順ではなく長い順にすることもできる. これを使うと最長路問題(longeest path problem)を解くことができる.この問題は $NP$-困難であるが, グラフが疎な場合には,列挙によって求めることができる. 例として,格子グラフに対する最長路を求めてみる.

count = 0
for p in paths.max_iter(weight):
    count += 1
    if count >= 2:
        break
    print(p)
[((0, 0), (0, 1)), ((0, 1), (0, 2)), ((0, 2), (0, 3)), ((0, 3), (0, 4)), ((0, 4), (1, 4)), ((1, 0), (1, 1)), ((1, 0), (2, 0)), ((1, 1), (2, 1)), ((1, 2), (1, 3)), ((1, 2), (2, 2)), ((1, 3), (2, 3)), ((1, 4), (2, 4)), ((2, 0), (3, 0)), ((2, 1), (3, 1)), ((2, 2), (3, 2)), ((2, 3), (3, 3)), ((2, 4), (3, 4)), ((3, 0), (4, 0)), ((3, 1), (3, 2)), ((3, 3), (3, 4)), ((4, 0), (4, 1)), ((4, 1), (4, 2)), ((4, 2), (4, 3)), ((4, 3), (4, 4))]
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()

閉路の列挙

Graphillionを使うと,パスだけでなく閉路も列挙できる.例として,平面上に10個の点を配置し,直線距離の長さをもつグラフを考える.

点 $0$ はデポ(トラックの初期位置)であり,4件の顧客(他の点)を巡回して再びデポに戻る閉路を列挙し,そこからchoiceメソッドでランダムに1つの閉路を選んで描画する.

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);

Hamilton閉路の列挙

cyclesの引数is_hamiltonをTrueにすることによって,すべての点をちょうど1回ずつ通過する閉路(Hamilton閉路)を列挙できる. Hamilton閉路は,全部で181440あることが分かる.また,最も移動距離の短いHamilton閉路を描画する.これは,巡回セールスマン問題の最適解になる.

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
);
181440

多目的最短路問題

費用と時間の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目的の最短路問題に適用してみよう.格子グラフの枝に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)
 from 0  to 1
 from 1  to 2
[0, 1, 2, 3]
 from 1  to 3
[0, 1, 3]
 from 0  to 2
 from 2  to 1
[0, 2, 1, 3]
 from 2  to 3
[0, 2, 3]
 from 0  to 3
[0, 3]