Python3 – 対数

Pythonで対数を出すには、math.logを使います。import mathで使えるようになります。

import math
print(math.log(2, 10))

0.30102999566398114

これは、10を何乗したら2になるか?です。log102です。
試しに、10を0.30102999566398114乗してみます。

print(pow(10, 0.30102999566398114))

1.9999999999999998

100の0.5乗は、100の平方根なので10です。

print(pow(100, 0.5))

10.0

27の0.3333333333333333333333333333333333333333333乗は、大体27の立方根なので、3に近いはずです。

print(pow(27, 0.3333333333333333333333333333333333333333333))

3.0

おージャスト3.0になった。

print(pow(27, 0.333333333333333))

2.9999999999999964

print(pow(27, 0.3333333333333333))

3.0

少数15桁までだと2.99…になりますが、16桁にするとジャスト3.0になる。

10を底とする対数を常用対数といいます。

対数関数と指数関数は逆関数の関係にあります。逆関数はy=xに対して対称です。指数関数と対数関数のグラフを書いてみます。
y = axの指数関数と、y = logaxの対数関数を書きます。

numpyを使いたいので、math.logではなくnp.logを使います。aを2として、np.log2を使います。例えばnp.log2(4)とすると2になります。2を底とした4の対数を出しています。

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(0.1, 3, 0.1)
y = x
plt.plot(x, y, label='y=x')
y = 2 ** x
plt.plot(x, y, label='y=2**x')
y = np.log2(x)
plt.plot(x, y, label='y=log2 x')
plt.grid(True)
plt.legend()
plt.show()

対称なのか全然これじゃ分からない。

対数関数は、
対数が1より小さいとき、
底が1より大きければ、真数は底より小さい。
底が1より小さければ、真数は底より大きくなる。

対数関数を、底の大きさのパターンで2つ書いてみる。

import math
import numpy as np
import matplotlib.pyplot as plt

#底は1ではない正の数
#底が1より大きい場合
a = 2
#底が1より小さい場合
b = 0.5
#xは真数
x_list = np.arange(0.1, 3, 0.1).tolist()
#yは対数
y_a = []
y_b = []
for x in x_list:
    y_a.append(math.log(x, a))
    y_b.append(math.log(x, b))
plt.plot(x_list, y_a)
plt.plot(x_list, y_b)
plt.grid(True)
plt.show()

2100の桁数を出すには、常用対数を使います。
log102100 = 100log102 = 100 * 0.30102999566398114 = 30.102999566398114
よって31桁です。10の2乗は100なので3桁です。10の3乗は1000で4桁です。2.1乗は100より大きく、1000より小さいので3桁です。30.10乗は31桁です。

Python3 – xのn乗のグラフ(matplotlibのsubplotとアニメーション)

GIFアニメーション

matplotlibのアニメーション作成は2つ種類があって、ArtistAnimationとFuncAnimationとがある。

参考:matplotlib でアニメーションを作る

ArtistAnimation は、あらかじめ全てのグラフを配列の形で用意しておき、それを1枚ずつ流すというものである。FuncAnimationは予め完成したグラフを渡すのではなく、アニメーションの1フレームごとに関数を実行する。データが巨大すぎる場合や潜在的に無限に続く場合に便利。

ArtistAnimation

ArtistAnimationを使って、xの2乗~10乗のグラフを配列に入れて、アニメーションをつくってみる。

サンプルコード

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def xpown(s, e):
    graphs = []
    fig = plt.figure()
    for n in range(s, e + 1):
        x = np.arange(-100, 100, 0.1)
        y = x ** n
        g = plt.plot(x, y, label='y = x ** {}'.format(n))
        graphs.append(g)
    plt.legend()
    ani = animation.ArtistAnimation(fig, graphs, interval=500)
    ani.save('pow1.gif', writer='imagemagick')

if __name__ == '__main__':
    xpown(2, 10)

結果

これだと、10乗のグラフのyが大きすぎて、他がみんな直線みたいになってしまった。

