spiketools

137def binary_to_spiketimes(binary: np.ndarray, time: Iterable[float]) -> np.ndarray:
138    """
139    Convert a binned spike matrix into canonical `spiketimes`.
140
141    Parameters
142    ----------
143    binary:
144        Array with shape `(n_trials, n_bins)`. Values larger than `1` are
145        interpreted as multiple spikes in the same bin.
146    time:
147        One time value per bin, typically the left bin edges in ms.
148
149    Returns
150    -------
151    np.ndarray
152        `float` array with shape `(2, n_spikes_or_markers)`.
153
154    Notes
155    -----
156    Trials with no spikes are retained via placeholder columns
157    `[nan, trial_id]`.
158
159    This function assumes that each non-zero entry denotes one or more spike
160    events in that time bin. It is not intended for persistent state matrices
161    where a neuron stays at `1` across consecutive bins until the next update.
162
163    Examples
164    --------
165    Convert a spike-count raster on a regular grid:
166
167    >>> binary_to_spiketimes(
168    ...     np.array([[1, 0, 2], [0, 0, 0]]),
169    ...     [0.0, 1.0, 2.0],
170    ... ).tolist()
171    [[0.0, 2.0, 2.0, nan], [0.0, 0.0, 0.0, 1.0]]
172    """
173    time = pylab.array(time)
174    spiketimes: list[list[float]] = [[], []]
175    max_count = binary.max()
176    spikes = binary.copy()
177    while max_count > 0:
178        if max_count == 1:
179            trial, t_index = spikes.nonzero()
180        else:
181            trial, t_index = pylab.where(spikes == max_count)
182        spiketimes[0] += time[t_index].tolist()
183        spiketimes[1] += trial.tolist()
184        spikes[spikes == max_count] -= 1
185        max_count -= 1
186    spiketimes_array = pylab.array(spiketimes)
187
188    all_trials = set(range(binary.shape[0]))
189    found_trials = set(spiketimes_array[1, :])
190    missing_trials = list(all_trials.difference(found_trials))
191    for mt in missing_trials:
192        spiketimes_array = pylab.append(
193            spiketimes_array, pylab.array([[pylab.nan], [mt]]), axis=1
194        )
195    return spiketimes_array.astype(float)

239def cut_spiketimes(spiketimes: np.ndarray, tlim: Sequence[float]) -> np.ndarray:
240    """
241    Restrict `spiketimes` to a time interval.
242
243    Parameters
244    ----------
245    spiketimes:
246        Canonical spike representation.
247    tlim:
248        Two-element sequence `[tmin, tmax]` in ms. The interval is inclusive at
249        the start and exclusive at the end.
250
251    Returns
252    -------
253    np.ndarray
254        Cropped spike array. Empty trials are still represented by
255        `[nan, trial_id]` markers.
256
257    Examples
258    --------
259    >>> cut_spiketimes(
260    ...     np.array([[1.0, 3.0, 5.0], [0.0, 0.0, 1.0]]),
261    ...     [2.0, 5.0],
262    ... ).tolist()
263    [[3.0, nan], [0.0, 1.0]]
264    """
265    alltrials = list(set(spiketimes[1, :]))
266    cut_spikes = spiketimes[:, pylab.isfinite(spiketimes[0])]
267    cut_spikes = cut_spikes[:, cut_spikes[0, :] >= tlim[0]]
268
269    if cut_spikes.shape[1] > 0:
270        cut_spikes = cut_spikes[:, cut_spikes[0, :] < tlim[1]]
271    for trial in alltrials:
272        if trial not in cut_spikes[1, :]:
273            cut_spikes = pylab.append(
274                cut_spikes, pylab.array([[pylab.nan], [trial]]), axis=1
275            )
276    return cut_spikes

