1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119

.. currentmodule:: brian
.. index:: plotting
Analysis and plotting
=====================
.. index::
pair: plotting; pylab
Most plotting should be done with the PyLab commands, all of
which are loaded when you import Brian. See:
http://matplotlib.sourceforge.net/matplotlib.pylab.html
for help on PyLab. The scientific library `Scipy <http://www.scipy.org>`__ is also automatically
imported by the instruction ``from brian import *``.
The most useful plotting instruction is the Pylab function ``plot``. A typical use with Brian is::
plot(t/ms,vm/mV)
where t is a vector of times with units ms and vm is a vector of voltage values with units mV.
To display the figures on the screen, the function ``show()`` must be called once (this should be the
last line of your script), except when using IPython with the Pylab mode (``ipython pylab``).
Brian currently defines just two plotting functions of its own,
:func:`raster_plot` and :func:`hist_plot`. In addition, the :class:`StateMonitor`
object has a :meth:`~StateMonitor.plot` method.
Raster plots

Spike trains recorded by a :class:`SpikeMonitor` can be displayed as raster plots::
S=SpikeMonitor(group)
...
raster_plot(S)
Usual options of the ``plot`` command can also be passed to :func:`raster_plot`. One may also pass
several spike monitors as arguments.
State variable plots

State values recorded by a :class:`StateMonitor` can also be plotted as follows::
M = StateMonitor(group, 'V', record=[0,1,2])
...
M.plot()
Realtime plotting

Both :func:`raster_plot` and :meth:`StateMonitor.plot` have realtime versions
which update as the simulation runs, for example::
G = NeuronGroup(...)
spikemon = SpikeMonitor(G)
statemon = StateMonitor(G, 'V', record=range(5))
ion()
subplot(211)
raster_plot(spikemon, refresh=10*ms, showlast=200*ms)
subplot(212)
statemon.plot(refresh=10*ms, showlast=200*ms)
run(1*second)
ioff()
show()
The ``ion()`` and ``ioff()`` command activate and deactivate Pylab's interactive
plotting mode. The ``refresh`` parameter specifies how often (in simulation
time) to refresh the plot  smaller values will slow down the simulation. The
``showlast`` option only plots the most recent values.
With some IDEs, you may need to do something like the following at the
beginning of your script to make interactive mode work::
import matplotlib
matplotlib.use('WXAgg')
This is because the default graphical backend can sometimes interact badly with
the IDE. Other options to try are ``GTKAgg``, ``QTAgg``, ``TkAgg``.
Statistics

Here are a few functions to analyse first and second order statistical properties
of spike trains, defined as ordered lists of spike times:
* Firing rate: ``firing_rate(spikes)`` where ``spikes`` is a spike train (list of spike times).
* Coefficient of variation: ``CV(spikes)``.
* Crosscorrelogram: ``correlogram(T1,T2,width=20*ms,bin=1*ms,T=None)`` returns
the crosscorrelogram of spike trains T1 and T2 with lag in [width,width] and given bin size.
T is the total duration (optional) and should be greater than the duration of T1 and T2.
The result the rate of coincidences in each bin, returned as an array.
* Autocorrelogram: ``autocorrelogram(T0,width=20*ms,bin=1*ms,T=None)`` is the same as
``correlogram(T0,T0,width=20*ms,bin=1*ms,T=None)``.
* Crosscorrelation function: ``CCF(T1,T2,width=20*ms,bin=1*ms,T=None)`` returns the crosscorrelation
function of T1 and T2, which is the same as the crosscorrelogram divided by the bin size (which makes the
result independent of the bin size).
* Autocorrelation function: ``ACF(T0,width=20*ms,bin=1*ms,T=None)``, same as
``CCF(T0,T0,width=20*ms,bin=1*ms,T=None)``.
* Crosscovariance function: ``CCVF(T1,T2,width=20*ms,bin=1*ms,T=None)`` is the crosscorrelation
function of T1 and T2 minus for the crosscorrelation of independent spike trains with the same rates
(product of rates).
* Autocovariance function: ``ACVF(T0,width=20*ms,bin=1*ms,T=None)`` is the same as
``CCVF(T0,T0,width=20*ms,bin=1*ms,T=None)``.
* Total correlation coefficient: ``total_correlation(T1,T2,width=20*ms,T=None)`` is
the integral of the crosscovariance function divided by the rate of T1, typically (but not
always) between 0 and 1.
* Vector strength: ``vector_strength(spikes,period)`` returns the vector strength of the given
spike train with respect to the period. If each spike time with phase phi is represented by
a vector with angle phi, then the vector strength is the length of the average vector.
It equals 1 for spikes with constant phase and 0 for homogeneous phase distributions.
* Gamma precision factor: ``gamma_factor(source, target, delta)`` returns the gamma precision
factor between source and target trains, with precision delta.
These functions return NaN (not a number) when a spike train is empty.
.. TODO: docs:histograms
Histograms

