Metrics

Scalar Metrics

tad.metrics.scalar.firing_rate(r, *, channels=None, tstart=None, tstop=None, inclusive_stop=False, per_channel=True)[source]

Compute firing rate FR = N / T, as in the references.

Parameters:

  • r (Raster) – Input raster.

  • channels (Sequence[int | str] | None) – Channels to include. If None, uses all channels.

  • tstart (float | None) – Optional time window. If None, inferred from data.

  • tstop (float | None) – Optional time window. If None, inferred from data.

  • inclusive_stop (bool) – If False, use [tstart, tstop); if True, use [tstart, tstop].

  • per_channel (bool) – If True, return FR per channel; otherwise return pooled FR.

Returns:

fr – Firing rate per channel (Hz) or pooled rate (Hz).

Return type:

ndarray or float

Raises:

ValueError – If window duration is non-positive.

tad.metrics.scalar.mean_firing_rate_across_channels(r, *, channels=None, tstart=None, tstop=None, inclusive_stop=False, active_threshold_hz=None)[source]

Compute mean firing rate across channels (optionally over “active channels”).

This matches the typical workflow described in the references where per-channel rates are computed and an “active channel” criterion may be applied.

Parameters:

  • r (Raster) – Input raster.

  • channels (Sequence[int | str] | None) – Channels to include.

  • tstart (float | None) – Optional time window.

  • tstop (float | None) – Optional time window.

  • inclusive_stop (bool) – If False, use [tstart, tstop); if True, use [tstart, tstop].

  • active_threshold_hz (float | None) – If provided, only channels with FR >= threshold are included in the mean.

Returns:

Mean firing rate (Hz) over selected channels.

Return type:

float

Raises:

ValueError – If no channels remain after thresholding.

tad.metrics.scalar.mean_inter_event_interval(r, *, channels=None, tstart=None, tstop=None, inclusive_stop=False)[source]

Compute the mean inter-event interval (IEI) of pooled activity.

This corresponds to the pooled IEI analysis described in Pasquale et al. 2008 (IEIs computed from the globally ordered spike times across electrodes).

Parameters:

  • r (Raster) – Input raster.

  • channels (Sequence[int | str] | None) – Channels to include. If None, uses all channels.

  • tstart (float | None) – Optional time window.

  • tstop (float | None) – Optional time window.

  • inclusive_stop (bool) – Window right-bound policy.

Returns:

Mean IEI. Returns np.nan if fewer than 2 pooled events exist.

Return type:

float

tad.metrics.scalar.percent_random_spiking(r, burst_intervals, *, channels=None, tstart=None, tstop=None, inclusive_stop=False)[source]

Compute percentage of random spiking activity (%RS), i.e. fraction of spikes outside bursts.

In Pasquale et al. 2008, “% random spiking activity” is defined as the fraction of spikes that do not belong to bursts.

Parameters:

  • r (Raster) – Input raster.

  • burst_intervals (Dict[int | str, ndarray]) – Mapping channel -> array of burst intervals with shape (B, 2), where each row is (burst_start, burst_end). Intervals are assumed to be within the analysis window.

  • channels (Sequence[int | str] | None) – Channels to include. If None, uses all raster channels.

  • tstart (float | None) – Optional analysis window.

  • tstop (float | None) – Optional analysis window.

  • inclusive_stop (bool) – Window right-bound policy.

Returns:

Percentage (0..1) of spikes outside bursts. Returns np.nan if no spikes exist.

Return type:

float

Notes

This function assumes bursts are already detected elsewhere (next step in your pipeline).

tad.metrics.scalar.save_scalar_metric(list_of_metrics, fname, output_folder='results/')[source]

Function to save the list_of_metrics passed as an argument here to a json file

Parameters:

  • list_of_metrics (Dict[str, int | float | List[float] | List[int]]) – A dictionary containing the metrics to be saved.

  • fname (str) – The name of the file to save the metrics to (without extension).

  • output_folder (str) – The folder to save the results in (default “results/”).

tad.metrics.scalar.spike_count(r, *, channels=None, tstart=None, tstop=None, inclusive_stop=False, per_channel=True)[source]

Count spikes in a window.

This is the base quantity used by firing rate metrics in the references (FR = N / T).

