graphillionとnetworkxを用いて路線図のN番目に短い経路を求める

はじめに

「鉄道路線データを可視化し、最短経路問題を解く」に、数え上げお姉さん問題を高速に解けることで有名なライブラリ「graphillion」を適用することで、以下のような路線図のN番目に短い経路を求めた。
ここでは東京メトロの西日暮里から西新宿までの経路のうち、3番目に短いものを赤線で表す。

ライブラリのインストール

google colabを開き、以下のコードをコピペしていく。
まずはmatplotlibの日本語化に必要なライブラリとgraphillionをインストールする。

!apt-get -y install fonts-ipafont-gothic
!pip install japanize_matplotlib

!pip install graphillion

ライブラリの読み込み

import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib
import numpy as np
import networkx as nx
from graphillion import GraphSet, tutorial
import itertools

　

駅データのダウンロードと読み込み

駅データ.jpに会員登録し、ダウンロード→データダウンロードから無料データをダウンロードする。
そして、ダウンロードしたファイルをpandasで読み込む。

#csvファイルからpandsa形式のテーブルを作成
station = pd.read_csv("/content/drive/MyDrive/Colab Notebooks/railway_data/station20210312free.csv")
join = pd.read_csv("/content/drive/MyDrive/Colab Notebooks/railway_data/join20210312.csv")
line = pd.read_csv("/content/drive/MyDrive/Colab Notebooks/railway_data/line20210312free.csv")
company = ("/content/drive/MyDrive/Colab Notebooks/railway_data/company20200619.csv")

東京メトロの路線図を表示

#全国の駅から東京メトロの駅のみ抽出する 東京メトロ...company_cd == 18
metro = station[["station_cd", "station_name", "line_cd", "lon", "lat"]]
metro = pd.merge(metro, line, on = 'line_cd')
metro = metro[metro["company_cd"] == 18]
metro = metro[["station_cd", "station_name", "line_cd", "lon_x", "lat_x", "line_name", "line_color_c", "line_color_t"]]
lon = metro["lon_x"]
lat = metro["lat_x"]
metro["lon"] = lon
metro["lat"] = lat
metro = metro[["station_cd", "station_name", "line_cd", "lon", "lat", "line_name"]]

#東京メトロの接続辺を抽出する 路線...line_cd == 28001---28010
metro_join = join[(join["line_cd"]==28001)|(join["line_cd"]==28002)|(join["line_cd"]==28003)|(join["line_cd"]==28004)|(join["line_cd"]==28005)|(join["line_cd"]==28006)|(join["line_cd"]==28007)|(join["line_cd"]==28008)|(join["line_cd"]==28009)|(join["line_cd"]==28010)]
metro_join = metro_join[["station_cd1", "station_cd2"]]


#グラフの宣言
G = nx.Graph()
#頂点を駅名にする
G.add_nodes_from(metro["station_name"])


#plotの座標を設定
pos={}
for i, j, k in zip(metro["station_name"], metro["lon"], metro["lat"]):
  pos[i] = (j, k)
#リストeにstation_nameとstation_cdを格納し、リンクさせる
e = []
for i, j in zip(metro["station_name"], metro["station_cd"]):
  e.append([i, j])
#グラフに辺情報を加える
for i, j in zip(metro_join["station_cd1"], metro_join["station_cd2"]):
    for k in e:
      if k[1] == i:
        for l in e:
          if l[1] == j:
            G.add_edge(k[0], l[0])

#グラフの出力の設定
plt.figure(figsize=(10,10),dpi=200)
plt.title('東京メトロ', fontsize=20)
plt.axes().set_aspect('equal', 'datalim')
nx.draw_networkx(G, pos, node_color='b', alpha=0.8, node_size=10, font_size=5, font_family='IPAexGothic')
plt.show()

駅同士の距離をグラフに与える

#辺に重みとして駅間の距離を持たせるためのデータ作成
dist = []
cd1_lat = []
cd1_lon = []
cd2_lat = []
cd2_lon = []
not_exist = [] #joinテーブルには存在しているが、stationテーブルには存在していないstation_cdを格納
for i, j in zip(metro_join["station_cd1"], metro_join["station_cd2"]):
  flag = True
  for k, l, m in zip(metro["station_cd"], metro["lat"], metro["lon"]):
    if i == k:
      cd1_x = l
      cd1_y = m
      cd1_lat.append(l)
      cd1_lon.append(m)
      flag = False
    if j == k:
      cd2_x = l
      cd2_y = m
      cd2_lat.append(l)
      cd2_lon.append(m)
  if flag:
    not_exist.append([i, j])
    #print(i, j)
  else:
    dist.append(((cd1_x-cd2_x)**2 + (cd1_y-cd2_y)**2)**0.5)
#このまま実行するとエラー ValueError: Length of values does not match length of index
#どうやら"station_cd" == 2800701と2800702はstationテーブルに存在しておらず、joinテーブルから不要なので削除
#以下2行は...metro_join = metro_join[metro_join["station_cd1"] != 2800701]...と等価
for i in range(len(not_exist)):
  metro_join = metro_join[metro_join["station_cd1"] != not_exist[i][0]]
#joinテーブルに列を追加
metro_join["cd1_lat"] = cd1_lat
metro_join["cd1_lon"] = cd1_lon
metro_join["cd2_lat"] = cd2_lat
metro_join["cd2_lon"] = cd2_lon
metro_join["distance"] = dist

#nodes is station_name
#グラフに辺の重みを与える
for i, j, m in zip(metro_join["station_cd1"], metro_join["station_cd2"], metro_join["distance"]):
    for k in e:
      if k[1] == i:
        for l in e:
          if l[1] == j:
            G.add_edge(k[0], l[0], weight=m)


edge_weights = {(e1, e2): data['weight'] for e1, e2, data in G.edges(data=True)}

graphillionでN番目に短い経路を求める

最後にgraphillionのGraphSetオブジェクトにグラフ情報を与えてN番目に短い経路を求める。

def overlap_subgraph(G, subgraph):
  # ノードとエッジに色をつける
  node_color = ['red' if node in list(itertools.chain.from_iterable(subgraph)) else 'blue' for node in G.nodes()]
  edge_color = ['red' if edge in subgraph or (edge[1], edge[0]) in subgraph else 'black' for edge in G.edges()]

  #グラフの出力の設定
  plt.figure(figsize=(10,10),dpi=200)
  plt.title('東京メトロ', fontsize=20)
  plt.axes().set_aspect('equal', 'datalim')
  nx.draw_networkx(G, pos, node_color=node_color, edge_color=edge_color, alpha=0.8, node_size=10, font_size=5, font_family='IPAexGothic')
  plt.show()


# GraphSetオブジェクトに路線図を読み込む
GraphSet.set_universe(G.edges())

# すべてのパスを数え上げる、詳しくは「数え上げお姉さん問題」で検索
paths = GraphSet.paths("西日暮里", "西新宿")


iterator = paths.min_iter(edge_weights)
# 経路が短いパスから順に5つ取り出す
for i in range(5):
  # 最短経路から順にイテレータ、ここでedge_weightsを読み込ませることで各駅同士の距離の和が小さい順に経路を返す
  subgraph = next(iterator)
  overlap_subgraph(G, subgraph)

西日暮里から西新宿までのN番目に短い経路は以下のようになった。

おわりに

一般的にグラフの最短経路は簡単に求まるが、N番目に短い経路を求めることは意外に難しい。
graphillionはZDDというグラフ圧縮アルゴリズムを裏で動かすことで、高速に数え上げ問題を解くことができ、今回のような問題に応用できる。
他にもSNSや電力ネットワークなど、現実のグラフデータで試したいけどなかなか見つからないので、詳しい方がいたら教えて下さい。