45def get_time_limits(spiketimes: np.ndarray) -> list[float]:
46    """
47    Infer time limits from a `spiketimes` array.
48
49    Parameters
50    ----------
51    spiketimes:
52        Canonical spike representation with shape `(2, n_spikes)`.
53
54    Returns
55    -------
56    list[float]
57        `[tmin, tmax]` with inclusive start and exclusive end. If inference
58        fails, `[0.0, 1.0]` is returned as a safe fallback.
59
60    Examples
61    --------
62    >>> spikes = np.array([[5.0, 8.0], [0.0, 1.0]])
63    >>> get_time_limits(spikes)
64    [5.0, 9.0]
65    """
66    try:
67        tlim = [
68            float(pylab.nanmin(spiketimes[0, :])),
69            float(pylab.nanmax(spiketimes[0, :]) + 1),
70        ]
71    except Exception:
72        tlim = [0.0, 1.0]
73    if pylab.isnan(tlim).any():
74        tlim = [0.0, 1.0]
75    return tlim

 78def spiketimes_to_binary(
 79    spiketimes: np.ndarray,
 80    tlim: Optional[Sequence[float]] = None,
 81    dt: float = 1.0,
 82) -> tuple[np.ndarray, np.ndarray]:
 83    """
 84    Convert `spiketimes` into a spike-count matrix on a regular time grid.
 85
 86    Parameters
 87    ----------
 88    spiketimes:
 89        Canonical spike representation. Row `0` stores times in ms, row `1`
 90        stores trial or unit indices.
 91    tlim:
 92        Two-element sequence defining [tmin, tmax] in ms. Defaults to inferred bounds.
 93    dt:
 94        Temporal resolution in ms.
 95
 96    Returns
 97    -------
 98    tuple[np.ndarray, np.ndarray]
 99        `(binary, time)` where `binary.shape == (n_trials, n_bins)` and
100        `time` contains the left bin edges in ms.
101
102    Notes
103    -----
104    The output is not strictly binary if multiple spikes land in the same bin.
105    In that case the matrix stores spike counts per bin.
106
107    Examples
108    --------
109    >>> spikes = np.array([[0.0, 1.0, 1.0], [0.0, 0.0, 1.0]])
110    >>> binary, time = spiketimes_to_binary(spikes, tlim=[0.0, 3.0], dt=1.0)
111    >>> binary.astype(int).tolist()
112    [[1, 1, 0], [0, 1, 0]]
113    >>> time.tolist()
114    [-0.5, 0.5, 1.5]
115    """
116    if tlim is None:
117        tlim = get_time_limits(spiketimes)
118
119    time = pylab.arange(tlim[0], tlim[1] + dt, dt).astype(float)
120    if dt <= 1:
121        time -= 0.5 * float(dt)
122
123    trials = pylab.array([-1] + list(range(int(spiketimes[1, :].max() + 1)))) + 0.5
124
125    tlim_spikes = cut_spiketimes(spiketimes, tlim)
126    tlim_spikes = tlim_spikes[:, pylab.isnan(tlim_spikes[0, :]) == False]
127
128    if tlim_spikes.shape[1] > 0:
129        binary = pylab.histogram2d(
130            tlim_spikes[0, :], tlim_spikes[1, :], [time, trials]
131        )[0].T
132    else:
133        binary = pylab.zeros((len(trials) - 1, len(time) - 1))
134    return binary, time[:-1]