Parameters:

  • r (Raster) – Input raster.

  • channels (Sequence[int | str] | None) – Channels to include. If None, uses all channels.

  • tstart (float | None) – Optional time window. If None, inferred from data.

  • tstop (float | None) – Optional time window. If None, inferred from data.

  • inclusive_stop (bool) – If False, use [tstart, tstop); if True, use [tstart, tstop].

  • per_channel (bool) – If True, return an array of counts per channel; otherwise return total count.

Returns:

counts – Spike counts per channel (shape (n_channels,)) or pooled total count.

Return type:

ndarray or int

Firing Rate

class tad.metrics.rates.FiringRateCurveResult(t, fr_ch, fr_pop, dt, tstart, tstop, channels, bin_edges)[source]

Bases: object

Result of a binned firing-rate curve computation.

Parameters:

  • t (numpy.ndarray) – Bin times according to time_mode (centers by default), shape (n_bins,).

  • fr_ch (numpy.ndarray) – Per-channel firing rate per bin (Hz), shape (n_channels, n_bins).

  • fr_pop (numpy.ndarray) – Population firing rate per bin (Hz), i.e. sum over channels, shape (n_bins,).

  • dt (float) – Bin width (seconds or your time unit).

  • tstart (float) – Analysis window start (actual binning start).

  • tstop (float) – Analysis window stop (actual binning stop = last edge).

  • channels (List[str | int]) – Channel IDs corresponding to fr_ch rows.

  • bin_edges (numpy.ndarray) – Bin edges used, shape (n_bins+1,).

bin_edges: ndarray
channels: List[str | int]
dt: float
fr_ch: ndarray
fr_pop: ndarray
t: ndarray
tstart: float
tstop: float
tad.metrics.rates.firing_rate_curve(r, dt, *, tstart=None, tstop=None, channels=None, inclusive_stop=False, time_mode='center')[source]

Compute binned firing-rate curves from a Raster.

This is the natural “time-resolved” extension of FR = N/T from the references: counts in bins divided by bin width.

Parameters:

  • r (Raster) – Input raster.

  • dt (float) – Bin width.

  • tstart (float | None) – Optional time window. If None, inferred by Raster.bin_counts.

  • tstop (float | None) – Optional time window. If None, inferred by Raster.bin_counts.

  • channels (Sequence[int | str] | None) – Optional subset/order of channels.

  • inclusive_stop (bool) – If True, include events exactly at tstop.

  • time_mode (str) – How to represent time points: - “left”: left bin edges (starts), shape (n_bins,) - “center”: bin centers, shape (n_bins,) (default)

Returns:

Contains per-channel and population firing-rate curves.

Return type:

FiringRateCurveResult

Raises:

ValueError – If time_mode is not recognized.

Avalanches

class tad.metrics.avalanches.AvalancheResult(sizes, lifetimes, dt, tstart, tstop, channels, size_definition, bin_edges, active_bins, intervals_bins)[source]

Bases: object

Result of avalanche extraction on a raster.

