Main code structure =================== Overview -------- Brian features can be broadly categorised into *construction* of the network, and *running* the network. Constructing the network ~~~~~~~~~~~~~~~~~~~~~~~~ The following objects need to be specified by the user explicitly or implicitly: * :class:`NeuronGroup` * :class:`Connection` * Monitors * :class:`Network` After that, the network needs to be *prepared*. Preparation of the network involves initialising objects' data structures appropriately, in particular compressing the :class:`Connection` matrices. Connection matrices are initially stored as instances of a :class:`ConstructionMatrix` class (sparse, dense, etc.), and then later *compressed* into an instance of a :class:`ConnectionMatrix` class. Two levels are necessary, because at construction time, all matrices have to be editable, whereas at runtime, for efficiency reasons, some matrix types are read-only or partially read-only. Data structures appropriate to the construction of a matrix, particularly sparse matrices, are not the most efficient for runtime access. Constructing the :class:`NeuronGroup` object is a rather complicated operation, involving the construction of many subsididary objects. The most complicated aspect is the creation, manipulation and analysis of an :class:`Equations` object. Running the network ~~~~~~~~~~~~~~~~~~~ The network is run by repeatedly evaluating the 'update schedule' and updating the clock or clocks. The 'update schedule' is user specifiable, but usually consists of the following sequence of operations (interspersed with optional user network operation calls): * Update state variables of :class:`NeuronGroup` * Call thresholding function * Push spikes into :class:`SpikeContainer` * Propagate spikes (possibly with delays) via :class:`Connection` * Update state variables of :class:`Synapses` (possibly includes updating the state of targeted :class:`NeuronGroup` objects * Call reset function on neurons which have spiked Details of network construction ------------------------------- Construction of :class:`NeuronGroup` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The :class:`NeuronGroup` object is responsible for storing the state variables of each of its neurons, for updating them each time step, generating spikes via a thresholding mechanism, storing spikes so that they can be accessed with a delay, and resetting state variables after spiking. State variable update is done by a ``StateUpdater`` class defined in ``brian/stateupdater.py``. Thresholding is done by a :class:`Threshold` class defined in ``brian/threshold.py`` and resetting is done by a :class:`Reset` class defined in ``brian/reset.py``. The ``__init__`` method of :class:`NeuronGroup` takes these objects as arguments, but it also has various additional keywords which can be used more intuitively. In this case, the appropriate object is selected automatically. For example, if you specify ``reset=0*mV`` in the keyword arguments, Brian generates a ``Reset(0*mV)`` object. The :meth:`NeuronGroup.__init__` method code is rather complicated and deals with many special cases. The most complicated aspect of this is the definition of the state variables and the update procedure. Typically, the user simply gives a list of differential equations, and Brian attempts to automatically extract the appropriate state variable definitions, and creates a differential equation solver appropriate to them (it needs some help in this at the moment, e.g. specifying the order or type of the solver). The main work in this is done by the :func:`magic_state_updater` function, which uses the :class:`Equations` object (see next section). Once the state variables are defined, various internal objects are created. The state variables are stored in the ``_S`` attribute of a :class:`NeuronGroup`. This is an MxN matrix where M is the number of variables and N is the number of neurons. The other major data structure generated is the ``LS`` attribute (last spikes). This is a ``SpikeContainer`` instance, a circular array used to contain spikes. See ``brian/utils/circular.py``. Finally note that the construction of many of these objects requires a :class:`Clock` object, which can either be specified explicitly, or is guessed by the :func:`guess_clock` function which searches for clocks using the magic module (see below). :class:`EventClock` objects are excluded from this guessing. The :func:`magic_state_updater` function and the :class:`Equations` object ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ This function returns an appropriate ``StateUpdater`` object and a list of the dynamic variables of an :class:`Equations` object. It uses methods of the :class:`Equations` object to do this (such as :meth:`Equations.is_linear`). The :class:`Equations` object can be constructed in various ways. It can be constructed from a multi-line string or by adding (concatenating) other :class:`Equations` objects. Construction by multi-line string is done by pattern matching. The four forms are: 1. ``dx/dt = f : unit`` (differential equation) 2. ``x = f : unit`` (equation) 3. ``x = y`` (alias) 4. ``x : unit`` (parameter) Differential equations and parameters are dynamic variables, equations and aliases are just computed and substituted when necessary. The ``f`` patterns in the forms above are stored for differential equations and parameters. For the solution of nonlinear equations, these ``f`` patterns are executed as Python code in the evaluation of the state update. For example, the equations ``dV/dt = -V*V/(10*ms) : 1`` and ``dW/dt = cos(W)/(20*ms) : 1`` are numerically evaluated with an Euler method as the following code (generated from the list of dynamic variables and their defining strings):: V, W = P._S V__tmp, W__tmp = P._dS V__tmp[:] = -V*V/(10*ms) W__tmp[:] = cos(W)/(20*ms) P._S += dt*P._dS This code generation is done by the :meth:`Equations.forward_euler` family of methods. In the case where the equations are linear, they are solved by a matrix exponential using the code in :func:`get_linear_equations` defined in ``brian/stateupdater.py``. Finally, note that equations strings can contain references to names of objects that are defined in the namespace of the string, and the :class:`Equations` object can pick these out. It does this by inspecting the call stack, extracting the namespace for the appropriate level (which has to be kept track of), and plucking out the appropriate name. The ``level=...`` keywords you see dotted around Brian's code are there to keep track of these levels for this reason. Construction of :class:`Connection` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ :class:`Connection` objects provide methods for storing weight matrices and propagating spikes. Spike propagation is done via the :meth:`Connection.do_propagate` and :meth:`Connection.propagate` methods. Weight matrices are stored in the ``W`` attribute. Initially, weight matrices are :class:`ConstructionMatrix` objects and are converted by the :meth:`Connection.compress` method, via the matrices' ``connection_matrix()`` methods to :class:`ConnectionMatrix` objects. The idea is to have two data structures, one appropriate to the construction of a matrix, supporting adding and removing new synapses, and one appropriate to runtime behaviour, focussing on fast row access above all else. There are three matrix structures, 'dense', 'sparse' and 'dynamic' (and 'computed' may be added later). The 'dense' matrix is just a full 2D array, and the matrix objects just reproduce the functionality of numpy arrays. The 'sparse' and 'dynamic' structures are sparse matrices. The first doesn't allow you to add or remove elements at runtime and is very optimised, the second does allow you to and is less optimised. Both are more optimised than scipy's sparse matrices, which are used as the basis for the construction phase. :class:`Connection` objects can handle homogeneous delays (all synapses with the same delay) by pulling spikes from :class:`NeuronGroup`'s ``LS`` object with a delay. Heterogeneous delays (each synapse with a different delay) are done by the :class:`DelayConnection` object, which stores a ``_delayvec`` attribute alongside the ``W`` attribute. The delay matrix is of the same type as the weight matrix and in the case of sparse matrices must have nonzero elements at the same places. The :class:`Connection` object automatically turns itself into a :class:`DelayConnection` object if heterogeneous delays are requested, so that the user doesn't even need to know of the existence of :class:`DelayConnection`. The :class:`Connection` object also provides methods for initialising the weight matrices either fully, randomly, etc. (See the user docs.) Construction of monitors ~~~~~~~~~~~~~~~~~~~~~~~~ The :class:`SpikeMonitor` class of monitors derive from :class:`Connection`. Rather than propagating spikes to another group, they store them in a list. The work is done in the :meth:`SpikeMonitor.propagate` method. The :class:`StateMonitor` class of monitors derive from :class:`NetworkOperation`. Network operations are called once every time step and execute arbitrary Python code. The :class:`StateMonitor` class of monitors record state variables each time step to a list. Construction of :class:`Network` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ When the user calls the function :func:`run`, a :class:`MagicNetwork` object is created, and the :meth:`Network.run` method is called. :class:`MagicNetwork` gathers, using the magic module, a list of all appropriate objects and runs them together. Alternatively, the user can specify their own list of objects using a :class:`Network` object. Each time an object is added to a :class:`Network` either via the initialiser or the :meth:`Network.add` method, it checks to see if it has an attribute ``contained_objects``, and if so it adds all the objects in that to the network too. This allows, for example, the :class:`STDP` object to contained :class:`NeuronGroup` and :class:`Connection` objects which the user doesn't see, but are used to implement the STDP functionality. The :meth:`Network.run` method calls the :meth:`Connection.compress` method on every :class:`Connection` object to convert construction matrices to connection matrices. It also builds an update schedule (see below). The magic module ~~~~~~~~~~~~~~~~ The magic module is used for tracking instances of objects. A class that derives from the :class:`magic.InstanceTracker` class can be tracked in this way, including :class:`NeuronGroup`, :class:`Connection`, :class:`NetworkOperation` and :class:`Clock`. The :func:`find_instances` function can be used to search for instances. Note that :func:`find_instances` will only return instances which were instantiated in the same execution frame as the :func:`find_instances` calling frame, or (if the ``level`` keyword is used) one of the frames higher up in the call stack. The idea is that you may have modular code with objects defined in different places, but that you don't want to use all objects that exist at all in the network. This system causes a bit of trouble, but seems unavoidable. See the user manual section :ref:`projects-with-multiple-files` for details on getting around this. Details of network running -------------------------- Update schedules ~~~~~~~~~~~~~~~~ An update schedule gives the sequence of operations to be carried out each time step of the simulation. Typically this is: state update and threshold; propagation; reset, although an option is available for switching propagation and reset around. Network operations can be weaved in anywhere amongst these basic steps. See the reference documentation for the :class:`Network` object for more details. Simulation proceeds by choosing the clock with the lowest current time, selecting all objects which have that clock as their clock, and performing the update schedule on those objects, before applying the :meth:`Clock.tick` method to increment the clock time by ``dt``. Network operations ~~~~~~~~~~~~~~~~~~ A :class:`NetworkOperation` object is called as if it were a function, i.e. the ``__call__()`` method is called with no arguments. The :func:`network_operation` decorator exists to convert a function into a suitable :class:`NetworkOperation` object. This technique is used for the internal functioning of many of Brian's features (such as :class:`StateMonitor`). :class:`NeuronGroup` update ~~~~~~~~~~~~~~~~~~~~~~~~~~~ The :meth:`NeuronGroup.update` method does the following three things. First of all it calls the ``StateUpdater`` to update the state variables. Then it calls its :class:`Threshold` object if one exists to extract the indices of the spiking neurons. Finally it pushes these into the ``LS`` attribute for extraction by any :class:`Connection` objects. :class:`NeuronGroup` reset ~~~~~~~~~~~~~~~~~~~~~~~~~~ The :meth:`Reset.__call__` method pulls spike from the :class:`NeuronGroup`'s ``LS`` attribute and then resets the state variables for those. Spike propagation ~~~~~~~~~~~~~~~~~ The :meth:`Connection.do_propagate` method does two things, it gets the spike indices to propagate (with homogeneous delays if chosen) from the ``LS`` attribute of the :class:`NeuronGroup` and then passes these to its :meth:`Connection.propagate` method. This method extracts a list of connection matrix rows using the :meth:`ConnectionMatrix.get_rows` method. This method returns a list of :class:`ConnectionVector` instances. There are two types of :class:`ConnectionVector`, dense and sparse. Dense ones are simply numpy arrays, sparse ones consist of two numpy arrays, an array of values and an array of corresponding column indices. The :class:`SparseConnectionVector` class has some methods which make this distinction seamless to the user in most instances, although developers need to be aware of it. Finally, the :meth:`Connection.propagate` method goes through this list applying the row vectors one by one. The pure Python version of this is straightforward, but there is also a ``C++`` accelerated version which uses the scipy Weave package if the user has a suitable compiler on their system. This version is much more efficient, but the code is rather dense and difficult to understand.