import networkx as nx
from networkx import DiGraph
import matplotlib.pyplot as plt
from numpy import array
import numpy as np
from heapq import heappush, heappop
from itertools import count
import requests
import pprint
from heapdict import heapdict
from cspy import BiDirectional
最短路問題
ここで考える最短路問題は多くの応用をもつ基本的なネットワーク最適化問題である．
（1対1の）最短路問題（shortest path problem）
$n$個の点から構成される点集合 $V$，$m$ 本の枝から構成される枝集合 $E$， $V$ および $E$ からなるグラフ $G=(V,E)$，枝上に定義される非負の費用関数 $c: E \rightarrow \mathbf{R}_+$，始点 $s \in V$ および終点 $t \in V$ が与えられたとき， 始点 $s$ から終点 $t$ までのパスで，パスに含まれる枝の費用の合計が最小のものを求めよ．
最短路問題の目的である「パスに含まれる枝の費用の合計が最小のもの」を最短路（shortest path）とよぶ．
枝 $e$ がパスに含まれるとき $1$ になる変数を使うと，最短路問題は以下のように定式化できる．
$$ \begin{array}{c l} minimize & \sum_{e \in E} c_{e} x_{e} \\ s.t. & \sum_{w: wv \in E} x_{wv} - \sum_{w: vw \in E} x_{vw} = \left\{ \begin{array}{c l} -1 & v=s \\ 0 & \forall v \in V \setminus \{s,t\} \\ 1 & v=t \end{array} \right. \\ & x_e \in \{ 0,1 \} \ \ \ \forall e \in E \end{array} $$$0$-$1$変数 $x$ を $0$ 以上に緩和した線形最適化を解くことによって，整数解が得られることが知られている．
実用的なアルゴリズムとして，Dijkstra法（Dijkstra's method）がある。 Dijkstra法は，最短路問題に対する最も有名なアルゴリズムであり， ほとんどのテキストでは最初にとりあげられるアルゴリズムである．
networkXには様々な最短路アルゴリズムのための関数が準備されているが，単一始点から他のすべての点に対する最短路を求めるものとしては，dijkstra_predecessor_and_distanceが使いやすい． この関数の引数と返値は，以下の通りである．
引数：
- source: 始点を表すラベル
- cutoff: この値を超える費用をもつパスは探索しない（カットオフする）
- weight: 枝(u,v)の費用がG.edges[u,v][weight]に格納されていると仮定する
返値： 以下の2つの辞書のタプル
- pred: 点のラベルをキーとし，先行点のリストを値とした辞書
- distance: 点のラベルをキーとし，その点までの最短路の費用を値とした辞書
networkXを使うと，有向・無向の両方のグラフに対する最短路問題を簡単に解くことができる． 例として，小さな問題例（有向グラフを仮定）で求解してみよう．
D = nx.DiGraph()
D.add_weighted_edges_from(
[("s", 1, 10), ("s", 2, 5), (2, 1, 3), (1, "t", 4), (2, 3, 2), (3, "t", 6)]
)
pred, distance = nx.dijkstra_predecessor_and_distance(D, source="s")
print(pred)
print(distance)
edge_labels = {}
pos = nx.kamada_kawai_layout(D)
for (i, j) in D.edges():
edge_labels[i, j] = f"{ D[i][j]['weight'] }"
nx.draw_networkx_edge_labels(D, pos=pos, edge_labels=edge_labels)
nx.draw(
D, node_size=300, pos=pos, node_color="y", edge_color="g", width=3, with_labels=True
)
plt.show()
始点と終点が決まっている場合には， A*探索を使うと（多少は）高速になると言われているが，問題例によっては遅くなる場合もある．
ランダムな格子グラフで試してみる． $m \times n$ の格子グラフの各枝に， $[lb,ub]$ の一様整数乱数の費用を設定する． 端にある点 $(0,0)$ から $(m-1,n-1)$ までの最短路を，色々なアルゴリズムで求めてみよう．
m, n = 1000, 1000
lb, ub = 1, 300
G = nx.grid_2d_graph(m, n)
for (i, j) in G.edges():
G[i][j]["weight"] = random.randint(lb, ub)
%%time
pred, distance = nx.dijkstra_predecessor_and_distance(G, source=(0,0))
標準的なdijkstra_predecessor_and_distanceを使うと， $1000 \times 1000$ の問題例が8秒程度で解ける． 点$(m-1,n-1)$ までの最小費用は辞書distanceに格納されている．また，最短路自身は，パスの前の点を保持する辞書predに格納されているので，以下のように復元できる． $1000 \times 1000$ は大きすぎるので， $10 \times 10$ で計算してみる．
m, n = 10, 10
lb, ub = 1, 300
G = nx.grid_2d_graph(m, n)
for (i, j) in G.edges():
G[i][j]["weight"] = random.randint(lb, ub)
pred, distance = nx.dijkstra_predecessor_and_distance(G, source=(0, 0))
print("minimum cost=", distance[m - 1, n - 1])
i = (m - 1, n - 1)
path = []
while i != (0, 0):
path.append(i)
i = pred[i][0]
path.reverse()
print("optimal path =", path)
始点と終点が決まっている場合には，$A^{*}$探索が使える．この場合には，あまり高速化しない．
%%time
path = nx.astar_path(G, (0,0), (m-1,n-1))
始点と終点の両方から探索するという関数もあったので試してみる． 高速化どころか，5倍以上の計算時間がかかるので，これは使うべきではない．
%%time
length, path = nx.bidirectional_dijkstra(G,(0,0), (m-1,n-1))
1対全の最短路問題
Dijkstra法は，1つの始点から他のすべての点に対する最短路を一度に求めることができる． これは，1対全の最短路問題と呼ばれる．
$n$個の点から構成される点集合 $V$，$m$ 本の枝から構成される枝集合 $E$， $V$ および $E$ からなる有向グラフ $G=(V,E)$，枝上に定義される非負の費用関数 $c: E \rightarrow \mathbf{R}_+$，始点 $s \in V$ が与えられたとき， 始点 $s$ から他のすべての点までの最小費用のパスを求めよ．
dijkstra_predecessor_and_distanceを用いて得られるパスの前の点を保持する辞書predは，すべての最短路の情報を保持している． 直前の点 $pred[v]$ から $v$ に有向枝を張ることによって，始点 $s$ を根にした有向木が得られる． この有向木を最短路木（shortest path tree）とよぶ．
$10 \times 10$ の格子グラフの最短路木（左下の点が始点）を描画してみる．
D = nx.DiGraph()
for i in pred:
if len(pred[i]) >= 1:
D.add_edge(pred[i][0], i)
pos = {(i, j): (i, j) for (i, j) in G.nodes()}
plt.figure()
nx.draw(D, pos=pos, width=2, edge_color="b", node_size=10)
plt.show()
枝の費用が負のものがある場合
Dijkstra法は効率的な解法であるが，枝の費用が負の場合には使えない．しかも，負の閉路があると（そこをくるくる回ると永遠に儲かるので）アルゴリズムが破綻する． そのような場合にはBellman-Ford法を使う． この解法は，動的最適化の一種であり，以下のBellman等式（Bellman's equation）に基づく．
$$ \begin{array}{lll} y_w = & \min_{v: vw \in E} \{ y_v +c_{vw} \} & \forall v \in V \setminus \{s\} \\ y_s = & 0 & \end{array} $$使用する関数は bellman_ford_predecessor_and_distanceである． どこかに負の閉路（そこを巡回すると費用はマイナス無限大に発散する）があると，NetworkXUnbounded例外を発生して終了する． 引数としてheuristicが追加されており，それをTrueに設定すると，負の閉路がある場合に高速に計算できる． goldberg_radzik関数も同様の手法であるが，より高速である．
$100 \times 100$ の有向格子グラフを生成し，枝の費用を $[-100,100]$ の一様整数乱数として解いてみる．（有向グラフにしたのは，無向グラフだと負の費用をもつ枝が存在するだけで，負の閉路になるためである．）
m, n = 10, 10
lb, ub = -100, 100
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)
%%time
pred, distance ={},{}
try:
pred, distance = nx.bellman_ford_predecessor_and_distance(D, source=(0,0), heuristic=False)
except nx.NetworkXUnbounded as e:
print(e)
%%time
try:
pred, distance = nx.goldberg_radzik(D, source=(0,0))
except nx.NetworkXUnbounded as e:
print(e)
m, n = 50, 50
lb, ub = 1, 100
G = nx.grid_2d_graph(m, n)
for (i, j) in G.edges():
G[i][j]["weight"] = random.randint(lb, ub)
%%time
try:
paths = nx.johnson(G)
except nx.NetworkXUnbounded as e:
print(e)
%%time
nx.floyd_warshall_numpy(G)
道路ネットワークの最短路と前処理による高速化
実際の道路ネットワークでの最短路の計算に，通常のDijkstra法を使うことは推奨できない． 遠い2点間の最短路の計算には，結構な時間がかかるからだ．
実際には適当な前処理（これには時間がかかる）をしておき，実際の始点と終点が与えられたときには，前処理の情報を用いて超高速に最短路を計算するというテクニックが使われる．
オープンソースの地図であるOpenStreetMapに，Highway Hierarchyという前処理法を適用したプログラムとしてOSRM http://project-osrm.org/ がある．
Highway Hierarchyにおける前処理を簡単に解説する．
まず高速ネットワークの階層の概念を導入する． ここで言う高速とは，有料道路のことではなく，遠くの地点に行くときに用いられる疎なネットワークのことである．
始点 $s$ から他の全点に対する最短路をDijkstra法によって求める際， 優先キューから出される順番をDijkstraランクと定義する． 点 $v$ のDijkstraランクを $r_s(v)$ と記す． 始点 $s$ のDijkstraランク $r_s(s)$ は $0$ であり， $s$ に最も近い点 $v$ のDijkstraランク $r_s(v)$ は $1$ である．
$H$ 番目の（言い換えれば $H=r_s(v)$ を満たす） 点 $v$ への距離以下の点集合を $s$ の $H$ 近傍とよび， $N_H(s)$ と記す． パラメータ $H$ を与えたとき，第1階層の高速ネットワーク $G_1=(V_1,E_1)$ は， オリジナルのグラフ $G=(V,E)$ から以下のように作成する． ある始点 $s$ からある終点 $t$ への最短路に含まれる枝 $(v,w)$ において， $w \not\in N_H(s)$ かつ $v \not\in N_H(t)$ を満たすものをすべて合わせた集合を $E_1$ とする． $E_1$ は遠くの地点間に行くときに使用する枝の集合であり，高速遠路と考えれば良い．
次に，$E_1$ から構成されるグラフ $G_1=(V_1,E_1)$ において，点の次数が $2$ 以上のもの（インターチェンジと考えれば良い）を抽出し，さらに次数 $2$ の点を除いて， その点に隣接する2つの点の間に枝をはる操作を繰り返すことによってを縮約高速ネットワーク $G'_1$ を生成する． 上の操作を第2，第3階層と順に繰り返すことによって，多階層の高速ネットワークを得る．
このネットワークを使って，始点 $s$ から終点 $t$ までの最短路を計算する際には，始点 $s$ と終点 $t$ からの両方向探索を行う． また，通常のDijkstra法で上の階層に移動できるときには， 必ず上の階層を優先して探索することによって，高速に最短路を求めることができる．
この手法を使うと，日本全国の地図で（メモリが十分にあれば）数時間で前処理ができ，地点間の最短路が高速に計算できるようになる．
例として，東京海洋大学から東京駅までの移動距離と移動時間を求めてみる．
requestsでOSRMのサーバーにクエリを渡すと，結果が返される． responseをJSON化すると， 3675メートル，324秒（車で）かかることが分かる．
response = requests.get(
"http://router.project-osrm.org/route/v1/driving/139.792429,35.667864;139.768525,35.681010"
)
pprint.pprint(response.json())
単なる最短路ではなく，与えた複数の座標間の全対全の最短路長と移動時間のテーブルも（高速に）計算できる．
response = requests.get(
"http://router.project-osrm.org/table/v1/driving/13.388860,52.517037;13.397634,52.529407;13.428555,52.523219"
+ "?annotations=distance,duration"
)
pprint.pprint(response.json())
ベンチマーク問題例の読み込みとDial法
最短路問題のベンチマーク問題例（道路ネットワーク）が，以下のサイトで入手できる．
http://users.diag.uniroma1.it/challenge9/download.shtml
以下の関数を用いてデータを読み込んで，道路ネットワークに強いDial法を実装して解いてみる．
def ReadDimacs(filename):
"""Read the data from DIMACS shortest path format
that can be downloaded from DIMACS site"""
f = open(filename + ".co", "r")
f2 = open(filename + ".gr", "r")
# read data from file
line = f.readline()
x_cord = []
y_cord = []
G = nx.Graph()
G.position = {}
while line.find("v 1") == -1:
line = f.readline()
(dum, i, x, y) = line.split()
x_cord.append(float(x))
y_cord.append(float(y))
while 1:
line = f.readline()
if line.find("v") == 0:
(dum, i, x, y) = line.split()
x_cord.append(float(x))
y_cord.append(float(y))
else:
break
# edge info
edge_info = []
line = f2.readline()
for i in range(6):
line = f2.readline()
while 1:
line = f2.readline()
if line.find("a") == 0:
(dum, i, j, length) = line.split()
edge_info.append((int(i) - 1, int(j) - 1, int(length)))
else:
break
minx = float(min(x_cord))
maxx = float(max(x_cord))
xwidth = maxx - minx + 1
miny = float(min(y_cord))
maxy = float(max(y_cord))
ywidth = maxy - miny + 1
# add nodes after scaling coordinates
n = len(x_cord)
for i in range(n):
G.add_node(i)
x = (x_cord[i] - minx) / xwidth
y = (y_cord[i] - miny) / ywidth
G.position[i] = (x, y)
print("number of nodes=", n)
# add edges
m = len(edge_info)
C = 0
for (i, j, length) in edge_info:
G.add_edge(i, j, weight=length)
if C < length:
C = length
C += 1
print("number of edges=", m)
f.close()
f2.close()
return G, n, m, C
def DialSub(n, G, settled, reached, dist, prev, C, source, sink):
"""Dial method subroutine"""
circular = [[-1] for i in range(C)] # -1 represents the list is empty
circular[0].append(source)
reached[source] = 1
prev[source] = -1
current = 0
while 1:
first = current
while 1:
v = circular[current][-1] # v is the node with minimum potential
if v >= 0:
break
# the current list is empty, so we increment it
current += 1
current = current % C
if current == first:
# the circular list is empty, so we quit
v = -1
break
if v == -1: # the circular list is empty, so we quit
break
circular[current].pop()
settled[v] = 1
if v == sink:
break
for w in G.neighbors(v): # scan operation
if settled[w] == 0:
temp = dist[v] + G[v][w]["weight"]
if reached[w] == 0: # w is not in the circular list
reached[w] = 1
dist[w] = temp
prev[w] = v
circular[temp % C].append(w)
else: # w is in the circular list
if temp < dist[w]:
circular[dist[w] % C].remove(w)
circular[temp % C].append(w)
dist[w] = temp
prev[w] = v
filename = "USA-road-d.NY"
G, n, m, C = ReadDimacs(filename)
print("number of nodes in the largest connected componet=", n, "max length+1=", C)
source = 0
sink = 50000 # if sink =-1 then find all shortest paths from source to other nodes
settled = np.zeros(
n, np.int0
) # =1 if node is settled (or has an eternal label or is scanned)
reached = np.zeros(n, np.int0) # =1 if node is heap or list (or has a finite potential)
dist = np.zeros(n, np.int0)
prev = np.zeros(n, np.int0)
settled = [0 for i in range(n)]
reached = [0 for i in range(n)]
dist = [0 for i in range(n)]
prev = [0 for i in range(n)]
%%time
DialSub(n,G,settled,reached,dist,prev,C,source,sink)
%%time
length=nx.dijkstra_path_length(G,source,sink)
dist[50000], length
networkXのパス長だけを算出する関数 dijkstra_path_lengthを用いた方が多少速い．以下で，結果を描画する．
path = [sink]
w = sink
while 1:
v = prev[w]
path.append(v)
if v == source:
break
w = v
path.reverse()
print("Path Length=", dist[sink])
print("Optimal Path", path)
G.clear()
nx.draw(
G,
G.position,
node_size=1000 / n,
with_labels=None,
width=1,
node_color="b",
edge_color="g",
)
nx.add_path(G, path)
nx.draw(
G,
G.position,
node_size=1000 / n,
with_labels=None,
width=3,
node_color="b",
edge_color="r",
)
%%time
source=0
target =50000 #if sink =-1 then find all shortest paths from source to other nodes
weight='weight'
get_weight = lambda u, v, data: data.get(weight, 1)
G_succ = G.succ if G.is_directed() else G.adj
pred = {source: []}
paths={source: [source]}
push = heappush
pop = heappop
dist = {} # dictionary of final distances
seen = {source: 0}
c = count()
fringe = [] # use heapq with (distance,label) tuples
push(fringe, (0, next(c), source))
while fringe:
(d, _, v) = pop(fringe)
if v in dist:
continue # already searched this node.
dist[v] = d
if v == target:
break
for u, e in G_succ[v].items():
cost = get_weight(v, u, e)
if cost is None:
continue
vu_dist = dist[v] + get_weight(v, u, e)
if u in dist:
if vu_dist < dist[u]:
raise ValueError('Contradictory paths found:',
'negative weights?')
elif u not in seen or vu_dist < seen[u]:
seen[u] = vu_dist
push(fringe, (vu_dist, next(c), u))
if paths is not None:
paths[u] = paths[v] + [u]
if pred is not None:
pred[u] = [v]
elif vu_dist == seen[u]:
if pred is not None:
pred[u].append(v)
dist[50000]
%%time
source=0
target =50000 #if sink =-1 then find all shortest paths from source to other nodes
weight='weight'
get_weight = lambda u, v, data: data.get(weight, 1)
G_succ = G.succ if G.is_directed() else G.adj
pred = {source: []}
paths={source: [source]}
dist = {} # dictionary of final distances
hd = heapdict()
seen = {source: 0}
c = count()
hd[source] = 0
dist
while hd:
(v, d) = hd.popitem()
# decrease-keyの場合には以下が不要
#if v in dist:
# continue # already searched this node.
dist[v] = d
if v == target:
break
for u, e in G_succ[v].items():
cost = get_weight(v, u, e)
if cost is None:
continue
vu_dist = dist[v] + get_weight(v, u, e)
if u in dist:
if vu_dist < dist[u]:
raise ValueError('Contradictory paths found:',
'negative weights?')
elif u not in seen or vu_dist < seen[u]:
seen[u] = vu_dist
#push(fringe, (vu_dist, next(c), u))
hd[u] = vu_dist
if paths is not None:
paths[u] = paths[v] + [u]
if pred is not None:
pred[u] = [v]
elif vu_dist == seen[u]:
if pred is not None:
pred[u].append(v)
dist[50000]
時刻依存の移動時間をもつ最短路問題
道路ネットワークに対する最短路（最短時間路）を求めたい場合には，道路の混雑などで移動時間が時刻によって変化する場合がある． ここでは，そのような拡張が簡単に解けることを示す．
地点 $i,j$ 間の移動時間が時刻 $t$ によって $c_{ij}(t)$ という関数で与えられている場合を考える．
任意の2つの時刻 $t \leq t'$ に対して，$t+c_{ij}(t) \leq t' + c_{ij}(t')$ が成立するとき，移動時間は先入れ・先出し（First In First Out: FIFO）であるとよばれる．
点の上で「待つ」ことを許す場合には，FIFOになるように移動時間を変形できる．待ちを許さない場合には，$NP$-困難になる．
FIFOの場合には，通常のDijkstra法において，枝の移動時間を計算する部分を，到着時刻関数に変更するだけで動作する．
interval = [0, 4, 8, 10] # time interval
v = [10, 20, 30, 20] # velocity
distance = 10
interval, b, a = arrival_func(interval, v, distance)
print(interval, b, a)
for t in range(interval[-1] + 5):
print(t, arrival(t, interval, b, a))
plot_arrival(interval, b, a)
G = nx.Graph()
G.add_weighted_edges_from(
[("s", 1, 50), ("s", 2, 30), (2, 1, 20), (1, "t", 80), (2, 3, 50), (3, "t", 60)]
)
interval = [0, 2, 4, 7] # time interval
velocity = [16, 20, 35, 20] # velocity
for (i, j) in G.edges():
G[i][j]["piecewise"] = arrival_func(interval, velocity, distance=G[i][j]["weight"])
for (i, j) in G.edges:
interval, b, a = G[i][j]["piecewise"]
plot_arrival(interval, b, a)
上で準備した例題に対して，FIFOを仮定した時刻依存最短路問題を解く．
source = "s"
target = "t"
G_succ = G.succ if G.is_directed() else G.adj
pred = {source: []}
paths = {source: [source]}
push = heappush
pop = heappop
dist = {} # dictionary of final distances
seen = {source: 0}
c = count()
fringe = [] # use heapq with (distance,label) tuples
push(fringe, (0, next(c), source))
while fringe:
(d, _, v) = pop(fringe)
if v in dist:
continue # already searched this node.
dist[v] = d
if v == target:
break
for u, e in G_succ[v].items():
interval, b, a = G[v][u]["piecewise"]
vu_dist = arrival(dist[v], interval, b, a)
if u in dist:
if vu_dist < dist[u]:
raise ValueError("Contradictory paths found:", "negative weights?")
elif u not in seen or vu_dist < seen[u]:
seen[u] = vu_dist
push(fringe, (vu_dist, next(c), u))
if paths is not None:
paths[u] = paths[v] + [u]
if pred is not None:
pred[u] = [v]
elif vu_dist == seen[u]:
if pred is not None:
pred[u].append(v)
print(pred)
print(dist)
資源制約付き最短路問題
ここでは，運搬車スジューリング問題や乗務員スケジューリング問題を列生成法を用いて解く際の子問題としてあらわれる資源制約付き最短路問題について考える． 列生成法では，主問題で線形最適化問題を解き，子問題では主問題の最適双対変数の情報を点の価値とみなした制約付きのパスを求める必要がある． 資源制約付き最短路問題は，付加制約を一般化した資源としてモデル化し，価値を負の費用として最短路を求める．
資源制約付き最短路問題(resource constrained shortest path problem)とは，以下の仮定をもつ問題である．
始点から終点までのパスを求める． これは人や運搬車の経路を抽象化したものである． 実際には，初等パス（同じ枝を2度以上通過しないパス）である必要があるが， この条件を外した（すなわち，同じ枝を2度以上通過することを許した）パスを考える場合もある．
資源の集合 $R$ が与えられており， 終点における資源量 $r$ の下限 $RLB_r$ と上限 $RUB_r$ が与えられている． これは，時刻，重量，容量などに関する諸制約を抽象化したものである．
枝 $(i,j)$ をパスが通過した際には，資源 $r$ が $t_{ij}^r$ だけ増加するものとする． 荷物の積み降ろしなどで重量や容量が減少する場合や，タスクを処理することによる利益を扱う場合には，対応する $t_{ij}^r$ を負の値に設定すれば良い．
基本モデルにおいては，点 $i$ から点 $j$ に移動したとき，資源 $r$ の量は定数 $t_{ij}^{r}$ だけ大きくなると仮定する． すなわち，点 $i$ における資源 $r$ の量を $T_i^r$ としたとき， 点 $i$ から点 $j$ のへの移動の際には， 以下の関係が成り立つ． $$ T_j^r =T_i^r + t_{ij}^r $$
しかし実際には，点 $i$ における幾つかの資源の量に依存して点 $j$ における資源量が決められる場合がある． ここでは，点 $i$ における資源量を表すベクトル $T_i$ を入力したときに点 $j$ における資源 $r$ の量を返す関数 資源拡張関数 (resource extension function） $f_{ij}^r: \mathbf{R}^{|R|} \rightarrow \mathbf{R}$ を 導入することによってモデルの一般化を行う．
資源拡張関数は，以下の2つの条件を満たす必要がある．
- 最初の資源は単調増加（たとえば使用した枝の本数）
- 逆写像が計算可能
以下のcspyパッケージを用いる．
https://github.com/torressa/cspy
使用法は以下の通り．
入力グラフは networkX の有向グラフ DiGraph のインスタンスとする．
入力グラフは， 資源数を表す n_res属性（グラフの生成時に指定する）をもたなければならない．
入力グラフは，単一の始点 'Source' と終点 'Sink'をもたなければならない．
枝上の属性に $t_{ij}^r$ を表すNumPyの配列型 res_cost と枝の費用を表す特別な資源を表す weightの属性をもたなければならない．
最適化は様々なクラスで行われるが，クラスを生成する際に，終点における資源の上限を表すリスト max_res と下限を表すリスト min_res を設定する．
使用できる最適化手法のクラスは，以下の通り．
BiDirectional: 双方向型の厳密解法
Tabu: タブーサーチ
GreedyElim: 貪欲削減法
GRASP: 貪欲ランダム化適応型探索法 (greedy randomized adaptive search procedure)
PSOLGENT: 蟻コロニー法の亜種 (Particle Swarm Optimization with combined Local and Global Expanding Neighbourhood Topology: PSOLGENT)
使用法は，クラスを生成（elementaryをTrueにすると初等パスを探索）し， runメソッドで求解する．パスはpath属性，総費用はtotal_cost属性，使用した資源量はcomsumed_resources属性で得ることができる．
max_res, min_res = [4, 20], [1, 0]
G = DiGraph(directed=True, n_res=2)
G.add_edge("Source", "A", res_cost=[1, 2], weight=0)
G.add_edge("A", "B", res_cost=[1, 0.3], weight=0)
G.add_edge("A", "C", res_cost=[1, 0.1], weight=0)
G.add_edge("B", "C", res_cost=[1, 3], weight=-10)
G.add_edge("B", "Sink", res_cost=[1, 2], weight=10)
G.add_edge("C", "Sink", res_cost=[1, 10], weight=0)
for (i, j) in G.edges():
G[i][j]["res_cost"] = np.array(G[i][j]["res_cost"])
bidirec = BiDirectional(G, max_res, min_res, elementary=True)
bidirec.run()
print("Path=", bidirec.path)
print("Cost=", bidirec.total_cost)
print("Resource=",bidirec.consumed_resources)
nx.draw(G, with_labels=True, node_color="y")
BOARD_ROWS = 3
BOARD_COLS = 4
WIN_STATE = (0, 3)
LOSE_STATE = (1, 3)
START = (2, 0)
DETERMINISTIC = True
MOVE_PROB = [0.8, 0.9, 1.0] # 確率的移動のときの確率の閾値 (UPのときUP=0.8,LEFU=0.1,RIGHT=0.1)
class State:
"""
状態を定義するクラス
"""
def __init__(self, state=START):
self.board = np.zeros([BOARD_ROWS, BOARD_COLS])
self.board[1, 1] = -1 # 壁
self.state = state # 現在の位置
self.isEnd = False # 終了か否か
self.determine = DETERMINISTIC # 確定的ならTrue
def giveReward(self):
if self.state == WIN_STATE:
return 1
elif self.state == LOSE_STATE:
return -1
else:
return 0
def isEndFunc(self):
if (self.state == WIN_STATE) or (self.state == LOSE_STATE):
self.isEnd = True
def nxtPosition(self, action):
"""
action: up, down, left, right
次の位置を返す
return next position
"""
if self.determine:
if action == "up":
nxtState = (self.state[0] - 1, self.state[1])
elif action == "down":
nxtState = (self.state[0] + 1, self.state[1])
elif action == "left":
nxtState = (self.state[0], self.state[1] - 1)
else:
nxtState = (self.state[0], self.state[1] + 1)
# 壁や外側ではない
if (nxtState[0] >= 0) and (nxtState[0] <= (BOARD_ROWS - 1)):
if (nxtState[1] >= 0) and (nxtState[1] <= (BOARD_COLS - 1)):
if nxtState != (1, 1): # 壁
return nxtState
# 違法な移動なのでそのまま
return self.state
else: # 確率的推移
rnd = np.random.random() # [0,1) unoform random
if action == "up":
if rnd < MOVE_PROB[0]:
nxtState = (self.state[0] - 1, self.state[1]) # up
elif rnd < MOVE_PROB[1]:
nxtState = (self.state[0], self.state[1] - 1) # left
else:
nxtState = (self.state[0], self.state[1] + 1) # right
elif action == "down":
if rnd < MOVE_PROB[0]:
nxtState = (self.state[0] + 1, self.state[1])
elif rnd < MOVE_PROB[1]:
nxtState = (self.state[0], self.state[1] + 1) # right
else:
nxtState = (self.state[0], self.state[1] - 1) # left
elif action == "left":
if rnd < MOVE_PROB[0]:
nxtState = (self.state[0], self.state[1] - 1)
elif rnd < MOVE_PROB[1]:
nxtState = (self.state[0] - 1, self.state[1]) # up
else:
nxtState = (self.state[0] + 1, self.state[1])
else: # right
if rnd < MOVE_PROB[0]:
nxtState = (self.state[0], self.state[1] + 1)
elif rnd < MOVE_PROB[1]:
nxtState = (self.state[0] + 1, self.state[1]) # down
else:
nxtState = (self.state[0] - 1, self.state[1]) # up
# 壁や外側ではない
if (nxtState[0] >= 0) and (nxtState[0] <= (BOARD_ROWS - 1)):
if (nxtState[1] >= 0) and (nxtState[1] <= (BOARD_COLS - 1)):
if nxtState != (1, 1): # 壁
return nxtState
# 違法な移動なのでそのまま
return self.state
def showBoard(self):
"""
盤面を返す
"""
self.board[self.state] = 1
for i in range(0, BOARD_ROWS):
print("-----------------")
out = "| "
for j in range(0, BOARD_COLS):
if self.board[i, j] == 1:
token = "*"
if self.board[i, j] == -1:
token = "z"
if self.board[i, j] == 0:
token = "0"
out += token + " | "
print(out)
print("-----------------")
# Agent of player
class Agent:
def __init__(self):
self.states = []
self.actions = ["up", "down", "left", "right"]
self.State = State()
self.lr = 0.2 #
self.exp_rate = 0.3 # ランダムに選択する確率
# initial state reward
self.state_values = {}
for i in range(BOARD_ROWS):
for j in range(BOARD_COLS):
self.state_values[(i, j)] = 0.0 # set initial value to 0
def chooseAction(self):
# choose action with most expected value
mx_nxt_reward = 0
action = ""
if np.random.uniform(0, 1) <= self.exp_rate:
action = np.random.choice(self.actions)
else:
# 貪欲法
for a in self.actions:
# 移動が確定的のとき
nxt_reward = self.state_values[self.State.nxtPosition(a)]
if nxt_reward >= mx_nxt_reward:
action = a
mx_nxt_reward = nxt_reward
return action
def takeAction(self, action):
position = self.State.nxtPosition(action)
return State(state=position)
def reset(self):
self.states = []
self.State = State()
def play(self, rounds=10, LOG=False):
i = 0
while i < rounds:
# to the end of game back propagate reward
if LOG:
self.State.showBoard()
if self.State.isEnd:
# back propagate
reward = self.State.giveReward()
# explicitly assign end state to reward values
self.state_values[self.State.state] = reward # this is optional
if LOG:
print("Game End Reward", reward)
for s in reversed(self.states):
reward = self.state_values[s] + self.lr * (
reward - self.state_values[s]
)
self.state_values[s] = round(reward, 3)
self.reset()
i += 1
if LOG:
self.showValues()
else:
action = self.chooseAction()
# append trace
self.states.append(self.State.nxtPosition(action))
if LOG:
print(
"current position {} action {}".format(self.State.state, action)
)
# by taking the action, it reaches the next state
self.State = self.takeAction(action)
# mark is end
self.State.isEndFunc()
if LOG:
print("nxt state", self.State.state)
print("---------------------")
def showValues(self):
for i in range(0, BOARD_ROWS):
print("----------------------------------")
out = "| "
for j in range(0, BOARD_COLS):
out += str(self.state_values[(i, j)]).ljust(6) + " | "
print(out)
print("----------------------------------")
ag = Agent()
ag.play(5000, LOG=False)
ag.showValues()