Parameters:

  • sizes (numpy.ndarray) – Avalanche sizes, one per avalanche. Definition depends on size_definition.

  • lifetimes (numpy.ndarray) – Avalanche lifetimes in number of bins (Δt units), one per avalanche.

  • dt (float) – Bin width used to define avalanches.

  • tstart (float) – Start time of analysis window used for binning.

  • tstop (float) – Stop time of analysis window used for binning.

  • channels (List[str | int]) – Channel IDs included in the analysis (row order of internal matrices).

  • size_definition (int) – Size definition used: - 1: sum over bins of (# active electrodes in bin) - 2: number of distinct electrodes active at least once in avalanche

  • bin_edges (numpy.ndarray) – Bin edges used, shape (n_bins+1,).

  • active_bins (numpy.ndarray) – Boolean array indicating which bins were active (at least one spike in any channel), shape (n_bins,).

  • intervals_bins (List[Tuple[int, int]]) – List of (start_bin, stop_bin) bin indices for each avalanche, where stop_bin is exclusive.

active_bins: ndarray
bin_edges: ndarray
channels: List[str | int]
dt: float
intervals_bins: List[Tuple[int, int]]
lifetimes: ndarray
size_definition: int
sizes: ndarray
tstart: float
tstop: float
tad.metrics.avalanches.extract_avalanches(r, dt, *, tstart=None, tstop=None, channels=None, inclusive_stop=False, size_definition=1)[source]

Extract neuronal avalanches (sizes and lifetimes) from a raster.

This follows the definitions in the provided references:

  • Bin time into windows of width Δt.

  • A bin is “active” if at least one electrode fires ≥1 spike in that bin; “blank” otherwise.

  • An avalanche is a maximal contiguous sequence of active bins bounded by blank bins.

  • Lifetime = number of consecutive active bins.

  • Size:

    1. sum over bins of number of active electrodes in that bin (multiple activations counted)

    2. number of distinct electrodes that were active at least once during the avalanche

Parameters:

  • r (Raster) – Input raster.

  • dt (float) – Bin width (Δt). Must be positive.

  • tstart (float | None) – Optional analysis window. If None, inferred from the data in the selected channels. (Inference follows your Raster.bin_counts behavior.)

  • tstop (float | None) – Optional analysis window. If None, inferred from the data in the selected channels. (Inference follows your Raster.bin_counts behavior.)

  • channels (Sequence[int | str] | None) – Channels to include. If None, uses all channels in the raster.

  • inclusive_stop (bool) – Whether to include events at exactly tstop in binning window.

  • size_definition (int) – 1 or 2, as above.

Returns:

Extracted avalanche sizes, lifetimes, and supporting metadata.

Return type:

AvalancheResult

Raises:

ValueError – If size_definition is not 1 or 2.

Bursts

class tad.metrics.burst.Burst(start, end, n_spikes, duration, intra_rate_hz)[source]

Bases: object

Single burst description.

Parameters:

  • start (float) – Burst start time (time of first spike in burst).

  • end (float) – Burst end time (time of last spike in burst).

  • n_spikes (int) – Number of spikes in the burst.

  • duration (float) – Burst duration defined as end - start.

  • intra_rate_hz (float) – Intra-burst firing rate (Hz), defined here as (n_spikes - 1) / duration if duration > 0, else np.inf.

duration: float
end: float
intra_rate_hz: float
n_spikes: int
start: float
class tad.metrics.burst.BurstChannelResult(channel, isi_th, bursts, diagnostics=None)[source]

Bases: object

Bursts for one channel.

Parameters:

  • channel (int | str) – Channel id.

  • isi_th (float) – ISI threshold used for this channel.

  • bursts (List[tad.metrics.burst.Burst]) – List of detected bursts.

  • diagnostics (tad.metrics.burst.LogISIThresholdDiagnostics | None) – Threshold-selection diagnostics if method=”logisih”, otherwise None. If threshold_scope=”pooled”, this may be the same object for all channels.

bursts: List[Burst]
channel: int | str
diagnostics: LogISIThresholdDiagnostics | None = None
isi_th: float
class tad.metrics.burst.BurstDetectionResult(per_channel, method, tstart, tstop, channels)[source]

Bases: object

Burst detection result across channels.

Parameters:

  • per_channel (Dict[int | str, tad.metrics.burst.BurstChannelResult]) – Mapping channel -> BurstChannelResult.

  • method (str) – A human-readable description of the method used.

  • tstart (float) – Time window used for detection.

  • tstop (float) – Time window used for detection.

  • channels (List[str | int]) – Channel IDs included (order).

channels: List[str | int]
method: str
per_channel: Dict[int | str, BurstChannelResult]
tstart: float
tstop: float
class tad.metrics.burst.LogISIThresholdDiagnostics(isi_th, centers_log10, hist, hist_smooth, peak1_idx, peak2_idx, valley_idx, method)[source]

Bases: object

Diagnostics for adaptive ISI threshold selection using the log-ISIH.

Parameters:

  • isi_th (float) – Chosen ISI threshold (seconds).

  • centers_log10 (numpy.ndarray) – Histogram bin centers in log10(ISI) domain.

  • hist (numpy.ndarray) – Raw histogram values (counts or density depending on call).

  • hist_smooth (numpy.ndarray) – Smoothed histogram used for peak/valley detection.

  • peak1_idx (int) – Index of the left peak (short-ISI regime) in centers_log10, or -1.

  • peak2_idx (int) – Index of the right peak (long-ISI regime) in centers_log10, or -1.

  • valley_idx (int) – Index of the valley between peaks used as threshold, or -1 when fallback used.

  • method (str) – Selection path, e.g. “valley_between_peaks” or “fallback_quantile_*”.

centers_log10: ndarray
hist: ndarray
hist_smooth: ndarray
isi_th: float
method: str
peak1_idx: int
peak2_idx: int
valley_idx: int
tad.metrics.burst.choose_isi_threshold_logisih(isi_values=None, r=None, channel=None, tstart=None, tstop=None, inclusive_stop=False, *, bins='pasquale', smooth_window='from_bins', density=False, min_isi=1e-06, fallback=0.1)[source]

Choose an ISI threshold from the log10(ISI) histogram (log-ISIH).

This is an adaptive threshold selection strategy commonly used for burst detection in MEA analyses: short ISIs (within bursts) and long ISIs (between bursts) can create a structured distribution in log10(ISI). When two regimes are detectable, the threshold is taken as the valley between the two main peaks.

Algorithm (robust, NumPy-only): 1) Compute y = log10(ISI) for ISI > 0. 2) Compute histogram of y with bins bins. 3) Smooth histogram with moving average window smooth_window using 20% of bins as window. 4) Detect local maxima (peaks) on smoothed histogram. 5) If two sufficiently separated peaks exist, pick the valley (minimum) between them

