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

• Could you post more context, perhaps something that can be run as a minimal example? Also what's quad and how does density_profile look like? Also is this with Python 2.7 or Python 3.x? – ferada Sep 8 '16 at 17:34
• quad is part of the scipy.integrate package to numerically integrate a function. An exampIe for the density profile would be rho = (1.0 / r**2) times a function of the parameters. I am using Python 2.7 – Sketos Sep 8 '16 at 18:00

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.
• I have tried your suggestions before but they don't really improve the speed (they are more of "how to write a better code in python"). The aim i think is to find a way (if there is one) to avoid computing all these integrals, basically improve the surface_denstiy_profile function somehow. Indeed one solution whould be to make the code parallel. – Sketos Sep 8 '16 at 17:56
• If you need the value of those integrals and there's no mathematical transformation that makes it easier then yes. However you can still make it fast by getting as many computations as possible out of the quad call, e.g. if I remember correctly the multiplication by a constant can be moved to the outside, maybe there's other integral transformations which you can use to split the integration into constant and variable parts. – ferada Sep 8 '16 at 18:30