spiketools.surrogates

  1"""Spike-train surrogate generation utilities.
  2
  3Examples
  4--------
  5Generate the shared documentation dataset:
  6
  7```python
  8import numpy as np
  9
 10from spiketools import gamma_spikes
 11
 12np.random.seed(0)
 13rates = 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)
 14orders = 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)
 15spiketimes = gamma_spikes(rates=rates, order=orders, tlim=[0.0, 5000.0], dt=1.0)
 16```
 17"""
 18
 19from __future__ import annotations
 20
 21import numpy as np
 22import pylab
 23
 24__all__ = ["gamma_spikes"]
 25
 26
 27def gamma_spikes(rates, order=[1], tlim=[0.0, 1000.0], dt=0.1):
 28    r"""Generate surrogate spike trains from homogeneous gamma-renewal processes.
 29
 30    Parameters
 31    ----------
 32    rates:
 33        Scalar or per-train firing rates in spikes/s.
 34    order:
 35        Gamma-process shape parameter per train. `1` corresponds to a Poisson
 36        process, values larger than `1` are more regular, and values below `1`
 37        are more bursty.
 38    tlim:
 39        Two-element time interval `[tmin, tmax]` in ms.
 40    dt:
 41        Output quantization step in ms. Generated spike times are rounded to
 42        this grid for backwards compatibility.
 43
 44    Returns
 45    -------
 46    np.ndarray
 47        Canonical `spiketimes` array.
 48
 49    Definition
 50    ----------
 51    For each spike train, consecutive inter-spike intervals are drawn
 52    independently as
 53
 54    $$
 55    I_n \sim \mathrm{Gamma}(k, \theta),
 56    \qquad
 57    k = \mathrm{order},
 58    \qquad
 59    \theta = \frac{1000}{\mathrm{rate}\,k}
 60    $$
 61
 62    so that
 63
 64    $$
 65    \mathrm{E}[I_n] = \frac{1000}{\mathrm{rate}},
 66    \qquad
 67    \mathrm{CV}^2 = \frac{1}{k}.
 68    $$
 69
 70    Examples
 71    --------
 72    >>> np.random.seed(0)
 73    >>> spikes = gamma_spikes([10.0], order=[0.5], tlim=[0.0, 10.0], dt=1.0)
 74    >>> spikes.shape[0]
 75    2
 76    """
 77    rates_arr = np.atleast_1d(np.asarray(rates, dtype=float))
 78    order_arr = np.atleast_1d(np.asarray(order, dtype=float))
 79
 80    if rates_arr.size == 1 and order_arr.size > 1:
 81        rates_arr = np.full(order_arr.shape, rates_arr.item(), dtype=float)
 82    elif order_arr.size == 1 and rates_arr.size > 1:
 83        order_arr = np.full(rates_arr.shape, order_arr.item(), dtype=float)
 84    elif rates_arr.size != order_arr.size:
 85        raise ValueError("`rates` and `order` must have matching lengths or be scalar.")
 86
 87    if np.any(rates_arr <= 0.0):
 88        raise ValueError("All firing rates must be positive.")
 89    if np.any(order_arr <= 0.0):
 90        raise ValueError("All gamma orders must be positive.")
 91
 92    tmin, tmax = float(tlim[0]), float(tlim[1])
 93    if tmax <= tmin:
 94        raise ValueError("`tlim` must satisfy tmax > tmin.")
 95    if dt <= 0.0:
 96        raise ValueError("`dt` must be positive.")
 97
 98    spiketimes = [[], []]
 99    for trial_id, (rate_hz, shape) in enumerate(zip(rates_arr, order_arr)):
100        scale_ms = 1000.0 / (rate_hz * shape)
101        t = tmin + np.random.gamma(shape=shape, scale=scale_ms)
102        generated = []
103        while t < tmax:
104            rounded = np.round((t - tmin) / dt) * dt + tmin
105            if rounded < tmax:
106                generated.append(float(rounded))
107            t += np.random.gamma(shape=shape, scale=scale_ms)
108
109        if generated:
110            spiketimes[0].extend(generated)
111            spiketimes[1].extend([float(trial_id)] * len(generated))
112        else:
113            spiketimes[0].append(pylab.nan)
114            spiketimes[1].append(float(trial_id))
115
116    return pylab.array(spiketimes, dtype=float)