and set ISI_th = 10**(center_at_valley).

  1. Otherwise fallback to a conservative quantile of the ISI distribution.

Parameters:

  • isi_values (ndarray | None) – 1D array of ISIs (>0).

  • bins (str | int) – Number of bins in log10(ISI) space, or method from np.histogram.

  • smooth_window (str | int) – Moving average smoothing window length, or method from number of bins.

  • density (bool) – If True, histogram is density; otherwise counts.

  • min_isi (float) – Minimum ISI value before log transform (prevents log(0)).

  • fallback (float) – ISI used if peak/valley detection fails.

  • r (Raster | None)

  • channel (int | str | None)

  • tstart (float | None)

  • tstop (float | None)

  • inclusive_stop (bool)

Returns:

Contains chosen ISI_th and intermediate arrays/indices for inspection.

Return type:

LogISIThresholdDiagnostics

tad.metrics.burst.detect_bursts(r, *, method='fixed', isi_th=0.1, min_spikes=3, channels=None, tstart=None, tstop=None, inclusive_stop=False, threshold_scope='per_channel', logisih_bins='auto', logisih_smooth_window='from_bins', fallback=0.1)[source]

Detect single-channel bursts in a Raster.

Three threshold-selection methods are supported:

  1. method=”fixed”

    Uses the constant threshold isi_th for all channels.

  2. method=”logisih”

    Derives the ISI threshold from the log10(ISI) histogram (log-ISIH) either: - per channel (threshold_scope=”per_channel”), or - from pooled ISIs across selected channels (threshold_scope=”pooled”).

  3. method=”fixed_from_logisih”

    Uses pre-calculated ISI thresholds from log-ISIH analysis. For this method, isi_th must be an array with one threshold value per channel in the same order as channels.

Burst definition (segmentation): A burst is a maximal consecutive sequence of spikes such that every consecutive inter-spike interval satisfies ISI <= ISI_th, and the sequence contains at least min_spikes spikes.

Parameters:

  • r (Raster) – Input raster.

  • method (Literal['fixed', 'logisih', 'fixed_from_logisih']) – Threshold selection method: “fixed”, “logisih”, or “fixed_from_logisih”.

  • isi_th (float | ndarray) – ISI threshold(s) in seconds. For “fixed” and “logisih”, a single float. For “fixed_from_logisih”, an array of floats with one per channel.

  • min_spikes (int) – Minimum number of spikes required to accept a burst.

  • channels (Sequence[int | str] | None) – Channels to include. If None, uses all channels in the raster.

  • tstart (float | None) – Optional time window. If None, inferred from data across selected channels.

  • tstop (float | None) – Optional time window. If None, inferred from data across selected channels.

  • inclusive_stop (bool) – If True, include events exactly at tstop; otherwise window is [tstart, tstop).

  • threshold_scope (Literal['per_channel', 'pooled']) – Only used when method=”logisih”: - “per_channel”: compute a separate ISI_th per channel from that channel’s ISIs - “pooled”: compute a single ISI_th from pooled ISIs across channels and reuse it

  • logisih_bins (str | int) – Number of bins in log10(ISI) histogram, or method for numpy.histogram (“auto”, “fd”, “doane”, “sqrt” etc) (method=”logisih”).

  • logisih_smooth_window (str | int) – Smoothing window length for histogram, or method to determine it from number of bins (method=”logisih”).

  • fallback (float) – ISI threshold value used when log-ISIH peak/valley detection fails.

