Распараллеливание/оптимизация интегралов в Python ⇐ Python

[Anonymous]

Мне нужно вычислить большой набор НЕЗАВИСИМЫХ интегралов, и сейчас я написал для этого следующий код на Python:

Код: Выделить всё

# We define the integrand
def integrand(tau,k,t,z,config, epsilon=1e-7):
u = np.sqrt(np.maximum(epsilon,tau**2 - z**2))
return np.sin(config.omega * (t - tau))*k/2, np.sin(config.omega * (t - tau)) * j1(k * u) / u

NON-VECTORISED
partial_integral = np.zeros((len(n_values),len(t_values),len(z_values)))

for n in range(1, len(n_values)): # We skip the n=0 case as it is trivially = 0
for j in range(1, len(z_values)): # We skip the z=0 case as it is trivially = 0
for i in range(1, len(t_values)): # We skip the t=0 case as it is trivially = 0
partial_integral[n,i,j],_ = quad(integrand, x_min[n,i,j], x_max[n,i,j], args=(k_n_values[n],t_values[i],z_values[j],config), limit=np.inf) # We use quad

Теперь все len(n_values),len(t_values),len(z_values) являются большими числами, поэтому я хотел бы максимально ускорить код. Есть предложения?

Помимо экспериментов с различными библиотеками интеграции (среди которых Scipy.quad кажется лучшей), я думал о векторизации кода:

Код: Выделить всё

# VECTORISED
def compute_integral(n,i,j):
quad(integrand, x_min[n,i,j], x_max[n,i,j], args=(k_n_values[n],t_values[i],z_values[j],config), limit=np.inf) # We use quad

# Use np.meshgrid to create the index grids for n, i, j (starting from 1 to avoid 0-index)
n_grid, i_grid, j_grid = np.meshgrid(np.arange(0, len(n_values)), np.arange(0, len(t_values)), np.arange(0, len(z_values)), indexing='ij')

# Flatten the grids to vectorize the loop over n, i, j
indices = np.vstack([n_grid.ravel(), i_grid.ravel(), j_grid.ravel()]).T

# Vectorize the integral computation using np.vectorize
vectorized_integral = np.vectorize(lambda n, i, j: compute_integral(n, i, j))

# Apply the vectorized function to all combinations of (n, i, j)
partial_integral = np.empty((len(n_values),len(t_values),len(z_values)))
partial_integral[tuple(indices.T)] = vectorized_integral(*indices.T)

Но не факт, что это меня сильно выиграет...
Я также пытался использовать Numba (с numba-scipy для j1 ) для JIT подынтегральной функции в надежде немного повысить производительность, и это действительно помогло мне получить производительность x5!
PS: В качестве дополнения я сейчас запускаю этот скрипт на моем компьютере, но я, вероятно, смогу запустить его на кластер.

Подробнее здесь: https://stackoverflow.com/questions/793 ... -in-python