198def spiketimes_to_list(spiketimes: np.ndarray) -> list[np.ndarray]:
199    """
200    Convert canonical `spiketimes` into one array per trial or unit.
201
202    Parameters
203    ----------
204    spiketimes:
205        Spike array with times in row `0` and indices in row `1`.
206
207    Returns
208    -------
209    list[np.ndarray]
210        `result[i]` contains the spike times of trial or unit `i`.
211
212    Examples
213    --------
214    >>> spiketimes_to_list(np.array([[1.0, 2.0, 4.0], [0.0, 1.0, 1.0]]))
215    [array([1.]), array([2., 4.])]
216    """
217    if spiketimes.shape[1] == 0:
218        return []
219    trials = range(int(spiketimes[1, :].max() + 1))
220    orderedspiketimes = spiketimes.copy()
221    orderedspiketimes = orderedspiketimes[
222        :, pylab.isnan(orderedspiketimes[0]) == False
223    ]
224    spike_order = pylab.argsort(orderedspiketimes[0], kind="mergesort")
225    orderedspiketimes = orderedspiketimes[:, spike_order]
226    trial_order = pylab.argsort(orderedspiketimes[1], kind="mergesort")
227    orderedspiketimes = orderedspiketimes[:, trial_order]
228
229    spikelist: list[Optional[np.ndarray]] = [None] * len(trials)
230    start = 0
231    for trial in trials:
232        end = bisect_right(orderedspiketimes[1], trial)
233        trialspikes = orderedspiketimes[0, start:end]
234        start = end
235        spikelist[trial] = trialspikes
236    return [spikes if spikes is not None else pylab.array([]) for spikes in spikelist]

 88def gaussian_kernel(sigma, dt=1.0, nstd=3.0):
 89    r"""Return a normalized Gaussian kernel for rate smoothing.
 90
 91    Parameters
 92    ----------
 93    sigma:
 94        Kernel width in ms.
 95    dt:
 96        Temporal resolution of the target grid in ms.
 97    nstd:
 98        Half-width of the kernel in units of `sigma`.
 99
100    Definition
101    ----------
102    On a grid $t_k = k\,dt$ centered at zero, this function samples the
103    truncated Gaussian
104
105    $$
106    G_{\sigma, n_\mathrm{std}}(t) =
107    \begin{cases}
108    \frac{\exp\left(-(t / \sigma)^2\right)}
109    {\int_{-n_\mathrm{std}\sigma}^{n_\mathrm{std}\sigma} \exp\left(-(u / \sigma)^2\right)\,du}
110    & \text{for } |t| \le n_\mathrm{std}\,\sigma, \\
111    0 & \text{otherwise,}
112    \end{cases}
113    $$
114
115    and rescales the sampled support so that `kernel.sum() * dt = 1`.
116
117    <img src="spiketools_assets/gaussian_kernel_example.png"
118         alt="Gaussian kernel example"
119         onerror="this.onerror=null;this.src='../spiketools_assets/gaussian_kernel_example.png';" />
120
121    Examples
122    --------
123    >>> kernel = gaussian_kernel(2.0, dt=1.0, nstd=1.0)
124    >>> round(kernel.sum(), 3)
125    1.0
126    """
127    t = pylab.arange(-nstd * sigma, nstd * sigma + dt, dt)
128    gauss = pylab.exp(-t ** 2 / sigma ** 2)
129    gauss /= gauss.sum() * dt
130    return gauss

