SpikeSliceStack

The SpikeSliceStack class stores event-aligned slices of raw spike train data. It holds a list of SpikeData objects (one per trial or event window) and provides methods for converting to raster arrays of shape (num_units, num_bins, num_slices) and computing trial-averaged statistics.

class spikelab.spikedata.spikeslicestack.SpikeSliceStack(data_obj=None, times_start_to_end=None, time_peaks=None, time_bounds=None, spike_stack=None, neuron_attributes=None, drop_slice_attributes=True)

Bases: object

A list of SpikeData objects, one per slice, with spike-based comparison capabilities.

U is units (neurons) and S is slices (bursts, events, etc). Construct from either a single SpikeData with time specifications, or directly from a pre-built list of SpikeData objects.

Parameters:

• data_obj (SpikeData or None) – A SpikeData object to slice. Provide either this or spike_stack, not both.

• times_start_to_end (list or None) – Each entry is a tuple (start, end) representing the start and end times of a desired slice. Each tuple must have the same duration.

• time_peaks (list or None) – List of times as int or float where there is a burst peak or stimulation event. Must be paired with time_bounds. Alternative to times_start_to_end.

• time_bounds (tuple or None) – Single tuple (left_bound, right_bound). For example, (250, 500) means 250 ms before peak and 500 ms after peak. Must be paired with time_peaks.

• spike_stack (list or None) – List of SpikeData objects, one per slice. All must have the same number of units. Spike times must be relative to the slice (0-based or event-centered via start_time), not absolute recording times. Provide either this or data_obj, not both.

• neuron_attributes (list or None) – List of attribute dicts, one per unit. If None, inherited from data_obj when available.

• drop_slice_attributes (bool) – If True (default), neuron_attributes are removed from individual SpikeData slices after construction. The shared copy is stored at neuron_attributes. This avoids duplicating large per-unit data (e.g. waveform templates) across every slice. Set to False to keep per-slice attributes.

spike_stack

List of SpikeData objects, one per slice. Spike times are relative to the slice window. For 0-based slices, times run from 0 to duration. For event-centered slices, times run from -pre_ms to +post_ms with t=0 at the event. Use self.times for absolute recording time positions.

Type:

list

times

List of (start, end) time bounds for each slice in absolute recording time, sorted chronologically. Length equals S. Example: [(100, 350), (500, 750), (1000, 1250)].

Type:

list

N

Number of units.

Type:

int

neuron_attributes

List of attribute dicts, one per unit. None if not provided. When drop_slice_attributes is True (default), this is the only copy and individual slices will have neuron_attributes set to None.

Type:

list or None

__init__(data_obj=None, times_start_to_end=None, time_peaks=None, time_bounds=None, spike_stack=None, neuron_attributes=None, drop_slice_attributes=True)

subslice(slices)

Extract a subset of slices from the spike stack.

Parameters:

slices (int or list) – Slice index or list of slice indices to extract.

Returns:

New SpikeSliceStack containing only the

specified slices. Shape changes from S to S_trimmed. All units and neuron_attributes are carried over.

Return type:

result (SpikeSliceStack)

subset(units, by=None, preserve_order=False)

Extract a subset of units from every slice in the spike stack.

Parameters:

• units (int, str, or list) – Unit indices to extract. If by is None, must be int(s). If by is set, values to match in neuron_attributes.

• by (str or None) – If set, select units by this neuron_attribute key instead of by index.

• preserve_order (bool) – When False (default), output is sorted ascending by index — consistent with the other SpikeLab data classes. When True, output respects the order of the input units list. Duplicates are deduplicated either way.

Returns:

New SpikeSliceStack containing only the

specified units across all slices. All slices and neuron_attributes are carried over.

Return type:

result (SpikeSliceStack)

Notes

• If IDs are not unique (when using by), every matching neuron is included.

subtime_by_index(start_idx, end_idx)

Trim each slice to a sub-window specified by millisecond indices.

Indices are measured from the start of each slice (1 index = 1 ms). Trims along the time axis while preserving all slices and units.

Parameters:

• start_idx (int) – Start index in ms from slice start (inclusive). Supports negative indexing.

• end_idx (int) – End index in ms from slice start (exclusive). Supports negative indexing.

Returns:

New SpikeSliceStack where each slice is

