[Python] Threading or multiprocess to run code in parallel seems slower

I've been searching with little success how to solve this problem. The script below is supposed to perform planet simulations. planet1_pars will define 1st planet parameters. set_grids_fakePlanet will create a grid for each of the parameters of a hypothetical planet put into the system. This function will return a generator not a list/array with tons of parameter values. planet2_pars will give me a set of parameters previously created in set_grids_fakePlanet, hence each time I execute planet2_pars it will give me a different set of parameters from the hypothetical planet. ComputeTTV will do some calculations and return what I need each time I execute run_rebound, which is my main function that will call all these mentioned functions above. Whenever I execute run_rebound, I need to give it the hypothetical planet parameter so it run the simulation.

def planet1_pars():
P_p1,m_p1,e_p1 = 0.7920639164 / 365.25, 29.32*3.0027e-6, 0.0 #P[yrs], m[solar],e[fixed]
inc_p1,omega_p1,M_p1 = 77.4041697839 * np.pi/180, 90., 0.
return P_p1,m_p1,e_p1,inc_p1,omega_p1,M_p1

def set_grids_fakePlanet(pars_p1):
P_p1,m_p1,e_p1,inc_p1,omega_p1,M_p1 = [*pars_p1]

#min max periods in which to put a planet
Pmin = P_p1 * 2.02 # Pmin ~ 0.0196/365.25 #shortest period found so far in exoplanet.eu
Pmax = P_p1 * 2.05
#set grids
P_grid = np.arange(Pmin, Pmax, P_p1 * 0.005)
m_p2_grid = np.arange(.5, 320, 1) * 3.0027e-6 # every 1 Earth mass to 1 Jupiter
e_grid = [0.0]#np.linspace(0,0.1, 10) # e=1 may cause code to blow up
inc_grid = [inc_p1]#np.linspace(60,90, 5)
omega_grid = [0.0]#np.linspace(0,360, 5)
M_grid = [0.0]#np.linspace(0,360, 5)

#store grid vals, each column is a parameter, last column TTV amplitude
#[n,m] n is max_size(P,e,inc,omega,M) ** m. m is # of orbital parameters + 1 ttv amp
size = len(P_grid) * len(m_p2_grid) * len(e_grid) * len(inc_grid) * len(omega_grid) * len(M_grid)
results = np.zeros([size,6+1]) * np.nan

peiom_grid = ((x,k,y,w,j,z) for x in P_grid for k in m_p2_grid for y in e_grid for w in inc_grid
for j in omega_grid for z in M_grid)
return peiom_grid

def planet2_pars():
for pars_p2 in peiom_grid:
return pars_p2

#2nd planet
# m_p2, P_p2, e_p2, inc_p2, omega_p2, M_p2 = system_parameters(n*m_p1, P_p2,e_p2,inc_p2,omega_p2,M_p2)


def computeTTVs(sim, P_p1, P_p2):
N=34
transittimes = np.zeros(int(N))
p = sim.particles
i = 0
while i<N:
y_old = p[1].y - p[0].y # (Thanks to David Martin for pointing out a bug in this line!)
t_old = sim.t
if P_p1 > P_p2:
sim.integrate(sim.t+ (P_p2 * 0.05)) # check for transits every 0.5 time units. Note that 0.5 is shorter than one orbit
else:
sim.integrate(sim.t+ (P_p1 * 0.05)) #5% of period ~ 1h which is shorter than Tdur=2h

t_new = sim.t
if y_old*(p[1].y-p[0].y)<0. and p[1].x-p[0].x>0.: # sign changed (y_old*y<0), planet in front of star (x>0)
while t_new-t_old>1e-7: # bisect until prec of 1e-5 reached
if y_old*(p[1].y-p[0].y)<0.:
t_new = sim.t
else:
t_old = sim.t
sim.integrate( (t_new+t_old)/2.)
transittimes = sim.t
i += 1
sim.integrate(sim.t+ P_p1 * 0.01) # integrate 0.05 to be past the transit

A = np.vstack([np.ones(N), range(N)]).T
c, m = np.linalg.lstsq(A, transittimes, rcond=-1)[0] # fits a linear model and get period m and t0 c
comp_t0s = c + m*np.array(range(N))

OC = transittimes-comp_t0s # in years
OC *= 365.25*24*60
amp = rms(OC)
# amp = np.diff([np.min(OC), np.max(OC)])[0]
return amp #in minutes

def run_rebound(pars_p2):
ms=1.02 #solar unit
P_p1,m_p1,e_p1,inc_p1,omega_p1,M_p1 = planet1_pars()
P_p2,m_p2,e_p2,inc_p2,omega_p2,M_p2 = [*pars_p2] #fake planet
#start simulation
sim = rebound.Simulation()
sim.G = 39.478 #AU^3 yr^-2 Ms^-1
sim.add(m=ms)
sim.add(m=m_p1, P=P_p1, e=e_p1, inc=inc_p1, omega=omega_p1, M=M_p1)
sim.add(m=m_p2, P=P_p2, e=e_p2, inc=inc_p2, omega=omega_p2, M=M_p2)
#put outcomes in a list
results = [P_p2,m_p2,e_p2,inc_p2*(180/np.pi),omega_p2,M_p2, computeTTVs(sim, P_p1, P_p2)]
return results


Question: I tried to make it parallel using the threading library as in:

peiom_grid = set_grids_fakePlanet(planet1_pars()) #make the fake planet grid as a generator variable

import threading

start = time.time()
for pars in peiom_grid:
t1 = threading.Thread(target=run_rebound, args=(pars,))

t1.start()
t1.join()


end = time.time()

print((end-start) /60, 'min')


In this manner, I see the 8 CPU I got is being used but at a rate which is less than 50%. And it takes ~ 1.2 min to run (the grids are small because I am testing, but ideally the grids should be lager so it may take days to run).

I also tried MultiProcessing

from multiprocessing import Process
start = time.time()

if __name__ == '__main__':
for pars in peiom_grid:

p = Process(target=run_rebound, args=(pars,))
p.start()
p.join()

end = time.time()
print((end-start) /60, 'min')


it takes ~ 1.7min

and without any parallelization

start = time.time()

for pars in peiom_grid:
run_rebound(pars)

end = time.time()
print((end-start)/60, 'min')


it takes ~ 1.34 min

I think I am not doing any parallelization because the difference between the runs above with/without parallelization isn't significant. I cannot find where the issue is. I followed a few examples and check several examples on stack overflow but nothing... Hope you guys can give me some feedback.

Continue reading...