spiketools.population

 1"""Population-level summary statistics for binned spike matrices.
 2
 3Examples
 4--------
 5Shared example setup used throughout the documentation:
 6
 7```python
 8import numpy as np
 9
10from spiketools import gamma_spikes, spiketimes_to_binary
11from spiketools.population import synchrony
12
13np.random.seed(0)
14rates = np.array([6.0] * 10 + [5.6, 6.3, 5.9, 6.5, 5.8, 6.1, 5.7, 6.4, 6.0, 5.5], dtype=float)
15orders = np.array([0.2] * 10 + [1.0, 2.0, 2.0, 3.0, 1.0, 2.0, 3.0, 2.0, 1.0, 3.0], dtype=float)
16spiketimes = gamma_spikes(rates=rates, order=orders, tlim=[0.0, 5000.0], dt=1.0)
17
18binary, _ = spiketimes_to_binary(spiketimes, tlim=[0.0, 5000.0], dt=50.0)
19sync = synchrony(binary)
20```
21"""
22
23from __future__ import annotations
24
25import pylab
26
27__all__ = ["synchrony"]
28
29
30def synchrony(spikes, ignore_zero_rows=True):
31    r"""
32    Calculate the Golomb & Hansel (2000) population synchrony measure.
33
34    Parameters
35    ----------
36    spikes:
37        Binned spike matrix with shape `(n_units, n_time_bins)` or a stack of
38        such matrices with trials along the first axis.
39    ignore_zero_rows:
40        If `True`, units with zero spikes are excluded from the statistic.
41
42    Returns
43    -------
44    float
45        Synchrony estimate between `0` and `1` for typical inputs.
46
47    Definition
48    ----------
49    If $x_i(t)$ is the binned activity of unit $i$ and $\langle \cdot \rangle_i$
50    denotes the population average, this implementation returns
51
52    $$
53    \chi =
54    \sqrt{
55    \frac{\mathrm{Var}_t\left[\langle x_i(t) \rangle_i\right]}
56    {\left\langle \mathrm{Var}_t[x_i(t)] \right\rangle_i}
57    }.
58    $$
59
60    Values near `0` indicate largely independent activity, while larger values
61    indicate that the population fluctuates together on the chosen time grid.
62
63    Notes
64    -----
65    This function expects a dense spike-count matrix, not canonical
66    `spiketimes`. Convert first with `spiketimes_to_binary(...)` if needed.
67
68    Examples
69    --------
70    >>> round(float(synchrony(pylab.array([[1, 0, 1], [0, 1, 0]]))), 3)
71    0.0
72    """
73    if len(spikes.shape) > 2:
74        return pylab.array([synchrony(s, ignore_zero_rows) for s in spikes]).mean()
75    if ignore_zero_rows:
76        mask = spikes.sum(axis=1) > 0
77        sync = spikes[mask].mean(axis=0).var() / spikes[mask].var(axis=1).mean()
78    else:
79        sync = spikes.mean(axis=0).var() / spikes.var(axis=1).mean()
80    return sync ** 0.5