ctypesによるarrayの受け渡し

計算処理時間がpythonを使っていると問題になることがあります。
pythonはスクリプト言語なので仕方ないといえばそれまでですが、一部分だけでもC言語で書いて実行形式で計算したいと思うことがあります。
ただ、そう思うのが、開発を始めた頃ではなく、開発がほとんど終盤に差し掛かった頃で、移行が簡単ではないことが多いです。

自作class やscipyやnumpyを使い、引数をやりとりしていると、簡単には置き換えしにくくなります。
今回はnumpyのctypesを使ってnumpy.arrayを引数として渡す方法についてまとめてみます。

手順は

1. c言語でダイナミックリンクライブラリを作ります
2. pythonプログラムでライブラリを読み込む設定を書き込みます
1. load_libraryで読み込み
2. 引数の型設定
3. 戻り値の型設定
4. テスト

という感じです。


１
C言語によるダイナミックリンクライブラリ作成
arrayで渡したい部分は配列にします。
たとえば
double[] hoge( double x[], int n ){
...
}

と言う感じです。
ファイル名はここでは、libhoge.soと言う名前でつくります。
ウィンドウズの場合は拡張子dllにする必要があります。
gccの場合、
　gcc -shared hoge.c -o libhoge.so
とします。


2.1 python側
ダイナミックリンク読み込み
dll = numpy.ctypeslib.load_library( 'libhoge.so', '.' )
通常ctypesというライブラリが標準でついているのそちらを使っていますが、numpyのctypeslibはarrayを使うときに便利です。
numpyの場合、実行ディレクトリを指定する必要があります。


2.2. 引数型の設定
ndtype = numpy.ctypeslib.ndpointer( dtype=numpy.float64, ndim=1, flags='C_CONTIGUOUS' )
dll.hoge.argtypes = [ ndtype, numpy.ctypeslib.ct.c_int,  ]


2.3.戻り値型指定
dll.hoge.restype = numpy.ctypeslib.ndpointer( dtype=numpy.float64, shape=(100,), flags='C_CONTIGUOUS' )

2.4　実行テスト

x = dll.hoge( numpy.array([ 0.1, 0.3, 0.5],dtype=float), 3 )

として実行してみます。

これでxがnumpy.arrayになっていればOKです。

2.2や2.3の型の設定がややこしいところです。

確率統計：モンテカルロ

モンテカルロ計算を考えてみます。

統計を含め、科学技術計算には、numpyが使えます。
行列操作ができ、乱数や確率に関するライブラリが豊富です。
numpyは数値計算用のライブラリであり、
中身の逆行列計算などのアルゴリズムはCなどで書かれています。
そのため、うまく書けばpythonとは言えども高速な計算が可能です。

モンテカルロと言うのはモナコにある町の名前です。
計算機科学の世界でいういわゆるモンテカルロ法というのは確率的なサンプリングを行うことによって、計算を行うことを指し、
フォンノイマンによって命名されたと言われています。

今回はnumpyで簡単なモンテカルロ計算を行います。

numpyをまずimportします。
import numpy

中心０分散１の正規乱数を発生させるには、

x = numpy.random.randn(10000)

とします。

同一正規分布に従う二つの独立な変数x1,x2があったとします。
常にx1<x2であるという制約があった場合、この二つはどのような分布となるでしょうか。

x1 = numpy.random.randn(10000)
x2 = numpy.random.randn(10000)
index = x1 < x2

x1とx2をプロットすると

この場合、切断正規分布という分布になります。
こうするとx1とx2は独立ではなくなり、p(x1,x2)と書くと、x1の確率は

p(x1) = ∫p(x1,x2) dX2

となります。
x1とx2のヒストグラムをとるとこんな感じです。

このようにサンプリングを通して、面積や分布などを調べることをモンテカルロ法といいます。
よく円の面積などの例が教科書などに載っていると思います。
モンテカルロは概念が簡単で容易に理解できますが、奥が深いです。


プロットはmatplotlibを使いました。
matplotlibについてはまたそのうちレポートします。

擬似unixコマンド１

特に仕事でウィンドウズ使っているとunixコマンドを使いたくなることがあります。

そんなの時のために擬似的に似たような動作をする簡易的なpythonプログラムを作っています。
例えばgrep
第1引数の文字列が含まれている行のみ表示するコマンドです。

grep x file

こんな風に使います。
つまり、この場合fileの中身でxが入っている行のみを表示します。
といってもgrepにはもっといろんな機能が付いていますが、簡易版なので、この機能のみを実装してみます。
pythonで書くとかなり簡単に書けます。


import sys

for line in sys.stdin:
    if sys.argv[1] in line:
        print line.rstrip()


処理が遅いと思うのならばリスト内包表記で

import sys
print '\n'.join( [ line.rstrip()   for line in sys.stdin if sys.argv[1] in line ] )

とか書くと少し処理が速くなります。
これをpgrep.pyとかいう名前で保存し、

pgrep.py import pgrep.py
などとコマンドラインで実行すると
文字列importのある行のみ抜き出して表示されます。
いまは作成したプログラムpgrep.pyに対して処理してみましたが、3番目のファイルはテキストならなんでも構いません。
書式は下のとおりです。

pgrep.py [検索文字列] [対象ファイル名]

はじめに

数年前からpythonを使って仕事しています。
科学技術計算のためのソフトウェアとして最近注目されてきているように思います。
ロボット分野だとそれなりにユーザがいて、
ライブラリも充実しています。
ロボット用だけではなくWEB関連もあり、画像処理や統計計算など種類も豊富です。

しかし、新しいことをやろうとすると日本語ドキュメントがないことも多いです。
なので、備忘録もかねて、pythonの使い方についてレポートします。
内容は専門的なことを書きます。

文法構造等の基本的なことは、もっとわかりやすく書ける人がいるので、他のサイトを参考にしてください。



pythonは
使っているライブラリの対応の関係から、version 2.6を使っています。