Python3 – アルゴリズム(重み付きグラフ・最短経路探索・ダイクストラのアルゴリズム)

参考:欲張り法 (greedy strategy)

重み付きグラフ

networkxでつくった重みグラフ

networkxのコードサンプル

上記は重み付きのグラフです。これはnetworkxで書きました。コードは下記です。

import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
g = [
    [ 0, 10,  0, 20,  0,  0,  0, 30],   # A (0)
    [10,  0, 10,  0,  0,  0,  0,  0],   # B (1)
    [ 0, 10,  0,  0, 20,  0,  0,  0],   # C (2)
    [20,  0,  0,  0,  0, 20,  0,  0],   # D (3)
    [ 0,  0, 20,  0,  0,  0,  0, 20],   # E (4)
    [ 0,  0,  0, 20,  0,  0, 10,  0],   # F (5)
    [ 0,  0,  0,  0,  0, 10,  0, 10],   # G (6)
    [30,  0,  0,  0, 20,  0, 10,  0]    # H (7)
]
nodes = np.array(['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h'])
edges = []
edge_labels = {}
for hi, hv  in enumerate(g):
    for wi, wv in enumerate(hv):
        if(wv):
            tpl = (nodes[hi], nodes[wi])
            edges.append(tpl)
            edge_labels[tpl] = wv
G = nx.Graph()
G.add_nodes_from(nodes)
G.add_edges_from(edges)
pos = nx.spring_layout(G)
nx.draw_networkx(G, pos, with_labels=True)
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels)
plt.axis("off")
plt.show()

この重み付きグラフは参考サイトにのってるやつです。参考サイトのダイクストラのアルゴリズムをやってみます。

ダイクストラのアルゴリズム

重み付きグラフの最短経路探索のアルゴリズムの一つです。効率は結構いいですが、各地点のコストに負の値があると使えません。

コードサンプル

def print_path(prev, cost):
    for i in range(len(prev)):
        print("%d, prev = %d, cost = %d" %  (i, prev[i], cost[i]))

def get_path(start, goal, prev):
    path = []
    now = goal
    path.append(now)
    while True:
        path.append(prev[now])
        if prev[now] == start: break
        now = prev[now]
    path.reverse()
    return path

def search(glaph, start, goal):
    MAX_VAL = 0x10000000
    g_size = len(glaph)
    visited = [False] * g_size
    cost = [MAX_VAL] * g_size
    prev = [None] * g_size
    cost[start] = 0
    prev[start] = start
    now = start
    while True:
        min = MAX_VAL
        next = -1
        visited[now] = True
        for i in range(g_size):
            if visited[i]: continue
            if glaph[now][i]:
                tmp_cost = glaph[now][i] + cost[now]
                if cost[i] > tmp_cost:
                    cost[i] = tmp_cost
                    prev[i] = now
            if min > cost[i]:
                min = cost[i]
                next = i
        if next == -1: break
        now = next
    print_path(prev, cost)
    return [get_path(start, goal, prev), cost[goal]]

g = [
    [ 0, 10,  0, 20,  0,  0,  0, 30],   # A (0)
    [10,  0, 10,  0,  0,  0,  0,  0],   # B (1)
    [ 0, 10,  0,  0, 20,  0,  0,  0],   # C (2)
    [20,  0,  0,  0,  0, 20,  0,  0],   # D (3)
    [ 0,  0, 20,  0,  0,  0,  0, 20],   # E (4)
    [ 0,  0,  0, 20,  0,  0, 10,  0],   # F (5)
    [ 0,  0,  0,  0,  0, 10,  0, 10],   # G (6)
    [30,  0,  0,  0, 20,  0, 10,  0]    # H (7)
]

path, cost = search(g, 5, 1)
print(path, 'cost: ', cost)

結果

0, prev = 3, cost = 40
1, prev = 0, cost = 50
2, prev = 4, cost = 60
3, prev = 5, cost = 20
4, prev = 7, cost = 40
5, prev = 5, cost = 0
6, prev = 5, cost = 10
7, prev = 6, cost = 20
[5, 3, 0, 1] cost:  50

上記はfからbまでの最短経路を探しております。どういう動きなのかを言葉で書くのが難しい。。

アルゴリズムの仕組み

最初にする主なこと

・スタート地点のコストを0にする。
・現在地をスタート地点にする。

ループでしている主なこと

