Error while running parallel standalone simulation in Brian2?

Support
,
1

Description of problem

I have a simulation that runs fine without the cpp_standalone parallelization. However, if I do that with parallelization described here, it throws error. What am I doing wrong ? :smiling_face_with_tear:

Minimal code to reproduce problem



#define overlap percentage
from IPython.display import clear_output
import numpy as np
import matplotlib.pyplot as plt
import os
import shutil
from brian2 import *
from tqdm import tqdm


# #delete the previous build directory if it exists
if os.path.exists('output'):
    shutil.rmtree('output')


set_device('cpp_standalone', build_on_run=False)

prefs.devices.cpp_standalone.openmp_threads = 20
n_ab = 0.025 # this percentage recieves both sine and pulse
n_a = (1 - n_ab)/2 # only sine
n_b = (1 - n_ab)/2 # only pulse


########## neuron inputs #############
start_scope()
defaultclock.dt = 0.05*ms
duration = 4*second

Cm = 1*uF  # /cm**2
gL = 0.1*msiemens
EL = -65*mV
ENa = 55*mV
EK = -90*mV
gNa = 35*msiemens
gK = 9*msiemens

# Synaptic variables
taue = 5*ms
taui = 10*ms
Ee = 0*mV
Ei = -80*mV
we = 15*uS  # excitatory synaptic weight
wi = 67*uS  # inhibitory synaptic weight

# Hodgkin-Huxley + synaptic + external currents
eqs = '''
dv/dt = (-gNa*m**3*h*(v-ENa) - gK*n**4*(v-EK) - gL*(v-EL)
         + Iamp0 + Iamp * sin(2*pi*Ifreq*t)
         + ge*(Ee-v) + gi*(Ei-v)) / Cm : volt

# gating vars
m = alpha_m/(alpha_m+beta_m) : 1
alpha_m = 0.1/mV*10*mV/exprel(-(v+35*mV)/(10*mV))/ms : Hz
beta_m = 4*exp(-(v+60*mV)/(18*mV))/ms : Hz

dh/dt = 5*(alpha_h*(1-h)-beta_h*h) : 1
alpha_h = 0.07*exp(-(v+58*mV)/(20*mV))/ms : Hz
beta_h = 1./(exp(-0.1/mV*(v+28*mV))+1)/ms : Hz

dn/dt = 5*(alpha_n*(1-n)-beta_n*n) : 1
alpha_n = 0.01/mV*10*mV/exprel(-(v+34*mV)/(10*mV))/ms : Hz
beta_n = 0.125*exp(-(v+44*mV)/(80*mV))/ms : Hz

# synaptic conductances
dge/dt = -ge*(1./taue) : siemens
dgi/dt = -gi*(1./taui) : siemens

# external input parameters
Iamp0 : amp   # constant or pulsed current
Iamp : amp    # sinusoidal amplitude
Ifreq : Hz    # sinusoidal frequency
'''

# Neuron group
N = 100
group = NeuronGroup(N, eqs, threshold='v>-20*mV', refractory=3*ms,
                    method='rk4')
group.v = -70*mV
group.h = 1


######### noise goes everywhere ##########
# Noise input
rate = 20*Hz
Noise_in = PoissonGroup(N, rates=rate)
S = Synapses(Noise_in, group, on_pre='ge += we')
S.connect(j='i')

# Example: set input currents in n_a neurons

n_a_arr = np.zeros(N)
# first n_a percentage becomes 1
n_a = n_a * N  # convert percentage to number of neurons
n_a_arr[:int(n_a)] = 1  # set first n_a neurons to 1  

# sinusoidal input
n_b_arr = np.zeros(N)
#last n_b percentage becomes 1
n_b = n_b * N  # convert percentage to number of neurons
n_b_arr[-int(n_b):] = 1  # set last n_b neurons to 1

n_ab_arr = n_a_arr + n_b_arr  # overlap neurons
#make 0s as 1s and 1s as 0s
n_ab_arr[n_ab_arr == 0] = 1
n_ab_arr[n_ab_arr == 1] = 0