trimmed to the corresponding absolute time window. Absolute spike times are preserved (not shifted). self.times is updated to reflect the new absolute time bounds.

Return type:

result (SpikeSliceStack)

Raises:

ValueError – If the underlying slice duration (times[0][1] -     times[0][0]) is not an integer number of milliseconds. Use SpikeData.subtime() with explicit ms bounds for non-integer windows.

Notes

• Indices are relative to each slice’s own start (index 0 = slice start ms). They are converted to absolute recording times internally before trimming.

• Original absolute timestamps are preserved. To get shifted-to-zero timestamps, create a new SpikeSliceStack.

• All slices and neuron_attributes are carried over from the original.

to_raster_array(bin_size=1.0, absolute_times=False)

Convert the spike stack into a 3D raster array of shape (N, T, S).

Each slice is rasterized with the given bin size, producing a spike count matrix where entry (n, t, s) is the number of spikes unit n fired in time bin t of slice s.

Parameters:

• bin_size (float) – Time bin size in ms (default 1.0).

• absolute_times (bool) – If False (default), time bin 0 corresponds to the start of each slice (0-based). If True, each slice’s spikes are offset by its absolute start time from self.times, so bin indices reflect the original recording position. The T dimension is sized to cover the full time span from the earliest slice start to the latest slice end. Caution: this can produce very large arrays when the recording span is long and bin_size is small.

Returns:

3D array of shape (N, T, S) with

non-negative integer spike counts. When absolute_times is True, T covers the full recording span and all slices share the same time axis.

Return type:

raster_stack (np.ndarray)

baseline_normalized_raster(bin_size, baseline_window_ms, *, mode='subtract')

Per-slice raster normalized against a per-slice baseline rate.

Wraps to_raster_array(bin_size) and converts each bin into a baseline-normalized response value. The baseline rate is computed from spikes inside baseline_window_ms (in milliseconds relative to each slice’s time origin) and projected to each bin via rate * bin_size. Output shape matches the raster: (U, T, S).

Parameters:

• bin_size (float) – Raster bin size in milliseconds. Passed to to_raster_array.

• baseline_window_ms (tuple[float, float]) – (start_ms, end_ms) window relative to slice origin used to estimate the per-slice baseline rate.

• mode (str) –

  Normalization mode: - "subtract" (default) — counts above baseline expectation. - "ratio" — counts / expected_counts (NaN where expected

  is 0).

  • "zscore" — (counts - expected) / sqrt(expected), the Poisson z-score (NaN where expected is 0).

Returns:

Float array of shape (U, T, S).

Return type:

normalized (np.ndarray)

Notes

• Baseline window is validated against each slice’s actual time range. ValueError if any slice doesn’t contain it.

• For uniform-bin response counts (no normalization), use to_raster_array(bin_size) directly; this method adds the per-slice baseline correction on top.

responsive_units(bin_size, baseline_window_ms, *, response_window_ms=None, z_threshold=2.0, aggregator='mean')

Identify units that show a significant evoked response.

Builds the Poisson-z-scored baseline-normalized raster, optionally restricts to a response time window, aggregates across slices (mean or max), and returns a unit mask where any time bin’s aggregated z-score exceeds z_threshold.

Parameters:

• bin_size (float) – Raster bin size in milliseconds.

• baseline_window_ms (tuple[float, float]) – Baseline window (start_ms, end_ms) relative to slice origin used to estimate the per-slice baseline rate.

• response_window_ms (tuple[float, float] or None) – Optional response window (start_ms, end_ms) relative to slice origin. When None (default), the full slice is searched.

• z_threshold (float) – Z-score threshold (default 2.0).

• aggregator (str) – How to combine z-scores across slices before thresholding. "mean" (default) or "max".

Returns:

Boolean array of shape (U,). True for

responsive units.

Return type:

mask (np.ndarray)

Notes

• Units with no baseline spikes in any slice are flagged non-responsive (z-scores are NaN there).

decode_slice_labels(labels, response_window_ms, *, bin_size, baseline_window_ms=None, classifier='ridge', cv='loo', classifier_kwargs=None, random_state=None)

Decode per-slice labels (e.g. stim identity) from population responses.

Builds an (S, U) feature matrix by summing per-unit spike counts in response_window_ms (optionally with baseline subtraction) and runs cross-validated classifier decoding via spikelab.spikedata.decoding.cross_validated_decode().

