Running multiple standalone simulations [Request for Feedback]

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 :wave:, 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:

Let me know what you think and whether it works for you!

1 Like

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 :wink: 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 :tada: