Fail to simulate a long run in cpp_standalone mode.

[NeoNeuron]

Issue statement

I fail to simulate an LIF network with 2 neurons for a long duration run (1e6*second) in cpp_standalone mode. NO warning or error appears, and no data of spike trains is generated.

However, if I switch to a shorter run duration, like 1e5*second, the code works fine. And also, if I split the whole run duration into multiple pieces, such as replacing run(duration) with

for _ in range(10):
    run(duration/10)

the code works fine. I'm wondering is there a maximum run duration value in Brian2 cpp_standalone mode, or is it due to the data type issue? Could you please give me some help?

I attach my code below, and the command to run the simulation case I mentioned is below.

python simulate_LIF.py -T 1e6 --ref 0 -v 20 -g 10 -f 0.4 -u 0.1 --verbose

Versions

• python: 3.7.7
• brian2: 2.3.0.2

# %%

if __name__ == '__main__':
    # %%
    import numpy as np
    import sys, os
    import subprocess as sp
    from utils import save2bin
    from brian2 import *
    import argparse
    import argparse
    parser = argparse.ArgumentParser()
    parser.add_argument("-T",       type=float, default=10000, help="simulation period, ms", )
    parser.add_argument("--ref",    type=float, default=2, help="refractory period", )
    parser.add_argument('-v', "--tauv", type=float, default=20, help="time constant of v", )
    parser.add_argument('-g', "--taug", type=float, default=5, help="time constant of g", )
    parser.add_argument("-f",       type=float, default=0.4, help="Poisson strength", )
    parser.add_argument("-u",       type=float, default=0.4, help="Poisson rate, in kHz", )
    parser.add_argument("--verbose", action='store_true', help="show progress bar", )
    args = parser.parse_args()

    class ProgressBar(object):
        def __init__(self, toolbar_width=40):
            self.toolbar_width = toolbar_width
            self.ticks = 0

        def __call__(self, elapsed, complete, start, duration):
            if complete == 0.0:
                # setup toolbar
                sys.stdout.write("[%s]" % (" " * self.toolbar_width))
                sys.stdout.flush()
                sys.stdout.write("\b" * (self.toolbar_width + 1)) # return to start of line, after '['
            else:
                ticks_needed = int(round(complete * self.toolbar_width))
                if self.ticks < ticks_needed:
                    sys.stdout.write("=" * (ticks_needed-self.ticks))
                    sys.stdout.flush()
                    self.ticks = ticks_needed
            if complete == 1.0:
                sys.stdout.write("\n")

    pid = os.getpid()
    directory = f"standalone{pid}"
    set_device('cpp_standalone', directory=directory, build_on_run=False)
    prefs.devices.cpp_standalone.extra_make_args_unix = ['-j1']
    N = 2
    Vr = -65*mV
    theta = -40*mV
    tau_v = args.tauv*ms
    tau_g = args.taug*ms
    delay = 0*ms
    ref = args.ref*ms
    duration = args.T*second
    J = 2e-2
    ffwd_J = args.f
    freq =   args.u*kHz
    muext = 25*mV
    sigmaext = 1*mV
    E_Na = 0*mV

    eqs = """
    dV/dt = (-(V-Vr) - g*(V-E_Na))/tau_v : volt (unless refractory)
    dg/dt = -g/tau_g : 1
    """
    # %%
    start_scope()
    group = NeuronGroup(N, eqs, threshold='V>theta',
                        reset='V=Vr', refractory=ref, method='euler')
    group.V = Vr
    group.g = 0.
    conn = Synapses(group, group, on_pre='g += J', delay=delay)
    for idx in [(0,1)]: #, (1,2), (0,3), (1,3), (0,4)]:
        conn.connect(i=idx[0], j=idx[1])
    adjmat = np.eye(2, k=1, dtype=float)
    P = PoissonGroup(N, freq, name='ffwd_Poisson')
    ffwd_conn = Synapses(P, group, on_pre='g += ffwd_J')
    ffwd_conn.connect('i==j')

    M = SpikeMonitor(group)
    # V = StateMonitor(group, 'V', record=0)

    report_func = '''
        int remaining = (int)((1-completed)/completed*elapsed+0.5);
        if (completed == 0.0)
        {
            std::cout << "Starting simulation at t=" << start << " s for duration " << duration << " s"<<std::flush;
        }
        else
        {
            int barWidth = 70;
            std::cout << "\\r[";
            int pos = barWidth * completed;
            for (int i = 0; i < barWidth; ++i) {
                    if (i < pos) std::cout << "=";
                    else if (i == pos) std::cout << ">";
                    else std::cout << " ";
            }
            std::cout << "] " << int(completed * 100.0) << "% completed. | "<<int(remaining) <<"s remaining"<<std::flush;
        }
    '''
    if args.verbose:
        # for _ in range(10):
        #    run(duration/10, report=report_func, report_period=1*second)
        run(duration, report=ProgressBar(), report_period=1*second)
        device.build(directory=directory, compile=True, run=True, debug=False,)
    else:
        # for _ in range(10):
        #    run(duration/10)
        run(duration)
        device.build(directory=directory, compile=True, run=True, debug=False,)
    try:
        assert M.i.shape[0] > 0
        print('\nMfr : ', M.i.shape[0]*1.0/N/duration)
    except AssertionError:
        print('\nNo spike is generated.')

    spk_data = np.vstack((M.t/ms, M.i)).T
    # %%
    data_out = np.vstack((M.t/ms, M.i)).T
    save2bin(
             f"LIFp=0.25s={J:.3f}" +\
             f"ref={ref/ms:.1f}" +\
             f"tauv={tau_v/ms:.1f}" +\
             f"taug={tau_g/ms:.1f}" +\
             f"f={ffwd_J:.3f}" +\
             f"u={freq/kHz:.3f}" +\
             "_spike_train.dat", data_out, verbose=True)

    save2bin(f"connect_matrix-p=0.250.dat", adjmat)

    # %%
    # delete tmp folder
    sp.run(['rm', '-rf', directory])

[mstimberg]

Hi @NeoNeuron. Thanks for the report. This is an actual bug in Brian's C++ standalone mode. Internally, Brian uses integer time steps to avoid rounding issues, (i.e. instead of doing t += dt every time step, we do i += 1; t = i*dt). When setting this up in the beginning of a run, the code by mistake casts the time step to a 32bit integer, and it can therefore not handle more than ~2·10⁹ time steps (your 10⁶s correspond to 10¹⁰ time steps with the default dt). Basically it wraps over and the "end timestep" becomes smaller than the "start time step", and nothing gets simulated.
I did actually fix this issue a couple of years ago as part of #1057, but progress stalled on this PR (which was originally about an unrelated issue about using multiple clocks with dts of different orders of magnitude). I am going to pull out the fix into an independent PR that we can hopefully merge soon.
In case you can't wait and want to fix things manually, you can manually replace all mentions of int by int64_t in brian2/devices/cpp_standalone/brianlib/clocks.h.

[NeoNeuron]

@mstimberg Thanks for quick reply. And your advice for temporal fixing the issue works for me. 👍