Cython: медленные массивы numpy

Question

Я пытаюсь ускорить мой код с помощью cython. После перевода кода на cython из python я вижу, что я не получил никакой скорости. Я думаю, что источником проблемы является плохая производительность, которую я получаю с помощью массивов numpy в cython.

я придумал очень простую программу, чтобы показать это:

############### test.pyx ################# 
import numpy as np 
cimport numpy as np 
cimport cython 

def func1(long N): 

    cdef double sum1,sum2,sum3 
    cdef long i 

    sum1 = 0.0 
    sum2 = 0.0 
    sum3 = 0.0 

    for i in range(N): 
     sum1 += i 
     sum2 += 2.0*i 
     sum3 += 3.0*i 

    return sum1,sum2,sum3 

def func2(long N): 

    cdef np.ndarray[np.float64_t,ndim=1] sum_arr 
    cdef long i 

    sum_arr = np.zeros(3,dtype=np.float64) 

    for i in range(N): 
     sum_arr[0] += i 
     sum_arr[1] += 2.0*i 
     sum_arr[2] += 3.0*i 

    return sum_arr 

def func3(long N): 

    cdef double sum_arr[3] 
    cdef long i 

    sum_arr[0] = 0.0 
    sum_arr[1] = 0.0 
    sum_arr[2] = 0.0 

    for i in range(N): 
     sum_arr[0] += i 
     sum_arr[1] += 2.0*i 
     sum_arr[2] += 3.0*i 

    return sum_arr 
########################################## 

################## test.py ############### 
import time 
import test as test 

N = 1000000000 

for i in xrange(10): 

    start = time.time() 
    sum1,sum2,sum3 = test.func1(N) 
    print 'Time taken = %.3f'%(time.time()-start) 

print '\n' 
for i in xrange(10): 
    start = time.time() 
    sum_arr = test.func2(N) 
    print 'Time taken = %.3f'%(time.time()-start) 

print '\n' 
for i in xrange(10): 
    start = time.time() 
    sum_arr = test.func3(N) 
    print 'Time taken = %.3f'%(time.time()-start) 
############################################ 

И из питона test.py я получаю:

Time taken = 1.445 
Time taken = 1.433 
Time taken = 1.434 
Time taken = 1.428 
Time taken = 1.449 
Time taken = 1.425 
Time taken = 1.421 
Time taken = 1.451 
Time taken = 1.483 
Time taken = 1.418 

Time taken = 2.623 
Time taken = 2.603 
Time taken = 2.977 
Time taken = 3.237 
Time taken = 2.748 
Time taken = 2.798 
Time taken = 2.811 
Time taken = 2.783 
Time taken = 2.585 
Time taken = 2.595 

Time taken = 1.503 
Time taken = 1.529 
Time taken = 1.509 
Time taken = 1.543 
Time taken = 1.427 
Time taken = 1.425 
Time taken = 1.423 
Time taken = 1.415 
Time taken = 1.414 
Time taken = 1.418 

Мой вопрос: почему func2 почти 2x медленнее чем func1 и func3?

Есть ли способ улучшить это?

Спасибо!

######## UPDATE

Моя настоящая проблема заключается в следующем. Я вызываю функцию, которая принимает 3D-массив (например, P [i, j, k]). Функция будет проходить через каждый элемент и вычислять несколько величин: величину, которая зависит от значения массива в этой позиции (например, A = f (P [i, j, k])) и других величин, которые зависят только от позиции самого массива (B = g (i, j, k)). Схематически все будет выглядеть так:

for i in xrange(N): 
    corr1 = h(i,val) 

    for j in xrange(N): 
     corr2 = h(j,val) 

     for k in xrange(N): 
      corr3 = h(k,val) 

      A = f(P[i,j,k]) 
      B = g(i,j,k) 
      Arr[B] += A*corr1*corr2*corr3 

где val является свойством 3D-массива, представленного рядом. Это число может быть различным для разных полей.

Поскольку я должен выполнять эту операцию по многим 3D-массивам, я бы подумал, что было бы лучше, если бы я создал новую процедуру, которая принимает множество различных входных 3D-массивов, оставляя количество неизвестных априорных массивов. Идея состоит в том, что, поскольку B будет точно таким же во всех массивах, я могу избежать вычисления его для каждого массива и только вычислить его один раз. Проблема заключается в том, что corr1, corr2, corr3 выше станет массивами:

Если у меня есть ряд 3D-массивы, равных num_3D_arrays я делаю что-то как:

for i in xrange(N): 
    for p in xrange(num_3D_arrays): 
     corr1[p] = h(i,val[p]) 

    for j in xrange(N): 
     for p in xrange(num_3D_arrays): 
      corr2[p] = h(j,val[p]) 

     for k in xrange(N): 
      for p in xrange(num_3D_arrays): 
       corr3[p] = h(k,val[p]) 

      B = g(i,j,k) 
      for p in xrange(num_3D_arrays): 
       A[p] = f(P[i,j,k]) 
       Arr[p,B] += A[p]*corr1[p]*corr2[p]*corr3[p] 

Так вал, что я меняюсь переменные corr1, corr2, corr3 и A от скаляра к массивам убивают производительность, которую я ожидал бы избежать большого цикла.

#

Answer 1

• почему func2 почти 2x медленнее, чем func1?

  Это потому, что индексирование вызывает косвенное отношение, поэтому вы удваиваете количество элементарных операций. Рассчитать сумму, как в func1, а затем повлиять на sum=array([sum1,sum2,sum3])

• Как ускорить код python?

  1. Numpy - первая хорошая идея, она поднимет почти скорость C без усилий.

  2. Numba может заполнить разрыв без усилий, и это очень просто.

  3. Cython для критических случаев.

