h.frecord_init() and Python Multiprocessing

[Robert Claar]

Hello,

I ran into an odd issue with my model in NEURON+Python that I think I have fixed but would appreciate further clarity to check that I have addressed it correctly.

I am simulating a large network of endocrine cells with lots of variables and need to run long simulations (15-30 minutes) due to slowly changing hormone concentrations that we are interested in studying. For this reason I have followed the advice here on the forum to break the simulations into chunks and have the following workflow:

Code: Select all

# Create a dictionary with recorded variable names as keys and values as h.Vectors
rec_dict = create_rec_dict(section_list)

# Chunk end times is a list of end times for each chunk
# For example if I was running a 5000 ms simulation in chunks of 1000 ms
# it would be the list [1000, 2000, 3000, 4000, 5000]
for chunk_end_time in chunk_end_times:
	
	# Run simulation until the end of the next time chunk
	h.continuerun(chunk_end_time)
	
	# Save chunk data... Omitting code for brevity 
	
	# Resize recording vectors so that max amount of memory consumed is close to
	# the amount required to store one time chunks worth of data
	h.frecord_init()

This works great, but my simulations were still very long and I wanted to increase the efficiency wherever I could. To do this I had the idea to use python multiprocessing to have my simulations run on my main process, have worker processes that process each of my recording h.Vectors, and then have those worker processes pass the processed data to another process through a multiprocessing.Manager.Queue that handles writing the data to disk. This way my simulations could continue running while data is being processed and written. This looked something like

Code: Select all

import multiprocessing as mp

# Create a dictionary with recorded variable names as keys and values as h.Vectors
rec_dict = create_rec_dict(section_list)

# Create multiprocessing manager to manage data queue
with mp.Manager() as manager:
                
	# Create data queue
	data_queue = manager.Queue()
	
	# Create writer process that will get processed data arrays from the data queue and write them to disk
	data_writer = mp.Process(target=data_writer_func, args=(data_writer_args))
	
	# Create multiprocessing pool to handle data_processor worker processes
	with mp.Pool() as pool:

		# Chunk end times is a list of end times for each chunk
		# For example if I was running a 5000 ms simulation in chunks of 1000 ms
		# it would be the list [1000, 2000, 3000, 4000, 5000]
		for chunk_end_time in chunk_end_times:
	
			# Run simulation until the end of the next time chunk
			h.continuerun(chunk_end_time)
			
			# Convert time vector to numpy array
			time_array = rec_dict["Time"].as_numpy()
			
			# Save chunk data
			for var, vec in rec_dict.items():
				
				if var == "Time":
					continue
					
				# Convert vector to numpy array
				var_array = vec.as_numpy()
				
				# Create worker process to process this data before passing it the data writer
				pool.apply_async(target=data_processor, args=(var, var_array, time_array))
	
			# Resize recording vectors so that max amount of memory consumed is close to
			# the amount required to store one time chunks worth of data
			h.frecord_init()

My simulations became much faster but I started to realized when visualizing my data that many of the variables appeared constant for multiple time chunks when I know they should be changing. Skipping past many tests I ran to figure out the issue, I believe that before all of the worker processes start h.frecord_init() is called which resizes them. Since numpy arrays created from using the as_numpy method on recording vectors map to the same place in memory as the recording vectors, any changes to the recording vectors will also be present in the numpy arrays. Therefore, the data passed to any worker process that hasn't started yet when h.frecord_init() is called has their data array reduced to the last value and this is why the data appears constant when I plot it.

I tested this by running a very small simulation (takes around 1 second to complete) with 2 time chunks and added the statement

Code: Select all

import time
# Sleep for 10 seconds
time.sleep(10)

right before my call to h.frecord_init(). This corrected the issue. I also tested this by creating deepcopies (copies of my recording vectors stored in a different place in memory than the original) of my numpy arrays and then passed those copies to my worker processes and this also corrected the issue.

Does all this make sense? And if so, does anyone have any recommendations for best practices on how to have multiple processes handle running simulations, processing, and writing data?

Any advice or thoughts are appreciated!

[ted]

