GRL_1_C - All Pairs Shortest Path / PythonでAizu Online Judgeを解く

Posted on 2018/10/11

概要

重み付きの有向グラフ \(G(V, E)\) が与えられます。任意の2頂点について最短経路の距離を求めてください。また、グラフが負閉路を持つ場合はそれを判定しましょう。

制約

  • \(1 \leq |V| \leq 100\)
  • \(0 \leq |E| \leq 9900\)
  • \(-2 \times 10^7 \leq d_i \leq 2 \times 10^7\)
  • グラフ \(G\) は平行辺を持たない
  • グラフ \(G\) は自己ループを持たない

入力

%python
def input(f):
    global v, e, s, t, d
    
    v, e = map(int, f.readline().split())
    s = []
    t = []
    d = []
    for _ in range(e):
        s_, t_, d_ = map(int, f.readline().split())
        s.append(s_)
        t.append(t_)
        d.append(d_)

サンプル入力

%sh
cat << EOF > /tmp/input1
4 6
0 1 1
0 2 5
1 2 2
1 3 4
2 3 1
3 2 7
EOF

cat << EOF > /tmp/input2
4 6
0 1 1
0 2 -5
1 2 2
1 3 4
2 3 1
3 2 7
EOF

cat << EOF > /tmp/input3
4 6
0 1 1
0 2 5
1 2 2
1 3 4
2 3 1
3 2 -7
EOF

Python: 入力をnetworkxで可視化してみる

%sh
apt update && apt install -y python3-pip
pip3 install networkx matplotlib
WARNING: apt does not have a stable CLI interface. Use with caution in scripts.

Hit:1 http://cran.rstudio.com/bin/linux/ubuntu xenial/ InRelease
Hit:2 http://security.ubuntu.com/ubuntu xenial-security InRelease
Hit:3 http://archive.ubuntu.com/ubuntu xenial InRelease
Hit:4 http://archive.ubuntu.com/ubuntu xenial-updates InRelease
Hit:5 http://archive.ubuntu.com/ubuntu xenial-backports InRelease
Reading package lists...
Building dependency tree...
Reading state information...
69 packages can be upgraded. Run 'apt list --upgradable' to see them.

WARNING: apt does not have a stable CLI interface. Use with caution in scripts.

Reading package lists...
Building dependency tree...
Reading state information...
python3-pip is already the newest version (8.1.1-2ubuntu0.4).
0 upgraded, 0 newly installed, 0 to remove and 69 not upgraded.
Requirement already satisfied (use --upgrade to upgrade): networkx in /usr/local/lib/python3.5/dist-packages
Collecting matplotlib
  Downloading https://files.pythonhosted.org/packages/7b/ca/8b55a66b7ce426329ab16419a7eee4eb35b5a3fbe0d002434b339a4a7b09/matplotlib-3.0.0-cp35-cp35m-manylinux1_x86_64.whl (12.8MB)
Requirement already satisfied (use --upgrade to upgrade): decorator>=4.3.0 in /usr/local/lib/python3.5/dist-packages (from networkx)
Collecting kiwisolver>=1.0.1 (from matplotlib)
  Downloading https://files.pythonhosted.org/packages/7e/31/d6fedd4fb2c94755cd101191e581af30e1650ccce7a35bddb7930fed6574/kiwisolver-1.0.1-cp35-cp35m-manylinux1_x86_64.whl (949kB)
Collecting numpy>=1.10.0 (from matplotlib)
  Downloading https://files.pythonhosted.org/packages/75/22/355e68c80802d6f488223788fbda75c1daab83c3ef609153676c1f17be5f/numpy-1.15.2-cp35-cp35m-manylinux1_x86_64.whl (13.8MB)
Collecting cycler>=0.10 (from matplotlib)
  Using cached https://files.pythonhosted.org/packages/f7/d2/e07d3ebb2bd7af696440ce7e754c59dd546ffe1bbe732c8ab68b9c834e61/cycler-0.10.0-py2.py3-none-any.whl
Collecting pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 (from matplotlib)
  Downloading https://files.pythonhosted.org/packages/2b/4a/f06b45ab9690d4c37641ec776f7ad691974f4cf6943a73267475b05cbfca/pyparsing-2.2.2-py2.py3-none-any.whl (57kB)
Collecting python-dateutil>=2.1 (from matplotlib)
  Using cached https://files.pythonhosted.org/packages/cf/f5/af2b09c957ace60dcfac112b669c45c8c97e32f94aa8b56da4c6d1682825/python_dateutil-2.7.3-py2.py3-none-any.whl
Requirement already satisfied (use --upgrade to upgrade): setuptools in /usr/lib/python3/dist-packages (from kiwisolver>=1.0.1->matplotlib)
Collecting six (from cycler>=0.10->matplotlib)
  Downloading https://files.pythonhosted.org/packages/67/4b/141a581104b1f6397bfa78ac9d43d8ad29a7ca43ea90a2d863fe3056e86a/six-1.11.0-py2.py3-none-any.whl
Installing collected packages: kiwisolver, numpy, six, cycler, pyparsing, python-dateutil, matplotlib
Successfully installed cycler-0.10.0 kiwisolver-1.0.1 matplotlib-3.0.0 numpy-1.15.2 pyparsing-2.2.2 python-dateutil-2.7.3 six-1.11.0
You are using pip version 8.1.1, however version 18.1 is available.
You should consider upgrading via the 'pip install --upgrade pip' command.

%python
import os