# Monitors
monitor = StateMonitor(group, 'v', record=True)
monitor_s = SpikeMonitor(group)
# syn_rate = PopulationRateMonitor(Noise_in)
# monitor_rate = PopulationRateMonitor(group)
# ge_monitor = StateMonitor(group, 'ge', record=True)
# input_monitor = SpikeMonitor(Noise_in)


n_sims = 100

for i in tqdm(range(n_sims)):    
    group.Iamp0 = 0*0.55*uA * (n_a_arr + n_ab_arr)   # baseline or pulse amplitude
    group.Iamp  = 0*0.35*uA * (n_b_arr + n_ab_arr)  # sinusoidal amplitude
    group.Ifreq = 8*Hz     # sinusoidal frequency

    run(duration/5.0)

    group.Iamp0 = 0.55*uA * (n_a_arr + n_ab_arr)   # baseline or pulse amplitude
    group.Iamp  = 0.35*uA * (n_b_arr + n_ab_arr)  # sinusoidal amplitude
    group.Ifreq = 8*Hz     # sinusoidal frequency

    run(3*duration/5.0)

    group.Iamp0 = 0*0.55*uA * (n_a_arr + n_ab_arr)   # baseline or pulse amplitude
    group.Iamp  = 0*0.35*uA * (n_b_arr + n_ab_arr)  # sinusoidal amplitude
    group.Ifreq = 8*Hz     # sinusoidal frequency   

    run(duration/5.0)
device.build()
print('Whole sim completed')