For embarrassingly parallel problems like this I have always used NEURON's bulletin-board style parallelization with MPI, so I can't give you any advice about python multiprocessing. That said, here's a suggestion that could reduce data storage requirements a bit: record the state variables that have very slow dynamics at long intervals instead of at every new time step. You probably already thought of this. Others who may read this might wonder "fine--given a particular signal, what should the sampling interval be?" Here's a useful reference: Acquiring an Analog Signal: Bandwidth, Nyquist Sampling Theorem, and Aliasing https://www.ni.com/en/shop/data-acquisi ... orem-.html

[ramcdougal]

I think your hypothesis is correct: apply_async does a task in parallel with whatever happens after that, which in your case is to reset the recording vectors. as_numpy is a O(1) method because all it does is create a new data structure referring to the same memory locations as the Vector.

Effectively the same thing as your deep copy solution: You could try making a copy with: vec_array = list(vec)

That should work, but it is O(n) in both memory and time.

[Robert Claar]

Thank you both for your responses.

Regarding Teds response, I have slowly changing hormone concentrations with periods of about 5 minutes that are dependent upon much faster changing electrical, calcium, and metabolic components that I need to record. Is there a way to specify different sampling frequencies for different variables? I am using CVODE variable time step and then resampling my recording vectors to a fixed time step.

[ted]

Time to (re-)read the documentation of the Vector class's record() method at nrn.readthedocs.io. I bet it has been a long time since you last read it. It's hard to retain all that stuff after a single reading, especially the parts that don't pertain to the problem at hand (which is usually simple data capture after each fadvance()). Also, there's a good chance that the documentation has been revised/updated/clarified since the last time you read it.

Note that sampling at fixed intervals should work even if you're using adaptive integration. If a sample time does not coincide with a solution time, the sampled value is supposed to be determined by linear interpolation between the two adjacent solution times. You'll want to verify that by testing on a toy problem. A quick, convincing visual test might be to plot the detailed time course of the solution as a line and decorate the line with markers at each solution time (short vertical line could work), then on the same graph plot the recorded samples with markers of a different color (could be dots of a different color).

[Robert Claar]

I actually created a post about this titled "Saving fixed timestep data using cvode variable timestep integration" which can be found here viewtopic.php?t=4668 and followed your (Ted) suggested solution to run my simulation using variable timestep integration and then resample my data using linear interpolation.

This post was earlier this year, but after looking at the documentation again I ran the following test and see that I can just supply a time h.Vector or a desired Dt as an argument for record and accomplish the same thing.

Code: Select all

from neuron import h
h.load_file("stdrun.hoc")

# Create test section
sec = h.Section(name="sec")

# Set simulation run time in ms
run_time = 500.0

# Set desired final dt (i.e., sampling rate) in ms
dt = 0.1

# Create time vector with desired sampling rate
time_vec = h.Vector().indgen(0, run_time, dt)

# Insert hh mechanism
sec.insert(h.hh)

# Enable variable timestep integration
cvode = h.CVode()
cvode.active(True)

# Create three time vectors, one using our prespecified time vector above, 
# one using dt, and another with variable timestep
t_time_vec = h.Vector().record(h._ref_t, time_vec)
t_dt = h.Vector().record(h._ref_t, dt)
t = h.Vector().record(h._ref_t)            

# Create three corresponding Vm recording vectors
v_time_vec = h.Vector().record(sec(0.5)._ref_v, time_vec)
v_dt = h.Vector().record(sec(0.5)._ref_v, dt)
v = h.Vector().record(sec(0.5)._ref_v)

# Run simulation until provided run time
h.finitialize()
h.continuerun(run_time)

After I printed out the sizes of each different vector and each of the different time vector sizes were equal to their corresponding Vm vector size. I also visualized their timeseries and the data looks the same. The Dt and time_vec approaches generated the same data besides some differing values due to floating point accuracy I'm guessing (i.e.,

Code: Select all

print(t_time_vec.to_python())

outputs
...
2.2,
2.3000000000000003,
2.4000000000000004,
2.5,
2.6,
2.7,
2.8000000000000003,
2.9000000000000004,
3.0,
3.1,
...

and

Code: Select all

print(t_dt.to_python())