if os.environ.get('ZEPPELIN_HOME'):
    import networkx as nx
    import matplotlib.pyplot as plt
    
    def visualize():
        G = nx.MultiDiGraph()
        for i in range(v):
            G.add_node(i, node_color='b')
        for i in range(e):
            G.add_edge(s[i], t[i], t=d[i])

        plt.figure()
        plt.style.use('seaborn')
        plt.axis("off")
        pos = nx.spring_layout(G)
        nodes = nx.draw_networkx_nodes(G, pos, node_color='w', linewidths=2)
        nodes.set_edgecolor('black')
        nx.draw_networkx_labels(G, pos)
        nx.draw_networkx_edges(G, pos)

        edge_labels = dict([((a, b), c['t']) for a, b, c in G.edges(data=True)])
        nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, label_pos=0.25, font_size=16)

        plt.show()
        z.show(plt)

%python
if os.environ.get('ZEPPELIN_HOME'):
    with open('/tmp/input1') as f:
        print('Example #1')
        input(f)
        visualize()
Example #1

%python
if os.environ.get('ZEPPELIN_HOME'):
    with open('/tmp/input2') as f:
        print('Example #2')
        input(f)
        visualize()
Example #2

%python
if os.environ.get('ZEPPELIN_HOME'):
    with open('/tmp/input3') as f:
        print('Example #3')
        input(f)
        visualize()
Example #3

頂点2と3の間を往復することでコストを下げ続けることができます。

考え方 - Warshall–Floyd Algorithm

100頂点くらいの規模感なグラフ上ですべての頂点のペアについて最短経路問題を解くとなると Warshall–Floyd 法を使え使えと言われているような気がします。このアルゴリズムでは2つの頂点について別のある頂点を経由したときに掛かるコストと現在の最小コストを比較して、別の頂点を経由したほうがコストが低ければ最小コストを更新していきます。

何となくで運用できてしまうアルゴリズムだと思うのですが、仕組みとしては動的計画法に基づいたものです。

参考: “ALGORITHM NOTE Warshall Floyd’s algorithm ワーシャル フロイド”
http://algorithms.blog55.fc2.com/blog-entry-52.html

$$
A^k[i, j] = 頂点 \{1, 2, 3, \dots, k\} のみを経由したときの頂点iと頂点jの最小コスト
$$

$$
P^k[i, j] = そのような最短経路の一つ
$$

とすると \(k=0\) のとき \(A^0[i, j]\) が最小コストであるのは明らかです。また \(k>0\) のとき、\(A^k\) について \(P^k[i, j]\) が頂点 k を経由しない場合は \(P^{k-1}[i, j]\) と等しくなり、頂点kを経由する場合は \(P^k[i, k]\) と \(P^k[k, j]\) の2つの部分路に分けることができてそれぞれ \(P^{k-1}[i, k]\) と \(P^{k-1}[k, j]\) に等しいものと考えることができます。したがって次のような漸化式を立てることができます。

$$
A^k[i,j] = min(A^{k-1}[i, j], A^{k-1}[i,k] + A^{k-1}[k,j])
$$

また、\(A^k[i, i] = 0\) より

$$
A^k[i,k] = min(A^{k-1}[i, k], A^{k-1}[i,k] + A^{k-1}[k,k]) = A^{k-1}[i,k]
$$

$$
A^k[k,j] = min(A^{k-1}[k, j], A^{k-1}[k,k] + A^{k-1}[k,j]) = A^{k-1}[k,j]
$$

となるので、実装においては各ステップで2次元配列の値を使い回すことができます。

実装

%python
import sys


def init():
    global inf
    inf = sys.maxsize

def generate_graph():
    global g
    g = [[inf for _ in range(v)] for _ in range(v)]
    for i in range(v):
        g[i][i] = 0
    for i in range(e):
        g[s[i]][t[i]] = d[i]
        

def solve():
    init()
    generate_graph()
    
    for k in range(v):
        for i in range(v):
            for j in range(v):
                if g[i][k] == inf or g[k][j] == inf:
                    continue
                g[i][j] = min([g[i][j], g[i][k] + g[k][j]])

    for i in range(v):
        if g[i][i] != 0:
            print('NEGATIVE CYCLE')
            return
        
    for i in range(v):
        for j in range(v):
            if j > 0:
                print(' ', end='')
            if g[i][j] == inf:
                print('INF', end='')
            else:
                print(g[i][j], end='')
        print('')

Example Output #1

%python
if os.environ.get('ZEPPELIN_HOME'):
    os.environ['INPUT'] = '/tmp/input1'

with open(os.environ.get('INPUT') or '/dev/stdin') as f:
    input(f)
    solve()
0 1 3 4
INF 0 2 3
INF INF 0 1
INF INF 7 0

Example Output #2

%python
if os.environ.get('ZEPPELIN_HOME'):
    os.environ['INPUT'] = '/tmp/input2'

    with open(os.environ.get('INPUT') or '/dev/stdin') as f:
        input(f)
        solve()
0 1 -5 -4
INF 0 2 3
INF INF 0 1
INF INF 7 0

Example Output #3

%python
if os.environ.get('ZEPPELIN_HOME'):
    os.environ['INPUT'] = '/tmp/input3'

    with open(os.environ.get('INPUT') or '/dev/stdin') as f:
        input(f)
        solve()
NEGATIVE CYCLE

計算量的には重いけれど実装はとてもシンプルで素早く問題を解決することができます。Warshall–Floyd 法はトレードオフという概念を初めて自分に突きつけてくれた思い入れのあるアルゴリズムの一つです。