札幌から那覇までの最短経路を求める（グラフ理論）

グラフ理論の勉強として、札幌から那覇までの最短経路を求める。条件としては以下の通り。

都市間の距離は直線距離を使用する。
各都市は最も近い3都市と接続する。

今回の計算は以下のページを参考にした。

まず初めに都市名、緯度、経度を取得する。（コードはほぼコピペ）

url = 'https://raw.githubusercontent.com/maskot1977/ipython_notebook/master/toydata/location.txt'

import urllib.request
import matplotlib.pyplot as plt
import math
import numpy as np

urllib.request.urlretrieve(url, 'location.txt')


col1 = [] # ０列目の数字を格納する予定のリスト
col2 = [] # １列目の数字を格納する予定のリスト
col3 = [] # ２列目の数字を格納する予定のリスト
for i, line in enumerate(open('location.txt')): # ファイルを開いて一行一行読み込む
    if i == 0: # ０番目の行の場合
        continue # 次の行に行く
    c = line.split(",") # 行をコンマで分割したものをcというリストに入れる
    col1.append(c[0]) # ０列目の単語col1に入れる
    col2.append(float(c[1])) # １列目の単語を実数に変換してcol2に入れる
    col3.append(float(c[2])) # ２列目の単語を実数に変換してcol3に入れる

以下のコードでプロット。

plt.figure(figsize=(10, 8))
plt.scatter(col3, col2, alpha=0.5)
plt.title("Prefectural capitals in Japan")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.grid(True)
for city, x, y in zip(col1, col3, col2):
    plt.plot(x, y)
    #plt.text(x, y, city, alpha=0.5, size=12)
plt.show()

以下のコードで都市間を接続する。各都市名をKeyとして、接続されている都市をリストとして持つ辞書(d)を作成する。注意点として、一方から接続されている場合、もう片方からも接続するようにする。これを行わないと例えば沖縄から接続はあるものの、沖縄へは接続がないため、沖縄にわたることができない。

def distance(x1, y1, x2, y2):
    return math.sqrt((x1 - x2)**2 + (y1 - y2)**2)

d = {}
tmp_dist_l = []
for city1, x1, y1 in zip(col1, col3, col2):
    tmp_dist = np.array([])
    for city2, x2, y2 in zip(col1, col3, col2):
        tmp_dist = np.append(tmp_dist, distance(x1, y1, x2, y2))
    d[city1] = {}
    d[city1]["idx"] = []
    tmp_dist_l.append(tmp_dist)
    # 隣接3都市を追加する
    for i in range(3):
        idx = np.argsort(tmp_dist)[i+1]
        d[city1]["idx"].append(idx)

# 一方から接続されている場合、もう片方からも接続
for idx, k in enumerate(d.keys()):
    for i in d[k]["idx"]:
        if not idx in d[col1[i]]["idx"]:
            d[col1[i]]["idx"].append(idx)

# プロット
plt.figure(figsize=(12, 8))
plt.title("Prefectural capitals in Japan")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.grid(True)
for city1, x1, y1 in zip(col1, col3, col2):
    for i in d[city1]["idx"]:
        plt.plot([x1, col3[i]], [y1, col2[i]], 'k-', alpha=0.3, lw=2)
plt.scatter(col3, col2)
plt.show()

うまくいけば以下のような画像が出来上がる。

本題の経路探索を行う。探索方法は幅優先探索で、アルゴリズムのイメージは次の図のような感じ。

これをスタートからゴールまで行う。コードは以下の通り。まずは札幌から、中間あたりの岐阜の経路を探索する。経路探索は候補が1000個見つかった時点で打ち切っている。（コードが早く動くようになっていないのと、そもそもPythonなので、全経路探索は難しい。limitを変更すればいいので、時間とメモリに余裕がある方はぜひ。）

start = 'Sapporo'
goal = 'Gifu'

