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

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

def brianreset():
clear(True,True)
gc.collect()
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]:
%%timeit
for i in range(1000):
for j in range(1000):
i*j

10 loops, best of 3: 42.6 ms per loop


In [25]:
%%timeit
for i in xrange(1000):
for j in xrange(1000):
i*j

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]:
%%timeit
for a, b in zip(A, B):
a*b

10 loops, best of 3: 63.9 ms per loop


In [36]:
import itertools

In [35]:
%%timeit
for a, b in itertools.izip(A, B):
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

In [26]:
import pickle
obj = ['a', 'list', 'of', 'words']
pickle.dump(obj, open('alistofwords.pkl', 'wb'), -1)
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
shelf.close()
newshelf = shelve.open('test.shelf')
print newshelf['meaning']

42


# 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]:
subplot(121)
imshow(getimg())
subplot(122)
imshow(getimg_vectorised())

Out[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


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
break
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]

3.0


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!

# joblib¶

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

@cache_mem.cache
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
5.29
Computing f(3.5)
--- I am currently computing x*x for x = 3.5
12.25
Computing f(2.3) again
5.29


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 *
set_global_preferences(
useweave=True,
weavecompiler='gcc',
gcc_options=['-ffast-math', '-march=native'],
usecodegen=True,
usecodegenweave=True,
usenewpropagate=True,
usecstdp=True,
)

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):
set_global_preferences(usecodegen=codegen)
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()
net.run(duration)
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)
gc.collect()
defaultclock.t = 0*second


Explanation:

• 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():
clear(True,True)
gc.collect()
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):
brianreset()

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
run(duration)
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

Out[17]:
[<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]:
%%px
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

Out[20]:
[<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):
brianreset()

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
run(duration)
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

Out[25]:
[<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 *

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
if dataman.itemcount()<N:
M = N-dataman.itemcount()
X, Y = zip(*dataman.values())
plot(X, Y, '.')
xlabel('k')
ylabel('Firing rate (Hz)')
show()

In [2]:
from IPython.core.display import Image

Out[2]:

# 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]:
brianreset()

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


# Some useful Brian features you may not know¶

### TimedArray¶

In [6]:
brianreset()

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)

run(duration)

subplot(211)
M['I'].plot()
ylabel('I (amp)')
subplot(212)
M['V'].plot()
ylabel('V (volt)')

Out[6]:
<matplotlib.text.Text at 0x9ac5630>

In [7]:
brianreset()

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)
M_nervefibres = MultiStateMonitor(nervefibres, record=True)

run(duration)

subplot(221)
M_haircells['input'].plot()
ylabel('haircell.input')
subplot(222)
M_haircells['y'].plot()
ylabel('haircell.y')
subplot(223)
M_nervefibres['I'].plot()
ylabel('nervefibres.I')
subplot(224)
M_nervefibres['V'].plot()
ylabel('nervefibres.V')

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

Out[7]:
<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]:
brianreset()

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()

run(1e10*second)

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)'))
client.stop()


# 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():
brianreset()

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


But this will not:

In [11]:
def getsol_numerical():
brianreset()

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)
run(100*ms)
return M[0]

In [12]:
exact = getsol_exact()
numerical = getsol_numerical()
figure(figsize=(10,4))
subplot(121)
title('Solutions')
plot(exact)
plot(numerical)
subplot(122)
title('Difference')
plot(exact-numerical)

Out[12]:
[<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]:
brianreset()
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)
run(500*ms)
figure(figsize=(10, 8))
M.plot()
ylim(-0.1, 3.0)
legend()

brian.experimental.new_c_propagate: WARNING  Using new C based propagation function.

Out[23]:
<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]:
brianreset()
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)
run(500*ms)
figure(figsize=(10, 8))
M.plot()
ylim(-0.1, 3.0)
legend()

brian.experimental.new_c_propagate: WARNING  Using new C based propagation function.

Out[25]:
<matplotlib.legend.Legend at 0x974ecf0>

Any questions?