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 法はトレードオフという概念を初めて自分に突きつけてくれた思い入れのあるアルゴリズムの一つです。