#plot the raster plot for neuron group
plt.figure(figsize=(10, 6), dpi=150)
plt.plot(monitor_s.t/ms, monitor_s.i, 'k|', markersize= 3)
plt.xlabel('Time (ms)', fontsize=12)
plt.ylabel('Neuron index', fontsize=12)
plt.title('Raster plot of neuron spikes', fontsize=14)
plt.tight_layout()
plt.xlim(0, 2*(duration/ms))
plt.show() ```

What you have already tried

I have already tried installing brian2 on conda environment with python3.13 and python3.12, using both pip and conda. I have tried this on my linux system and macOS. Same error each time.

EDIT: It runs fine if you reduce the n_sims to 100 –> 10 . Or when the size is reduced.

Expected output (if relevant)

if you use the same code after commenting the standalone activation lines (below), it works just fine but takes a long time. Here is the output attached also, that looks fine.



#define overlap percentage
from IPython.display import clear_output
import numpy as np
import matplotlib.pyplot as plt
import os
import shutil
from brian2 import *
from tqdm import tqdm


# #delete the previous build directory if it exists
# if os.path.exists('output'):
#     shutil.rmtree('output')


# set_device('cpp_standalone', build_on_run=False)

# prefs.devices.cpp_standalone.openmp_threads = 20
n_ab = 0.025 # this percentage recieves both sine and pulse
n_a = (1 - n_ab)/2 # only sine
n_b = (1 - n_ab)/2 # only pulse


########## neuron inputs #############
start_scope()
defaultclock.dt = 0.05*ms
duration = 4*second

Cm = 1*uF  # /cm**2
gL = 0.1*msiemens
EL = -65*mV
ENa = 55*mV
EK = -90*mV
gNa = 35*msiemens
gK = 9*msiemens

# Synaptic variables
taue = 5*ms
taui = 10*ms
Ee = 0*mV
Ei = -80*mV
we = 15*uS  # excitatory synaptic weight
wi = 67*uS  # inhibitory synaptic weight

# Hodgkin-Huxley + synaptic + external currents
eqs = '''
dv/dt = (-gNa*m**3*h*(v-ENa) - gK*n**4*(v-EK) - gL*(v-EL)
         + Iamp0 + Iamp * sin(2*pi*Ifreq*t)
         + ge*(Ee-v) + gi*(Ei-v)) / Cm : volt

# gating vars
m = alpha_m/(alpha_m+beta_m) : 1
alpha_m = 0.1/mV*10*mV/exprel(-(v+35*mV)/(10*mV))/ms : Hz
beta_m = 4*exp(-(v+60*mV)/(18*mV))/ms : Hz

dh/dt = 5*(alpha_h*(1-h)-beta_h*h) : 1
alpha_h = 0.07*exp(-(v+58*mV)/(20*mV))/ms : Hz
beta_h = 1./(exp(-0.1/mV*(v+28*mV))+1)/ms : Hz

dn/dt = 5*(alpha_n*(1-n)-beta_n*n) : 1
alpha_n = 0.01/mV*10*mV/exprel(-(v+34*mV)/(10*mV))/ms : Hz
beta_n = 0.125*exp(-(v+44*mV)/(80*mV))/ms : Hz

# synaptic conductances
dge/dt = -ge*(1./taue) : siemens
dgi/dt = -gi*(1./taui) : siemens

# external input parameters
Iamp0 : amp   # constant or pulsed current
Iamp : amp    # sinusoidal amplitude
Ifreq : Hz    # sinusoidal frequency
'''

# Neuron group
N = 100
group = NeuronGroup(N, eqs, threshold='v>-20*mV', refractory=3*ms,
                    method='rk4')
group.v = -70*mV
group.h = 1


######### noise goes everywhere ##########
# Noise input
rate = 20*Hz
Noise_in = PoissonGroup(N, rates=rate)
S = Synapses(Noise_in, group, on_pre='ge += we')
S.connect(j='i')

# Example: set input currents in n_a neurons

n_a_arr = np.zeros(N)
# first n_a percentage becomes 1
n_a = n_a * N  # convert percentage to number of neurons
n_a_arr[:int(n_a)] = 1  # set first n_a neurons to 1  

# sinusoidal input
n_b_arr = np.zeros(N)
#last n_b percentage becomes 1
n_b = n_b * N  # convert percentage to number of neurons
n_b_arr[-int(n_b):] = 1  # set last n_b neurons to 1

n_ab_arr = n_a_arr + n_b_arr  # overlap neurons
#make 0s as 1s and 1s as 0s
n_ab_arr[n_ab_arr == 0] = 1
n_ab_arr[n_ab_arr == 1] = 0

# Monitors
monitor = StateMonitor(group, 'v', record=True)
monitor_s = SpikeMonitor(group)
# syn_rate = PopulationRateMonitor(Noise_in)
# monitor_rate = PopulationRateMonitor(group)
# ge_monitor = StateMonitor(group, 'ge', record=True)
# input_monitor = SpikeMonitor(Noise_in)


n_sims = 100

for i in tqdm(range(n_sims)):    
    group.Iamp0 = 0*0.55*uA * (n_a_arr + n_ab_arr)   # baseline or pulse amplitude
    group.Iamp  = 0*0.35*uA * (n_b_arr + n_ab_arr)  # sinusoidal amplitude
    group.Ifreq = 8*Hz     # sinusoidal frequency

    run(duration/5.0)

    group.Iamp0 = 0.55*uA * (n_a_arr + n_ab_arr)   # baseline or pulse amplitude
    group.Iamp  = 0.35*uA * (n_b_arr + n_ab_arr)  # sinusoidal amplitude
    group.Ifreq = 8*Hz     # sinusoidal frequency

    run(3*duration/5.0)

    group.Iamp0 = 0*0.55*uA * (n_a_arr + n_ab_arr)   # baseline or pulse amplitude
    group.Iamp  = 0*0.35*uA * (n_b_arr + n_ab_arr)  # sinusoidal amplitude
    group.Ifreq = 8*Hz     # sinusoidal frequency   

    run(duration/5.0)
# device.build()
print('Whole sim completed')


#plot the raster plot for neuron group
plt.figure(figsize=(10, 6), dpi=150)
plt.plot(monitor_s.t/ms, monitor_s.i, 'k|', markersize= 3)
plt.xlabel('Time (ms)', fontsize=12)
plt.ylabel('Neuron index', fontsize=12)
plt.title('Raster plot of neuron spikes', fontsize=14)
plt.tight_layout()
plt.xlim(0, 2*(duration/ms))
plt.show() 

```