171def kernel_rate(spiketimes, kernel, tlim=None, dt=1.0, pool=True):
172    r"""Estimate smoothed firing rates from `spiketimes`.
173
174    Parameters
175    ----------
176    spiketimes:
177        Canonical spike representation with times in ms.
178    kernel:
179        One-dimensional kernel normalized to unit integral in seconds.
180    tlim:
181        Optional `[tmin, tmax]` interval in ms.
182    dt:
183        Bin width in ms used for discretization before convolution.
184    pool:
185        If `True`, average across trials or units before convolution and return
186        a single population rate trace. If `False`, return one trace per row.
187
188    Returns
189    -------
190    tuple[np.ndarray, np.ndarray]
191        `(rates, time)` where `rates` is in spikes/s.
192
193    Definition
194    ----------
195    If `b_i[k]` is the discretized spike train on row `i`, then the returned
196    estimate is the discrete kernel convolution
197
198    $$
199    r_i[n] = 1000 \sum_k b_i[k]\,K[n-k]
200    $$
201
202    where $K$ is normalized so that $\sum_k K[k]\,dt = 1$. If `pool=True`, the
203    mean spike train across rows is convolved before returning the result.
204
205    <img src="spiketools_assets/kernel_rate_example.png"
206         alt="Kernel rate example"
207         onerror="this.onerror=null;this.src='../spiketools_assets/kernel_rate_example.png';" />
208
209    Examples
210    --------
211    >>> spikes = np.array([[0.0, 2.0], [0.0, 0.0]])
212    >>> kernel = gaussian_kernel(1.0, dt=1.0, nstd=1.0)
213    >>> rates, time = kernel_rate(spikes, kernel, tlim=[0.0, 4.0], dt=1.0, pool=False)
214    >>> rates.shape, time.shape
215    ((1, 2), (2,))
216    """
217    if tlim is None:
218        tlim = get_time_limits(spiketimes)
219
220    binary, time = spiketimes_to_binary(spiketimes, tlim, dt)
221
222    if pool:
223        binary = binary.mean(axis=0)[pylab.newaxis, :]
224
225    rates = convolve2d(binary, kernel[pylab.newaxis, :], "same")
226    kwidth = len(kernel)
227    rates = rates[:, int(kwidth / 2) : -int(kwidth / 2)]
228    time = time[int(kwidth / 2) : -int(kwidth / 2)]
229    return rates * 1000.0, time

281def rate_integral(rate, dt):
282    r"""Integrate a rate trace in spikes/s to expected spike counts.
283
284    Parameters
285    ----------
286    rate:
287        One-dimensional rate trace in spikes/s.
288    dt:
289        Sampling interval in ms.
290
291    Definition
292    ----------
293    The cumulative expected spike count is
294
295    $$
296    I[n] = \sum_{k \le n} r[k]\,dt / 1000
297    $$
298
299    which is the discrete-time version of integrating the rate over time.
300
301    <img src="spiketools_assets/rate_integral_example.png"
302         alt="Rate integral example"
303         onerror="this.onerror=null;this.src='../spiketools_assets/rate_integral_example.png';" />
304
305    Examples
306    --------
307    >>> rate_integral(np.array([500.0, 500.0]), dt=1.0).tolist()
308    [0.5, 1.0]
309    """
310    return pylab.cumsum(rate / 1000.0) * dt

