2018-08-20

並列計算の効果検証用プログラム

並列化の試行やJITの効果テストなどのために，Forをたくさん回す計算例が必要になる。そのような例としてモンテカルロ法を置いておく。

ToDo 並列計算のテスト用としては，モンテカルロ法のような単純に並列化される(データ共有が無い)ものだけでなく，適時データ共有が必要となる例を別途置いておく。

モンテカルロ法による円周率計算

モンテカルロ法の例題でよく見る計算。[0,1]の範囲の乱数を生成して，その点が半径1の円の中に入るか否か？(点が円の中か否かは，点の距離dがsqrt(x²+y²)で計算できるから，円周率piを知らなくても判定できる)をランダムに試行して，N回サンプリングすると，円の中に入る比率が1:pi/4になることから推定する(円の面積がpirrなのは知っている)

# monte_carlo.py
import random

def calc_pi(iter):
    cnt = 0
    for _ in range(iter):
        x, y = random.random(), random.random()
        if x*x + y*y < 1:
            cnt += 1
    return (cnt*4.0)/iter

全く同じ書き方でC言語でも。

#include <stdio.h>
#include <stdlib.h>

double calc_pi(int iter)
{
  int cnt=0;
  double x, y;
  
  for (int i=0; i<iter; i++) {
    x = (double)rand()/RAND_MAX;
    y = (double)rand()/RAND_MAX;
    if ((x*x + y*y) < 1)
      cnt++;
  }
  return cnt*4.0/iter;
}

int main(int argc, char** argv) {
  int iter;
  double pi;
  
  iter = atoi(argv[1]);
  pi = calc_pi(iter);
  printf("pi:%f", pi);

  return 0;
}

ちなみに，生のPythonは20倍くらい遅い。

# Python
$ time python3 monte_carlo.py 10000000 # 10M
real 0m3.838s                          # 0.4us/iter

# C(Clang, -O3)
$ time ./a.out 10000000                # 10M
real 0m0.214s                          # 0.02us/iter

2018-08-18

モンテカルロ法1 (モンテカルロ積分)

統計・機械学習 python

Rによるモンテカルロ法入門のメモ，その１。

イントロ

統計的推論では最適化と積分の計算，特に実用上は数値計算が必要になる。
数値積分の手法としてモンテカルロ積分がある。シンプソン法などの区分法は多次元の場合に区分を細かくしようとすると計算量が増加してしまうので適用が困難である。そのような多次元の場合に対してモンテカルロ積分は有用。

モンテカルロ 積分

モンテカルロ積分の考え方は非常に単純で，下記の期待値計算の積分を対象の確率分布からのサンプリングによる平均値で近似する方法。
例えば，f(x)が確率密度関数で，xがf(x)に従う時のh(x)の期待値を計算したいときに，

$E_f(h(x)) = \int_{x}f(x)dx$

を

$\frac{1}{m}\sum_{i=1}^{m}h(x_j)$

で近似する。

同様に考えて，定積分も対象の範囲に対する一様乱数を生成して，その時のh(x)の値を求めて平均取れば，その区間の平均値となるので，積分区間を掛けて上げれば積分値が得られる。(b-a倍するのを忘れないように)

import numpy as np
from scipy import integrate
import math

f = lambda x: math.sin(x)/x  # 被積分関数

def mc_integrate(f, a, b, n=10000):
    def scaling(x):
        '''
        np.random.randで生成する[0,1)の一様乱数を[a,b)に変換する関数
        '''
        return (b-a)*x + a

    xs = np.random.rand(n)
    xs = np.vectorize(scaling)(xs)
    return (b-a) * np.vectorize(f)(xs).sum()/n


def main():
    v = mc_integrate(f, 0, 0.0, 2.0*math.pi)
    
    # Scipyを使って検算
    v_val = integrate.quad(f, 0.0, 2.0*math.pi)[0]
    print('Scipy, v={:.5}, MC, v={:.5}'.format(v_val, v))

if __name__ == '__main__':
    main()

2018-08-16

パフォーマンス計測 in Python

python

timeモジュール(time.perf_counter)

