Recording with Monitors

Between two calls to simulate(), all neural and synaptic variables can be accessed through the generated attributes. The evolution of neural or synaptic variables during a simulation phase can be selectively recorded using Monitor objects.

The Monitor object can be created at any time (before or after compile()) to record any variable of a Population, PopulationView, Dendrite or Projection.

Important

The value of each variable is stored for every simulation step in the RAM. For huge networks and long simulations, this can very rapidly fill up the available memory and lead to cache defaults, thereby degrading strongly the performance. It is the user’s responsability to record only the needed variables and to regularly save the values in a file.

Neural variables

The Monitor object takes four arguments:

  • obj: the object to monitor. It can be a population, a population view (a slice of a population or an individual neuron), a dendrite (the synapses of a projection which reach a single post-synaptic neuron) or a projection.
  • variables: a (list of) variable name(s) which should be recorded. They should be variables of the neuron/synapse model of the corresponding object. Although it generally makes no sense, you can also record parameters of an object. By definition a parameter is constant throughout a simulation, but it maybe useful when tracking externally-set inputs, for example. You can know which attributes are recordable by checking the attributes attribute of the object (pop.attributes or proj.attributes).
  • period: the period in ms at which recordings should be made. By default, recording is done after each simulation step (dt), but this may be overkill in long simulations.
  • start: boolean value stating if the recordings should start immediately after the creation of the monitor (default), or if it should be started later.

Some examples:

m = ann.Monitor(pop, 'r') # record r in all neurons of pop
m = ann.Monitor(pop, ['r', 'v']) # record r and v of all neurons
m = ann.Monitor(pop[:100], 'r', period=10.0) # record r in the first 100 neurons of pop, every 10 ms
m = ann.Monitor(pop, 'r', start=False) # record r in all neurons, but do not start recording

Spiking networks additionally allow to record the spike events in a population (see later). You also can record conductances (e.g. g_exc) and weighted sums of inputs in rate-coded networks (sum(exc)) the same way:

m = ann.Monitor(pop, ['spike', 'g_exc', 'g_inh'])
m = ann.Monitor(pop, ['r', 'sum(exc)', 'sum(inh)'])

Starting the recordings

If start is set to False, recordings can be started later by calling the start() method:

m = ann.Monitor(pop, 'r', start=False)
ann.simulate(100.)
m.start()
ann.simulate(100.)

In this case, only the last 100 ms of the simulation will be recorded. Otherwise, recording would start immediately after the creation of the object.

Pausing/resuming the recordings

If you are interested in recording only specific periods of the simulation, you can ause and resume recordings:

m = ann.Monitor(pop, 'r')
ann.simulate(100.)
m.pause()
ann.simulate(1000.)
m.resume()
ann.simulate(100.)

In this example, only the first and last 100 ms of the simulation are recorded.

Retrieving the recordings

The recorded values are obtained through the get() method. If no argument is passed, a dictionary is returned with one element per recorded variable. If the name of a variable is passed (for example get('r')), the recorded values for this variable are directly returned:

m = ann.Monitor(pop, ['r', 'v'])
ann.simulate(100.)
data = m.get()
ann.simulate(100.)
r = m.get('r')
v = m.get('v')

In the example above, data is a dictionary with two keys 'r' and 'v', while r and v are directly the recorded arrays.

The recorded values are Numpy arrays with two dimensions, the first one representing time, the second one representing the ranks of the recorded neurons.

For example, the time course of the firing rate of the neuron of rank 15 is accessed through:

data['r'][:, 15]

The firing rates of the whole population after 50 ms of simulation are accessed with:

data['r'][50, :]
Note

Once you call get(), the internal data is erased, so calling it immediately afterwards will return an empty recording data. You need to simulate again in order to retrieve new values.

Representation of time

The time indices are in simulation steps (integers), not in real time (ms). If dt is different from 1.0, this indices must be multiplied by dt() in order to plot real times:

ann.setup(dt=0.1)
# ...
m = ann.Monitor(pop, 'r')
ann.simulate(100.)
r = m.get('r')
plt.plot(dt()*np.arange(100), r[:, 15])

If recordings used the pause() and resume() methods, get() returns only one array with all recordings concatenated. You can access the steps at which the recording started or paused with the times() method:

m = ann.Monitor(pop, 'r')
ann.simulate(100.)
m.pause()
ann.simulate(1000.)
m.resume()
ann.simulate(100.)
r = m.get('r') # A (200, N) Numpy array
print(m.times()) # {'start': [0, 1100], 'stop': [100, 1200]}

Special case for spiking neurons

Any variable defined in the neuron type can be recorded. An exception for spiking neurons is the spike variable itself, which is never explicitely defined in the neuron type but can be recorded:

m = ann.Monitor(pop, ['v', 'spike'])

Unlike other variables, the binary value of spike is not recorded at each time step, which would lead to very sparse matrices, but only the times (in steps, not milliseconds) at which spikes actually occur.

As each neuron fires differently (so each neuron will have recorded spikes of different lengths), get() in this case does not return a Numpy array, but a dictionary associating to each recorded neuron a list of spike times:

m = ann.Monitor(pop, ['v', 'spike'])
ann.simulate(100.0)
data = m.get('spike')
print(data[0]) # [23, 76, 98]