Parameters:

• labels (array-like) – Per-slice labels of length S (e.g. stim electrode index, treatment category).

• response_window_ms (tuple[float, float]) – Window relative to slice origin over which response counts are summed.

• bin_size (float) – Raster bin size in ms.

• baseline_window_ms (tuple[float, float] or None) – Optional per-slice baseline window; when provided, counts are baseline-subtracted via baseline_normalized_raster.

• classifier (str) – "ridge" (default), "mlp", or "random_forest".

• cv (str or int) – "loo" (default) or int >= 2.

• classifier_kwargs (dict or None) – Forwarded to the sklearn classifier constructor.

• random_state (int or None) – Reproducibility seed.

Returns:

Same shape as

spikelab.spikedata.decoding.cross_validated_decode() — accuracy, predictions, true_labels, confusion_matrix, per_fold_accuracy, classes, classifier_name.

Return type:

result (dict)

Notes

• Requires scikit-learn (optional dependency).

• For decoding from the full (U, T) raster (not just summed counts), call decoding.cross_validated_decode directly on self.to_raster_array(bin_size).reshape(U * T, S).T.

group_pair_similarity(stim_labels, *, metric='cosine', bin_size=1.0, slice_indices=None)

Pairwise similarity between mean response vectors for each stim class.

For each unique stimulus label, averages the per-slice (U, T) raster across all slices that share that label, then computes a (K, K) similarity matrix between the resulting per-stim mean response vectors. Lets you ask: “how distinguishable are responses to different stims?”.

Parameters:

• stim_labels (array-like) – Per-slice stim label of length S.

• metric (str) – "cosine" (default), "pearson", "euclidean" (distance), or "cross_entropy".

• bin_size (float) – Raster bin size in ms (default 1.0).

• slice_indices (array-like or None) – Optional subset of slice indices to use (e.g. an “early-cycle” or “late-cycle” window). When None (default), uses all slices.

Returns:

(K, K) similarity matrix with

labels set to the unique stim classes (sorted).

Return type:

sim (PairwiseCompMatrix)

Notes

• Stim classes that have no slices in slice_indices are dropped from the output.

responsive_units_per_group(group_labels, bin_size, baseline_window_ms, *, response_window_ms=None, z_threshold=2.0, aggregator='mean')

Per-cycle responsive-unit mask for tracking responsiveness over time.

For each unique cycle, runs responsive_units on the slices belonging to that cycle and returns a (U, n_cycles) boolean matrix. Use the per-cycle masks to compute gained / lost / preserved responsive units across cycle groups, or to correlate responsiveness changes with intrinsic activity changes per unit.

Parameters:

• group_labels (array-like) – Per-slice cycle index of length S.

• bin_size (float) – Raster bin size in ms.

• baseline_window_ms (tuple[float, float]) – Baseline window for Poisson z-score normalization.

• response_window_ms (tuple[float, float] or None) – Optional response window (default: full slice).

• z_threshold (float) – Per-cycle z-threshold (default 2.0).

• aggregator (str) – "mean" (default) or "max" across slices within each cycle.

Returns:

• cycles (np.ndarray): Sorted unique cycle indices.

• mask (np.ndarray): (U, n_cycles) boolean responsiveness mask.

• responsive_count (np.ndarray): Per-cycle responsive unit count, shape (n_cycles,).

Return type:

result (dict)

responsiveness_change(group_labels, early_groups, late_groups, bin_size, baseline_window_ms, *, response_window_ms=None, z_threshold=2.0, aggregator='mean')

Gained / lost / preserved responsive units between two cycle groups.

Computes responsive-unit masks separately for slices in early_groups and late_groups (any iterables of cycle indices), and reports which units become responsive (“gained”), stop being responsive (“lost”), or stay responsive (“preserved”).

Parameters:

• group_labels (array-like) – Per-slice cycle index of length S.

• early_groups (array-like) – Cycle indices for the early group.

• late_groups (array-like) – Cycle indices for the late group.

• bin_size (float) – Raster bin size in ms.

• baseline_window_ms (tuple[float, float]) – Baseline window.

• response_window_ms (tuple[float, float] or None) – Response window (default: full slice).

• z_threshold (float) – Per-group z-threshold.

• aggregator (str) – "mean" (default) or "max".

Returns:

• early_mask (np.ndarray (U,) bool): Responsive in early.