使いどころは特定のコードブロックに対する計測。精度はマイクロ秒程度。精度を求める場合はこれ。(逆にcProfile/line_profileなどのプロファイリングは基本的にクロック程度の精度,1/100秒程度しかなく，ボトルネックの発見が目的と考えてよいのかもしれない。)
使い方は開始位置と終了位置でtime.time()ならエポック秒(UNIX時間)を，time.perf_counter()ならパフォーマンスカウンタの値を取得して差分として計算する。
Pythonのドキュメントに依ると，time.time()システムによっては1秒程度しか精度が無い場合もある。しかし，UNIXにおいてはgettimeofday()を使って時刻を取得するようなので，マイクロ秒単位程度の精度があるようだ。 # Web上にはtime.time()が1/60~1/100秒程度の精度の記述が見られるが，Python3のドキュメントにはその記述は見つけられなかった。
より明確にパフォーマンスカウンタ値を返す関数としてtime.perf_counter()がある。time.time()と同じような気もするけど，使い方は同じであるし，無難にこちらを使っておく。 # 実際に，time.get_clock_info('time')とtime.get_clock_info('perf_counter')のresolutionはどちらも1e-09。

import time
start_time = time.perf_counter()
# some heavy calculation
# ...
end_time =  time.perf_counter()
print('Elapsed time:{:.7}'.format(end_time - start_time))

cProfile

使いどころは関数ごとの処理時間やコール回数を計測したい場合。精度は1/100秒程度。
使い方はスクリプト実行時に，cProfileをインポートして実行するだけ。
表示の意味をまとめておく。
- ncalls: 呼び出し回数。再帰の場合が別途表示されていて，"再帰含む場合"/”再帰を含まない場合”の２つが出力される。
- tottime: その関数で消費された時間。その関数が呼び出す関数の消費時間は含まない。
- cumtime: その関数で消費された時間。その関数が呼び出す関数の消費時間は含む。
'-s'オプションでソートして出力することがほとんど。
- time: tottimeでソート
- cumulative: cumtimeでソート
- calls: ncallsでソート

# tak.py (Takeuchi Function) is tested as an example.
hoge@foo:~/$ python3 -m cProfile -s time tak.py 

         12604868 function calls (8 primitive calls) in 3.457 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
12604861/1    3.457    0.000    3.457    3.457 tak.py:3(tak)
        1    0.000    0.000    0.000    0.000 {built-in method builtins.print}
        1    0.000    0.000    3.457    3.457 tak.py:15(main)
        1    0.000    0.000    3.457    3.457 tak.py:1(<module>)
        2    0.000    0.000    0.000    0.000 {built-in method time.perf_counter}
        1    0.000    0.000    3.457    3.457 {built-in method builtins.exec}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

line_profiler

使いどころは行単位のプロファイルを取りたい場合。cProfileであたりを関数単位で付けて，計測したい関数をline_profileで指定(デコレータ)する。
いくつかの使い方があるが，コードを修正してしまうので面倒ではあるが単純なのは関数に@profileデコレータを付けて，kerfnprof.py(line_profileと一緒にインストールされる)を使う。

# hoge.py
@profile
def f():
    # heavy calculation    
    ...
    ...
    ...
    return

pipでインストールして，かつline_profilerをgit cloneしていて例があるが，kernprof.py自体はローカルにあるので不要。

$ pip3 install line_profiler --user # ~/.local/lib/python3.6/site-packages/にインストールされる。
$ python3 ~/.local/lib/python3.6/site-packages/kernprof.py -l -v hoge.py

2018-07-28

訓練・検証・テストの３つにデータを分割する理由

統計・機械学習

MLの教科書とか記事を読んでいると，データを訓練，検証, テストの３つに分けている場合が多い。訓練と検証の２つじゃだめなの？とついつい勘違いしてしまうのでメモ。

テストはなんのため？

ここの勘違いが全てなんだと思う。
テストの目的は実運用時の性能を推定すること。勘違いしがちなのは，モデル選択にテストデータでの評価結果を使うこと。
繰り返すが，テストは性能を推定するだけ，モデル選択するためじゃない。

テストデータと検証データの使い方

