性能比较Fortran,Numpy,Cython和Numexpr

Question

我有以下功能:

def get_denom(n_comp,qs,x,cp,cs):
'''
len(n_comp) = 1 # number of proteins
len(cp) = n_comp # protein concentration
len(qp) = n_comp # protein capacity
len(x) = 3*n_comp + 1 # fit parameters
len(cs) = 1

'''
    k = x[0:n_comp]
    sigma = x[n_comp:2*n_comp]
    z = x[2*n_comp:3*n_comp]

    a = (sigma + z)*( k*(qs/cs)**(z-1) )*cp
    denom = np.sum(a) + cs
    return denom

我将它与Fortran实现(我的第一个Fortran函数)进行比较:

subroutine get_denom (qs,x,cp,cs,n_comp,denom)

! Calculates the denominator in the SMA model (Brooks and Cramer 1992)
! The function is called at a specific salt concentration and isotherm point
! I loops over the number of components

implicit none

! declaration of input variables
integer, intent(in) :: n_comp ! number of components
double precision, intent(in) :: cs,qs ! salt concentration, free ligand concentration
double precision, dimension(n_comp), INTENT(IN) ::cp ! protein concentration
double precision, dimension(3*n_comp + 1), INTENT(IN) :: x ! parameters

! declaration of local variables
double precision, dimension(n_comp) :: k,sigma,z
double precision :: a
integer :: i

! declaration of outpur variables
double precision, intent(out) :: denom

k = x(1:n_comp) ! equlibrium constant
sigma = x(n_comp+1:2*n_comp) ! steric hindrance factor
z = x(2*n_comp+1:3*n_comp) ! charge of protein

a = 0.
do i = 1,n_comp
    a = a + (sigma(i) + z(i))*(k(i)*(qs/cs)**(z(i)-1.))*cp(i)
end do

denom = a + cs

end subroutine get_denom

我使用以下命令编译.f95文件:

1) f2py -c -m get_denom get_denom.f95 --fcompiler=gfortran

2)f2py -c -m get_denom_vec get_denom.f95 --fcompiler=gfortran --f90flags='-msse2'(最后一个选项应该打开自动矢量化)

我通过以下方式测试函数:

import numpy as np
import get_denom as fort_denom
import get_denom_vec as fort_denom_vec
from matplotlib import pyplot as plt
%matplotlib inline

def get_denom(n_comp,qs,x,cp,cs):
    k = x[0:n_comp]
    sigma = x[n_comp:2*n_comp]
    z = x[2*n_comp:3*n_comp]
    # calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992)
    a = (sigma + z)*( k*(qs/cs)**(z-1) )*cp
    denom = np.sum(a) + cs
    return denom

n_comp = 100
cp = np.tile(1.243,n_comp)
cs = 100.
qs = np.tile(1100.,n_comp)
x= np.random.rand(3*n_comp+1)
denom = np.empty(1)
%timeit get_denom(n_comp,qs,x,cp,cs)
%timeit fort_denom.get_denom(qs,x,cp,cs,n_comp)
%timeit fort_denom_vec.get_denom(qs,x,cp,cs,n_comp)

我添加了以下Cython代码:

import cython
# import both numpy and the Cython declarations for numpy
import numpy as np
cimport numpy as np

@cython.boundscheck(False)
@cython.wraparound(False)
def get_denom(int n_comp,np.ndarray[double, ndim=1, mode='c'] qs, np.ndarray[double, ndim=1, mode='c'] x,np.ndarray[double, ndim=1, mode='c'] cp, double cs):

    cdef int i
    cdef double a
    cdef double denom   
    cdef double[:] k = x[0:n_comp]
    cdef double[:] sigma = x[n_comp:2*n_comp]
    cdef double[:] z = x[2*n_comp:3*n_comp]
    # calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992)
    a = 0.
    for i in range(n_comp):
    #a += (sigma[i] + z[i])*( pow( k[i]*(qs[i]/cs), (z[i]-1) ) )*cp[i]
        a += (sigma[i] + z[i])*( k[i]*(qs[i]/cs)**(z[i]-1) )*cp[i]

    denom = a + cs

    return denom

编辑:

使用一个线程添加Numexpr:

def get_denom_numexp(n_comp,qs,x,cp,cs):
    k = x[0:n_comp]
    sigma = x[n_comp:2*n_comp]
    z = x[2*n_comp:3*n_comp]
    # calculates the denominator in Equ 14a - 14c (Brooks & Cramer 1992)
    a = ne.evaluate('(sigma + z)*( k*(qs/cs)**(z-1) )*cp' )
    return cs + np.sum(a)

