Есть ли способ улучшить скорость кода на Cython? - PullRequest
Новая
       10

Есть ли способ улучшить скорость кода на Cython?

0 голосов
/ 11 апреля 2019

Я ищу, чтобы ускорить следующие коды питона: 

def fun_np(m,data):
a, b, c = data[:,0], data[:,1], data[:,2] 

M = len(data[:,0]) 
n = round((m+1)*(m+2)*(m+3)/6) 
u =np.zeros((M,n))

C = 0
for i in range(0,m+1):
    for j in range(0,i+1):
        for k in range(0,j+1):
            if ((i-j)!=0):
                u[:,C] = (j-k)*(a)**(i-j)*(b)**(j-k-1)*(c)**k

        C=C+1  
return u

соответствующие коды Cython: 

%%cython 
import numpy as np
cimport numpy as np
from cython import wraparound, boundscheck, nonecheck

@boundscheck(False)
@wraparound(False)
@nonecheck(False)

cpdef fun_cyt(int m,np.ndarray[np.float64_t, ndim=2] data):

cdef:
    np.ndarray[np.float64_t, ndim=1] a = data[:,0]
    np.ndarray[np.float64_t, ndim=1] b = data[:,1]
    np.ndarray[np.float64_t, ndim=1] c = data[:,2]
    int M, n
    Py_ssize_t i, j, k, s
M = len(data[:,0]) 
n = round((m+1)*(m+2)*(m+3)/6)   
cdef np.ndarray[np.float64_t, ndim=2]  u = np.zeros((M,n), dtype=np.float64)

cdef int C = 0
for i in range(m+1): #range(0,m+1):
    for j in range(i+1):
        for k in range(j+1):
            for s in range(M):
                if (i-j)!=0:
                    u[s,C] = (j-k)*(a[s])**(i-j)*(b[s])**(j-k-1)*(c[s])**k

            C=C+1
return u

Вот время 

z = np.random.randn(6000, 3); m=20;

%timeit fun_np(m,z);

результат: 1,97 с ± 11,2 мс на цикл (среднее ± стандартное отклонение из 7 циклов, по 1 циклу каждый) 

%timeit fun_cyt(m,z);

результат: 1,91 с ± 12,7 мс на цикл (среднее ± стандартное отклонение из 7 циклов, по 1 циклу каждый)

Как видите, между кодами numpy и cython нет значимости. Буду признателен, если вы сможете помочь оптимизировать коды Cython, если это возможно.

Аннотированный HTML кодов Cython HTML

Ответы [ 2 ]

1 голос
/ 17 апреля 2019

Очень интересный пример!Большинство операций выполняются на 6000 элементных векторах.Cython не может быть быстрее, чем numpy, когда дело касается большой векторной мощности, умножения и сложения.Вы можете быть такими же быстрыми, как NumPy, внедрив это в Cython, и, возможно, даже получить от 10% до 20%, удалив некоторые накладные расходы NUMPY.

Однако есть и другие способы ускорить этот расчет.Операции с вектором - это операции с тремя столбцами вашего вектора данных, и вы записываете в столбцы выходного вектора.По умолчанию у пустых массивов есть порядок следования строк, то есть в памяти строки расположены в памяти рядом.Для операций, сделанных здесь, это плохо.Далее читаем: https://en.wikipedia.org/wiki/Row-_and_column-major_order.

Две функции в основном одинаковы, они будут идентичны, если создание выходного вектора будет происходить вне функции.

Обратите внимание на следующее: я заменилu [:, C] = ... с u [:, C] + =, потому что в противном случае результат определяется только k = j и, следовательно, всегда 0. Я не знаю, в чем смысл этих вычислений, но чтоВероятно, это не так. 

import numpy as np
def fun_np(m,data):
    a, b, c = data[:,0], data[:,1], data[:,2] 

    M = len(data[:,0]) 
    n = round((m+1)*(m+2)*(m+3)/6) 
    u = np.zeros((M,n))

    C = 0
    for i in range(0,m+1):
        for j in range(0,i+1):
            for k in range(0,j+1):
                if ((i-j)!=0):
                    u[:,C] += (j-k)*(a)**(i-j)*(b)**(j-k-1)*(c)**k

            C=C+1  
    return u

def fun_npF(m,data):
    a, b, c = data[:,0], data[:,1], data[:,2] 

    M = len(data[:,0]) 
    n = round((m+1)*(m+2)*(m+3)/6) 
    u = np.zeros((M,n),order='F')

    C = 0
    for i in range(0,m+1):
        for j in range(0,i+1):
            for k in range(0,j+1):
                if ((i-j)!=0):
                    u[:,C] += (j-k)*(a)**(i-j)*(b)**(j-k-1)*(c)**k

            C=C+1  
    return u

z = np.random.randn(6000, 3); m=20;
print("Numpy Row-major")
%timeit fun_np(m,z);

# Fortran order, because vector operations on columns
print("Numpy Column-major")
zF = np.asarray(z.copy(),order='F')
%timeit fun_npT(m,zF);

# Check if output the same
diff = (max(np.ravel(abs(fun_np(m,z)-fun_npF(m,zF)))))
max_rm = (max(np.ravel(abs(fun_np(m,z)))))
max_cm = (max(np.ravel(abs(fun_npF(m,zF)))))
print("Dffference: %f, Max value Row-major: %f, Max value Column-major: %f"%(diff, max_rm, max_cm))

Это дает мне 

Numpy Row-major
1.64 s ± 12.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
Numpy Column-major
16 ms ± 345 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Dffference: 0.000000, Max value Row-major: 196526643123.792450, Max value Column-major: 196526643123.792450

Вы можете получить еще больше, если подумаете о том, где поставить «если» и объедините это с Cython, но опять же толькопримерно на 10-20% я бы догадался.

1 голос
/ 12 апреля 2019

Как уже упоминалось в комментариях, вы можете попробовать это с помощью Numba.Я бы порекомендовал дополнительно распараллелить циклы: 

from numba import prange, jit

@jit(nopython=True, parallel=True)
def fun_numba(m,data):
    a, b, c = data[:,0], data[:,1], data[:,2] 

    M = len(data[:,0]) 
    n = round((m+1)*(m+2)*(m+3)/6) 
    u = np.zeros((M,n))

    C = 0
    for i in range(0,m+1):
        for j in range(0,i+1):
            for k in prange(0,j+1):
                if ((i-j)!=0):
                    u[:,C] = (j-k)*(a)**(i-j)*(b)**(j-k-1)*(c)**k

            C=C+1  
    return u

дает мне на моей машине: 

In [11]: %timeit fun_np(m,z)                                                                         
642 ms ± 4.13 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [12]: %timeit fun_numba(m,z)                                                                      
101 ms ± 7.15 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
...