How to compute the time-average values on the fly?

Support
1

Description of problem

For the time being, Brian2 is able to compute the population averages. But what I want is the average values over time. The simplest way is to monitor the states every timestep. However, in my simulation, I have to run for a long time range. Storing the net states at every timestep will consume a lot of time, memory, and disk space. So I prefer to compute the time average on the fly. And I monitor the states every 5 timesteps to keep the store file not so huge, which made me inable to get the state at every timestep.
In this case, I implemented my average with the following code:

#  v and all v_prev have been defined in equations
neuron.run_regularly('''
    v_prev3 = v_prev2
    v_prev2 = v_prev1
    v_prev1 = v
''', when='before_groups')
run(1*sec)
vs = np.vstack((neuron.v, neuron.v_prev1, neuron.v_prev2, neuron.v_prev3))
v_mean = vs.mean(axis=0)

This is feasible, but very awkward. If I want to average over 100 timesteps, I have to define 99 v_prevs, extract them respectively, which is very clumsy.

So I tried another way. I defined a deque with a fixed length. In each timestep, I compute the average of the deque:

from collections import deque
vs = deque([0] * 5, maxlen=5)

@check_units(v=mV, result=mV)
def var_avg(v):
    global deque
    deque.append(v)
    return sum(v) / 5

eqs = """
...
v_avg = var_avg(v) : volt
"""

Then I ran into the error:

SyntaxError: The subexpression ‘v_avg’ refers to a stateful function (e.g. rand()). Such expressions should only be evaluated once per timestep, add the ‘constant over dt’ flag.

All right, then I added the flag:

eqs = """
...
v_avg = var_avg(v) : volt (constant over dt)
"""

Another error occured this time:

NotImplementedError: Cannot use function ‘var_avg’: "No implementation available for target ‘cython’. Available implementations: "

So I wonder if there is a better way to obtain the time average values.

Minimal code to reproduce problem

from brian2 import *

start_scope()

N = 1000
tau = 10 * ms
vr = -70 * mV
vt0 = -50 * mV
delta_vt0 = 5 * mV
tau_t = 100 * ms
sigma = 0.5 * (vt0 - vr)
v_drive = 2 * (vt0 - vr)
duration = 20 * ms

from collections import deque

vs = deque([0] * 5, maxlen=5)


@check_units(v=mV, result=mV)
def var_avg(v):
    global deque
    deque.append(v)
    return sum(v) / 5


eqs = '''
# dv/dt = (v_drive+vr-v)/tau : volt
dv/dt = (v_drive+vr-v)/tau + sigma*xi*tau**-0.5 : volt 
dvt/dt = (vt0-vt)/tau_t : volt 
v_avg = var_avg(v) : volt (constant over dt)
'''
reset = '''
 v = vr
 vt += delta_vt0  
'''

G = NeuronGroup(N, eqs, threshold='v>vt',
                reset=reset, refractory=5 * ms, method='euler')
G.v = 'rand()*(vt0-vr)+vr'
G.vt = vt0
run(1*ms)

What you have aready tried

  1. average v, v_prev1, …, v_prevN (clusmy)
  2. use a deque (error encountered)
    as described in Description of problem
2

Hi @Jasmine9691. Just to make sure I understand; you want to record every X timesteps, time-averaged value of the membrane potential of the last X timesteps, for each neuron, right? The “clumsy” solution requires you to do many short runs to recover the values, which is quite inefficient (in addition to the clumsiness). Your dequeue solution would work (but I think the code you posted is a bit wrong, since it does not refer to vs at all), but if you provide a function implementation, it needs to be in the target language of the code generation. By default, this is Cython, so implementing the function in Python is not enough. You could either: 1) switch your code generation target to Python (prefs.codegen.target = 'numpy'), 2) provide an implementation for Cython (not trivial), or 3) use a network_operation instead (see Input stimuli — Brian 2 2.8.0 documentation), which can run arbitrary Python code even if the target language is Cython.
But all of these solutions are a bit wasteful, since you do not need to store the previous values just to calculate the average, just the sum should be enough.

So something like the following should do the trick:

eqs = '''...
v_avg : volt
'''
G = NeuronGroup(...)
# Average over 100 timesteps
G.run_regularly("v_avg += v/100.0", when='after_groups')
G.run_regularly("v_avg = 0*mV", dt=100*defaultclock.dt, when='end')
avg_mon = StateMonitor(G, "v_avg", record=True, dt=100*defaultclock.dt)

If you want to avoid a weird value at the very first recording, you can initialize v_avg with something like vr.

Hope that work for you!

[Edit: The StateMonitor above was missing a record=True]

1 Like
3

Thanks a lot @mstimberg . I have tried network operation and run_regularly. They both worked for me. But I prefer the former ( network operation) because I can get the time-average values at any timestep. For example, I can average v(1), v(2) and v(3) to obtain v_avg(3); and average v(2), v(3), v(4) to obtain v_avg(4). With run_regularly, I can only get average values at specified intervals, which is a bit limited for me.

but I think the code you posted is a bit wrong, since it does not refer to vs at all

Yes, it has some typos. Also, the deque only applies to a single neuron. For multiple neurons, I think a numpy array is necessary. My implementations of time-average are provided below in case anybody else needs it.

  1. run_regularly (nearly the same as @mstimberg 's solution)
eqs = '''
v_avg : volt
'''
G = NeuronGroup(...)
G.run_regularly('''
    v_avg += v / n_callback
''', when='after_groups')
G.run_regularly('v_avg = 0', dt=n_callback * defaultclock.dt, when='start')
  1. network_operation
eqs = '''
v_avg : volt
'''
G = NeuronGroup(n_pop, eqs)
vs = zeros((n_pop, n_callback), dtype=float)

@network_operation(dt=defaultclock.dt, when='after_groups')
def update_v_avg():
    nonlocal vs
    vs = hstack((vs[:, 1:], G.v[:, newaxis]))
    G.v_avg = vs.mean(axis=-1)
4

Hmm, but if you record time-averaged values at every timestep, you don’t save any memory, do you? If that’s what you are looking for, you can also simply record the values at every timestep without any network_operation or run_regularly, and then calculate the rolling average afterwards. There are various ways to do this efficiently, e.g. the following should work:

w = 10  # average over 10 time steps
average = np.cumsum(state_monitor.v[:], axis=1)
average[:, w:] = average[:, w:] - average[:, :-w]
average = average[:, w-1:] / w
5

I think monitoring at each timestep is different from using network operation.
For example, I run 10000 steps, and set the average window length as 4. If I monitor at each timestep, I would get a numpy array with shape of (n_pop, 10000). But if I use network operation above, the numpy array is very small, with shape of (n_pop, 4). I only need to save four values for each neuron.

6

Hmm, no, I don’t think so. Your vs array would be small, but if you record G.v_avg at every time step, you’d have as many values as when you record G.v itself. But you are right that the two approaches are not perfectly identical: with your network_operation, you have the option to e.g. average over the last 100 time steps, but to record the value every, say, 10 time steps and get meaningful values. The disadvantage of the network_operation method is that it can slow down a simulation significantly, but of course whether this matters depends how complex your simulation is in general. Just to mention it: another common approach to have a smooth version of the membrane potential, very similar to the moving average is to use a low-pass filter in the differential equations, which can also be interpreted as exponential smoothing (aka exponential moving average). You’d add

dv_avg/dt = (v - v_avg)/tau_average : volt

to your equation, and set tau_average according to your desired smoothing. You can then also decide how often you want to record these values, independent of the tau_average value.

7

I just carried out a performance test. Two cases are tested.

Case 1: net.run(1*second)

Network operation: 1. s (100%) simulated in 8s
run_regularly: 1. s (100%) simulated in 5s

Case 2: net.run(0.4*ms) for 100 times in a loop

Network operation: 1 min 2 sec
run_regularly: 1 min 3 sec

In my project, I’m coupling brian2 with another package (which is the driver code) in the following way:

# callback function of another package
def callback(...):
    ... # get stim from another package
    net['stim'] = stim
    net.run(dt)
    v_avg = net['G'].v_avg
    return v_avg

This is more like Case 2, i.e., I must run brian2 at regular intervals. Although run_regularly is faster than network_operation in Case 1, they are almost as fast as each other in Case 2.

Exponential smoothing is indeed a faster approach to get the average. But it would assign bigger weights to newer values. In my simulation, I want to assign the same weight to the values to be averaged. So arithmetic mean is chosen.

Anyway, thanks for your kind reminder!

8

Just as some final remarks on your use case: in your example, the details of the implementation do indeed not matter, since the overhead from calling run repeatedly dominates everything.
If performance matters, you could use Brian’s “C++ standalone mode”, and in particular its run_args feature, as described here: Computational methods and efficiency — Brian 2 2.8.0 documentation This would probably run faster than Case 1, even when using the approach from Case 2 (not counting the time that the other package takes to give you the stimulus). If you want to try this and run into issues, happy to help (but probably best in a new discussion thread)!

1 Like