・現在地を訪問済み(チェック済み)にする。
・現在地から訪問可能な地点の、各コスト(現在地までのコスト+現在地から訪問可能地点までのコスト)を出して、今まで記憶しているコストより低ければ更新する。
・訪問済みでない地点の記憶されているコストで最も低いコストの地点を現在地にする。
・以降、上記3つを繰り返す。
・訪問済みの地点がなくなったらループを終了させる。

ループ終了後どうなっているか?

・スタート地点からの各地点までの最短コストと、パスを記憶している。上記コードの場合、cost配列にコストが入っていて、prev配列に各地点の直前の地点が入っている。

注意点

・各地点のコストがマイナスのものがあるとこのアルゴリズムは使えない。

改善点

・ヒープを使うとより効率的になるらしいので、次にヒープを使ってつくってみる。
・最短経路をnetworkxで表示させるようにしたい。

Python3 – ヒープをnetworkxではなくgraphvizで書いてみる

networkxだときれいに表示できない。posに表示スタイルを設定するようですが、spring_layout、spectral_layout、shell_layoutとかありますが、全部キレイな木を表示してくれる感じじゃない。外部のdotファイルとかいうのを読み込むとできるようですが、graphvizというのがあるので、そっちをやってみます。

Networkxの場合

import networkx as nx
import matplotlib.pyplot as plt

h = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
G = nx.Graph()
G.add_nodes_from(h)
for i, _ in enumerate(h):
    if(i - 1) // 2 >= 0:
        G.add_edge(h[(i - 1) // 2], h[i])
pos = nx.spring_layout(G)
nx.draw_networkx(G, pos, with_labels=True)
plt.axis("off")
plt.show()

Graphvizの場合

インストール

ここからGraphvizをダウンロードしてインストールします。
その後、下記のようにpipインストールします。

$ pip install graphviz

あと、windows10ですが、pathを通す必要がありました。インストールしたフォルダのbinをフォルダのパスを通してください。

サンプルコード

from graphviz import Graph

h = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
G = Graph(format='png')
G.attr('node', shape='circle')
for i in h: G.node(str(i))
for i, _ in enumerate(h):
    if(i - 1) // 2 >= 0:
        G.edge(str(h[(i - 1) // 2]), str(h[i]))
print(G)
G.render('heap')

参考:graphvizを使ってPython3で木構造を描く

Python3 – ヒープ

ヒープはノードに常に子供が2個あって、子供より親は小さい値を持ちます。ヒープは配列で表せます。配列の0番目は一番上の根ノードです。あとは、k番目のノードの子供は、左が2k+1、右が2k+2番目になります。挿入は最後に挿入して、親ノードと比較して必要に応じて交換することを繰り返します。

Python3にはヒープのデータ構造がデフォルトで用意されているようなので、使ってみます。heapqです。これはヒープキュー、優先度キュー、プライオリティキューとかいわれるものみたいで、「優先度が高い順に取り出せるデータ構造」だそうです。(参考:優先度キュー

普通の配列をヒープにする

コード

h = [10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0]
hq.heapify(h)
print(h)

結果

[0, 1, 4, 2, 6, 5, 8, 3, 7, 9, 10]

コード

h = [10, 5, 12, 43, 1, 30, 45, 12, 14, 15, 43, 53]
hq.heapify(h)
print(h)

結果

[1, 5, 12, 12, 10, 30, 45, 43, 14, 15, 43, 53]

ヒープに挿入する

import heapq as hq

h = [
    (5, 'a'),
    (6, 'b'),
    (3, 'c'),
    (1, 'd'),
    (9, 'e')
]
hq.heapify(h)
hq.heappush(h, (7, 'f'))
hq.heappush(h, (2, 'g'))
hq.heappush(h, (2, 'h'))
print(h)

結果

[(1, 'd'), (2, 'h'), (2, 'g'), (5, 'a'), (9, 'e'), (7, 'f'), (3, 'c'), (6, 'b')]

ヒープからpopする

import heapq as hq

h = [
    (5, 'a'),
    (6, 'b'),
    (3, 'c'),
    (1, 'd'),
    (9, 'e')
]
hq.heapify(h)
hq.heappush(h, (7, 'f'))
hq.heappush(h, (2, 'g'))
hq.heappush(h, (4, 'h'))
print(h)
print('--- heappop ---')
print(hq.heappop(h))
print(h)
print('--- heappushpop ---')
print(hq.heappushpop(h, (1, 'i')))
print(h)
print('--- heapreplace ---')
print(hq.heapreplace(h, (1, 'j')))
print(h)

結果

[(1, 'd'), (4, 'h'), (2, 'g'), (5, 'a'), (9, 'e'), (7, 'f'), (3, 'c'), (6, 'b')]
--- heappop ---
(1, 'd')
[(2, 'g'), (4, 'h'), (3, 'c'), (5, 'a'), (9, 'e'), (7, 'f'), (6, 'b')]
--- heappushpop ---
(1, 'i')
[(2, 'g'), (4, 'h'), (3, 'c'), (5, 'a'), (9, 'e'), (7, 'f'), (6, 'b')]
--- heapreplace ---
(2, 'g')
[(1, 'j'), (4, 'h'), (3, 'c'), (5, 'a'), (9, 'e'), (7, 'f'), (6, 'b')]

heappushpopはpushしてから、popする。heappushしてからheappopするよりも効率的だそうです。pushしてからpopするので、pushしたものも含めて最小値をpopします。heapreplaceはpopしてからpushするので、pushした要素が最小値でもpopされません。

Python3 – アルゴリズム(グラフ)

参考:http://www.geocities.jp/m_Hiroi/light/pyalgo05.html

グラフをプログラムする場合、よく使われる方法に「隣接行列」と「隣接リスト」があります。隣接行列は 2 次元配列で頂点の連結を表す方法です。頂点が N 個ある場合、隣接行列は N 行 N 列の行列で表すことができます。

グラフの例

グラフをプログラムで書く

隣接行列

縦横がそれぞれa-gまでのノードを表す。エッジがある場合は1にする。g[0][1]は、a-b間のエッジがあることを示す。

g = np.array([
    [0, 1, 1, 0, 0, 0, 0],
    [1, 0, 1, 1, 0, 0, 0],
    [1, 1, 0, 0, 1, 0, 0],
    [0, 1, 0, 0, 1, 1, 0],
    [0, 0, 1, 1, 0, 0, 1],
    [0, 0, 0, 1, 0, 0, 0],
    [0, 0, 0, 0, 1, 0, 0]
])

隣接行列の欠点は、辺の数が少ない場合でも N 行 N 列の行列が必要になることです。つまり、ほとんどの要素が 0 になってしまい、メモリを浪費してしまうのです。この欠点を補う方法に隣接リストがあります。これはつながっている頂点を格納する方法です。

隣接リスト

つながってる頂点を配列にいれる。

g = np.array([
    [1, 2],    # A
    [0, 2, 3], # B
    [0, 1, 4], # C
    [1, 4, 5], # D
    [2, 3, 6], # E
    [3],       # F
    [4]        # G
])

探索

深さ優先探索

def search(goal, path):
    n = path[-1]
    if n == goal:
        print(path)
    else:
        for x in g[n]:
            if x not in path:
                path.append(x)
                search(goal, path)
                path.pop()

使い方

search(6, [0])

結果

[0, 1, 2, 4, 6]
[0, 1, 3, 4, 6]
[0, 2, 1, 3, 4, 6]
[0, 2, 4, 6]

再帰をうまく使うと芸術的だな。分かりづらいけど。goalは到達したい場所。pathはスタート地点と通過したい場所を配列で入れる。
AからGに行きたい場合は、search(6, [0])という風にやります。最初にpathの最後の要素を取り出しそれがゴール地点だったら、pathつまりスタート地点からゴール地点までの順番を表示しております。ゴール地点でなければ、現在位置からいける地点を順に取り出し、pathに加えます。その際に、戻らないように今まで通った場所は、次にいける地点から除外します。そして、またsearch(goal, path)を実行します。再帰的に繰り返すことで、Goalに到達可能であれば、いずれGoalが次にいける地点に現れます。これによって、goalへの道順を発見できたことになります。最後のsearchの実行で、pathがprintされます。これで最後のsearchは終了します。もしゴールに到達できなければ、最終的にはg[n]が空になるか、if x not in pathにマッチしないものだけになるか、になるので、これもまた最後のsearchは終了します。そうすると1個前のsearchの後のpopによってpathが1つ前の状態に戻り、forが残ってたら残りのforを回すし、なければ終わりになって、もう1つ前のsearchの実行に戻りまっす。まあ文章で書いてもめちゃくちゃになるけど。再帰を自分でパッと思いつけるようになるのは結構大変だろうな。

これが深さ優先探索といわれるのは、到達可能な深い場所にとりあえずどんどんいって、また戻っては次にいける深いところにどんどん行くことを繰り返すからだと思います。最短経路探索と言われているものをしようとしても、全部やってみて一番短いやつを選択することになるので、いけるところは全部見る必要があるのかなと思います。

幅優先探索

幅優先探索は、最短経路探索に適していて深さ優先探索より効率的だそうです。ただし、深さ優先の場合は諸突猛進型でとりあえず一番深いところまで行って終わったら、戻ってきて、もう今までのことは忘れるのですが、幅優先の場合は、全ての経路を並行して深堀っていきます。そのため、最短経路が発見されるまで全ての経路の道順を記憶しておく必要があり、メモリの必要量が深さ優先に比べて圧倒的に多くなる可能性があります。でも、ものすごい深いグラフだと深さ優先グラフでもメモリが膨大になる可能性があります。

def search(start, goal):
    q = [[start]]
    while len(q) > 0:
        path = q.pop(0)
        n = path[-1]
        if n == goal:
            print(path)
        else:
            for x in g[n]:
                if x not in path:
                    new_path = path[:]
                    new_path.append(x)
                    q.append(new_path)

使い方

search(0, 6)

結果

[0, 2, 4, 6]
[0, 1, 2, 4, 6]
[0, 1, 3, 4, 6]
[0, 2, 1, 3, 4, 6]

qをキューとして使っている。先入れ先出ししている。先入れ先出し。このqで全部の経路を記憶しており、最初に記憶したものは階層が低いパスになるので、最初に記憶したものから順番に取り出している。qがからになるまでwhileを続けている。popで一番古いパスを取り出して、そのパスの最後尾がゴールだったらprintで経路を書き出している。ゴールでない場合、現在地から次にいける地点を取り出して、順番に処理している。処理内容は、単純に今までの経路に取り出した次にいける地点を追加したものを、qに加えているだけ。

このキューをスタックにすると、深さ優先探索になるらしいのでやってみる。

深さ優先探索のwhileバージョン

def search(start, goal):
    q = [[start]]
    while len(q) > 0:
        path = q.pop()
        n = path[-1]
        if n == goal:
            print(path)
        else:
            for x in g[n]:
                if x not in path:
                    new_path = path[:]
                    new_path.append(x)
                    q.append(new_path)

こういうことでしょうか?キューは先入れ先出し、スタックは後入れ先出。qに入っているものをpop()で末尾から取り出せば、後入れ先出になる。とりあえず全部探索できるし、結果もあってるんだけど、最初に再帰でやったやつとは出力順序が異なる。これは、再帰でやったやつは、各配列の先頭から深堀していくのに対して、whileの場合は、配列の末尾から深堀していくから。メモリ的には再帰でやった方がいいんじゃないかと思う。再帰の場合記憶しているのは一つの経路だけだから。

反復深化

反復深化というのは、深さ優先の深堀する最大階層を決めるやつ。1階層目まででまず深さ優先をやって、ダメだったら2階層目まででやるという感じで進めていく。必要となるメモリ量は深さ優先と同じになるので、幅優先と比べるとメモリ使用量は少ない。ただし、1階層目でだめだったら、2階層目までの制限で再度一からやり直すことになり、同じ探索を何度もすることになるため、時間的には無駄が多く効率が悪い。

def id_search(limit, goal, path):
    n = len(path)
    m = path[-1]
    if n == limit:
        if m == goal: print(path)
    else:
        for x in g[m]:
            if x not in path:
                path.append(x)
                id_search(limit, goal, path)
                path.pop()

for x in range(1, g.shape[0] + 1):
    print(x, 'moves')
    id_search(x, 6, [0])

結果

1 moves
2 moves
3 moves
4 moves
[0, 2, 4, 6]
5 moves
[0, 1, 2, 4, 6]
[0, 1, 3, 4, 6]
6 moves
[0, 2, 1, 3, 4, 6]
7 moves

Python3 – NetworkXでグラフを表示する

https://networkx.github.io/

インストール

https://github.com/networkx/networkx/

$ pip install networkx
$ pip install decorator

Decoratorというのも必要っぽい。両方ともすでに入ってた。

サンプルコード

import networkx as nx

G=nx.Graph()
G.add_node("spam")
G.add_edge(1,2)
print(list(G.nodes()))
print(list(G.edges()))

結果

[1, 2, 'spam']
[(1, 2)]

使い方

matplotlibでグラフを表示する

https://networkx.readthedocs.io/en/stable/tutorial/tutorial.html

import matplotlib.pyplot as plt
nx.draw(G)
plt.show()

サンプルコード

import numpy as np
import matplotlib.pyplot as plt
import networkx as nx

#隣接行列
g = np.array([
    [0, 1, 1, 0, 0, 0, 0],
    [1, 0, 1, 1, 0, 0, 0],
    [1, 1, 0, 0, 1, 0, 0],
    [0, 1, 0, 0, 1, 1, 0],
    [0, 0, 1, 1, 0, 0, 1],
    [0, 0, 0, 1, 0, 0, 0],
    [0, 0, 0, 0, 1, 0, 0]
])
nodes = np.array(['a', 'b', 'c', 'd', 'e', 'f', 'g'])
G = nx.Graph()
G.add_nodes_from(nodes)
edges = []
for hi, hv  in enumerate(g):
    for wi, wv in enumerate(hv):
        if(wv): edges.append((nodes[hi], nodes[wi]))
G.add_edges_from(edges)
pos = nx.spring_layout(G)
nx.draw_networkx(G, pos, with_labels=True)
plt.axis("off")
plt.show()

参考:
[Python]NetworkXでQiitaのタグ関係図を描く
Pythonで迷路を解く – アルゴリズムクイックリファレンス6章の補足 –
NetworkXでグラフを描いた(最短経路他)

Python3 – アルゴリズム(累乗)・計算量

累乗

pow0は遅いやつ。計算量はループをn回回すので、O(n)という。pow1は速いやつ。計算量はO(log n)らしい。

def pow0(x, n):
    value = 1
    for i in range(n): value *= x
    return value

def pow1(x, n):
    if n == 0: return 1
    value = pow1(x, int(n / 2))
    value *= value
    if n % 2 == 1: value *= x
    return value

ちなみに、自分のPCで速度を計ったら、下記だった。(python3内蔵のpow関数は、pow1と同じ速度だった)
print(pow0(2, 1000000))は17.3秒
print(pow1(2, 1000000))は1.3秒

なんで、pow1の計算量はO(log n)なのか?

pow1を呼び出す度に、nを半分にしている(2のマイナス1乗を掛けてる)。nを2のk乗とおくと、pow1の呼び出し回数は、k回となる。よってO(k)となる、n = 2のk乗だから、k = log2 n。O()表記の場合、対数の底は無視するので、log n。

参考:
再帰定義 (recursive definition)
[初心者向け] プログラムの計算量を求める方法

greedyアルゴリズム(貪欲法)

greedyアルゴリズムは、全部をN回試して、報酬の平均が最も高いやつを選択するというアルゴリズムです。

当たりか外れがでる機械が4個あって、どれが一番当たり率が高いか分からないのでgreedyアルゴリズムでやってみる想定にします。機械はa~dまであってそれぞれ当たり率は下記のとおりとします。

a b c d
0.3 0.6 0.8 0.95

当たりが出たら1点もらえて、外れたら0点とします。

import numpy as np
import random

m = np.array([0.3, 0.6, 0.8, 0.95])
n = 3
total = 50

m_r = np.zeros([m.shape[0]])
m_c = np.zeros([m.shape[0]])
m_h = np.zeros([total])
cost = 0
reward = 0

for i in range(m.shape[0]):
    for _ in range(n):
        m_h[cost] = i
        if random.random() <= m[i]:
            reward += 1
            m_r[i] += 1
        cost += 1
        m_c[i] += 1
for i in range(total - m.shape[0] * n):
    good_m = np.argmax(m_r / m_c)
    m_h[cost] = good_m
    if random.random() <= m[good_m]:
        reward += 1
        m_r[good_m] += 1
    cost += 1
    m_c[good_m] += 1
print(cost)
print(reward)
print(reward / cost)
print(m_h)

結果

50
45
0.9
[ 0.  0.  0.  1.  1.  1.  2.  2.  2.  3.  3.  3.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  3.  3.  3.  3.  3.
  3.  3.  3.  3.  3.  3.  3.  3.  3.  3.  3.  3.  3.  3.]