NumPy を扱う

普通の Python と全く同様に、Cython からも NumPy を操作できます。ただし、 そのままでは、相当な高速化の可能性を捨てていることになります。というの も、Cython には NumPy のアレイに高速にアクセスする機能があるからです。 その仕組みを、簡単な例でみていきましょう。

下のコードは、ある画像とフィルタを、二次元の離散畳込みしています (もち ろん、もっといい方法があります!デモのためということでご容赦ください) このコードは、 Python のコードとしても、 Cython のコードとしても有効で す。 Python 版は convolve_py.py, Cython 版は convolve1.pyx というファイルで参照することにします – Cython のファイル拡張子は ”.pyx” です。

from __future__ import division
import numpy as np
def naive_convolve(f, g):
    # f は画像で、 (v, w) でインデクスする
    # g はフィルタカーネルで、 (s, t) でインデクスする
    #   ディメンジョンは奇数でなくてはならない
    # h は出力画像で、 (x, y) でインデクスする
    #   クロップしない
    if g.shape[0] % 2 != 1 or g.shape[1] % 2 != 1:
        raise ValueError("Only odd dimensions on filter supported")
    # smid と tmid は中心ピクセルからエッジまでのピクセル数、
    # つまり 5x5 のフィルタならどちらも 2
    #
    # 出力画像のサイズは入力画像の両端に smid と tmid を足した
    # 値となる
    vmax = f.shape[0]
    wmax = f.shape[1]
    smax = g.shape[0]
    tmax = g.shape[1]
    smid = smax // 2
    tmid = tmax // 2
    xmax = vmax + 2*smid
    ymax = wmax + 2*tmid
    # 出力画像の確保
    h = np.zeros([xmax, ymax], dtype=f.dtype)
    # コンボリューションの演算
    for x in range(xmax):
        for y in range(ymax):
            # (x,y) におけるピクセル値 h を計算。
            # フィルタ g の各ピクセル (s, t) に対するコンポーネン
            # トを加算
            s_from = max(smid - x, -smid)
            s_to = min((xmax - x) - smid, smid + 1)
            t_from = max(tmid - y, -tmid)
            t_to = min((ymax - y) - tmid, tmid + 1)
            value = 0
            for s in range(s_from, s_to):
                for t in range(t_from, t_to):
                    v = x - smid + s
                    w = y - tmid + t
                    value += g[smid - s, tmid - t] * f[v, w]
            h[x, y] = value
    return h

このコードを Cython でコンパイルすると、 (Linux システムでは) yourmod.so が生成されます。Python セッションを起動して、 Python 版 (.py ファイルから import) と、コンパイルした Cython モ ジュールの両方をテストしてみます。

In [1]: import numpy as np
In [2]: import convolve_py
In [3]: convolve_py.naive_convolve(np.array([[1, 1, 1]], dtype=np.int),
...     np.array([[1],[2],[1]], dtype=np.int))
Out [3]:
array([[1, 1, 1],
    [2, 2, 2],
    [1, 1, 1]])
In [4]: import convolve1
In [4]: convolve1.naive_convolve(np.array([[1, 1, 1]], dtype=np.int),
...     np.array([[1],[2],[1]], dtype=np.int))
Out [4]:
array([[1, 1, 1],
    [2, 2, 2],
    [1, 1, 1]])
In [11]: N = 100
In [12]: f = np.arange(N*N, dtype=np.int).reshape((N,N))
In [13]: g = np.arange(81, dtype=np.int).reshape((9, 9))
In [19]: %timeit -n2 -r3 convolve_py.naive_convolve(f, g)
2 loops, best of 3: 1.86 s per loop
In [20]: %timeit -n2 -r3 convolve1.naive_convolve(f, g)
2 loops, best of 3: 1.41 s per loop

ご覧のとおり、大した違いはありません; というのも、 C のコードが実行し ている内容が、(例えば、数値を使うたびにオブジェクトを生成するといった 点で) 実際にPython インタプリタがやっているのと全く同じ処理だからです。 -a オプションで Cython が生成する html ファイルを見て何をしなけれ ばならないかを考えれば、最も単純な実行文でさえ問題をすぐに見つけられる でしょう。 Cython にはもっと沢山の情報を教えてやらねばなりません; まずは、型を追加するところからです。

型を追加する

型を追加するには、 Cython の構文を使わねばなりません。そのため、この時 点で Python とのソース互換性は失われてしまいます。以下のようなコードを 考えましょう (コメントを読んでくださいね!)

