3

Cython を使用して、異なる次元の配列に対して機能する高速な一般関数を作成する方法はありますか? たとえば、関数のエイリアスを解除するこの単純なケースの場合:

import numpy as np
cimport numpy as np

ctypedef np.uint8_t DTYPEb_t
ctypedef np.complex128_t DTYPEc_t


def dealiasing1D(DTYPEc_t[:, :] data, 
                 DTYPEb_t[:] where_dealiased):
    """Dealiasing data for 1D solvers."""
    cdef Py_ssize_t ik, i0, nk, n0

    nk = data.shape[0]
    n0 = data.shape[1]

    for ik in range(nk):
        for i0 in range(n0):
            if where_dealiased[i0]:
                data[ik, i0] = 0.


def dealiasing2D(DTYPEc_t[:, :, :] data, 
                 DTYPEb_t[:, :] where_dealiased):
    """Dealiasing data for 2D solvers."""
    cdef Py_ssize_t ik, i0, i1, nk, n0, n1

    nk = data.shape[0]
    n0 = data.shape[1]
    n1 = data.shape[2]

    for ik in range(nk):
        for i0 in range(n0):
            for i1 in range(n1):
                if where_dealiased[i0, i1]:
                    data[ik, i0, i1] = 0.


def dealiasing3D(DTYPEc_t[:, :, :, :] data, 
                 DTYPEb_t[:, :, :] where_dealiased):
    """Dealiasing data for 3D solvers."""
    cdef Py_ssize_t ik, i0, i1, i2, nk, n0, n1, n2

    nk = data.shape[0]
    n0 = data.shape[1]
    n1 = data.shape[2]
    n2 = data.shape[3]

    for ik in range(nk):
        for i0 in range(n0):
            for i1 in range(n1):
                for i2 in range(n2):
                    if where_dealiased[i0, i1, i2]:
                        data[ik, i0, i1, i2] = 0.

ここでは、1 次元、2 次元、3 次元の場合に 3 つの関数が必要です。すべての (妥当な) ディメンションに対して機能する関数を作成する良い方法はありますか?

PS:ここでは、メモリビューを使用しようとしましたが、これが正しい方法かどうかはわかりません。if where_dealiased[i0]: data[ik, i0] = 0.コマンドによって生成された注釈付きのhtmlで線が白くないことに驚いていcython -aます。何か問題がありますか？

4

2 に答える 2

2

私が最初に言うことは、3 つの関数を保持したい理由があるということです。より一般的な関数を使用すると、cython コンパイラと c コンパイラの両方からの最適化を逃す可能性があります。

これら 3 つの関数をラップする 1 つの関数を作成することは非常に実行可能です。2 つの配列を Python オブジェクトとして取得し、形状をチェックして、関連する他の関数を呼び出すだけです。

しかし、これを試みる場合は、最高次元の関数を作成し、新しい軸表記を使用して低次元の配列を高次元の配列として再キャストすることを試みます。

cdef np.uint8_t [:] a1d = np.zeros((256, ), np.uint8) # 1d
cdef np.uint8_t [:, :] a2d = a1d[None, :]             # 2d
cdef np.uint8_t [:, :, :] a3d = a1d[None, None, :]    # 3d
a2d[0, 100] = 42
a3d[0, 0, 200] = 108
print(a1d[100], a1d[200])
# (42, 108)

cdef np.uint8_t [:, :] data2d = np.zeros((128, 256), np.uint8) #2d
cdef np.uint8_t [:, :, :, :] data4d = data2d[None, None, :, :] #4d
data4d[0, 0, 42, 108] = 64
print(data2d[42, 108])
# 64

ご覧のとおり、メモリ ビューは高次元にキャストでき、それを使用して元のデータを変更できます。新しいビューを最高次元の関数に渡す前に、これらのトリックを実行するラッパー関数を書きたいと思うでしょう。このトリックはあなたのケースでは非常にうまくいくと思いますが、データでやりたいことができるかどうかを知るために試してみる必要があります.

あなたの PS: には、非常に簡単な説明があります。「余分なコード」は、インデックス エラー、タイプ エラーを生成し、[-1] を使用して配列の先頭ではなく末尾からインデックス付けできるようにするコードです (ラップアラウンド)。これらの余分な python 機能を無効にして、コンパイラ ディレクティブを使用して C 配列機能に減らすことができます。たとえば、ファイル全体からこの余分なコードを削除するには、ファイルの先頭にコメントを含めることができます。

# cython: boundscheck=False, wraparound=False, nonecheck=False

コンパイラ ディレクティブは、デコレータを使用して関数レベルで適用することもできます。ドクターが説明します。

于 2014-10-06T01:53:48.133 に答える
1

位置が次のように決定されるように、オブジェクトの属性を使用numpy.ndindex()してフラット化された配列を一般的な方法で読み取ることができます。stridednp.ndarray

indices[0]*strides[0] + indices[1]*strides[1] + ... + indices[n]*strides[n]

が1 次元配列である(strides*indices).sum()場合、これは簡単に実行できます。strides以下のコードは、実際の例を作成する方法を示しています。

#cython profile=True
#blacython wraparound=False
#blacython boundscheck=False
#blacython nonecheck=False
#blacython cdivision=True
cimport numpy as np
import numpy as np

def readNDArray(x):
    if not isinstance(x, np.ndarray):
        raise ValueError('x must be a valid np.ndarray object')
    if x.itemsize != 8:
        raise ValueError('x.dtype must be float64')
    cdef np.ndarray[double, ndim=1] v # view of x
    cdef np.ndarray[int, ndim=1] strides
    cdef int pos

    shape = list(x.shape)
    strides = np.array([s//x.itemsize for s in x.strides], dtype=np.int32)
    v = x.ravel()
    for indices in np.ndindex(*shape):
        pos = (strides*indices).sum()
        v[pos] = 2.
    return np.reshape(v, newshape=shape)

このアルゴリズムは、元の配列が C contiguous の場合、コピーしません。

def main():
    # case 1
    x = np.array(np.random.random((3,4,5,6)), order='F')
    y = readNDArray(x)
    print(np.may_share_memory(x, y))
    # case 2
    x = np.array(np.random.random((3,4,5,6)), order='C')
    y = readNDArray(x)
    print np.may_share_memory(x, y)
    return 0

結果:

False
True
于 2014-10-06T06:25:17.957 に答える