Surface Density Profile 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])
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.