MLのモデルはハイパーパラメタがある場合がほとんど。
よくある勘違いは訓練データで種々のモデル(モデル自体やハイパーパラメタが違うモデル)を訓練し，テストデータで性能を測定し，最も性能が良いものを採用してしまう方法。これはあくまでも，そのテストデータに対して最も良いモデルを選んだだけであって，その値を持って実運用時の性能を推定することはできない。
くどいけれど，あくまでも，テストはテストであって，選択基準にしてはならない。
適切な使い方は，訓練データでパラメタを学習して，検証データでモデル選択をして，テストデータで性能を測定する。
だから，名前のままなんだ。訓練データで訓練して，検証データで検証してみて，テストデータでテストする。
そうか，検証とテストという言葉が似すぎているのか！検証は自分側に裁量がまだあるが，テストはもはや裁量権を失っているようなイメージなのかな。モデルを検証しましたので，テストお願いします，みたいな感じなのかな。

補足(クロスバリデーション)

データを分割すると，貴重なデータ数が目減りして見える。
より有効的に訓練と検証をする方法がクロスバリデーション。
データを訓練データとテストデータに分けて，訓練データに対して交差検証をする方法。これはデータをうまく使いつつ，良いモデルを選択できる。

2018-06-16

TkinterによるPythonのGUIプログラミング

python

趣旨

たまに使おうとするとウェブを漁りまくることになるので，注意点だけまとめておく。

ドキュメント

tkとttk

テーマを設定すれば一括してウィジットの見た目を変更できるのがttk。
多くのウィジットが同名のウィジットがあるけれど，canvas, menu, messageはttkにはないので注意する。

ttkテーマの設定

style = ttk.Style()
available_themes = style.theme_names() # 利用可能テーマの確認
style.theme_use('clam')  # テーマを指定する。

ウィジット配置

３つのマネージャーが用意されているが，座標指定のplaceはほぼ使わない。１行または1列に並ぶだけの場合にはpackマネージャーがが有用。ただし，多くの場合はpackでは済まないので，gridを使うことが基本になる。グリッドと言っても，columnspan/rowspanで複数行列に跨ることは可能なので，大概は問題ない。
適切な配置をするためにexpand/anchor/fill/padxy/ipadxyのオプションについて理解しておく。
- packのオプション
  - expand:親ウィジットの拡大縮小に応じてリサイズするかどうか。
  - anchor:左寄せ，上寄せなど，スペースに空きがある場合にどこに寄せるか。東西南北(EWSN)により8方向の中から指定する。例えば，右上は'NE'
  - fill:スペースがある場合にウィジットを拡大する。x, y, bothでどちらに拡大するかを指定する。
  - padxy:外側のスペース。
  - ipadxy:内側のスペース。
- gridのオプション
  - row/column:行列を指定。0から開始。
  - row/columnspan:ぶち抜く場合には幅を指定する。
  - padxy,ipadxy:packと同じ。
  - sticky:packのfill+anchorのオプション。sticky=tk.Eだと右寄せ。fillする場合は'+'でつなぐ。例えば，横に伸ばす場合にはsticky=tk.W+tk.E，全体はsticky=tk.W+tk.E+tk.N+tk.S，みたいな要領。

ウィンドウの拡大縮小

親フレームをexpand, fill, anchorでrootと連動して拡大縮小するようにpackしておく。
子ウィジットはgridで設置する際に，拡大縮小したいウィジットをのcolumnconfigure, rowconfigureのweightを1にする。デフォルトでは0になっている。よって，例えば，画面サイズに依存せず同じ大きさで良いボタンなどはデフォルトのまま(0)にしておいて，textやentryなどを1にしておくと，そのウィジットだけ拡大縮小される。

parent_frame.columnconfigure(1, weight=1)  # 列方向は1列目が拡大縮小される。
parent_frame.rowconfigure(2, weight=1)     # 行方向は2行目が拡大縮小される。

値の取得方法

textやentryなどの値は，ウィジット生成時にtextvariableにtk.Variableを指定しておけば自動的にその変数に値が格納される。
tk.Variableは型ごとにIntVar, StringVar, BooleanVarなどが用意されていて，変換もしてくれる。
値を取得する場合にはget()，設定する場合はset()。生の変数(StringVarなど)はそのままでは値ではないので注意する。

entry_var = Tk.StringVar()
entry = Tk.Entry(self, textvariable=entry_var)