Returns:

Per-channel burst lists and metadata. For method=”logisih”, each channel result includes a diagnostics object that can be plotted/inspected; when threshold_scope=”pooled”, this diagnostic may be shared across channels.

Return type:

BurstDetectionResult

Raises:

ValueError – If parameters are inconsistent.

ISI Analysis

class tad.metrics.isi.ISIResult(isi, mode, tstart, tstop, channels)[source]

Bases: object

Result container for ISI extraction.

Parameters:

  • isi (Dict[int | str, numpy.ndarray] | numpy.ndarray) – ISIs. If mode=”per_channel”, this is a dict {channel_id: isi_array}. If mode=”pooled”, this is a single 1D array.

  • mode (Literal['per_channel', 'pooled']) – “per_channel” or “pooled”.

  • tstart (float) – Window used for extraction.

  • tstop (float) – Window used for extraction.

  • channels (List[str | int]) – Channel IDs included in extraction (order).

channels: List[str | int]
isi: Dict[int | str, ndarray] | ndarray
mode: Literal['per_channel', 'pooled']
tstart: float
tstop: float
tad.metrics.isi.isi(r, *, channels=None, tstart=None, tstop=None, inclusive_stop=False, mode='per_channel')[source]

Compute inter-spike intervals (ISIs) from a raster.

Parameters:

  • r (Raster) – Input raster.

  • channels (Sequence[int | str] | None) – Channels to include. If None, uses all channels.

  • tstart (float | None) – Optional analysis window. If None, inferred from data.

  • tstop (float | None) – Optional analysis window. If None, inferred from data.

  • inclusive_stop (bool) – Right-bound policy: [tstart,tstop) if False, [tstart,tstop] if True.

  • mode (Literal['per_channel', 'pooled']) –

    • “per_channel”: compute ISIs separately for each channel

    • ”pooled”: pool spike times across channels (global ordering) then compute ISIs

Returns:

ISIs and metadata.

Return type:

ISIResult

Notes

  • Per-channel ISIs are computed from each channel’s sorted spike times.

  • Pooled ISIs correspond to the IEI concept used in the references when pooling across electrodes (global event train).

tad.metrics.isi.isih(isi_values, *, bins=50, density=False, log=False)[source]

Compute an ISI histogram (ISIH), optionally in log domain.

Parameters:

  • isi_values (ndarray) – 1D array of ISIs (>0).

  • bins (int | ndarray) – Number of bins or explicit bin edges (passed to np.histogram).

  • density (bool) – If True, return a probability density (area = 1). If False, return counts.

  • log (bool) – If True, histogram log10(ISI) rather than ISI.

Returns:

  • centers (ndarray) – Bin centers (in ISI units if log=False; in log10 units if log=True).

  • hist (ndarray) – Histogram values (counts or density).

Return type:

Tuple[ndarray, ndarray]

Synchrony

class tad.metrics.synchrony.PearsonSynchronyResult(corr, channels, method, n_time, masked_channels)[source]

Bases: object

Pearson correlation-based synchrony result computed from firing-rate curves.

Parameters:

  • corr (numpy.ndarray) – Pearson correlation matrix, shape (n_channels, n_channels).

  • channels (List) – Channel IDs corresponding to corr rows/cols (same order as input FR).

  • method (str) – Description of computation.

  • n_time (int) – Number of time bins used.

  • masked_channels (List[Tuple[int, object]]) – Channels removed due to zero variance (if drop_constant=True). Each entry is (index_in_input, channel_id).