rastor
rastor1500Ă—900 21.3 KB

Actual output (if relevant)

i am uploading output.txt for a neat post.

output.txt (918.7 KB)

Full traceback of error (if relevant)

I am uploading the brain2 log files

brian_stderr_shwjw3wh.log (125 Bytes)

brian_stdout_17isr1mk.log (916.9 KB)

2

Hi @Jeevan_codes8580 ,
In cpp_standalone mode, the run() commands don’t execute immediately; they add instructions to a queue. The device.build() at the end takes this entire queue (300 run calls in your case) and builds a single C++ project to execute everything in one go. The SpikeMonitor and StateMonitor are therefore trying to allocate memory for all 400 seconds of simulation time, which is likely causing your computer to run out of memory and crash. The reason it works for 10 simulations (a 40s total duration) but not 100 is that the memory requirement is 10 times smaller.

To fix this, you need to run truly independent trials. Since the store/restore mechanism is not available in C++ standalone, the correct approach is to manually reset the state of your simulation before each trial. This allows you to compile the C++ code just once and then re-run it with a clean slate for each trial.

WHat you’ll need to do is :

  1. Single Compilation: Create your core simulation objects (NeuronGroup, Synapses) once, outside the main loop. Brian2 will compile the code for these objects on the first run() call.
  2. Manual State Reset: Inside your for loop, at the beginning of each trial, explicitly reset all the state variables (group.v, group.h, group.ge, etc.) back to their initial values.
  3. Fresh Monitors: To get clean recordings for each trial, create new SpikeMonitor and StateMonitor objects inside the loop.

code would look like :

#define overlap percentage
from IPython.display import clear_output
import numpy as np
import matplotlib.pyplot as plt
import os
import shutil
from brian2 import *
from tqdm import tqdm

# --- 1. SET UP THE DEVICE AND BUILD DIRECTORY ---
set_device('cpp_standalone', build_on_run=True)
prefs.devices.cpp_standalone.openmp_threads = 20

if os.path.exists('output'):
    shutil.rmtree('output')


# --- 2. DEFINE PARAMETERS AND EQUATIONS (outside the loop) ---
start_scope()
defaultclock.dt = 0.05*ms
duration = 4*second
n_ab = 0.025
n_a = (1 - n_ab)/2
n_b = (1 - n_ab)/2

Cm = 1*uF
gL = 0.1*msiemens
EL = -65*mV
ENa = 55*mV
EK = -90*mV
gNa = 35*msiemens
gK = 9*msiemens
taue = 5*ms
taui = 10*ms
Ee = 0*mV
Ei = -80*mV
we = 15*uS
wi = 67*uS

# Hodgkin-Huxley + synaptic + external currents
eqs = '''
dv/dt = (-gNa*m**3*h*(v-ENa) - gK*n**4*(v-EK) - gL*(v-EL)
         + Iamp0 + Iamp * sin(2*pi*Ifreq*t)
         + ge*(Ee-v) + gi*(Ei-v)) / Cm : volt
m = alpha_m/(alpha_m+beta_m) : 1
alpha_m = 0.1/mV*10*mV/exprel(-(v+35*mV)/(10*mV))/ms : Hz
beta_m = 4*exp(-(v+60*mV)/(18*mV))/ms : Hz
dh/dt = 5*(alpha_h*(1-h)-beta_h*h) : 1
alpha_h = 0.07*exp(-(v+58*mV)/(20*mV))/ms : Hz
beta_h = 1./(exp(-0.1/mV*(v+28*mV))+1)/ms : Hz
dn/dt = 5*(alpha_n*(1-n)-beta_n*n) : 1
alpha_n = 0.01/mV*10*mV/exprel(-(v+34*mV)/(10*mV))/ms : Hz
beta_n = 0.125*exp(-(v+44*mV)/(80*mV))/ms : Hz
dge/dt = -ge*(1./taue) : siemens
dgi/dt = -gi*(1./taui) : siemens
Iamp0 : amp
Iamp : amp
Ifreq : Hz
'''