コールバックの設定(command/bind)について

ボタンを押した場合などの処理をcommandにて設定することが可能である。
注意としては，commandに渡す関数は引数を取れないことである。よって，関数生成時に引数をクロージャで包み込んでおいて，それをcommandに渡すようにする。

textウィジット

文字列の書き込みはinsertを使う。第一引数で書き込み開始位置を，第２引数で文字列を渡す。
現在の最後尾がtk.ENDで指定できる。
書き込むだけだと，表示位置は動いてくれないので，seeで表示位置も動かす。

txt_wdgt.insert(tk.END, (''hogehoge))  # 文字列はなぜかカッコで囲む。
txt_wdgt.see(tk.END)

イベント処理

キーイベント

ウィジット毎にバインドできる。当然rootウィンドウにもバインドできる。
callbackは１引数としてevent情報を取る関数。(クラス内ならselfも入れて2引数)
イベントには下記のプロパティが設定されている。
- x,y: マウス位置
- num: マウスボタンの番号(1:左，2:真ん中，3:右)
- time: イベント時刻
- char: キー文字
- keysym: キーに対する名前

# widget.bind(<Event Sequence>, callback)の書式
# Event Sequence = <modifier>-<modifier>-<type>-<detail>の書式
#  i.e.) <Control-Shift-A>, <Buttun-1-KeyPress>, <Button-1-KeyPress-Motion>
root.bind('<Right>', move_right)  # arrow key:[Right, Left, Up, Down]
root.bind('<Escape>', move_right)  # arrow key:[Right, Left, Up, Down]

def move_right(e):
    print(e.x)

ウィジット情報の取得と設定

cgetで取得し，configureで設定する。ただし，canvasの項に記載の通り，ウィジットの大きさに関してはcgetは初期値の設定情報のようで，リアルタイムな値はwinfo_xxxを確認する。
ウィジットはそれぞれwinfo_xxx(widget information)を持っている。
xxxとしてはgeometory, width, height, x, y, rootx, rootyがある。ｘ，ｙは親ウィジット内での位置，rootx,yはディスプレイ上の位置。
ウィジットが生成された直後に出力してもwinfo_xxxが正しく取れない場合があるので注意。その場合はupdate_idletasksをして実行させるか，他の処理で時間を消費してからアクセスする。

Canvasウィジット

色の指定

'#RGB'で各色16進(0〜F)で指定する。各色２バイトにも出来る(0〜FF)。

オブジェクト操作

Canvas内のオブジェクト(Rectangle)などはそれぞれIDを持っており，move, deleteなどIDを介して操作が可能。
scale関数もあるので，ウィンドウのリサイズに応じてオブジェクトをリサイズすることも可能。

ウィンドウのリサイズに伴うオブジェクトのリサイズ

Canvasウィジット自体はpackならexpandに，gridならcolumnconfigureなどで自動でリサイズできる。
ただし，Canvas内の描画オブジェクトはリサイズされないので別途対処が必要。
ウィンドウ(rootウィジット)のりサイズはにて検出する。そして，検出イベントのコールバックとして各図形をリサイズする。
キャンバスのりサイズ後の大きさはcget(configure get)では得られない。cgetでは初期設定値が返ってくる。代わりに，winfo_width(), winfo_height()にて取得可能。
ウィジットが生成された直後に出力してもwinfo_xxxが正しく取れない場合があるので注意。その場合はupdate_idletasksをして実行させるか，他の処理で時間を消費してからアクセスする。

root.bind('<Configure>', rescale_objects)

クラスとして書く

フラットに書かれている例が多いが，クラス化しておくと見通しも良いので，クラスの書き方テンプレートをもとに拡張する。
フラットに書いてあるものと基本は同じ。

import tkinter as tk

class Application(tk.Frame):
    def __init__(self, master=None):
        tk.Frame.__init__(self, master)
        master.title("My application")
        self.pack(expand=1, fill=tk.BOTH, anchor=tk.NW)
        self.create_widgets()

    def create_widgets(self):
        self.label = tk.Label(self, text='Input file')
        self.label.grid(column=0, row=0)
        
root = tk.Tk()
app = Application(master=root)
app.mainloop()

