我大体上同意 @chepner 和 @juanpa.arrivillaga 在 cmets 中提出的建议。 Numpy 是一个高性能库，它执行的底层计算确实是用 C 语言编写的。此外，语法很简洁，并且在 numpy 数组的所有元素上应用标量操作很简单。
但是，如果我们使用以下假设（并且可以容忍丑陋的代码），由于您的特定算法的结构方式，实际上有一种方法可以使用 cython 显着提高代码的性能：
- 您的数组都是一维的，因此迭代数组中的每个项目非常简单。我们不需要替换更困难的 numpy 函数，例如
numpy.dot，因为您代码中的所有操作仅将标量与矩阵结合起来。
- 虽然在 python 中使用
for 循环是不可想象的，但在 cython 中迭代每个索引是非常可行的。此外，最终输出中的每个项目仅取决于对应于该项目索引的输入（即第 0 个项目使用u[0]、PorosityProfile[0] 等）。
- 您对任何中间数组都不感兴趣，只对
compute_python 函数中返回的最终结果感兴趣。因此，为什么要浪费时间为所有这些中间 numpy 数组分配内存？
- 使用
x**y 语法出奇地慢。我使用gcc 编译器选项--ffast-math 来显着改善这一点。我还使用了几个 cython 编译器指令来避免 python 检查和开销。
- 创建 numpy 数组本身可能会产生 python 开销，因此我使用类型化的 memoryviews（无论如何都是首选的、更新的语法）和 malloc-ed 指针来创建输出数组，而无需与 python 进行太多交互（只有两行，得到输出大小和返回语句显示了显着的 python 交互，如 cython 注释文件中所示）。
考虑到所有这些因素，这里是修改后的代码。它的执行速度比我笔记本电脑上的朴素 python 版本快了近一个数量级。
sublimation.pyx
from libc.stdlib cimport malloc, free
def compute_cython(float[:] u, float[:] porosity_profile,
float[:] density_ice_profile, float[:] density_dust_profile,
float[:] density_profile):
cdef:
float dust_j, dust_f, dust_g, dust_h, dust_i
float ice_i, ice_c, ice_d, ice_e, ice_f, ice_g, ice_h
int size, i
float dt, result_dust, x, dust
float result_ice_numer, result_ice_denom, result_ice, ice
float* out
dust_j, dust_f, dust_g, dust_h, dust_i = \
250.0, 633.0, 2.513, -2.2e-3, -2.8e-6
ice_i, ice_c, ice_d, ice_e, ice_f, ice_g, ice_h = \
273.16, 1.843e5, 1.6357e8, 3.5519e9, 1.6670e2, 6.4650e4, 1.6935e6
size = len(u)
out = <float *>malloc(size * sizeof(float))
for i in range(size):
dt = u[i] - dust_j
result_dust = dust_f + (dust_g*dt) + (dust_h*dt**2) + (dust_i*dt**3)
x = u[i] / ice_i
result_ice_numer = x**3*(ice_c + ice_d*x**2 + ice_e*x**6)
result_ice_denom = 1 + ice_f*x**2 + ice_g*x**4 + ice_h*x**8
result_ice = result_ice_numer / result_ice_denom
ice = density_ice_profile[i]*result_ice
dust = density_dust_profile[i]*result_dust
out[i] = (dust + ice)/density_profile[i]
return <float[:size]>out
setup.py
from distutils.core import setup
from Cython.Build import cythonize
from distutils.core import Extension
def create_extension(ext_name):
global language, libs, args, link_args
path_parts = ext_name.split(".")
path = "./{0}.pyx".format("/".join(path_parts))
ext = Extension(ext_name, sources=[path], libraries=libs, language=language,
extra_compile_args=args, extra_link_args=link_args)
return ext
if __name__ == "__main__":
libs = []#no external c libraries in this case
language = "c"#chooses c rather than c++ since no c++ features were used
args = ["-w", "-O3", "-ffast-math"]#assumes gcc is the compiler
link_args = []#none here, could use -fopenmp for parallel code
annotate = True#autogenerates .html files per .pyx
directives = {#saves typing @cython decorators and applies them globally
"boundscheck": False,
"wraparound": False,
"initializedcheck": False,
"cdivision": True,
"nonecheck": False,
}
ext_names = [
"sublimation",
]
extensions = [create_extension(ext_name) for ext_name in ext_names]
setup(ext_modules = cythonize(
extensions,
annotate=annotate,
compiler_directives=directives,
)
)
main.py
import numpy as np
import sublimation as sub
def compute_python(u, PorosityProfile, DensityIceProfile, DensityDustProfile, DensityProfile):
DustJ, DustF, DustG, DustH, DustI = 250.0, 633.0, 2.513, -2.2e-3, -2.8e-6
IceI, IceC, IceD, IceE, IceF, IceG, IceH = 273.16, 1.843e5, 1.6357e8, 3.5519e9, 1.6670e2, 6.4650e4, 1.6935e6
delta = u-DustJ
result_dust = DustF+DustG*delta+DustH*delta**2+DustI*(delta**3)
x = u/IceI
result_ice = (x**3)*(IceC+IceD*(x**2)+IceE*(x**6))/(1+IceF*(x**2)+IceG*(x**4)+IceH*(x**8))
return (DensityIceProfile*result_ice+DensityDustProfile*result_dust)/DensityProfile
size = 100
u = np.random.rand(size).astype(np.float32)
porosity = np.random.rand(size).astype(np.float32)
ice = np.random.rand(size).astype(np.float32)
dust = np.random.rand(size).astype(np.float32)
density = np.random.rand(size).astype(np.float32)
"""
Run these from the terminal to out the performance!
python3 -m timeit -s "from main import compute_python, u, porosity, ice, dust, density" "compute_python(u, porosity, ice, dust, density)"
python3 -m timeit -s "from main import sub, u, porosity, ice, dust, density" "sub.compute_cython(u, porosity, ice, dust, density)"
python3 -m timeit -s "import numpy as np; from main import sub, u, porosity, ice, dust, density" "np.asarray(sub.compute_cython(u, porosity, ice, dust, density))"
The first command tests the python version. (10000 loops, best of 3: 45.5 usec per loop)
The second command tests the cython version, but returns just a memoryview object. (100000 loops, best of 3: 4.63 usec per loop)
The third command tests the cython version, but converts the result to a ndarray (slower). (100000 loops, best of 3: 6.3 usec per loop)
"""
如果我对这个答案的工作原理的解释中有任何不清楚的部分，请告诉我，希望对您有所帮助！
更新 1：
不幸的是，我无法让 MSYS2 和 numba（它依赖于 LLVM）相互配合，所以我无法进行任何直接比较。但是，按照@max9111 的建议，我将-march=native 添加到我的setup.py 文件中的args 列表中；但是，时间安排与以前没有显着差异。
从this great answer 看来，numpy 数组和类型化内存视图之间的自动转换似乎存在一些开销，这在初始函数调用中都会进行（如果将结果转换回，则在 return 语句中也是如此）。恢复使用这样的函数签名：
ctypedef np.float32_t DTYPE_t
def compute_cython_np(
np.ndarray[DTYPE_t, ndim=1] u,
np.ndarray[DTYPE_t, ndim=1] porosity_profile,
np.ndarray[DTYPE_t, ndim=1] density_ice_profile,
np.ndarray[DTYPE_t, ndim=1] density_dust_profile,
np.ndarray[DTYPE_t, ndim=1] density_profile):
每次调用为我节省了大约 1us，将它减少到大约 3.6us 而不是 4.6us，这有点重要，特别是如果要多次调用该函数。 当然，如果您打算多次调用该函数，则直接传入二维 numpy 数组可能会更有效，从而节省大量 Python 函数调用开销并摊销numpy array -> typed memoryview 转换的成本。 此外，使用 numpy 结构化数组可能会很有趣，它可以在 cython 中转换为结构的类型化内存视图，因为这可以将缓存中的所有数据更紧密地放在一起并加快内存访问时间。
作为之前在 cmets 中承诺的最后一点，这里有一个使用 prange 的版本，它利用了并行处理的优势。请注意，这只能与类型化的内存视图一起使用，因为 python 的 GIL 必须在 prange 循环内释放（并使用 args 和 link_args 的 -fopenmp 标志编译：
from cython.parallel import prange
from libc.stdlib cimport malloc, free
def compute_cython_p(float[:] u, float[:] porosity_profile,
float[:] density_ice_profile, float[:] density_dust_profile,
float[:] density_profile):
cdef:
float dust_j, dust_f, dust_g, dust_h, dust_i
float ice_i, ice_c, ice_d, ice_e, ice_f, ice_g, ice_h
int size, i
float dt, result_dust, x, dust
float result_ice_numer, result_ice_denom, result_ice, ice
float* out
dust_j, dust_f, dust_g, dust_h, dust_i = \
250.0, 633.0, 2.513, -2.2e-3, -2.8e-6
ice_i, ice_c, ice_d, ice_e, ice_f, ice_g, ice_h = \
273.16, 1.843e5, 1.6357e8, 3.5519e9, 1.6670e2, 6.4650e4, 1.6935e6
size = len(u)
out = <float *>malloc(size * sizeof(float))
for i in prange(size, nogil=True):
dt = u[i] - dust_j
result_dust = dust_f + (dust_g*dt) + (dust_h*dt**2) + (dust_i*dt**3)
x = u[i] / ice_i
result_ice_numer = x**3*(ice_c + ice_d*x**2 + ice_e*x**6)
result_ice_denom = 1 + ice_f*x**2 + ice_g*x**4 + ice_h*x**8
result_ice = result_ice_numer / result_ice_denom
ice = density_ice_profile[i]*result_ice
dust = density_dust_profile[i]*result_dust
out[i] = (dust + ice)/density_profile[i]
return <float[:size]>out
更新 2：
根据 cmets 中来自 @max9111 的非常有用的附加建议，我将代码中的所有 float[:] 声明切换为 float[::1]。这样做的意义在于它允许连续存储数据，而 cython 不需要担心元素之间存在跨步。这允许 SIMD 矢量化，从而显着进一步优化代码。以下是使用以下命令生成的更新时序：
python3 -m timeit -s "from main import compute_python, u, porosity, ice, dust, density" "compute_python(u, porosity, ice, dust, density)"
python3 -m timeit -s "import numpy as np; from main import sub, u, porosity, ice, dust, density" "np.asarray(sub.compute_cython(u, porosity, ice, dust, density))"
python3 -m timeit -s "import numpy as np; from main import sub, u, porosity, ice, dust, density" "np.asarray(sub.compute_cython_p(u, porosity, ice, dust, density))"
size = 100
python: 44.7 usec per loop
cython serial: 4.44 usec per loop
cython parallel: 111 usec per loop
cython serial contiguous: 3.83 usec per loop
cython parallel contiguous: 116 usec per loop
size = 1000
python: 167 usec per loop
cython serial: 16.4 usec per loop
cython parallel: 115 usec per loop
cython serial contiguous: 8.24 usec per loop
cython parallel contiguous: 111 usec per loop
size = 10000
python: 1.32 msec per loop
cython serial: 128 usec per loop
cython parallel: 142 usec per loop
cython serial contiguous: 55.5 usec per loop
cython parallel contiguous: 150 usec per loop
size = 100000
python: 19.5 msec per loop
cython serial: 1.21 msec per loop
cython parallel: 691 usec per loop
cython serial contiguous: 473 usec per loop
cython parallel contiguous: 274 usec per loop
size = 1000000
python: 211 msec per loop
cython serial: 12.3 msec per loop
cython parallel: 5.74 msec per loop
cython serial contiguous: 4.82 msec per loop
cython parallel contiguous: 1.99 msec per loop