Surface Density Profile

Question

I am trying to compute the surface density profile, given the spherical density profile in 3D for different parameters in order to interpolate and have it as a function of them to use later on for some gravitational lensing modelling. However, looping through all parameters and evaluating the integral over the 2D Radius array R is extremely time consuming (many hours for complicated profiles).

The code that follows describe the procedure I am following in order to do that. Does anyone have a better suggestion to improve the performance of the following code in Python?

def density_profile(r, par_a, par_b, par_c):
    ...
    return ...

def surface_denstiy_profile(R, density_profile, *args):

    surface_density_array = []

    for R_i in R:

        surface_density_array.append( quad(lambda r : 2.0 * density_profile(r , *args, ) * r / np.sqrt(r**2. - R_i**2.), R_i, np.inf, epsabs=1.49e-03, epsrel=1.49e-03)[0] ) )

    return surface_density_array


surface_density = np.zeros((len(par_a_array), len(par_b_array), len(par_b_array)), dtype = object)

for par_a in range(0, len(par_a_array)) :

    for par_b in range(0, len(par_b_array)):

        for par_c in range(0, len(par_c_array)):

            surface_density[par_a, par_b, par_c] = surface_denstiy_profile(np.logspace(-2., 2., 100), density_profile, par_a_array[par_a], par_b_array[par_b], par_c_array[par_c])

Answer 1

Firstly try and do as few computation as possible - e.g. the np.logspace call is essentially constant, no reason to call it over and over for every single element.

Also the dtype of object for surface_density makes using NumPy somewhat useless for that array, try instead one dimension more and the right datatype / dimensions of the array to actually save some space / make things faster.

With NumPy you really want to avoid these loops and try to find some operations that work on the whole array instead of per-element. If there's no way to do that easily, consider some of the apply functions to at least make it look better and then look into processing parts of the array in multiple threads / processes to speed things up as it looks like to be a highly parallel computation. Also Cython, but it's not clear to me if it would help that much here, even though trying it is not very complicated, so you could just do that and see.

Also avoid slower language features like the variable arguments in often called functions, e.g. here pass a fixed number of arguments instead of doing the *args thing.

---

Unrelated to NumPy and won't help here, but related to general Python, just in case for your next scripts:

• There's also a typo in surface_denstiy_profile, it should be density.
• Don't use the range(..., len(...)) pattern, better would be to have enumerate since you get both the index and value delivered directly.