Inaccurate updates when having multiple differential equations in neuron groups

[jesusgf96]

Hi there! In Brian2GeNN, when having multiple differential equations in a neuron group, some are not updated continually. Probably, for efficiency reasons, they are not constantly updated when not needed in the neurons. However, when they are needed in the synapses, this can be a problem.
Below, you can find a toy script that shows this possible problem. Two equations for two variables ('v' and 'r') are included in a neuron group. The second variable is not needed to calculate the updates in the first one. However, it is needed to update the synapses. If we place a monitor in this second variable, it surprisingly updates correctly (as it is explicitly needed in every time step). When the monitor is removed, it updates wrongly. This behaviour can be tested by trying different values for the variable 'monitor' (True/False) in different executions. Also, it doesn't occur when using Brian2 alone (no GeNN backend). You can test it by changing the value of the variable 'GPU' to False.
I attached some images with the resulting distribution of the synaptic values for these different situations. The values only change when not tracking this second variable (monitor=False) using the GeNN backend (GPU=True). In this toy example, there's only a slight difference. But when using large networks in long simulations, the changes are immense.

ToyExample_Brian2GeNN.zip

[mstimberg]

Umm, many thanks for reporting this, this looks indeed like a serious issue! I will look into this in more detail at the end of the week.

[mstimberg]

I started looking into this, but I did not yet get completely to the bottom of it. The discrepancy between Brian2GeNN and C++ standalone is partly expected for your model, but I'm still figuring out what is going when you add the monitor (in my simplified example, the monitor does not seem to make any difference anymore, somehow it also seems to require a run_regularly operation 🤷‍♂️ ).
Why is the discrepancy expected? In your model, you have a differential equation updating r, a reset operation updating r after a spike, a differential equation updating w based on r of the pre and post-synaptic cell, and finally a synaptic on_pre operation updating the post-synaptic v based on the current synaptic weight w. The problem is that all these things happen during a single time step. In particular, a pre-synaptic spike triggers both an update of r in the pre-synaptic cell, and a propagation of the current value of w to the post-synaptic cell. In a way, this is "ill-defined" – should we first update r, then update w on the current value of r_pre and r_post and then use this w for the propagation? Or should we e.g. first transmit the spike by using the current w, then update w based on its equation and finally "reset" r to its new value?
In Brian, we try to make all possible orderings possible and provide tools to define (and check) what exactly is going on, see the scheduling documentation and also this related blog post. In contrast to this, Brian2GeNN/GeNN use a fixed schedule which is not exactly the same as Brian's default schedule (this is why you get an info message about the changed "default_schedule" preference) – unfortunately, we don't make this as clear as we should.
In Brian, by default we first update all state variables (differential equations), then check thresholds, propagate any resulting spikes, and finally apply the reset. In Brian2GeNN, we first update synaptic state variables, propagate all pending spikes, update neuronal state variables, check thresholds, and apply reset statements.
If you change the global schedule to

magic_network.schedule = ['start', 'synapses', 'groups', 'thresholds', 'resets', 'end']

and also move the state update of the w variable to an earlier time via:

S.state_updater.when = 'before_synapses'

then you should have the same basic scheduling in Brian as in Brian2GeNN.

As I said earlier, there is probably still some bug in Brian2GeNN (I'm suspecting related to moving memory from/to the GPU), since obviously results should not change with the addition of a StateMonitor. I will try to have another look next week.

[mstimberg]

It took me a long time to get back to this, but I finally figured out what is going on. As so often, the issue was very puzzling but the final explanation (and the fix) is quite simple... Here's a quick explanation: Brian2GeNN mostly runs things on the GPU (that's the whole idea, of course), but for simplicity, run_regulary operations and monitors are run on the CPU, reusing standard Brian C++ code. In order for this to work correctly, the relevant variables have to be copied between the GPU and the CPU: we need GPU → CPU copies for the monitors, GPU → CPU copies for all the variables that are read by the run_regulary operation, and finally CPU → GPU copies for all variables that were changed by the run_regularly operation. The PR #143 by @kernfel optimized the GPU → CPU copies by only copying the required variables, instead of the full state of a neuron/synapse. In that case, it was only a matter of performance, though. What I did not realize before, is that the CPU → GPU copy for a run_regularly operation also copies the full state of a neuron/synapse, instead of only the required variable. This is wasteful performance-wise, but since most run_regulary operations are executed with a large dt, this shouldn't matter much. What is more important, is that this can actually introduce bugs, as in your example: if a variable has been updated on the GPU, and is not affected by the run_regulary operation (like the variable r in your example), it will be overwritten with its initial value (the only value in CPU memory) after every run_regulary operation. This gets "fixed" when a monitor is introduced that records this variable: now, the variable (again, r in your case) gets copied over to CPU memory every time step, and the unnecessary copy from CPU → GPU therefore does not have any effect.
The fix is quite simple: instead of copying over the full state after every run_regulary execution, only copy over the changed variable (I in your example). I'll prepare a proper fix next week (of course, if e.g. @kernfel wants to have a go at it, please go ahead 😄). One minor complication: the functions to push neuronal variables to the GPU have the word Current in it (e.g. pushCurrentIneurongroupToDevice to push transfer the values of I), while the same functions for synapses do not (e.g. pushwsynapsesToDevice to push the variable w).