outputs
...
2.2000000000000006,
2.3000000000000007,
2.400000000000001,
2.500000000000001,
2.600000000000001,
2.700000000000001,
2.800000000000001,
2.9000000000000012,
3.0000000000000013,
3.1000000000000014,
...

[ted]

Yes, that's the imprecision of fixed precision floating point arithmetic. Then I remember that few parameters of any biological system are known to more than two significant figures and shudder

after looking at the documentation again I ran the following test and see that I can just supply a time h.Vector or a desired Dt as an argument for record and accomplish the same thing.

Absolutely right. Funny thing about documentation: it sometimes has implications that aren't apparent on the first reading. Or second. Or . . . time to stop before I run out of toes . . .

[Robert Claar]

Another follow-up here as providing a Dt value inside the Vector().record() method has now broken another part of my code base.

In a previous post (viewtopic.php?t=4605) I got some advice on how to ensure that, when using CVode and running my simulations in chunks, h.t is equal to my prespecified chunk end time at the end of each time chunk. This is desired because otherwise I will have random time chunks where my recording vectors are one element longer or shorter which raises errors. The methods that were recommended to me are combined and summarized as

Code: Select all

from neuron import h
h.load_file("stdrun.hoc")

# Create dummy variables representing values set in config file not present here
use_cvode = True
num_chunks = 5
chunk_size = 5000 # in ms

# Create Section
soma = h.Section("soma")

# Create recording vector
t = h.Vector().record(h._ref_t)
h.CVode().active(use_cvode)

h.finitialize(-65)

# Run simulation in chunks of size chunk_size
for chunk_num in range(num_chunks):
	
	# Get start and stop times for this chunk
	chunk_start_time = chunk_num * chunk_size
	chunk_stop_time = (chunk_num + 1) * chunk_size
	
	# If NEURONs h.t not equal to our chunk start time, then make them equal
	if use_cvode:
		if h.t != chunk_start_time:
			h.t = chunk_start_time
			
			# Re-initialize
			h.CVode().re_init()
		
		# Create event so simulation stops at chunk_end_time
		h.CVode().event(chunk_stop_time)
	
	# Run simulation until chunk stop time
	h.continuerun(chunk_stop_time)
	
	# Use custom function to resample recording vectors to desired Dt and save them
	save_chunk_data()
	
	# Resize recording vectors
	h.frecord_init()

In the setup above I was not specifying Dt as an argument for the record method and I would perform the linear interpolation on my recording vectors myself. Now, after realizing I was able to do so, I provide the Dt argument for all of my recording vectors (leaving all else the same) and started to get the following error which would result in a segmentation fault and core dump.

Code: Select all

python: /home/conda/feedstock_root/build_artifacts/neuron_1717659302460/work/src/nrncvode/netcvode.cpp:2269: int NetCvode::global_microstep(): Assertion `tdiff == 0.0 || (gcv_->tstop_begin_ <= tt && tt <= gcv_->tstop_end_)' failed.
Aborted (core dumped)

I didn't really understand this error but looking at my logs where I would output the value of h.t at the beginning and end of each time chunk I realized that:

1.) Despite creating an event at the chunk end times, the chunks were not finishing at the chunk end times and
2.) I would only get this segmentation fault when h.t was less than Dt away from my chunk end time


I then changed my condition for setting h.t and chunk_start_time equal to

Code: Select all

...
	# If the difference between h.t and chunk_start_time is greater than Dt, set them equal to each other.
	if use_cvode:
		if abs(h.t - chunk_start_time) > Dt:
			h.t = chunk_start_time
			
			# Re-initialize
			h.CVode().re_init()
		
		# Create event so simulation stops at chunk_end_time
		h.CVode().event(chunk_stop_time)
...

This fixed the core dumps, but creating the event to ensure my simulation chunks end at the chunk end times still did not make that so. So now I realized an issue that was present the whole time, but did not generate errors before because I was doing the linear interpolation myself. The error being that every now and then my time chunk recording vectors will have one more or less value than they should. I could always just add a few lines of code to catch this and create some small fix, but I feel like there should be a "right" way to do this.

Any help with this is appreciated. I also understand that this question/issue might be more related to my previous post linked above, so if I should post this there as well let me know.