channels: List
corr: ndarray
masked_channels: List[Tuple[int, object]]
method: str
n_time: int
tad.metrics.synchrony.pearson_corr_firing_rate(fr, *, channels=None, zscore=True, drop_constant=True, ddof=0)[source]

Compute Pearson correlation matrix between channels using firing-rate curves.

This is intended as an exploratory synchrony metric: - Take per-channel rate signals r_i(t) (from binning) - Optionally z-score each channel across time - Compute corr(i,j) = corr(r_i(t), r_j(t))

Parameters:

  • fr (FiringRateCurveResult) – Firing rate curves result (from tad.metrics.rates.firing_rate_curve).

  • channels (Sequence[int] | None) – Optional subset of channels specified as row indices into fr.fr_ch. If None, uses all channels in stored order.

  • zscore (bool) – If True, z-score each channel across time before correlation. If False, correlation is still Pearson (mean-centered, variance-normalized), but keeping this option helps debugging.

  • drop_constant (bool) – If True, drop channels with (near-)zero temporal variance to avoid NaNs. If False, those channels will produce NaNs in correlation.

  • ddof (int) – Degrees of freedom for std calculation in z-scoring (0 is fine for signals).

Returns:

Correlation matrix and metadata.

Return type:

PearsonSynchronyResult

Notes

  • Correlation is computed as (Z @ Z.T) / T where Z is per-channel standardized data.

  • If zscore=True, Z has mean 0 and std 1 (per channel), up to numerical tolerance.

Cross-Correlation

class tad.metrics.crosscorrelation.CrossCorrelationResult[source]

Bases: object

tad.metrics.crosscorrelation.compute_crosscorrelation()[source]

Power Law Fitting

class tad.metrics.powerlaw.PowerLawFitResult(variable, a, b, rmse, x_all, p_all, keep_mask, x_used, p_used, p_hat_used, exclude_unit, tail_fraction)[source]

Bases: object

Power-law fit result for discrete avalanche distributions.

Parameters:

  • variable (Literal['size', 'lifetime']) – “size” or “lifetime”.

  • a (float) – Prefactor in P(x) = a * x^b (in probability space).

  • b (float) – Exponent in P(x) = a * x^b.

  • rmse (float) – Root mean squared error between empirical P(x) and fitted P_hat(x), computed over the retained bins (after exclusions).

  • x_all (numpy.ndarray) – All support values (sorted unique x) from the empirical PMF.

  • p_all (numpy.ndarray) – Empirical PMF values corresponding to x_all.

  • keep_mask (numpy.ndarray) – Boolean mask over x_all/p_all indicating which points were used for fit/RMSE.

  • x_used (numpy.ndarray) – x values used for fit.

  • p_used (numpy.ndarray) – empirical P(x) used for fit.

  • p_hat_used (numpy.ndarray) – fitted P_hat(x) evaluated at x_used.

  • exclude_unit (bool) – Whether x=1 was excluded.

  • tail_fraction (float) – Tail cutoff fraction relative to max(P). Points with P < tail_fraction * max(P) were excluded.

a: float
b: float
exclude_unit: bool
keep_mask: ndarray
p_all: ndarray
p_hat_used: ndarray
p_used: ndarray
rmse: float
tail_fraction: float
variable: Literal['size', 'lifetime']
x_all: ndarray
x_used: ndarray
tad.metrics.powerlaw.fit_avalanche_powerlaw(aval, *, variable='size', exclude_unit=True, tail_fraction=0.01, min_points=5)[source]

Fit a power-law model P(x) = a x^b to an avalanche distribution and compute RMSE.

This implements the paper-style protocol used in your references: - Build a discrete histogram / PMF of sizes or lifetimes - Exclude unit events (x=1) (first bin) - Exclude tail bins whose probability is < (tail_fraction * max_probability)

(paper uses 1% of the maximum => tail_fraction=0.01)

  • Fit P(x) = a x^b

  • Compute RMSE between empirical P(x) and fitted P_hat(x) on retained points