FuncAnimation

FuncAnimationを使ってやってみる。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

x = np.arange(-100, 100, 0.1)
n = 1
max = 10
fig = plt.figure()

def plot_pow(data):
    plt.cla()
    global n
    n = n + 1 if n < max else 2
    y = x ** n
    im = plt.plot(x, y, label='y = x ** {}'.format(n))
    plt.legend()

if __name__ == '__main__':
    ani = animation.FuncAnimation(fig, plot_pow, interval=400, frames=9)
    ani.save("pow2.gif", writer="imagemagick")

結果

指数が偶数ならy軸に対称のグラフになり、奇数なら原点に対象のグラフになります。原点に対象っていうのか。

Subplot

matplotlibで複数のグラフを表示したいときにsubplotが使えます。
参考:[Python]Matplotlibで複数のグラフを描画する方法

サンプルコード

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(-100, 100, 0.1)
fig = plt.figure()
ax = []
for i in range(9):
    y = x ** (i ++ 2)
    ax.append(plt.subplot2grid((3, 3), (i // 3, i % 3)))
    ax[i].set_title('y = x ** {}'.format(i + 2))
    ax[i].plot(x, y)
plt.show()

結果

Python3 – 変数のスコープについて

下記コードのとき結果は1と表示されます。

n = 1
def hoge():
    print(n)
hoge()

下記コードのときエラーになります。

n = 1
def hoge():
    n += 1
    print(n)
hoge()

UnboundLocalError: local variable ‘n’ referenced before assignment

参考:なぜ変数に値があるのに UnboundLocalError が出るのですか?

これは、あるスコープの中で変数に代入を行うとき、その変数はそのスコープに対してローカルになり、外のスコープにある同じ名前の変数を隠すからです。foo の最後の文が x に新しい値を代入しているので、コンパイラはこれをローカル変数であると認識します。その結果、先の print(x) が初期化されていないローカル変数を表示しようとして結果はエラーとなります。

globalだよと宣言すると大丈夫になります。

n = 1
def hoge():
    global n
    n += 1
    print(n)
hoge()

ネストされたスコープだと、nonlocalが使えます。

def hoge():
    n = 1
    def page():
        nonlocal n
        n += 1
        print(n)
    page()
hoge()

Python3 – matplotlibでアニメーションGIFをつくる

matplotlib.animatioを使うとアニメーションがつくれます。自分の環境ではgifで保存しようとしたら、imagemagickがなくてエラーになりました。imagemagickのインストールとかはここに書いてありました。

Imagemagick

http://www.imagemagick.org/script/binary-releases.php#windows

Imagemagickのパスをmatplotlibに設定する必要がある。下記のようにするとmatplotlibが見ている設定ファイルの場所がわかる。

import matplotlib
print(matplotlib.matplotlib_fname())

自分はユーザディレクトリ内のAnaconda3の中の下記だった。

Anaconda3\lib\site-packages\matplotlib\mpl-data\matplotlibrc

matplotlibrcの最後の行らへんに、下記行があるので、ここにImagemagickのパスを入れる。

#animation.convert_path: 'convert' # Path to ImageMagick's convert binary.

コメントを外して、下記のように入力します。入力の際、パスを「’」で囲むと下記のようなエラーになります。

animation.convert_path: D:\ImageMagick-7.0.4-Q16\magick.exe


C:\Users\hoge\Anaconda3\lib\site-packages\matplotlib\animation.py:782: UserWarning: MovieWriter imagemagick unavailable
warnings.warn(“MovieWriter %s unavailable” % writer)
Traceback (most recent call last):
File “C:/Users/hoge/hoge.py”, line 15, in
ani.save(‘sample.gif’, writer=’imagemagick’)
File “C:\Users\hoge\Anaconda3\lib\site-packages\matplotlib\animation.py”, line 810, in save
writer.grab_frame(**savefig_kwargs)
File “C:\Users\hoge\Anaconda3\lib\contextlib.py”, line 66, in __exit__
next(self.gen)
File “C:\Users\hoge\Anaconda3\lib\site-packages\matplotlib\animation.py”, line 196, in saving
self.finish()
File “C:\Users\hoge\Anaconda3\lib\site-packages\matplotlib\animation.py”, line 389, in finish
+ ‘ Try running with –verbose-debug’)
RuntimeError: Error creating movie, return code: 1 Try running with –verbose-debug

サンプルコード

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def anim():
    graphs = []
    fig = plt.figure()
    x = np.arange(-50, 50, 0.1)
    y = x ** 2
    for i in range(50):
        plt.plot(x, y)
        g = plt.plot(i, i ** 2, 'o')
        graphs.append(g)
    ani = animation.ArtistAnimation(fig, graphs, interval=50)
    plt.show()
    ani.save('anim.gif', writer='imagemagick')

if __name__ == '__main__':
    anim()

Python3 – 素因数分解

素因数分解

正の整数 n を素因数分解するための最も単純な方法は、2 から順に √n までの素数で割っていく方法である(Trial division(英語版))。しかし、n が大きくなると、この方法では困難である。

結果
[3, 79, 519507173]

参考:Python Finding Prime Factors

Python3 – timeitを使ってコマンドライン上で関数の実行速度を計る

コマンドラインで実行できるので、プログラムに計測用のコードを書かなくてよくて便利。

timeit — 小さなコード断片の実行時間計測

計測方法

自分の環境だと、実行したいpythonファイル(hoge.py)がおいてあるディレクトリで、下記のようにやると動きます。(hoge.pyのtest関数の実行速度を計測したい場合)

$ python -m timeit 'import hoge' 'hoge.test()'

計測結果の見方

試しに、ここで作成した素数を出すコードの速度を計測してみます。ファイル名は、eratos.pyです。

$ python -m timeit 'import eratos' 'eratos.prime(100000)'

結果

100 loops, best of 3: 17.5 msec per loop

結果の見方がややこしい。
参考:ライブラリ:timeit

-nは、1試行あたりの実行回数。-rは試行回数らしい。下記のように指定できる。

$ python -m timeit -n 10 -r 5 'import eratos' 'eratos.prime(100000)'
10 loops, best of 5: 17.4 msec per loop

上記の場合は、処理を10回連続で実行したときの時間を5回計測し、そのうちの最小時間を返してくれているらしい。デフォルトだと連続で100回実行して、ランダムなタイミング(?)で3回だけ時間を計測して、3回のうち一番時間が短かった結果を返してくれているということかな?まあ何しろこれで実行時間が簡単に計測できるので便利です。

時間の単位は、下記になります。

単位名 単位
ナノ秒 1,000,000,000 nsec(ns) 10の-9乗 s
マイクロ秒 1,000,000 usec(µs) 10の-6乗 s
ミリ秒 1,000 msec(ms) 10の-3乗 s
1 sec(s) 1s

Python3 – 素数

エラトステネスの篩

エラトステネスの篩

結果
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997]

参考:素数と仲良しになる3つのプログラム

フェルマーの小定理

フェルマーの小定理で素数か判定できる。

Pythonを使って高速素数判定をしてみる」に詳しく書いてあった。フェルマーテストというので出せるけど、素数じゃないものも出しちゃう可能性があって、より高性能な方法がある。それもこのサイトに書いてある。

Python3 – Seleniumの使い方

Selenium
http://www.seleniumhq.org/

作業環境

私の環境は、下記です。
Windows10
Python3.5
Firefox 51.0.1
Chrome 56.0
ChromeDriver 2.27

Selenium WebDriver

SeleniumRCはJavascriptでブラウザを操作してたけどセキュリティ的に難しくなってSelenium WebDriverが主流になったそうです。WebDriverはブラウザの拡張機能やOSのネイティブ機能などを利用して操作するそうです。ブラウザのドライバに対して、JSON Wire ProtocolでHTTPリクエストを行うとブラウザが操作できるます。ドライバは各ブラウザ毎に異なります。

ドライバに対するリクエストは、一般的には各種プログラミング言語で作成します。使える言語は、公式だと、Java、Ruby、Python、C#、Javascript、PHP、Perlだそうです。その他非公式であったりするそうです。私はPythonでやってみようかなと思います。

Python用Seleniumモジュールをインストールする

pipでインストールできます。

$ pip install -U selenium

Python用のseleniumのAPIのドキュメントはここです。

chromeドライバのインストールと設定

Chrome Driverはこのページでダウンロードできます。ChromeDriver 2.27のchromedriver_win32.zipをダウンロードしました。展開して、chromedriver.exeを任意のフォルダに配置し、そのフォルダにパスを通しておきます。

Fireboxドライバのインストールと設定

Firefoxはブラウザをインストールしてあればドライバが使えると色々なところに書いてはあったのですが、Firefoxの新しいバージョンではそうでもないらしく、実際自分もエラーになってしまいました。下記のようなエラーが出る場合、ドライバが有効になっていない可能性があります。

AttributeError: 'Service' object has no attribute 'process'

上記のような場合は、ここからドライバをダウンロードしてみます。私は、geckodriver-v0.14.0-win64.zipをダウンロードして展開して、上記のchromeと同様に任意のフォルダにgeckodriver.exeを配置し、そのフォルダにパスを通しました。これで、Firefoxも動作するようになりました。

PythonでSeleniumを操作してみる

ブラウザを起動してURLを開く

from selenium import webdriver

#Firefox
fire = webdriver.Firefox()
fire.get('http://google.com')

#Chrome
chro = webdriver.Chrome()
chro.get('http://google.com')

要素の選択、入力、クリック

Googleに検索ワードを入力して検索ボタンを押すことで検索結果ページに遷移させてみます。

from selenium import webdriver

#Firefox
fire = webdriver.Firefox()
fire.get('http://google.com')
input = fire.find_element_by_id('lst-ib')
input.send_keys('Selenium Python3')
btn = fire.find_element_by_name('btnK')
btn.click()

#Chrome
chro = webdriver.Chrome()
chro.get('http://google.com')
input = chro.find_element_by_id('lst-ib')
input.send_keys('Selenium Python3')

chromeの場合は、検索ワードを入力したら勝手に検索結果のページに遷移したので、btnを取得・クリックしていません。

ログイン

上記と同じですが、ログインもしてみます。Twitterにログインします。

from selenium import webdriver

tw_user = 'hoge@hoge.com'
tw_pass = 'passwordpassword'

#Firefox
fire = webdriver.Firefox()
fire.get('https://twitter.com')
login = fire.find_element_by_class_name("StreamsLogin")
login.click()
input_user = fire.find_element_by_name('session[username_or_email]')
input_pass = fire.find_element_by_name('session[password]')
input_user.send_keys(tw_user)
input_pass.send_keys(tw_pass)
btn = fire.find_element_by_class_name('submit')
btn.click()

#Chrome
chro = webdriver.Chrome()
chro.get('https://twitter.com')
login = chro.find_element_by_class_name("StreamsLogin")
login.click()
input_user = chro.find_element_by_name('session[username_or_email]')
input_pass = chro.find_element_by_name('session[password]')
input_user.send_keys(tw_user)
input_pass.send_keys(tw_pass)
btn = chro.find_element_by_class_name('submit')
btn.click()

問題なくログインできるのですが、getでURLを開くと全て読み込むまで次に進まないようで、ログインボタンをクリックするまでにやたらと待たせれる点がいやです。解決方法があるか今度確認します。

参考:
第1回 Seleniumの仕組み
Seleniumとは―読み方、意味、WebDriverの使い方
Selenium Web DriverのPythonバインディングを使ってブラウザを操作する
5分でわかるSelenium IDEの使い方
Python で Selenium WebDriver を使ったブラウザテストをする方法
Webブラウザの自動操作 (Selenium with Rubyの実例集)

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)
[初心者向け] プログラムの計算量を求める方法