RateSliceStack

The RateSliceStack class stores event-aligned slices of firing rate data. It holds a 3-D array of shape (num_units, num_bins, num_slices) where each slice corresponds to one trial or event window. This structure supports trial-averaging, per-unit statistics, and visualization of peri-event firing rate profiles.

class spikelab.spikedata.rateslicestack.RateSliceStack(data_obj=None, times_start_to_end=None, time_peaks=None, time_bounds=None, sigma_ms=10, event_matrix=None, step_size=None, neuron_attributes=None)

Bases: object

A 3D firing rate matrix of shape (U, T, S) with correlation and similarity capabilities.

U is units (neurons), T is time bins, and S is slices (bursts, events, etc). Construct from either a data_obj (SpikeData or RateData) with time specifications, or directly from a pre-built event_matrix. The instance variables are the same regardless of input option.

Parameters:

• data_obj (SpikeData or RateData) – A data object to slice. Provide either this or event_matrix, 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.

• sigma_ms (float) – Smoothing factor for computing ISI if you input a SpikeData object. Otherwise not used.

• event_matrix (np.ndarray or None) – A 3D array of shape (U, T, S). Provide either this or data_obj, not both.

• step_size (float or None) – Time resolution in milliseconds between consecutive time bins. If None, defaults to 1.0.

• neuron_attributes (list or None) – List of attribute objects, one per unit, containing arbitrary metadata about each neuron.

event_stack

3D array of shape (U, T, S) where U is the number of units, T is the number of time bins, and S is the number of slices.

Type:

np.ndarray

times

List of (start, end) time bounds for each slice, sorted chronologically. Length equals S. Example: [(100, 200), (500, 600), (1000, 1100)].

Type:

list

step_size

Time resolution in milliseconds between consecutive time bins. Inferred from input data. For SpikeData input, defaults to 1.0 ms. Example: 1.0 means time bins are at [100, 101, 102, …] ms.

Type:

float

neuron_attributes

List of attribute objects, one per unit, containing arbitrary metadata about each neuron.

Type:

list or None

__init__(data_obj=None, times_start_to_end=None, time_peaks=None, time_bounds=None, sigma_ms=10, event_matrix=None, step_size=None, neuron_attributes=None)

order_units_across_slices(agg_func, MIN_RATE_THRESHOLD=0.1, MIN_FRAC_ACTIVE=0.0, frac_active=None, timing_matrix=None)

Reorder units from earliest to latest peak firing rate across slices.

Parameters:

• agg_func (str) – Either "median" or "mean". Used for calculating the median/mean time when each unit has peak firing rate.

• MIN_RATE_THRESHOLD (float) – Minimum peak firing rate for a slice to be included in the ordering calculation. Slices where a unit’s max rate is below this threshold are excluded from that unit’s typical peak time calculation. Ignored when timing_matrix is provided.

• MIN_FRAC_ACTIVE (float) – Minimum fraction of slices a unit must be active in to be placed in the highly-active group. Default 0.0 means all units are in the first group, so the second array in each output tuple will be empty.

• frac_active (np.ndarray or None) – Optional pre-computed fraction-active array of shape (U,) to use for the group split instead of the rate-based calculation. 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, MIN_RATE_THRESHOLD is ignored and this matrix is used directly.

Returns:

Tuple of two 3D

arrays from event_stack with the U dimension reordered temporally. The first array is the highly-active group and the second is the lower-activity group.

unit_ids_in_order (tuple of arrays): Two arrays of original unit

IDs in temporal order (highly-active, low-activity). For example, [3, 1, 0, 2] means unit 3 fires first. Use this to map back to original unit IDs.

unit_std_indices (tuple of arrays): Two arrays of standard

deviation of peak firing rate times (highly-active, low-activity). Lower values indicate more consistent timing across slices.

unit_peak_times (tuple of arrays): Two arrays of median/mean peak

firing time bin (highly-active, low-activity).

unit_frac_active (tuple of arrays): Two arrays of the fraction of

slices each unit was active in (highly-active, low-activity).

Return type:

reordered_slice_matrices (tuple of arrays)

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).

get_slice_to_slice_unit_corr_from_stack(compare_func=<function compute_cross_correlation_with_lag>, MIN_RATE_THRESHOLD=0.1, MIN_FRAC=0.3, max_lag=10, frac_active=None, n_jobs=-1)

Compute slice-to-slice similarity along the unit axis of event_stack (U, T, S).

Output is a PairwiseCompMatrixStack of shape (S, S, U).

Parameters:

• compare_func (callable) – Comparison function from utils. Specify cross-correlation or cosine similarity. The default is cross correlation. See utils.py for details.

• MIN_RATE_THRESHOLD (float) – Minimum mean firing rate to consider a slice valid for that neuron.

• MIN_FRAC (float) – Maximum fraction of slices that can be skipped before a unit is deemed invalid (default 0.3 = 30%).

• max_lag (int) – Maximum lag in frames to search for similarity. If None, lag is set to 0.

• frac_active (np.ndarray or None) – Optional pre-computed fraction-active array of shape (U,) to override the internal rate-based validity check for computing averages. When provided, a unit’s average is set to NaN if frac_active[u] < (1 - MIN_FRAC). MIN_RATE_THRESHOLD 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.

Notes

