Реализация гаммальн в Нумбе
Это может быть довольно трудоемкой реализацией некоторых часто используемых функций, не только для достижения производительности, но и для получения четко определенного уровня точности.
Таким образом, прямой способ будет просто обернуть рабочую реализацию .
В случае gammaln scipy- вызывает C-реализацию этой функции. Поэтому скорость реализации scipy также зависит от флагов компилятора и компилятора, используемых при компиляции зависимостей scipy.
Также не очень удивительно, что результаты производительности для одного значения могут сильно отличаться от результатов больших массивов.
В первом случае преобладают накладные расходы на вызовы (включая преобразования типов, проверку ввода, ...), во втором - производительность
реализация становится все более и более важной.
Улучшение вашей реализации
- Написать явные циклы. В Numba векторизованные операции расширяются до циклов, и после этого Numba пытается присоединиться к циклам. Часто лучше выписать и соединить эти циклы вручную.
- Подумайте о различиях в основных арифметических реализациях. Python всегда проверяет деление на 0 и вызывает исключение в таком случае, что очень дорого. Numba также использует это поведение по умолчанию, но вы также можете переключиться на проверку ошибок Numpy. В этом случае деление на 0 приводит к NaN. Способ обработки NaN и Inf -0 / + 0 в дальнейших вычислениях также зависит от флага быстрой математики.
Код
import numpy as np
from numpy import log
from scipy.special import gammaln
from numba import njit
import numba as nb
@njit(fastmath=True,error_model='numpy')
def gammaln_nr(z):
"""Numerical Recipes 6.1"""
#Don't use global variables.. (They only can be changed if you recompile the function)
coefs = np.array([
57.1562356658629235, -59.5979603554754912,
14.1360979747417471, -0.491913816097620199,
.339946499848118887e-4, .465236289270485756e-4,
-.983744753048795646e-4, .158088703224912494e-3,
-.210264441724104883e-3, .217439618115212643e-3,
-.164318106536763890e-3, .844182239838527433e-4,
-.261908384015814087e-4, .368991826595316234e-5])
out=np.empty(z.shape[0])
for i in range(z.shape[0]):
y = z[i]
tmp = z[i] + 5.24218750000000000
tmp = (z[i] + 0.5) * np.log(tmp) - tmp
ser = 0.999999999999997092
n = coefs.shape[0]
for j in range(n):
y = y + 1.
ser = ser + coefs[j] / y
out[i] = tmp + log(2.5066282746310005 * ser / z[i])
return out
@njit(fastmath=True,error_model='numpy',parallel=True)
def gammaln_nr_p(z):
"""Numerical Recipes 6.1"""
#Don't use global variables.. (They only can be changed if you recompile the function)
coefs = np.array([
57.1562356658629235, -59.5979603554754912,
14.1360979747417471, -0.491913816097620199,
.339946499848118887e-4, .465236289270485756e-4,
-.983744753048795646e-4, .158088703224912494e-3,
-.210264441724104883e-3, .217439618115212643e-3,
-.164318106536763890e-3, .844182239838527433e-4,
-.261908384015814087e-4, .368991826595316234e-5])
out=np.empty(z.shape[0])
for i in nb.prange(z.shape[0]):
y = z[i]
tmp = z[i] + 5.24218750000000000
tmp = (z[i] + 0.5) * np.log(tmp) - tmp
ser = 0.999999999999997092
n = coefs.shape[0]
for j in range(n):
y = y + 1.
ser = ser + coefs[j] / y
out[i] = tmp + log(2.5066282746310005 * ser / z[i])
return out
import matplotlib.pyplot as plt
import seaborn as sns
import time
n_trials = 8
scipy_times = np.zeros(n_trials)
fastats_times = np.zeros(n_trials)
fastats_times_p = np.zeros(n_trials)
for i in range(n_trials):
zs = np.linspace(0.001, 100, 10**i) # evaluate gammaln over this range
# dont take first timing - this is just compilation
start = time.time()
arr_1=gammaln_nr(zs)
end = time.time()
start = time.time()
arr_1=gammaln_nr(zs)
end = time.time()
fastats_times[i] = end - start
start = time.time()
arr_3=gammaln_nr_p(zs)
end = time.time()
fastats_times_p[i] = end - start
start = time.time()
start = time.time()
arr_3=gammaln_nr_p(zs)
end = time.time()
fastats_times_p[i] = end - start
start = time.time()
arr_2=gammaln(zs)
end = time.time()
scipy_times[i] = end - start
print(np.allclose(arr_1,arr_2))
print(np.allclose(arr_1,arr_3))
fig, ax = plt.subplots(figsize=(12,8))
sns.lineplot(np.logspace(0, n_trials-1, n_trials), fastats_times, label="numba");
sns.lineplot(np.logspace(0, n_trials-1, n_trials), fastats_times_p, label="numba_parallel");
sns.lineplot(np.logspace(0, n_trials-1, n_trials), scipy_times, label="scipy");
ax.set(xscale="log");
ax.set_xlabel("Array Size", fontsize=15);
ax.set_ylabel("Execution Time (s)", fontsize=15);
ax.set_title("Execution Time of Log Gamma");
fig.show()