def way_search(start, goal, print_log=True, limit=1000):
    ans = []
    candidate = [[col1.index(start), i] for i in d[start]["idx"]]
    offset = 0
    while len(candidate) > 0:
        if len(ans) > limit:
            break
        print("next", len(candidate), len(ans))
        for i, arr in enumerate(candidate):
            for ii in d[col1[arr[-1]]]["idx"]:
                # 同じ都市を通る場合、スキップ
                if ii in candidate[i] or col1[ii] == start:
                    continue
                # ゴールに到達した場合、ansに追加する
                if col1[ii] == goal:
                    tmp = candidate[i].copy()
                    tmp.append(ii)
                    if print_log:
                        print("ans", tmp)
                    ans.append(tmp)
                    continue
                # 上記以外で探索を継続する場合、候補に追加する
                else:
                    tmp = candidate[i].copy()
                    tmp.append(ii)
                    candidate.append(tmp)
            candidate.pop(i)
    return ans

ans = way_search(start, goal, print_log=False)

経路がすべて接続されているかを確認する。NGと表示されなければOK。

for arr in ans:
    for i, loc in enumerate(arr[:-1]):
        if arr[i+1] in d[col1[loc]]["idx"]:
            continue
        else:
            print("NG", arr, loc)

経路候補の中に同じものがないかチェック。hit!が同じインデックスであればOK。

for i1, arr1 in enumerate(ans):
    for i2, arr2 in enumerate(ans):
        if arr1 == arr2:
            print("hit!", i1, i2)

最後に各経路の総距離を計算して、最短経路を表示する。

def calc_dist(ans):
    ans_dist = []
    
    for arr in ans:
        s = arr[0]
        tmp = 0
        for i in arr[1:]:
            tmp += distance(col3[s], col2[s], col3[i], col2[i])
            s = i
        ans_dist.append(tmp)
    return ans_dist

ans_dist = calc_dist(ans)
min_idx = ans_dist.index(min(ans_dist))

print(ans[min_idx])
for i in ans[min_idx]:
    print(col1[i])

結果は以下の通り。北陸経由が早いらしい。

[0, 1, 4, 5, 14, 19, 15, 16, 20]
Sapporo
Aomori
Akita
Yamagata
Niigata
Nagano
Toyama
Kanazawa
Gifu

次に本題の札幌から那覇までの最短経路を求める。今回は打ち切り件数を10万件にしておく。

start = 'Sapporo'
goal = 'Naha'

ans = way_search(start, goal, print_log=False, limit=100000)

ans_dist = calc_dist(ans)
min_idx = ans_dist.index(min(ans_dist))

print(ans[min_idx])
for i in ans[min_idx]:
    print(col1[i])

結果は以下の通り。17県通ればいけるらしい。

[0, 1, 4, 5, 14, 19, 15, 17, 24, 26, 29, 35, 38, 37, 43, 44, 46]
Sapporo
Aomori
Akita
Yamagata
Niigata
Nagano
Toyama
Fukui
Otsu
Osaka
Wakayama
Tokushima
Kochi
Matsuyama
Oita
Miyazaki
Naha

おまけとして、最長距離の経路も求められる。

max_idx = ans_dist.index(max(ans_dist))

print(ans[max_idx])
for i in ans[max_idx]:
    print(col1[i])

結果は以下の通り。40県とほぼすべての県を回っている。

[0, 4, 1, 2, 3, 6, 8, 9, 10, 7, 11, 12, 21, 18, 19, 15, 17, 20, 23, 22, 24, 28, 26, 27, 29, 35, 38, 36, 30, 32, 31, 33, 34, 43, 39, 42, 44, 45, 41, 46]
Sapporo
Akita
Aomori
Morioka
Sendai
Fukushima
Utsunomiya
Maebashi
Saitama
Mito
Chiba
Tokyo
Shizuoka
Kofu
Nagano
Toyama
Fukui
Gifu
Tsu
Nagoya
Otsu
Nara
Osaka
Kobe
Wakayama
Tokushima
Kochi
Takamatsu
Tottori
Okayama
Matsue
Hiroshima
Yamaguchi
Oita
Fukuoka
Kumamoto
Miyazaki
Kagoshima
Nagasaki
Naha