Running Multiple Simulations in Parallel¶
This tutorial will go over the basics of running multiple fiber simulations in parallel. We are often curious about not just a single fiber but an array of fibers with various model types, geometric parameters, and/or electrical parameters. Or, we might want to run the same fiber through multiple simulations. We can leverage parallelism to run multiple simulations simultaneously, each on a separate processor core.
Create a function to parallelize¶
First, let’s create a function that we can call in parallel. The function should create a Fiber instance and solve for its activation threshold. We will use the fiber model and stimulation parameters from the simulation tutorial. Instead of a single fiber diameter, we will create a function which takes a fiber diameter as an argument, then returns the activation threshold of the fiber.
def create_and_run_sim(diam=5.7, temp=37):
"""Create a fiber and determine activate threshold.
:param diam: diameter of fiber (um).
:param temp: fiber temperature (C)
:return: returns activation threshold (mA)
"""
from pyfibers import build_fiber, FiberModel, ScaledStim
from scipy.interpolate import interp1d
# Create fiber object
fiber = build_fiber(
FiberModel.MRG_INTERPOLATION, diameter=diam, n_sections=265, temperature=temp
)
fiber.potentials = fiber.point_source_potentials(0, 250, fiber.length / 2, 1, 10)
# Setup for simulation
time_step = 0.001
time_stop = 20
waveform = interp1d([0, 0.2, 0.4, time_stop], [1, -1, 0, 0], kind="previous")
# Create stimulation object
stimulation = ScaledStim(waveform=waveform, dt=time_step, tstop=time_stop)
amp, ap = stimulation.find_threshold(fiber, silent=True) # Find threshold
return amp
Parallelization with multiprocess¶
The multiprocess package provides a way to create and manage multiple processes in Python, similar to how the
threading module handles threads. The Pool object creates a pool of processes which can be used to parallelize our
fiber jobs. See the multiprocess documentation.
Determine available cpus¶
Before submitting any jobs, first use the multiprocess package to see the number of cpus available on your machine.
import multiprocess
cpus = multiprocess.cpu_count() - 1
print(cpus)
3
Parallelize fiber jobs for a list of fibers¶
Now, create an instance of the multiprocess.Pool class. Finally, we can use the Pool.starmap() method in the Pool class to submit our jobs to the process pool. The Pool.starmap() method allows us to pass in a function with multiple arguments to simultaneously submit jobs. For this tutorial, we will demonstrate submitting local parallel jobs to find the activation threshold for a list of fibers, each with a unique diameter.
Note
Note, you must place the Pool.starmap() call inside of an if __name__ == "__main__": statement, as shown below, otherwise your Python code will generate an infinite loop. Besides function definitions, all other functionality you use should be under this statement as well.
from multiprocess import Pool
if __name__ == "__main__":
fiber_diams = [2.0, 5.7, 8.7, 11.5, 14.0]
temp = 37
params = [(diam, temp) for diam in fiber_diams]
with Pool(cpus) as p:
results = p.starmap(create_and_run_sim, params)
Running:
-1
Running:
-1
Running:
-1
N aps: 1, time 0.346
Running:
-0.01
N aps: 1, time 0.338
Running:
-0.01
N aps: 0, time None
Running:
-0.01
N aps: 0, time None
Running:
-0.505
N aps: 0, time None
Running:
-0.505
N aps: 0, time None
Running:
-0.7525
N aps: 0, time None
Running:
-0.7525
N aps: 1, time 0.394
Running:
-0.62875
N aps: 0, time None
Running:
-1.1
N aps: 0, time None
Running:
-0.87625
N aps: 1, time 0.377
Running:
-0.814375
N aps: 0, time None
Running:
-0.690625
N aps: 0, time None
Running:
-0.721562
N aps: 1, time 0.45
Running:
-0.706094
N aps: 0, time None
Running:
-0.845312
N aps: 1, time 0.415
Running:
-0.829844
N aps: 0, time None
Running:
-0.713828
N aps: 0, time None
Running:
-1.21
N aps: 1, time 0.569
Running:
-0.709961
N aps: 0, time None
Running:
-0.837578
N aps: 1, time 0.462
Running:
-0.833711
N aps: 1, time 0.55
Running:
-0.833711
N aps: 0, time None
Running:
-0.711895
N aps: 0, time None
Running:
-0.713828
N aps: 0, time None
Running:
-1.331
N aps: 1, time 0.55
Running:
-1
N aps: 1, time 0.318
Running:
-0.01
N aps: 0, time None
Running:
-0.505
N aps: 1, time 0.569
Running:
-1
N aps: 0, time None
Running:
-0.7525
N aps: 0, time None
Running:
-1.4641
N aps: 1, time 0.369
Running:
-0.62875
N aps: 1, time 0.297
Running:
-0.01
N aps: 1, time 0.401
Running:
-1.39755
N aps: 1, time 0.5
Running:
-1.364275
N aps: 0, time None
Running:
-0.690625
N aps: 1, time 0.432
Running:
-0.659688
N aps: 0, time None
Running:
-1.380913
N aps: 0, time None
Running:
-0.505
N aps: 0, time None
Running:
-0.675156
N aps: 0, time None
Running:
-1.389231
N aps: 0, time None
Running:
-0.7525
N aps: 1, time 0.349
Running:
-0.62875
N aps: 0, time None
Running:
-0.682891
N aps: 0, time None
Running:
-1.393391
N aps: 0, time None
Running:
-0.690625
N aps: 1, time 0.517
Running:
-0.679023
N aps: 1, time 0.383
Running:
-0.659688
N aps: 1, time 0.838
Running:
-1.393391
N aps: 0, time None
Running:
-0.675156
N aps: 0, time None
Running:
-0.680957
N aps: 1, time 0.584
Running:
-0.680957
N aps: 1, time 0.454
Running:
-0.667422
N aps: 0, time None
Running:
-0.671289
N aps: 1, time 0.838
N aps: 1, time 0.522
Running:
-0.669355
N aps: 1, time 0.625
Running:
-0.669355
N aps: 1, time 0.584
N aps: 1, time 0.625
Let’s plot the activation threshold vs. fiber diameter to see if a relationship between the two exists.
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
sns.set(font_scale=1.5, style='whitegrid', palette='colorblind')
plt.figure()
plt.plot(fiber_diams, -np.array(results), marker='o')
plt.xlabel('Diameter (microns)')
plt.ylabel('Activation threshold (mA)')
plt.show()
Parallelization is not just limited to running multiple fiber diameters. You could also test the same fiber with different stimulation parameters, or different numbers of sections. Let’s do another example, except this time, let’s vary the number of sections. Again, let’s visualize the data to see if a relationship exists between fiber length and activation threshold.
if __name__ == "__main__":
diam = 5.7
temps = [20, 26, 32, 37]
params = [(diam, temp) for temp in temps]
with Pool(cpus) as p:
results = p.starmap(create_and_run_sim, params)
Running:
-1
Running:
Running:
-1
-1
N aps: 1, time 0.414
Running:
-0.01
N aps: 0, time None
Running:
-0.01
N aps: 0, time None
Running:
-0.01
N aps: 0, time None
Running:
-0.505
N aps: 0, time None
Running:
-1.1
N aps: 0, time None
Running:
-0.7525
N aps: 0, time None
Running:
-0.87625
N aps: 0, time None
Running:
-1.21
N aps: 1, time 0.476
Running:
-1.155
N aps: 1, time 0.671
Running:
-1.1275
N aps: 0, time None
Running:
-1.14125
N aps: 0, time None
Running:
-0.938125
N aps: 0, time None
Running:
-1.1
N aps: 0, time None
Running:
-1.148125
N aps: 1, time 0.865
Running:
-1.144688
N aps: 0, time None
Running:
-0.969062
N aps: 1, time 0.449
Running:
-0.953594
N aps: 0, time None
Running:
-1.148125
N aps: 1, time 0.532
Running:
-0.945859
N aps: 0, time None
Running:
-0.949727
N aps: 0, time None
Running:
-1.21
N aps: 1, time 0.611
Running:
-0.949727
N aps: 1, time 0.865
Running:
-1
N aps: 1, time 0.346
Running:
-0.01
N aps: 1, time 0.611
N aps: 0, time None
Running:
-1.331
N aps: 0, time None
Running:
-0.505
N aps: 0, time None
Running:
-0.7525
N aps: 0, time None
Running:
-1.4641
N aps: 1, time 0.611
Running:
-1.39755
N aps: 0, time None
Running:
-0.87625
N aps: 1, time 0.377
Running:
-0.814375
N aps: 0, time None
Running:
-1.430825
N aps: 1, time 0.808
Running:
-1.414188
N aps: 1, time 1.276
Running:
-1.405869
N aps: 0, time None
Running:
-0.845312
N aps: 1, time 0.415
Running:
-0.829844
N aps: 0, time None
Running:
-1.410028
N aps: 0, time None
Running:
-0.837578
N aps: 1, time 0.462
Running:
-0.833711
N aps: 1, time 0.55
Running:
-0.833711
N aps: 0, time None
Running:
-1.414188
N aps: 1, time 0.55
N aps: 1, time 1.276
plt.figure()
plt.plot(temps, -np.array(results), marker='o')
plt.xlabel('Temperature (C)')
plt.ylabel('Activation threshold (mA)')
plt.show()
