Make `SpikeGeneratorGroup` work for `run_args`

[mstimberg]

The run_args feature is very useful for simulations where you want to present different stimuli to a network, but currently this only works if you can represent the stimulus as the initial value of a state variable (e.g. the rates of a PoissonGroup), or as the values in a TimedArray. You cannot change the values of a SpikeGeneratorGroup, since it currently checks and prepares the data structures in Python code. If we moved this preparation into the templates, it would become possible to switch out the spikes indices/times for every run via run_args. See this related discussion: https://brian.discourse.group/t/dynamically-loading-spiketimes-in-a-spikegeneratorgroup-on-c-standalone-mode/1369/3

[alikhere]

Hi @mstimberg I’m Md Khurshid Alam, a final-year B.Tech student skilled in C++, Python, JavaScript, algorithms, data structures, and open-source development. I would love to work on this!

[mstimberg]

Hi @alikhere, glad to hear that! But please note that this task is not as trivial as it might look, as it requires several steps to get things to work.

The first step towards this would be to move the Python code in SpikeGeneratorGroup.before_run

brian2/brian2/input/spikegeneratorgroup.py

Lines 207 to 283 in 83847e3

	def before_run(self, run_namespace):
	# Do some checks on the period vs. dt
	dt = self.dt_[:] # make a copy
	period = self.period_
	if period < np.inf and period != 0:
	if period < dt:
	raise ValueError(
	f"The period of '{self.name}' is {self.period[:]!s}, "
	f"which is smaller than its dt of {dt*second!s}."
	)
	if self._spikes_changed:
	current_t = self.variables["t"].get_value().item()
	timesteps = timestep(self._spike_time, dt)
	current_step = timestep(current_t, dt)
	in_the_past = np.nonzero(timesteps < current_step)[0]
	if len(in_the_past):
	logger.warn(
	"The SpikeGeneratorGroup contains spike times "
	"earlier than the start time of the current run "
	f"(t = {current_t*second!s}), these spikes will be "
	"ignored.",
	name_suffix="ignored_spikes",
	)
	self.variables["_lastindex"].set_value(in_the_past[-1] + 1)
	else:
	self.variables["_lastindex"].set_value(0)
	# Check that we don't have more than one spike per neuron in a time bin
	if self._previous_dt is None or dt != self._previous_dt or self._spikes_changed:
	# We shift all the spikes by a tiny amount to make sure that spikes
	# at exact multiples of dt do not end up in the previous time bin
	# This shift has to be quite significant relative to machine
	# epsilon, we use 1e-3 of the dt here
	shift = 1e-3 * dt
	timebins = np.asarray(
	np.asarray(self._spike_time + shift) / dt, dtype=np.int32
	)
	# time is already in sorted order, so it's enough to check if the condition
	# that timebins[i]==timebins[i+1] and self._neuron_index[i]==self._neuron_index[i+1]
	# is ever both true
	if (
	np.logical_and(np.diff(timebins) == 0, np.diff(self._neuron_index) == 0)
	).any():
	raise ValueError(
	f"Using a dt of {self.dt!s}, some neurons of "
	f"SpikeGeneratorGroup '{self.name}' spike more than "
	"once during a time step."
	)
	self.variables["_timebins"].set_value(timebins)
	period_bins = np.round(period / dt)
	max_int = np.iinfo(np.int32).max
	if period_bins > max_int:
	logger.warn(
	f"Periods longer than {max_int} timesteps "
	f"(={max_int*dt*second!s}) are not "
	"supported, the period will therefore be "
	"considered infinite. Set the period to 0*second "
	"to avoid this "
	"warning.",
	"spikegenerator_long_period",
	)
	period = period_bins = 0
	if np.abs(period_bins * dt - period) > period * np.finfo(dt.dtype).eps:
	raise NotImplementedError(
	f"The period of '{self.name}' is "
	f"{self.period[:]!s}, which is "
	"not an integer multiple of its dt "
	f"of {dt*second!s}."
	)
	self.variables["_period_bins"].set_value(period_bins)
	self._previous_dt = dt
	self._spikes_changed = False
	super().before_run(run_namespace=run_namespace)

into the before_code blocks of the numpy, Cython, and C++ templates. I'd ignore Cython for now, and only focus on numpy and C++. I'd start with numpy, since this is mostly a copy of the existing code, except for replacing expressions such as self.variables["t"].get_value().item() by {{t}}[0] – I know that dealing with these Jinja2 templates is not obvious, though 😬 It might be helpful to print out the generated code to get an idea of whether it is what you expect. Also, I'd ignore 1) all the error/warning checks for now, and only focus on getting the "correct case" to work and 2) ignore the checks about whether it is actually necessary to re-calculate things – simply always recalculate things. If all this works, try converting the code to C++ in the corresponding C++ template. Once you get this far, we can discuss the following steps 😊

[mahipalimkar]

Hi @alikhere ! I came across this issue, and wanted to confirm if you are still working on it. If you are not, I would love to take it forward. Let me know. Thanks !😊

[alikhere]

You can work.

[mahipalimkar]

Hi @mstimberg! I have completed moving the python code into the before_code blocks of numpy and c++ as you explained above. What are the next steps?

[alikhere]

Hi @mstimberg! I have completed moving the python code into the before_code blocks of numpy and c++ as you explained above. What are the next steps?

Raise a PR i guess.

[mahipalimkar]

I don't think this is all that has to be done to solve the issue.

[mstimberg]

Agreed, please open a draft PR so I can have a look !

[mahipalimkar]

Okay @mstimberg I will do that at the earliest

[maverick4code]

@mstimberg,
I have some doubts,

1. What exactly needs to be calculated in the template(like time_bins, period_bins) like what’s the core job here?
2. How can I check if the template is working as expected?

[mstimberg]

@maverick4code

What exactly needs to be calculated in the template(like time_bins, period_bins) like what’s the core job here?

The core job is converting the floating point times and period into their integer equivalents, based on the current time step. Later, it should also do some related checks (e.g. whether the time step size means that multiple events fall into the same bin), but for a first version, I'd assume that everything is correct.

How can I check if the template is working as expected?

I'd either write a simple script that uses a SpikeGeneratorGroup, or use an existing example, e.g. https://brian2.readthedocs.io/en/stable/examples/frompapers.Diesmann_et_al_1999.html

To check whether this works with the Python template (assuming you also deleted the previous SpikeGeneratorGroup.before_run function), you'd add prefs.codegen.target = 'numpy' to the beginning . To check it with the C++ standalone template, you'd instead use set_device('cpp_standalone').