その他

バージョン確認

Python3.6付属のTkinter8.0以降は見た目がきれいになっている。

$ python -c "import tkinter;print(tkinter.TkVersion)"

ウィジットのプロパティへのアクセス

各ウィジットのプロパティ一覧はconfigure()にて取得できる。

w = tk.WidgetXXX()
print(w.configure())

スクロールについて

textなどにおいて，スクロールバーは付けなくてもどんどん延長は自動でしてくれる。記述が煩雑で見た目もきれいではないので付けなくても良い気がする。

2018-06-12

MatplotlibのimshowでRGBデータを表示

python

imshowは２Dデータだけでなく，RGBデータを含む(H, W, CH)のデータを表示することができる。
その際に，floatかintかによって，値の範囲が異なるようだ(調べきれていないので自信が無いけど)
float32のときは，[0,1]の実数値，int32などの整数値のときは[0,255]の256値でやるとうまく行く。
vmin, vmaxを設定してもうまく認識されなかった。

# uintの場合
import numpy as np
import matplotlib.pyplot as plt
img = np.array([ [[255,0,0], [0, 255,0], [0, 0, 255]],
                 [[255,255,0], [0, 255,255], [255, 0, 255]],
                 [[0,0,0], [128, 255,128], [255, 255, 255]]])
plt.imshow(img, vmin=0, vmax=255, interpolation='none')
plt.show()

f:id:nobUnaga:20180612234017p:plain

# floatの場合
import numpy as np
import matplotlib.pyplot as plt
img = np.random.randint(0, 2, (5,5,3)).astype(np.float32)
plt.imshow(img, vmin=0, vmax=1, interpolation='none')
plt.show()

f:id:nobUnaga:20180612234009p:plain

2018-06-12

Pythonのパッケージ構成とimport文

python

パッケージについて

PyPIとかへの登録とかはおいておいて，自作パッケージをimportしたり，他のプログラムから使おうとするときの方法を整理しておく。下記に従えば，最悪，mypackage以下を開発中のパッケージにコピーすれば，そのまま使える。

ディレクトリ構成

git管理ディレクトリMyPackageに以下の構成にてディレクトリを作成する。
パッケージ名はmypackage(全部小文字)とし，すべてMyPackage/mypackage以下に関連ソースコード，データを配置する。基本的にはこのディレクトリがパッケージのトップディレクトリであると考える。
パッケージの中で静的データを使う場合(例えばゲームのマップファイルなど)は，mypackage以下に置く。間違ってMyPackage以下に置かない。繰り返すが，トップディレクトリはmypackage。
- サンプル実行コードてテストコード，ドキュメントはMyPackage以下に配置する。これらはpackageのインストール自体では入らない。

MyPackage
    - mypakcage  # ソースコードをインストールするディレクトリ
        - hoge.py
        - foo
            - bar.py
        - data
            - data1.dat
            - data2.dat
    - sample
        - sample.py
    - doc
    - test

各ファイルの注意

パッケージ内のimport文

基本的に相対インポートは書かない。絶対インポートを書く。(相対インポートよりもわかりやすいと思う。ディレクトリ構成がコロコロ変わると厄介だけど)
例えば，bar.pyをhoge.pyがインポートする場合は下記のように書く。

# hoge.py
from mypackage.foo import bar  # Good
from foo import bar # <= Bad

sample/test内のimport文

sample/sample.pyは兄弟ディレクトリのmypakcage以下のファイルをそのままではインポートできない。よって，pathに追加する。
その際にはfileがそのファイルのパスを表すこと活用して，親ディレクトリをパスに追加する。os.path.dirnameを２回適用することで，親ディレクトリの絶対パスを求めている。".."を使うと実行ディレクトリに依存するのでその記述は避ける(Pythonは実行ディレクトリをカレントディレクトリとする)

import os
import sys
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))

静的データのパス

マップデータやイメージなどの静的データをロードする場合も同様に，fileがそのファイルのパスを表すことを利用して，パッケージが配置される絶対パスからパスを指定する。

# hoge.pyからdata/data1.datをロードする場合
import os
dir_path = os.path.dirname(os.path.abspath(__file__))
data = dir_path + '/data/data1.dat'