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

Python3 – 動的計画法(フィボナッチ数列)

動的計画法は、分割統治法メモ化を合わせた方法のことらしい。分割統治法は、問題を細分化して、細かい部分を順に解いていくことで全体を解明するようなことの総称らしい。

分割統治法は、コード的には下記のようになり、再帰することになる。

function conquer(x)
  if xは簡単にconquerできる then
    return conquerの結果
  else
    (x1, x2, ...) := divide(x)     // 複数の小さな副問題に分割
    (y1, y2, ...) := (conquer(x1), conquer(x2), ...)  // 再帰
    return synthesize(y1, y2, ...)  // 各副問題の解を合成(最大値を選ぶ、等)

メモ化とは、再帰の際に同じ副問題を複数回解いてしまうことによる効率低下を防ぐために、副問題の結果を保存して再利用することらしい。キャッシュみたいな感じ。

下記のフィボナッチ数列を表示するプログラムは、動的計画法の具体的な例らしい。

int fib(unsigned int n) {
    int memo[1000] = {0, 1}, i;
    for (i = 2; i <= n; i++) {
        memo[i] = memo[i - 1] + memo[i - 2];
    }
    return memo[n];
}

Python3だと下記のような感じ。

def fib(n):
    memo = [0] * n
    memo[0:2] = [0, 1]
    for i in range(2, n):
        memo[i] = memo[i - 1] + memo[i - 2]
    return memo[0:n]
print(fib(100))

ちなみに、ネットでフィボナッチ数列のPython版しらべたら下記がでてきたけど、これはまさしく動的計画法じゃない例だと思った。毎回メモらずに計算しているようだ。

def fib2(n, a=1, b=0):
    return b if n < 1 else fib2(n - 1, a + b, a)
print([fib2(i) for i in range(100)])

せっかくだから時間図ってみる。

import time

def fib(n):
    memo = [0] * n
    memo[0:2] = [0, 1]
    for i in range(2, n):
        memo[i] = memo[i - 1] + memo[i - 2]
    return memo[0:n]
start = time.time()
print(fib(900))
print(time.time() - start)

def fib2(n, a=1, b=0):
    return b if n < 1 else fib2(n - 1, a + b, a)
start = time.time()
print([fib2(i) for i in range(900)])
print(time.time() - start)

結果