Parameters:

  • aval (AvalancheResult) – AvalancheResult containing extracted avalanche sizes/lifetimes.

  • variable (Literal['size', 'lifetime']) – Choose which distribution to fit: “size” or “lifetime”.

  • exclude_unit (bool) – If True, exclude x=1.

  • tail_fraction (float) – Tail probability cutoff relative to max(P). Default 0.01.

  • min_points (int) – Minimum number of retained points required to fit. If fewer, returns NaNs.

Returns:

Fit parameters and diagnostics.

Return type:

PowerLawFitResult

Notes

The fit is performed by linear least squares on log space:

log P = log a + b log x

Then a is recovered by exponentiation.

RMSE is computed in probability space:

RMSE = sqrt(mean((P - P_hat)^2))

over the retained points, matching the “goodness-of-fit via RMSE” style described.

Evoked Responses

class tad.metrics.evoked.EvokedPeakResult(channels, baseline_hz, peak_hz, peak_latency_s, peak_minus_baseline_hz, auc_above_baseline)[source]

Bases: object

Scalar evoked-response metrics derived from a PSTH.

Parameters:

  • channels (List[str | int]) – Channel IDs corresponding to rows.

  • baseline_hz (numpy.ndarray) – Baseline rate in Hz (mean over baseline window).

  • peak_hz (numpy.ndarray) – Peak rate in Hz (max over response window).

  • peak_latency_s (numpy.ndarray) – Time of peak in seconds (relative to stimulus).

  • peak_minus_baseline_hz (numpy.ndarray) – Peak above baseline (Hz).

  • auc_above_baseline (numpy.ndarray) – Area above baseline over response window (spikes per stimulus), computed as sum(max(rate-baseline,0))*dt.

auc_above_baseline: ndarray
baseline_hz: ndarray
channels: List[str | int]
peak_hz: ndarray
peak_latency_s: ndarray
peak_minus_baseline_hz: ndarray
tad.metrics.evoked.evoked_peak_metrics(psth, *, baseline_window=(-0.05, 0.0), response_window=(0.0, 0.05), rectify=True)[source]

Compute per-channel evoked-response scalar metrics from PSTH.

Parameters:

  • psth (PSTHResult) – PSTHResult as returned by compute_psth.

  • baseline_window (Tuple[float, float]) – (t0, t1) window for baseline mean rate (seconds, relative to stimulus).

  • response_window (Tuple[float, float]) – (t0, t1) window where peak/latency/AUC are computed.

  • rectify (bool) – If True, AUC uses max(rate-baseline, 0). If False, uses signed (rate-baseline).

Returns:

Per-channel scalar metrics.

Return type:

EvokedPeakResult

Raises:

ValueError – If windows are invalid or contain no bins.

tad.metrics.evoked.response_probability(r, stim_times, *, t0=0.0, t1=0.05, channels=None, tstart=None, tstop=None, inclusive_stop=False, inclusive_zero=True)[source]

Estimate per-channel response probability to a stimulus.

Definition used here: For each channel c and each stimulus i, define an indicator:

I_{c,i} = 1 if channel c has >=1 spike in [stim_i + t0, stim_i + t1)

else 0

Response probability per channel:

p_c = (1/N) * sum_i I_{c,i}

This complements PSTH peak metrics and is a standard evoked-response quantity.

Parameters:

  • r (Raster) – Input raster.

  • stim_times (Sequence[float]) – Stimulus onset times.

  • t0 (float) – Response window bounds relative to stimulus (seconds), with [t0, t1).

  • t1 (float) – Response window bounds relative to stimulus (seconds), with [t0, t1).

  • channels (Sequence[int | str] | None) – Optional channels subset/order.

  • tstart (float | None) – Optional admissible stimulus window; stimuli are kept only if the full [stim+t0, stim+t1) lies within [tstart, tstop) (or inclusive_stop).

  • tstop (float | None) – Optional admissible stimulus window; stimuli are kept only if the full [stim+t0, stim+t1) lies within [tstart, tstop) (or inclusive_stop).

  • inclusive_stop (bool) – Right-bound policy for admissible stimulus window.

  • inclusive_zero (bool) – If True, include spikes exactly at stim+t0 when t0==0 (uses searchsorted side).

Returns:

channels_list is list of ChannelId, p is ndarray shape (n_channels,) with probabilities.

Return type:

channels_list, p

PSTH (Post-Stimulus Time Histogram)

