Я в целом согласен с советами, представленными @chepner и @ juanpa.arrivillaga в комментариях.Numpy является производительной библиотекой, и верно, что выполняемые ею базовые вычисления написаны на языке C. Кроме того, синтаксис чистый и тривиально применять скалярные операции ко всем элементам массива numpy.
Однакона самом деле есть способ значительно улучшить производительность вашего кода с помощью Cython благодаря структуре вашего конкретного алгоритма, если мы используем следующие допущения (и можем допускать некрасивый код):
- Ваши массивывсе одномерные, что делает итерации по каждому элементу в массиве чрезвычайно тривиальным.Нам не нужно заменять более сложные простые функции, такие как, например,
numpy.dot, поскольку все операции в вашем коде объединяют скаляры только с матрицами. - Хотя использование цикла
for в python было бы немыслимо, повторять всеИндекс очень выполним в Cython.Кроме того, каждый элемент в конечном выводе зависит только от входных данных, которые соответствуют индексу этого элемента (т. Е. 0-й элемент использует u[0], PorosityProfile[0] и т. Д.). - Вас не интересует ни одно изпромежуточные массивы, только в конечном результате, возвращенном в вашей функции
compute_python.Следовательно, зачем тратить время на выделение памяти для всех этих промежуточных массивов numpy? - Использование синтаксиса
x**y на удивление медленное.Я использую опцию компилятора gcc, --ffast-math, чтобы значительно улучшить это.Я также использую несколько директив компилятора Cython, чтобы избежать проверок и накладных расходов на Python. - Создание самих numpy массивов может иметь накладные расходы на Python, поэтому я использую комбинацию типизированных представлений памяти (в любом случае предпочтительный, более новый синтаксис) и указателей malloc-edсоздать выходной массив без особого взаимодействия с python (только две строки, получение выходного размера и оператор return показывают значительное взаимодействие с 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 в список args в моем файле setup.py;Однако сроки не значительно отличались от ранее.
Из этого замечательного ответа выясняется, что при автоматическом преобразовании между массивными массивами и типизированными представлениями памяти возникают некоторые накладные расходы, которые происходят как при первоначальном вызове функции (так и в операторе 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, что несколько важно, особенно если функция должна бытьзвонил много раз. Конечно, если вы планируете вызывать функцию много раз, вместо этого может быть более эффективным вместо этого просто передавать двумерные массивы, а не просто сохранять значительную нагрузку на вызовы функций Python и амортизировать стоимость numpy array -> typed memoryview преобразования. Кроме того, было бы интересно использовать пустые структурированные массивы, которые в cython могут быть преобразованы в типизированное представление памяти структур, поскольку это может сблизить все данные в кэше и сократить время доступа к памяти.
В качестве заключительного замечания, как было обещано в комментариях ранее, вот версия, использующая prange, которая использует преимущества параллельной обработки. Обратите внимание, что это может использоваться только с типизированными представлениями памяти, поскольку GIL Python должен быть освобожден в цикле prange (и скомпилирован с использованием флага -fopenmp для args и link_args:
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:
Следуя очень полезному дополнительному совету @ max9111 в комментариях, я переключил все объявления float[:] в своем коде на float[::1]. Значение этого заключается в том, что он позволяет хранить данные непрерывно, и Cython не нужно беспокоиться о наличии шага между элементами. Это позволяет использовать SIMD Vectorization, что значительно оптимизирует код. Ниже приведены обновленные тайминги, которые генерируются с помощью следующих команд:
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