# --- 3. CREATE THE CORE SIMULATION OBJECTS (outside the loop) ---
N = 100
group = NeuronGroup(N, eqs, threshold='v>-20*mV', refractory=3*ms, method='rk4')
Noise_in = PoissonGroup(N, rates=20*Hz)
S = Synapses(Noise_in, group, on_pre='ge += we')
S.connect(j='i')


n_a_arr = np.zeros(N)
n_a_arr[:int(n_a * N)] = 1
n_b_arr = np.zeros(N)
n_b_arr[-int(n_b * N):] = 1
n_ab_arr = np.zeros(N)
n_ab_arr[int(n_a * N):-int(n_b * N)] = 1

# --- 4. STORE INITIAL VALUES (outside the loop) ---
v_init = -70*mV
h_init = 1.0

# --- 5. RUN MULTIPLE INDEPENDENT TRIALS ---
n_sims = 100
all_spike_times = []
all_spike_indices = []


for i in tqdm(range(n_sims)):
    # ** Manually reset all state variables to their initial conditions **
    group.v = v_init
    group.h = h_init
    group.ge = 0 * siemens
    group.gi = 0 * siemens
    # 'n' will start at its steady-state value, which is fine.

    # ** Create new monitors for each run to get clean data **
    monitor_s = SpikeMonitor(group)

    # Now what you do is create a new Network for this trial, including the new monitor
    net = Network(group, Noise_in, S, monitor_s)

    # Set the current injection patterns for this trial
    group.Iamp0 = 0*uA
    group.Iamp  = 0*uA
    group.Ifreq = 8*Hz
    net.run(duration/5.0)

    group.Iamp0 = 0.55*uA * (n_a_arr + n_ab_arr)
    group.Iamp  = 0.35*uA * (n_b_arr + n_ab_arr)
    net.run(3*duration/5.0)

    group.Iamp0 = 0*uA
    group.Iamp  = 0*uA
    net.run(duration/5.0)

    # Store the results from this trial's monitor
    all_spike_times.append(monitor_s.t)
    all_spike_indices.append(monitor_s.i)


# --- plotting ---
plt.figure(figsize=(10, 6), dpi=150)
plt.plot(all_spike_times[-1]/ms, all_spike_indices[-1], 'k|', markersize=3)
plt.xlabel('Time (ms)', fontsize=12)
plt.ylabel('Neuron index', fontsize=12)
plt.title('Raster plot of neuron spikes (last trial)', fontsize=14)
plt.tight_layout()
plt.xlim(0, duration/ms)
plt.show()

Let me know if this worked :smile:

3

Hi. Thanks for your reply :smiley: . I think I understand what you are trying to do. I do not know about net.run (). I will try to read its usage. However, as of now the code you provided shows me some error :smiling_face_with_tear: . I will paste it in a .txt file and attach it. Let me try this approach of running a particular network and saving it’s output multiple times in a list. I was trying the same with Joblib but it freezes my systems more often than not.

output.txt (1.6 KB)

4

ahh .. shoot , let me look more into this … the problem is I cannot try it out nicely locally as openmp does not play well on macOS for me `:(

5

Hi all. @Mrigesh_Thakur’s diagnosis is basically correct: the original code creates a huge simulation, since the run calls only add things to the “code to be generated”. On my machine, I do not even get to running it, since the compilation is trying to compile too many files in parallel. Actually, in your error log it even mentions that it cannot compile things at all, since the command line (which states all the filenames) is too long. Do you need one long simulation, or would it also work for you to start each simulation (not necessarily each run, it could be e.g. the three runs within your loop body together) from scratch? If the latter is an option, you can use the run_args feature which only compiles the network once, and then runs it with different parameters: Computational methods and efficiency — Brian 2 0.0.post128 documentation