Source code for gridcells.analysis.fields

'''
===============================================================
:mod:`gridcells.analysis.fields` - grid field related analysis
===============================================================

The :mod:`~gridcells.analysis.fields` module contains routines to analyse
spiking data either from experiments involoving a rodent running in an arena or
simulations involving an animat running in a simulated arena.

Functions
---------
.. autosummary::

    gridnessScore
    occupancy_prob_dist
    spatialAutoCorrelation
    spatialRateMap

'''
from __future__ import absolute_import, division, print_function

import os

import numpy    as np
import numpy.ma as ma

from scipy.integrate             import trapz
from scipy.signal                import correlate2d
from scipy.ndimage.interpolation import rotate

# Do not import when in RDT environment
on_rtd = os.environ.get('READTHEDOCS', None) == 'True'
if not on_rtd:
    from . import _fields
from ..core import Pair2D


[docs]def spatialRateMap(spikeTimes, positions, arena, sigma): '''Compute spatial rate map for spikes of a given neuron. Preprocess neuron spike times into a smoothed spatial rate map, given arena parameters. Both spike times and positional data must be aligned in time! The rate map will be smoothed by a gaussian kernel. Parameters ---------- spikeTimes : np.ndarray Spike times for a given neuron. positions : gridcells.core.Position2D Positional data for these spikes. The timing must be aligned with ``spikeTimes`` arena : gridcells.core.Arena The specification of the arena in which movement was carried out. sigma : float Standard deviation of the Gaussian smoothing kernel. Returns ------- rateMap : np.ma.MaskedArray The 2D spatial firing rate map. The shape will be determined by the arena type. ''' spikeTimes = np.asarray(spikeTimes, dtype=np.double) edges = arena.getDiscretisation() rateMap = _fields.spatialRateMap(spikeTimes, positions.x, positions.y, positions.dt, edges.x, edges.y, sigma) # Mask values which are outside the arena rateMap = np.ma.MaskedArray(rateMap, mask=arena.getMask(), copy=False) return rateMap.T
[docs]def spatialAutoCorrelation(rateMap, arenaDiam, h): '''Compute autocorrelation function of the spatial firing rate map. This function assumes that the arena is a circle and masks all values of the autocorrelation that are outside the `arenaDiam`. .. warning:: This function will undergo serious interface changes in the future. Parameters ---------- rateMap : np.ndarray Spatial firing rate map (2D). The shape should be `(arenadiam/h+1, arenadiam/2+1)`. arenaDiam : float Diameter of the arena. h : float Precision of the spatial firing rate map. Returns ------- corr : np.ndarray The autocorrelation function, of shape `(arenadiam/h*2+1, arenaDiam/h*2+1)` xedges, yedges : np.ndarray Values of the spatial lags for the correlation function. The same shape as `corr.shape[0]`. ''' precision = arenaDiam/h xedges = np.linspace(-arenaDiam, arenaDiam, precision*2 + 1) yedges = np.linspace(-arenaDiam, arenaDiam, precision*2 + 1) X, Y = np.meshgrid(xedges, yedges) corr = ma.masked_array(correlate2d(rateMap, rateMap), mask = np.sqrt(X**2 + Y**2) > arenaDiam) return corr, xedges, yedges
[docs]def gridnessScore(rateMap, arenaDiam, h, corr_cutRmin): '''Calculate gridness score of a spatial firing rate map. Parameters ---------- rateMap : np.ndarray Spatial firing rate map. arenaDiam : float The diameter of the arena. h : float Precision of the spatial firing rate map. Returns ------- G : float Gridness score. crossCorr : np.ndarray An array containing cross correlation values of the rotated autocorrelations, with the original autocorrelation. angles : np.ndarray An array of angles corresponding to the `crossCorr` array. Notes ----- This function computes gridness score accoring to [1]_. The auto correlation of the firing rate map is rotated in 3 degree steps. The resulting gridness score is the difference between a minimum of cross correlations at 60 and 90 degrees, and a maximum of cross correlations at 30, 90 and 150 degrees. The center of the auto correlation map (given by corr_cutRmin) is removed from the map. References ---------- .. [1] Hafting, T. et al., 2005. Microstructure of a spatial map in the entorhinal cortex. Nature, 436(7052), pp.801-806. ''' rateMap_mean = rateMap - np.mean(np.reshape(rateMap, (1, rateMap.size))) autoCorr, autoC_xedges, autoC_yedges = SNAutoCorr(rateMap_mean, arenaDiam, h) # Remove the center point and X, Y = np.meshgrid(autoC_xedges, autoC_yedges) autoCorr[np.sqrt(X**2 + Y**2) < corr_cutRmin] = 0 da = 3 angles = list(range(0, 180+da, da)) crossCorr = [] # Rotate and compute correlation coefficient for angle in angles: autoCorrRot = rotate(autoCorr, angle, reshape=False) C = np.corrcoef(np.reshape(autoCorr, (1, autoCorr.size)), np.reshape(autoCorrRot, (1, autoCorrRot.size))) crossCorr.append(C[0, 1]) max_angles_i = np.array([30, 90, 150]) / da min_angles_i = np.array([60, 120]) / da maxima = np.max(np.array(crossCorr)[max_angles_i]) minima = np.min(np.array(crossCorr)[min_angles_i]) G = minima - maxima return G, np.array(crossCorr), angles
def extractSpikePositions(spikeTimes, positions): spikeIdx = spikeTimes / positions.dt pos_x = _fields.extractSpikePos(spikeIdx, positions.x) pos_y = _fields.extractSpikePos(spikeIdx, positions.y) return Pair2D(pos_x, pos_y), np.max(spikeIdx)
[docs]def occupancy_prob_dist(arena, pos): '''Calculate a probability distribution for animal positions in an arena. Parameters ---------- arena : :class:`~gridcells.core.arena.Arena` Arena the animal was running in. pos : :class:`~gridcells.core.common.Position2D` Positions of the animal. Returns ------- dist : numpy.ndarray Probability distribution for the positional data, given the discretisation of the arena. The first dimension is the y axis, the second dimension is the x axis. The shape of the distribution is equal to the number of items in the discretised edges of the arena. ''' edges = arena.getDiscretisation() dx = arena.getDiscretisationSteps() xedges = np.hstack((edges.x, [edges.x[-1] + dx.x])) yedges = np.hstack((edges.y, [edges.y[-1] + dx.y])) H, _, _ = np.histogram2d(pos.x, pos.y, bins=[xedges, yedges], normed=False) return (H / len(pos)).T