• late_mask (np.ndarray (U,) bool): Responsive in late.

• gained (np.ndarray (U,) bool): NOT responsive in early AND responsive in late.

• lost (np.ndarray (U,) bool): Responsive in early AND NOT responsive in late.

• preserved (np.ndarray (U,) bool): Responsive in BOTH.

• early_count, late_count, gained_count, lost_count, preserved_count (int).

Return type:

result (dict)

Notes

• Pair this with intrinsic-activity changes per unit (e.g. cv_isi differences) and correlate via stat_utils.linear_regression to ask whether responsiveness changes track changes in baseline activity.

slice_to_slice_similarity(metric='cosine', *, bin_size=1.0)

Pairwise similarity between slice-wise population response vectors.

Each slice is converted to a (U * T) flat vector via to_raster_array(bin_size).reshape(U*T, S).T and a square (S, S) similarity matrix is computed using the requested metric.

Parameters:

• metric (str) – One of "cosine" (default), "pearson", "euclidean" (distance), or "cross_entropy" (symmetric KL on normalized bin distributions).

• bin_size (float) – Raster bin size in ms (default 1.0).

Returns:

(S, S) similarity matrix. For

cosine and pearson, higher = more similar (diagonal ~1.0); for euclidean and cross_entropy, lower = more similar (diagonal 0).

Return type:

sim (PairwiseCompMatrix)

Notes

• cosine and pearson return values in [-1, 1].

• euclidean returns raw L2 distance.

• cross_entropy returns symmetric KL divergence (i.e. (KL(p||q) + KL(q||p)) / 2) between bin distributions normalized to sum to 1.

• Use PairwiseCompMatrix.extract_lower_triangle() for feature extraction.

per_unit_response_regression(bin_size, response_window_ms, *, x_values=None, baseline_window_ms=None, min_valid_slices=3)

Per-unit OLS regression of evoked response amplitude across slices.

For each slice and unit, computes response amplitude as the sum of spike counts in response_window_ms — optionally with a per-slice baseline subtraction. Then fits a linear regression of amplitude against x_values for every unit. Use this to detect facilitation / depression of the evoked response across cycles or stimulus intensities.

Parameters:

• bin_size (float) – Raster bin size in milliseconds (passed to to_raster_array).

• response_window_ms (tuple[float, float]) – (start_ms, end_ms) window relative to slice origin over which response counts are summed.

• x_values (array-like or None) – Per-slice x values for the regression (e.g. cycle index, stimulus intensity). Length must equal the number of slices S. When None (default), uses np.arange(S).

• baseline_window_ms (tuple[float, float] or None) – Optional baseline window for per-slice subtraction. When None (default), uses raw response counts; otherwise subtracts the expected count per bin (baseline_rate * bin_size) before summing.

• min_valid_slices (int) – Minimum number of valid (non-NaN) (x, y) pairs required to fit a regression. Units with fewer return NaN for all coefficients. Default 3.

Returns:

Dictionary with keys:

• slope (np.ndarray (U,)): Slope per unit.

• intercept (np.ndarray (U,)): Intercept per unit.

• r_squared (np.ndarray (U,)): R² per unit.

• p_value (np.ndarray (U,)): Two-sided p-value of the slope per unit.

• stderr (np.ndarray (U,)): Standard error of the slope per unit.

• amplitudes (np.ndarray (U, S)): Per-slice response amplitude (raw or baseline-subtracted).

• x_values (np.ndarray (S,)): The x values used.

Return type:

result (dict)

Notes

• Requires scipy (optional dependency); raises ImportError with installation instructions if missing.

• Units with constant amplitudes get r_squared = 0, slope = 0 and p_value = 1.0.

compute_frac_active(min_spikes=2)

Compute the fraction of slices each unit is active in.

A unit counts as active in a slice if it has at least min_spikes spikes within that slice’s time window.

Parameters:

min_spikes (int) – Minimum number of spikes for a unit to count as active in a slice (default: 2).

Returns:

1-D array of shape (U,) with the

fraction of slices each unit is active in (values in [0, 1]).

Return type:

frac_active (np.ndarray)

Notes

• The returned array can be passed as frac_active to RateSliceStack.order_units_across_slices, RateSliceStack.get_slice_to_slice_unit_corr_from_stack, SpikeSliceStack.order_units_across_slices, or SpikeSliceStack.get_slice_to_slice_unit_comparison to override their internal activity calculation.