class tad.metrics.psth.PSTHResult(t, bin_edges, counts, rate_hz, dt, t_pre, t_post, stim_times_used, channels)[source]

Bases: object

Post-Stimulus Time Histogram (PSTH) result.

The PSTH is computed by aligning spikes to each stimulus time, binning the relative spike times in a window [-t_pre, t_post), and averaging across stimuli.

Parameters:

  • t (numpy.ndarray) – Bin times (centers if time_mode=”center”, else left edges), shape (n_bins,).

  • bin_edges (numpy.ndarray) – Bin edges in relative time (seconds), shape (n_bins+1,).

  • counts (numpy.ndarray) – Mean spike counts per bin per stimulus, shape (n_channels, n_bins).

  • rate_hz (numpy.ndarray) – Mean firing rate per bin (Hz), i.e. counts / dt, shape (n_channels, n_bins).

  • dt (float) – Bin width (seconds).

  • t_pre (float) – Pre-stimulus window length (seconds).

  • t_post (float) – Post-stimulus window length (seconds).

  • stim_times_used (numpy.ndarray) – Stimulus times actually used (after filtering), shape (n_stim_used,).

  • channels (List[str | int]) – Channel IDs corresponding to counts rows.

bin_edges: ndarray
channels: List[str | int]
counts: ndarray
dt: float
property population_counts: ndarray

Population PSTH in counts/bin/stimulus (sum across channels).

Return type:

ndarray, shape (n_bins,)

property population_rate_hz: ndarray

Population PSTH in Hz (sum across channels).

Return type:

ndarray, shape (n_bins,)

rate_hz: ndarray
stim_times_used: ndarray
t: ndarray
t_post: float
t_pre: float
tad.metrics.psth.compute_psth(r, stim_times, *, dt, t_pre, t_post, channels=None, tstart=None, tstop=None, inclusive_stop=False, inclusive_zero=True, time_mode='center', dtype=dtype('float64'))[source]

Compute PSTH aligned to stimulus onsets.

This implements the paper definition:

P(t) = (1/N) * sum_i ST(t - t_i^stim)

in discrete form by binning relative spike times per stimulus and averaging.

Parameters:

  • r (Raster) – Input raster.

  • stim_times (Sequence[float]) – Iterable of stimulus onset times (seconds).

  • dt (float) – Bin width (seconds).

  • t_pre (float) – Pre-stimulus window length (seconds). Window starts at -t_pre.

  • t_post (float) – Post-stimulus window length (seconds). Window ends at +t_post (exclusive).

  • channels (Sequence[int | str] | None) – Optional subset/order of channels to include.

  • tstart (float | None) – Optional global time window for admissible stimuli. If None, inferred from raster. A stimulus is used only if its full PSTH window lies within [tstart, tstop) (or [tstart, tstop] if inclusive_stop=True).

  • tstop (float | None) – Optional global time window for admissible stimuli. If None, inferred from raster. A stimulus is used only if its full PSTH window lies within [tstart, tstop) (or [tstart, tstop] if inclusive_stop=True).

  • inclusive_stop (bool) – Right-bound policy for the admissible stimulus window.

  • inclusive_zero (bool) – If True, spikes exactly at stimulus time (relative time 0) are included. This affects the left/right searchsorted side used for slicing.

  • time_mode (str) –

    • “center”: return bin centers in t

    • ”left”: return left bin edges in t

  • dtype (dtype) – Floating dtype for output arrays.

Returns:

PSTH counts and rates per channel.

Return type:

PSTHResult

Raises:

ValueError – If parameters are invalid.

Utilities

tad.metrics.utils.pooled_spike_times(r, *, channels=None, tstart=None, tstop=None, inclusive_stop=False)[source]

Pool spike times across channels into a single sorted array.

Parameters:

  • r (Raster) – Input raster.

  • channels (Sequence[int | str] | None) – Channels to include. If None, uses all channels.

  • tstart (float | None) – Optional time window.

  • tstop (float | None) – Optional time window.

  • inclusive_stop (bool) – Window right-bound policy: [tstart,tstop) if False, [tstart,tstop] if True.

Returns:

Sorted pooled times in the requested window.

Return type:

ndarray, shape (N,)