Сравнение ускорителей Python (Cython, Numba, f2py) с einsum Numpy

Я сравниваю ускорители Python (Numba, Cython, f2py) с простыми циклами For и einsum Numpy для конкретной задачи (см. ниже). На данный момент Numpy является самым быстрым для этой проблемы (в 6 раз быстрее), но я хотел получить отзывы, если есть дополнительные оптимизации, которые я должен попробовать, или если я делаю что-то не так. Этот простой код основан на более крупном коде, в котором есть несколько вызовов einsum, но нет явных циклов for. Я проверяю, может ли какой-либо из этих ускорителей работать лучше.

Время выполнено с Python 2.7.9 в Mac OS X Yosemite с установленным gcc-5.3.0 (--with-fortran --without-multilib) из Homebrew. Также звонил %timeit; эти тайминги одиночных вызовов довольно точны.

In [1]: %run -i test_numba.py
test_numpy: 0.0805640220642
Matches Numpy output: True

test_dumb: 1.43043899536
Matches Numpy output: True

test_numba: 0.464295864105
Matches Numpy output: True

test_cython: 0.627640008926
Matches Numpy output: True

test_f2py: 5.01890516281
Matches Numpy output: True

test_f2py_order: 2.31424307823
Matches Numpy output: True

test_f2py_reorder: 0.507861852646
Matches Numpy output: True

Основной код:

import numpy as np
import numba
import time
import test_f2py as tf2py
import pyximport
pyximport.install(setup_args={'include_dirs':np.get_include()})
import test_cython as tcyth

def test_dumb(f,b):
    fnew = np.empty((f.shape[1],f.shape[2]))
    for i in range(f.shape[0]):
        for l in range(f.shape[3]):
            fnew += f[i,:,:,l] * b[i,l]
    return fnew


def test_dumber(f,b):
    fnew = np.empty((f.shape[1],f.shape[2]))
    for i in range(f.shape[0]):
        for j in range(f.shape[1]):
            for k in range(f.shape[2]):
                for l in range(f.shape[3]):
                    fnew[j,k] += f[i,j,k,l] * b[i,l]
    return fnew

@numba.jit(nopython=True)
def test_numba(f,b):
    fnew = np.zeros((f.shape[1],f.shape[2])) #NOTE: can't be empty, gives errors
    for i in range(f.shape[0]):
        for j in range(f.shape[1]):
            for k in range(f.shape[2]):
                for l in range(f.shape[3]):
                    fnew[j,k] += f[i,j,k,l] * b[i,l]
    return fnew

def test_numpy(f,b):
    return np.einsum('i...k,ik->...',f,b)

def test_f2py(f,b):
    return tf2py.test_f2py(f,b)

def test_f2py_order(f,b):
    return tf2py.test_f2py(f,b)

def test_f2py_reorder(f,b):
    return tf2py.test_f2py_reorder(f,b)

def test_cython(f,b):
    return tcyth.test_cython(f,b)

if __name__ == '__main__':

    #goal is to create: fnew = sum f*b over dim 0 and 3.
    f = np.random.rand(32,33,2000,64)
    b = np.random.rand(32,64)

    f1 = np.asfortranarray(f)
    b1 = np.asfortranarray(b)

    f2 = np.asfortranarray(np.transpose(f,[1,2,0,3]))

    funcs = [test_dumb,test_numba, test_cython, \
            test_f2py,test_f2py_order,test_f2py_reorder]

    tstart = time.time()    
    fnew_numpy= test_numpy(f,b)
    tstop = time.time()
    print test_numpy.__name__+': '+str(tstop-tstart)
    print 'Matches Numpy output: '+str(np.allclose(fnew_numpy,fnew_numpy))
    print ''

    for func in funcs:
        tstart = time.time()
        if func.__name__ == 'test_f2py_order':
            fnew = func(f1,b1)
        elif func.__name__ == 'test_f2py_reorder':
            fnew = func(f2,b1)
        else:
            fnew = func(f,b)
        tstop = time.time()
        print func.__name__+': '+str(tstop-tstart)
        print 'Matches Numpy output: '+str(np.allclose(fnew,fnew_numpy))
        print ''

Файл f2py (скомпилированный с помощью f2py -c -m test_f2py test_f2py.F90):

!file: test_f2py
subroutine test_f2py(f,b,fnew,n1,n2,n3,n4)

integer  :: n1,n2,n3,n4
real(8), dimension(n1,n2,n3,n4) :: f
real(8), dimension(n1,n4) :: b
real(8), dimension(n2,n3) :: fnew
!f2py intent(in) f
!f2py intent(in) b
!f2py intent(out) fnew
!f2py intent(in) n1
!f2py intent(in) n2
!f2py intent(in) n3
!f2py intent(in) n4

integer :: i1,i2,i3,i4

do i1=1,n1
    do i2=1,n2
        do i3=1,n3
            do i4=1,n4
                fnew(i2,i3) = fnew(i2,i3) + f(i1,i2,i3,i4)*b(i1,i4)
            enddo
        enddo
    enddo
enddo

end subroutine test_f2py

subroutine test_f2py_reorder(f,b,fnew,n1,n2,n3,n4)

integer  :: n1,n2,n3,n4
real(8), dimension(n1,n2,n3,n4) :: f
real(8), dimension(n3,n4) :: b
real(8), dimension(n1,n2) :: fnew
!f2py intent(in) f
!f2py intent(in) b
!f2py intent(out) fnew
!f2py intent(in) n1
!f2py intent(in) n2
!f2py intent(in) n3
!f2py intent(in) n4

integer :: i1,i2,i3,i4

do i3=1,n3
    do i4=1,n4
        do i1=1,n1
            do i2=1,n2
                fnew(i1,i2) = fnew(i1,i2) + f(i1,i2,i3,i4)*b(i3,i4)
            enddo
        enddo
    enddo
enddo

end subroutine test_f2py_reorder

И файл Cython .pyx (скомпилированный с помощью pyximport в основной процедуре):

#/usr/bin python
import numpy as np
cimport numpy as np

def test_cython(np.ndarray[np.float64_t,ndim=4] f, np.ndarray[np.float64_t,ndim=2] b):
    # cdef np.ndarray[np.float64_t,ndim=4] f
    # cdef np.ndarray[np.float64_t,ndim=2] b
    cdef np.ndarray[np.float64_t,ndim=2] fnew = np.empty((f.shape[1],f.shape[2]),dtype=np.float64)
    cdef int i,j,k,l
    cdef int Ni = f.shape[0]
    cdef int Nj = f.shape[1]
    cdef int Nk = f.shape[2]
    cdef int Nl = f.shape[3]

    for i in range(Ni):
        for j in range(Nj):
            for k in range(Nk):
                for l in range(Nl):
                    fnew[j,k] += f[i,j,k,l] * b[i,l]
    return fnew