"""Spike-history basis objects and discrete-time filter banks.
This module mirrors MATLAB ``History.m``. History bases turn a vector
of window boundary times into a set of piecewise-constant basis
functions that GLMs use to capture refractoriness and bursting.
Exported symbols
----------------
- :class:`History` — primary spike-history basis (Matlab ``History``).
- :class:`HistoryFilter` — immutable discrete-time MATLAB-style
transfer-function pair ``(numerator, denominator)``.
- :class:`HistoryFilterBank` — collection of :class:`HistoryFilter`
objects, returned by :meth:`History.computeHistory`.
Time vectors and window boundaries are in **seconds**; sample rates in
**Hz**. A history window with bounds ``[w_i, w_{i+1}]`` (in seconds)
contributes a basis column equal to ``1`` for time lags in that interval
and ``0`` elsewhere — exactly the MATLAB convention.
"""
from __future__ import annotations
from dataclasses import dataclass
from collections.abc import Sequence
from typing import Any
import matplotlib.pyplot as plt
import numpy as np
from .core import Covariate, SignalObj, nspikeTrain
@dataclass(frozen=True)
class HistoryFilter:
"""Discrete-time MATLAB-style transfer function numerator/denominator pair."""
numerator: np.ndarray
denominator: np.ndarray
delta: float
variable: str = "z^-1"
@dataclass(frozen=True)
class HistoryFilterBank:
"""Matrix-like collection of discrete history-window filters."""
numerators: tuple[np.ndarray, ...]
denominators: tuple[np.ndarray, ...]
delta: float
variable: str = "z^-1"
@property
def shape(self) -> tuple[int, int]:
return (len(self.numerators), 1)
@property
def numFilters(self) -> int:
return len(self.numerators)
def __len__(self) -> int:
return self.numFilters
def __getitem__(self, index: int) -> HistoryFilter:
return HistoryFilter(
numerator=np.asarray(self.numerators[index], dtype=float).copy(),
denominator=np.asarray(self.denominators[index], dtype=float).copy(),
delta=float(self.delta),
variable=self.variable,
)
def combine(self, coefficients) -> HistoryFilter:
coeffs = np.asarray(coefficients, dtype=float).reshape(-1)
if coeffs.size != self.numFilters:
raise ValueError("Number of coefficients must match the number of history filters.")
max_len = max(len(numerator) for numerator in self.numerators)
padded = np.zeros((self.numFilters, max_len), dtype=float)
for idx, numerator in enumerate(self.numerators):
arr = np.asarray(numerator, dtype=float).reshape(-1)
padded[idx, : arr.size] = arr
numerator = coeffs @ padded
denominator = np.zeros(max_len, dtype=float)
denominator[0] = 1.0
return HistoryFilter(numerator=np.asarray(numerator, dtype=float), denominator=denominator, delta=float(self.delta), variable=self.variable)
def __rmatmul__(self, coefficients) -> HistoryFilter:
return self.combine(coefficients)
[docs]
class History:
"""Spike-history basis defined by a set of window boundaries.
Given a vector of *N* window-times, *N − 1* piecewise-constant basis
functions are created. Window *i* spans from ``windowTimes[i]`` to
``windowTimes[i+1]``.
Parameters
----------
windowTimes : array_like
Strictly increasing vector of window boundary times (seconds).
Must contain at least two entries.
minTime : float or None
Optional lower time bound for the computed history covariate.
maxTime : float or None
Optional upper time bound for the computed history covariate.
name : str, default ``'History'``
Display name.
Examples
--------
>>> h = History([0, 0.001, 0.01, 0.05, 0.1])
>>> h.numWindows
4
See Also
--------
TrialConfig, Analysis.computeHistLag
"""
[docs]
def __init__(self, windowTimes, minTime: float | None = None, maxTime: float | None = None, name: str = "History") -> None:
times = np.asarray(windowTimes, dtype=float).reshape(-1)
if times.size <= 1:
raise ValueError("At least two times points must be specified to determine a window")
if np.any(np.diff(times) <= 0):
raise ValueError("windowTimes must be strictly increasing")
self.windowTimes = times
self.minTime = None if minTime is None else float(minTime)
self.maxTime = None if maxTime is None else float(maxTime)
self.name = str(name)
# parity: python-extension — convenience accessor (no MATLAB counterpart)
@property
def lags(self) -> np.ndarray:
return np.asarray(self.windowTimes[1:], dtype=float).copy()
# parity: python-extension — convenience accessor (no MATLAB counterpart)
@property
def numWindows(self) -> int:
return int(self.windowTimes.size - 1)
[docs]
def computeHistory(self, trains, historyIndex: int | None = None, time_grid=None):
"""Matlab-facing alias for :meth:`compute_history`."""
return self.compute_history(trains, historyIndex, time_grid=time_grid)
[docs]
def setWindow(self, windowTimes) -> None:
"""Replace the window-times vector.
Parameters
----------
windowTimes : array_like
New window boundary times (must be strictly increasing,
length ≥ 2).
"""
replacement = History(windowTimes, self.minTime, self.maxTime, self.name)
self.windowTimes = replacement.windowTimes
self.minTime = replacement.minTime
self.maxTime = replacement.maxTime
[docs]
def plot(self, *_, handle=None, **__):
"""Plot each of the history windows as step functions.
Parameters
----------
handle : Axes or None
Matplotlib axes to draw into. If ``None``, a new figure is
created.
Returns
-------
Axes or list
The axes or plot handles.
"""
tmin = np.asarray(self.windowTimes[:-1], dtype=float)
tmax = np.asarray(self.windowTimes[1:], dtype=float)
sampleRate = 1000.0
num_samples = max(1, int(round((float(np.max(tmax)) - float(np.min(tmin))) * sampleRate)))
data = np.zeros((num_samples, tmax.size), dtype=float)
dataLabels: list[str] = []
for index, (start, stop) in enumerate(zip(tmin, tmax)):
indMin = max(1, int(round((float(start) - float(np.min(tmin))) * sampleRate)))
indMax = int(round((float(stop) - float(np.min(tmin))) * sampleRate))
if indMax >= indMin:
data[indMin - 1 : indMax, index] = 1.0
dataLabels.append(f"[{start:.3g},{stop:.3g}]")
time = np.linspace(float(np.min(tmin)), float(np.max(tmax)), num_samples)
signal = SignalObj(time, data, "History", "time", "s", "", dataLabels)
created_ax = handle is None
ax = handle if handle is not None else plt.subplots(1, 1, figsize=(6.0, 2.2))[1]
plot_handles = signal.plot(handle=ax)
return ax if created_ax else plot_handles
[docs]
def toFilter(self, delta: float) -> HistoryFilterBank:
"""Convert the history windows to a discrete-time filter bank.
Parameters
----------
delta : float
Sample period (seconds). Must be positive.
Returns
-------
HistoryFilterBank
Bank of FIR filters (one per window) in the ``z⁻¹`` domain.
"""
delta = float(delta)
if delta <= 0:
raise ValueError("delta must be positive")
tmin = np.asarray(self.windowTimes[:-1], dtype=float)
tmax = np.asarray(self.windowTimes[1:], dtype=float)
numerators: list[np.ndarray] = []
denominators: list[np.ndarray] = []
for row, (window_start, window_stop) in enumerate(zip(tmin, tmax)):
num_samples = int(np.ceil(float(window_stop) / delta))
start_sample = int(np.ceil(float(window_start) / delta)) + 1
del row
# ``toFilter`` numerator is ``size = num_samples + 1`` with
# a leading zero — MATLAB-faithful per ``History.m:165-168``
# ``b(i,(StartSample:NumSamples)+1)=1; %delay by 1``. The
# ``computeHistory`` / ``_compute_single_history`` path uses
# ``size = num_samples`` and applies its own ``[0, 1]`` delay
# filter externally, so the two construction paths produce
# the same final convolution with 1-sample shift offsets
# that cancel in the output. Both match the MATLAB gold
# fixture ``history_exactness.mat``. Audit finding H2:
# verified harmless; do NOT "unify" by dropping the leading
# zero — that breaks ``test_history_matches_matlab_gold_fixture``.
numerator = np.zeros(num_samples + 1, dtype=float)
denominator = np.zeros(num_samples + 1, dtype=float)
denominator[0] = 1.0
numerator[start_sample : num_samples + 1] = 1.0
numerators.append(numerator)
denominators.append(denominator)
return HistoryFilterBank(numerators=tuple(numerators), denominators=tuple(denominators), delta=delta)
def _compute_single_history(self, train: nspikeTrain, historyIndex: int | None = None, time_grid=None) -> Covariate:
sigrep = train.getSigRep() if time_grid is None else train.getSigRep(None, float(np.min(time_grid)), float(np.max(time_grid)))
tmin = np.asarray(self.windowTimes[:-1], dtype=float)
tmax = np.asarray(self.windowTimes[1:], dtype=float)
data_columns: list[np.ndarray] = []
data_labels: list[str] = []
for window_start, window_stop in zip(tmin, tmax, strict=False):
num_samples = int(np.ceil(float(window_stop) * float(train.sampleRate)))
numerator = np.zeros(max(num_samples, 0), dtype=float)
start_sample = int(np.ceil(float(window_start) * float(train.sampleRate))) + 1
if num_samples > 0 and start_sample <= num_samples:
numerator[max(start_sample - 1, 0) : num_samples] = 1.0
filtered = sigrep.filter(numerator if numerator.size else [0.0], [1.0])
delayed = filtered.filter([0.0, 1.0], [1.0])
data_columns.append(np.asarray(delayed.dataToMatrix(), dtype=float))
if historyIndex is None:
data_labels.append(f"[{window_start:.3g},{window_stop:.3g}]")
else:
data_labels.append(f"[{window_start:.3g},{window_stop:.3g}]_{historyIndex}")
data = np.hstack(data_columns) if data_columns else np.zeros((sigrep.time.size, 0), dtype=float)
name = "History" if not getattr(train, "name", "") else f"History {train.name}"
cov = Covariate(sigrep.time, data, name, sigrep.xlabelval, sigrep.xunits, sigrep.yunits, data_labels)
if time_grid is not None:
return cov
if data.size == 0:
return Covariate([], data, name, sigrep.xlabelval, sigrep.xunits, sigrep.yunits, data_labels)
if (self.minTime is not None or self.maxTime is not None) and round(float(cov.sampleRate), 9) != round(float(train.sampleRate), 9):
cov.resampleMe(float(train.sampleRate))
min_time = float(cov.minTime) if self.minTime is None else float(self.minTime)
max_time = float(cov.maxTime) if self.maxTime is None else float(self.maxTime)
windowed = cov.getSigInTimeWindow(min_time, max_time)
windowed.setMinTime(float(train.minTime))
windowed.setMaxTime(float(train.maxTime))
windowed.minTime = float(train.minTime)
windowed.maxTime = float(train.maxTime)
return windowed
# parity: python-extension — snake_case implementation; MATLAB API is `computeHistory`
[docs]
def compute_history(self, trains, historyIndex: int | None = None, time_grid=None):
"""Compute the history covariate(s) for one or more spike trains.
Parameters
----------
trains : nspikeTrain or nstColl or sequence of nspikeTrain
Spike train(s) whose history to compute.
historyIndex : int or None
Optional label index appended to the data labels.
time_grid : array_like or None
Optional external time grid. If ``None``, the spike train's
own time grid is used.
Returns
-------
CovariateCollection
One :class:`Covariate` per spike train, each with
``numWindows`` columns corresponding to the history windows.
"""
from .trial import CovariateCollection
if isinstance(trains, nspikeTrain):
cov = self._compute_single_history(trains, historyIndex, time_grid=time_grid)
if historyIndex is not None:
cov.name = f"History #{historyIndex} for {trains.name}"
return CovariateCollection([cov])
if hasattr(trains, "getNST") and hasattr(trains, "numSpikeTrains"):
covariates = [
self._compute_single_history(trains.getNST(index), historyIndex, time_grid=time_grid)
for index in range(int(trains.numSpikeTrains))
]
return CovariateCollection(covariates)
if isinstance(trains, Sequence) and not isinstance(trains, (str, bytes, np.ndarray)):
covariates = [self._compute_single_history(train, historyIndex, time_grid=time_grid) for train in trains]
return CovariateCollection(covariates)
raise TypeError("History can only be computed from nspikeTrain, nstColl, or sequences of nspikeTrain")
[docs]
def toStructure(self) -> dict[str, Any]:
"""Serialise the History to a plain dictionary."""
return {
"windowTimes": self.windowTimes.tolist(),
"minTime": self.minTime,
"maxTime": self.maxTime,
"name": self.name,
}
[docs]
@staticmethod
def fromStructure(structure: dict[str, Any] | None) -> "History" | None:
"""Reconstruct a History from a dictionary (inverse of :meth:`toStructure`)."""
if structure is None:
return None
if "windowTimes" in structure:
windowTimes = structure["windowTimes"]
elif "lags" in structure:
lags = np.asarray(structure["lags"], dtype=float).reshape(-1)
windowTimes = np.concatenate([[0.0], lags])
else:
windowTimes = [0.0, 1.0]
return History(
windowTimes,
minTime=structure.get("minTime"),
maxTime=structure.get("maxTime"),
name=structure.get("name", "History"),
)
[docs]
@classmethod
def raisedCosine(cls, K: int, tMin: float | None = None, tMax: float | None = None) -> "History":
"""Construct a ``History`` with ``K`` log-spaced windows in ``[tMin, tMax]``.
Windowed-basis analog of the Pillow et al. 2008 raised-cosine basis on
log time — log-spaced lag scales capture refractory effects at short
lags (small bins) and slow adaptation at long lags (large bins).
Parameters
----------
K : int
Number of history windows (typical: 5–10). Must be ``>= 2``.
tMin : float, optional
Shortest lag in seconds. Default ``0.002`` (2 ms — absolute
refractory edge).
tMax : float, optional
Longest lag in seconds. Default ``0.100`` (100 ms — typical
slow-adaptation horizon).
Returns
-------
History
History object whose ``windowTimes`` is ``[0, peaks...]`` with
``K`` log-spaced peaks in ``[tMin, tMax]``. The leading ``0``
defines the first window from spike-onset to ``tMin``.
Notes
-----
Returns the WINDOWED-basis approximation of the Pillow 2008
raised-cosine basis. The true Pillow basis is a set of smooth cosine
bumps; reproducing those requires a separate ``RaisedCosineBasis``
class with overloaded ``computeHistory``. For PP-GLM history-coefficient
fitting where the basis is collapsed into bin counts the log-spaced
windowed approximation is sufficient.
Mirrors MATLAB ``History.raisedCosine`` (History.m:197).
"""
if tMin is None:
tMin = 0.002
if tMax is None:
tMax = 0.100
K = int(K)
if K < 2:
raise ValueError(f"K must be >= 2 (need at least 2 log-spaced peaks); got {K}.")
if tMin <= 0:
raise ValueError(f"tMin must be > 0 for log spacing (got {tMin}).")
if tMax <= tMin:
raise ValueError(f"tMax ({tMax}) must exceed tMin ({tMin}).")
peaks = np.logspace(np.log10(float(tMin)), np.log10(float(tMax)), K)
windowTimes = np.concatenate(([0.0], peaks))
return cls(windowTimes)
HistoryBasis = History
__all__ = ["History", "HistoryBasis", "HistoryFilter", "HistoryFilterBank"]