Вот некоторые иллюстрации, что:

# python way 
def func1(N): 
    sum1 = 0.0 
    sum2 = 0.0 
    sum3 = 0.0 

    for i in range(N): 
     sum1 += i 
     sum2 += 2.0*i 
     sum3 += 3.0*i 

    return sum1,sum2,sum3 

# numpy way 
def func2(N): 
    aran=arange(float(N)) 
    sum1=aran.sum() 
    sum2=(2.0*aran).sum() 
    sum3=(3.0*aran).sum() 
    return sum1,sum2,sum3 

#numba way 
import numba  
func3 =numba.njit(func1) 

""" 
In [609]: %timeit func1(10**6) 
1 loop, best of 3: 710 ms per loop 

In [610]: %timeit func2(1e6) 
100 loops, best of 3: 22.2 ms per loop 

In [611]: %timeit func3(10e6) 
100 loops, best of 3: 2.87 ms per loop 
""" 

Answer 2

Посмотрите на html производимого cython -a ...pyx.

Для func1, то sum1 += i линия расширяется:

+15:   sum1 += i 
    __pyx_v_sum1 = (__pyx_v_sum1 + __pyx_v_i); 

для func3, с массивом C

+45:   sum_arr[0] += i 
    __pyx_t_3 = 0; 
    (__pyx_v_sum_arr[__pyx_t_3]) = ((__pyx_v_sum_arr[__pyx_t_3]) + __pyx_v_i); 

Немного сложнее, но прямо вперед c.

Но func2:

+29:   sum_arr[0] += i 
    __pyx_t_12 = 0; 
    __pyx_t_6 = -1; 
    if (__pyx_t_12 < 0) { 
     __pyx_t_12 += __pyx_pybuffernd_sum_arr.diminfo[0].shape; 
     if (unlikely(__pyx_t_12 < 0)) __pyx_t_6 = 0; 
    } else if (unlikely(__pyx_t_12 >= __pyx_pybuffernd_sum_arr.diminfo[0].shape)) __pyx_t_6 = 0; 
    if (unlikely(__pyx_t_6 != -1)) { 
     __Pyx_RaiseBufferIndexError(__pyx_t_6); 
     __PYX_ERR(0, 29, __pyx_L1_error) 
    } 
    *__Pyx_BufPtrStrided1d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_sum_arr.rcbuffer->pybuffer.buf, __pyx_t_12, __pyx_pybuffernd_sum_arr.diminfo[0].strides) += __pyx_v_i; 

Гораздо сложнее со ссылками на numpy функций (например Pyx_BUfPtrStrided1d). Даже инициализируешь массив сложно:

+26:  sum_arr = np.zeros(3,dtype=np.float64) 
    __pyx_t_1 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 26, __pyx_L1_error) 
    __Pyx_GOTREF(__pyx_t_1); 
    .... 

Я ожидаю, что перемещение sum_arr создания вызывающему Python, и передать его в качестве аргумента func2 бы сэкономить время.

Вы читали это руководство для использования memoryviews:

http://cython.readthedocs.io/en/latest/src/userguide/memoryviews.html

Вы получите лучший cython производительность, если сосредоточиться на написании операции низкого уровня, чтобы они переводят на простой c. В

for k in xrange(N): 
     corr3 = h(k,val) 

     A = f(P[i,j,k]) 
     B = g(i,j,k) 
     Arr[B] += A*corr1*corr2*corr3 

Это не петли на i,j,k, который замедляет работу. Он оценивает h, f и g каждый раз, а также Arr[B] +=.... Эти функции должны быть плотно закодированы cython, а не общие функции Python. Посмотрите на скомпилированную простоту функции sum3d в руководстве memoryview.

Answer 3

Есть несколько вещей, которые вы можете сделать, чтобы ускорить индексацию массива в Cython:

• Поворот bounds checking and wraparound.
• Использование typed memoryviews.
• Объявить массив contiguous.

Так что для вашей функции:

@cython.boundscheck(False) 
@cython.wraparound(False) 
def func2(long N): 

    cdef np.float64_t[::1] sum_arr 
    cdef long i 

    sum_arr = np.zeros(3,dtype=np.float64) 

    for i in range(N): 
     sum_arr[0] += i 
     sum_arr[1] += 2.0*i 
     sum_arr[2] += 3.0*i 

    return sum_arr 

Для исходного кода Cython производится следующий код C для линии sum_arr[0] += i:

__pyx_t_12 = 0; 
__pyx_t_6 = -1; 
if (__pyx_t_12 < 0) { 
    __pyx_t_12 += __pyx_pybuffernd_sum_arr.diminfo[0].shape; 
    if (unlikely(__pyx_t_12 < 0)) __pyx_t_6 = 0; 
} else if (unlikely(__pyx_t_12 >= __pyx_pybuffernd_sum_arr.diminfo[0].shape)) __pyx_t_6 = 0; 
if (unlikely(__pyx_t_6 != -1)) { 
    __Pyx_RaiseBufferIndexError(__pyx_t_6); 
    {__pyx_filename = __pyx_f[0]; __pyx_lineno = 13; __pyx_clineno = __LINE__; goto __pyx_L1_error;} 
} 
*__Pyx_BufPtrStrided1d(__pyx_t_5numpy_float64_t *, __pyx_pybuffernd_sum_arr.rcbuffer->pybuffer.buf, __pyx_t_12, __pyx_pybuffernd_sum_arr.diminfo[0].strides) += __pyx_v_i; 

С улучшениями выше:

__pyx_t_8 = 0; 
*((double *) (/* dim=0 */ ((char *) (((double *) __pyx_v_sum_arr.data) + __pyx_t_8)))) += __pyx_v_i;