Ignore this, just a few imports needed for the tutorial.

In [ ]:
from brian import *
import time
import pickle
import joblib
import gc

def brianreset():
    defaultclock.t = 0*second

Python data structures

Python has some great datatypes, it's worth understanding and using them.

  • List
    • Appending is fast, inserting anywhere else is slow
    • Indexing x[i] is fast
    • Uses a 'dynamic array' structure
  • Dictionary
    • X[key] = value
    • Key can be (almost) anything
    • Fast!
  • Set
    • Fast to add or subtract elements
    • Fast for testing if an item is in the set or not

Python iteration

Don't do this:

In [24]:
for i in range(1000):
    for j in range(1000):
10 loops, best of 3: 42.6 ms per loop

Do this instead:

In [25]:
for i in xrange(1000):
    for j in xrange(1000):
10 loops, best of 3: 38.3 ms per loop

Similarly, don't do this:

In [33]:
A = rand(100000)
B = rand(100000)
In [34]:
for a, b in zip(A, B):
10 loops, best of 3: 63.9 ms per loop

Do this instead:

In [36]:
import itertools
In [35]:
for a, b in itertools.izip(A, B):
10 loops, best of 3: 51.5 ms per loop
  • range creates a list, whereas xrange can only be used for iteration.
  • There are many other nice tools for iteration in the itertools module of Python
  • Best of all is to use numpy vectorisation instead of loops

Saving and loading data in Python

Saving and loading is easy in Python with pickle

In [26]:
import pickle
obj = ['a', 'list', 'of', 'words']
pickle.dump(obj, open('alistofwords.pkl', 'wb'), -1)
newobj = pickle.load(open('alistofwords.pkl', 'rb'))
print ' '.join(newobj)
a list of words

Shelves are a dict-like database

In [28]:
import shelve
shelf = shelve.open('test.shelf')
shelf['meaning'] = 42
newshelf = shelve.open('test.shelf')
print newshelf['meaning']

There are many other databases in Python

A Python example

Compute firing rates for a two parameter family of experimental data recorded in random order with an unknown number of trials per condition.

First we generate some fake data - ignore this.

In [40]:
from collections import namedtuple
paramA = linspace(0, 1, 10)
paramB = linspace(0, 1, 10)
trials = []
Trial = namedtuple("Trial", ['a', 'b', 'spikes'])
for i in xrange(100):
    a = paramA[randint(len(paramA))]
    b = paramB[randint(len(paramB))]
    r = (a-0.5)**2+(b-0.5)**2
    spikes = rand(int(r*100))
    trials.append(Trial(a=a, b=b, spikes=spikes))

Now we compute the firing rates and plot them. Using a dictionary means that we don't need to precompute the set of which parameters exist and which don't, etc., but just immediately start computing what we need.

In [42]:
from collections import defaultdict

counts = defaultdict(int)
numtrials = defaultdict(int)

for trial in trials:
    counts[trial.a, trial.b] += len(trial.spikes)
    numtrials[trial.a, trial.b] += 1
rates = dict()

for a, b in counts.keys():
    rate = counts[a,b]*1.0/numtrials[a,b]
    plot([a], [b], 'o', ms=rate, c='k')
    rates[a, b] = rates

pickle.dump(rates, open('rates.pkl', 'wb'), -1)

Python: final thoughts

There are several other useful toolboxes that are worth finding out about.

  • Multiprocessing for running simulations simultaneously on multiple CPUs or computers (including over the internet). There are also other, simpler solutions for this.
  • GUI toolboxes are surprisingly easy to use and can be useful for interactively exploring your data/simulations.

Vectorising with numpy

Don't do this:

In [47]:
a = linspace(0, 1, 100)
b = linspace(0, 1, 100)
In [54]:
def getimg():
    img = zeros((100, 100))
    for i in xrange(100):
        for j in xrange(100):
            img[i, j] = (a[i]-0.5)**2+(b[j]-0.5)**2
    return img

Do this:

In [56]:
def getimg_vectorised():
    return (a[:, newaxis]-0.5)**2+(b[newaxis, :]-0.5)**2

They're equivalent:

In [58]:
<matplotlib.image.AxesImage at 0x4b8f590>

And the vectorised version is much faster!

In [59]:
%timeit getimg()
%timeit getimg_vectorised()
10 loops, best of 3: 39.3 ms per loop
10000 loops, best of 3: 40.9 us per loop

Advanced vectorisation

With a bit of ingenuity, almost everything can be vectorised, and numpy includes some nice functions to help with this. For the cases that can't, we'll come back to that later.

Here's a more complicated example, computing the coincidence factor between two spike trains.

In [61]:
delta = 0.005
train1 = sort(rand(50))
train2 = sort(rand(55))

Don't do this:

In [62]:
def cf_loop():
    num_coinc = 0
    for t1 in train1:
        for t2 in train2:
            if abs(t1-t2)<delta:
                num_coinc += 1
    return num_coinc

Do this:

In [63]:
def cf_vectorised():
    ends = 0.5*(train2[1:]+train2[:-1])
    i = digitize(train1, ends)
    return sum(abs(train1-train2[i])<delta)

As before, they're equivalent but the vectorised version is an order of magnitude faster.

In [64]:
print cf_loop(), cf_vectorised()
19 19
In [65]:
%timeit cf_loop()
%timeit cf_vectorised()
100 loops, best of 3: 2.22 ms per loop
100000 loops, best of 3: 18.5 us per loop

It can be a little difficult to understand the vectorised algorithms sometimes. Talk to me later if you want to know how this one works.

When vectorisation doesn't work

You can run C++ code inline in various simple ways:

Scipy has a package weave which allows inline C++ code:

