はじめに

友だちPython シリーズのご紹介

友だちPython シリーズは、Python の小ネタを短文でお届けします。
小さなエピソードを重ねてPythonと仲良しになれたら、と願ってシリーズ名を付けました。

話題

Pythonによる実務で役立つ最適化問題100+ (1)

書籍「Pythonによる実務で役立つ最適化問題100+ (1)」（朝倉書店）をご紹介いたします。

この書籍シリーズは、いわゆる「最適化問題」について、コードを写経しつつ、基礎から学べる良書です。
数式とPythonコード中心の硬派な一冊です。
対象の書籍は３巻構成の第１巻（初版第１刷）です。

この記事は networkx と 私の戯れの記録です。

記事中のイラストは「いらすとや」さんからお借りしました。
ありがとうございます！

今回の友だち

NetworkX 公式サイトの Gallery ページをご覧ください。

話題の概要

第２章「最短路問題」の「2.2.2 １対全の最適化問題」で networkx の数理的な使い方を学びました。
２つの遊びポイントを描きます。

実行環境
・networkx ver. 2.8.8

遊んだ内容１

networkx を用いて、Dijkstra法による「１つの始点から他の全点に対する最短路」を算出し、迷路のような図を描画します。

■ 遊び方
書籍の28～31ページのコードを混ぜ混ぜして次のサンプルコードを作りました。
コメントだらけになってしまいました＼(^o^)／

### 5×7の格子グラフの最短路木を描画

import networkx as nx
from networkx import DiGraph
import matplotlib.pyplot as plt
import random

### 乱数を用いて2次元格子グラフを作成
# 初期値設定
x, y = 5, 7             # 格子の数を設定(x, y)
lower, upper = 1, 30    # 重み：コストの範囲を設定
random.seed(123)        # 乱数シード
# networkxの2次元格子グラフGridの初期化
Grid = nx.grid_2d_graph(x, y)
# Gridの辺に重み：コストを設定（ランダムな整数値）
for (i, j) in Grid.edges():
    Grid[i][j]['weight'] = random.randint(lower, upper)

### Gridの最短絡とコストを計算: 始点の座標は(0, 0)
pred, distance = nx.dijkstra_predecessor_and_distance(Grid, source=(0, 0))
print('predの一部 '
      f'{[{k: v} for k, v in pred.items() if k in ((0, 0), (1, 0), (2, 0))]}')

### 最短路木の描画
# networkxの有向グラフDirectedの初期化
Directed = nx.DiGraph()
# predに保持する最短路をDirectedの辺に追加
for i in pred:
    if len(pred[i]) >= 1:
        Directed.add_edge(pred[i][0], i)

# Gridの頂点の座標情報の取得
pos = {(i, j): (i, j) for (i, j) in Grid.nodes()}

# Gridの辺の重み：コスト情報の取得
edge_labels = {}
for (i, j) in Grid.edges():
    edge_labels[i, j] = f'{ Grid[i][j]["weight"] }'

## Directedを用いて最短路木を描画
plt.figure(figsize=(5, 7))
# 有向グラフの描画
nx.draw(
    Directed, pos=pos, width=2,
    node_color='tomato', edge_color='tab:blue', node_size=50,
    with_labels=True, font_size=10,        # 座標ラベルを表示する
)
# 辺の重み：コストの描画
nx.draw_networkx_edge_labels(
    Directed, pos=pos,
    edge_labels=edge_labels, font_size=10, # 重みラベルを表示する
)
plt.show()

【実行結果】

【図の見方】
左下の 頂点 (0,0) からスタートして、青い線（辺です）の上の数字（コストです）が小さくなるように（最小化です）、全頂点を青い矢印線で結んでいます。

【経路は pred に格納】
最小化の計算を networkx の dijkstra_predecessor_and_distance() が行って、経路情報を pred に格納しています。
実行結果の最初の pred の一部に注目しましょう。

始点 (0, 0) の前ステップには何もなし（[]）です。
頂点 (1, 0) の前ステップは (0, 0)、つまり始点です。
右下の (0, 0) から (1, 0) へ青い線（辺です）が矢印付きで（有向です）繋がっています。
同じく、頂点 (2, 0) の前ステップは (1, 0) です。
(1, 0) から (2, 0) へ青い線（辺です）が矢印付きで（有向です）繋がっています。
「predに保持する最短路をDirectedの辺に追加」して、図を描いています！