When max_lag is 0 or None, the inner S x S loop is replaced by a vectorized matrix multiplication, which is significantly faster for large S. For non-zero max_lag, the S x S comparisons are computed in a serial loop per unit (parallelised across units). This can be slow for large S (e.g. S > 100).

Returns:

Pairwise

correlation scores between all slice pairs for each unit. Shape is (S, S, U) in the stack attribute.

av_slice_corr_scores (np.ndarray): Average correlation per neuron

across all valid slice pairs. Shape is (U,).

Return type:

all_slice_corr_scores (PairwiseCompMatrixStack)

get_slice_to_slice_time_corr_from_stack(compare_func=<function compute_cosine_similarity_with_lag>, max_lag=0, n_jobs=-1)

Compute slice-to-slice similarity along the time axis of event_stack (U, T, S).

Output is a PairwiseCompMatrixStack of shape (S, S, T).

Parameters:

• compare_func (callable) – Comparison function from utils. Specify cross-correlation or cosine similarity. The default is cosine similarity. See utils.py for details.

• max_lag (int) – Maximum lag in frames to search for similarity. If None, lag is set to 0.

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

Returns:

Pairwise

correlation scores between all slice pairs for each time bin. Shape is (S, S, T) in the stack attribute.

av_slice_corr_scores (np.ndarray): Average correlation per time

bin across all valid slice pairs. Shape is (T,).

Return type:

all_slice_corr_scores (PairwiseCompMatrixStack)

convert_to_list_of_RateData()

Create a list of RateData objects from the 3D event_stack.

Returns:

List of RateData objects. Length equals S.

Return type:

output (list)

unit_to_unit_correlation(compare_func=<function compute_cross_correlation_with_lag>, max_lag=10, n_jobs=-1)

Compute unit-to-unit similarity along the slice axis of event_stack (U, T, S).

Output is a PairwiseCompMatrixStack of shape (U, U, S).

Parameters:

• compare_func (callable) – Comparison function from utils. Specify cross-correlation or cosine similarity. The default is cross correlation. See utils.py for details.

• max_lag (int) – Maximum lag in frames to search for similarity. If None, lag is set to 0.

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

Returns:

Pairwise correlation

scores between all unit pairs for each slice. Shape is (U, U, S) in the stack attribute.

max_corr_lag_array (PairwiseCompMatrixStack): Lag where

correlation between each pair is at its maximum. Shape is (U, U, S) in the stack attribute.

av_max_corr (np.ndarray): Average correlation per slice across

all valid unit pairs. Shape is (S,).

av_max_corr_lag (np.ndarray): Average lag where correlation

between each pair is at its maximum. Shape is (S,).

Return type:

max_corr_array (PairwiseCompMatrixStack)

subset(units, by=None)

Extract a subset of units/neurons from the rate slice stack.

Parameters:

• units (list or array) – Unit indices to extract. If by is None, this should always be a list of ints. If by is not None, the list can contain ints or strings.

• by (str or None) – Neuron attribute key to match against. Only use this if you initialized the object with neuron_attributes. Set to the key that contains neuron_id values. None selects by index (default).

Returns:

New RateSliceStack object containing

only the specified units.

Return type:

result (RateSliceStack)

subtime_by_index(start_idx, end_idx)

Extract a subset of time bins from every slice using index values.

Trims along the time axis (T dimension) while preserving all slices (S dimension).

Parameters:

• start_idx (int) – Starting time bin index (inclusive). Supports negative indexing.

• end_idx (int) – Ending time bin index (exclusive). Supports negative indexing.

Returns:

New RateSliceStack where each slice

contains only the specified time bins. Shape changes from (U, T, S) to (U, T_trimmed, S).

Return type:

result (RateSliceStack)

Notes

• Original timestamps are preserved (not shifted to zero). To get shifted-to-zero timestamps, create a new RateSliceStack.

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

subslice(slices)

Extract a subset of slices from the event stack using index values.

Trims along the slice axis (S dimension) while preserving all time bins (T dimension).

Parameters:

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

Returns:

New RateSliceStack containing only the

specified slices. Shape changes from (U, T, S) to (U, T, S_trimmed).

Return type:

result (RateSliceStack)

Notes

• All units, neuron_attributes, and step_size are carried over from the original.

get_unit_timing_per_slice(MIN_RATE_THRESHOLD=0.1)

Compute the peak firing rate time bin for each unit in each slice.

Returns a (U, S) matrix where entry [u, s] is the time bin index of the peak firing rate for unit u in slice s. Units whose peak rate falls below MIN_RATE_THRESHOLD in a slice are marked NaN.

Parameters:

MIN_RATE_THRESHOLD (float) – Minimum peak firing rate for a unit to count as active in a slice (default: 0.1).

Returns:

Array of shape (U, S) with peak

time bin indices (integers cast to float for NaN support). These can be used to index directly into the (U, T, S) event stack. NaN where the unit is inactive.

Return type:

timing_matrix (np.ndarray)

Notes

• Values are bin indices, not milliseconds. This differs from SpikeSliceStack.get_unit_timing_per_slice which returns milliseconds. 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.

rank_order_correlation(timing_matrix=None, MIN_RATE_THRESHOLD=0.1, min_overlap=3, min_overlap_frac=None, n_shuffles=100, 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 MIN_RATE_THRESHOLD.

• MIN_RATE_THRESHOLD (float) – Minimum peak firing rate threshold (default: 0.1). 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)