I have a variable ‘k’ set equal to time t. In the threshold argument, I want brian to register a spike every x milliseconds, (where x is an integer, 10ms in example below), for which I check k’s remainder with x. However, brian is missing the condition sometimes (attached plot). The StateMonitor on the other hand shows the values of k modulo x as expected, as shown in the output.
Minimal code to reproduce problem
from brian2 import *
eqs = '''
k = t : second
'''
G = NeuronGroup(1, eqs, method = 'euler', threshold = 'k/ms % 10. == 0.')
S = StateMonitor(G, variables = True, record = True)
M = SpikeMonitor(source = G, variables = 'k')
run(120*ms)
plot(M.t/ms, M.i, 'k.')
show()
for i in S.k[0]:
if i/ms % 10. == 0.:
print(i)
What I have already tried
-
Different values of x (all miss the condition at some point. there seems to be no pattern at which the condition is missed, but the time points at which condition is missed seems to be the same for given x)
-
running on a different machine, and a colab instance : gives same output.
-
using ‘exact’ instead of ‘euler’ as method passed in NeuronGroup : gives same output.
-
no difference on changing to different backends
-
using a differential equation instead : something like dk/dt = 10/(1*ms) : this works perfectly and as expected. However, the question is still valid, because in theory k is taking the same values.
Expected output
Spike should be registered every time the condition is true
Actual output
StateMonitor values whenever same condition is true: (from the for-loop)
0. s
10. ms
20. ms
30. ms
40. ms
50. ms
60. ms
70. ms
80. ms
100. ms
110. ms
note:
(StateMonitor misses 90ms, but that is happening because of rounding errors - 90.0ms is not exactly reached by our variable k. The problem is that 30.0 and 60.0 ms exist as values that k takes but SpikeMonitor isn’t picking them up)
SpikeMonitor plot:
(t = 30, 60 ms are missed)
Sorry, I am not familiar with the internals of the way threshold is implemented, except that it evaluates a boolean. What is the cause of this behavior then? Is it a bug, or a result of the way threshold works? Maybe I’m missing something obvious, but I haven’t been able to solve this. Thank you for your time!
Hi. This comes up every once in a while, and we should probably document this better. At its core, this is not a Brian issue, but due to the way floating point numbers are represented in binary. A number like 0.1 cannot be represented exactly as a floating point number, the actual value is something like 0.1000000000006. Most of the time this is not much of an issue, but if you try an exact comparison with ==, this will often fail. The simplest example I know of is to try 0.1 * 3 == 0.3 in Python which will return False. Or equivalently, and directly relevant for your use case, (0.1 * 0.3) % 0.1 == 0 will return False as well.
To fix this issue, there are two options:
- Instead of an exact comparison, you use something like
abs(k % 10*ms) < 0.00001*ms (this would be the common approach to compare floating point numbers)
- You use Brian’s built-in
timestep function, which is used internally to handle this issues in situations where we need to translate times into an exact number of time steps (e.g. refractoriness conditions). It takes a time and a dt value, and returns an integer time step. In your case, you could use something like timestep(k, dt) % timestep(10*ms, dt) == 0.
All that said: I’m not sure I quite understand why you define k = t : second – why don’t you use t directly in the threshold? Also, if you want to have a regular spike train, the easiest is to use a SpikeGeneratorGroup instead of a NeuronGroup. The following group will emit a spike every 10ms:
G = SpikeGeneratorGroup(1, [0], [0*ms], period=10*ms)
1 Like
Hello Marcel,
Thank you for your reply! I understand that numbers are not represented exactly, but I was confused because the same condition gets satisfied by StateMonitor’s recording, but not satisfied by the threshold as a spike (for the same value ; as shown in actual output)
Thank you for pointing to timestep. that will be quite helpful.
As for my use case, I was trying to see what would happen to a network if any neuron can spike on its own if it hasn’t spiked in last x seconds. So my approach was to put something like this in the threshold argument : 'v > v_thresh or k % x == 0' where k would act like that neuron’s internal clock.
I got around the issue of remainder-check by using :
eqs = '''
dk/dt = 10/(1*ms) : 1
'''
Essentially StateMonitor for k in this case would yield array([ 0., 1., 2., 3., ... ]) if k is starting from zero. This way I can set the initial starting point of k for each neuron, and so I don’t use t directly in the threshold.
And though I am not 100% sure why (My guess is that it is because k is dimensionless now), but doing this always correctly satisfies the remainder condition.(will not work if dk/dt was, say = 1/(1*ms)).
That said I think using timestep would be best here. Please let me know if there is a better way to approach this.
Thank you very much!
Aaradhya
Hi again, sorry for having explained the obvious to you (but maybe not obvious to everyone else), I hope it did not come over as patronizing… I should have read your post a bit more closely, your observation about the StateMonitor / SpikeMonitor is indeed a bit surprising. The reason is that in equations, threshold conditions, etc., Brian may re-arrange terms for optimization purposes (most importantly to only evaluate terms once that are shared across neurons). These operations of course leave the mathematical meaning intact, but they might lead to slightly different rounding which may affect exact comparisons like yours (which should be avoided in general). Concretely, your threshold k/ms % 10. == 0. was actually evaluated as (1/ms * k) % 10. == 0. – if you plug this into your loop, you should get the same behaviour as with the SpikeMonitor.
For your use case, I think I wouldn’t use modulo operations and an exact comparison, but rather a simple >=, comparing the current time to the time of the last spike. This is very similar to what the refractoriness mechanism does. If you had a refractoriness mechanism, you’d automatically have access to the lastspike variable that you could use for that purpose, but without refractoriness you can easily create this variable yourself. I’d formulate your example as something like this
group = NeuronGroup(..., '''
dv/dt = ...
lastspike : second''',
threshold='v > v_thres or (t - lastspike) >= 10*ms',
reset='lastspike = t; v = v_reset')
This way the threshold is very self-explanatory. Now, this formulation above might still run into similar issues because it will trigger a new spike sometimes exactly 10ms after the last one, and sometimes ± 1 timestep later. If you don’t want that, either use a value that is not exactly a multiple of your simulation dt (e.g. replace 10*ms above with 9.99*ms), or use the timestep function: ... or timestep(t - lastspike, dt) >= timestep(10*ms, dt).
Happy modeling!
1 Like
Sorry for the late reply, but this answered both the questions; about the difference in SpikeMonitor/StateMonitor handling the condition, and my specific use case. Thank you very much!
I should have worded my question’s title better, it probably had caused the misunderstanding.
Thanks again,
Aaradhya
1 Like