In the example above, the neuron of rank 0 has spiked 3 times (at t = 23, 76 and 98 ms if dt = 1.0) during the first 100 ms of the simulation.

Raster plots

In order to easily display raster plots, the method raster_plot() is provided to transform this data into an easily plottable format:

spike_times, ranks = m.raster_plot(data)
plt.plot(spike_times, ranks, '.')

raster_plot() returns two Numpy arrays, whose length is the total number of spikes emitted during the simulation. The first array contains the spike times (ín ms) while the second contains the ranks of the neurons who fired. They can be directly used t produce the raster plot with Matplotlib.

An example of the use of raster_plot() can be seen in the Izhikevich pulse network section.

Inter-spike interval (ISI) and coefficient of variation (ISI CV)

In addition to a raster plot, the distribution of inter-spike intervals could be considered for evaluation. The inter-spike interval (short ISI) is defined as the time in milliseconds between two consecutive spike events. The method inter_spike_interval() transforms the recorded spike events into a list of ISI across all neurons (default) or for indivdual neurons (add per_neuron = True to argument list) which could be fed into a histogram method provided by matplotlib:

pop_isi = m.inter_spike_interval(data)
plt.hist(pop_isi)

The coefficient of variation is a measure often reported together with the inter-spike interval. These values can be easily obtained using another function of the Monitor object:

pop_isi_cv = m.coefficient_of_variation(data)
plt.hist(pop_isi_cv)

An example of the use of inter_spike_interval() and coefficient_of_variation() can be seen in the COBA network section.

Mean firing rate

The mean firing rate in the population can be easily calculated using the length of the arrays returned by raster_plot:

N = 1000 # number of neurons
duration = 500. # duration of the simulation
data = m.get('spike')
spike_times, ranks = m.raster_plot(data)
print('Mean firing rate:', len(spike_times)/float(N)/duration*1000., 'Hz.')

For convenience, this value is returned by the mean_fr() method, which has access to the number of recorded neurons and the duration of the recordings:

print('Mean firing rate:', m.mean_fr(data), 'Hz.')

Firing rates

Another useful method is smoothed_rate(). It allows to display the instantaneous firing rate of each neuron based on the spike recordings:

rates = m.smoothed_rate(data)
plt.imshow(rates, aspect='auto')

For each neuron, it returns an array with the instantaneous firing rate during the whole simulation. The instantaneous firing rate is computed by inverting the inter-spike interval (ISI) between two consecutive spikes, and assigning it to all simulation steps between the two spikes.

As this value can be quite fluctuating, a smooth argument in milliseconds can be passed to smoothed_rate() to apply a low-pass filter on the firing rates:

rates = m.smoothed_rate(data, smooth=200.0)
plt.imshow(rates, aspect='auto')

A smoothed firing rate for the whole population is also accessible through population_rate():

fr = m.population_rate(data, smooth=200.0)

Histogram

histogram() allows to count the spikes emitted in the whole population during successive bins of the recording duration:

histo = m.histogram(data, bins=1.0)
plt.plot(histo)

bins represents the size of each bin, here 1 ms. By default, the bin size is dt.

Note : the methods to analyse the spike patterns are also available outside the monitors. For example if you save the spike recordings into a file using numpy:

spikes = m.get('spike')

np.save('spikes.npy', spikes)

you can analyze them in a separate file like this:

# Load the data
spikes = np.load('spikes.npy').item()

# Compute the raster plot
t, n = ann.raster_plot(spikes)

# Compute the population firing rate
fr = ann.histogram(spikes, bins=1.)

# Smoothed firing rate
sr = ann.smoothed_rate(spikes, smooth=10.0)

# Population firing rate
pr = ann.population_rate(spikes, smooth=10.0)

# Global firing rate
mfr = ann.mean_fr(spikes)

Synaptic variables

Recording of synaptic variables such as weights w during learning is also possible using the monitor object. However, it can very easily lead to important memory consumption. Let’s suppose we have a network composed of two populations of 1000 neurons each, fully connected: each neuron of the second population receives 1000 synapses. This makes a total of 1 million synapses for the projection and, supposing the weights w use the double floating precision, requires 4 MB of memory. If you record w during a simulation of 1 second (1000 steps, with dt=1.0), the total added memory consumption would already be around 4GB.

To avoid fast memory fills, you should either 1) record the projection variables infrequently (by setting the period argument of the Monitor), or 2) selectively record particular dendrites. The corresponding dendrite should be simply passed to the monitor:

dendrite = proj.dendrite(12) # or simply proj[12]
m = ann.Monitor(dendrite, 'w')
ann.simulate(1000.0)
data = m.get('w')

The Monitor object has the same start(), pause(), resume() and get() methods as for populations. get() returns also 2D Numpy arrays, the first index being time, the second being the index of the synapse. To know to which pre-synaptic neuron it corresponds, use the pre_ranks attribute of the dendrite:

dendrite.pre_ranks # [0, 3, 12, ..]

To record a complete projection, simply pass it to the Monitor:

m = ann.Monitor(proj, 'w', period=1000.)
ann.simulate(10000.0)
data = m.get('w')

One last time, do not record all weights of a projection at each time step!

Warning

Recording synaptic variables with CUDA is not available.