n-D 배열을 반복하는 일반 함수

Question

Cython을 사용하면 다른 차원의 배열에서도 사용할 수있는 빠른 일반 함수를 작성할 수 있습니까? dealiasing 기능이 간단한 경우에 대한 예를 들어 : 여기서

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. 

, I는 일차원 2 차원 및 3 차원의 경우 세 가지 기능을 필요로한다. 모든 (합리적인) 차원에서 작업을 수행하는 함수를 작성하는 좋은 방법이 있습니까?

추신 : 여기, 나는 memoryviews를 사용하려고 시도했지만, 이것이 올바른 방법이라고 확신하지 못합니다. 나는 if where_dealiased[i0]: data[ik, i0] = 0.이라는 줄이 cython -a 명령에 의해 작성된 주석이 달린 HTML에서 흰색이 아니라는 사실에 놀랐다. 뭔가 잘못 됐니?

Answer 1

가장 먼저 말할 것은 3 가지 기능을 유지하고자하는 이유가 있습니다.보다 일반적인 기능을 사용하면 cython 컴파일러와 c 컴파일러 모두 최적화를 놓칠 수 있습니다.

3 개의 함수를 래핑하는 하나의 함수를 만드는 것은 매우 어렵습니다. 두 개의 배열을 파이썬 객체로 가져 와서 모양을 확인하고 관련 다른 함수를 호출합니다.

그러나

new axis 표기법을 사용하여 더 높은 차원의 배열로 개주 난 그냥 낮은 차원의 배열 한 후 가장 높은 차원의 함수를 작성하고, 시도 할 것이 무엇 다음이 시도하려고 한 경우 :

을

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]을 사용하여 배열의 끝에서부터 시작 (줄 바꿈) 대신 색인 할 수 있습니다.

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

: 파일의 시작에 주석을 포함 할 수 당신은 이러한 추가 파이썬 기능을 해제 할 수 있습니다 예를 전체 파일이 추가 코드를 제거하기 위해, compiler directives의 사용을 통해 C 배열 기능에 감소

컴파일러 지시문은 데코레이터를 사용하여 함수 수준에서 적용 할 수도 있습니다. 이 문서에 설명되어 있습니다.

Answer 2

당신은 일반적인 방식으로 읽을 수있는 np.ndarray 객체 strided 특성과 numpy.ndindex()을 사용하여 평면화 어레이와 같은 위치에 의해 결정된다 : strides가 쉽게 (strides*indices).sum() 하 이루어진다

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

, 일차원 배열

#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 경우 가 일본어 배열을 복사하지이 알고리즘 : 아래의 코드는 동작 예를 구성하는 방법을 도시

False 
True 

:에

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 

결과