ne.set_num_threads(1)  # using just 1 thread
%timeit get_denom_numexp(n_comp,qs,x,cp,cs)

结果是(越小越好):

为什么随着阵列数量的增加,Fortran的速度越来越接近Numpy？我怎么能加速Cython？使用指针？

Answer 1

暂停了

好的，最后，我们被允许在我们的其中一个盒子上安装Numpy等，这可能对您的原始帖子进行了全面的解释。

简短的答案是，从某种意义上说，您的原始问题是“错误的问题”。此外，其中一个贡献者还存在许多令人不安的滥用和错误信息，这些错误和虚假行为值得关注，以免任何人误以为是错误的信念和代价。

另外，出于以下原因和适当性，我决定将此答案作为单独的答案提交，而不是编辑4月14日的答案。

A部分：OP的答案

首先，处理原始帖子中的问题：您可能还记得我只能在Fortran方面发表评论，因为我们的政策严格限制可以安装哪些软件以及在计算机上的位置，并且我们没有Python等。 （直到现在）。我也曾多次指出，就我们所称的弯曲特征或“凹度”而言，结果的特征很有趣。

此外，仅根据“相对”结果（因为您没有发布绝对时间，并且那时我没有Numpy），我曾几次表示某些重要信息可能会潜伏其中。

确实是这样。

首先，我想确保可以重现您的结果，因为我们通常不使用Python / F2py，所以您的结果中隐含了什么编译器设置等并不明显，因此我进行了各种测试以确保我们一对一交谈（我4月14日的回答证明Debug与Release / O2的区别很大）。

图1显示了我在以下三种情况下的Python结果：您的原始Python / Numpy内部子程序（称为P，我刚刚剪切/粘贴了您的原始），您已发布的基于Do的原始Fortran s / r（称为此） FDo，我只是像以前一样复制/粘贴了您的原始文件，而我之前建议的一种变体就是依靠Array Sections，因此需要Sum（）（调用此FSAS，是通过编辑原始FDo创建的）。图1显示了通过timeit的绝对时序。

图2显示了基于您的相对策略除以Python / Numpy（P）时序的相对结果。仅显示两个（相对）Fortran变体。

显然，这些内容可以在您的原始情节中再现角色，并且我们可能相信我的测试似乎与您的测试一致。

现在，您最初的问题是“为什么随着数组大小的增加，Fortran的速度越来越接近Numpy？”。

其实不是。纯粹是纯粹依靠可能给人以“相对”时机的伪影/失真。

从图1可以看到，绝对的时间安排，很明显Numpy和Fortran正在分歧。因此，实际上，如果您愿意，Fortran的结果正在偏离Numpy，反之亦然。

一个更好的问题（为什么我首先将这些向上弯曲（例如，比较线性）？我以前的仅Fortran结果显示相对性能比率（几乎是平坦的）（4月14日图表/答案中的黄线），因此我推测Python方面发生了一些有趣的事情，建议进行检查。

显示这一点的一种方法是采用另一种相对度量。我将每个（绝对）序列除以自己的最大值（即n_comp = 10k），以了解这种“内部相对”性能如何发展（这些值称为10k值，代表n_comp = 10,000的时间） 。图3将P，FDo和FSAS的这些结果显示为P / P10k，FDo / FDo10k和FSAS / FSAS10k。为了清楚起见，y轴具有对数刻度。显然，随着P结果的n_comp减少（例如，带红色圆圈的注释部分），Fortran变体的预成型相对更好。

换句话说，Fortran在减小数组大小方面非常（非线性）更为高效。不确定到底为什么Python会随着n_comp的减少而变得更糟...但是确实存在，并且可能是内部开销/设置等问题以及解释器与编译器之间的差异等问题。

因此，并不是说Fortran正在使用Python“追赶”，恰恰相反，它正在与Python保持距离（见图1）。但是，Python和Fortran之间的非线性差异（随着n_comp的减少）会产生具有明显反直觉和非线性特征的“相对”时序。

因此，随着n_comp的增加和每种方法“稳定”到或多或少的线性模式，曲线趋于平坦，表明它们的差异呈线性增长，并且相对比率稳定为近似常数（忽略内存争用，smp问题等）。 ）...这很容易看出是否允许n_comp> 10k，但是我的4月14日答案中的黄线已经显示了仅Fortran的s / r。

另外：我的偏好是创建自己的计时例程/功能。timeit似乎很方便，但是“黑匣子”内部正在发生很多事情。设置自己的循环和结构，并确定定时功能的属性/分辨率对于正确评估很重要。