Parallelisierung/Optimierung von Integralen in Python ⇐ Python

[Guest]

Ich muss eine große Menge UNABHÄNGIGER Integrale berechnen und habe derzeit den folgenden Code dafür in Python geschrieben:

Code: Select all

# 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

Jetzt sind alle len(n_values),len(t_values),len(z_values) große Zahlen, daher möchte ich den Code so weit wie möglich beschleunigen. Irgendwelche Vorschläge?

Neben dem Experimentieren mit verschiedenen Integrationsbibliotheken (von denen Scipy.quad die beste zu sein scheint) habe ich darüber nachgedacht, den Code zu vektorisieren:

Code: Select all

# 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)

Aber es ist nicht klar, ob mir das viel bringt...
Ich habe auch versucht, Numba zu verwenden (mit numba-scipy für die j1 ) zum JIT des Integranden in der Hoffnung, ein wenig Leistung zu gewinnen, und das hat mir tatsächlich geholfen, eine x5-Leistung zu erzielen!
PS: Als Extra führe ich derzeit dieses Skript aus in meinem PC, aber ich könnte es wahrscheinlich auf einem ausführen Cluster.