run_regularly() with large dt in standalone mode

[felix11h]

I'm running into an issue with using run_regularly() in standalone mode that I believe could possibly be a bug.

Below is the network that reproduces the issue. The simulation segment 'T2', although very long, is completely skipped. For example I get an output of

Starting simulation at t=0 s for duration 1 s 1 s (100%) simulated in 0.23194 s Starting simulation at t=1 s for duration 100000 s 100000 s (100%) simulated in 1.35702e-94 s Starting simulation at t=100001 s for duration 1 s 1 s (100%) simulated in 0.228624 s

Two things fix the problem:

1. using a shorter simluation time, for example t2=100*second, or,
2. disabling the run_regularly() step

Brian2 version '2.2.1'. Can anyone reproduce this? Thanks!

import numpy as np
from brian2 import *

N_e, N_i = 50, 50

tau = 20*ms               
tau_e = 3*ms              
tau_i = 5*ms              
El = -60*mV               
Ee = 0*mV                 
Ei = -80*mV               
mu_e = 9.0*mV
mu_i = 8.5*mV
sigma_e = 0.5**0.5 *mV    
sigma_i = 0.5**0.5 *mV

Vr_e = -60*mV
Vr_i = -60*mV
Vt_e = -57.5*mV
Vt_i = -58*mV

a_ee = 0.002

taupre = 15*ms
taupost = 30*ms
Aplus = 0.0015/10
Aminus = -0.00075/10
amax = 0.8
ATotalMax = 0.15

T1 = 1*second
T2 = 100000*second
T3 = 1*second

insert_P = 0.001
p_inactivate = 0.25
p_ee = 0.1

condlif_memnoise = '''
              dV/dt = (El-V + (gfwd+ge)*(Ee-V) + gi*(Ei-V))/tau +  mu/tau + (sigma * xi) / (tau **.5) : volt
              Vt : volt 
              dge /dt = -ge/tau_e : 1
              dgfwd /dt = -gfwd/tau_e : 1
              dgi /dt = -gi/tau_i : 1

              Asum : 1

              sigma: volt (constant)
              mu : volt (constant)
              '''

nrnEE_thrshld = 'V > Vt'

nrnEE_reset = 'V = Vr_e'

synEE_mod = '''
            a : 1
            syn_active : integer

            dApre  /dt = -Apre/taupre  : 1 (event-driven)
            dApost /dt = -Apost/taupost : 1 (event-driven)

            Asum_post = a : 1 (summed)         
            insert_P : 1 (shared) 
            p_inactivate : 1 (shared)
            '''

synEE_p_activate = '''
                   r = rand()
                   syn_active = int(r < p_ee)
                   a = syn_active*a
                   '''

synEE_pre = '''
            ge_post += syn_active*a
            Apre += syn_active*Aplus
            '''

synEE_post = '''
             Apost+= syn_active*Aminus
             '''

set_device('cpp_standalone', directory='./builds/0000',
           build_on_run=False)

T = T1 + T2 + T3

neuron_model = condlif_memnoise

GExc = NeuronGroup(N=N_e, model=neuron_model,
                   threshold=nrnEE_thrshld,
                   reset=nrnEE_reset)
GInh = NeuronGroup(N=N_i, model=neuron_model,
                   threshold ='V > Vt',
                   reset='V=Vr_i')

GExc.mu, GInh.mu = mu_e, mu_i
GExc.sigma, GInh.sigma = sigma_e, sigma_i

GExc.Vt, GInh.Vt = Vt_e, Vt_i
GExc.V , GInh.V  = np.random.uniform(Vr_e/mV, Vt_e/mV,
                                     size=N_e)*mV, \
                   np.random.uniform(Vr_i/mV, Vt_i/mV,
                                     size=N_i)*mV

synEE_pre_mod = synEE_pre
synEE_post_mod = synEE_post

SynEE = Synapses(target=GExc, source=GExc, model=synEE_mod,
                 on_pre=synEE_pre_mod, on_post=synEE_post_mod)

SynEE.connect(p=1)

SynEE.a = a_ee        
SynEE.insert_P = insert_P
SynEE.p_inactivate = p_inactivate

# this causes problems when T2 is large!
# compare T2 = 100000*second with T2= 100*second
SynEE.run_regularly(synEE_p_activate, dt=T, when='start',
                    order=-100)

net = Network(GExc, GInh, SynEE)

net.run(T1, report='text')

net.run(T2, report='text')

net.run(T3, report='text')

device.build(directory='builds/0000', clean=True,
             compile=True, run=True, debug=False)

[mstimberg]

Hi Felix, thanks for the report, I can reproduce the issue with your example script. It seems to be specific to C++ standalone mode, I could imagine that there is some conversion of time to integer time steps at some point that goes over 32bit bounds, or something along these lines. I'll have a closer look soon.

[felix11h]

Thanks for the fast reply! I'm very curious to learn what's happening here. I wasn't successful in reproducing the behaviour in a minimal example network, so I've cut my current simulation down to what I hope is a manageable size to track down the issue.

The initial idea was to set the dt of the run_regularly() to the duration of the simulation in order to do something once at the beginning of the simulation (activate some units). I was able to do this directly now without using a model equation. I still wanted to report the issue though as it might be important!

[mstimberg]

I looked into this a bit, and I no longer think it has to do with 32bit integers. It is rather about some rounding tests where we consider things to fall into the same time step if they differ by less than a 10000th of dt – in your example, the first run time is less than that (for the dt of your run_regularly) and I think this messes things up. I also noticed that this is not only a standalone problem, with runtime targets you can get similar issues. I'll need to look into this more thoroughly, but I think I have a general idea what the problem is.

[felix11h]

As mentioned in #1057 , removing the run_regularly statement but simply increasing the simulation time to T2 = 500000*second produces similar behaviour:

Starting simulation at t=0 s for duration 1 s 1 s (100%) simulated in 0.322394 s Starting simulation at t=1 s for duration 500000 s 500000 s (100%) simulated in 0 s Starting simulation at t=500001 s for duration 1 s 1 s (100%) simulated in 0 s

This issue does not seem to be fixed by https://github.com/brian-team/brian2/commit/a59121462dafabf59aefbb74a4cc0cea966c6499.