[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89,...
0.001961946487426758
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89,...
0.0782785415649414

おー40倍速いってことか。さすが動的計画法だ。

下記のような方法もある。

#include <stdbool.h>

int memo[1000] = {0, 1};
bool in_memo[1000] = {true, true};

int fib(unsigned int n) {
    if (!in_memo[n]) {
        memo[n] = fib(n - 1) + fib(n - 2);
        in_memo[n] = true;
    }
    return memo[n];
}

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

Python3 – random

randomを試してみます。

下記をやってみます。100回randintを0~100まででやってみます。

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

arr = np.zeros([100])
for _ in range(100):
    a = random.randint(1, 100)
    arr[a - 1] += 1
x = np.arange(0, 100, 1)
plt.plot(x, arr)
plt.savefig('./img/rand.jpg')

結果

10000回やってみます。

100000回やってみます。

100万回やってみます。

seed設定してみます。

random.seed(123)
print(random.random())
print(random.randint(0, 100))

何回やっても、下記でした。

0.052363598850944326
11
random.seed(123)
print(random.random())
print(random.randint(0, 100))

a = random.random()
print(a)
b = 0
cnt = 0
for i in range(5):
    b = random.random()
    if a != b:
        cnt += 1
        print(b)
print(cnt)


a = random.random()
print(a)
b = 0
cnt = 0
for i in range(5):
    b = random.random()
    if a != b:
        cnt += 1
        print(b)
print(cnt)

何回やっても下記でした。

0.052363598850944326
11
0.7689563885870707
0.26655381253138355
0.8385035164743577
0.8748257461052661
0.3791243324890884
0.5623187149479814
5
0.34079822787707015
0.05190165459858176
0.1350574593038607
0.5609620089366762
0.7016897890356433
0.16377684475236043
5

seedが同じだと、プログラムを実行する度に、まったく同じ数字が同じ順番ででてくる。

60%の確率で正解を出す機械をつくってみます。
random.random()が0.6以内だったら当たりにします。

per = 0.6
cnt = 5
results = np.zeros([100])

for i in range(100):
    ok = 0
    for _ in range(cnt):
        if random.random() <= per: ok += 1
    results[i] = ok / cnt
print(results)
print(results.mean())

結果

[ 0.2  0.6  0.8  0.8  0.6  0.4  0.8  0.6  1.   0.4  1.   0.4  0.8  0.6  0.4
  0.8  0.6  0.6  0.6  1.   0.2  0.4  0.2  0.6  0.6  0.6  0.8  0.2  0.6  0.8
  0.8  1.   0.6  0.6  0.4  0.4  0.4  0.6  0.6  0.8  1.   0.8  0.8  0.6  0.8
  1.   0.4  0.6  0.6  0.8  0.8  0.6  0.8  1.   0.6  0.2  0.2  0.6  0.8  0.6
  0.4  0.6  0.2  0.8  0.6  0.4  0.4  0.4  0.6  1.   0.4  0.8  0.4  0.6  0.6
  0.6  0.4  0.4  0.6  0.6  0.8  0.4  0.6  0.2  0.6  0.6  0.8  0.8  0.6  0.4
  1.   0.2  0.6  0.6  0.2  0.8  0.4  0.6  0.   0.6]
0.594

手書き文字を作れるJavascriptをつくってTensorFlowで予測させてみた(2)

この前、「手書き文字を作れるJavascriptをつくってTensorFlowで予測させてみた」という投稿でブラウザ上で手書きした文字画像を、MNISTで訓練したモデルで予測してみましたが、ものすごく精度が悪かったです。今回改めて、CNNを使ってやってみたらかなり精度が上がりました。何%か測ったりしてませんが、自分の手書きだと90%は超える感じでした。やっぱりCNNはすごいなーと思いました。でももしかしたら前回のものにミスがあり、CNNではなくても精度は本当はもっと高い可能性はあります。

もうちょっとやるとしたら、文字を画像の中心に適度な大きさで書く必要があり、例えば右上に小さく2と書いても認識されません。あとは、現在はMNISTに合わせて、手書き文字画像も背景黒、文字色白で作成するように固定していますが、これらの色を変えても認識するようにしたいです。今度やってみます。

Github

https://github.com/endoyuta/mnist_test

index.html

<html>
<head>
<title>MNIST TEST</title>
</head>
<body>
<h1>MNIST TEST</h1>
<canvas id="canvas1" width="400" height="400" style="border: 1px solid #999;"></canvas><br><br>
<input id="clear" type="button" value="Clear" onclick="canvasClear();">
<input id="submit" type="button" value="Submit" onclick="saveImg();"><br><br>
<img id="preview"><span id="answer"></span>
<script src="http://code.jquery.com/jquery-3.1.1.min.js"></script>
<script>
var url = 'http://127.0.0.1:8000/cgi-bin/mnist.py';
var lineWidth = 40;
var lineColor = '#ffffff';
var imgW = imgH = 28;

var canvas = document.getElementById('canvas1');
var ctx = canvas.getContext('2d');
var cleft = canvas.getBoundingClientRect().left;
var ctop = canvas.getBoundingClientRect().top;
var mouseX = mouseY = null;

canvasClear();
canvas.addEventListener('mousemove', mmove, false);
canvas.addEventListener('mousedown', mdown, false);
canvas.addEventListener('mouseup', mouseInit, false);
canvas.addEventListener('mouseout', mouseInit, false); 

function mmove(e){
    if (e.buttons == 1 || e.witch == 1) {
        draw(e.clientX - cleft, e.clientY - ctop);
    };
}

function mdown(e){
    draw(e.clientX - cleft, e.clientY - ctop);
}

function draw(x, y){
    ctx.beginPath();
    if(mouseX === null) ctx.moveTo(x, y);
    else ctx.moveTo(mouseX, mouseY);
    ctx.lineTo(x, y);
    ctx.lineCap = "round";
    ctx.lineWidth = lineWidth;
    ctx.strokeStyle = lineColor;
    ctx.stroke();
    mouseX = x;
    mouseY = y;
}

function mouseInit(){
    mouseX = mouseY = null;
}

function canvasClear(){
    ctx.clearRect(0, 0, canvas.width, canvas.height);
    ctx.fillStyle = '#000000';
    ctx.fillRect(0, 0, canvas.width, canvas.height);
    $('#preview').attr('src', '');
    $('#answer').empty();
}

function toImg(){
    var tmp = document.createElement('canvas');
    tmp.width = imgW;
    tmp.height = imgH;
    var tmpCtx = tmp.getContext('2d');
    tmpCtx.drawImage(canvas, 0, 0, canvas.width, canvas.height, 0, 0, imgW, imgH);
    var img = tmp.toDataURL('image/jpeg');
    $('#preview').attr('src', img)
    return img;
}

function saveImg(){
    var img = toImg();
    console.log(img);
    $.ajax({
        url: url,
        type: 'POST',
        data: {img: img},
        dataType: 'json',
        success: function(data){
            if(data.status){
                $('#answer').html('は、' + data.num + 'です');
            }else{
                $('#answer').html('は、分かりません');
            }
        },
    });
}
</script>
</body>
</html>

cgi-bin/mnist.py

#!/usr/bin/env python

import sys
import os
import cgi
import json
import cgitb
cgitb.enable()

from PIL import Image
import numpy as np
from io import BytesIO
from binascii import a2b_base64

import mytensor

print('Content-Type: text/json; charset=utf-8')
print()

if os.environ['REQUEST_METHOD'] == 'POST':
    data = cgi.FieldStorage()
    img_str = data.getvalue('img', None)
    if img_str:
        b64_str = img_str.split(',')[1]
        img = Image.open(BytesIO(a2b_base64(b64_str))).convert('L')
        img_arr = np.array(img).reshape(1, -1)
        img_arr = img_arr / 255
        result = mytensor.predict(img_arr)
        print(json.dumps({'status': True, 'num': result}))
        sys.exit()
print(json.dumps({'status': False, 'num': False}))

cgi-bin/mytensor.py

import tensorflow as tf
import numpy as np
from tensorflow.examples.tutorials.mnist import input_data

def _weight_variable(shape):
    initial = tf.truncated_normal(shape, stddev=0.1)
    return tf.Variable(initial)

def _bias_variable(shape):
    initial = tf.constant(0.1, shape=shape)
    return tf.Variable(initial)

def _conv2d(x, W):
    return tf.nn.conv2d(x, W, strides=[1, 1, 1, 1], padding='SAME')

def interface():
    x = tf.placeholder(tf.float32, shape=[None, 784])
    y_ = tf.placeholder(tf.float32, shape=[None, 10])

    x_image = tf.reshape(x, [-1, 28, 28, 1])
    W_conv1 = _weight_variable([5, 5, 1, 32])
    b_conv1 = _bias_variable([32])
    h_conv1 = tf.nn.relu(_conv2d(x_image, W_conv1) + b_conv1)
    h_pool1 = tf.nn.max_pool(h_conv1, ksize=[1, 2, 2, 1], strides=[1, 2, 2, 1], padding='SAME')

    W_conv2 = _weight_variable([5, 5, 32, 64])
    b_conv2 = _bias_variable([64])
    h_conv2 = tf.nn.relu(_conv2d(h_pool1, W_conv2) + b_conv2)
    h_pool2 = tf.nn.max_pool(h_conv2, ksize=[1, 2, 2, 1], strides=[1, 2, 2, 1], padding='SAME')

    W_fc1 = _weight_variable([7 * 7 * 64, 1024])
    b_fc1 = _bias_variable([1024])
    h_pool2_flat = tf.reshape(h_pool2, [-1, 7 * 7 * 64])
    h_fc1 = tf.nn.relu(tf.matmul(h_pool2_flat, W_fc1) + b_fc1)

    keep_prob = tf.placeholder(tf.float32)
    h_fc1_drop = tf.nn.dropout(h_fc1, keep_prob)

    W_fc2 = _weight_variable([1024, 10])
    b_fc2 = _bias_variable([10])
    y_conv = tf.matmul(h_fc1_drop, W_fc2) + b_fc2

    cross_entropy = tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits(y_conv, y_))
    train_step = tf.train.AdamOptimizer(1e-4).minimize(cross_entropy)
    correct_prediction = tf.equal(tf.argmax(y_conv, 1), tf.argmax(y_, 1))
    accuracy = tf.reduce_mean(tf.cast(correct_prediction, tf.float32))
    saver = tf.train.Saver()

    class CNNModel():
        pass
    model = CNNModel()
    model.x = x
    model.y_ = y_
    model.keep_prob = keep_prob
    model.y_conv = y_conv
    model.train_step = train_step
    model.accuracy = accuracy
    model.saver = saver
    return model

def predict(img):
    ckpt = tf.train.get_checkpoint_state('./cgi-bin/ckpt')
    if not ckpt: return False
    m = interface()
    with tf.Session() as sess:
        m.saver.restore(sess, ckpt.model_checkpoint_path)
        result = sess.run(m.y_conv, feed_dict={m.x: img, m.keep_prob:1.0})
        return int(np.argmax(result))

def train():
    if tf.train.get_checkpoint_state('./ckpt'):
        print('train ok')
        return
    mnist = input_data.read_data_sets('./mnist', one_hot=True, dtype=tf.float32)
    m = interface()
    with tf.Session() as sess:
        sess.run(tf.global_variables_initializer())
        for i in range(20000):
            batch = mnist.train.next_batch(100)
            m.train_step.run(feed_dict={m.x: batch[0], m.y_: batch[1], m.keep_prob: 0.5})
            if i % 100 == 0:
                train_accuracy = m.accuracy.eval(feed_dict={m.x:batch[0], m.y_: batch[1], m.keep_prob: 1.0})
                print("step %d, training accuracy %g"%(i, train_accuracy))
        m.saver.save(sess, './ckpt/model.ckpt')
    print('train ok')

if __name__ == '__main__':
    train()