Recording the maximum of a state variable with Brian 2

brian2_max_mon.py

from brian2 import *
# Does also work in standalone mode, enable with:
# set_device('cpp_standalone')

N = 10
G = NeuronGroup(N,
                '''
                dv/dt = -v / (10*ms) + 1*mV*ms**-0.5*xi : volt
                v_max  : volt
                ''')
G.run_regularly('v_max = population_max(t, v)', when='before_start')

@implementation('cpp',
'''
double population_max(double t, double v)
{
    static double last_update = -1.;
    static double highest_value;
    if (t > last_update)
    {
        highest_value = -1e9;  // should be lower than the lowest possible value
        last_update = t;
    }
    if (v > highest_value)
        highest_value = v;
    return highest_value;
}
''')
@check_units(t=second, v=volt, result=volt)
def population_max(t, v):
    # For numpy, this will only be called once with an array of all values
    return max(v)

state_mon = StateMonitor(G, 'v', record=True)
# Only the last entry of v_max contains the actual maximum!
max_mon = StateMonitor(G, 'v_max', record=N-1)
run(100*ms)

plot(state_mon.t/ms, state_mon.v.T/mV, 'k')
plot(state_mon.t/ms, np.max(state_mon.v.T/mV, axis=1), 'r')
plot(state_mon.t/ms, max_mon.v_max[0]/mV, 'bx')

show()