In [5]:
from scipy import weave
N = 100
x = zeros(N)
code = '''
for(int i=0; i<N; i++)
    x[i] = i;
namespace = {'x': x, 'N': N}
weave.inline(code, namespace.keys(), namespace, compiler='gcc')
print x[3]

There are also other options:

  • Cython allows you to write code in a Python-like syntax with a few extra keywords to define types, and then generates C++ code.
  • Swig makes it straightforward to write a separate C++ library and link it to Python
  • numba is a very new option, not yet complete, by the people who designed numpy. It allows you to write 'naive' Python code with loops, and it JIT compiles to LLVM and then machine code behind the scenes. In principle, this may turn out to be the best solution, and something we are considering for Brian 2.

Numpy, scipy and other packages

Don't forget to investigate features in numpy, scipy and other third party Python packages. Check out the list of related software on the scipy website. No need to reinvent the wheel!


One of my favourite packages is joblib which I use in every project for cacheing computations to disk rather than managing this manually.

In [69]:
import joblib

cache_mem = joblib.Memory(cachedir='.', verbose=0)
cache_mem.clear() # don't do this normally, it's just for the example

def f(x):
    print "--- I am currently computing x*x for x =", str(x)
    return x*x

print "Computing f(2.3)"
print f(2.3)
print "Computing f(3.5)"
print f(3.5)
print "Computing f(2.3) again"
print f(2.3)
WARNING:root:[Memory(cachedir='.\\joblib')]: Flushing completely the cache
Computing f(2.3)
--- I am currently computing x*x for x = 2.3
Computing f(3.5)
--- I am currently computing x*x for x = 3.5
Computing f(2.3) again

The same thing works with arbitrarily complicated arguments including lists, dicts, complex objects, etc.

You can use this to write your code as if everything was being computed from scratch each time, but on subsequent runs it will complete fast.

Efficient Brian code

There are several key elements to writing efficient Brian code:

  • Using the C++ extensions
  • Using code generation
  • Knowing when to use Connection and Synapses
  • Handling looping over several simulation runs correctly

Brian's C++ extensions

Brian is designed to run in pure Python. However, if you do have a C++ compiler (if you are running on Linux, or Windows with MinGW installed) then you can build several extensions which replace their pure Python equivalents. These can make a big difference to speed. The latest release of Brian (1.4.1) will try to detect if you have a compiler installed and build these automatically, but if it fails or if you're running an older version of Brian you can do it manually (see the docs for details).

  • ccircular package
  • cspikequeue package

Using code generation

Brian has support for automatically generating C++ code in the background, compiling it and then using this rather than working in pure Python. Sometimes this can make a big speed difference.

Brian has a preference mechanism, you can create a file brian_global_config.py anywhere on your Python path, and set preference values in it. The following is the recommended settings for most computers if you have gcc installed (on linux or Windows):

In [1]:
from brian.globalprefs import *
    gcc_options=['-ffast-math', '-march=native'],
C:\Users\goodmad\programming\brian\brian\synapses\spikequeue.py:490: UserWarning: Using C++ SpikeQueue
  warnings.warn('Using C++ SpikeQueue')

In Brian 1.x code generation doesn't work for all models, so you might need to disable it for some simulations (which you can do in the script rather than the preferences).

Here's an example of the speed difference:

In [2]:
from brian import *
import time
def runsim(N, duration=1*second, codegen=False):
    clk = Clock(dt=0.1*ms)
    G = NeuronGroup(N, 'dv/dt=(2-v**2)/(10*ms):1', reset=0, threshold=1, clock=clk)
    G.v = rand(len(G))
    net = Network(G)
    start = time.time()
    end = time.time()
    return end-start
In [3]:
N = 10
runsim(N, codegen=True, duration=1*ms)
t_without = runsim(N, codegen=False)
t_with = runsim(N, codegen=True)
print 'N =', N
print 'Without code generation:', t_without 
print 'With code generation:   ', t_with
print 'Speedup:                ', t_without/t_with
brian.stateupdater: WARNING  Using codegen CStateUpdater
brian.stateupdater: WARNING  Using codegen CStateUpdater
N = 10
Without code generation: 0.461999893188
With code generation:    0.175000190735
Speedup:                 2.63999651228

Connection and Synapses

Brian 1.4 introduced the Synapses object, which is much more flexible than the older Connection object.

However, Synapses is currently much slower than Connection (we're working on this for Brian 2).

Conclusion: only use Synapses if you cannot do what you need with Connection (e.g. nonlinear synapses).

Working with large Connections

Some users have run out of memory when working with very large numbers of synapses in Connection objects. This is because:

  • The underlying matrix object uses two data structures and converting between them is expensive.
  • When building the network we use a data structure optimised for inserting elements.
  • When running the simulation we use a data structure optimised for fast extraction of rows and columns.
  • When converting between them a large amount of memory is used.

To solve this problem, break your Connection into multiple sub-connections:

  • Instead of one connection from neurons [0, ..., N-1] to [0, ..., N-1]
  • Use two connections:
    • from [0, ..., N/2-1] to [0, ..., N-1]
    • from [N/2-1, ..., N-1] to [0, ..., N-1]

Looping over simulation runs

Here's how to do a loop over several parameters:

In [4]:
import gc
for params in []:
    clear(True, True)
    defaultclock.t = 0*second
    # your code here


  • Each time you call run() it collects all existing Brian objects
  • Because Python doesn't always immediately deallocate objects, sometimes objects from the previous run still exist for the next run
  • This means the second time you call run() you might be doing twice as much work as the first time.
  • Calling clear(True, True) deletes all memory from currently existing Brian objects.
  • Calling gc.collect() forces Python to deallocate these objects.
  • Don't forget to reset the global clock time to 0!

All of this will work much better in Brian 2, the above is a workaround for Brian 1.x. For now we use the following function to abbreviate this:

In [9]:
def brianreset():
    defaultclock.t = 0*second

Another way to run a simulation for several parameters is in parallel...

Running simulations in parallel

There are many options. Most of them rely on you being able to pickle your Brian objects. This usually works, but sometimes objects cannot be pickled and you need to find workarounds. For example f = lambda x: x**2 is a definition of the function f(x)=x**2, but lambda functions cannot be pickled. Simply use a def instead of lambda.

We'll use this I-F curve function as an example:

In [14]:
def IF(v0):
    eqs = '''
    dv/dt=(v0-v)/(10*ms) : volt
    group = NeuronGroup(1, eqs, threshold='v>10*mV', reset='v = 0*mV')   
    monitor = SpikeMonitor(group)
    duration = 1*second
    return monitor.nspikes/duration

