Ignore this, just a few imports needed for the tutorial.
from brian import *
import time
import pickle
import joblib
import gc
def brianreset():
clear(True,True)
gc.collect()
defaultclock.t = 0*second
Python has some great datatypes, it's worth understanding and using them.
X[key] = value
Don't do this:
%%timeit
for i in range(1000):
for j in range(1000):
i*j
10 loops, best of 3: 42.6 ms per loop
Do this instead:
%%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:
A = rand(100000)
B = rand(100000)
%%timeit
for a, b in zip(A, B):
a*b
10 loops, best of 3: 63.9 ms per loop
Do this instead:
import itertools
%%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.itertools
module of Pythonimport 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
import shelve
shelf = shelve.open('test.shelf')
shelf['meaning'] = 42
shelf.close()
newshelf = shelve.open('test.shelf')
print newshelf['meaning']
42
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.
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.
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)
There are several other useful toolboxes that are worth finding out about.
Don't do this:
a = linspace(0, 1, 100)
b = linspace(0, 1, 100)
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:
def getimg_vectorised():
return (a[:, newaxis]-0.5)**2+(b[newaxis, :]-0.5)**2
They're equivalent:
subplot(121)
imshow(getimg())
subplot(122)
imshow(getimg_vectorised())
<matplotlib.image.AxesImage at 0x4b8f590>
And the vectorised version is much faster!
%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.
delta = 0.005
train1 = sort(rand(50))
train2 = sort(rand(55))
Don't do this:
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:
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.
print cf_loop(), cf_vectorised()
19 19
%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.
You can run C++ code inline in various simple ways:
Scipy has a package weave which allows inline C++ code:
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:
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.
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.
There are several key elements to writing efficient Brian code:
Connection
and Synapses
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
packagecspikequeue
packageBrian 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):
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:
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
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
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).
Some users have run out of memory when working with very large numbers of synapses in Connection
objects. This is because:
To solve this problem, break your Connection
into multiple sub-connections:
[0, ..., N-1]
to [0, ..., N-1]
[0, ..., N/2-1]
to [0, ..., N-1]
[N/2-1, ..., N-1]
to [0, ..., N-1]
Here's how to do a loop over several parameters:
import gc
for params in []:
clear(True, True)
gc.collect()
defaultclock.t = 0*second
# your code here
Explanation:
run()
it collects all existing Brian objectsrun()
you might be doing twice as much work as the first time.clear(True, True)
deletes all memory from currently existing Brian objects.gc.collect()
forces Python to deallocate these objects.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:
def brianreset():
clear(True,True)
gc.collect()
defaultclock.t = 0*second
Another way to run a simulation for several parameters is 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:
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:
%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):
if 0:
import multiprocessing
pool = multiprocessing.Pool()
plot(V0, pool.map(IF, V0))
Using IPython:
from IPython import parallel
rc = parallel.Client()
view = rc[:]
%%px
from brian import *
import gc
%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:
run_tasks
function, but we'll come back to that later.Can you just do something like this?
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
[<matplotlib.lines.Line2D at 0x93c20d0>]
Good options include:
Brian has a specific solution based on shelf files, but which allows for the use of multiple CPUs and machines in parallel.
DataManager
objectrun_tasks
functionmap
DataManager
objectif 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, '.')
xlabel('k')
ylabel('Firing rate (Hz)')
show()
from IPython.core.display import Image
Image(filename='taskcontrol.jpg')
When you call run()
, what happens is:
Network()
objectNetwork.run()
on thisThe way Network.run()
works is:
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.
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, :]);
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)')
<matplotlib.text.Text at 0x9ac5630>
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)
nervefibres.I = linked_var(haircells, 'y')
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
<matplotlib.text.Text at 0xa005050>
Allows us to control a running Brian script.
Run the script with a RemoteControlServer
object, like this:
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):
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()
When you create a NeuronGroup
in Brian, it automatically selects what it thinks is the best way to integrate the equations.
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:
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:
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]
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)
[<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).
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:
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.
<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:
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.
<matplotlib.legend.Legend at 0x974ecf0>
Any questions?