Ok, understood. But note (as @rth mentioned above), that you can manually store and restore variable values (membrane potential, weights, etc.), if you do not use synaptic delays.
Hi everyone , just to let you know that the “change parameters between standalone runs” feature has now been merged into master and is ready for more extensive testing. It will be a few months before the next release, so please understand that some details of the feature might still change, but these changes should be minor (unless someone finds a major issue, of course). To test, you can install Brian’s latest development version as explained in the installation docs (the version that introduced the feature is 2.5.4.post54
), or of course run Brian directly from its git repository. There is quite a bit of documentation for the new feature, and I added/updated examples to use it as well (make sure you are on the “latest” version of the documentation, the “stable” version won’t show it). Please note that – thanks to the suggestions in this discussion thread – this feature not only works for neuronal variables, but also for synapses and TimedArray
.
A few links:
- user documentation for the feature
- developer documentation for the feature
- Bisection example from Stimberg et al. 2019, using the feature
- Adapted example for efficient multiple runs of standalone
- Using standalone with multiple processes
Let me know what you think and whether it works for you!
Hello! I have been struggling adapting my code to the example that you link to above (Using standalone with multiple processes) but I cannot get it to work. It keeps complaining with AttributeError: Can't pickle local object
.
Then I went back to your example and tried running it and I found that it also gave the same error. I also get the error if I try to run it directly using the launch-binder
button that you have on the website. How can I get that code to work?
Thanks!
Albert
Hi @acompte. Could it be that you are running the latest released version (2.5.4) instead of a development version (>=2.5.4.post54)? Unfortunately, we still did not do a release since my announcement. The binder link in the “latest” documentation will also use the released version instead of the development version, so unfortunately it isn’t functional for examples that rely on unreleased features.
Thank you for your quick reply. I am running the latest development version (2.5.4.post111) of Brian2. I wanted to try it on the 2.5.4.post54 version but this command
python -m pip install git+https://github.com/brian-team/brian2.git@2.5.4.post54
is giving the error:
pathspec ‘2.5.4.post54’ did not match any file(s) known to git
Hi again, I finally managed to install development version 2.5.4.post54 and it runs the code as in the documentation. If I switch to my environment with the 2.5.4.post111 Brian2 version, however, the code does not run.
You could install the post54
version from test.pypi.org
, to install it from github would use a slightly different syntax (using a commit hash instead of a version number). But both versions should work. Sorry to insist on this, but are you sure you are using the .post111
version when you run your code (e.g. can you check that the output of import brian2; print(brian2.__version__)
is 2.5.4.post111
)?
When launching the example with its launch-binder
button, I can successfully run it after installing the development version with:
%pip install --upgrade --pre -i https://test.pypi.org/simple/ Brian2
Sorry, I did not see your new post before answering: could it be that there is something wrong with the installation in your .post111
environment, e.g. do you have two versions of Brian 2 installed? Try adding
import brian2
print(brian2.__version__)
to your script to make sure it uses the correct version.
Yes, you are right. My environment lists 2.5.4.post111 but from within python the version reported is 2.5.4
Sorry about this mess. I will clean up my python environment
This happened to me more than once, so I’m in no position to complain Hope this makes it work for you!
Dear @mstimberg . I am still stuck with my standalone multiple-simulation code. I managed to get your example to work using the right Brian2 version, but when I write my code following that recipe I do not get it to run with multiprocessing. It complains about KeyError: ‘Stim’ or AttributeError: Can’t pickle local object ‘check_units..do_check_units..new_f’
If I remove the multiprocessing loop and instead I do a regular loop so that simulations are run in series (see last commented section in the code), everything works fine. Below is my code, which really follows the structure of your example quite literally. I am at a loss with this. Any help will be most appreciated.
from brian2 import *
import sys
import time
import numpy as np
import os
# import brian2cuda
# set_device("cuda_standalone")
set_device("cpp_standalone", build_on_run=False)
def decode(firing_rate, N_e):
R = np.sum(np.dot(firing_rate,np.exp(1j*neurspaceE)))/(2*N_e)
return np.angle(R)
class SimWrapper:
def __init__(self):
self.net = Network()
## Input currents
#equations for each neuron
eqs = '''
dV/dt = (gea + gen + gi - V + Ix + Inp*Stim(t))/tau : volt (unless refractory)
dgea/dt = -gea/taua : volt
dgen/dt = -gen/taun : volt
dgi/dt = -gi/taug : volt
tau : second
Ix : volt
Inp : volt
theta : 1
stimat : 1 (shared, constant)
spike_count : integer
'''
eqsSTP = '''
dx/dt=(1-x)/taud : 1 (event-driven)
du/dt=(U-u)/tauf : 1 (event-driven)
wA : volt
wN : volt
'''
onpreSTP = '''
gea_post += u*x*wA
gen_post += u*x*wN
x = x*(1 - u)
u = u + U*(1 - u)
'''
ring = NeuronGroup(NE+NI,model=eqs,
threshold='V > 20*mV',
reset='''V = -3.33*mV
spike_count+=1''',
refractory=refE, method='exact', name='ring')
ring.run_regularly('spike_count = 0',
dt=500*ms)
networkE = ring[:NE]
networkE.theta = neurspaceE
networkE.tau=tauE
networkE.Ix=1.66*sqrt(KE)*mV
networkE.V = Vr + rand(NE)*(Vt - Vr)
networkI= ring[NE:]
networkI.theta = neurspaceI
networkI.tau=tauI
networkI.Ix=1.85*0.83*sqrt(KE)*mV
networkI.V = Vr + rand(NI)*(Vt - Vr)
C1=Synapses(networkE, networkE, eqsSTP, on_pre=onpreSTP, delay=1.*ms, name='EE')
C3=Synapses(networkE, networkI, '''wA : volt
wN : volt''',
on_pre = '''gea +=wA
gen +=wN''', name='EI')
C5=Synapses(networkI, networkE, 'w: volt', on_pre='gi +=w', name='IE')
C6=Synapses(networkI, networkI, 'w: volt', on_pre='gi +=w', name='II')
fE = float(KE)/float(NE)
fI = float(KI)/float(NI)
C1.connect(p='fE*exp(-0.5*(arccos(cos(theta_pre - theta_post))/(sigEE*2*pi))**2)/sqrt(2*pi)/sigEE')
C1.wA = gEEA
C1.wN = gEEN
C1.x = 1
C1.u = U
C3.connect(p='fE*exp(-0.5*(arccos(cos(theta_pre - theta_post))/(sigEI*2*pi))**2)/sqrt(2*pi)/sigEI')
C3.wA = gEIA
C3.wN = gEIN
C5.connect(p='fI*exp(-0.5*(arccos(cos(theta_pre - theta_post))/(sigIE*2*pi))**2)/sqrt(2*pi)/sigIE')
C5.w = gIE
C6.connect(p='fI*exp(-0.5*(arccos(cos(theta_pre - theta_post))/(sigII*2*pi))**2)/sqrt(2*pi)/sigII')
C6.w = gII
self.device = get_device()
self.device.apply_run_args()
networkI.Inp = 0
networkE.Inp = 'stimE*exp(-0.5*(arccos(cos(theta - stimat))/(epsE*2*pi))**2)'
self.net.add([ring, C1, C3, C5, C6])
self.net.run(runtime)
self.device.build(run=False, directory=None)
def do_run(self, theta_i):
from brian2.devices import device_module
device_module.active_device = self.device
result_dir = f'result_{theta_i}'
self.device.run(run_args={self.net['ring'].stimat: theta_i},
results_directory=result_dir)
dec = decode(self.net['ring'].spike_count[:NE], NE)
return ([theta_i, dec])
if __name__ == "__main__":
start_time = time.time()
# number of neurons and connections
N=20000 # total number of neurons
K=500 # total number of inputs
NE=int(ceil(0.8*N)) # number of excitatory neurons
NI=int(floor(0.2*N)) # number of inhibitory neurons
KE=int(ceil(0.8*K)) # number of excitatory inputs
KI=int(floor(0.2*K)) # number of inhibitory inputs
# neuron labels
neurspaceE = np.linspace(-np.pi, np.pi, NE, endpoint=False)
neurspaceI = np.linspace(-np.pi, np.pi, NI, endpoint=False)
# general simulation
nstims=1 #50 # number of simulations in the series
stim_on=2000*ms # onset of stimulus
stim_off=3000*ms # offset of stimulus
runtime=5*second #7*second # total time of each simulation run
# stimulus input parameters
dtm = 10*ms
ntms = int(ceil(runtime/dtm))
values = np.ones(ntms)
values[:int(stim_on/dtm)] *= 0
values[int(stim_off/dtm):] *= 0
Stim = TimedArray(values, dt=dtm)
stimE=0.24*mV # stimulus input
epsE=60
# cellular and synaptic parameters
tauE=20*ms
tauI=10*ms
taua = 3*ms # AMPA synapse decay
taun = 50*ms # NMDA synapse decay
taug = 4*ms # GABA synapse decay
Vt = 20*mV # spike threshold ##CHANGE THEM IN EQS
Vr = -3.33*mV # reset value
refE= 0*ms # refractory periods
refI= 0*ms # refractory periods
# connectivity parameters
gEEA=533.3*mV*ms
gEEN=0.92*533.3*mV*ms
gEIA=67.2*mV*ms
gEIN=1*7.4*mV*ms
gIE=-138.6*mV*ms
gII=-90.6*mV*ms
sigmaEE=30 # E-to-E footprint in degrees
sigmaEI=35 #70 # E-to-E footprint in degrees
sigmaIE=30 #60 # E-to-E footprint in degrees
sigmaII=30 #60 # E-to-E footprint in degrees
# parameters of STP
taud = 200*ms
tauf = 450*ms
U=0.03
#parameter scalings
sigEE=sigmaEE/360.0
sigEI=sigmaEI/360.0
sigIE=sigmaIE/360.0
sigII=sigmaII/360.0
epsE = epsE/360.0
gEEA=gEEA/sqrt(KE)/taua
gEEN=gEEN/sqrt(KE)/taun
gEIA=gEIA/sqrt(KE)/taua
gEIN=gEIN/sqrt(KE)/taun
gIE=gIE/sqrt(KI)/taug
gII=gII/sqrt(KI)/taug
stimE=stimE*sqrt(KE)
## random stimulus location in each simulation run
stimuli = np.random.randint(low=0, high=32, size=nstims)/32.*2*np.pi - np.pi
sim = SimWrapper()
## run the simulations in parallel
from multiprocessing import Pool
with Pool(nstims) as pool:
output = pool.map(sim.do_run, stimuli)
## run the simulations in series
# output = []
# for stim in stimuli:
# out = sim.do_run(stim)
# output.append(out)
print(f"Done in {time.time() - start_time}")
np.savetxt("ringattractor_series_results.csv", output, delimiter=",")
Dear @acompte. Apologies for this, from a cursory look I’d say it is because we don’t handle “non-standard” functions (in your case, a TimedArray
) well in a multiprocessing context. The underlying issue is that Python uses its pickle
framework to convert the network into a single byte string that is sent to the individual processes and unwrapped again into a network there. This mechanism works reasonably well for simple objects, but “local” functions (functions not defined in a module) are often problematic. I think we can fix this on Brian’s side, but until then, you can try replacing the use of Python’s built-in multiprocessing
module by pathos
(pathos · PyPI). This package is much more powerful in general (e.g. it could launch the processes on separate machines via ssh
), but in particular uses the dill
(dill · PyPI) package for pickling, which handles things like local functions much better. After installing pathos
, you’d only need to change the multiprocessing
lines into the very similar:
from pathos.multiprocessing import ProcessingPool
with ProcessingPool(nstims) as pool:
output = pool.map(sim.do_run, stimuli)
Hopefully this fixes your issue, but as I said, I’ll also have a look into fixing this from our side (otherwise we should probably make our example use pathos
in the first place…).
Thank you @mstimberg . I just tried pathos
and it also seems not to like the TimedArray. Instead, I removed that TimedArray object from my code and I defined the input within the equations of the NeuronGroup as
Stim = int(t > stim_on) * int(stim_off > t) : 1 (shared, constant over dt)
and now it apparently works well.
Finally, no chance I can get this to work with cuda_standalone
, right?
Thank you for your help!
Albert
That’s odd, I could run your code with it on my machine, but these things can also depend on OS and other things. Glad that it works with the equation-based approach for you (which is probably preferred in this case, anyway).
Currently, cuda_standalone
does not work with the new run_args
argument, but we hope to make it work soon after the new Brian release. This feature (running simulations without re-compilation) would make a lot of sense for CUDA, where compilation times are very long. Running things in parallel with multiprocessing on the GPU is much trickier, though, and I am not sure that it would be useful in many cases. For small networks (small enough so that you can fit several of them on the GPU at once), the approach started here by @denisalevi’s student seems more fruitful: GitHub - denisalevi/brian2-network-multiplier: Run multiple versions of a `brian2.network` in a single simulation.
Dear @mstimberg , the code is now working nice in cpp_standalone, but I was just wondering if this implementation makes each run of the simulation use the same connectivity implementation, or if the probabilistic connection matrix is generated anew in each new call to do_run
. If so, how could I make all simulations share the same connectivity structure?
Dear @acompte. By default, all random values (including connectivity) will be independent for each run. If you want to have the same values each time, you can set a seed like in this example:
from brian2 import *
set_device('cpp_standalone', build_on_run=False)
class SimWrapper:
def __init__(self):
seed(123) # Same connections every time
self.net = Network()
G = NeuronGroup(5, "")
S = Synapses(G, G, "")
S.connect(p=0.5)
device.build(run=False, directory=None)
self.net.add([S, G])
self.net.run(1 * ms)
self.device = get_device()
self.device.build(run=False, directory=None) # compile the code, but don't run it yet
def do_run(self, tau_i):
# Workaround to set the device globally in this context
from brian2.devices import device_module
device_module.active_device = self.device
result_dir = f'result_{tau_i}'
self.device.run(results_directory=result_dir)
return list(zip(self.net["synapses"].i[:], self.net["synapses"].j[:]))
if __name__ == "__main__":
sim = SimWrapper()
import multiprocessing
with multiprocessing.Pool(5) as pool:
results = pool.map(sim.do_run, range(5))
for idx, result in enumerate(results):
print(idx, ":", result)
If you remove the seed, you will see that it generates a different connectivity matrix each time.
Cf. Random numbers — Brian 2 2.5.4 documentation
For anyone interested in the feature described in this thread, please note that it is included in the new Brian 2.6 release