Vectorizing Python code to improve performance

I am writing scientific code in python to calculate the energy of a system. Here is my function: cte1, cte2, cte3, cte4 - these are the constants previously calculated; pii - np.pi (calculated in advance since it slows down the loop otherwise). I calculate the 3 components of the total energy, then add them up.

def calc_energy(diam): 
    Energy1 = cte2*((pii*diam**2/4)*t)
    Energy2 = cte4*(pii*diam)*t
    d=diam/t
    u=np.sqrt((d)**2/(1+d**2))
    cc= u**2
    E = sp.special.ellipe(cc) 
    K = sp.special.ellipk(cc) 
    Id=cte3*d*(d**2+(1-d**2)*E/u-K/u)
    Energy3 = cte*t**3*Id
    total_energy = Energy1+Energy2+Energy3
    return (total_energy,Energy1)

      

My first idea was to just loop over all the diameter values:

start_diam, stop_diam, step_diam = 1e-10, 500e-6, 1e-9 #Diametre
diametres = np.arange(start_diam,stop_diam,step_diam)

for d in diametres:  
    res1,res2 = calc_energy(d)
    totalEnergy.append(res1)
    Energy1.append(res2)

      

In an attempt to speed up calculations, I decided to use numpy for vectorization as shown below:

diams = diametres.reshape(-1,1) #If not reshaped, calculations won't run
r1 = np.apply_along_axis(calc_energy,1,diams)

      

However, the "vectorized" solution does not work as expected. When syncing, I get 5 seconds for the first solution and 18 seconds for the second.

I think I am doing something wrong, but I cannot figure out what.

+3


source to share


1 answer


With your current approach, you are applying a Python function to each element of your array, which incurs additional overhead. Instead, you can pass the entire array to your function and get an array of responses. Your existing function works fine without any changes.



import numpy as np
from scipy import special
cte = 2
cte1 = 2
cte2 = 2
cte3 = 2
cte4 = 2
pii = np.pi

t = 2

def calc_energy(diam): 
    Energy1 = cte2*((pii*diam**2/4)*t)
    Energy2 = cte4*(pii*diam)*t
    d=diam/t
    u=np.sqrt((d)**2/(1+d**2))
    cc= u**2
    E = special.ellipe(cc) 
    K = special.ellipk(cc) 
    Id=cte3*d*(d**2+(1-d**2)*E/u-K/u)
    Energy3 = cte*t**3*Id
    total_energy = Energy1+Energy2+Energy3
    return (total_energy,Energy1)

start_diam, stop_diam, step_diam = 1e-10, 500e-6, 1e-9 #Diametre
diametres = np.arange(start_diam,stop_diam,step_diam)

a = calc_energy(diametres) # Pass the whole array

      

+2


source







All Articles