I'm trying to numerically integrate a 2D complex-valued function multiple times. The problem with the current setup which I have which works fine but it is understandably slow.
The slow part of the below code is ssht.create_ylm(theta, phi, self.L) which basically finds the spherical harmonics up to a band-limit L for the angles in question (can be arrays rather than floats). At the moment I'm calling self.f twice for the real and imaginary parts but this, of course, means that I'm being very inefficient. However, I'm not sure how to about saving intermediate values.
I looked into quadpy which seems to be able to handle complex-valued functions but not sure how to go about this.
The below example is a sub-set of a class:
import numpy as np
from scipy import integrate
import pyssht as ssht # user-defined cython code
def f(self, theta, phi, i, j):
ylm = ssht.create_ylm(theta, phi, self.L)
ylm = ylm.reshape(ylm.size)
f = ylm[i] * np.conj(ylm[j]) * np.sin(theta)
return f
def real_func(self, theta, phi, i, j):
return self.f(theta, phi, i, j).real
def imag_func(self, theta, phi, i, j):
return self.f(theta, phi, i, j).imag
def integral(self, f, i, j):
F = integrate.dblquad(
f,
self.phi_min,
self.phi_max,
lambda t: self.theta_min,
lambda t: self.theta_max,
args=(i, j),
)[0]
return F
def D_integral(self, i, j):
F_real = self.integral(self.real_func, i, j)
F_imag = self.integral(self.imag_func, i, j)
return F_real + 1j * F_imag
def D_matrix(self):
D = np.zeros((self.side.size, self.side.size), dtype=complex)
for i in range(self.side.size):
for j in range(i, self.side.size):
integral = self.D_integral(self.side[i], self.side[j])
D[i][j] = integral
if i != j:
D[j][i] = np.conj(integral)
return D
Ideally, I'd like to be able to do this in a more clever way, and if not at least in a way which avoids the need to call ssht.create_ylm twice per one integration step.
