Вы можете использовать Numba или низкоуровневую
Почти ваш пример
Я просто передаю функцию непосредственно в scipy.integrate.dblquad вместо вашего метода, использующего лямбды для генерации функций.
import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import time
q = np.linspace(0.03, 1.0, 1000)
start = time.time()
def f(t, z, q):
return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * (1 / (np.sqrt(2 * np.pi) * 2)) * np.exp(
-0.5 * ((z - 40) / 2) ** 2)
def lower_inner(z):
return 10.
def upper_inner(z):
return 60.
y = np.empty(len(q))
for n in range(len(q)):
y[n] = integrate.dblquad(f, 0, 50, lower_inner, upper_inner,args=(q[n],))[0]
end = time.time()
print(end - start)
#143.73969149589539
Это уже чуть-чуть быстрее (143 против 151 с), но единственное использование - простой пример оптимизации.
Простая компиляция функции, использующие Numba
Чтобы запустить это приложение, вам необходимо дополнительно Numba и numba-scipy . Назначение numba-scipy - предоставить упакованные функции из scipy.special.
import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import time
import numba as nb
q = np.linspace(0.03, 1.0, 1000)
start = time.time()
#error_model="numpy" -> Don't check for division by zero
@nb.njit(error_model="numpy",fastmath=True)
def f(t, z, q):
return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * (1 / (np.sqrt(2 * np.pi) * 2)) * np.exp(
-0.5 * ((z - 40) / 2) ** 2)
def lower_inner(z):
return 10.
def upper_inner(z):
return 60.
y = np.empty(len(q))
for n in range(len(q)):
y[n] = integrate.dblquad(f, 0, 50, lower_inner, upper_inner,args=(q[n],))[0]
end = time.time()
print(end - start)
#8.636585235595703
Использование низкоуровневого вызова
Функции scipy.integrate также предоставляют возможность передать функцию C -callback вместо функции Python. Эти функции могут быть написаны, например, в C, Cython или Numba, которые я использую в этом примере. Основным преимуществом является то, что при вызове функции не требуется Python взаимодействие с интерпретатором.
Отличный ответ @Jacques Gaudin показывает простой способ сделать это, включая дополнительные аргументы.
import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import time
import numba as nb
from numba import cfunc
from numba.types import intc, CPointer, float64
from scipy import LowLevelCallable
q = np.linspace(0.03, 1.0, 1000)
start = time.time()
def jit_integrand_function(integrand_function):
jitted_function = nb.njit(integrand_function, nopython=True)
#error_model="numpy" -> Don't check for division by zero
@cfunc(float64(intc, CPointer(float64)),error_model="numpy",fastmath=True)
def wrapped(n, xx):
ar = nb.carray(xx, n)
return jitted_function(ar[0], ar[1], ar[2])
return LowLevelCallable(wrapped.ctypes)
@jit_integrand_function
def f(t, z, q):
return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * (1 / (np.sqrt(2 * np.pi) * 2)) * np.exp(
-0.5 * ((z - 40) / 2) ** 2)
def lower_inner(z):
return 10.
def upper_inner(z):
return 60.
y = np.empty(len(q))
for n in range(len(q)):
y[n] = integrate.dblquad(f, 0, 50, lower_inner, upper_inner,args=(q[n],))[0]
end = time.time()
print(end - start)
#3.2645838260650635