こんにちは、Gunosyのエンジニアの粟飯原です。Gunosyでは、主に広告配信サーバー全般の開発運用インフラを行いつつ、データ集計や分析等を行なっています。
今回は、自分が開発業務や分析業務で日常的に利用しているIPython Notebookを便利に使う方法を紹介させて頂きます。
基本的に、pipでライブラリがインストールできる環境とnumpy、scipyの環境が揃っていることが前提で進めます。windows環境であれば、ライブラリのインストールは以下のURLのパッケージ群を利用すると快適です。とはいえこのページで紹介しているライブラリはwindowsでは動かないものもあります。
今回はIPython Notebookを既に活用している人向けのエントリーになっていますが、逐次更新してPandasとStatsmodelsとScikitsを用いたデータ処理や、統計解析に関する様々なユースケースをこのIPython Notebook上に追記してまとめていきたいと考えています。
IPython Notebookは、pythonのインタラクティブシェルであるipythonをブラウザ上から利用する為のインターフェイスです。ipythonを起動するときに、オプションにnotebookとつけることで、インタプリタがブラウザ上に起動します。(依存ライブラリなどはこちらを参照 > http://ipython.org/ipython-doc/dev/install/install.html)
ipython notebook
ipythonの機能も含みますが、以下の様な機能と特徴を持っています。
分散処理周りは自前で作っているceleryやpyro4の環境があるので基本使っていないのですが、アルゴリズムの確認や、アドホックな集計分析依頼に対応する際には、ブラウザ上でグラフや出力結果を確認しながらコードも快適に記述できるということは集計・解析の作業速度を大幅に早めることになります。
データの分析や、アルゴリズム実装のテストを行うにあたって、IPython NotebookとPythonの分析系のエコシステムは非常に協力なツールです。特に、ビジネスサイドからのアドホックな要求に対して、データの取得、分析処理、可視化をインタラクティブに実行でき、そのすべてのフェーズにおいてPythonの様々なライブラリを利用することができるのがRやMatlabとの大きな差分になります。
Gunosyでは、解析用のログ情報をAWS上のRedshiftに蓄積しています。定常的なKPIの取得や可視化等はバッチで集計した上で社内向けのWebアプリケーション上で閲覧することができるようにしています。仮説検証やアルゴリズムの構築検証を行う際に、IPython Notebookを活用しています。VPN経由か、SSH Portforwardingを用いて手元のIPython NotebookからRedshiftを叩いてデータの分析などを行っています。Redshiftから取得したデータをPandasのDataFrameにして処理を加えたあと、インタラクティブに可視化を行いつつ、データの分析やアルゴリズムの開発を日々行っています。
Pythonでデータの可視化を行うにあたって、標準的にmatplotlibというmatlabやRのplotライブラリのようなインターフェイスで利用可能な協力なplottingライブラリがあり、こちらもIPython Notebook上で気軽に利用することが出来ます。IPython Notebookで利用する場合は、matplotlibをインストールの後に、IPython Notebookを起動する際に、以下のようなオプションをつけることで、そのままブラウザ上で画像を表示することが出来ます。
matplotlibのインストールは、windowsの場合はインストーラで、macやlinuxの場合はportやyum等のパッケージ管理システムでインストールするのが楽ですが、pipでも入ります。
% ipython notebook --pylab inline
IPython Notebookを起動後にplotを呼び出すと、グラフが画像の形でブラウザ上に表示されます。matplotlibの使い方を知っていれば、通常のものと基本的には利用方法は変わらないため、特に問題なくグラフが描画できると思います。subplotで複数のグラフを表示することも可能です。
クロススペクトルを表示するサンプルをIPython Notebook上で実行するとそのまま、IPython Notebookのセルの下にグラフが描画されます。
import numpy as np
dt = 0.01
t = np.arange(0, 30, dt)
nse1 = np.random.randn(len(t)) # white noise 1
nse2 = np.random.randn(len(t)) # white noise 2
r = np.exp(-t/0.05)
cnse1 = np.convolve(nse1, r, mode='same')*dt # colored noise 1
cnse2 = np.convolve(nse2, r, mode='same')*dt # colored noise 2
# two signals with a coherent part and a random part
s1 = 0.01*np.sin(2*np.pi*10*t) + cnse1
s2 = 0.01*np.sin(2*np.pi*10*t) + cnse2
plt.subplot(211)
plt.xlim(0,5)
plt.plot(t, s1, 'b-', t, s2, 'g-')
plt.xlabel('time')
plt.ylabel('s1 and s2')
plt.grid(True)
plt.subplot(212)
cxy, f = plt.csd(s1, s2, 256, 1./dt)
plt.ylabel('CSD (db)')
plt.show()
IPython Notebookの中でmatplotlibのplotを呼ぶと、そのままpngの画像を生成されてHTMLの中に埋め込まれるのですが、元々の機能としてあった、拡大縮小や移動・回転などがインタラクティブにできなくなってしまいます。画像の生成時にサイズや位置などはいじれますが、インタラクティブにグラフを見ることが出来ないのは不便です。散布図を描画した時など、データ点が密集していると細かいところが見えなかったりするので、拡大縮小はデータ分析時にはとても欲しい機能です。
そこで、IPython Notebook上でインタラクティブにグラフを操作することができるJavascriptベースのプロッティングライブラリがBokehです。
pipでインストール出来ますが、Python自体をAnaconda Scientific Python Distributionに変えて、それのパッケージ管理システムであるcondaを利用すると他のライブラリの導入なども楽にできるようになります。
$ conda install bokeh or $ pip install bokeh
output_notebookという関数をインポートして実行することで、IPython Notebook上でbokehのグラフ描画が可能になります。
from bokeh.plotting import output_notebook
output_notebook()
Configuring embedded BokehJS mode.
グラフにlineプロットやscatterプロットをする前に、figure関数で描画面を作成します。この際、toolsに拡大縮小や保存等のオプションを付与することで、インタラクティブに操作が可能なグラフを描画できるようになります。描画面を作成後、hold関数を呼ぶことで、ひとつのfigureに複数のプロットを描画することが出来ます。
from bokeh.plotting import line, show, figure, xaxis, hold, yaxis
from bokeh.objects import Range1d
import numpy as np
dt = 0.01
t = np.arange(0, 30, dt)
nse1 = np.random.randn(len(t)) # white noise 1
nse2 = np.random.randn(len(t)) # white noise 2
r = np.exp(-t/0.05)
cnse1 = np.convolve(nse1, r, mode='same')*dt # colored noise 1
cnse2 = np.convolve(nse2, r, mode='same')*dt # colored noise 2
# two signals with a coherent part and a random part
s1 = 0.01*np.sin(2*np.pi*10*t) + cnse1
s2 = 0.01*np.sin(2*np.pi*10*t) + cnse2
hold()
figure(tools="pan,wheel_zoom,box_zoom,reset,previewsave", title="Signals", plot_height=400)
line(t[:500], s1[:500], color="blue", legend='s1', x_range = Range1d(start=0, end=5))
line(t[:500], s2[:500], color='green', legend='s2', x_range = Range1d(start=0, end=5))
yaxis().axis_label='s1 and s2'
yaxis().axis_label_text_font_size = "12pt"
hold(False)
show()
hold()
figure(tools="pan,wheel_zoom,box_zoom,reset,previewsave", title="CSD", plot_height=200)
csd = matplotlib.mlab.csd(s1, s2, 256, 1./dt)
line(csd[1], np.log10(np.abs(csd[0]))*10, x_range = Range1d(start=0, end=50))
yaxis().axis_label='CSD(db)'
xaxis().axis_label='Frequency'
xaxis().axis_label_text_font_size = "12pt"
yaxis().axis_label_text_font_size = "12pt"
show()
拡大縮小や移動が出来る上、グラフも非常にfancyです。ほとんどドキュメントが無いのがたまにキズですが、普通に使うだけならそこそこ綺麗に可視化されるので、ビジネスサイドの人間と議論しながら集計して可視化するとなかなかに受けが良いです。
データに対する集計を行う場合、numpyとpandasをうまく活用することでそこそこの速度で実行して、IPython Notebookで可視化出来ますが、それでも動的計画法などの多重Loopが必要になる場合などはPythonの実行速度の遅さが気になります。大抵の場合、C-Extensionを作ることになるのですが、ビルドが面倒だったり、setup.pyを作りこんだりと結構めんどいです。また、アルゴリズムの実装などであれば、インタプリタ上でC-Extensionを書きながらテストできると実装の効率も良くなります。 pythonのコードの中にC++のコードを埋め込むscipy.waeveというものも有りますが、結局PythonのC拡張を書くのと変わらないため、PythonのC言語レベルのAPIの知識が必要になります。PythonのJitコンパイラであるPyPyは、Pythonのコードをそのまま早くしてくれるという意味では便利なのですが、Pythonの数値計算系のライブラリで対応しているものも少なく、そもそもPyPy向けのnumpyも現在開発中で完全に利用できるわけではないので、アドホック分析に利用するというにはまだまだ環境が整っていないという現状が有ります。
そこで、IPython Notebookで気軽に動かせ、既存のPythonのデータ分析系の環境と親和性の高い高速化ソリューションとして、numbaとcythonを紹介します。
numbaはAnaconda Scientific Python Distributionという数値計算向けのPython Distributionを作っているContinuum Analyticsが作成している、Numpy AwareなPythonコードのJitコンパイラです。
llvmをバックエンドにしており、Pythonのコードをネイティブコード変換してくれたり、Cuda向けのGPGPUコードにコンパイルしてくれる優れ物です(とはえい、numbaでコンパイルできるようにコードを書き換えたり、既存のNumpyを利用して書いた複雑なコードをコンパイルしてもきちんとチューニングしないと性能が出なかったりしますが……)。
pypyと異なり、CPythonの環境の中でJitコンパイルを利用することができ、他のPythonライブラリや既存環境にそのまま組み込むことができます。以前は有料版のNumbaProのみがGPGPUに対応していたのですが、現在は部分的にではありますが、オープンソース版のNumbaでもGPGPU向けにコンパイルしてくれるようです。
インストールはpipでも出来ますが、Anaconda Pythonに付いてくるライブラリ管理システムcondaでインストールするのが楽です。Anaconda PythonはPython2.7ベースで普通にpipなどでライブラリもインストールできるので、いっその事皆さん手元のPython環境自体Anaconda Pythonに置き換えると良いでしょう。Anaconda Python自体は無料で、有料のパッケージであるAnaconda Accelerateを購入するとIntel MKLをリンクしてビルドされていたり、GPUに対応したNumpyが使えて手軽に高速化が出来ますし、科学計算用ライブラリも結構整備されてます(ステマ)。もちろんLLMM を利用しているので、LLVMが動く環境で無いと動作しません。
インストールはコマンドラインで一発です。そのままllvmpyも入れてくれます
$ conda install numba or $ pip install numba
公式ドキュメントにサンプルのマンデルブロー集合を描画するコードをIPython Notebookで動かしてみましょう。
基本的に、numbaを使う場合は、numba.jitというデコレータを使ってPythonの関数をJitコンパイルします。
import numpy as np
def mandel(x, y, max_iters):
"""
Given the real and imaginary parts of a complex number,
determine if it is a candidate for membership in the Mandelbrot
set given a fixed number of iterations.
"""
i = 0
c = complex(x,y)
z = 0.0j
for i in range(max_iters):
z = z*z + c
if (z.real*z.real + z.imag*z.imag) >= 4:
return i
return 255
def create_fractal(min_x, max_x, min_y, max_y, image, iters):
height = image.shape[0]
width = image.shape[1]
pixel_size_x = (max_x - min_x) / width
pixel_size_y = (max_y - min_y) / height
for x in range(width):
real = min_x + x * pixel_size_x
for y in range(height):
imag = min_y + y * pixel_size_y
color = mandel(real, imag, iters)
image[y, x] = color
return image
とりあえずJitコンパイルを利用しない生のPythonだと、1ループ5秒程度。
image = np.zeros((1000, 1300), dtype=np.uint8)
%timeit create_fractal(-2.0, 1.0, -1.0, 1.0, image, 20)
ret = create_fractal(-2.0, 1.0, -1.0, 1.0, image, 20)
imshow(ret)
1 loops, best of 3: 5.57 s per loop
<matplotlib.image.AxesImage at 0x10aa4f910>
関数定義に対してjitデコレータをつけて実行してみましょう。とりあえずデコレータをつけるのではなく、関数定義に対してnumba.jitを適用してみます。普通にやるなら、関数定義分の上に@numba.jitとつけるだけです。
import numba
mandel = numba.jit(mandel)
create_fractal = numba.jit(create_fractal)
先ほどと同様に実行してみると、一周50ms程度になっています。爆速ですね。このように多重ループが大量に有るようなコードであれば比較的簡単に高速化出来ます。
image = np.zeros((1000, 1300), dtype=np.uint8)
%timeit create_fractal(-2.0, 1.0, -1.0, 1.0, image, 20)
ret = create_fractal(-2.0, 1.0, -1.0, 1.0, image, 20)
1 loops, best of 3: 48.9 ms per loop
jitコンパイル時にエラーが出たり、うまく実行されない場合は、jitコンパイル済みの関数に対して、inspect_types関数を呼ぶことで、numbaがどのように解釈したかを表示することが出来ます。その際には、関数の引数と返り値の型をjitデコレータにわたして上げる必要があります。
型の表記方法については以下を参照してください。
@numba.jit("f8[:,:](i8,i8,i8,i8,f8[:,:],i8)")
def create_fractal2(min_x, max_x, min_y, max_y, image, iters):
height = image.shape[0]
width = image.shape[1]
pixel_size_x = (max_x - min_x) / width
pixel_size_y = (max_y - min_y) / height
for x in range(width):
real = min_x + x * pixel_size_x
for y in range(height):
imag = min_y + y * pixel_size_y
color = mandel(real, imag, iters)
image[y, x] = color
return image
create_fractal2.inspect_types()
create_fractal2 (int64, int64, int64, int64, array(float64, 2d, A), int64) -------------------------------------------------------------------------------- # File: <ipython-input-49-ddaab3774099> # --- LINE 1 --- @numba.jit("f8[:,:](i8,i8,i8,i8,f8[:,:],i8)") # --- LINE 2 --- def create_fractal2(min_x, max_x, min_y, max_y, image, iters): # --- LINE 3 --- # label 0 # min_x.1 = min_x :: int64 # max_x.1 = max_x :: int64 # min_y.1 = min_y :: int64 # max_y.1 = max_y :: int64 # image.1 = image :: array(float64, 2d, A) # iters.1 = iters :: int64 # $0.1 = getattr(attr=shape, value=image.1) :: (int64 x 2) # $0.2 = const(<type 'int'>, 0) :: int32 # $0.3 = getitem(index=$0.2, target=$0.1) :: int64 # height = $0.3 :: int64 height = image.shape[0] # --- LINE 4 --- # $0.4 = getattr(attr=shape, value=image.1) :: (int64 x 2) # $0.5 = const(<type 'int'>, 1) :: int32 # $0.6 = getitem(index=$0.5, target=$0.4) :: int64 # width = $0.6 :: int64 width = image.shape[1] # --- LINE 5 --- # --- LINE 6 --- # $0.7 = max_x.1 - min_x.1 :: int64 # $0.8 = $0.7 /? width :: int64 # pixel_size_x = $0.8 :: int64 pixel_size_x = (max_x - min_x) / width # --- LINE 7 --- # $0.9 = max_y.1 - min_y.1 :: int64 # $0.10 = $0.9 /? height :: int64 # pixel_size_y = $0.10 :: int64 pixel_size_y = (max_y - min_y) / height # --- LINE 8 --- # jump 54 # label 67 # $67.1 = iternext(value=$54.3) :: int64 # $67.2 = itervalid(value=$54.3) :: bool # branch $67.2, 70, 161 # label 70 # $70.1 = $67.1 :: int64 # x = $70.1 :: int64 # label 54 # $54.1 = global(range: <built-in function range>) :: range # $54.2 = call $54.1(width, ) :: (int64,) -> range_state64 # $54.3 = getiter(value=$54.2) :: range_iter64 # jump 67 for x in range(width): # --- LINE 9 --- # $70.2 = x * pixel_size_x :: int64 # $70.3 = min_x.1 + $70.2 :: int64 # real = $70.3 :: int64 real = min_x + x * pixel_size_x # --- LINE 10 --- # label 100 # $100.1 = iternext(value=$87.3) :: int64 # $100.2 = itervalid(value=$87.3) :: bool # branch $100.2, 103, 157 # jump 87 # label 103 # $103.1 = $100.1 :: int64 # y = $103.1 :: int64 # label 87 # $87.1 = global(range: <built-in function range>) :: range # $87.2 = call $87.1(height, ) :: (int64,) -> range_state64 # $87.3 = getiter(value=$87.2) :: range_iter64 # jump 100 for y in range(height): # --- LINE 11 --- # $103.2 = y * pixel_size_y :: int64 # $103.3 = min_y.1 + $103.2 :: int64 # imag = $103.3 :: int64 imag = min_y + y * pixel_size_y # --- LINE 12 --- # $103.4 = global(mandel: CPUOverloaded(<function mandel at 0x10ae4a938>)) :: Dispatcher(CPUOverloaded(<function mandel at 0x10ae4a938>)) # $103.5 = call $103.4(real, imag, iters.1, ) :: (int64, int64, int64) -> int64 # color = $103.5 :: int64 color = mandel(real, imag, iters) # --- LINE 13 --- # label 161 # del $54.3 # $103.6 = build_tuple(items=[Var(y, <ipython-input-49-ddaab3774099> (10)), Var(x, <ipython-input-49-ddaab3774099> (8))]) :: (int64 x 2) # image.1[$103.6] = color :: (array(float64, 2d, A), (int64 x 2), float64) -> none # jump 100 # label 157 # del $87.3 # jump 158 # label 158 # jump 67 image[y, x] = color # --- LINE 14 --- # --- LINE 15 --- # jump 162 # label 162 # return image.1 return image ================================================================================
numbaはこのようにPythonのコードをかなりお手軽に高速化してくれますが、既存のnumpyの関数やscipyの関数を用いたコードをそのまま早くできるかというと実際はそうではありません。numpyの関数によっては、numba.jitの中で利用することが出来ず、高速化どころか実行すら出来ないようになってしまうものがあります(sumなどの返り値の型やshapeが引数によって様々に変わってしまうもの等)。numpy、scipyを使い倒したような既存の実装をnumbaで高速化するには大幅な書き換えが必要になります。とはいえ、numpy、scipyをそれほど使用していないコードに関してはPythonで書いたコードをかなり高速化してくれるので、集計時のボトルネックつぶしなどには非常に有用だと言えるでしょう。
numba以外で処理を高速化をするためのライブラリとしては、Cythonがあります。
Cythonは、numbaのようにPythonのコードをjitコンパイルしてくれるものではなく、Pythonに型情報を付けたようなDSLで書いたコードをC or C++にコンパイルしてpythonから利用できる形のsoファイルを作成してくれます。numbaと違ってJITコンパイラではなく、Cのコードに変換してC拡張を作成するものであるため、Cython言語という型付のPython Likeなコードを書く必要があります。PythonのコードをそのままCythonでコンパイルしても、型宣言をしないとそのままPythonのオブジェクトをC上で動かすことになるため高速化されません。Cへの変換器であるので、Cのコンパイラも必要です。
以前は自前でCythonからCへの変換コマンドを叩いたり、setup.pyを書いたりといろいろ面倒でしたが、最近はpyximportという同梱されているモジュールを利用すると、Cythonのコード(*.pyx)をimport時に自動でコンパイルして読み込んでくれたりと便利になっています。さらに、IPython Notebook上でcythonを利用する場合は、便利なmagic functionが用意されており、IPython Notebook上でコードを書いてコンパイルして読み込むことが出来ます。
magic functionをloadするには、まず以下のコマンドを実行します
%load_ext cythonmagic
まず、pythonの式をそのままcythonでコンパイルして実行するにはcython_inlineコマンドを利用します。ほとんど使うことは無いと思いますが、dynamicに式をコンパイルして実行することができます。
%%cython_inline
cdef int a = 10
cdef int b = 5
return a*b
Compiling /Users/aihara/Library/Caches/cython/inline/_cython_inline_c29d81f5f7b0e71e9ade91a50f5bb972.pyx because it changed. Cythonizing /Users/aihara/Library/Caches/cython/inline/_cython_inline_c29d81f5f7b0e71e9ade91a50f5bb972.pyx
50
IPython Notebook上で、cythonで書いた関数をコンパイルして利用する場合は、cython_pyximportコマンドを利用します。比較のためにnumbaのsampleをCythonに書き直してみましょう。numbaの時のように実行時に型情報を認識してくれないので、型宣言を多く追加してあります。
%%cython_pyximport cython_mandelbrot
cdef unsigned char mandel(double x, double y, int max_iters):
"""
Given the real and imaginary parts of a complex number,
determine if it is a candidate for membership in the Mandelbrot
set given a fixed number of iterations.
"""
cdef unsigned char i = 0
cdef double complex c = complex(x,y)
cdef double complex z = 0.0j
cdef double zreal, zimag
for i in range(max_iters):
z = z*z + c
if (z.real*z.real + z.imag*z.imag) >= 4:
return i
return 255
def create_fractal(double min_x, double max_x, double min_y, double max_y, unsigned char[:, ::1] image, int iters):
cdef int height = image.shape[0]
cdef int width = image.shape[1]
cdef double pixel_size_x = (max_x - min_x) / width
cdef double pixel_size_y = (max_y - min_y) / height
cdef int x, y,
cdef unsigned char color
cdef double real, imag
for x in range(width):
real = min_x + x * pixel_size_x
for y in range(height):
imag = min_y + y * pixel_size_y
color = mandel(real, imag, iters)
image[y, x] = color
return image
cython_pyximportを記載したIPython Notebookのcellを実行するとコンパイルされます。それでは実行してみましょう。cython_pyximportの引数に与えた名前がモジュール名になるので、そのモジュールからオブジェクトをインポートして実行します。
from cython_mandelbrot import create_fractal
image = np.zeros((1000, 1300), dtype=np.uint8)
%timeit create_fractal(-2.0, 1.0, -1.0, 1.0, image, 20)
ret = create_fractal(-2.0, 1.0, -1.0, 1.0, image, 20)
imshow(ret)
1 loops, best of 3: 249 ms per loop
<matplotlib.image.AxesImage at 0x10b3cbdd0>
型宣言をいっぱい追加した割には、245msとnumbaのJITコンパイルを利用したものより5倍程度遅い結果になってしまいました。
これだけ見るとCythonは型宣言をする必要があったりと面倒な上遅いしいいところがないようにも見えますが、numbaはまだ発展段階であり、複雑なコードだと高速化どころか遅くなってしまうことも多いです。numba向けのコードを最初から書くというのも有りですが、Cythonはチューニングポイントに関するノウハウも多く、 細かいチューニングポイントさえ抑えてCython向けのコードを書くことでかなり複雑なアルゴリズム高速化することが出来ます。また、既存のC/C++のコードをラップしたりする場合にはとても有用です。
Cythonの使い方を勉強するなら、以下のライブラリのソースコードを見てみると良いかと思います。
CythonやC-Extensionとnumpyを組み合わせた画像処理、音声信号処理、自然言語処理向けのライブラリを複数公開しています。勉強中に作ったようなコードも有りますが、Cythonやnumpyの活用のサンプルコードとして活用していただければと思います。