V0 = [v0*mV for v0 in range(20)]

log_level_error() # suppress a few warning messages

Using a loop:

In [17]:
%time plot(V0, [IF(v0) for v0 in V0])
CPU times: user 102.20 s, sys: 0.00 s, total: 102.20 s
Wall time: 102.18 s
[<matplotlib.lines.Line2D at 0x8083130>]

Using multiprocessing (doesn't work in IPython):

In [21]:
if 0:
    import multiprocessing
    pool = multiprocessing.Pool()
    plot(V0, pool.map(IF, V0))

Using IPython:

In [18]:
from IPython import parallel
rc = parallel.Client()
view = rc[:]
In [19]:
from brian import *
import gc
In [20]:
%time plot(V0, view.map_sync(IF, V0))
CPU times: user 19.93 s, sys: 0.00 s, total: 19.93 s
Wall time: 19.92 s
[<matplotlib.lines.Line2D at 0x80bbc90>]

There are also many other options for parallel execution of code, including:

  • playdoh written by Brian's very own Cyrille Rossant. This library also allows you to easily apply global optimisation algorithms running in parallel over multiple CPUs and machines, for example for fitting models to data.
  • Brian's run_tasks function, but we'll come back to that later.

But before you start using multiple CPUs...

Can you just do something like this?

In [25]:
def IF_vectorised(V0):
    eqs = '''
    dv/dt=(v0-v)/(10*ms) : volt
    v0 : volt
    group = NeuronGroup(len(V0), eqs, threshold='v>10*mV', reset='v = 0*mV')   
    group.v0 = V0
    monitor = SpikeCounter(group)
    duration = 1*second
    return monitor.count/duration

%time plot(V0, IF_vectorised(V0))
CPU times: user 5.72 s, sys: 0.00 s, total: 5.72 s
Wall time: 5.72 s
[<matplotlib.lines.Line2D at 0x93c20d0>]

Managing simulation data

Good options include:

  • Using joblib and writing all your code as if it was running from scratch each time, but automatically cache the results.
  • Writing data to a shelf file or database.

Brian has a specific solution based on shelf files, but which allows for the use of multiple CPUs and machines in parallel.

  • DataManager object
    • uses multiple shelf files (one per machine)
    • controls access via a single process (because the databases underlying shelf are not process safe)
  • run_tasks function
    • like a parallel map
    • the data is written into a DataManager object
    • has an optional GUI allowing control over running simulations
In [ ]:
if 0: # doesn't work in IPython
    from brian import *
    from brian.tools.datamanager import *
    from brian.tools.taskfarm import *
    def find_rate(k, report):
        eqs = '''
        dV/dt = (k-V)/(10*ms) : 1
        G = NeuronGroup(1000, eqs, reset=0, threshold=1)
        M = SpikeCounter(G)
        run(30*second, report=report)
        return (k, mean(M.count)/30)
    if __name__=='__main__':
        N = 20
        dataman = DataManager('taskfarmexample')
        if dataman.itemcount()<N:
            M = N-dataman.itemcount()
            run_tasks(dataman, find_rate, rand(M)*19+1)
        X, Y = zip(*dataman.values())
        plot(X, Y, '.')
        ylabel('Firing rate (Hz)')
In [2]:
from IPython.core.display import Image

Understanding Brian's Network

When you call run(), what happens is:

  • Brian creates a Network() object
  • It adds all existing Brian objects to this network
    • Actually, the rule is slightly more complex than this, see the discussion of 'magic' in the docs for more details.
  • It calls Network.run() on this

The way Network.run() works is:

  • It builds an 'update schedule'.
    • This is the order in which objects are updated
    • Normally, this is
      • Numerical integration step
      • Thresholding operation
      • Synaptic propagation
      • Resets
    • You can change this.
    • You can choose to have certain objects update at specific points in this schedule
      • always at the start
      • always at the end
      • always just before synaptic propagation
      • etc.
  • At each time step:
    • Choose a clock to process
      • The one with the lowest current time
      • If several have the same time, they can be ordered for more precise control
    • Process all objects with this clock
    • Increment that clock's time
  • This is repeated until all clocks have run for the specified duration.

Network operations

The NetworkOperation object, and network_operation decorator allow you to run arbitrary code at every timestep. This is one of the main mechanisms for extending Brian to do new things that we didn't think of.

In [26]:

N = 5
G = NeuronGroup(N, 'dV/dt = -V/(10*ms) : 1')
def random_input():
    i = randint(N)
    G.V[i] += 1
M = StateMonitor(G, 'V', record=True)
plot(M.times, M.values.T+3*arange(N)[newaxis, :]);

Some useful Brian features you may not know


In [6]:

N = 5
duration = 100 * ms
Vr = -60 * mV
Vt = -50 * mV
tau = 10 * ms
Rmin = 1 * Mohm
Rmax = 10 * Mohm
freq = 50 * Hz
k = 10 * nA

eqs = '''
dV/dt = (-(V-Vr)+R*I)/tau : volt
R : ohm
I : amp

G = NeuronGroup(N, eqs, reset='V=Vr', threshold='V>Vt')
G.R = linspace(Rmin, Rmax, N)

t = linspace(0 * second, duration, int(duration / defaultclock.dt))
I = clip(k * sin(2 * pi * freq * t), 0, Inf)
G.I = TimedArray(I)

M = MultiStateMonitor(G, record=True)


ylabel('I (amp)')
ylabel('V (volt)')
<matplotlib.text.Text at 0x9ac5630>

Linked variables

In [7]:

N = 5
f = 50 * Hz
a_min = 1.0
a_max = 100.0
tau_haircell = 50 * ms
tau = 10 * ms
duration = 100 * ms

eqs_haircells = '''
input = a*sin(2*pi*f*t) : 1
x = clip(input, 0, Inf)**(1.0/3.0) : 1
a : 1
dy/dt = (x-y)/tau_haircell : 1 

haircells = NeuronGroup(N, eqs_haircells)
haircells.a = linspace(a_min, a_max, N)
M_haircells = MultiStateMonitor(haircells, vars=('input', 'y'), record=True)

eqs_nervefibres = '''
dV/dt = (I-V)/tau : 1
I : 1
nervefibres = NeuronGroup(N, eqs_nervefibres, reset=0, threshold=1)
nervefibres.I = linked_var(haircells, 'y')
M_nervefibres = MultiStateMonitor(nervefibres, record=True)


brian.equations   : WARNING  Equation variable t also exists in the namespace
brian.equations   : WARNING  Equation variable x also exists in the namespace
brian.equations   : WARNING  Equation variable I also exists in the namespace
<matplotlib.text.Text at 0xa005050>

Remote control

Allows us to control a running Brian script.

Run the script with a RemoteControlServer object, like this:

In [8]:

eqs = '''
dV/dt = (I-V)/(10*ms)+0.1*xi*(2/(10*ms))**.5 : 1
I : 1

G = NeuronGroup(3, eqs, reset=0, threshold=1)
M = RecentStateMonitor(G, 'V', duration=50*ms)

server = RemoteControlServer()

brian.equations   : WARNING  Equation variable I also exists in the namespace

Now you can control it like this (we'll run these commands in a separate IPython console):

In [9]:
if 0:
    from brian import *
    client = RemoteControlClient()
    clf(); plot(*client.evaluate('(M.times, M.values)'))
    client.execute('G.I = 1.1')
    clf(); plot(*client.evaluate('(M.times, M.values)'))

Linear, nonlinear and multilinear equations

When you create a NeuronGroup in Brian, it automatically selects what it thinks is the best way to integrate the equations.

  • If the equations are linear, they can be solved exactly.
  • If the equations are nonlinear, a numerical method is used.
  • By default, the numerical method is forward Euler, but alternative methods can be selected.

There is a limitation with linear equations though, which will not be present in Brian 2, namely that the equations have to be identical for all neurons. For example, this will be integrated exactly:

In [10]:
def getsol_exact():
    tau = 10*ms
    eqs = '''
    dv/dt = -v/tau : 1
    G = NeuronGroup(1, eqs)
    G.v = 1
    M = StateMonitor(G, 'v', record=True)
    return M[0]

But this will not:

In [11]:
def getsol_numerical():
    eqs = '''
    dv/dt = -v/tau : 1
    tau : second
    G = NeuronGroup(1, eqs)
    G.v = 1
    G.tau = 10*ms
    M = StateMonitor(G, 'v', record=True)
    return M[0]
In [12]:
exact = getsol_exact()
numerical = getsol_numerical()
[<matplotlib.lines.Line2D at 0x9c85d30>]

The reason is that the second form with tau as a parameter would allow each neuron to be a different linear equation, and so it has to treat it as nonlinear. There is a solution though: brian.experimental.multilinearstateupdater.MultiLinearNeuronGroup (see the docs).

Deep Brian voodoo: custom resets and refractoriness

The system for resets and refractoriness in Brian is slightly limited and this will be corrected in Brian 2. For the moment though, there are some workarounds:

  • SimpleCustomRefractoriness
  • CustomRefractoriness

The way the reset and refractoriness conditions work in Brian is that you can specify an arbitrary reset condition, and during the refractory period this will be called each time step on the neurons that are refractory.

This would not allow you to, for example, combine an adaptive threshold with a refractory period. Here's how to do that:

In [23]:
eqs = '''
dv/dt = -v/(10*ms) : 1
dvt/dt = (1-vt)/(100*ms) : 1
def reset(G, spikes):
    G.vt[spikes] += 1
    G.v[spikes] = 0
inp = PoissonGroup(1, 100*Hz)
refrac = SimpleCustomRefractoriness(reset, 10*ms, 'v')
G = NeuronGroup(1, eqs, threshold='v>vt', reset=refrac)
G.vt = 1
C = Connection(inp, G, 'v', weight=0.5)
M = MultiStateMonitor(G, record=True)
figure(figsize=(10, 8))
ylim(-0.1, 3.0)
brian.experimental.new_c_propagate: WARNING  Using new C based propagation function.
<matplotlib.legend.Legend at 0x9c7ed10>

For even more control, you can choose a separate reset and refractoriness function to be called at the beginning and during the refractory period, respectively:

In [25]:
eqs = '''
dv/dt = -v/(10*ms) : 1
dvt/dt = (1-vt)/(100*ms) : 1
savevt : 1
def resetfunc(G, spikes):
    G.vt[spikes] += 1
    G.savevt[spikes] = G.vt[spikes]
    G.v[spikes] = 0
def refracfunc(G, indices):
    G.vt[indices] = G.savevt[indices]
    G.v[indices] = 0
inp = PoissonGroup(1, 100*Hz)
refrac = CustomRefractoriness(resetfunc, 10*ms, refracfunc)
G = NeuronGroup(1, eqs, threshold='v>vt', reset=refrac)
G.vt = 1
C = Connection(inp, G, 'v', weight=0.5)
M = MultiStateMonitor(G, record=True)
figure(figsize=(10, 8))
ylim(-0.1, 3.0)
brian.experimental.new_c_propagate: WARNING  Using new C based propagation function.
<matplotlib.legend.Legend at 0x974ecf0>

And finally... (if there is any time)

Any questions?