from __future__ import division
import numpy as np
# "cimport" は、 numpy モジュールを使うコードのコンパイル時に必要
# な情報を import するのに使います。
# (この情報は numpy.pxd に入っています。現状、 numpy.pxd は Cython
# ディストリビューションに入っています)
cimport numpy as np
# 変数 DTYPE を使って、アレイのデータ型に手を加えます。DTYPE には、
# NumPy のランタイム型情報オブジェクト (runtime type info object)
# を代入します。
DTYPE = np.int
# "ctypedef" は、コンパイル時に決定した DTYPE_t で参照できるよ
# うにします。numpy モジュールのデータ型には、すべて _t という
# サフィックスのついたコンパイル時用の型があります。
ctypedef np.int_t DTYPE_t
# "def" は引数の型指定を行えますが、戻り値の型を指定できません。
# "def" 関数型の引数の型は、関数に入るときに動的にチェックされます。
#
# アレイ f, g, h は "np.ndarray" インスタンス型に型付けされていま
# す。その効果は a) 関数の引数が本当に NumPy アレイかどうかチェッ
# クが入ることと、 b) f.shape[0] のようなアトリビュートアクセスの
# 一部が格段に効率的になること、だけです (この例の中では、どちらも
# さして意味はありません)。
def naive_convolve(np.ndarray f, np.ndarray g):
    if g.shape[0] % 2 != 1 or g.shape[1] % 2 != 1:
        raise ValueError("Only odd dimensions on filter supported")
    assert f.dtype == DTYPE and g.dtype == DTYPE
    # "cdef" キーワードは、関数の中で変数の型を宣言するのにも使い
    # ます。 "cdef" は、関数内のトップインデントレベルでしか使えま
    # せん (他の場所で使えるようにするのは大したことではありません。
    # もしいい考えがあったら提案してください)。
    #
    # インデクスには、 "int" 型を使います。この型は C の int 型に
    # 対応しています。 ("unsigned int" のような) 他の型も使えます。
    # アレイのインデクスとして適切な型を使いたい純粋主義の人は、
    # "Py_ssize_t" を使ってもかまいません。
    cdef int vmax = f.shape[0]
    cdef int wmax = f.shape[1]
    cdef int smax = g.shape[0]
    cdef int tmax = g.shape[1]
    cdef int smid = smax // 2
    cdef int tmid = tmax // 2
    cdef int xmax = vmax + 2*smid
    cdef int ymax = wmax + 2*tmid
    cdef np.ndarray h = np.zeros([xmax, ymax], dtype=DTYPE)
    cdef int x, y, s, t, v, w
    # 変数全てについて型を定義するのがとても大事です。型を定義し忘
    # れても何の警告もでませんが、(変数が暗黙のうちに Python オブ
    # ジェクトに片付けされるので) コードは極端に遅くなります。
    cdef int s_from, s_to, t_from, t_to
    # 変数 value に対しては、アレイに保存されいているのと同じデー
    # タ型を使いたいので、上で定義した "DTYPE_t" を使います。
    # 注! この操作には、重大な副作用があります。 "value" がデータ
    # 型の定義域をオーバフローすると、 Python のように例外が送出さ
    # れるのではなく、 C のときと同様に単なる桁落ち (wrap around)
    # を起こします。
    cdef DTYPE_t value
    for x in range(xmax):
        for y in range(ymax):
            s_from = max(smid - x, -smid)
            s_to = min((xmax - x) - smid, smid + 1)
            t_from = max(tmid - y, -tmid)
            t_to = min((ymax - y) - tmid, tmid + 1)
            value = 0
            for s in range(s_from, s_to):
                for t in range(t_from, t_to):
                    v = x - smid + s
                    w = y - tmid + t
                    value += g[smid - s, tmid - t] * f[v, w]
            h[x, y] = value
    return h

このモジュールをビルドして、 (とてもざっくりした) ベンチマークを実行す ると、手元ではこんな結果になりました:

In [21]: import convolve2
In [22]: %timeit -n2 -r3 convolve2.naive_convolve(f, g)
2 loops, best of 3: 828 ms per loop

インデクシングの効率化

プログラムには、まだパフォーマンスを低下させるボトルネックがあります。 それは、アレイの要素に対するルックアップと代入です。 [] オペレータ は、まだ完全に Python 的な操作を使っています – これを止めるには、デー タバッファに C のスピードで直接アクセスします。

そこで、次に必要になってくるのが ndarray オブジェクトの中身に対 する型付けです。型付けには、特殊な “buffer” 構文を使います。この構文は データ型 (最初の引数) と、ディメンジョンの数 (“ndim” キーワード引数。 省略すると、一次元とみなされます) を指定せねばなりません。

以下の変更が必要です:

...
def naive_convolve(np.ndarray[DTYPE_t, ndim=2] f, np.ndarray[DTYPE_t, ndim=2] g):
...
cdef np.ndarray[DTYPE_t, ndim=2] h = ...

以下のように使います:

In [18]: import convolve3
In [19]: %timeit -n3 -r100 convolve3.naive_convolve(f, g)
3 loops, best of 100: 11.6 ms per loop

この変更の重要性が理解できたでしょう。

落とし穴: インデクシングを高速化できるのは、特定のインデクス操作だけ、 中でも、 ndim に指定した数ぶん、整数に型づけしたインデクスを伴う操 作だけです。従って、例えば v を型付けしていなければ、 f[v, w] のルックアップは最適化されません。ただ、見方を変えれば、依然として、 あたかも型付けされていないアレイを操作するかのように、Python オブジェ クトを使って、洗練された書き方で動的なスライスなどを実行できるともいえ ます。

インデクシングをもっとチューニングする

まだ、以下に述べる二つの要素が、アレイのルックアップを遅くしています:

  1. 境界チェック (bounds checking) を行なっている。
  2. 負のインデクスをチェックしていて、それに対する処理を適切に行なって いる。

上のコードでは、負のインデクスを使わないように明にコードを書いているの で、常に (おそらく) 境界内の値にしかアクセスしていないはずです。 境界チェックは、デコレータを追加して無効にできます:

...
cimport cython
@cython.boundscheck(False) # 関数全体で境界チェックを無効化
def naive_convolve(np.ndarray[DTYPE_t, ndim=2] f, np.ndarray[DTYPE_t, ndim=2] g):
...

これで、境界チェックを行わなくなります (そして、副作用として、もし境界 の外のデータに「アクセスしてしまった」ら、少なくともプログラムがクラッ シュ、最悪の場合データが破損することでしょう)。境界チェックのスイッチ にはいくつか方法があります。詳しくは [:docs/compilerdirectives: コンパイラディレクティブ] を参照してくださ い。

負のインデクスは、 Cython に変数の値を unsigned integer にキャストさせ、 インデクスを正の値にして処理します (変数の値が負だと、キャストによって 巨大な正の値ができるので、やはり境界外へのアクセスが起き得ます)。 値のキャストは <> 構文で行います。以下のコードでは、適宜変数を unsigned int にして、キャストを施しています:

...
cdef int s, t                                                                            # changed
cdef unsigned int x, y, v, w                                                             # changed
cdef int s_from, s_to, t_from, t_to
cdef DTYPE_t value
for x in range(xmax):
    for y in range(ymax):
        s_from = max(smid - x, -smid)
        s_to = min((xmax - x) - smid, smid + 1)
        t_from = max(tmid - y, -tmid)
        t_to = min((ymax - y) - tmid, tmid + 1)
        value = 0
        for s in range(s_from, s_to):
            for t in range(t_from, t_to):
                v = <unsigned int>(x - smid + s)                                         # changed
                w = <unsigned int>(y - tmid + t)                                         # changed
                value += g[<unsigned int>(smid - s), <unsigned int>(tmid - t)] * f[v, w] # changed
        h[x, y] = value
...

関数呼び出しのオーバヘッドがそろそろ効きはじめてくるので、後者の二つの 例を大きめの N で実行して比較してみます:

In [11]: %timeit -n3 -r100 convolve4.naive_convolve(f, g)
3 loops, best of 100: 5.97 ms per loop
In [12]: N = 1000
In [13]: f = np.arange(N*N, dtype=np.int).reshape((N,N))
In [14]: g = np.arange(81, dtype=np.int).reshape((9, 9))
In [17]: %timeit -n1 -r10 convolve3.naive_convolve(f, g)
1 loops, best of 10: 1.16 s per loop
In [18]: %timeit -n1 -r10 convolve4.naive_convolve(f, g)
1 loops, best of 10: 597 ms per loop

(ただし、このベンチマークでも、アレイを関数呼び出しの中でアロケーショ ンしているので、純粋なアルゴリズムの速度比較ではありません)

Warning

スピードには代償がつきものです。殊に、型づけしたオブジェクト (このサンプルでは f, g, h) の値を None にするのは 危険です。値が None であること自体は全く問題はないのですが、オ ブジェクトには None かどうかのチェックしかできません。それ以外の操 作 (アトリビュートのルックアップやインデクシング) を行うと、 (Python では例外を送出するところが) segfault やデータの破損を引き 起こします。

実際のルールを突き詰めると少々ややこしいのですが、総体として言える ことははっきりしています: 値が None にセットされないことを知ってい るオブジェクト以外は、型づけするのをやめましょう。

より汎用的なコード

こんな書き方もできます:

def naive_convolve(object[DTYPE_t, ndim=2] f, ...):

つまり、 np.ndarray の代わりに object を使う書き方です。 Python 3.0 では、この書き方によって、バッファインタフェースを備えたラ イブラリなら何でも上記のアルゴリズムと組み合わせられます。興味のある人 がいれば、 Python Imaging Library などへのサポートを、 Python 2.x でも 簡単に追加できるでしょう。

ただし、この書き方にはちょっとだけ速度上のペナルティがあります (データ 型が np.ndarray にセットされていると、コンパイル時にいくつか条 件が仮定されます。とりわけ、データは完全にひとかたまり (pure stride) で、 indirect モードで格納されていないと仮定します)。