前回
PythonにはNumpy（とScipy）というMatlabに似たインターフェースを持つ数値計算に非常に便利なパッケージがあるのはご存知でしょう。というか、Weaveについて調べてココにたどり着いたなら、知っていて当然ですね。例えばこんな感じ。
import numpy as np # numpyのインポート
# 乱数配列の生成 2行3列
a = np.random.rand(2,3)
b = np.random.rand(2,3)
# 各要素ごとの掛け算 Matlabだと a .* b
c = a * b
# 行列としての掛け算 Matlabだと a' * b
d = np.dot(a.T, b) # C = A^T * B （Cは2行2列になります）
# 固有値分解
[E, V] = np.linalg.eig(d)
Matlabでできることはほぼ全部できます。しかも、無料で。このページは参考になります。
Numpy for Matlab Users
行列演算や要素同士、固有値計算などのメソッドなどはCで実装されているので高速に実行できます。しかし、こういった汎用的な演算でない場合はPython上でループを回したりして配列にランダムアクセスするようなコードを書く必要がでてきます。
Pythonのようなスクリプト言語でそういったコードを書くと激遅になりますよね。例えば、ある配列の行のインデックスと列のインデックスを足した値が3の倍数だったら0に書き換えるというコードをPythonで書くと次のようになります（あまり良い例ではないですが、お許しください！）。
def numpy_mult3_0_pure(npArr):
for r in xrange(npArr.shape[0]):
for c in xrange(npArr.shape[1]):
if (r + c) % 3 == 0:
npArr[r, c] = 0
これを次のように呼び出す。
>>> side = 1000
>>> npArr = np.ones((side, side))
>>> numpy_mult3_0_pure(npArr)
>>> print npArr
結果は、以下。
[[ 0. 1. 1. ..., 1. 1. 0.]
[ 1. 1. 0. ..., 1. 0. 1.]
[ 1. 0. 1. ..., 0. 1. 1.]
...,
[ 1. 1. 0. ..., 1. 0. 1.]
[ 1. 0. 1. ..., 0. 1. 1.]
[ 0. 1. 1. ..., 1. 1. 0.]]
時間を計測すると、0.384463787079[秒]でした。
これをWeaveを使ってC化して高速化します。この例に移る前にNumpyをWeaveに使う場合の基本から見ていきます。
weave.inlineはPython変数を渡すと自動的にC内で適当な型に変換してくれます（converterを指定することも可能）。Numpyのndarrayを渡した場合も自動で色々と変換を施してくれます。次の例を見てみます。
def numpy_converted_vars():
npArr = np.zeros((2, 3, 4), dtype=np.uint16)
npArr[1, 1, 3] = 1
print '(py) ndim:', npArr.ndim
print '(py) shape:', npArr.shape
print '(py) strides:', npArr.strides
code = """
// Number of dimension
std::cout << "ndim: " << DnpArr << std::endl;
// Shape
std::cout << "shape:";
for (int i = 0; i < DnpArr; i++) std::cout << " " << NnpArr[i];
std::cout << std::endl;
// Strides
std::cout << "strides:";
for (int i = 0; i < DnpArr; i++) std::cout << " " << SnpArr[i];
std::cout << std::endl;
// Access to data as a three dimensional array
std::cout << "npArr[1, 1, 3]: " << NPARR3(1, 1, 3) << std::endl;
"""
weave.inline(code, [u'npArr'])
このメソッド内では、まず2×3×4の形のすべての要素が0の３次元配列を作ります。確認の為に[1,1,3]の要素だけ1を代入しておきます。
データ型がuint16=2[bytes]で、3次元、C-contiguousのメモリ上のデータ配置で、2 * 3 * 4の形なので、上の3つのprintの結果は
(py) ndim: 3
(py) shape: (2, 3, 4)
(py) strides: (24, 8, 2)
となります。stridesのところとかよくわからない人はこの画像を見て理解するかググってください。
Cのinlineコードの中でこれらにアクセスする方法を見ていきます。ndim, shape, stridesはinlineメソッドの内部でN, D, Sという接頭辞をつけてNnpArr（int）, DnpArr（npy_int*）, SnpArr（npy_int*）という名前で使用可能になります。npy_int型はint64かint32（環境に依存）のtypedefだと思います。実際にCの内部でどのように変換されているかを生成された.cファイルを見てみます。npArrのオブジェクトはCコードの中では、PyArrayObjectというNumpyのC内部での構造体として解釈し、npArr_arrayという名前で使えるようにしてくれています。
PyArrayObject* npArr_array = convert_to_numpy(py_npArr,"npArr");
NumpyのC APIについては各自ぐぐって欲しいのですが、PyArrayObject構造体のメンバにstridesなどがあるので、
npy_intp* NnpArr = npArr_array->dimensions;
npy_intp* SnpArr = npArr_array->strides;
int DnpArr = npArr_array->nd;
で値を取り出してくれています。inlineコードのndim, shape, stridesの出力は次のようになり、python側での出力と一致しています。
ndim: 3
shape: 2 3 4
strides: 24 8 2
肝心のnpArrという名前は以下のように配列のデータへの先頭ポインタのアドレスを取り出したものとして宣言しています。
double* npArr = (double*) npArr_array->data;
一次元配列としてアクセスする場合はnpArr[0]などとアクセスすれば値を取り出せますが、この例では3次元配列なので、stridesを使って要素にアクセス必要があります。たとえばnpArr[i, j, k]にアクセスしたい場合は、以下のようにする必要があります（PyArrayObject*->dataはchar*型）。
*((double*)(npArr_array->data + (i * SnpArr[0] + j * SnpArr[1] + k * SnpArr[2])))
しかし、これはめんどくさい。実はこれをやるためのマクロ関数を最高4次元配列の場合まで定義してくれています。
#define NPARR1(i) (*((double*)(npArr_array->data + (i)*SnpArr[0])))
#define NPARR2(i,j) (*((double*)(npArr_array->data + (i)*SnpArr[0] + (j)*SnpArr[1])))
#define NPARR3(i,j,k) (*((double*)(npArr_array->data + (i)*SnpArr[0] + (j)*SnpArr[1] + (k)*SnpArr[2])))
#define NPARR4(i,j,k,l) (*((double*)(npArr_array->data + (i)*SnpArr[0] + (j)*SnpArr[1] + (k)*SnpArr[2] + (l)*SnpArr[3])
3次元配列の場合はnpArrの名前を大文字にして後ろに3をつけた名前が付けられていますので、上のようにpython側でnpArr[1,1,3]=1で変更した値にNPARR3(1,1,3)でアクセスできています。
npArr[1, 1, 3]: 1
今回はNumpyのweave.inlineへの受け渡しの基本を紹介しました。
次回に実際計算してみます。