• SpikeData.get_frac_active produces a compatible (U,) array based on burst edges and can be used in the same way.

order_units_across_slices(agg_func='median', timing='median', min_spikes=2, min_frac_active=0.0, frac_active=None, timing_matrix=None)

Reorder units by their typical spike timing across slices.

For each unit in each slice, computes a representative spike time (median, mean, or first spike) relative to the slice’s time origin. These per-slice values are aggregated across slices to obtain a single typical timing per unit. Units are then sorted by this value from earliest to latest and optionally split into a highly-active group and a low-activity group.

Parameters:

• agg_func (str) – How to aggregate per-slice timing values across slices. "median" (default) or "mean".

• timing (str) – Which spike time to extract per unit per slice. "median" — median spike time within the slice (default). "mean" — mean spike time within the slice. "first" — first spike time (onset latency). Ignored when timing_matrix is provided.

• min_spikes (int) – Minimum number of spikes for a unit to count as active in a slice (default: 2). Ignored when timing_matrix is provided.

• min_frac_active (float or None) – Minimum fraction of slices a unit must be active in to be placed in the highly-active group. 0.0 or None (default: 0.0) skips the split entirely and places all units in the highly-active group without computing activity fractions.

• frac_active (np.ndarray or None) – Optional pre-computed fraction-active array of shape (U,) to override the internal calculation for the group split. Only used when min_frac_active > 0. Compatible sources: SpikeSliceStack.compute_frac_active and SpikeData.get_frac_active (frac_per_unit output).

• timing_matrix (np.ndarray or None) – Optional pre-computed (U, S) timing matrix from get_unit_timing_per_slice. When provided, timing and min_spikes are ignored and this matrix is used directly.

Returns:

Two SpikeSliceStack objects

(highly_active, low_active) with units reordered by typical timing. The low-activity stack is None when the group is empty.

unit_ids_in_order (tuple): Two ndarray

(highly_active, low_active) of original unit indices in the reordered sequence.

unit_std (tuple): Two ndarray (highly_active, low_active)

of standard deviation of per-slice timing values. Lower values indicate more consistent timing across slices.

unit_peak_times_ms (tuple): Two ndarray

(highly_active, low_active) of the aggregated typical timing in milliseconds relative to slice start. NaN for units with no active slices.

unit_frac_active (tuple): Two ndarray

(highly_active, low_active) of the fraction of slices each unit was active in.

Return type:

reordered_stacks (tuple)

Notes

• Call get_unit_timing_per_slice first to pre-compute the timing matrix if you want to reuse it across multiple calls (e.g. rank_order_correlation and this method).

• When frac_active is None and min_frac_active > 0, activity fraction is computed via compute_frac_active.

• Analogous to RateSliceStack.order_units_across_slices but operates on raw spike trains instead of firing rate curves.

apply(func, *args, **kwargs)

Apply a function to each SpikeData in the stack and return stacked results.

Calls func(sd, *args, **kwargs) on every slice and stacks the outputs into a single numpy array with a new leading axis of size S (number of slices).

Parameters:

• func (callable) – Function that accepts a SpikeData as its first argument and returns a numeric value (scalar, 1-D, or 2-D array). Output shape must be consistent across all slices.

• *args – Additional positional arguments forwarded to func.

• **kwargs – Additional keyword arguments forwarded to func.

Returns:

Stacked results with shape (S, ...).

Return type:

result (np.ndarray)

Notes

• Intended for use with stacks built by SpikeData.frames, SpikeData.align_to_events, SpikeData.spike_shuffle_stack, or SpikeData.subset_stack. Pair with shuffle_z_score, shuffle_percentile, slice_trend, or slice_stability from utils to interpret the results.

unit_to_unit_comparison(metric='ccg', delt=20.0, bin_size=1.0, max_lag=350, n_jobs=-1)

Compute pairwise unit-to-unit similarity within each slice using spike-based metrics.

For each slice, computes a (U, U) similarity matrix between all unit pairs, then stacks the results into a PairwiseCompMatrixStack (U, U, S).

Parameters:

• metric (str) – Similarity metric to use. "ccg" for cross-correlogram on binned rasters (default), "sttc" for spike time tiling coefficient.

