How to write efficient Brian code¶
There are a few keys to writing fast and efficient Brian code. The first is to use Brian itself efficiently. The second is to write good vectorised code, which is using Python and NumPy efficiently. For more performance tips, see also Compiled code.
Brian specifics¶
You can switch off Brian’s entire unit checking module by including the line:
import brian_no_units
before importing Brian itself. Good practice is to leave unit checking on most of the time when developing and debugging a model, but switching it off for long runs once the basic model is stable.
Another way to speed up code is to store references to arrays rather
than extracting them from Brian objects each time you need them. For
example, if you know the custom reset object in the code above is
only ever applied to a group custom_group
say, then you could
do something like this:
def myreset(P, spikes):
custom_group_V_[spikes] = 0*mV
custom_group_Vt_[spikes] = 2*mV
custom_group = ...
custom_group_V_ = custom_group.V_
custom_group_Vt_ = custom_group.Vt_
In this case, the speed increase will be quite small, and probably
not worth doing because it makes it less readable, but in more
complicated examples where you repeatedly refer to custom_group.V_
it could add up.
Vectorisation¶
Python is a fast language, but each line of Python code has an associated overhead attached to it. Sometimes you can get considerable increases in speed by writing a vectorised version of it. A good guide to this in general is the Performance Python page. Here we will do a single worked example in Brian.
Suppose you wanted to multiplicatively depress the connection strengths every time step by some amount, you might do something like this:
C = Connection(G1, G2, 'V', structure='dense')
...
@network_operation(when='end')
def depress_C():
for i in range(len(G1)):
for j in range(len(G2)):
C[i,j] = C[i,j]*depression_factor
This will work, but it will be very, very slow.
The first thing to note is that the Python expression range(N)
actually constructs a list [0,1,2,...,N-1]
each time it is called,
which is not really necessary if you are only iterating over the list.
Instead, use the xrange
iterator which doesn’t construct the list
explicitly:
for i in xrange(len(G1)):
for j in xrange(len(G2)):
C[i,j] = C[i,j]*depression_factor
The next thing to note is that when you call C[i,j] you are doing an
operation on the Connection
object, not directly on the underlying
matrix. Instead, do something like this:
C = Connection(G1, G2, 'V', structure='dense')
C_matrix = asarray(C.W)
...
@network_operation(when='end')
def depress_C():
for i in xrange(len(G1)):
for j in xrange(len(G2)):
C_matrix[i,j] *= depression_factor
What’s going on here? First of all, C.W
refers to the ConnectionMatrix
object, which is a 2D NumPy array with some extra stuff - we don’t need the extra
stuff so we convert it to a straight NumPy array asarray(C.W)
. We also store
a copy of this as the variable C_matrix
so we don’t need to do this every
time step. The other thing we do is to use the *=
operator instead of the *
operator.
The most important step of all though is to vectorise the entire operation. You
don’t need to loop over i
and j
at all, you can manipulate the array
object with a single NumPy expression:
C = Connection(G1, G2, 'V', structure='dense')
C_matrix = asarray(C.W)
...
@network_operation(when='end')
def depress_C():
C_matrix *= depression_factor
This final version will probably be hundreds if not thousands of times faster than the original. It’s usually possible to work out a way using NumPy expressions only to do whatever you want in a vectorised way, but in some very rare instances it might be necessary to have a loop. In this case, if this loop is slowing your code down, you might want to try writing that loop in inline C++ using the SciPy Weave package. See the documentation at that link for more details, but as an example we could rewrite the code above using inline C++ as follows:
from scipy import weave
...
C = Connection(G1, G2, 'V', structure='dense')
C_matrix = asarray(C.W)
...
@network_operation(when='end')
def depress_C():
n = len(G1)
m = len(G2)
code = '''
for(int i=0;i<n;i++)
for(int j=0;j<m;j++)
C_matrix(i,j) *= depression_factor
'''
weave.inline(code,
['C_matrix', 'n', 'm', 'depression_factor'],
type_converters=weave.converters.blitz,
compiler='gcc',
extra_compile_args=['-O3'])
The first time you run this it will be slower because it compiles the
C++ code and stores a copy, but the second time will be much faster as
it just loads the saved copy. The way it works is that Weave converts
the listed Python and NumPy variables (C_matrix
, n
, m
and depression_factor
) into C++ compatible data types. n
and
m
are turned into int``s, ``depression_factor
is turned into
a double
, and C_matrix
is turned into a Weave
Array
class. The only thing you need to know about this is that
elements of a Weave array are referenced with parentheses rather than
brackets, i.e. C_matrix(i,j)
rather than C_matrix[i,j]
. In
this example, I have used the gcc
compiler and added the optimisation
flag -O3
for maximum optimisations. Again, in this case it’s much
simpler to just use the C_matrix *= depression_factor
NumPy expression,
but in some cases using inline C++ might be necessary, and as you can see
above, it’s very easy to do this with Weave, and the C++ code for a
snippet like this is often almost as simple as the Python code would be.