.. currentmodule:: brian

.. index::
	single: efficient code

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 :ref:`compiled-code`.

Brian specifics
---------------

.. removed for Brian 1.0

	The slowest part of Brian is the :class:`QuantityArray` class.
	This is an extension of the NumPy ``ndarray`` array class with
	added support for physical units, unfortunately it's very slow
	in comparison to NumPy. For writing code which runs once, these
	arrays with units are a good idea because they check that you
	are using the right units, thus avoiding bugs, but for code that
	runs repeatedly (in a loop for example), it should be avoided.
	
	If you have a :class:`QuantityArray` object ``x``, you can convert it
	to a NumPy array by writing ``y=asarray(x)``. If you do this only
	once and from then on use ``y`` instead of ``x`` you will see a
	significant increase in speed.
	
	Also note that most Brian code that returns arrays of data supports
	a mechanism for returning a version without units. For example,
	if ``G`` is a :class:`NeuronGroup` with a variable ``V`` you would
	write ``G.V`` to get the array of values of ``V`` for each neuron
	in the group. This array will be a :class:`QuantityArray` and so
	will be slow. For code that is running often or in a loop, you can
	write ``G.V_`` to return the NumPy array version, which will be
	much faster. This same mechanism works for any :class:`NeuronGroup`
	neuron variable, just add an underscore ``_`` character to the end
	of the variable name to refer to the NumPy version without units.
	
	As an example, suppose you were writing a custom reset function,
	you might do something like this::
	
		def myreset(P, spikes):
			P.V[spikes] = 0*mV
			P.Vt[spikes] += 2*mV
	
	This will work, but it will be very slow, so you should write this
	instead::
	
		def myreset(P, spikes):
			P.V_[spikes] = 0*mV
			P.Vt_[spikes] += 2*mV
	
	Note that the constants can still be given with units, if you set
	a NumPy array to a value with units it just ignores the units. It's
	much better to write ``2*mV`` than ``0.002`` even though they're
	equivalent.

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.

.. index::
	pair: efficient code; vectorisation
	single: vectorisation
	single: speed; vectorisation

.. _efficiency-vectorisation:

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 <http://www.scipy.org/PerformancePython>`__
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 :class:`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 :class:`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 <http://www.scipy.org/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.