232def sliding_counts(spiketimes, window, dt=1.0, tlim=None):
233    r"""Count spikes inside a sliding window.
234
235    Parameters
236    ----------
237    spiketimes:
238        Canonical spike representation.
239    window:
240        Window width in ms.
241    dt:
242        Step size of the discretized binary representation in ms.
243    tlim:
244        Optional `[tmin, tmax]` interval in ms.
245
246    Definition
247    ----------
248    With $W = \text{window} / dt$ bins and binned spike train `b_i[k]`, this function
249    returns
250
251    $$
252    C_i[n] = \sum_{k = 0}^{W - 1} b_i[n + k]
253    $$
254
255    for each row `i`.
256
257    <img src="spiketools_assets/sliding_counts_example.png"
258         alt="Sliding counts example"
259         onerror="this.onerror=null;this.src='../spiketools_assets/sliding_counts_example.png';" />
260
261    Examples
262    --------
263    >>> spikes = np.array([[0.0, 2.0], [0.0, 0.0]])
264    >>> counts, time = sliding_counts(spikes, window=2.0, dt=1.0, tlim=[0.0, 4.0])
265    >>> counts.astype(int).tolist()
266    [[1, 1, 1]]
267    """
268    if tlim is None:
269        tlim = get_time_limits(spiketimes)
270    binary, time = spiketimes_to_binary(spiketimes, dt=dt, tlim=tlim)
271
272    kernel = pylab.ones((1, int(window // dt)))
273    counts = convolve2d(binary, kernel, "valid")
274
275    dif = time.shape[0] - counts.shape[1]
276    time = time[int(np.ceil(dif / 2)) : int(-dif / 2)]
277
278    return counts, time

133def triangular_kernel(sigma, dt=1):
134    r"""Return a normalized triangular smoothing kernel.
135
136    Parameters
137    ----------
138    sigma:
139        Target width parameter in ms.
140    dt:
141        Temporal resolution of the target grid in ms.
142
143    Definition
144    ----------
145    Let $a = \sigma \sqrt{6}$. The underlying continuous kernel is
146
147    $$
148    T_\sigma(t) = \frac{\max\left(1 - |t| / a, 0\right)}
149    {\int_{-\infty}^{\infty} \max\left(1 - |u| / a, 0\right)\,du}
150    $$
151
152    and the discrete samples are normalized so that `kernel.sum() * dt = 1`.
153
154    <img src="spiketools_assets/triangular_kernel_example.png"
155         alt="Triangular kernel example"
156         onerror="this.onerror=null;this.src='../spiketools_assets/triangular_kernel_example.png';" />
157
158    Examples
159    --------
160    >>> kernel = triangular_kernel(1.0, dt=1.0)
161    >>> round(kernel.sum(), 3)
162    1.0
163    """
164    half_base = int(pylab.around(sigma * pylab.sqrt(6)))
165    half_kernel = pylab.linspace(0.0, 1.0, half_base + 1)
166    kernel = pylab.append(half_kernel, half_kernel[:-1][::-1])
167    kernel /= dt * kernel.sum()
168    return kernel

192def rate_warped_analysis(
193    spiketimes,
194    window,
195    step=1.0,
196    tlim=None,
197    rate=None,
198    func=lambda x: x.shape[1] / x[1].max(),
199    kwargs=None,
200    rate_kernel=None,
201    dt=1.0,
202):
203    r"""
204    Apply a function after warping time according to instantaneous rate.
205
206    Parameters
207    ----------
208    spiketimes:
209        Canonical spike representation.
210    window:
211        Window size in warped time. Use `"full"` to analyze the complete warped
212        recording at once.
213    step:
214        Step size in warped time.
215    tlim:
216        Optional `[tmin, tmax]` interval in ms.
217    rate:
218        Optional precomputed `(rate, time)` tuple as returned by
219        `kernel_rate(...)`.
220    func:
221        Statistic to evaluate on the warped spike train.
222    kwargs:
223        Optional keyword arguments forwarded to `func`.
224    rate_kernel:
225        Kernel used when `rate` is not provided.
226    dt:
227        Sampling interval in ms used for the rate estimate.
228
229    Definition
230    ----------
231    The idea is to replace real time $t$ by warped time $\tau(t)$,
232
233    $$
234    \tau(t) = \int_{t_0}^{t} \frac{r(u)}{1000}\,du,
235    $$
236
237    where $r(u)$ is an estimated firing rate in spikes/s. Equal distances on
238    the warped axis therefore correspond to equal expected spike counts. High-
239    rate epochs are stretched, low-rate epochs are compressed, and the chosen
240    statistic is then evaluated on this transformed spike train. If windowed
241    output is requested, the resulting window centers are mapped back to real
242    time before returning.
243
244    Examples
245    --------
246    >>> spikes = pylab.array([[0.0, 3.0, 6.0], [0.0, 0.0, 0.0]])
247    >>> result, warped_duration = rate_warped_analysis(
248    ...     spikes,
249    ...     window="full",
250    ...     func=lambda s: s.shape[1],
251    ...     rate=(pylab.array([[500.0, 500.0, 500.0]]), pylab.array([0.0, 1.0, 2.0])),
252    ...     dt=1.0,
253    ... )
254    >>> int(result), round(float(warped_duration), 3)
255    (3, 1.5)
256    """
257    if kwargs is None:
258        kwargs = {}
259    if rate_kernel is None:
260        rate_kernel = gaussian_kernel(50.0)
261    if rate is None:
262        rate = kernel_rate(spiketimes, rate_kernel, tlim)
263
264    rate, trate = rate
265    rate_tlim = [trate.min(), trate.max() + 1]
266    spiketimes = cut_spiketimes(spiketimes, rate_tlim)
267
268    ot = pylab.cumsum(rate) / 1000.0 * dt
269    w_spiketimes = spiketimes.copy()
270    w_spiketimes[0, :] = time_warp(spiketimes[0, :], trate, ot)
271
272    if window == "full":
273        return func(w_spiketimes, **kwargs), ot.max()
274
275    result, tresult = time_resolved(
276        w_spiketimes, window, func, kwargs, tstep=step
277    )
278    tresult = time_warp(tresult, ot, trate)
279    return result, tresult

 52def time_resolved(
 53    spiketimes,
 54    window,
 55    func: Callable,
 56    kwargs=None,
 57    tlim: Optional[Sequence[float]] = None,
 58    tstep=1.0,
 59):
 60    r"""
 61    Apply a function to successive windows of a spike train.
 62
 63    Parameters
 64    ----------
 65    spiketimes:
 66        Canonical spike representation.
 67    window:
 68        Window size in ms.
 69    func:
 70        Callable receiving the cropped `spiketimes` of each window.
 71    kwargs:
 72        Optional keyword arguments passed to `func`.
 73    tlim:
 74        Optional analysis interval `[tmin, tmax]` in ms.
 75    tstep:
 76        Window step size in ms.
 77
 78    Returns
 79    -------
 80    tuple[list, np.ndarray]
 81        Function outputs and window-center times.
 82
 83    Definition
 84    ----------
 85    If $S$ denotes the spike train and $f$ the analysis function, then
 86
 87    $$
 88    y_n = f\left(S \cap [t_n, t_n + W)\right),
 89    \qquad
 90    t_n^\mathrm{center} = t_n + \frac{W}{2},
 91    $$
 92
 93    where $W$ is the window width and consecutive windows are offset by
 94    `tstep`.
 95
 96    Examples
 97    --------
 98    >>> spikes = pylab.array([[0.0, 2.0, 4.0], [0.0, 0.0, 0.0]])
 99    >>> values, times = time_resolved(spikes, window=2.0, func=lambda s: s.shape[1], tlim=[0.0, 5.0], tstep=1.0)
100    >>> values
101    [1, 1, 1, 1]
102    >>> times.tolist()
103    [0.5, 1.5, 2.5, 3.5]
104    """
105    if kwargs is None:
106        kwargs = {}
107    if tlim is None:
108        tlim = get_time_limits(spiketimes)
109    cut_spikes = cut_spiketimes(spiketimes, tlim)
110
111    tmin = tlim[0]
112    tmax = tmin + window
113    tcenter = tmin + 0.5 * (tmax - tstep - tmin)
114
115    time = []
116    func_out = []
117    while tmax <= tlim[1]:
118        windowspikes = cut_spiketimes(cut_spikes, [tmin, tmax])
119        time.append(tcenter)
120        func_out.append(func(windowspikes, **kwargs))
121
122        tmin += tstep
123        tmax += tstep
124        tcenter += tstep
125
126    return func_out, pylab.array(time)

129def time_resolved_new(
130    spiketimes,
131    window,
132    func: Callable,
133    kwargs=None,
134    tlim: Optional[Sequence[float]] = None,
135    tstep=1.0,
136):
137    r"""
138    Windowed evaluation using an intermediate binary representation.
139
140    This variant can be faster for repeated window extraction because the spike
141    train is binned once before the loop.
142
143    Definition
144    ----------
145    This function evaluates the same windowed statistic as `time_resolved(...)`,
146    but it first bins the spike train into a count matrix $B$ and reconstructs
147    each window from the relevant slice:
148
149    Let $k_{n} = \lfloor t_{n} \rfloor$ be the left window index on the binned
150    grid and let $K = \lfloor W \rfloor$ be the corresponding window width in
151    bins. Let $B_{n}$ denote the corresponding slice of the count matrix and let
152    $\mathcal{C}$ denote the conversion from a binary window back to spike
153    times. The statistic is then
154
155    $$
156    y_n = f\left(\mathcal{C}(B_{n})\right).
157    $$
158
159    Examples
160    --------
161    >>> spikes = pylab.array([[0.0, 2.0, 4.0], [0.0, 0.0, 0.0]])
162    >>> values, times = time_resolved_new(spikes, window=2.0, func=lambda s: s.shape[1], tlim=[0.0, 5.0], tstep=1.0)
163    >>> values
164    [1, 1, 1, 1]
165    """
166    if kwargs is None:
167        kwargs = {}
168    if tlim is None:
169        tlim = get_time_limits(spiketimes)
170
171    binary, btime = spiketimes_to_binary(spiketimes, tlim=tlim)
172
173    tmin = tlim[0]
174    tmax = tmin + window
175    tcenter = 0.5 * (tmax + tmin)
176
177    time = []
178    func_out = []
179    while tmax <= tlim[1]:
180        windowspikes = binary_to_spiketimes(
181            binary[:, int(tmin) : int(tmax)], btime[int(tmin) : int(tmax)]
182        )
183        time.append(tcenter)
184        func_out.append(func(windowspikes, **kwargs))
185
186        tmin += tstep
187        tmax += tstep
188        tcenter += tstep
189    return func_out, pylab.array(time)

181def cv2(
182    spiketimes,
183    pool=True,
184    return_all=False,
185    bessel_correction=False,
186    minvals=0,
187):
188    r"""Compute the coefficient-of-variation statistic from inter-spike intervals.
189
190    Parameters
191    ----------
192    spiketimes:
193        Canonical spike representation.
194    pool:
195        If `True`, pool intervals across trials or units before computing the
196        statistic. Otherwise compute per row and average.
197    return_all:
198        If `True`, return per-row values instead of their mean.
199    bessel_correction:
200        Use `ddof=1` in the interval variance.
201    minvals:
202        Minimum number of finite intervals required for a per-row value.
203
204    Definition
205    ----------
206    Despite the historical function name, this returns the squared coefficient
207    of variation of the inter-spike intervals $I_n$:
208
209    $$
210    \mathrm{CV}^2 = \frac{\mathrm{Var}[I_n]}{\mathrm{E}[I_n]^2}
211    $$
212
213    Examples
214    --------
215    >>> round(float(cv2(np.array([[0.0, 2.0, 5.0], [0.0, 0.0, 0.0]]))), 3)
216    0.04
217    """
218    if spiketimes.shape[1] < 3:
219        if return_all:
220            return pylab.array([pylab.nan])
221        return pylab.nan
222    ddof = 1 if bessel_correction else 0
223    spikelist = spiketimes_to_list(spiketimes)
224    maxlen = max(len(sl) for sl in spikelist)
225    if maxlen < 3:
226        return pylab.nan
227    spikearray = pylab.zeros((len(spikelist), maxlen)) * pylab.nan
228    spike_counts = []
229    for i, sl in enumerate(spikelist):
230        spikearray[i, : len(sl)] = sl
231        spike_counts.append(len(sl))
232    spike_counts = pylab.array(spike_counts)
233    spikearray = spikearray[spike_counts > 2]
234    intervals = pylab.diff(spikearray, axis=1)
235
236    if pool:
237        var = pylab.nanvar(intervals, ddof=ddof)
238        mean_squared = pylab.nanmean(intervals) ** 2
239        return var / mean_squared
240
241    intervals[pylab.isfinite(intervals).sum(axis=1) < minvals, :] = pylab.nan
242    var = pylab.nanvar(intervals, axis=1, ddof=ddof)
243    mean_squared = pylab.nanmean(intervals, axis=1) ** 2
244    if return_all:
245        return var / mean_squared
246    return pylab.nanmean(var / mean_squared)

266def cv_two(spiketimes, min_vals=20):
267    r"""Compute the local Cv2 measure from consecutive inter-spike intervals.
268
269    Definition
270    ----------
271    For consecutive inter-spike intervals $I_n$ and $I_{n+1}$, the local Cv2
272    statistic is
273
274    $$
275    \mathrm{Cv2} =
276    \left\langle
277    \frac{2\,|I_{n+1} - I_n|}{I_{n+1} + I_n}
278    \right\rangle_n
279    $$
280
281    where the average is taken across all finite neighboring interval pairs.
282
283    Examples
284    --------
285    >>> round(float(cv_two(np.array([[0.0, 2.0, 5.0, 9.0], [0.0, 0.0, 0.0, 0.0]]), min_vals=2)), 3)
286    0.343
287    """
288    consecutive_isis = _consecutive_intervals(spiketimes)
289    ms = (
290        2
291        * pylab.absolute(consecutive_isis[:, 0] - consecutive_isis[:, 1])
292        / (consecutive_isis[:, 0] + consecutive_isis[:, 1])
293    )
294    if len(ms) >= min_vals:
295        return pylab.mean(ms)
296    return pylab.nan

 99def ff(spiketimes, mintrials=None, tlim=None):
100    r"""Compute the Fano factor across trials or units.
101
102    Parameters
103    ----------
104    spiketimes:
105        Canonical spike representation.
106    mintrials:
107        Optional minimum number of rows required for a valid estimate.
108    tlim:
109        Optional `[tmin, tmax]` counting window in ms.
110
111    Definition
112    ----------
113    If `N` is the spike count in the analysis window, the Fano factor is
114
115    $$
116    \mathrm{FF} = \frac{\mathrm{Var}[N]}{\mathrm{E}[N]}
117    $$
118
119    Examples
120    --------
121    >>> ff(np.array([[0.0, 1.0, 0.0, 1.0], [0.0, 0.0, 1.0, 1.0]]), tlim=[0.0, 2.0])
122    0.0
123    """
124    if tlim is None:
125        tlim = get_time_limits(spiketimes)
126    dt = tlim[1] - tlim[0]
127    binary, _ = spiketimes_to_binary(spiketimes, tlim=tlim, dt=dt)
128    counts = binary.sum(axis=1)
129    if (counts == 0).all():
130        return pylab.nan
131    if mintrials is not None and len(counts) < mintrials:
132        return pylab.nan
133    return counts.var() / counts.mean()

136def kernel_fano(spiketimes, window, dt=1.0, tlim=None, components=False):
137    """Estimate a sliding-window Fano factor from binned spike counts.
138
139    Parameters
140    ----------
141    spiketimes:
142        Canonical spike representation.
143    window:
144        Counting window in ms.
145    dt:
146        Bin width in ms.
147    tlim:
148        Optional analysis interval in ms.
149    components:
150        If `True`, return variance and mean separately instead of their ratio.
151
152    Examples
153    --------
154    >>> fano, time = kernel_fano(
155    ...     np.array([[0.0, 2.0, 0.0, 2.0], [0.0, 0.0, 1.0, 1.0]]),
156    ...     window=2.0,
157    ...     dt=1.0,
158    ...     tlim=[0.0, 4.0],
159    ... )
160    >>> fano.shape, time.shape
161    ((3,), (3,))
162    """
163    if tlim is None:
164        tlim = get_time_limits(spiketimes)
165    if tlim[1] - tlim[0] == window:
166        binary, time = spiketimes_to_binary(spiketimes, tlim=tlim)
167        counts = binary.sum(axis=1)
168        return pylab.array([counts.var() / counts.mean()]), pylab.array(time.mean())
169
170    counts, time = sliding_counts(spiketimes, window, dt, tlim)
171
172    vars_ = counts.var(axis=0)
173    means = counts.mean(axis=0)
174    fano = pylab.zeros_like(vars_) * pylab.nan
175    if components:
176        return vars_, means, time
177    fano[means > 0] = vars_[means > 0] / means[means > 0]
178    return fano, time