From d1e40f6a59537d200a633419cc2d553eefd8c046 Mon Sep 17 00:00:00 2001 From: Chadwick Boulay Date: Wed, 21 Oct 2026 13:50:43 -0400 Subject: [PATCH] Blackrock: compare NEV spike and event times as seconds, select segments by id At the 1 ns timestamp resolution of PTP files, spikes and events at the edge of a segment were kept or dropped by rounding: the segment bounds went through float seconds back to ticks, and np.searchsorted compared uint64 timestamps with a Python int as float64. spike_count and event_count used the same window and could disagree with the timestamps returned. Without t_start/t_stop, items are selected by segment id alone, as the getters already did, and the counts select the same way. Explicit bounds are compared with the times neo reports (timestamps as float64 seconds), so a bound taken from neo, such as segment_t_start, selects the item it came from. t_stop stays exclusive; documented, with the advice to pass t_stop + 1e-5 to include an item at t_stop. --- neo/rawio/blackrockrawio.py | 76 ++++++++++------------- neo/test/rawiotest/test_blackrockrawio.py | 29 +++++++++ 2 files changed, 63 insertions(+), 42 deletions(-) diff --git a/neo/rawio/blackrockrawio.py b/neo/rawio/blackrockrawio.py index 83d1a2771..cc41bda25 100644 --- a/neo/rawio/blackrockrawio.py +++ b/neo/rawio/blackrockrawio.py @@ -65,7 +65,6 @@ import os import re import warnings -import math import numpy as np import quantities as pq @@ -131,6 +130,10 @@ class BlackrockRawIO(BaseRawIO): 2.2 and 2.3 are automatically extracted and the data set is split into different segments. + * The times of spikes and events are compared with `t_start` and `t_stop` as float64 seconds, and + `t_stop` is exclusive. To include a spike at `t_stop`, such as the last of a segment, pass + `t_stop + 1e-5`. + * The Blackrock data format consists not of a single file, but a set of different files. This constructor associates itself with a set of files that constitute a common data set. By default, all files belonging to @@ -826,18 +829,15 @@ def _get_blackrock_timestamps(self, block_index, seg_index, i_start, i_stop, str def _spike_count(self, block_index, seg_index, unit_index): channel_id, unit_id = self.internal_unit_ids[unit_index] - all_spikes = self.nev_data["Spikes"][0] - mask = (all_spikes["packet_id"] == channel_id) & (all_spikes["unit_class_nb"] == unit_id) - if self._nb_segment == 1: - # very fast - nb = int(np.sum(mask)) - else: - # must clip in time time range - timestamp = all_spikes[mask]["timestamp"] - sl = self._get_timestamp_slice(timestamp, seg_index, None, None) - timestamp = timestamp[sl] - nb = timestamp.size - return nb + all_spikes, event_segment_ids = self.nev_data["Spikes"] + # Selected by segment id, as _get_spike_timestamps does, so the count always equals + # the number of timestamps it returns. + mask = ( + (all_spikes["packet_id"] == channel_id) + & (all_spikes["unit_class_nb"] == unit_id) + & (event_segment_ids == seg_index) + ) + return int(np.sum(mask)) def _get_spike_timestamps(self, block_index, seg_index, unit_index, t_start, t_stop): channel_id, unit_id = self.internal_unit_ids[unit_index] @@ -859,27 +859,27 @@ def _get_spike_timestamps(self, block_index, seg_index, unit_index, t_start, t_s return timestamp def _get_timestamp_slice(self, timestamp, seg_index, t_start, t_stop): - if self._nb_segment > 1: - # we must clip event in seg time limits - if t_start is None: - t_start = self._seg_t_starts[seg_index] - if t_stop is None: - t_stop = self._seg_t_stops[seg_index] + 1 / float(self._nev_basic_header["timestamp_resolution"]) - - if t_start is None: - ind_start = None - else: - ts = math.ceil(t_start * self._nev_basic_header["timestamp_resolution"]) - ind_start = np.searchsorted(timestamp, ts) + """ + Slice of `timestamp` (one segment's, in file order) within [t_start, t_stop), in seconds. - if t_stop is None: - ind_stop = None - else: - ts = int(t_stop * self._nev_basic_header["timestamp_resolution"]) - ind_stop = np.searchsorted(timestamp, ts) # +1 + The segment itself is already selected by segment id, so None means no limit. - sl = slice(ind_start, ind_stop) - return sl + The bounds are compared with the times neo reports (the timestamps in float64 seconds), so a + time obtained from neo, such as `segment_t_start`, selects the item it came from. Converting the + bounds to ticks instead does not work at the 1 ns resolution of PTP files: a float64 number of + seconds (about 1.75e9) is only precise to about 240 ns, so the converted bound can fall a few + ticks past the item and drop it, and np.searchsorted compares uint64 with a Python int as float64. + + `t_stop` is exclusive. To include an item at `t_stop`, such as the last of a segment, + pass `t_stop + 1e-5` (10 us, shorter than the interval between samples at 30 kHz and longer + than the precision of float64 seconds). + """ + if t_start is None and t_stop is None: + return slice(None) + times = timestamp.astype("float64") / self._nev_basic_header["timestamp_resolution"] + ind_start = None if t_start is None else np.searchsorted(times, t_start, side="left") + ind_stop = None if t_stop is None else np.searchsorted(times, t_stop, side="left") + return slice(ind_start, ind_stop) def _rescale_spike_timestamp(self, spike_timestamps, dtype): spike_times = spike_timestamps.astype(dtype) @@ -924,17 +924,9 @@ def _event_count(self, block_index, seg_index, event_channel_index): else: events_data, event_segment_ids = self.nev_data["NonNeural"] ev_dict = self._nonneural_evdicts[self._nev_spec](events_data)[name] + # Selected by segment id, as _get_event_timestamps does. mask = ev_dict["mask"] & (event_segment_ids == seg_index) - if self._nb_segment == 1: - # very fast - nb = int(np.sum(mask)) - else: - # must clip in time time range - timestamp = events_data[ev_dict["mask"]]["timestamp"] - sl = self._get_timestamp_slice(timestamp, seg_index, None, None) - timestamp = timestamp[sl] - nb = timestamp.size - return nb + return int(np.sum(mask)) def _get_event_timestamps(self, block_index, seg_index, event_channel_index, t_start, t_stop): name = self.header["event_channels"]["name"][event_channel_index] diff --git a/neo/test/rawiotest/test_blackrockrawio.py b/neo/test/rawiotest/test_blackrockrawio.py index 03c3481ae..409866265 100644 --- a/neo/test/rawiotest/test_blackrockrawio.py +++ b/neo/test/rawiotest/test_blackrockrawio.py @@ -374,6 +374,35 @@ def test_gap_tolerance_ms_parameter(self): segments_strict = reader_strict.segment_count(0) self.assertEqual(segments_strict, 3) # + def test_nev_spikes_at_segment_edges_ptp(self): + """ + With 1 ns timestamps (PTP files), converting the segment bounds to ticks or comparing uint64 with + floats drops spikes at the edge of a segment. The NEV alone of this file has 4 segments. + """ + reader = BlackrockRawIO( + filename=self.get_local_path( + "blackrock/blackrock_ptp_with_missing_samples/Hub1-NWBtestfile_neural_wspikes" + ), + nsx_to_load=[], + ) + reader.parse_header() + all_spikes, segment_ids = reader.nev_data["Spikes"] + for seg_index in range(reader.segment_count(0)): + t_start, t_stop = reader.segment_t_start(0, seg_index), reader.segment_t_stop(0, seg_index) + for unit_index, (channel_id, unit_id) in enumerate(reader.internal_unit_ids): + in_segment = all_spikes["timestamp"][ + (all_spikes["packet_id"] == channel_id) + & (all_spikes["unit_class_nb"] == unit_id) + & (segment_ids == seg_index) + ] + for bounds in ((None, None), (t_start, t_stop + 1e-5)): + spikes = reader.get_spike_timestamps(0, seg_index, unit_index, *bounds) + assert_equal(spikes, in_segment) + self.assertEqual(reader.spike_count(0, seg_index, unit_index), len(in_segment)) + # t_stop is exclusive + exclusive = reader.get_spike_timestamps(0, seg_index, unit_index, t_start, t_stop) + self.assertLessEqual(len(in_segment) - len(exclusive), 1) + if __name__ == "__main__": unittest.main()