【感想】
networkx はグラフィックツールだと思って今まで使ってきました。
この書籍を学んだことによって、最適化問題を解けるのだと知りました。
networkx は凄いです！！！

遊んだ内容２

「2.4 ベンチマーク問題例の読み込みとDial法」で紹介されている道路ネットワークのデータ取得方法について、お知らせしたいことがあります。
テキスト掲載のURLではデータ取得サイトにアクセスできないようです。

次のリンクでダウンロードサイトに飛んでくださいね🚀

https://www.diag.uniroma1.it/challenge9/download.shtml

ダウンロードするファイルは下図の「２箇所の赤枠のファイル」です！
「NY のデータ」です。

ダウンロードした gzip ファイルを解凍して、読み込みフォルダ（作業フォルダ等）に適切に配置することで、書籍39ページの コードが正しく作動するようになります。

次の２つのコードは、Python でローカルPCにダウンロードしてファイルを解凍するコードです。参考まで。

### ファイルのダウンロード

# 設定：ダウンロードファイルの格納先パス＋ファイル名
filename1 = r'./USA-road-d.NY.gr'
filename2 = r'./USA-road-d.NY.co'
ext = '.gz'

# ダウンロードの実行 filename1
url1 = r'https://www.diag.uniroma1.it/challenge9/data/USA-road-d/USA-road-d.NY.gr.gz'
res1 = requests.get(url1)
with open(filename1 + ext, mode='wb') as f:
  f.write(res1.content)

# ダウンロードの実行 filename2
url2 = r'https://www.diag.uniroma1.it/challenge9/data/USA-road-d/USA-road-d.NY.co.gz'
res2 = requests.get(url2)
with open(filename2 + ext, mode='wb') as f:
  f.write(res2.content)

### gzファイルの解凍
# 参考サイト: https://wannabe-data-engineer.net/python-gunzip/

# ライブラリのインポート
import gzip
import shutil

# ファイルの解凍 filename1
with gzip.open(filename1 + ext, mode='rb') as gz_file:
    with open(filename1, mode='wb') as file:
        shutil.copyfileobj(gz_file, file)

# ファイルの解凍 filename2
with gzip.open(filename2 + ext, mode='rb') as gz_file:
    with open(filename2, mode='wb') as file:
        shutil.copyfileobj(gz_file, file)

続いて、もう一つお知らせがあります。
書籍39ページの「データを読み込むコード」でワーニングが出ます。

ワーニングの原因はNumPyの「np.int0」が廃止予定入りしたためです。

ワーニング回避コードを置いておきます！

### 実行

filename = 'USA-road-d.NY'
G, n, m, C = ReadDimacs(filename)
print(f'num of nodes in the largest connected componet={n}, max length+1={C}')

source = 0      # 始点
sink = 50000    # 終点

settled = np.zeros(n, np.intp) # ワーニング回避
reached = np.zeros(n, np.intp) # ワーニング回避
dist = np.zeros(n, np.intp)    # ワーニング回避
prev = np.zeros(n, np.intp)    # ワーニング回避

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

これで無事に書籍39～40ページのコードを動かせます！
40ページのコードで描画した networkx の「NYの最短路のパス図」をご覧ください！

みなさんも最適化問題を存分に堪能してください！

おわり

ブログの紹介

noteで５つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計
統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。

２．実験！たのしいベイズモデリングをPyMC Ver.5で
書籍「たのしいベイズモデリング」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で
書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析トピックを PythonとPyMC Ver.5で取り組みます。
豊富なテーマ（お題）を実践することによって、PythonとPyMCの基礎体力づくりにつながる、そう信じています。
日々、Web検索に勤しみ、時系列モデルの理解、Pythonパッケージの把握、R・Stanコードの翻訳に励んでいます！
このシリーズがPython時系列分析の入門者の参考になれば幸いです🍀

４．Python機械学習プログラミング実践記
書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learn と PyTorch の教科書です。
よかったらぜひ、お試しくださいませ。

５．データサイエンスっぽいことを綴る
統計、データ分析、AI、機械学習、Python のコラムを不定期に綴っています。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

最後までお読みいただきまして、ありがとうございました。