From 6898e21e98aac3ad37cf21ebf7cde3ecd74c943d Mon Sep 17 00:00:00 2001 From: Kamran Hussain Date: Sun, 9 Nov 2025 13:46:03 -0800 Subject: [PATCH 1/5] add neuron to channel raster and frame level raster (mask of spikes over frames for each channel --- spikedata/spikedata.py | 260 ++++++++++++++++++++++++++++++++++++----- test_spikedata.py | 166 ++++++++++++++++++++++++++ 2 files changed, 398 insertions(+), 28 deletions(-) diff --git a/spikedata/spikedata.py b/spikedata/spikedata.py index 511bb94..cf4795a 100644 --- a/spikedata/spikedata.py +++ b/spikedata/spikedata.py @@ -5,7 +5,7 @@ import warnings from collections import namedtuple from dataclasses import dataclass -from typing import Literal +from typing import Literal, Union import numpy as np from numpy.typing import NDArray @@ -191,7 +191,7 @@ def from_thresholding( data: NDArray, fs_Hz=20e3, threshold_sigma=5.0, - filter: dict | bool = True, + filter: Union[dict, bool] = True, hysteresis=True, direction: Literal["both", "up", "down"] = "both", ): @@ -234,7 +234,7 @@ def __init__( neuron_attributes=None, metadata={}, raw_data=None, - raw_time: NDArray | float | None = None, + raw_time: Union[NDArray, float, None] = None, ): """ Initialize a SpikeData object using a list of spike trains, each a @@ -524,6 +524,211 @@ def raster(self, bin_size=20.0): """ return self.sparse_raster(bin_size).toarray() + def frame_raster(self, frame_rate_hz=30.0): + """ + Bin all spike times at a specific frame rate to create a frame-level raster. + + This is useful for aligning spike data with video frame rates or other + frame-based data sources. The bin size is automatically calculated from the + frame rate (bin_size = 1000 / frame_rate_hz ms). + + Parameters + ---------- + frame_rate_hz : float, default=30.0 + Frame rate in Hz (frames per second). The bin size will be 1000/frame_rate_hz ms. + + Returns + ------- + numpy.ndarray + Dense array of shape (N, T) where N is the number of neurons and T is the + number of frames. Entry (i, j) is the number of times neuron i fired in frame j. + """ + bin_size_ms = 1000.0 / frame_rate_hz + return self.raster(bin_size_ms) + + def sparse_frame_raster(self, frame_rate_hz=30.0): + """ + Bin all spike times at a specific frame rate to create a sparse frame-level raster. + + This is useful for aligning spike data with video frame rates or other + frame-based data sources when memory efficiency is important. + + Parameters + ---------- + frame_rate_hz : float, default=30.0 + Frame rate in Hz (frames per second). The bin size will be 1000/frame_rate_hz ms. + + Returns + ------- + scipy.sparse.csr_array + Sparse array of shape (N, T) where N is the number of neurons and T is the + number of frames. Entry (i, j) is the number of times neuron i fired in frame j. + """ + bin_size_ms = 1000.0 / frame_rate_hz + return self.sparse_raster(bin_size_ms) + + def channel_raster( + self, channel_map, bin_size=20.0, sparse_output=False, binary=False + ): + """ + Create a raster organized by channels instead of individual neurons. + + Aggregates spikes from multiple neurons that belong to the same channel into + a single channel-level representation. This is useful when neurons are recorded + from the same physical channel or when you want to analyze activity at the + channel level rather than the neuron level. + + Parameters + ---------- + channel_map : array-like, dict, or str + Mapping from neuron indices to channel indices. Can be: + - Array-like of length N: channel_map[i] is the channel index for neuron i + - Dict: channel_map[i] is the channel index for neuron i + - String: name of attribute in neuron_attributes (e.g., 'channel_id') + bin_size : float, default=20.0 + Bin size in milliseconds for the raster. + sparse_output : bool, default=False + If True, return a sparse array. If False, return a dense array. + binary : bool, default=False + If True, return a binary mask (0/1) indicating presence of any spike. + If False, return spike counts per channel per bin. + + Returns + ------- + numpy.ndarray or scipy.sparse.csr_array + Array of shape (C, T) where C is the number of channels and T is the + number of time bins. Entry (c, t) indicates spike activity for channel c + in time bin t (either count or binary mask depending on `binary` parameter). + + Examples + -------- + >>> # Map neurons 0,1,2 to channel 0, neuron 3 to channel 1 + >>> channel_map = [0, 0, 0, 1] + >>> channel_raster = sd.channel_raster(channel_map, bin_size=10) + """ + # Resolve channel mapping + if isinstance(channel_map, str): + if self.neuron_attributes is None: + raise ValueError( + f"Cannot use channel_map='{channel_map}' without neuron_attributes" + ) + channel_map = [ + getattr(attr, channel_map) for attr in self.neuron_attributes + ] + elif isinstance(channel_map, dict): + channel_map = [channel_map.get(i, -1) for i in range(self.N)] + + channel_map = np.asarray(channel_map) + if len(channel_map) != self.N: + raise ValueError( + f"channel_map length ({len(channel_map)}) must match number of neurons ({self.N})" + ) + + # Get unique channels and create mapping + unique_channels = np.unique(channel_map) + unique_channels = unique_channels[unique_channels >= 0] # Remove -1 if present + n_channels = len(unique_channels) + channel_to_idx = {ch: idx for idx, ch in enumerate(unique_channels)} + + # Get neuron-level raster + neuron_raster = self.sparse_raster(bin_size) + + # Aggregate by channel + if sparse.issparse(neuron_raster): + # Convert to COO format for easier manipulation + neuron_raster_coo = neuron_raster.tocoo() + # Map neuron indices to channel indices + channel_indices = np.array( + [ + channel_to_idx.get(channel_map[neuron_idx], -1) + for neuron_idx in neuron_raster_coo.row + ] + ) + # Filter out invalid channels + valid_mask = channel_indices >= 0 + channel_indices = channel_indices[valid_mask] + time_indices = neuron_raster_coo.col[valid_mask] + values = neuron_raster_coo.data[valid_mask] + + # Aggregate spikes for the same (channel, time) pair + n_bins = neuron_raster.shape[1] + if binary: + # Binary mask: just mark presence + # Use a temporary dense array to aggregate, then convert to sparse + temp_dense = np.zeros((n_channels, n_bins), dtype=int) + for ch_idx, t_idx in zip(channel_indices, time_indices): + temp_dense[ch_idx, t_idx] = 1 + channel_raster = sparse.csr_array(temp_dense) + else: + # Count spikes - need to sum values for same (channel, time) pairs + # Use a temporary dense array to aggregate properly + temp_dense = np.zeros((n_channels, n_bins), dtype=float) + for ch_idx, t_idx, val in zip(channel_indices, time_indices, values): + temp_dense[ch_idx, t_idx] += val + channel_raster = sparse.csr_array(temp_dense) + else: + # Dense case + n_bins = neuron_raster.shape[1] + channel_raster = np.zeros((n_channels, n_bins), dtype=int) + for neuron_idx in range(self.N): + channel_idx = channel_to_idx.get(channel_map[neuron_idx], -1) + if channel_idx >= 0: + if binary: + channel_raster[channel_idx] = ( + channel_raster[channel_idx] + | (neuron_raster[neuron_idx] > 0) + ).astype(int) + else: + channel_raster[channel_idx] += neuron_raster[neuron_idx] + + if not sparse_output and sparse.issparse(channel_raster): + return channel_raster.toarray() + elif sparse_output and not sparse.issparse(channel_raster): + return sparse.csr_array(channel_raster) + else: + return channel_raster + + def sparse_channel_raster(self, channel_map, bin_size=20.0, binary=False): + """ + Create a sparse raster organized by channels instead of individual neurons. + + This is a convenience method that calls channel_raster with sparse_output=True. + See channel_raster() for detailed documentation. + """ + return self.channel_raster( + channel_map, bin_size, sparse_output=True, binary=binary + ) + + def frame_channel_raster( + self, channel_map, frame_rate_hz=30.0, sparse_output=False, binary=False + ): + """ + Create a frame-level raster organized by channels. + + Combines frame-level binning with channel aggregation. Useful for aligning + channel-level spike activity with video frames or other frame-based data. + + Parameters + ---------- + channel_map : array-like, dict, or str + Mapping from neuron indices to channel indices. See channel_raster() for details. + frame_rate_hz : float, default=30.0 + Frame rate in Hz (frames per second). + sparse_output : bool, default=False + If True, return a sparse array. + binary : bool, default=False + If True, return a binary mask (0/1) indicating presence of any spike. + + Returns + ------- + numpy.ndarray or scipy.sparse.csr_array + Array of shape (C, F) where C is the number of channels and F is the + number of frames. Entry (c, f) indicates spike activity for channel c + in frame f. + """ + bin_size_ms = 1000.0 / frame_rate_hz + return self.channel_raster(channel_map, bin_size_ms, sparse_output, binary) + def interspike_intervals(self): "Produce a list of arrays of interspike intervals per unit." return [np.diff(ts) for ts in self.train] @@ -799,7 +1004,7 @@ def population_firing_rate(self, bin_size=10, w=5, average=False): def population_firing_rate( - trains: list[NDArray] | NDArray, + trains: Union[list[NDArray], NDArray], rec_length=None, bin_size=10, w=5, @@ -836,7 +1041,7 @@ def population_firing_rate( return bins, fr_pop -def spike_time_tiling(tA, tB, delt=20.0, length: float | None = None): +def spike_time_tiling(tA, tB, delt=20.0, length: Union[float, None] = None): """ Calculate the spike time tiling coefficient [1] between two spike trains. STTC is a metric for correlation between spike trains with some improved intuitive properties @@ -870,7 +1075,7 @@ def _spike_time_tiling(tA, tB, TA, TB, delt): return (aa + bb) / 2 -def best_effort_sample(counts, M, rng: np.random.Generator | None = None): +def best_effort_sample(counts, M, rng: Union[np.random.Generator, None] = None): """ Given a discrete distribution over the integers 0...N-1 in the form of an array of N counts, sample M elements from the distribution without replacement if possible. If @@ -893,7 +1098,7 @@ def best_effort_sample(counts, M, rng: np.random.Generator | None = None): return ret -def randomize_raster(raster, seed: int | None = None, method="poprate_greedy"): +def randomize_raster(raster, seed: Union[int, None] = None, method="poprate_greedy"): """ Generate a randomized version of a spike raster by reassigning all its spike times to new bins, while preserving the total number of spikes for each unit. @@ -926,7 +1131,7 @@ def _okun_swap(ar, idxs, rng): return True -def randomize_raster_okun(raster, seed: int | None = None, swap_per_spike=5): +def randomize_raster_okun(raster, seed: Union[int, None] = None, swap_per_spike=5): """ Generate a randomized version of a spike raster, preserving population rate. The input raster MUST have at most one spike per neuron per bin! @@ -958,7 +1163,7 @@ def randomize_raster_okun(raster, seed: int | None = None, swap_per_spike=5): return raster -def randomize_raster_greedy(raster, seed: int | None = None): +def randomize_raster_greedy(raster, seed: Union[int, None] = None): """ Generate a randomized version of a spike raster, preserving population rate. @@ -1212,8 +1417,8 @@ def burst_detection(spike_times, burst_threshold, spike_num_thr=3): def butter_filter( data, - lowcut: float | None = None, - highcut: float | None = None, + lowcut: Union[float, None] = None, + highcut: Union[float, None] = None, fs=20000.0, order=5, ): @@ -1234,23 +1439,22 @@ def butter_filter( Returns: The filtered output with the same shape as data """ - match (lowcut, highcut): - case None, None: - raise ValueError( - "Need at least a low cutoff (lowcut) or high cutoff (highcut) frequency!" - ) - case None, _: - filter_type = "lowpass" - Wn = highcut / fs * 2 - case _, None: - filter_type = "highpass" - Wn = lowcut / fs * 2 - case _, _: - if lowcut >= highcut: - raise ValueError("lowcut must be smaller than highcut") - filter_type = "bandpass" - band = [lowcut, highcut] - Wn = [e / fs * 2 for e in band] + if lowcut is None and highcut is None: + raise ValueError( + "Need at least a low cutoff (lowcut) or high cutoff (highcut) frequency!" + ) + elif lowcut is None: + filter_type = "lowpass" + Wn = highcut / fs * 2 + elif highcut is None: + filter_type = "highpass" + Wn = lowcut / fs * 2 + else: + if lowcut >= highcut: + raise ValueError("lowcut must be smaller than highcut") + filter_type = "bandpass" + band = [lowcut, highcut] + Wn = [e / fs * 2 for e in band] filter_coeff = signal.iirfilter( order, Wn, analog=False, btype=filter_type, output="sos" diff --git a/test_spikedata.py b/test_spikedata.py index dfa4506..5766794 100644 --- a/test_spikedata.py +++ b/test_spikedata.py @@ -640,3 +640,169 @@ def test_randomization_issue_13(self): rr = randomize_raster_greedy(r) self.assertAll(r.sum(0) == rr.sum(0)) self.assertAll(r.sum(1) == rr.sum(1)) + + def test_frame_raster(self): + # Test frame-level binning at different frame rates + # Create spike data with known spike times + sd = SpikeData([[0, 33.33, 66.67, 100]]) # Spikes at 0, 33.33, 66.67, 100 ms + + # At 30 fps (33.33 ms per frame), binning uses ceil(length/bin_size) bins + # length=100, bin_size=33.33, so ceil(100/33.33)=3 bins + # Spikes: 0->bin 0, 33.33->bin 0, 66.67->bin 2, 100->bin 2 + # (Note: ceil(66.67/33.33)-1 = ceil(2.0)-1 = 2-1 = 1, but 66.67 falls in bin 2 due to right-closed) + frame_raster_30 = sd.frame_raster(frame_rate_hz=30.0) + self.assertEqual(frame_raster_30.shape[0], 1) # 1 neuron + self.assertEqual(frame_raster_30.shape[1], 3) # 3 frames + self.assertAll( + frame_raster_30[0] == [2, 0, 2] + ) # 2 spikes in bin 0, 0 in bin 1, 2 in bin 2 + + # At 60 fps (16.67 ms per frame), should get more frames + frame_raster_60 = sd.frame_raster(frame_rate_hz=60.0) + self.assertGreater(frame_raster_60.shape[1], 4) + + # Test sparse version + sparse_frame_raster = sd.sparse_frame_raster(frame_rate_hz=30.0) + self.assertTrue(sparse.issparse(sparse_frame_raster)) + self.assertAll(sparse_frame_raster.toarray() == frame_raster_30) + + # Test with multiple neurons + sd_multi = SpikeData([[0, 50], [25, 75]], length=100) + frame_raster_multi = sd_multi.frame_raster( + frame_rate_hz=20.0 + ) # 50 ms per frame + self.assertEqual(frame_raster_multi.shape[0], 2) # 2 neurons + self.assertEqual(frame_raster_multi.shape[1], 2) # 2 frames (0-50, 50-100) + # Neuron 0: spikes at 0 and 50 both go to bin 0 (right-closed), neuron 1: spike at 25->bin 0, 75->bin 1 + self.assertAll(frame_raster_multi[0] == [2, 0]) # Neuron 0: 2 spikes in bin 0 + self.assertAll(frame_raster_multi[1] == [1, 1]) # Neuron 1: 1 spike in each bin + + def test_channel_raster(self): + # Create spike data with 4 neurons + # Neurons 0,1 map to channel 0; neurons 2,3 map to channel 1 + sd = SpikeData( + [[0, 20, 40], [10, 30], [5, 25], [15, 35]], length=50 + ) # 4 neurons + + # Test with array channel map + channel_map = [0, 0, 1, 1] + channel_raster = sd.channel_raster(channel_map, bin_size=10) + self.assertEqual(channel_raster.shape[0], 2) # 2 channels + self.assertEqual(channel_raster.shape[1], 5) # 5 bins (0-10, 10-20, ..., 40-50) + + # Channel 0 (neurons 0,1): spikes at 0,10,20,30 + # Channel 1 (neurons 2,3): spikes at 5,15,25,35 + # Bins are left-open, right-closed except first bin captures t=0 + # Bin 0: (0, 10] captures spikes at 0, 5, 10 -> Channel 0: 2 (0,10), Channel 1: 1 (5) + # Bin 1: (10, 20] captures spikes at 15, 20 -> Channel 0: 1 (20), Channel 1: 1 (15) + # Bin 2: (20, 30] captures spikes at 25, 30 -> Channel 0: 1 (30), Channel 1: 1 (25) + # Bin 3: (30, 40] captures spikes at 35, 40 -> Channel 0: 1 (40), Channel 1: 1 (35) + # Bin 4: (40, 50] captures no spikes -> Channel 0: 0, Channel 1: 0 + self.assertEqual( + channel_raster[0, 0], 2 + ) # Neuron 0 spike at 0, Neuron 1 spike at 10 + self.assertEqual(channel_raster[0, 1], 1) # Neuron 0 spike at 20 + self.assertEqual(channel_raster[0, 2], 1) # Neuron 1 spike at 30 + self.assertEqual(channel_raster[0, 3], 1) # Neuron 0 spike at 40 + self.assertEqual(channel_raster[1, 0], 1) # Neuron 2 spike at 5 + self.assertEqual(channel_raster[1, 1], 1) # Neuron 3 spike at 15 + self.assertEqual(channel_raster[1, 2], 1) # Neuron 2 spike at 25 + self.assertEqual(channel_raster[1, 3], 1) # Neuron 3 spike at 35 + + # Test binary mode + channel_raster_binary = sd.channel_raster(channel_map, bin_size=10, binary=True) + self.assertAll((channel_raster_binary == 0) | (channel_raster_binary == 1)) + # Should have same pattern but with max value 1 + self.assertAll(channel_raster_binary <= 1) + + # Test sparse version + sparse_channel_raster = sd.sparse_channel_raster(channel_map, bin_size=10) + self.assertTrue(sparse.issparse(sparse_channel_raster)) + self.assertAll(sparse_channel_raster.toarray() == channel_raster) + + # Test with dict channel map + channel_map_dict = {0: 0, 1: 0, 2: 1, 3: 1} + channel_raster_dict = sd.channel_raster(channel_map_dict, bin_size=10) + self.assertAll(channel_raster_dict == channel_raster) + + # Test with neuron_attributes + @dataclass + class ChannelAttributes: + channel_id: int + + attrs = [ + ChannelAttributes(0), + ChannelAttributes(0), + ChannelAttributes(1), + ChannelAttributes(1), + ] + sd_with_attrs = SpikeData( + [[0, 20, 40], [10, 30], [5, 25], [15, 35]], + length=50, + neuron_attributes=attrs, + ) + channel_raster_attr = sd_with_attrs.channel_raster("channel_id", bin_size=10) + self.assertAll(channel_raster_attr == channel_raster) + + # Test error cases + with self.assertRaises(ValueError): + sd.channel_raster([0, 1], bin_size=10) # Wrong length + + with self.assertRaises(ValueError): + sd.channel_raster("channel_id", bin_size=10) # No neuron_attributes + + def test_frame_channel_raster(self): + # Test combining frame-level binning with channel aggregation + sd = SpikeData( + [[0, 33.33, 66.67], [16.67, 50], [8.33, 41.67]], length=100 + ) # 3 neurons + + # Map neurons 0,1 to channel 0; neuron 2 to channel 1 + channel_map = [0, 0, 1] + + # At 30 fps, should get frame-level channel raster + # Bin size = 1000/30 = 33.33 ms, length=100, so ceil(100/33.33) = 3 bins + frame_channel_raster = sd.frame_channel_raster(channel_map, frame_rate_hz=30.0) + self.assertEqual(frame_channel_raster.shape[0], 2) # 2 channels + self.assertEqual(frame_channel_raster.shape[1], 3) # 3 frames at 30 fps + + # Channel 0 (neurons 0,1) should have spikes in multiple frames + # Channel 1 (neuron 2) should have spikes in multiple frames + self.assertGreater(frame_channel_raster[0].sum(), 0) + self.assertGreater(frame_channel_raster[1].sum(), 0) + + # Test binary mode + frame_channel_binary = sd.frame_channel_raster( + channel_map, frame_rate_hz=30.0, binary=True + ) + self.assertAll((frame_channel_binary == 0) | (frame_channel_binary == 1)) + + # Test sparse version + sparse_frame_channel = sd.frame_channel_raster( + channel_map, frame_rate_hz=30.0, sparse_output=True + ) + self.assertTrue(sparse.issparse(sparse_frame_channel)) + self.assertAll(sparse_frame_channel.toarray() == frame_channel_raster) + + def test_channel_raster_aggregation(self): + # Test that multiple neurons on the same channel properly aggregate + # Create data where multiple neurons spike in the same bin + sd = SpikeData([[0, 10], [5, 15], [8, 18]], length=20) + + # All neurons map to channel 0 + channel_map = [0, 0, 0] + channel_raster = sd.channel_raster(channel_map, bin_size=10) + + self.assertEqual(channel_raster.shape[0], 1) # 1 channel + self.assertEqual(channel_raster.shape[1], 2) # 2 bins + + # First bin (0, 10]: neurons 0,1,2 spike at 0,5,8,10 -> 4 spikes total (10 is included due to right-closed) + self.assertEqual(channel_raster[0, 0], 4) + + # Second bin (10, 20]: neurons 0,1,2 spike at 15,18 -> 2 spikes total + self.assertEqual(channel_raster[0, 1], 2) + + # Test binary mode - should still be 1 even with multiple spikes + channel_raster_binary = sd.channel_raster(channel_map, bin_size=10, binary=True) + self.assertEqual(channel_raster_binary[0, 0], 1) + self.assertEqual(channel_raster_binary[0, 1], 1) From ffd5b6f0df30ed644bd0b279154e9c637f840293 Mon Sep 17 00:00:00 2001 From: Kamran Hussain Date: Sun, 9 Nov 2025 13:51:15 -0800 Subject: [PATCH 2/5] add automatic channel mapping resolution --- spikedata/spikedata.py | 203 +++++++++++++++++++++++++++++++++++------ test_spikedata.py | 162 +++++++++++++++++++++++--------- 2 files changed, 296 insertions(+), 69 deletions(-) diff --git a/spikedata/spikedata.py b/spikedata/spikedata.py index cf4795a..a8e590d 100644 --- a/spikedata/spikedata.py +++ b/spikedata/spikedata.py @@ -568,7 +568,12 @@ def sparse_frame_raster(self, frame_rate_hz=30.0): return self.sparse_raster(bin_size_ms) def channel_raster( - self, channel_map, bin_size=20.0, sparse_output=False, binary=False + self, + bin_size=20.0, + sparse_output=False, + binary=False, + attribute_name=None, + from_raw_data=False, ): """ Create a raster organized by channels instead of individual neurons. @@ -578,13 +583,12 @@ def channel_raster( from the same physical channel or when you want to analyze activity at the channel level rather than the neuron level. + The channel mapping is automatically derived using get_channel_map(), which + attempts to find channel information from neuron_attributes, raw_data shape, + or metadata. + Parameters ---------- - channel_map : array-like, dict, or str - Mapping from neuron indices to channel indices. Can be: - - Array-like of length N: channel_map[i] is the channel index for neuron i - - Dict: channel_map[i] is the channel index for neuron i - - String: name of attribute in neuron_attributes (e.g., 'channel_id') bin_size : float, default=20.0 Bin size in milliseconds for the raster. sparse_output : bool, default=False @@ -592,6 +596,13 @@ def channel_raster( binary : bool, default=False If True, return a binary mask (0/1) indicating presence of any spike. If False, return spike counts per channel per bin. + attribute_name : str, optional + Name of the attribute in neuron_attributes that contains channel information. + If None, will try common names: 'channel_id', 'channel', 'ch', 'channel_idx'. + Passed to get_channel_map(). + from_raw_data : bool, default=False + If True and raw_data is available, derive mapping from raw_data shape. + Passed to get_channel_map(). Returns ------- @@ -600,23 +611,32 @@ def channel_raster( number of time bins. Entry (c, t) indicates spike activity for channel c in time bin t (either count or binary mask depending on `binary` parameter). + Raises + ------ + ValueError + If no channel mapping can be automatically determined from the available data. + Examples -------- - >>> # Map neurons 0,1,2 to channel 0, neuron 3 to channel 1 - >>> channel_map = [0, 0, 0, 1] - >>> channel_raster = sd.channel_raster(channel_map, bin_size=10) - """ - # Resolve channel mapping - if isinstance(channel_map, str): - if self.neuron_attributes is None: - raise ValueError( - f"Cannot use channel_map='{channel_map}' without neuron_attributes" - ) - channel_map = [ - getattr(attr, channel_map) for attr in self.neuron_attributes - ] - elif isinstance(channel_map, dict): - channel_map = [channel_map.get(i, -1) for i in range(self.N)] + >>> # Auto-detect channel mapping from neuron_attributes + >>> channel_raster = sd.channel_raster(bin_size=10) + >>> + >>> # Specify custom attribute name + >>> channel_raster = sd.channel_raster(bin_size=10, attribute_name='electrode_id') + >>> + >>> # Derive from raw_data shape + >>> channel_raster = sd.channel_raster(bin_size=10, from_raw_data=True) + """ + # Automatically derive channel mapping + channel_map = self.get_channel_map( + attribute_name=attribute_name, from_raw_data=from_raw_data + ) + if channel_map is None: + raise ValueError( + "Could not automatically derive channel mapping. " + "Ensure neuron_attributes contains channel information (e.g., 'channel_id'), " + "or provide raw_data with from_raw_data=True, or store 'channel_map' in metadata." + ) channel_map = np.asarray(channel_map) if len(channel_map) != self.N: @@ -688,7 +708,9 @@ def channel_raster( else: return channel_raster - def sparse_channel_raster(self, channel_map, bin_size=20.0, binary=False): + def sparse_channel_raster( + self, bin_size=20.0, binary=False, attribute_name=None, from_raw_data=False + ): """ Create a sparse raster organized by channels instead of individual neurons. @@ -696,11 +718,20 @@ def sparse_channel_raster(self, channel_map, bin_size=20.0, binary=False): See channel_raster() for detailed documentation. """ return self.channel_raster( - channel_map, bin_size, sparse_output=True, binary=binary + bin_size, + sparse_output=True, + binary=binary, + attribute_name=attribute_name, + from_raw_data=from_raw_data, ) def frame_channel_raster( - self, channel_map, frame_rate_hz=30.0, sparse_output=False, binary=False + self, + frame_rate_hz=30.0, + sparse_output=False, + binary=False, + attribute_name=None, + from_raw_data=False, ): """ Create a frame-level raster organized by channels. @@ -708,16 +739,22 @@ def frame_channel_raster( Combines frame-level binning with channel aggregation. Useful for aligning channel-level spike activity with video frames or other frame-based data. + The channel mapping is automatically derived using get_channel_map(). + Parameters ---------- - channel_map : array-like, dict, or str - Mapping from neuron indices to channel indices. See channel_raster() for details. frame_rate_hz : float, default=30.0 Frame rate in Hz (frames per second). sparse_output : bool, default=False If True, return a sparse array. binary : bool, default=False If True, return a binary mask (0/1) indicating presence of any spike. + attribute_name : str, optional + Name of the attribute in neuron_attributes that contains channel information. + Passed to get_channel_map(). + from_raw_data : bool, default=False + If True and raw_data is available, derive mapping from raw_data shape. + Passed to get_channel_map(). Returns ------- @@ -725,9 +762,121 @@ def frame_channel_raster( Array of shape (C, F) where C is the number of channels and F is the number of frames. Entry (c, f) indicates spike activity for channel c in frame f. + + Raises + ------ + ValueError + If no channel mapping can be automatically determined from the available data. """ bin_size_ms = 1000.0 / frame_rate_hz - return self.channel_raster(channel_map, bin_size_ms, sparse_output, binary) + return self.channel_raster( + bin_size_ms, + sparse_output, + binary, + attribute_name=attribute_name, + from_raw_data=from_raw_data, + ) + + def get_channel_map(self, attribute_name=None, from_raw_data=False): + """ + Derive or extract the channel-to-neuron mapping from available data. + + This method attempts to automatically determine the channel mapping using + multiple strategies: + 1. From neuron_attributes if an attribute name is provided or common names exist + 2. From raw_data shape if it has channel structure + 3. Returns None if no mapping can be determined + + Parameters + ---------- + attribute_name : str, optional + Name of the attribute in neuron_attributes that contains channel information. + If None, will try common names: 'channel_id', 'channel', 'ch', 'channel_idx'. + from_raw_data : bool, default=False + If True and raw_data is available, derive mapping from raw_data shape. + Assumes raw_data has shape (channels, time) and neurons correspond to channels. + + Returns + ------- + numpy.ndarray or None + Array of length N where entry i is the channel index for neuron i. + Returns None if no mapping can be determined. + + Examples + -------- + >>> # Auto-detect from neuron_attributes with common attribute names + >>> channel_map = sd.get_channel_map() + >>> + >>> # Specify custom attribute name + >>> channel_map = sd.get_channel_map(attribute_name='electrode_id') + >>> + >>> # Derive from raw_data shape + >>> channel_map = sd.get_channel_map(from_raw_data=True) + """ + # Strategy 1: Try neuron_attributes + if self.neuron_attributes is not None: + # Try provided attribute name or common names + candidates = [] + if attribute_name: + candidates.append(attribute_name) + else: + candidates = [ + "channel_id", + "channel", + "ch", + "channel_idx", + "electrode_id", + "electrode", + ] + + for attr_name in candidates: + try: + channel_map = [ + getattr(attr, attr_name) for attr in self.neuron_attributes + ] + # Check if all values are valid (not None, numeric) + if all( + v is not None + and isinstance(v, (int, float, np.integer, np.floating)) + for v in channel_map + ): + return np.asarray(channel_map) + except AttributeError: + continue + + # Strategy 2: Derive from raw_data shape + if from_raw_data and hasattr(self, "raw_data") and self.raw_data.size > 0: + raw_shape = self.raw_data.shape + # If raw_data has shape (channels, time) or (channels, ...), + # and number of channels matches number of neurons + if len(raw_shape) >= 2: + n_channels = raw_shape[0] + if n_channels == self.N: + # Each neuron corresponds to one channel + return np.arange(n_channels) + elif n_channels < self.N: + # Multiple neurons per channel - need to map them + # This is a heuristic: assume neurons are grouped sequentially + neurons_per_channel = self.N // n_channels + channel_map = [] + for ch in range(n_channels): + channel_map.extend([ch] * neurons_per_channel) + # Handle remainder + remainder = self.N % n_channels + if remainder > 0: + channel_map.extend([n_channels - 1] * remainder) + return np.asarray(channel_map) + + # Strategy 3: Check metadata + if hasattr(self, "metadata") and isinstance(self.metadata, dict): + if "channel_map" in self.metadata: + channel_map = self.metadata["channel_map"] + if isinstance(channel_map, (list, np.ndarray)): + channel_map = np.asarray(channel_map) + if len(channel_map) == self.N: + return channel_map + + return None def interspike_intervals(self): "Produce a list of arrays of interspike intervals per unit." diff --git a/test_spikedata.py b/test_spikedata.py index 5766794..4f729d4 100644 --- a/test_spikedata.py +++ b/test_spikedata.py @@ -680,13 +680,24 @@ def test_frame_raster(self): def test_channel_raster(self): # Create spike data with 4 neurons # Neurons 0,1 map to channel 0; neurons 2,3 map to channel 1 + @dataclass + class ChannelAttributes: + channel_id: int + + attrs = [ + ChannelAttributes(0), + ChannelAttributes(0), + ChannelAttributes(1), + ChannelAttributes(1), + ] sd = SpikeData( - [[0, 20, 40], [10, 30], [5, 25], [15, 35]], length=50 + [[0, 20, 40], [10, 30], [5, 25], [15, 35]], + length=50, + neuron_attributes=attrs, ) # 4 neurons - # Test with array channel map - channel_map = [0, 0, 1, 1] - channel_raster = sd.channel_raster(channel_map, bin_size=10) + # Test auto-detection from neuron_attributes + channel_raster = sd.channel_raster(bin_size=10) self.assertEqual(channel_raster.shape[0], 2) # 2 channels self.assertEqual(channel_raster.shape[1], 5) # 5 bins (0-10, 10-20, ..., 40-50) @@ -710,22 +721,29 @@ def test_channel_raster(self): self.assertEqual(channel_raster[1, 3], 1) # Neuron 3 spike at 35 # Test binary mode - channel_raster_binary = sd.channel_raster(channel_map, bin_size=10, binary=True) + channel_raster_binary = sd.channel_raster(bin_size=10, binary=True) self.assertAll((channel_raster_binary == 0) | (channel_raster_binary == 1)) # Should have same pattern but with max value 1 self.assertAll(channel_raster_binary <= 1) # Test sparse version - sparse_channel_raster = sd.sparse_channel_raster(channel_map, bin_size=10) + sparse_channel_raster = sd.sparse_channel_raster(bin_size=10) self.assertTrue(sparse.issparse(sparse_channel_raster)) self.assertAll(sparse_channel_raster.toarray() == channel_raster) - # Test with dict channel map - channel_map_dict = {0: 0, 1: 0, 2: 1, 3: 1} - channel_raster_dict = sd.channel_raster(channel_map_dict, bin_size=10) - self.assertAll(channel_raster_dict == channel_raster) + # Test with custom attribute name + channel_raster_attr = sd.channel_raster( + bin_size=10, attribute_name="channel_id" + ) + self.assertAll(channel_raster_attr == channel_raster) - # Test with neuron_attributes + # Test error case - no channel mapping available + sd_no_map = SpikeData([[0, 20, 40], [10, 30], [5, 25], [15, 35]], length=50) + with self.assertRaises(ValueError): + sd_no_map.channel_raster(bin_size=10) # No channel mapping available + + def test_frame_channel_raster(self): + # Test combining frame-level binning with channel aggregation @dataclass class ChannelAttributes: channel_id: int @@ -734,35 +752,16 @@ class ChannelAttributes: ChannelAttributes(0), ChannelAttributes(0), ChannelAttributes(1), - ChannelAttributes(1), ] - sd_with_attrs = SpikeData( - [[0, 20, 40], [10, 30], [5, 25], [15, 35]], - length=50, - neuron_attributes=attrs, - ) - channel_raster_attr = sd_with_attrs.channel_raster("channel_id", bin_size=10) - self.assertAll(channel_raster_attr == channel_raster) - - # Test error cases - with self.assertRaises(ValueError): - sd.channel_raster([0, 1], bin_size=10) # Wrong length - - with self.assertRaises(ValueError): - sd.channel_raster("channel_id", bin_size=10) # No neuron_attributes - - def test_frame_channel_raster(self): - # Test combining frame-level binning with channel aggregation sd = SpikeData( - [[0, 33.33, 66.67], [16.67, 50], [8.33, 41.67]], length=100 + [[0, 33.33, 66.67], [16.67, 50], [8.33, 41.67]], + length=100, + neuron_attributes=attrs, ) # 3 neurons - # Map neurons 0,1 to channel 0; neuron 2 to channel 1 - channel_map = [0, 0, 1] - # At 30 fps, should get frame-level channel raster # Bin size = 1000/30 = 33.33 ms, length=100, so ceil(100/33.33) = 3 bins - frame_channel_raster = sd.frame_channel_raster(channel_map, frame_rate_hz=30.0) + frame_channel_raster = sd.frame_channel_raster(frame_rate_hz=30.0) self.assertEqual(frame_channel_raster.shape[0], 2) # 2 channels self.assertEqual(frame_channel_raster.shape[1], 3) # 3 frames at 30 fps @@ -772,14 +771,12 @@ def test_frame_channel_raster(self): self.assertGreater(frame_channel_raster[1].sum(), 0) # Test binary mode - frame_channel_binary = sd.frame_channel_raster( - channel_map, frame_rate_hz=30.0, binary=True - ) + frame_channel_binary = sd.frame_channel_raster(frame_rate_hz=30.0, binary=True) self.assertAll((frame_channel_binary == 0) | (frame_channel_binary == 1)) # Test sparse version sparse_frame_channel = sd.frame_channel_raster( - channel_map, frame_rate_hz=30.0, sparse_output=True + frame_rate_hz=30.0, sparse_output=True ) self.assertTrue(sparse.issparse(sparse_frame_channel)) self.assertAll(sparse_frame_channel.toarray() == frame_channel_raster) @@ -787,11 +784,19 @@ def test_frame_channel_raster(self): def test_channel_raster_aggregation(self): # Test that multiple neurons on the same channel properly aggregate # Create data where multiple neurons spike in the same bin - sd = SpikeData([[0, 10], [5, 15], [8, 18]], length=20) + @dataclass + class ChannelAttributes: + channel_id: int + + attrs = [ + ChannelAttributes(0), + ChannelAttributes(0), + ChannelAttributes(0), + ] + sd = SpikeData([[0, 10], [5, 15], [8, 18]], length=20, neuron_attributes=attrs) # All neurons map to channel 0 - channel_map = [0, 0, 0] - channel_raster = sd.channel_raster(channel_map, bin_size=10) + channel_raster = sd.channel_raster(bin_size=10) self.assertEqual(channel_raster.shape[0], 1) # 1 channel self.assertEqual(channel_raster.shape[1], 2) # 2 bins @@ -803,6 +808,79 @@ def test_channel_raster_aggregation(self): self.assertEqual(channel_raster[0, 1], 2) # Test binary mode - should still be 1 even with multiple spikes - channel_raster_binary = sd.channel_raster(channel_map, bin_size=10, binary=True) + channel_raster_binary = sd.channel_raster(bin_size=10, binary=True) self.assertEqual(channel_raster_binary[0, 0], 1) self.assertEqual(channel_raster_binary[0, 1], 1) + + def test_get_channel_map(self): + # Test auto-detection from neuron_attributes + @dataclass + class ChannelAttributes: + channel_id: int + + attrs = [ + ChannelAttributes(0), + ChannelAttributes(0), + ChannelAttributes(1), + ChannelAttributes(1), + ] + sd = SpikeData( + [[0, 20], [10, 30], [5, 25], [15, 35]], + length=50, + neuron_attributes=attrs, + ) + + # Auto-detect using common attribute name + channel_map = sd.get_channel_map() + self.assertIsNotNone(channel_map) + self.assertAll(channel_map == [0, 0, 1, 1]) + + # Test with custom attribute name + @dataclass + class CustomAttributes: + electrode_id: int + + attrs_custom = [ + CustomAttributes(2), + CustomAttributes(2), + CustomAttributes(3), + ] + sd_custom = SpikeData( + [[0], [10], [5]], length=50, neuron_attributes=attrs_custom + ) + channel_map_custom = sd_custom.get_channel_map(attribute_name="electrode_id") + self.assertIsNotNone(channel_map_custom) + self.assertAll(channel_map_custom == [2, 2, 3]) + + # Test with raw_data + raw_data = np.random.rand(2, 100) # 2 channels, 100 time points + sd_raw = SpikeData( + [[0, 20], [10, 30]], + length=50, + raw_data=raw_data, + raw_time=np.arange(100) / 2.0, + ) + channel_map_raw = sd_raw.get_channel_map(from_raw_data=True) + self.assertIsNotNone(channel_map_raw) + self.assertAll(channel_map_raw == [0, 1]) # Each neuron = one channel + + # Test with metadata + sd_meta = SpikeData( + [[0, 20], [10, 30]], length=50, metadata={"channel_map": [0, 1]} + ) + channel_map_meta = sd_meta.get_channel_map() + self.assertIsNotNone(channel_map_meta) + self.assertAll(channel_map_meta == [0, 1]) + + # Test when no mapping can be determined + sd_no_map = SpikeData([[0, 20], [10, 30]], length=50) + channel_map_none = sd_no_map.get_channel_map() + self.assertIsNone(channel_map_none) + + # Test that channel_raster automatically uses get_channel_map + channel_raster_auto = sd.channel_raster(bin_size=10) + self.assertEqual(channel_raster_auto.shape[0], 2) # 2 channels + + # Test that it raises error when no mapping available + with self.assertRaises(ValueError): + sd_no_map.channel_raster(bin_size=10) From eefe1eb014706efca4e9ff55396caeb42afe2a5f Mon Sep 17 00:00:00 2001 From: Kamran Hussain Date: Sun, 9 Nov 2025 22:32:48 -0800 Subject: [PATCH 3/5] update channel raster --- spikedata/spikedata.py | 108 ++++++++++++++++++++++++++++++++--------- 1 file changed, 86 insertions(+), 22 deletions(-) diff --git a/spikedata/spikedata.py b/spikedata/spikedata.py index a8e590d..319453e 100644 --- a/spikedata/spikedata.py +++ b/spikedata/spikedata.py @@ -574,6 +574,7 @@ def channel_raster( binary=False, attribute_name=None, from_raw_data=False, + expected_num_channels=None, ): """ Create a raster organized by channels instead of individual neurons. @@ -603,13 +604,19 @@ def channel_raster( from_raw_data : bool, default=False If True and raw_data is available, derive mapping from raw_data shape. Passed to get_channel_map(). + expected_num_channels : int, optional + If provided, the output raster will be padded or trimmed to have exactly + this many channels. If the raster has fewer channels, zeros are appended. + If it has more channels, excess channels are trimmed. This is useful when + you need the raster to match a specific channel count from external data. Returns ------- numpy.ndarray or scipy.sparse.csr_array - Array of shape (C, T) where C is the number of channels and T is the - number of time bins. Entry (c, t) indicates spike activity for channel c - in time bin t (either count or binary mask depending on `binary` parameter). + Array of shape (C, T) where C is the number of channels (or expected_num_channels + if provided) and T is the number of time bins. Entry (c, t) indicates spike + activity for channel c in time bin t (either count or binary mask depending on + `binary` parameter). Raises ------ @@ -626,6 +633,9 @@ def channel_raster( >>> >>> # Derive from raw_data shape >>> channel_raster = sd.channel_raster(bin_size=10, from_raw_data=True) + >>> + >>> # Ensure output has exactly 128 channels (pad or trim as needed) + >>> channel_raster = sd.channel_raster(bin_size=10, expected_num_channels=128) """ # Automatically derive channel mapping channel_map = self.get_channel_map( @@ -701,6 +711,41 @@ def channel_raster( else: channel_raster[channel_idx] += neuron_raster[neuron_idx] + # Handle expected_num_channels: pad or trim to match + if expected_num_channels is not None: + current_channels = channel_raster.shape[0] + if current_channels < expected_num_channels: + # Pad with zeros + if sparse.issparse(channel_raster): + # Convert to dense, pad, then convert back + channel_raster_dense = channel_raster.toarray() + padding_shape = ( + expected_num_channels - current_channels, + channel_raster_dense.shape[1], + ) + padding = np.zeros(padding_shape, dtype=channel_raster_dense.dtype) + channel_raster_dense = np.concatenate( + [channel_raster_dense, padding], axis=0 + ) + channel_raster = ( + sparse.csr_array(channel_raster_dense) + if sparse_output + else channel_raster_dense + ) + else: + padding_shape = ( + expected_num_channels - current_channels, + channel_raster.shape[1], + ) + padding = np.zeros(padding_shape, dtype=channel_raster.dtype) + channel_raster = np.concatenate([channel_raster, padding], axis=0) + elif current_channels > expected_num_channels: + # Trim excess channels + if sparse.issparse(channel_raster): + channel_raster = channel_raster[:expected_num_channels, :] + else: + channel_raster = channel_raster[:expected_num_channels, :] + if not sparse_output and sparse.issparse(channel_raster): return channel_raster.toarray() elif sparse_output and not sparse.issparse(channel_raster): @@ -709,7 +754,12 @@ def channel_raster( return channel_raster def sparse_channel_raster( - self, bin_size=20.0, binary=False, attribute_name=None, from_raw_data=False + self, + bin_size=20.0, + binary=False, + attribute_name=None, + from_raw_data=False, + expected_num_channels=None, ): """ Create a sparse raster organized by channels instead of individual neurons. @@ -723,6 +773,7 @@ def sparse_channel_raster( binary=binary, attribute_name=attribute_name, from_raw_data=from_raw_data, + expected_num_channels=expected_num_channels, ) def frame_channel_raster( @@ -732,6 +783,7 @@ def frame_channel_raster( binary=False, attribute_name=None, from_raw_data=False, + expected_num_channels=None, ): """ Create a frame-level raster organized by channels. @@ -755,13 +807,18 @@ def frame_channel_raster( from_raw_data : bool, default=False If True and raw_data is available, derive mapping from raw_data shape. Passed to get_channel_map(). + expected_num_channels : int, optional + If provided, the output raster will be padded or trimmed to have exactly + this many channels. If the raster has fewer channels, zeros are appended. + If it has more channels, excess channels are trimmed. This is useful when + you need the raster to match a specific channel count from external data. Returns ------- numpy.ndarray or scipy.sparse.csr_array - Array of shape (C, F) where C is the number of channels and F is the - number of frames. Entry (c, f) indicates spike activity for channel c - in frame f. + Array of shape (C, F) where C is the number of channels (or expected_num_channels + if provided) and F is the number of frames. Entry (c, f) indicates spike + activity for channel c in frame f. Raises ------ @@ -775,6 +832,7 @@ def frame_channel_raster( binary, attribute_name=attribute_name, from_raw_data=from_raw_data, + expected_num_channels=expected_num_channels, ) def get_channel_map(self, attribute_name=None, from_raw_data=False): @@ -782,10 +840,11 @@ def get_channel_map(self, attribute_name=None, from_raw_data=False): Derive or extract the channel-to-neuron mapping from available data. This method attempts to automatically determine the channel mapping using - multiple strategies: - 1. From neuron_attributes if an attribute name is provided or common names exist - 2. From raw_data shape if it has channel structure - 3. Returns None if no mapping can be determined + multiple strategies (in order of priority): + 1. From metadata if 'channel_map' is stored (highest priority, most reliable) + 2. From neuron_attributes if an attribute name is provided or common names exist + 3. From raw_data shape if it has channel structure + 4. Returns None if no mapping can be determined Parameters ---------- @@ -812,8 +871,22 @@ def get_channel_map(self, attribute_name=None, from_raw_data=False): >>> >>> # Derive from raw_data shape >>> channel_map = sd.get_channel_map(from_raw_data=True) + >>> + >>> # Use channel_map from metadata (highest priority) + >>> sd.metadata['channel_map'] = [0, 0, 1, 1, 2, 2] + >>> channel_map = sd.get_channel_map() """ - # Strategy 1: Try neuron_attributes + # Strategy 1: Check metadata first (highest priority, most reliable) + # This allows external code to explicitly set the channel map + if hasattr(self, "metadata") and isinstance(self.metadata, dict): + if "channel_map" in self.metadata: + channel_map = self.metadata["channel_map"] + if isinstance(channel_map, (list, np.ndarray)): + channel_map = np.asarray(channel_map) + if len(channel_map) == self.N: + return channel_map + + # Strategy 2: Try neuron_attributes if self.neuron_attributes is not None: # Try provided attribute name or common names candidates = [] @@ -844,7 +917,7 @@ def get_channel_map(self, attribute_name=None, from_raw_data=False): except AttributeError: continue - # Strategy 2: Derive from raw_data shape + # Strategy 3: Derive from raw_data shape if from_raw_data and hasattr(self, "raw_data") and self.raw_data.size > 0: raw_shape = self.raw_data.shape # If raw_data has shape (channels, time) or (channels, ...), @@ -867,15 +940,6 @@ def get_channel_map(self, attribute_name=None, from_raw_data=False): channel_map.extend([n_channels - 1] * remainder) return np.asarray(channel_map) - # Strategy 3: Check metadata - if hasattr(self, "metadata") and isinstance(self.metadata, dict): - if "channel_map" in self.metadata: - channel_map = self.metadata["channel_map"] - if isinstance(channel_map, (list, np.ndarray)): - channel_map = np.asarray(channel_map) - if len(channel_map) == self.N: - return channel_map - return None def interspike_intervals(self): From dce7ba1c05bd2bd57ea8b554b7fb2a89a74aacd7 Mon Sep 17 00:00:00 2001 From: Kamran Hussain Date: Mon, 10 Nov 2025 00:14:45 -0800 Subject: [PATCH 4/5] update raster --- spikedata/spikedata.py | 28 ++-- test_spikedata.py | 310 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 324 insertions(+), 14 deletions(-) diff --git a/spikedata/spikedata.py b/spikedata/spikedata.py index 319453e..f59c6a8 100644 --- a/spikedata/spikedata.py +++ b/spikedata/spikedata.py @@ -717,22 +717,22 @@ def channel_raster( if current_channels < expected_num_channels: # Pad with zeros if sparse.issparse(channel_raster): - # Convert to dense, pad, then convert back - channel_raster_dense = channel_raster.toarray() - padding_shape = ( - expected_num_channels - current_channels, - channel_raster_dense.shape[1], - ) - padding = np.zeros(padding_shape, dtype=channel_raster_dense.dtype) - channel_raster_dense = np.concatenate( - [channel_raster_dense, padding], axis=0 - ) - channel_raster = ( - sparse.csr_array(channel_raster_dense) - if sparse_output - else channel_raster_dense + # Use sparse operations to avoid memory-intensive dense conversion + # Create a sparse padding matrix (all zeros, so very memory efficient) + n_bins = channel_raster.shape[1] + padding_channels = expected_num_channels - current_channels + # Create empty sparse matrix for padding (no data stored since all zeros) + padding = sparse.csr_array( + (padding_channels, n_bins), dtype=channel_raster.dtype ) + # Vertically stack the original raster with padding + # Use vstack to combine sparse matrices efficiently + channel_raster = sparse.vstack([channel_raster, padding]) + # Convert to dense only if sparse_output is False + if not sparse_output: + channel_raster = channel_raster.toarray() else: + # Dense case: pad with zeros padding_shape = ( expected_num_channels - current_channels, channel_raster.shape[1], diff --git a/test_spikedata.py b/test_spikedata.py index 4f729d4..4e260a1 100644 --- a/test_spikedata.py +++ b/test_spikedata.py @@ -884,3 +884,313 @@ class CustomAttributes: # Test that it raises error when no mapping available with self.assertRaises(ValueError): sd_no_map.channel_raster(bin_size=10) + + def test_expected_num_channels_padding(self): + """Test expected_num_channels parameter with padding (fewer channels -> more channels)""" + @dataclass + class ChannelAttributes: + channel_id: int + + # Create data with 3 neurons mapping to 2 channels + attrs = [ + ChannelAttributes(0), + ChannelAttributes(0), + ChannelAttributes(1), + ] + sd = SpikeData( + [[0, 10], [5, 15], [8, 18]], length=20, neuron_attributes=attrs + ) + + # Request 5 channels (padding needed) + channel_raster = sd.channel_raster( + bin_size=10, expected_num_channels=5 + ) + self.assertEqual(channel_raster.shape[0], 5) # Should be padded to 5 channels + self.assertEqual(channel_raster.shape[1], 2) # 2 bins + + # First 2 channels should have data, last 3 should be zeros + self.assertGreater(channel_raster[0].sum(), 0) + self.assertGreater(channel_raster[1].sum(), 0) + self.assertAll(channel_raster[2:] == 0) # Padded channels should be zero + + def test_expected_num_channels_trimming(self): + """Test expected_num_channels parameter with trimming (more channels -> fewer channels)""" + @dataclass + class ChannelAttributes: + channel_id: int + + # Create data with 5 neurons mapping to 5 channels + attrs = [ + ChannelAttributes(0), + ChannelAttributes(1), + ChannelAttributes(2), + ChannelAttributes(3), + ChannelAttributes(4), + ] + sd = SpikeData( + [[0], [5], [10], [15], [20]], length=25, neuron_attributes=attrs + ) + + # Request 3 channels (trimming needed) + channel_raster = sd.channel_raster( + bin_size=10, expected_num_channels=3 + ) + self.assertEqual(channel_raster.shape[0], 3) # Should be trimmed to 3 channels + self.assertEqual(channel_raster.shape[1], 3) # 3 bins + + # Should only have data from first 3 channels + self.assertGreater(channel_raster[0].sum(), 0) + self.assertGreater(channel_raster[1].sum(), 0) + self.assertGreater(channel_raster[2].sum(), 0) + + def test_expected_num_channels_exact_match(self): + """Test expected_num_channels when it matches the actual number of channels""" + @dataclass + class ChannelAttributes: + channel_id: int + + attrs = [ + ChannelAttributes(0), + ChannelAttributes(1), + ChannelAttributes(2), + ] + sd = SpikeData( + [[0], [5], [10]], length=15, neuron_attributes=attrs + ) + + # Request exactly 3 channels (no padding/trimming needed) + channel_raster = sd.channel_raster( + bin_size=10, expected_num_channels=3 + ) + self.assertEqual(channel_raster.shape[0], 3) + # Should be identical to without expected_num_channels + channel_raster_no_param = sd.channel_raster(bin_size=10) + self.assertAll(channel_raster == channel_raster_no_param) + + def test_large_channel_count_sparse(self): + """Test with large channel counts (like Maxwell 26.4k channels) using sparse operations""" + # Simulate a large number of channels (use smaller number for testing but test the logic) + num_channels = 1000 # Use 1000 for testing, but logic should scale to 26.4k + num_neurons = 100 # Fewer neurons than channels + + # Create channel map: map neurons to first few channels + channel_map = np.arange(num_neurons) % (num_channels // 10) + metadata = {"channel_map": channel_map.tolist()} + + # Create spike data with neurons + train = [[i * 10 + j * 0.1 for j in range(5)] for i in range(num_neurons)] + sd = SpikeData(train, length=1000, metadata=metadata) + + # Test sparse channel raster with large expected_num_channels + # This should use sparse operations and be memory efficient + channel_raster_sparse = sd.channel_raster( + bin_size=10, + sparse_output=True, + expected_num_channels=num_channels, + ) + self.assertTrue(sparse.issparse(channel_raster_sparse)) + self.assertEqual(channel_raster_sparse.shape[0], num_channels) + # Should have data in first few channels, zeros in rest + self.assertGreater(channel_raster_sparse[: num_channels // 10].sum(), 0) + # Rest should be zeros (sparse, so sum should be efficient) + self.assertEqual(channel_raster_sparse[num_channels // 10 :].sum(), 0) + + def test_large_channel_count_dense(self): + """Test with large channel counts using dense output""" + num_channels = 500 # Use smaller number for dense testing + num_neurons = 50 + + channel_map = np.arange(num_neurons) % (num_channels // 5) + metadata = {"channel_map": channel_map.tolist()} + + train = [[i * 10 + j * 0.1 for j in range(3)] for i in range(num_neurons)] + sd = SpikeData(train, length=600, metadata=metadata) + + # Test dense channel raster with large expected_num_channels + channel_raster_dense = sd.channel_raster( + bin_size=10, + sparse_output=False, + expected_num_channels=num_channels, + ) + self.assertFalse(sparse.issparse(channel_raster_dense)) + self.assertEqual(channel_raster_dense.shape[0], num_channels) + # First channels should have data + self.assertGreater(channel_raster_dense[: num_channels // 5].sum(), 0) + # Padded channels should be zeros + self.assertAll(channel_raster_dense[num_channels // 5 :] == 0) + + def test_frame_channel_raster_expected_num_channels(self): + """Test frame_channel_raster with expected_num_channels parameter""" + @dataclass + class ChannelAttributes: + channel_id: int + + attrs = [ + ChannelAttributes(0), + ChannelAttributes(1), + ] + sd = SpikeData( + [[0, 33.33], [16.67, 50]], length=100, neuron_attributes=attrs + ) + + # Request 5 channels (padding needed) + frame_raster = sd.frame_channel_raster( + frame_rate_hz=30.0, expected_num_channels=5 + ) + self.assertEqual(frame_raster.shape[0], 5) + self.assertEqual(frame_raster.shape[1], 3) # 3 frames at 30 fps + + # First 2 channels should have data, last 3 should be zeros + self.assertGreater(frame_raster[0].sum(), 0) + self.assertGreater(frame_raster[1].sum(), 0) + self.assertAll(frame_raster[2:] == 0) + + def test_expected_num_channels_binary_mode(self): + """Test expected_num_channels with binary mode""" + @dataclass + class ChannelAttributes: + channel_id: int + + attrs = [ + ChannelAttributes(0), + ChannelAttributes(1), + ] + sd = SpikeData( + [[0, 10, 20], [5, 15]], length=25, neuron_attributes=attrs + ) + + # Test binary mode with padding + channel_raster_binary = sd.channel_raster( + bin_size=10, binary=True, expected_num_channels=5 + ) + self.assertEqual(channel_raster_binary.shape[0], 5) + # Binary mode: values should only be 0 or 1 + self.assertAll((channel_raster_binary == 0) | (channel_raster_binary == 1)) + # Padded channels should be zeros + self.assertAll(channel_raster_binary[2:] == 0) + + def test_expected_num_channels_sparse_channel_raster(self): + """Test sparse_channel_raster with expected_num_channels""" + @dataclass + class ChannelAttributes: + channel_id: int + + attrs = [ + ChannelAttributes(0), + ChannelAttributes(1), + ] + sd = SpikeData( + [[0, 10], [5, 15]], length=20, neuron_attributes=attrs + ) + + # Test sparse version with padding + sparse_raster = sd.sparse_channel_raster( + bin_size=10, expected_num_channels=5 + ) + self.assertTrue(sparse.issparse(sparse_raster)) + self.assertEqual(sparse_raster.shape[0], 5) + # Convert to dense to check padding + dense_raster = sparse_raster.toarray() + self.assertAll(dense_raster[2:] == 0) + + def test_expected_num_channels_zero_padding(self): + """Test that padding with zeros works correctly""" + @dataclass + class ChannelAttributes: + channel_id: int + + attrs = [ChannelAttributes(0)] + sd = SpikeData([[0, 10, 20]], length=30, neuron_attributes=attrs) + + # Request 10 channels (1 real + 9 padded) + channel_raster = sd.channel_raster( + bin_size=10, expected_num_channels=10 + ) + self.assertEqual(channel_raster.shape[0], 10) + # First channel should have spikes + self.assertGreater(channel_raster[0].sum(), 0) + # All other channels should be zeros + self.assertAll(channel_raster[1:] == 0) + + def test_expected_num_channels_metadata_priority(self): + """Test that metadata channel_map is used correctly with expected_num_channels""" + # Create spike data with metadata channel_map + metadata = {"channel_map": [0, 0, 1, 1, 2, 2]} # 6 neurons -> 3 channels + train = [[i * 10] for i in range(6)] + sd = SpikeData(train, length=60, metadata=metadata) + + # Request 5 channels (3 real + 2 padded) + channel_raster = sd.channel_raster( + bin_size=10, expected_num_channels=5 + ) + self.assertEqual(channel_raster.shape[0], 5) + # First 3 channels should have data + self.assertGreater(channel_raster[:3].sum(), 0) + # Last 2 should be zeros + self.assertAll(channel_raster[3:] == 0) + + def test_expected_num_channels_edge_cases(self): + """Test edge cases for expected_num_channels""" + @dataclass + class ChannelAttributes: + channel_id: int + + attrs = [ + ChannelAttributes(0), + ChannelAttributes(1), + ] + sd = SpikeData( + [[0, 10], [5, 15]], length=20, neuron_attributes=attrs + ) + + # Test with expected_num_channels=1 (trimming to single channel) + channel_raster_1 = sd.channel_raster( + bin_size=10, expected_num_channels=1 + ) + self.assertEqual(channel_raster_1.shape[0], 1) + + # Test with expected_num_channels=None (should work normally) + channel_raster_none = sd.channel_raster( + bin_size=10, expected_num_channels=None + ) + channel_raster_default = sd.channel_raster(bin_size=10) + self.assertAll(channel_raster_none == channel_raster_default) + + def test_large_channel_count_memory_efficiency(self): + """Test that large channel counts use sparse operations efficiently""" + # Simulate Maxwell-like scenario: many neurons but fewer unique channels + num_neurons = 1000 + num_unique_channels = 100 + num_expected_channels = 26400 # Maxwell-like channel count + + # Map neurons to channels (many neurons per channel) + channel_map = np.arange(num_neurons) % num_unique_channels + metadata = {"channel_map": channel_map.tolist()} + + # Create spike trains + train = [ + [i * 0.1 + j * 0.01 for j in range(10)] for i in range(num_neurons) + ] + sd = SpikeData(train, length=100, metadata=metadata) + + # Test sparse output - should be memory efficient + sparse_raster = sd.channel_raster( + bin_size=1.0, + sparse_output=True, + expected_num_channels=num_expected_channels, + ) + self.assertTrue(sparse.issparse(sparse_raster)) + self.assertEqual(sparse_raster.shape[0], num_expected_channels) + # Should only store non-zero values (first 100 channels have data) + self.assertGreater(sparse_raster[:num_unique_channels].sum(), 0) + self.assertEqual(sparse_raster[num_unique_channels:].sum(), 0) + + # Test that converting to dense works (but may be memory intensive) + # Only test with smaller expected_num_channels to avoid memory issues + dense_raster_small = sd.channel_raster( + bin_size=1.0, + sparse_output=False, + expected_num_channels=200, # Smaller for testing + ) + self.assertEqual(dense_raster_small.shape[0], 200) + self.assertAll(dense_raster_small[num_unique_channels:] == 0) From a7f20ae40384c725154ab06193551362770ddf10 Mon Sep 17 00:00:00 2001 From: Kamran Hussain Date: Sat, 13 Dec 2025 11:30:09 -0800 Subject: [PATCH 5/5] address review comments, remove sparse data check --- spikedata/spikedata.py | 111 +++++++++++++++-------------------------- 1 file changed, 41 insertions(+), 70 deletions(-) diff --git a/spikedata/spikedata.py b/spikedata/spikedata.py index f59c6a8..0e30847 100644 --- a/spikedata/spikedata.py +++ b/spikedata/spikedata.py @@ -660,85 +660,57 @@ def channel_raster( n_channels = len(unique_channels) channel_to_idx = {ch: idx for idx, ch in enumerate(unique_channels)} - # Get neuron-level raster + # Get neuron-level raster (always sparse) neuron_raster = self.sparse_raster(bin_size) # Aggregate by channel - if sparse.issparse(neuron_raster): - # Convert to COO format for easier manipulation - neuron_raster_coo = neuron_raster.tocoo() - # Map neuron indices to channel indices - channel_indices = np.array( - [ - channel_to_idx.get(channel_map[neuron_idx], -1) - for neuron_idx in neuron_raster_coo.row - ] - ) - # Filter out invalid channels - valid_mask = channel_indices >= 0 - channel_indices = channel_indices[valid_mask] - time_indices = neuron_raster_coo.col[valid_mask] - values = neuron_raster_coo.data[valid_mask] - - # Aggregate spikes for the same (channel, time) pair - n_bins = neuron_raster.shape[1] - if binary: - # Binary mask: just mark presence - # Use a temporary dense array to aggregate, then convert to sparse - temp_dense = np.zeros((n_channels, n_bins), dtype=int) - for ch_idx, t_idx in zip(channel_indices, time_indices): - temp_dense[ch_idx, t_idx] = 1 - channel_raster = sparse.csr_array(temp_dense) - else: - # Count spikes - need to sum values for same (channel, time) pairs - # Use a temporary dense array to aggregate properly - temp_dense = np.zeros((n_channels, n_bins), dtype=float) - for ch_idx, t_idx, val in zip(channel_indices, time_indices, values): - temp_dense[ch_idx, t_idx] += val - channel_raster = sparse.csr_array(temp_dense) + # Convert to COO format for easier manipulation + neuron_raster_coo = neuron_raster.tocoo() + # Map neuron indices to channel indices + channel_indices = np.array( + [ + channel_to_idx.get(channel_map[neuron_idx], -1) + for neuron_idx in neuron_raster_coo.row + ] + ) + # Filter out invalid channels + valid_mask = channel_indices >= 0 + channel_indices = channel_indices[valid_mask] + time_indices = neuron_raster_coo.col[valid_mask] + values = neuron_raster_coo.data[valid_mask] + + # Aggregate spikes for the same (channel, time) pair + n_bins = neuron_raster.shape[1] + if binary: + # Binary mask: just mark presence + # Use a temporary dense array to aggregate, then convert to sparse + temp_dense = np.zeros((n_channels, n_bins), dtype=int) + for ch_idx, t_idx in zip(channel_indices, time_indices): + temp_dense[ch_idx, t_idx] = 1 + channel_raster = sparse.csr_array(temp_dense) else: - # Dense case - n_bins = neuron_raster.shape[1] - channel_raster = np.zeros((n_channels, n_bins), dtype=int) - for neuron_idx in range(self.N): - channel_idx = channel_to_idx.get(channel_map[neuron_idx], -1) - if channel_idx >= 0: - if binary: - channel_raster[channel_idx] = ( - channel_raster[channel_idx] - | (neuron_raster[neuron_idx] > 0) - ).astype(int) - else: - channel_raster[channel_idx] += neuron_raster[neuron_idx] + # Count spikes - need to sum values for same (channel, time) pairs + # Use a temporary dense array to aggregate properly + temp_dense = np.zeros((n_channels, n_bins), dtype=float) + for ch_idx, t_idx, val in zip(channel_indices, time_indices, values): + temp_dense[ch_idx, t_idx] += val + channel_raster = sparse.csr_array(temp_dense) # Handle expected_num_channels: pad or trim to match if expected_num_channels is not None: current_channels = channel_raster.shape[0] if current_channels < expected_num_channels: # Pad with zeros + # Since we already create a dense temp array during aggregation, + # we can work with the dense representation directly if sparse.issparse(channel_raster): - # Use sparse operations to avoid memory-intensive dense conversion - # Create a sparse padding matrix (all zeros, so very memory efficient) - n_bins = channel_raster.shape[1] - padding_channels = expected_num_channels - current_channels - # Create empty sparse matrix for padding (no data stored since all zeros) - padding = sparse.csr_array( - (padding_channels, n_bins), dtype=channel_raster.dtype - ) - # Vertically stack the original raster with padding - # Use vstack to combine sparse matrices efficiently - channel_raster = sparse.vstack([channel_raster, padding]) - # Convert to dense only if sparse_output is False - if not sparse_output: - channel_raster = channel_raster.toarray() - else: - # Dense case: pad with zeros - padding_shape = ( - expected_num_channels - current_channels, - channel_raster.shape[1], - ) - padding = np.zeros(padding_shape, dtype=channel_raster.dtype) - channel_raster = np.concatenate([channel_raster, padding], axis=0) + channel_raster = channel_raster.toarray() + padding_shape = ( + expected_num_channels - current_channels, + channel_raster.shape[1], + ) + padding = np.zeros(padding_shape, dtype=channel_raster.dtype) + channel_raster = np.concatenate([channel_raster, padding], axis=0) elif current_channels > expected_num_channels: # Trim excess channels if sparse.issparse(channel_raster): @@ -1666,8 +1638,7 @@ def butter_filter( if lowcut >= highcut: raise ValueError("lowcut must be smaller than highcut") filter_type = "bandpass" - band = [lowcut, highcut] - Wn = [e / fs * 2 for e in band] + Wn = [lowcut / fs * 2, highcut / fs * 2] filter_coeff = signal.iirfilter( order, Wn, analog=False, btype=filter_type, output="sos"