• delt (float) – STTC time window in milliseconds (default: 20.0). Only used when metric is "sttc".

• bin_size (float) – Bin size in milliseconds for the binary raster (default: 1.0). Only used when metric is "ccg".

• max_lag (float) – Maximum lag in milliseconds to search for the peak correlation (default: 350). Only used when metric is "ccg".

Returns:

Pairwise similarity scores between

all unit pairs for each slice. Shape is (U, U, S).

lag_stack (PairwiseCompMatrixStack or None): Lag at which maximum

similarity occurs for each pair per slice. Shape is (U, U, S). None when metric is "sttc" (STTC has no lag).

av_corr (np.ndarray): Average similarity per slice across all unit

pairs in the lower triangle. Shape is (S,).

av_lag (np.ndarray or None): Average lag per slice. Shape is (S,).

None when metric is "sttc".

Return type:

corr_stack (PairwiseCompMatrixStack)

Notes

• Analogous to RateSliceStack.unit_to_unit_correlation but operates on raw spike trains instead of firing rate time series.

get_slice_to_slice_unit_comparison(metric='ccg', delt=20.0, bin_size=1.0, max_lag=350, min_spikes=2, min_frac=0.3, frac_active=None, n_jobs=-1)

Compute slice-to-slice similarity for each unit using spike-based metrics.

For each unit independently, compares its spike train across every pair of slices. Asks: “Does unit X fire in the same temporal pattern across repeated events?” Returns a PairwiseCompMatrixStack (S, S, U).

Parameters:

• metric (str) – Similarity metric to use. "ccg" for cross-correlogram on binned rasters (default), "sttc" for spike time tiling coefficient.

• delt (float) – STTC time window in milliseconds (default: 20.0). Only used when metric is "sttc".

• bin_size (float) – Bin size in milliseconds for the binary raster (default: 1.0). Only used when metric is "ccg".

• max_lag (float) – Maximum lag in milliseconds to search for the peak correlation (default: 350). Only used when metric is "ccg".

• min_spikes (int) – Minimum number of spikes in a slice for a unit to be considered valid in that slice (default: 2).

• min_frac (float) – Maximum fraction of slices that can be invalid before a unit’s average is set to NaN (default: 0.3).

• frac_active (np.ndarray or None) – Optional pre-computed fraction-active array of shape (U,) to override the internal per-unit validity check for computing averages. When provided, a unit’s average is set to NaN if frac_active[u] < (1 - min_frac). min_spikes still controls which individual slice pairs are computed. Compatible sources: SpikeSliceStack.compute_frac_active and SpikeData.get_frac_active (frac_per_unit output).

• n_jobs (int) – Number of threads for parallel computation. -1 uses all cores (default), 1 disables parallelism, None is serial.

Returns:

Pairwise similarity between all

slice pairs for each unit. Shape is (S, S, U).

all_lag (PairwiseCompMatrixStack or None): Lag at which maximum

similarity occurs for each slice pair per unit. Shape is (S, S, U). None when metric is "sttc".

av_corr (np.ndarray): Average similarity per unit across all valid

slice pairs. Shape is (U,).

av_lag (np.ndarray or None): Average lag per unit. Shape is (U,).

None when metric is "sttc".

Return type:

all_corr (PairwiseCompMatrixStack)

Notes

• Analogous to RateSliceStack.get_slice_to_slice_unit_corr_from_stack but operates on raw spike trains.

• Spike times within each slice are relative to the slice time origin (0-based or event-centered) for aligned comparison.

get_unit_timing_per_slice(timing='median', min_spikes=2)

Compute a representative spike time for each unit in each slice.

Returns a (U, S) matrix where entry [u, s] is the timing value (in milliseconds relative to the slice’s time origin) for unit u in slice s. For event-centered slices, t=0 is the event. Units with fewer than min_spikes spikes in a slice are marked NaN.

Parameters:

• timing (str) – Which spike time to extract per unit per slice. "median" (default), "mean", or "first" (onset latency).

• min_spikes (int) – Minimum number of spikes for a unit to count as active in a slice (default: 2).

Returns:

Array of shape (U, S) with timing

values in milliseconds relative to each slice’s time origin. NaN where the unit is inactive.

Return type:

timing_matrix (np.ndarray)

Notes

• Values are in milliseconds, not bin indices. This differs from RateSliceStack.get_unit_timing_per_slice which returns bin indices (suitable for direct indexing into the event stack). Both representations preserve rank order, so rank_order_correlation produces identical results either way.

• The returned matrix can be passed to rank_order_correlation to compute Spearman rank correlations between slice pairs, or used as input to order_units_across_slices for manual inspection of per-slice timing values.

rank_order_correlation(timing_matrix=None, timing='median', min_spikes=2, min_overlap=3, n_shuffles=100, min_overlap_frac=None, seed=1, n_jobs=-1)

Compute Spearman rank-order correlation of unit timing between all slice pairs.

For each pair of slices, only units active in both slices (non-NaN in both columns of the timing matrix) are included. If the overlap falls below the required minimum, the pair is set to NaN.

When n_shuffles > 0, the rank orders are shuffled n_shuffles times for each pair to build a null distribution, and the raw correlation is z-score normalised against it.

Parameters:

• timing_matrix (np.ndarray or None) – Array of shape (U, S) with timing values per unit per slice. NaN entries mark inactive units. Typically produced by get_unit_timing_per_slice. When None, computed automatically using timing and min_spikes.

• timing (str) – Which spike time to extract per unit per slice. "median" (default), "mean", or "first". Only used when timing_matrix is None.

• min_spikes (int) – Minimum spikes for activity (default: 2). Only used when timing_matrix is None.

• min_overlap (int) – Minimum number of units that must be active in both slices (default: 3).

• min_overlap_frac (float or None) – Minimum fraction of total units that must be active in both slices (default: None). When provided, the effective threshold is max(min_overlap, ceil(min_overlap_frac * U)).

• n_shuffles (int) – Number of shuffle iterations for z-scoring (default: 100). Set to 0 to return raw Spearman correlations. Values between 1 and 4 are rejected (minimum 5 required for a meaningful null distribution).

• seed (int or None) – Random seed for reproducibility of the shuffle (default: 1).

• n_jobs (int) – Number of threads for parallel computation. -1 uses all cores (default), 1 disables parallelism, None is serial.

Returns:

Spearman correlation matrix of

shape (S, S). When n_shuffles > 0, values are z-scores. When n_shuffles == 0, values are raw Spearman correlations.

av_corr (float): Average correlation (or z-score) across all valid

lower-triangle pairs.

overlap_matrix (PairwiseCompMatrix): Matrix of shape (S, S)

with fraction of units active in both slices.

Return type:

corr_matrix (PairwiseCompMatrix)

plot_aligned_slice_single_unit(unit_idx, ax=None, color_vals=None, color_label='', cmap='viridis', time_offset=0, xlabel='Rel. time (ms)', ylabel='Burst', x_range=None, vlines=None, show_colorbar=True, marker_size=20, font_size=None, style='scatter', invert_y=False, linewidths=0.5)

Plot a single unit’s spike times across all slices as a raster.

Extracts the spike train for unit_idx from every slice and delegates to plot_aligned_slice_single_unit().

Parameters:

• unit_idx (int) – Index of the unit to plot.

• ax (matplotlib.axes.Axes or None) – Target axes. If None, a new figure and axes are created.

• color_vals (np.ndarray or None) – Per-slice colour values.

• color_label (str) – Colorbar label.

• cmap (str) – Matplotlib colormap name.

• time_offset (float) – Value subtracted from every spike time before plotting. Slices from align_to_events are already event-centered (spike times in [-pre_ms, +post_ms]), so use the default time_offset=0. Only set a non-zero value when spike times are not already centered on the event.

• xlabel (str) – X-axis label.

• ylabel (str) – Y-axis label.

• x_range (tuple or None) – (xmin, xmax) for the x-axis.

• vlines (list[dict] or None) – Vertical reference lines. Each dict must contain 'x' and may optionally include 'color', 'linestyle', 'linewidth'.

• show_colorbar (bool) – Add a colorbar when color_vals is provided.

• marker_size (float) – Scatter marker size.

• font_size (int or None) – Font size for labels/ticks.

• style (str) – "scatter" for dot markers, "eventplot" for vertical line markers.

• invert_y (bool) – If True, first slice at top, last at bottom.

• linewidths (float) – Line width for eventplot markers.

Returns:

(fig, ax, sc) when ax is None, otherwise just sc.

sc is the scatter PathCollection (or None if no colour coding).

Return type:

result