From d62f28dc1102f688ececc5a03196a5829ce6424b Mon Sep 17 00:00:00 2001 From: Heberto Mayorquin Date: Mon, 16 Feb 2026 22:32:35 -0600 Subject: [PATCH 1/9] Ad gap logic to nsx formats 2.1 - 3.0 --- neo/rawio/blackrockrawio.py | 679 +++++++++++++--------- neo/test/iotest/test_blackrockio.py | 8 +- neo/test/rawiotest/test_blackrockrawio.py | 70 ++- 3 files changed, 448 insertions(+), 309 deletions(-) diff --git a/neo/rawio/blackrockrawio.py b/neo/rawio/blackrockrawio.py index fddd5386b..6ab438503 100644 --- a/neo/rawio/blackrockrawio.py +++ b/neo/rawio/blackrockrawio.py @@ -310,7 +310,6 @@ def _parse_header(self): self._nsx_spec = {} self._nsx_basic_header = {} self._nsx_ext_header = {} - self._nsx_data_header = {} self._nsx_sampling_frequency = {} # Read headers @@ -325,20 +324,6 @@ def _parse_header(self): sampling_rate = 30_000.0 / nsx_period self._nsx_sampling_frequency[nsx_nb] = float(sampling_rate) - # Parase data packages - for nsx_nb in self._avail_nsx: - - # The only way to know if it is the Precision Time Protocol of file spec 3.0 - # is to check for nanosecond timestamp resolution. - is_ptp_variant = ( - "timestamp_resolution" in self._nsx_basic_header[nsx_nb].dtype.names - and self._nsx_basic_header[nsx_nb]["timestamp_resolution"] == 1_000_000_000 - ) - if is_ptp_variant: - data_header_spec = "3.0-ptp" - else: - data_header_spec = spec_version - # nsx_to_load can be either int, list, 'max', 'all' (aka None) # here make a list only if self.nsx_to_load is None or self.nsx_to_load == "all": @@ -375,43 +360,36 @@ def _parse_header(self): # Remove if raw loading becomes possible # raise IOError("For loading Blackrock file version 2.1 .nev files are required!") - self.nsx_datas = {} - # Keep public attribute for backward compatibility but let's use the private one and maybe deprecate this at some point - self.sig_sampling_rates = { - nsx_number: self._nsx_sampling_frequency[nsx_number] for nsx_number in self.nsx_to_load - } + # Compute session-level data spec (all nsx files share the same spec). + # The PTP distinction comes from timestamp_resolution; if one is PTP, all are. if len(self.nsx_to_load) > 0: - for nsx_nb in self.nsx_to_load: - basic_header = self._nsx_basic_header[nsx_nb] - spec_version = self._nsx_spec[nsx_nb] - # The only way to know if it is the Precision Time Protocol of file spec 3.0 - # is to check for nanosecond timestamp resolution. - is_ptp_variant = ( - "timestamp_resolution" in basic_header.dtype.names - and basic_header["timestamp_resolution"] == 1_000_000_000 - ) - if is_ptp_variant: - data_spec = "3.0-ptp" - else: - data_spec = spec_version - - # Parse data blocks (creates memmap, extracts data+timestamps) - data_blocks = self._parse_nsx_data(data_spec, nsx_nb) + first_nsx = self.nsx_to_load[0] + basic_header = self._nsx_basic_header[first_nsx] + is_ptp = ( + "timestamp_resolution" in basic_header.dtype.names + and basic_header["timestamp_resolution"] == 1_000_000_000 + ) + self._nsx_data_spec = "3.0-ptp" if is_ptp else self._nsx_spec[first_nsx] + else: + # No nsx files to load (nev-only mode) + self._nsx_data_spec = None - # Segment the data (analyzes gaps, reports issues) - segments = self._segment_nsx_data(data_blocks, nsx_nb) + self._segmented_data_headers = {} + if len(self.nsx_to_load) > 0: + for nsx_nb in self.nsx_to_load: + # Parse data headers (file offsets, sample counts, timestamps) + parsed_data_headers = self._parse_nsx_data(self._nsx_data_spec, nsx_nb) - # Store in existing structures for backward compatibility - self._nsx_data_header[nsx_nb] = { - seg_idx: {k: v for k, v in seg.items() if k != "data"} for seg_idx, seg in segments.items() - } - self.nsx_datas[nsx_nb] = {seg_idx: seg["data"] for seg_idx, seg in segments.items()} + # Segment the data (gap detection, groups headers into segments) + segmented_data_headers = self._segment_nsx_data(parsed_data_headers, nsx_nb, self.gap_tolerance_ms) + self._segmented_data_headers[nsx_nb] = segmented_data_headers # Match NSX and NEV segments for v2.3 if self._avail_files["nev"]: self._match_nsx_and_nev_segment_ids(nsx_nb) sr = self._nsx_sampling_frequency[nsx_nb] + spec_version = self._nsx_spec[nsx_nb] if spec_version in ["2.2", "2.3", "3.0"]: ext_header = self._nsx_ext_header[nsx_nb] @@ -449,7 +427,7 @@ def _parse_header(self): signal_channels.append((ch_name, ch_id, sr, sig_dtype, units, gain, offset, stream_id, buffer_id)) # check nb segment per nsx - nb_segments_for_nsx = [len(self.nsx_datas[nsx_nb]) for nsx_nb in self.nsx_to_load] + nb_segments_for_nsx = [len(self._segmented_data_headers[nsx_nb]) for nsx_nb in self.nsx_to_load] if not all(nb == nb_segments_for_nsx[0] for nb in nb_segments_for_nsx): raise NeoReadWriteError("Segment nb not consistent across nsX files") self._nb_segment = nb_segments_for_nsx[0] @@ -457,35 +435,33 @@ def _parse_header(self): self._delete_empty_segments() # t_start/t_stop for segment are given by nsx limits or nev limits - self._sigs_t_starts = {nsx_nb: [] for nsx_nb in self.nsx_to_load} self._seg_t_starts, self._seg_t_stops = [], [] for data_bl in range(self._nb_segment): + t_start = float("inf") t_stop = 0.0 for nsx_nb in self.nsx_to_load: - spec = self._nsx_spec[nsx_nb] - if "timestamp_resolution" in self._nsx_basic_header[nsx_nb].dtype.names: - ts_res = self._nsx_basic_header[nsx_nb]["timestamp_resolution"] - elif spec == "2.1": - ts_res = 30_000 # v2.1 always uses 30kHz timestamp resolution - else: - ts_res = 30_000 - period = self._nsx_basic_header[nsx_nb]["period"] - sec_per_samp = period / 30_000 # Maybe 30_000 should be ['sample_resolution'] - length = self.nsx_datas[nsx_nb][data_bl].shape[0] - timestamps = self._nsx_data_header[nsx_nb][data_bl]["timestamp"] - if timestamps is None: - # V2.1 format has no timestamps - t_start = 0.0 - t_stop = max(t_stop, length / self._nsx_sampling_frequency[nsx_nb]) - elif hasattr(timestamps, "size") and timestamps.size == length: - # FileSpec 3.0 with PTP -- use the per-sample timestamps - t_start = timestamps[0] / ts_res - t_stop = max(t_stop, timestamps[-1] / ts_res + sec_per_samp) + seg = self._segmented_data_headers[nsx_nb][data_bl] + t_start = min(t_start, seg["t_start"]) + nb_pts = seg["nb_data_points"] + sr = self._nsx_sampling_frequency[nsx_nb] + if self._nsx_data_spec == "3.0-ptp": + # PTP: read actual last timestamp (jitter makes t_start + n/sr imprecise) + channels = int(self._nsx_basic_header[nsx_nb]["channel_count"]) + filename = f"{self._filenames['nsx']}.ns{nsx_nb}" + ptp_dt = NSX_DATA_HEADER_TYPES["3.0-ptp"](channels) + memmap = np.memmap( + filename, dtype=ptp_dt, mode="r", + offset=seg["ptp_data_offset"], + shape=(seg["nb_data_points"],), + ) + ts_res = float(self._nsx_basic_header[nsx_nb]["timestamp_resolution"]) + last_ts = float(memmap["timestamps"][-1]) / ts_res + del memmap + seg_t_stop = last_ts + 1.0 / sr else: - # Standard format with scalar timestamp - t_start = timestamps / ts_res - t_stop = max(t_stop, t_start + length / self._nsx_sampling_frequency[nsx_nb]) - self._sigs_t_starts[nsx_nb].append(t_start) + # Standard/v2.1: exact from t_start + nb_pts / sr + seg_t_stop = seg["t_start"] + nb_pts / sr + t_stop = max(t_stop, seg_t_stop) if self._avail_files["nev"]: max_nev_time = 0 @@ -529,7 +505,6 @@ def _parse_header(self): self._seg_t_starts = [v / float(resolution) for k, v in sorted(min_nev_times.items())] self._seg_t_stops = [v / float(resolution) for k, v in sorted(max_nev_times.items())] self._nb_segment = len(self._seg_t_starts) - self._sigs_t_starts = [None] * self._nb_segment # finalize header spike_channels = np.array(spike_channels, dtype=_spike_channel_dtype) @@ -650,22 +625,105 @@ def _segment_t_stop(self, block_index, seg_index): def _get_signal_size(self, block_index, seg_index, stream_index): stream_id = self.header["signal_streams"][stream_index]["id"] nsx_nb = int(stream_id) - memmap_data = self.nsx_datas[nsx_nb][seg_index] - return memmap_data.shape[0] + return self._segmented_data_headers[nsx_nb][seg_index]["nb_data_points"] def _get_signal_t_start(self, block_index, seg_index, stream_index): stream_id = self.header["signal_streams"][stream_index]["id"] nsx_nb = int(stream_id) - return self._sigs_t_starts[nsx_nb][seg_index] + return self._segmented_data_headers[nsx_nb][seg_index]["t_start"] def _get_analogsignal_chunk(self, block_index, seg_index, i_start, i_stop, stream_index, channel_indexes): stream_id = self.header["signal_streams"][stream_index]["id"] nsx_nb = int(stream_id) - memmap_data = self.nsx_datas[nsx_nb][seg_index] + seg = self._segmented_data_headers[nsx_nb][seg_index] + channels = int(self._nsx_basic_header[nsx_nb]["channel_count"]) + filename = f"{self._filenames['nsx']}.ns{nsx_nb}" if channel_indexes is None: channel_indexes = slice(None) - sig_chunk = memmap_data[i_start:i_stop, channel_indexes] - return sig_chunk + + if self._nsx_data_spec == "3.0-ptp": + ptp_dt = NSX_DATA_HEADER_TYPES["3.0-ptp"](channels) + memmap = np.memmap( + filename, dtype=ptp_dt, mode="r", + offset=seg["ptp_data_offset"], + shape=(seg["nb_data_points"],), + ) + return memmap["samples"][i_start:i_stop, channel_indexes] + else: + specs = seg["memmap_specs"] + if len(specs) == 1: + data = np.memmap( + filename, dtype="int16", mode="r", + offset=specs[0]["offset"], + shape=(specs[0]["num_samples"], channels), + ) + return data[i_start:i_stop, channel_indexes] + else: + return self._read_multi_block_chunk( + filename, specs, channels, i_start, i_stop, channel_indexes, + ) + + def _read_multi_block_chunk(self, filename, memmap_specs, channels, + i_start, i_stop, channel_indexes): + """ + Read a chunk of analog signal data that spans multiple data blocks + merged into a single segment. + + In the standard format (v2.2, v2.3, v3.0), consecutive data blocks + without significant gaps are merged into one segment by the segmenter. + Each block is stored at a different file offset, so reading a contiguous + sample range [i_start, i_stop) may require stitching data from several + blocks. This method creates temporary memmaps only for the blocks that + overlap the requested range, slices each one, and concatenates the + results. + + Parameters + ---------- + filename : str + Path to the NSx file. + memmap_specs : list[dict] + Each dict has "offset" (byte offset in file) and "num_samples" + (number of samples in that block). Ordered sequentially within + the segment. + channels : int + Number of channels (columns in the int16 data matrix). + i_start : int or None + First sample index (segment-relative). None means 0. + i_stop : int or None + Stop sample index (exclusive, segment-relative). None means end. + channel_indexes : slice or array-like + Which channels to return. + + Returns + ------- + np.ndarray + Signal data of shape (i_stop - i_start, len(channel_indexes)), + dtype int16. Returns a memmap view when only one block is touched, + otherwise a copied concatenation. + """ + total_samples = sum(spec["num_samples"] for spec in memmap_specs) + if i_start is None: + i_start = 0 + if i_stop is None: + i_stop = total_samples + pieces = [] + cumulative = 0 + for spec in memmap_specs: + block_start = cumulative + block_end = cumulative + spec["num_samples"] + if block_end > i_start and block_start < i_stop: + local_start = max(0, i_start - block_start) + local_stop = min(spec["num_samples"], i_stop - block_start) + data = np.memmap( + filename, dtype="int16", mode="r", + offset=spec["offset"], + shape=(spec["num_samples"], channels), + ) + pieces.append(data[local_start:local_stop, channel_indexes]) + cumulative = block_end + if len(pieces) == 1: + return pieces[0] + return np.concatenate(pieces, axis=0) def _get_blackrock_timestamps(self, block_index, seg_index, i_start, i_stop, stream_index): """ @@ -710,24 +768,26 @@ def _get_blackrock_timestamps(self, block_index, seg_index, i_start, i_stop, str """ stream_id = self.header["signal_streams"][stream_index]["id"] nsx_nb = int(stream_id) + seg = self._segmented_data_headers[nsx_nb][seg_index] - # Resolve None to concrete indices - size = self.nsx_datas[nsx_nb][seg_index].shape[0] + size = seg["nb_data_points"] i_start = i_start if i_start is not None else 0 i_stop = i_stop if i_stop is not None else size - # Check if this segment has per-sample timestamps (PTP format) - raw_timestamps = self._nsx_data_header[nsx_nb][seg_index]["timestamp"] - - if isinstance(raw_timestamps, np.ndarray) and raw_timestamps.size == size: - # PTP: real hardware timestamps + if self._nsx_data_spec == "3.0-ptp": + channels = int(self._nsx_basic_header[nsx_nb]["channel_count"]) + filename = f"{self._filenames['nsx']}.ns{nsx_nb}" + ptp_dt = NSX_DATA_HEADER_TYPES["3.0-ptp"](channels) + memmap = np.memmap( + filename, dtype=ptp_dt, mode="r", + offset=seg["ptp_data_offset"], + shape=(seg["nb_data_points"],), + ) ts_res = float(self._nsx_basic_header[nsx_nb]["timestamp_resolution"]) - return raw_timestamps[i_start:i_stop].astype("float64") / ts_res + return memmap["timestamps"][i_start:i_stop].astype("float64") / ts_res else: - # Non-PTP: reconstruct from t_start + index / sampling_rate - t_start = self._sigs_t_starts[nsx_nb][seg_index] sr = self._nsx_sampling_frequency[nsx_nb] - return t_start + np.arange(i_start, i_stop, dtype="float64") / sr + return seg["t_start"] + np.arange(i_start, i_stop, dtype="float64") / sr def _spike_count(self, block_index, seg_index, unit_index): channel_id, unit_id = self.internal_unit_ids[unit_index] @@ -1073,28 +1133,20 @@ def _parse_nsx_data(self, spec, nsx_nb): Returns ------- - dict - Dictionary mapping block index to block information: - { - block_idx: { - "data": np.ndarray, - View into memory-mapped file with shape (samples, channels) - "timestamps": scalar, np.ndarray, or None, - - Standard format: scalar (one timestamp per block) - - PTP format: array (one timestamp per sample) - - v2.1 format: None (no timestamps) - # Additional metadata as needed - }, - ... - } + list[dict] + Each dict contains: + - "memmap_kwargs": {"offset": int, "num_samples": int} + Byte offset and sample count for creating memmaps at read time. + - "nsx_block_timestamp": scalar (standard only) + Raw timestamp tick from the block header. + - "ptp_timestamps": np.ndarray of uint64 (PTP only) + Per-sample timestamps, temporary (used for gap detection then discarded). Notes ----- - - This function creates the file memmap internally - - Data views are created using np.ndarray with buffer parameter (memory efficient) - - Returned data is NOT YET SEGMENTED (segmentation happens in a separate step) - - For standard format, each block from the file is one dict entry - - For PTP format, all data is in a single block (block_idx=0) + - Returned data is NOT YET SEGMENTED (segmentation happens in _segment_nsx_data) + - For standard format, each data block from the file is one list entry + - For PTP format, the list has a single entry """ if spec == "2.1": return self._parse_nsx_data_v21(nsx_nb) @@ -1114,31 +1166,17 @@ def _parse_nsx_data_v21(self, nsx_nb): Returns ------- - dict - {0: {"data": np.ndarray, "timestamps": None}} + list[dict] + [{"memmap_kwargs": {"offset": int, "num_samples": int}}] """ filename = f"{self._filenames['nsx']}.ns{nsx_nb}" - - # Create file memmap - file_memmap = np.memmap(filename, dtype="uint8", mode="r") - - # Calculate header size and data points for v2.1 channels = int(self._nsx_basic_header[nsx_nb]["channel_count"]) bytes_in_headers = ( self._nsx_basic_header[nsx_nb].dtype.itemsize + self._nsx_ext_header[nsx_nb].dtype.itemsize * channels ) filesize = self._get_file_size(filename) num_samples = int((filesize - bytes_in_headers) / (2 * channels) - 1) - offset = bytes_in_headers - # Create data view into memmap - data = np.ndarray(shape=(num_samples, channels), dtype="int16", buffer=file_memmap, offset=offset) - - return { - 0: { - "data": data, - "timestamps": None, - } - } + return [{"memmap_kwargs": {"offset": bytes_in_headers, "num_samples": num_samples}}] def _parse_nsx_data_v22_v30(self, spec, nsx_nb): """ @@ -1151,25 +1189,18 @@ def _parse_nsx_data_v22_v30(self, spec, nsx_nb): Returns ------- - dict - {block_idx: {"data": np.ndarray, "timestamps": scalar}, ...} + list[dict] + [{"memmap_kwargs": {"offset": int, "num_samples": int}, + "nsx_block_timestamp": scalar}, ...] """ filename = f"{self._filenames['nsx']}.ns{nsx_nb}" - - # Create file memmap - file_memmap = np.memmap(filename, dtype="uint8", mode="r") - - # Get file parameters filesize = self._get_file_size(filename) channels = int(self._nsx_basic_header[nsx_nb]["channel_count"]) current_offset = int(self._nsx_basic_header[nsx_nb]["bytes_in_headers"]) - data_blocks = {} - block_idx = 0 + parsed_data_headers = [] - # Loop through file, reading block headers while current_offset < filesize: - # Read header at current position header = self._read_nsx_dataheader(spec, nsx_nb, current_offset) if header["header_flag"] != 1: @@ -1180,22 +1211,16 @@ def _parse_nsx_data_v22_v30(self, spec, nsx_nb): num_samples = int(header["nb_data_points"]) data_offset = current_offset + header.dtype.itemsize - timestamp = header["timestamp"] - - # Create data view into memmap for this block - data = np.ndarray(shape=(num_samples, channels), dtype="int16", buffer=file_memmap, offset=data_offset) - data_blocks[block_idx] = { - "data": data, - "timestamps": timestamp, - } + parsed_data_headers.append({ + "memmap_kwargs": {"offset": data_offset, "num_samples": num_samples}, + "nsx_block_timestamp": header["timestamp"], + }) - # Jump to next block data_size_bytes = num_samples * channels * 2 # int16 = 2 bytes current_offset = data_offset + data_size_bytes - block_idx += 1 - return data_blocks + return parsed_data_headers def _parse_nsx_data_v30_ptp(self, nsx_nb): """ @@ -1208,36 +1233,35 @@ def _parse_nsx_data_v30_ptp(self, nsx_nb): Returns ------- - dict - {0: {"data": np.ndarray, "timestamps": np.ndarray}} + list[dict] + [{"memmap_kwargs": {"offset": int, "num_samples": int}, + "ptp_timestamps": np.ndarray (uint64)}] """ filename = f"{self._filenames['nsx']}.ns{nsx_nb}" - # Get file parameters filesize = self._get_file_size(filename) header_size = int(self._nsx_basic_header[nsx_nb]["bytes_in_headers"]) channel_count = int(self._nsx_basic_header[nsx_nb]["channel_count"]) - # Create structured memmap (timestamp + samples per packet) ptp_dt = NSX_DATA_HEADER_TYPES["3.0-ptp"](channel_count) npackets = int((filesize - header_size) / np.dtype(ptp_dt).itemsize) - file_memmap = np.memmap(filename, dtype=ptp_dt, shape=npackets, offset=header_size, mode="r") - # Verify this is truly PTP (all packets should have 1 sample) - if not np.all(file_memmap["num_data_points"] == 1): - # Not actually PTP! Fall back to standard format + # Temporary memmap for verification only + temp_memmap = np.memmap(filename, dtype=ptp_dt, shape=npackets, offset=header_size, mode="r") + + if not np.all(temp_memmap["num_data_points"] == 1): + # Not actually PTP -- fall back to standard v3.0 block parsing + del temp_memmap return self._parse_nsx_data_v22_v30("3.0", nsx_nb) - # Extract data and timestamps from structured array - data = file_memmap["samples"] - timestamps = file_memmap["timestamps"] + # Copy timestamps for gap detection (released after segmentation) + ptp_timestamps = np.array(temp_memmap["timestamps"]) + del temp_memmap - return { - 0: { - "data": data, - "timestamps": timestamps, - } - } + return [{ + "memmap_kwargs": {"offset": header_size, "num_samples": npackets}, + "ptp_timestamps": ptp_timestamps, + }] def _format_gap_report(self, gap_indices, timestamps_in_seconds, time_differences, nsx_nb): """ @@ -1281,117 +1305,208 @@ def _format_gap_report(self, gap_indices, timestamps_in_seconds, time_difference + "+-----------------+-----------------------+-----------------------+\n" ) - def _segment_nsx_data(self, data_blocks_dict, nsx_nb): + def _classify_gaps(self, deviations, sampling_rate, gap_tolerance_ms, nsx_nb, + timestamps_for_report, intervals_for_report): + """ + Classify deviations from expected sample intervals as gaps. + + Used by both standard (block-level) and PTP (sample-level) segmentation. + Forward gaps (pauses, dropped samples) are filtered by gap_tolerance_ms. + Backward jumps (abnormal condition) always create segment boundaries. + + Parameters + ---------- + deviations : np.ndarray + Difference between actual inter-sample/inter-block intervals and the + expected interval (1/sampling_rate). Positive = forward gap, + negative = backward jump. + sampling_rate : float + Sampling rate in Hz. + gap_tolerance_ms : float or None + If None, raises ValueError when gaps are detected. + Otherwise, forward gaps smaller than this (in ms) are ignored. + nsx_nb : int + NSX file number (for error reporting). + timestamps_for_report : np.ndarray + Timestamps in seconds (for gap report formatting). + intervals_for_report : np.ndarray + Time differences between consecutive timestamps (for gap report). + + Returns + ------- + np.ndarray + Indices into deviations where significant gaps occur. """ - Segment NSX data based on timestamp gaps. + half_sample_period = 0.5 / sampling_rate - Takes the data blocks returned by _parse_nsx_data() and creates segments. - Segmentation logic depends on the file format: + forward_mask = deviations > half_sample_period + backward_mask = deviations < -half_sample_period - - Standard format (multiple blocks): Each block IS a segment - - PTP format (single block with timestamp array): Detect gaps in timestamps - - V2.1 format (no timestamps): Single segment + if not np.any(forward_mask | backward_mask): + return np.array([], dtype=np.intp) + + all_gap_indices = np.flatnonzero(forward_mask | backward_mask) + gap_report = self._format_gap_report( + all_gap_indices, timestamps_for_report, intervals_for_report, nsx_nb, + ) + + if gap_tolerance_ms is None: + raise ValueError( + f"Detected {len(all_gap_indices)} timestamp gaps in ns{nsx_nb}.\n" + f"{gap_report}\n" + f"Provide gap_tolerance_ms to segment at gaps." + ) + + gap_tolerance_s = gap_tolerance_ms / 1000.0 + forward_indices = np.flatnonzero(forward_mask) + backward_indices = np.flatnonzero(backward_mask) + significant_forward = forward_indices[deviations[forward_indices] > gap_tolerance_s] + return np.union1d(significant_forward, backward_indices) + + def _segment_nsx_data(self, parsed_data_headers, nsx_nb, gap_tolerance_ms): + """ + Dispatch to the appropriate spec-specific segmentation function. Parameters ---------- - data_blocks_dict : dict - Dictionary from _parse_nsx_data(): - {block_idx: {"data": np.ndarray, "timestamps": scalar/array/None}} + parsed_data_headers : list[dict] + Data headers from _parse_nsx_data(). nsx_nb : int - NSX file number + NSX file number. + gap_tolerance_ms : float or None + Gap tolerance in milliseconds. Returns ------- - dict + list[dict] + Each dict contains: + - "nb_data_points": int + - "t_start": float (seconds) + - "memmap_specs": list of {"offset": int, "num_samples": int} + (standard/v2.1 only) + - "ptp_data_offset": int (PTP only) + - "nsx_block_timestamp": scalar raw ticks (standard only) + """ + data_spec = self._nsx_data_spec + if data_spec == "2.1": + return self._segment_nsx_v21(parsed_data_headers) + elif data_spec in ("2.2", "2.3", "3.0"): + return self._segment_nsx_v22_v30(parsed_data_headers, nsx_nb, gap_tolerance_ms) + elif data_spec == "3.0-ptp": + return self._segment_nsx_ptp(parsed_data_headers, nsx_nb, gap_tolerance_ms) + + def _segment_nsx_v21(self, parsed_data_headers): + """ + Segment v2.1 data. Single entry, no timestamps, single segment. + """ + header = parsed_data_headers[0] + return [ { - seg_idx: { - "data": np.ndarray, - "timestamps": scalar, array, or None, - "nb_data_points": int, - "header": int or None, - "offset_to_data_block": None (deprecated but kept for compatibility) - }, - ... + "memmap_specs": [header["memmap_kwargs"]], + "nb_data_points": header["memmap_kwargs"]["num_samples"], + "t_start": 0.0, } + ] + + def _segment_nsx_v22_v30(self, parsed_data_headers, nsx_nb, gap_tolerance_ms): + """ + Segment standard format data (v2.2, v2.3, v3.0) using block-level gap detection. + + Operates on N block timestamps (typically 1-5 scalars). Consecutive blocks + without significant gaps are merged into a single segment. """ - segments = {} - - # Case 1: Multiple blocks (Standard format) - each block is a segment - if len(data_blocks_dict) > 1: - for block_idx, block_info in data_blocks_dict.items(): - segments[block_idx] = { - "data": block_info["data"], - "timestamp": block_info["timestamps"], # Use singular for backward compatibility - "nb_data_points": block_info["data"].shape[0], - "header": 1, # Standard format has headers - "offset_to_data_block": None, # Not needed (have data directly) + ts_res = float(self._nsx_basic_header[nsx_nb]["timestamp_resolution"]) + sr = self._nsx_sampling_frequency[nsx_nb] + headers = parsed_data_headers + + if len(headers) == 1: + h = headers[0] + return [ + { + "memmap_specs": [h["memmap_kwargs"]], + "nb_data_points": h["memmap_kwargs"]["num_samples"], + "t_start": float(h["nsx_block_timestamp"]) / ts_res, + "nsx_block_timestamp": h["nsx_block_timestamp"], } + ] - # Case 2: Single block - check if PTP (array timestamps) or simple (no timestamps) - elif len(data_blocks_dict) == 1: - block_info = data_blocks_dict[0] - data = block_info["data"] - timestamps = block_info["timestamps"] - - # PTP format: array of timestamps - need to detect gaps - if isinstance(timestamps, np.ndarray): - # Analyze timestamp gaps - sampling_rate = self._nsx_sampling_frequency[nsx_nb] - - # Detection threshold: use strict 2x sampling period to find ALL gaps - detection_threshold = 2.0 / sampling_rate - - timestamps_sampling_rate = self._nsx_basic_header[nsx_nb]["timestamp_resolution"] - timestamps_in_seconds = timestamps / timestamps_sampling_rate - - time_differences = np.diff(timestamps_in_seconds) - gap_indices = np.argwhere(time_differences > detection_threshold).flatten() - - # If gaps found, check user's tolerance - if len(gap_indices) > 0: - gap_report = self._format_gap_report(gap_indices, timestamps_in_seconds, time_differences, nsx_nb) - - # Error by default - user must opt-in to segmentation - if self.gap_tolerance_ms is None: - raise ValueError( - f"Detected {len(gap_indices)} timestamp gaps in ns{nsx_nb} file.\n" - f"{gap_report}\n" - f"To load this data, provide gap_tolerance_ms parameter to automatically " - f"segment at gaps larger than the specified tolerance." - ) + block_t_starts = np.array([float(h["nsx_block_timestamp"]) / ts_res for h in headers]) + block_sizes = np.array([h["memmap_kwargs"]["num_samples"] for h in headers]) + block_t_ends = block_t_starts + block_sizes / sr - # User provided tolerance - filter gaps and segment - gap_tolerance_s = self.gap_tolerance_ms / 1000.0 - significant_gap_mask = time_differences[gap_indices] > gap_tolerance_s - significant_gap_indices = gap_indices[significant_gap_mask] + inter_block_intervals = block_t_starts[1:] - block_t_ends[:-1] + expected_interval = 1.0 / sr + deviations = inter_block_intervals - expected_interval - # Use significant gaps for segmentation (no warning - user opted in) - gap_indices = significant_gap_indices + gap_boundary_indices = self._classify_gaps( + deviations, sr, gap_tolerance_ms, nsx_nb, + block_t_starts, inter_block_intervals, + ) - # Create segments based on gaps - segment_starts = np.hstack((0, gap_indices + 1)) - segment_boundaries = list(segment_starts) + [len(data)] + # Group consecutive blocks into segments + segments = [] + seg_start = 0 + + for gap_after in sorted(gap_boundary_indices): + seg_headers = headers[seg_start : gap_after + 1] + segments.append({ + "memmap_specs": [h["memmap_kwargs"] for h in seg_headers], + "nb_data_points": sum(h["memmap_kwargs"]["num_samples"] for h in seg_headers), + "t_start": block_t_starts[seg_start], + "nsx_block_timestamp": headers[seg_start]["nsx_block_timestamp"], + }) + seg_start = gap_after + 1 + + # Last segment (after last gap, or all blocks if no gaps) + seg_headers = headers[seg_start:] + segments.append({ + "memmap_specs": [h["memmap_kwargs"] for h in seg_headers], + "nb_data_points": sum(h["memmap_kwargs"]["num_samples"] for h in seg_headers), + "t_start": block_t_starts[seg_start], + "nsx_block_timestamp": headers[seg_start]["nsx_block_timestamp"], + }) - for seg_idx, start in enumerate(segment_starts): - end = segment_boundaries[seg_idx + 1] + return segments - segments[seg_idx] = { - "data": data[start:end], - "timestamp": timestamps[start:end], # Use singular for backward compatibility - "nb_data_points": end - start, - "header": None, # PTP has no headers - "offset_to_data_block": None, - } + def _segment_nsx_ptp(self, parsed_data_headers, nsx_nb, gap_tolerance_ms): + """ + Segment PTP format data (v3.0-ptp) using sample-level gap detection. - # V2.1 or single block standard format: no segmentation needed - else: - segments[0] = { - "data": data, - "timestamp": timestamps, # Use singular for backward compatibility - "nb_data_points": data.shape[0], - "header": None, - "offset_to_data_block": None, - } + Operates on per-sample uint64 timestamps. The ptp_timestamps array from + the parser is consumed here and not stored on segments. + """ + ts_res = float(self._nsx_basic_header[nsx_nb]["timestamp_resolution"]) + sr = self._nsx_sampling_frequency[nsx_nb] + channel_count = int(self._nsx_basic_header[nsx_nb]["channel_count"]) + ptp_dt = NSX_DATA_HEADER_TYPES["3.0-ptp"](channel_count) + packet_size = np.dtype(ptp_dt).itemsize + + header = parsed_data_headers[0] + base_offset = header["memmap_kwargs"]["offset"] + num_samples = header["memmap_kwargs"]["num_samples"] + ptp_timestamps = header["ptp_timestamps"] # temporary, consumed here + + all_timestamps = ptp_timestamps.astype("float64") / ts_res + time_diffs = np.diff(all_timestamps) + deviations = time_diffs - 1.0 / sr + + gap_indices = self._classify_gaps( + deviations, sr, gap_tolerance_ms, nsx_nb, + all_timestamps, time_diffs, + ) + + segment_starts = np.concatenate(([0], gap_indices + 1)) + segment_ends = np.concatenate((gap_indices + 1, [num_samples])) + + segments = [] + for seg_index in range(len(segment_starts)): + start = int(segment_starts[seg_index]) + end = int(segment_ends[seg_index]) + segments.append({ + "ptp_data_offset": base_offset + start * packet_size, + "nb_data_points": end - start, + "t_start": float(all_timestamps[start]), + }) return segments @@ -1667,7 +1782,7 @@ def _match_nsx_and_nev_segment_ids(self, nsx_nb): # Only needs to be done for nev version 2.3 if self._nev_spec == "2.3": - nsx_offset = self._nsx_data_header[nsx_nb][0]["timestamp"] + nsx_offset = self._segmented_data_headers[nsx_nb][0]["nsx_block_timestamp"] # Multiples of 1/30.000s that pass between two nsX samples nsx_period = self._nsx_basic_header[nsx_nb]["period"] # NSX segments needed as dict and list @@ -1678,7 +1793,7 @@ def _match_nsx_and_nev_segment_ids(self, nsx_nb): # Nonempty segments are those containing at least 2 samples # These have to be able to be mapped to nev - for k, v in sorted(self._nsx_data_header[nsx_nb].items()): + for k, v in enumerate(self._segmented_data_headers[nsx_nb]): if v["nb_data_points"] > 1: nonempty_nsx_segments[k] = v list_nonempty_nsx_segments.append(v) @@ -1695,7 +1810,7 @@ def _match_nsx_and_nev_segment_ids(self, nsx_nb): # Last timestamp in this nsX segment # Not subtracting nsX offset from end because spike extraction might continue end_of_current_nsx_seg = ( - seg["timestamp"] + seg["nb_data_points"] * self._nsx_basic_header[nsx_nb]["period"] + seg["nsx_block_timestamp"] + seg["nb_data_points"] * self._nsx_basic_header[nsx_nb]["period"] ) mask_after_seg = (ev_ids == i) & (data["timestamp"] > end_of_current_nsx_seg + nsx_period) @@ -1703,7 +1818,7 @@ def _match_nsx_and_nev_segment_ids(self, nsx_nb): # Show warning if spikes do not fit any segment (+- 1 sampling 'tick') # Spike should belong to segment before mask_outside = (ev_ids == i) & ( - data["timestamp"] < int(seg["timestamp"]) - int(nsx_offset) - int(nsx_period) + data["timestamp"] < int(seg["nsx_block_timestamp"]) - int(nsx_offset) - int(nsx_period) ) if len(data[mask_outside]) > 0: @@ -1732,7 +1847,7 @@ def _match_nsx_and_nev_segment_ids(self, nsx_nb): # If reset and no segment detected in nev, then these segments cannot be # distinguished in nev, which is a big problem # XXX 96 is an arbitrary number based on observations in available files - elif list_nonempty_nsx_segments[i + 1]["timestamp"] - nsx_offset <= 96: + elif list_nonempty_nsx_segments[i + 1]["nsx_block_timestamp"] - nsx_offset <= 96: # If not all definitely belong to the next segment, # then it cannot be distinguished where some belong if len(data[ev_ids == i]) != len(data[mask_after_seg]): @@ -1991,43 +2106,35 @@ def _delete_empty_segments(self): """ If there are empty segments (e.g. due to a reset or clock synchronization across two systems), these can be discarded. - If this is done, all the data and data_headers need to be remapped to become a range - starting from 0 again. Nev data are mapped accordingly to stay with their corresponding + Nev data are mapped accordingly to stay with their corresponding segment in the nsX data. """ - # Discard empty segments - removed_seg = [] - for data_bl in range(self._nb_segment): - keep_seg = True - for nsx_nb in self.nsx_to_load: - length = self.nsx_datas[nsx_nb][data_bl].shape[0] - keep_seg = keep_seg and (length >= 2) + # Find which segments to keep (must have >= 2 samples in all nsx files) + keep_mask = [] + for seg_index in range(self._nb_segment): + keep = all( + self._segmented_data_headers[nsx_nb][seg_index]["nb_data_points"] >= 2 + for nsx_nb in self.nsx_to_load + ) + keep_mask.append(keep) - if not keep_seg: - removed_seg.append(data_bl) - for nsx_nb in self.nsx_to_load: - self.nsx_datas[nsx_nb].pop(data_bl) - self._nsx_data_header[nsx_nb].pop(data_bl) - - # Keys need to be increasing from 0 to maximum in steps of 1 - # To ensure this after removing empty segments, some keys need to be re mapped - for i in removed_seg[::-1]: - for j in range(i + 1, self._nb_segment): - # remap nsx seg index - for nsx_nb in self.nsx_to_load: - data = self.nsx_datas[nsx_nb].pop(j) - self.nsx_datas[nsx_nb][j - 1] = data + if all(keep_mask): + return - data_header = self._nsx_data_header[nsx_nb].pop(j) - self._nsx_data_header[nsx_nb][j - 1] = data_header + # Remap nev segment ids: shift down to fill gaps left by removed segments + if self._avail_files["nev"]: + for removed_index in [i for i, keep in enumerate(keep_mask) if not keep][::-1]: + for _, (_, ev_ids) in self.nev_data.items(): + ev_ids[ev_ids > removed_index] -= 1 - # Also remap nev data, ev_ids are the equivalent to keys above - if self._avail_files["nev"]: - for k, (data, ev_ids) in self.nev_data.items(): - ev_ids[ev_ids == j] -= 1 + # Filter segments + for nsx_nb in self.nsx_to_load: + self._segmented_data_headers[nsx_nb] = [ + seg for seg, keep in zip(self._segmented_data_headers[nsx_nb], keep_mask) if keep + ] - self._nb_segment -= 1 + self._nb_segment = sum(keep_mask) def _get_nonneural_evdicts_spec_v23(self, data): """ diff --git a/neo/test/iotest/test_blackrockio.py b/neo/test/iotest/test_blackrockio.py index 91fc2748b..3adccfd1d 100644 --- a/neo/test/iotest/test_blackrockio.py +++ b/neo/test/iotest/test_blackrockio.py @@ -289,12 +289,12 @@ def test_segment_detection_reset(self): # This fails, because in the nev there is no way to separate two segments with self.assertRaises(NeoReadWriteError): - reader = BlackrockIO(filename=filename, nsx_to_load=2, nev_override=filename_nev_fail) + reader = BlackrockIO(filename=filename, nsx_to_load=2, nev_override=filename_nev_fail, gap_tolerance_ms=0) # The correct file will issue a warning because a reset has occurred # and could be detected, but was not explicitly documented in the file with warnings.catch_warnings(record=True) as w: - reader = BlackrockIO(filename=filename, nsx_to_load=2) + reader = BlackrockIO(filename=filename, nsx_to_load=2, gap_tolerance_ms=0) self.assertGreaterEqual(len(w), 1) messages = [str(warning.message) for warning in w if warning.category == UserWarning] self.assertIn("Detected 1 undocumented segments within nev data after " "timestamps [5451].", messages) @@ -343,7 +343,7 @@ def test_segment_detection_pause(self): # And another one because there are spikes between segments with warnings.catch_warnings(record=True) as w: warnings.simplefilter("always") - reader = BlackrockIO(filename=filename, nsx_to_load=2, nev_override=filename_nev_outside_seg) + reader = BlackrockIO(filename=filename, nsx_to_load=2, nev_override=filename_nev_outside_seg, gap_tolerance_ms=0) self.assertGreaterEqual(len(w), 2) # Check that warnings are correct @@ -387,7 +387,7 @@ def test_segment_detection_pause(self): self.assertEqual(len(block.segments[1].analogsignals[0][:]), 4000) # This case is correct, no spikes outside segment or anything - reader = BlackrockIO(filename=filename, nsx_to_load=2) + reader = BlackrockIO(filename=filename, nsx_to_load=2, gap_tolerance_ms=0) block = reader.read_block(load_waveforms=False, signal_group_mode="split-all") # 2 segments diff --git a/neo/test/rawiotest/test_blackrockrawio.py b/neo/test/rawiotest/test_blackrockrawio.py index fa16c3e29..c8780b05e 100644 --- a/neo/test/rawiotest/test_blackrockrawio.py +++ b/neo/test/rawiotest/test_blackrockrawio.py @@ -267,20 +267,27 @@ def test_get_blackrock_timestamps_ptp(self): np.testing.assert_array_equal(timestamps_ns6[:5], expected_ns6) def test_get_blackrock_timestamps_ptp_with_gaps(self): - """Test _get_blackrock_timestamps for PTP format with gaps and multiple segments.""" + """Test _get_blackrock_timestamps for PTP format with gaps and multiple segments. + + This file has 3 timestamp discontinuities: + - Index 631: forward gap of ~0.97ms (dropped samples) + - Index 661: backward jump of ~1.9ms (time reversal) + - Index 689: forward gap of ~1.03ms (dropped samples) + Creating 4 segments: [0:632], [632:662], [662:690], [690:1000] + """ dirname = self.get_local_path("blackrock/blackrock_ptp_with_missing_samples/Hub1-NWBtestfile_neural_wspikes") gap_tolerance_ms = 0.5 reader = BlackrockRawIO(filename=dirname, nsx_to_load=6, gap_tolerance_ms=gap_tolerance_ms) reader.parse_header() n_segments = reader.segment_count(0) - self.assertEqual(n_segments, 3) + self.assertEqual(n_segments, 4) nanoseconds_per_second = 1_000_000_000.0 # This file has a single stream: ns6 (30 kHz, 1 channel) stream_index = 0 - expected_sizes = [632, 58, 310] + expected_sizes = [632, 30, 28, 310] # First 5 PTP clock values (nanoseconds since Unix epoch) per segment, read directly from file and hardcoded here for testing first_5_ptp_clock_ns_per_segment = [ @@ -288,12 +295,14 @@ def test_get_blackrock_timestamps_ptp_with_gaps(self): 1752531864717843077, 1752531864717876277], dtype="uint64"), np.array([1752531864739742985, 1752531864739776305, 1752531864739809665, 1752531864739842945, 1752531864739876385], dtype="uint64"), + np.array([1752531864738809624, 1752531864738842984, 1752531864738876304, + 1752531864738909664, 1752531864738942984], dtype="uint64"), np.array([1752531864740742986, 1752531864740776306, 1752531864740809626, 1752531864740842946, 1752531864740876346], dtype="uint64"), ] # Last PTP clock value (nanoseconds since Unix epoch) of each segment, hardcoded here for testing last_ptp_clock_ns_per_segment = np.array( - [1752531864738776304, 1752531864739709625, 1752531864751042999], dtype="uint64", + [1752531864738776304, 1752531864740709666, 1752531864739709625, 1752531864751042999], dtype="uint64", ) for seg_index in range(n_segments): @@ -308,13 +317,6 @@ def test_get_blackrock_timestamps_ptp_with_gaps(self): expected_last = last_ptp_clock_ns_per_segment[seg_index].astype("float64") / nanoseconds_per_second np.testing.assert_array_equal(timestamps[-1], expected_last) - # Verify the gaps between segments exceed the tolerance - for seg_index in range(n_segments - 1): - last_ts = last_ptp_clock_ns_per_segment[seg_index].astype("float64") / nanoseconds_per_second - first_ts_next = first_5_ptp_clock_ns_per_segment[seg_index + 1][0].astype("float64") / nanoseconds_per_second - gap_ms = (first_ts_next - last_ts) * 1000 - self.assertGreater(gap_ms, gap_tolerance_ms) - def test_gap_tolerance_ms_parameter(self): """ Test gap_tolerance_ms parameter for gap handling with files that have actual gaps. @@ -322,30 +324,60 @@ def test_gap_tolerance_ms_parameter(self): Tests the error-by-default behavior where files with timestamp gaps raise ValueError unless the user explicitly opts in with gap_tolerance_ms parameter. - See PR #1769 for the gap details on the example file used here. + This file has 2 forward gaps (~0.97ms and ~1.03ms) and 1 backward jump (~1.9ms). + Backward jumps always create segment boundaries regardless of tolerance. """ # Use stubbed files with missing samples (timestamp gaps) from SimulatedSpikes data dirname = self.get_local_path("blackrock/blackrock_ptp_with_missing_samples/Hub1-NWBtestfile_neural_wspikes") # Test 1: Default behavior (None) raises ValueError for files with gaps # This is the error-by-default behavior to ensure users are aware of data issues - with self.assertRaises(ValueError) as context: + with self.assertRaises(ValueError): reader = BlackrockRawIO(filename=dirname, nsx_to_load=6) reader.parse_header() - # Test 2: Explicit tolerance allows loading files with gaps - # User opts in by providing gap_tolerance_ms + # Test 2: Large tolerance filters forward gaps but backward jumps always split + # Forward gaps (~1ms) are under 10ms tolerance, but backward jump always creates boundary reader_with_tolerance = BlackrockRawIO(filename=dirname, nsx_to_load=6, gap_tolerance_ms=10.0) reader_with_tolerance.parse_header() segments_with_tolerance = reader_with_tolerance.segment_count(0) - self.assertEqual(1, segments_with_tolerance) # Gaps < 10ms are ignored + self.assertEqual(2, segments_with_tolerance) - # Test 3: Stricter tolerance creates 3 segments - # With strict tolerance (0.5ms), gaps > 0.5ms will create new segments + # Test 3: Strict tolerance creates 4 segments (2 forward gaps + 1 backward jump) reader_strict = BlackrockRawIO(filename=dirname, nsx_to_load=6, gap_tolerance_ms=0.5) reader_strict.parse_header() segments_strict = reader_strict.segment_count(0) - self.assertEqual(segments_strict, 3) # + self.assertEqual(segments_strict, 4) + + def test_standard_format_gap_error_by_default(self): + """Multi-block standard files raise ValueError when gap_tolerance_ms is None.""" + dirname = self.get_local_path("blackrock/segment/PauseCorrect/pause_correct") + with self.assertRaises(ValueError): + reader = BlackrockRawIO(filename=dirname, nsx_to_load=2) + reader.parse_header() + + def test_standard_format_gap_segmentation(self): + """Multi-block standard files create segments at block boundaries with gap_tolerance_ms=0.""" + dirname = self.get_local_path("blackrock/segment/PauseCorrect/pause_correct") + reader = BlackrockRawIO(filename=dirname, nsx_to_load=2, gap_tolerance_ms=0) + reader.parse_header() + self.assertEqual(reader.segment_count(0), 2) + + def test_standard_format_gap_merge(self): + """Multi-block standard files merge into 1 segment with very large gap_tolerance_ms.""" + dirname = self.get_local_path("blackrock/segment/PauseCorrect/pause_correct") + reader = BlackrockRawIO(filename=dirname, nsx_to_load=2, gap_tolerance_ms=100_000) + reader.parse_header() + self.assertEqual(reader.segment_count(0), 1) + + # Merged segment should have the sum of samples from both blocks + stream_index = 0 + total_size = reader.get_signal_size(0, 0, stream_index) + self.assertEqual(total_size, 8000) # 4000 + 4000 + + # Should be able to read the full signal + raw_sigs = reader.get_analogsignal_chunk(stream_index=stream_index) + self.assertEqual(raw_sigs.shape[0], 8000) if __name__ == "__main__": From 233de8556a3b8ddb3c7ee53e7554a7bd153b00bb Mon Sep 17 00:00:00 2001 From: Heberto Mayorquin Date: Mon, 16 Feb 2026 22:53:39 -0600 Subject: [PATCH 2/9] second refactor --- neo/rawio/blackrockrawio.py | 179 +++++++++++++++++++++--------------- 1 file changed, 104 insertions(+), 75 deletions(-) diff --git a/neo/rawio/blackrockrawio.py b/neo/rawio/blackrockrawio.py index 6ab438503..462f56645 100644 --- a/neo/rawio/blackrockrawio.py +++ b/neo/rawio/blackrockrawio.py @@ -440,27 +440,25 @@ def _parse_header(self): t_start = float("inf") t_stop = 0.0 for nsx_nb in self.nsx_to_load: - seg = self._segmented_data_headers[nsx_nb][data_bl] - t_start = min(t_start, seg["t_start"]) - nb_pts = seg["nb_data_points"] + segment_header = self._segmented_data_headers[nsx_nb][data_bl] + t_start = min(t_start, segment_header["t_start"]) + nb_pts = segment_header["nb_data_points"] sr = self._nsx_sampling_frequency[nsx_nb] - if self._nsx_data_spec == "3.0-ptp": + if "timestamps_memmap_kwargs" in segment_header: # PTP: read actual last timestamp (jitter makes t_start + n/sr imprecise) - channels = int(self._nsx_basic_header[nsx_nb]["channel_count"]) - filename = f"{self._filenames['nsx']}.ns{nsx_nb}" - ptp_dt = NSX_DATA_HEADER_TYPES["3.0-ptp"](channels) + ts_kw = segment_header["timestamps_memmap_kwargs"] memmap = np.memmap( - filename, dtype=ptp_dt, mode="r", - offset=seg["ptp_data_offset"], - shape=(seg["nb_data_points"],), + ts_kw["filename"], dtype=ts_kw["dtype"], mode="r", + offset=ts_kw["offset"], + shape=(ts_kw["num_samples"],), ) ts_res = float(self._nsx_basic_header[nsx_nb]["timestamp_resolution"]) - last_ts = float(memmap["timestamps"][-1]) / ts_res + last_ts = float(memmap[ts_kw["field"]][-1]) / ts_res del memmap seg_t_stop = last_ts + 1.0 / sr else: # Standard/v2.1: exact from t_start + nb_pts / sr - seg_t_stop = seg["t_start"] + nb_pts / sr + seg_t_stop = segment_header["t_start"] + nb_pts / sr t_stop = max(t_stop, seg_t_stop) if self._avail_files["nev"]: @@ -635,35 +633,37 @@ def _get_signal_t_start(self, block_index, seg_index, stream_index): def _get_analogsignal_chunk(self, block_index, seg_index, i_start, i_stop, stream_index, channel_indexes): stream_id = self.header["signal_streams"][stream_index]["id"] nsx_nb = int(stream_id) - seg = self._segmented_data_headers[nsx_nb][seg_index] + segment_header = self._segmented_data_headers[nsx_nb][seg_index] channels = int(self._nsx_basic_header[nsx_nb]["channel_count"]) - filename = f"{self._filenames['nsx']}.ns{nsx_nb}" if channel_indexes is None: channel_indexes = slice(None) - if self._nsx_data_spec == "3.0-ptp": - ptp_dt = NSX_DATA_HEADER_TYPES["3.0-ptp"](channels) - memmap = np.memmap( - filename, dtype=ptp_dt, mode="r", - offset=seg["ptp_data_offset"], - shape=(seg["nb_data_points"],), - ) - return memmap["samples"][i_start:i_stop, channel_indexes] - else: - specs = seg["memmap_specs"] - if len(specs) == 1: - data = np.memmap( - filename, dtype="int16", mode="r", - offset=specs[0]["offset"], - shape=(specs[0]["num_samples"], channels), + specs = segment_header["memmap_specs"] + if len(specs) == 1: + spec = specs[0] + field = spec.get("field") + if field is not None: + # PTP: structured dtype, extract named field (e.g. "samples") + memmap = np.memmap( + spec["filename"], dtype=spec["dtype"], mode="r", + offset=spec["offset"], + shape=(spec["num_samples"],), ) - return data[i_start:i_stop, channel_indexes] + return memmap[field][i_start:i_stop, channel_indexes] else: - return self._read_multi_block_chunk( - filename, specs, channels, i_start, i_stop, channel_indexes, + # Standard/v2.1: flat dtype (e.g. int16) + data = np.memmap( + spec["filename"], dtype=spec["dtype"], mode="r", + offset=spec["offset"], + shape=(spec["num_samples"], channels), ) + return data[i_start:i_stop, channel_indexes] + else: + return self._read_multi_block_chunk( + specs, channels, i_start, i_stop, channel_indexes, + ) - def _read_multi_block_chunk(self, filename, memmap_specs, channels, + def _read_multi_block_chunk(self, memmap_specs, channels, i_start, i_stop, channel_indexes): """ Read a chunk of analog signal data that spans multiple data blocks @@ -679,12 +679,10 @@ def _read_multi_block_chunk(self, filename, memmap_specs, channels, Parameters ---------- - filename : str - Path to the NSx file. memmap_specs : list[dict] - Each dict has "offset" (byte offset in file) and "num_samples" - (number of samples in that block). Ordered sequentially within - the segment. + Each dict has "filename" (path to file), "offset" (byte offset in + file) and "num_samples" (number of samples in that block). Ordered + sequentially within the segment. channels : int Number of channels (columns in the int16 data matrix). i_start : int or None @@ -715,7 +713,7 @@ def _read_multi_block_chunk(self, filename, memmap_specs, channels, local_start = max(0, i_start - block_start) local_stop = min(spec["num_samples"], i_stop - block_start) data = np.memmap( - filename, dtype="int16", mode="r", + spec["filename"], dtype=spec["dtype"], mode="r", offset=spec["offset"], shape=(spec["num_samples"], channels), ) @@ -768,26 +766,24 @@ def _get_blackrock_timestamps(self, block_index, seg_index, i_start, i_stop, str """ stream_id = self.header["signal_streams"][stream_index]["id"] nsx_nb = int(stream_id) - seg = self._segmented_data_headers[nsx_nb][seg_index] + segment_header = self._segmented_data_headers[nsx_nb][seg_index] - size = seg["nb_data_points"] + size = segment_header["nb_data_points"] i_start = i_start if i_start is not None else 0 i_stop = i_stop if i_stop is not None else size - if self._nsx_data_spec == "3.0-ptp": - channels = int(self._nsx_basic_header[nsx_nb]["channel_count"]) - filename = f"{self._filenames['nsx']}.ns{nsx_nb}" - ptp_dt = NSX_DATA_HEADER_TYPES["3.0-ptp"](channels) + if "timestamps_memmap_kwargs" in segment_header: + ts_kw = segment_header["timestamps_memmap_kwargs"] memmap = np.memmap( - filename, dtype=ptp_dt, mode="r", - offset=seg["ptp_data_offset"], - shape=(seg["nb_data_points"],), + ts_kw["filename"], dtype=ts_kw["dtype"], mode="r", + offset=ts_kw["offset"], + shape=(ts_kw["num_samples"],), ) ts_res = float(self._nsx_basic_header[nsx_nb]["timestamp_resolution"]) - return memmap["timestamps"][i_start:i_stop].astype("float64") / ts_res + return memmap[ts_kw["field"]][i_start:i_stop].astype("float64") / ts_res else: sr = self._nsx_sampling_frequency[nsx_nb] - return seg["t_start"] + np.arange(i_start, i_stop, dtype="float64") / sr + return segment_header["t_start"] + np.arange(i_start, i_stop, dtype="float64") / sr def _spike_count(self, block_index, seg_index, unit_index): channel_id, unit_id = self.internal_unit_ids[unit_index] @@ -1135,12 +1131,15 @@ def _parse_nsx_data(self, spec, nsx_nb): ------- list[dict] Each dict contains: - - "memmap_kwargs": {"offset": int, "num_samples": int} - Byte offset and sample count for creating memmaps at read time. - - "nsx_block_timestamp": scalar (standard only) + - "memmap_kwargs": {"filename": str, "dtype": str or np.dtype, + "offset": int, "num_samples": int} + Recipe for creating memmaps at read time. + - "nsx_block_timestamp": scalar (v2.2, v2.3, v3.0 only) Raw timestamp tick from the block header. - - "ptp_timestamps": np.ndarray of uint64 (PTP only) - Per-sample timestamps, temporary (used for gap detection then discarded). + - "ptp_timestamps_memmap_kwargs": {"filename": str, "dtype": np.dtype, + "offset": int, "num_samples": int} (v3.0-ptp only) + Recipe for creating a temporary memmap to read per-sample timestamps + during segmentation. Notes ----- @@ -1176,7 +1175,7 @@ def _parse_nsx_data_v21(self, nsx_nb): ) filesize = self._get_file_size(filename) num_samples = int((filesize - bytes_in_headers) / (2 * channels) - 1) - return [{"memmap_kwargs": {"offset": bytes_in_headers, "num_samples": num_samples}}] + return [{"memmap_kwargs": {"filename": filename, "dtype": "int16", "offset": bytes_in_headers, "num_samples": num_samples}}] def _parse_nsx_data_v22_v30(self, spec, nsx_nb): """ @@ -1213,7 +1212,7 @@ def _parse_nsx_data_v22_v30(self, spec, nsx_nb): data_offset = current_offset + header.dtype.itemsize parsed_data_headers.append({ - "memmap_kwargs": {"offset": data_offset, "num_samples": num_samples}, + "memmap_kwargs": {"filename": filename, "dtype": "int16", "offset": data_offset, "num_samples": num_samples}, "nsx_block_timestamp": header["timestamp"], }) @@ -1235,7 +1234,8 @@ def _parse_nsx_data_v30_ptp(self, nsx_nb): ------- list[dict] [{"memmap_kwargs": {"offset": int, "num_samples": int}, - "ptp_timestamps": np.ndarray (uint64)}] + "ptp_timestamps_memmap_kwargs": {"offset": int, "dtype": np.dtype, + "num_samples": int}}] """ filename = f"{self._filenames['nsx']}.ns{nsx_nb}" @@ -1246,7 +1246,7 @@ def _parse_nsx_data_v30_ptp(self, nsx_nb): ptp_dt = NSX_DATA_HEADER_TYPES["3.0-ptp"](channel_count) npackets = int((filesize - header_size) / np.dtype(ptp_dt).itemsize) - # Temporary memmap for verification only + # Temporary memmap to verify this is truly PTP (every packet has exactly 1 sample) temp_memmap = np.memmap(filename, dtype=ptp_dt, shape=npackets, offset=header_size, mode="r") if not np.all(temp_memmap["num_data_points"] == 1): @@ -1254,13 +1254,16 @@ def _parse_nsx_data_v30_ptp(self, nsx_nb): del temp_memmap return self._parse_nsx_data_v22_v30("3.0", nsx_nb) - # Copy timestamps for gap detection (released after segmentation) - ptp_timestamps = np.array(temp_memmap["timestamps"]) del temp_memmap return [{ - "memmap_kwargs": {"offset": header_size, "num_samples": npackets}, - "ptp_timestamps": ptp_timestamps, + "memmap_kwargs": {"filename": filename, "offset": header_size, "num_samples": npackets}, + "ptp_timestamps_memmap_kwargs": { + "filename": filename, + "offset": header_size, + "dtype": ptp_dt, + "num_samples": npackets, + }, }] def _format_gap_report(self, gap_indices, timestamps_in_seconds, time_differences, nsx_nb): @@ -1380,12 +1383,14 @@ def _segment_nsx_data(self, parsed_data_headers, nsx_nb, gap_tolerance_ms): ------- list[dict] Each dict contains: + - "memmap_specs": list of {"filename": str, "dtype": str or np.dtype, + "offset": int, "num_samples": int}, with optional "field": str + for structured dtypes (all formats) - "nb_data_points": int - "t_start": float (seconds) - - "memmap_specs": list of {"offset": int, "num_samples": int} - (standard/v2.1 only) - - "ptp_data_offset": int (PTP only) - - "nsx_block_timestamp": scalar raw ticks (standard only) + - "timestamps_memmap_kwargs": {"filename": str, "dtype": np.dtype, + "offset": int, "num_samples": int, "field": str} (v3.0-ptp only) + - "nsx_block_timestamp": scalar raw ticks (v2.2, v2.3, v3.0 only) """ data_spec = self._nsx_data_spec if data_spec == "2.1": @@ -1472,21 +1477,30 @@ def _segment_nsx_ptp(self, parsed_data_headers, nsx_nb, gap_tolerance_ms): """ Segment PTP format data (v3.0-ptp) using sample-level gap detection. - Operates on per-sample uint64 timestamps. The ptp_timestamps array from - the parser is consumed here and not stored on segments. + Creates a temporary memmap from ptp_timestamps_memmap_kwargs to read + per-sample uint64 timestamps for gap detection. The memmap is released + after segmentation; only byte offsets are stored on segments. """ ts_res = float(self._nsx_basic_header[nsx_nb]["timestamp_resolution"]) sr = self._nsx_sampling_frequency[nsx_nb] - channel_count = int(self._nsx_basic_header[nsx_nb]["channel_count"]) - ptp_dt = NSX_DATA_HEADER_TYPES["3.0-ptp"](channel_count) - packet_size = np.dtype(ptp_dt).itemsize header = parsed_data_headers[0] base_offset = header["memmap_kwargs"]["offset"] num_samples = header["memmap_kwargs"]["num_samples"] - ptp_timestamps = header["ptp_timestamps"] # temporary, consumed here + filename = header["memmap_kwargs"]["filename"] + + # Create temporary memmap to read timestamps for gap detection + ts_kwargs = header["ptp_timestamps_memmap_kwargs"] + temp_memmap = np.memmap( + ts_kwargs["filename"], dtype=ts_kwargs["dtype"], mode="r", + offset=ts_kwargs["offset"], shape=(ts_kwargs["num_samples"],), + ) + ptp_timestamps = temp_memmap["timestamps"] + ptp_dt = ts_kwargs["dtype"] + packet_size = np.dtype(ptp_dt).itemsize all_timestamps = ptp_timestamps.astype("float64") / ts_res + del temp_memmap # release file handle; all_timestamps is an independent copy time_diffs = np.diff(all_timestamps) deviations = time_diffs - 1.0 / sr @@ -1502,9 +1516,24 @@ def _segment_nsx_ptp(self, parsed_data_headers, nsx_nb, gap_tolerance_ms): for seg_index in range(len(segment_starts)): start = int(segment_starts[seg_index]) end = int(segment_ends[seg_index]) + seg_num_samples = end - start + seg_offset = base_offset + start * packet_size segments.append({ - "ptp_data_offset": base_offset + start * packet_size, - "nb_data_points": end - start, + "memmap_specs": [{ + "filename": filename, + "offset": seg_offset, + "num_samples": seg_num_samples, + "dtype": ptp_dt, + "field": "samples", + }], + "timestamps_memmap_kwargs": { + "filename": filename, + "offset": seg_offset, + "num_samples": seg_num_samples, + "dtype": ptp_dt, + "field": "timestamps", + }, + "nb_data_points": seg_num_samples, "t_start": float(all_timestamps[start]), }) From 50929c0c93b4cb5deecab4a97f3ab6d420de9d9b Mon Sep 17 00:00:00 2001 From: Heberto Mayorquin Date: Mon, 16 Feb 2026 23:17:29 -0600 Subject: [PATCH 3/9] Use buffer API pattern for blackrock --- neo/rawio/blackrockrawio.py | 253 ++++++++++++++++++++++++++---------- 1 file changed, 187 insertions(+), 66 deletions(-) diff --git a/neo/rawio/blackrockrawio.py b/neo/rawio/blackrockrawio.py index 462f56645..2239d1e7e 100644 --- a/neo/rawio/blackrockrawio.py +++ b/neo/rawio/blackrockrawio.py @@ -61,6 +61,7 @@ """ import datetime +import mmap import os import re import warnings @@ -447,14 +448,18 @@ def _parse_header(self): if "timestamps_memmap_kwargs" in segment_header: # PTP: read actual last timestamp (jitter makes t_start + n/sr imprecise) ts_kw = segment_header["timestamps_memmap_kwargs"] - memmap = np.memmap( - ts_kw["filename"], dtype=ts_kw["dtype"], mode="r", + fid = self._get_nsx_fid(nsx_nb) + timestamps = self._create_mmap_view( + fid=fid, + dtype=ts_kw["dtype"], offset=ts_kw["offset"], - shape=(ts_kw["num_samples"],), + num_samples=ts_kw["num_samples"], + packet_size=ts_kw.get("packet_size"), + item_offset=ts_kw.get("item_offset", 0), ) ts_res = float(self._nsx_basic_header[nsx_nb]["timestamp_resolution"]) - last_ts = float(memmap[ts_kw["field"]][-1]) / ts_res - del memmap + last_ts = float(timestamps[-1]) / ts_res + del timestamps seg_t_stop = last_ts + 1.0 / sr else: # Standard/v2.1: exact from t_start + nb_pts / sr @@ -630,6 +635,93 @@ def _get_signal_t_start(self, block_index, seg_index, stream_index): nsx_nb = int(stream_id) return self._segmented_data_headers[nsx_nb][seg_index]["t_start"] + @staticmethod + def _create_mmap_view(fid, dtype, offset, num_samples, num_channels=None, + packet_size=None, item_offset=0): + """ + Create an np.ndarray view over a raw mmap buffer from an open file. + + When packet_size is None, creates a standard contiguous memmap (for + standard/v2.1 formats where samples are stored contiguously as int16). + + When packet_size is provided, creates a strided view that extracts + interleaved fields from PTP packets. Each packet has a fixed size, + and the target field starts at item_offset bytes into the packet. + The stride between rows equals packet_size, allowing the view to + skip over other fields (timestamps, reserved bytes, etc.). + + Parameters + ---------- + fid : file-like + Open file object (must support .fileno()). + dtype : str or np.dtype + Data type of the target field (e.g. "int16" for samples, + "uint64" for timestamps). + offset : int + Byte offset in the file where the data region starts. + num_samples : int + Number of samples (rows) to read. + num_channels : int or None + Number of channels (columns). None for 1D arrays (e.g. timestamps). + packet_size : int or None + Stride between consecutive rows in bytes. None for contiguous data. + item_offset : int + Byte offset of the target field within each packet. + + Returns + ------- + np.ndarray + View into the memory-mapped file. Shape is (num_samples, num_channels) + when num_channels is provided, or (num_samples,) otherwise. + """ + dtype = np.dtype(dtype) + + if num_channels is not None: + shape = (num_samples, num_channels) + else: + shape = (num_samples,) + + if packet_size is None: + # Contiguous data (standard/v2.1) + bytes_per_sample = dtype.itemsize * (num_channels if num_channels is not None else 1) + start_byte = offset + length = num_samples * bytes_per_sample + else: + # Strided access (PTP) + start_byte = offset + item_offset + length = (num_samples - 1) * packet_size + ( + dtype.itemsize * num_channels if num_channels is not None else dtype.itemsize + ) + + # mmap offset must be aligned to ALLOCATIONGRANULARITY + mmap_offset, start_remainder = divmod(start_byte, mmap.ALLOCATIONGRANULARITY) + mmap_offset *= mmap.ALLOCATIONGRANULARITY + length += start_remainder + + raw_mmap = mmap.mmap(fid.fileno(), length=length, access=mmap.ACCESS_READ, offset=mmap_offset) + + if packet_size is not None: + strides = (packet_size, dtype.itemsize) if num_channels is not None else (packet_size,) + else: + strides = None # default contiguous strides + + return np.ndarray( + shape=shape, + dtype=dtype, + buffer=raw_mmap, + offset=start_remainder, + strides=strides, + ) + + def _get_nsx_fid(self, nsx_nb): + """Open NSX file on demand and cache the file descriptor for reuse.""" + if not hasattr(self, "_nsx_fids"): + self._nsx_fids = {} + if nsx_nb not in self._nsx_fids: + filename = f"{self._filenames['nsx']}.ns{nsx_nb}" + self._nsx_fids[nsx_nb] = open(filename, "rb") + return self._nsx_fids[nsx_nb] + def _get_analogsignal_chunk(self, block_index, seg_index, i_start, i_stop, stream_index, channel_indexes): stream_id = self.header["signal_streams"][stream_index]["id"] nsx_nb = int(stream_id) @@ -639,31 +731,25 @@ def _get_analogsignal_chunk(self, block_index, seg_index, i_start, i_stop, strea channel_indexes = slice(None) specs = segment_header["memmap_specs"] + fid = self._get_nsx_fid(nsx_nb) if len(specs) == 1: spec = specs[0] - field = spec.get("field") - if field is not None: - # PTP: structured dtype, extract named field (e.g. "samples") - memmap = np.memmap( - spec["filename"], dtype=spec["dtype"], mode="r", - offset=spec["offset"], - shape=(spec["num_samples"],), - ) - return memmap[field][i_start:i_stop, channel_indexes] - else: - # Standard/v2.1: flat dtype (e.g. int16) - data = np.memmap( - spec["filename"], dtype=spec["dtype"], mode="r", - offset=spec["offset"], - shape=(spec["num_samples"], channels), - ) - return data[i_start:i_stop, channel_indexes] + data = self._create_mmap_view( + fid=fid, + dtype=spec["dtype"], + offset=spec["offset"], + num_samples=spec["num_samples"], + num_channels=spec.get("num_channels", channels), + packet_size=spec.get("packet_size"), + item_offset=spec.get("item_offset", 0), + ) + return data[i_start:i_stop, channel_indexes] else: return self._read_multi_block_chunk( - specs, channels, i_start, i_stop, channel_indexes, + fid, specs, channels, i_start, i_stop, channel_indexes, ) - def _read_multi_block_chunk(self, memmap_specs, channels, + def _read_multi_block_chunk(self, fid, memmap_specs, channels, i_start, i_stop, channel_indexes): """ Read a chunk of analog signal data that spans multiple data blocks @@ -673,15 +759,17 @@ def _read_multi_block_chunk(self, memmap_specs, channels, without significant gaps are merged into one segment by the segmenter. Each block is stored at a different file offset, so reading a contiguous sample range [i_start, i_stop) may require stitching data from several - blocks. This method creates temporary memmaps only for the blocks that + blocks. This method creates mmap views only for the blocks that overlap the requested range, slices each one, and concatenates the results. Parameters ---------- + fid : file-like + Open file object (must support .fileno()). memmap_specs : list[dict] - Each dict has "filename" (path to file), "offset" (byte offset in - file) and "num_samples" (number of samples in that block). Ordered + Each dict has "offset" (byte offset in file), "dtype" (data type), + and "num_samples" (number of samples in that block). Ordered sequentially within the segment. channels : int Number of channels (columns in the int16 data matrix). @@ -696,7 +784,7 @@ def _read_multi_block_chunk(self, memmap_specs, channels, ------- np.ndarray Signal data of shape (i_stop - i_start, len(channel_indexes)), - dtype int16. Returns a memmap view when only one block is touched, + dtype int16. Returns a mmap view when only one block is touched, otherwise a copied concatenation. """ total_samples = sum(spec["num_samples"] for spec in memmap_specs) @@ -712,10 +800,12 @@ def _read_multi_block_chunk(self, memmap_specs, channels, if block_end > i_start and block_start < i_stop: local_start = max(0, i_start - block_start) local_stop = min(spec["num_samples"], i_stop - block_start) - data = np.memmap( - spec["filename"], dtype=spec["dtype"], mode="r", + data = self._create_mmap_view( + fid=fid, + dtype=spec["dtype"], offset=spec["offset"], - shape=(spec["num_samples"], channels), + num_samples=spec["num_samples"], + num_channels=channels, ) pieces.append(data[local_start:local_stop, channel_indexes]) cumulative = block_end @@ -774,13 +864,17 @@ def _get_blackrock_timestamps(self, block_index, seg_index, i_start, i_stop, str if "timestamps_memmap_kwargs" in segment_header: ts_kw = segment_header["timestamps_memmap_kwargs"] - memmap = np.memmap( - ts_kw["filename"], dtype=ts_kw["dtype"], mode="r", + fid = self._get_nsx_fid(nsx_nb) + timestamps = self._create_mmap_view( + fid=fid, + dtype=ts_kw["dtype"], offset=ts_kw["offset"], - shape=(ts_kw["num_samples"],), + num_samples=ts_kw["num_samples"], + packet_size=ts_kw.get("packet_size"), + item_offset=ts_kw.get("item_offset", 0), ) ts_res = float(self._nsx_basic_header[nsx_nb]["timestamp_resolution"]) - return memmap[ts_kw["field"]][i_start:i_stop].astype("float64") / ts_res + return timestamps[i_start:i_stop].astype("float64") / ts_res else: sr = self._nsx_sampling_frequency[nsx_nb] return segment_header["t_start"] + np.arange(i_start, i_stop, dtype="float64") / sr @@ -1233,9 +1327,12 @@ def _parse_nsx_data_v30_ptp(self, nsx_nb): Returns ------- list[dict] - [{"memmap_kwargs": {"offset": int, "num_samples": int}, - "ptp_timestamps_memmap_kwargs": {"offset": int, "dtype": np.dtype, - "num_samples": int}}] + [{"memmap_kwargs": {strided spec for samples}, + "ptp_timestamps_memmap_kwargs": {strided spec for timestamps}}] + + Strided specs contain "packet_size" and "item_offset" which allow + consumers to create np.ndarray views with custom strides over a raw + mmap buffer, extracting interleaved fields without structured dtypes. """ filename = f"{self._filenames['nsx']}.ns{nsx_nb}" @@ -1244,7 +1341,8 @@ def _parse_nsx_data_v30_ptp(self, nsx_nb): channel_count = int(self._nsx_basic_header[nsx_nb]["channel_count"]) ptp_dt = NSX_DATA_HEADER_TYPES["3.0-ptp"](channel_count) - npackets = int((filesize - header_size) / np.dtype(ptp_dt).itemsize) + ptp_dtype = np.dtype(ptp_dt) + npackets = int((filesize - header_size) / ptp_dtype.itemsize) # Temporary memmap to verify this is truly PTP (every packet has exactly 1 sample) temp_memmap = np.memmap(filename, dtype=ptp_dt, shape=npackets, offset=header_size, mode="r") @@ -1256,13 +1354,27 @@ def _parse_nsx_data_v30_ptp(self, nsx_nb): del temp_memmap + packet_size = ptp_dtype.itemsize + samples_item_offset = ptp_dtype.fields["samples"][1] + timestamps_item_offset = ptp_dtype.fields["timestamps"][1] + return [{ - "memmap_kwargs": {"filename": filename, "offset": header_size, "num_samples": npackets}, + "memmap_kwargs": { + "filename": filename, + "dtype": "int16", + "offset": header_size, + "num_samples": npackets, + "num_channels": channel_count, + "packet_size": packet_size, + "item_offset": samples_item_offset, + }, "ptp_timestamps_memmap_kwargs": { "filename": filename, + "dtype": "uint64", "offset": header_size, - "dtype": ptp_dt, "num_samples": npackets, + "packet_size": packet_size, + "item_offset": timestamps_item_offset, }, }] @@ -1383,13 +1495,15 @@ def _segment_nsx_data(self, parsed_data_headers, nsx_nb, gap_tolerance_ms): ------- list[dict] Each dict contains: - - "memmap_specs": list of {"filename": str, "dtype": str or np.dtype, - "offset": int, "num_samples": int}, with optional "field": str - for structured dtypes (all formats) + - "memmap_specs": list of {"filename": str, "dtype": str, + "offset": int, "num_samples": int}. PTP specs additionally + include "num_channels", "packet_size", and "item_offset" for + strided access (all formats) - "nb_data_points": int - "t_start": float (seconds) - - "timestamps_memmap_kwargs": {"filename": str, "dtype": np.dtype, - "offset": int, "num_samples": int, "field": str} (v3.0-ptp only) + - "timestamps_memmap_kwargs": {"filename": str, "dtype": str, + "offset": int, "num_samples": int, "packet_size": int, + "item_offset": int} (v3.0-ptp only) - "nsx_block_timestamp": scalar raw ticks (v2.2, v2.3, v3.0 only) """ data_spec = self._nsx_data_spec @@ -1477,30 +1591,34 @@ def _segment_nsx_ptp(self, parsed_data_headers, nsx_nb, gap_tolerance_ms): """ Segment PTP format data (v3.0-ptp) using sample-level gap detection. - Creates a temporary memmap from ptp_timestamps_memmap_kwargs to read - per-sample uint64 timestamps for gap detection. The memmap is released - after segmentation; only byte offsets are stored on segments. + Creates a temporary strided view over the raw file to read per-sample + uint64 timestamps for gap detection. The view is released after + segmentation; only strided memmap specs are stored on segments. """ ts_res = float(self._nsx_basic_header[nsx_nb]["timestamp_resolution"]) sr = self._nsx_sampling_frequency[nsx_nb] header = parsed_data_headers[0] - base_offset = header["memmap_kwargs"]["offset"] - num_samples = header["memmap_kwargs"]["num_samples"] - filename = header["memmap_kwargs"]["filename"] - - # Create temporary memmap to read timestamps for gap detection + samples_kwargs = header["memmap_kwargs"] ts_kwargs = header["ptp_timestamps_memmap_kwargs"] - temp_memmap = np.memmap( - ts_kwargs["filename"], dtype=ts_kwargs["dtype"], mode="r", - offset=ts_kwargs["offset"], shape=(ts_kwargs["num_samples"],), - ) - ptp_timestamps = temp_memmap["timestamps"] - ptp_dt = ts_kwargs["dtype"] - packet_size = np.dtype(ptp_dt).itemsize + base_offset = samples_kwargs["offset"] + num_samples = samples_kwargs["num_samples"] + filename = samples_kwargs["filename"] + packet_size = samples_kwargs["packet_size"] + + # Create temporary strided view to read timestamps for gap detection + fid = self._get_nsx_fid(nsx_nb) + ptp_timestamps = self._create_mmap_view( + fid=fid, + dtype=ts_kwargs["dtype"], + offset=base_offset, + num_samples=num_samples, + packet_size=packet_size, + item_offset=ts_kwargs["item_offset"], + ) all_timestamps = ptp_timestamps.astype("float64") / ts_res - del temp_memmap # release file handle; all_timestamps is an independent copy + del ptp_timestamps # release mmap buffer; all_timestamps is an independent copy time_diffs = np.diff(all_timestamps) deviations = time_diffs - 1.0 / sr @@ -1521,17 +1639,20 @@ def _segment_nsx_ptp(self, parsed_data_headers, nsx_nb, gap_tolerance_ms): segments.append({ "memmap_specs": [{ "filename": filename, + "dtype": samples_kwargs["dtype"], "offset": seg_offset, "num_samples": seg_num_samples, - "dtype": ptp_dt, - "field": "samples", + "num_channels": samples_kwargs["num_channels"], + "packet_size": packet_size, + "item_offset": samples_kwargs["item_offset"], }], "timestamps_memmap_kwargs": { "filename": filename, + "dtype": ts_kwargs["dtype"], "offset": seg_offset, "num_samples": seg_num_samples, - "dtype": ptp_dt, - "field": "timestamps", + "packet_size": packet_size, + "item_offset": ts_kwargs["item_offset"], }, "nb_data_points": seg_num_samples, "t_start": float(all_timestamps[start]), From a689ad25f9a146f7d44a48b22e080d888c7eb0fb Mon Sep 17 00:00:00 2001 From: Heberto Mayorquin Date: Thu, 16 Jul 2026 12:18:21 -0600 Subject: [PATCH 4/9] Extend gap_tolerance_ms to standard nsx formats (v2.2, v2.3, v3.0) Master detects gaps only for PTP files, where every sample carries its own timestamp. Standard-format files were segmented by making each data block a segment, with no gap detection and no way to merge blocks that a pause never really separated. Split the segmentation into a dispatcher plus one segmenter per format, and give the standard path block-level gap detection through the same _classify_gaps step the PTP path uses. Both hand it the interval between consecutive samples, so the two paths share one threshold and one report. Merging blocks means a segment can now span several of them, at different file offsets, so segments carry a list of memmap specs rather than a single one and _read_multi_block_chunk stitches reads that straddle a boundary. Backward jumps now split as well. A clock reset makes the timestamps run backwards rather than leaving a gap, so there is no duration to compare against the tolerance and the blocks either side must never merge. --- neo/rawio/blackrockrawio.py | 375 +++++++++++++++------- neo/test/iotest/test_blackrockio.py | 4 +- neo/test/rawiotest/test_blackrockrawio.py | 124 +++++-- 3 files changed, 372 insertions(+), 131 deletions(-) diff --git a/neo/rawio/blackrockrawio.py b/neo/rawio/blackrockrawio.py index ccf40bb8c..f4517de4e 100644 --- a/neo/rawio/blackrockrawio.py +++ b/neo/rawio/blackrockrawio.py @@ -124,6 +124,10 @@ class BlackrockRawIO(BaseRawIO): - PTP format (v3.0-ptp): Gaps in per-sample timestamps - Standard format (v2.2/2.3/3.0): Gaps between data blocks + Backward jumps, where the clock runs backwards from one sample to the next, + are not gaps and have no duration to compare against this tolerance. They + always create a segment boundary, whatever the value given here. + Notes ----- * Note: This routine will handle files according to specification 2.1, 2.2, @@ -739,8 +743,15 @@ def _get_analogsignal_chunk(self, block_index, seg_index, i_start, i_stop, strea nsx_nb = int(stream_id) seg = self._nsx_data_header[nsx_nb][seg_index] fid = self._get_nsx_fid(nsx_nb) - memmap_kwargs = seg["memmap_kwargs"] + memmap_specs = seg["memmap_specs"] channels = int(self._nsx_basic_header[nsx_nb]["channel_count"]) + if channel_indexes is None: + channel_indexes = slice(None) + + if len(memmap_specs) > 1: + return self._read_multi_block_chunk(fid, memmap_specs, channels, i_start, i_stop, channel_indexes) + + memmap_kwargs = memmap_specs[0] data = self._create_mmap_view( fid=fid, dtype=memmap_kwargs["dtype"], @@ -750,10 +761,66 @@ def _get_analogsignal_chunk(self, block_index, seg_index, i_start, i_stop, strea packet_size=memmap_kwargs.get("packet_size"), item_offset=memmap_kwargs.get("item_offset", 0), ) - if channel_indexes is None: - channel_indexes = slice(None) return data[i_start:i_stop, channel_indexes] + def _read_multi_block_chunk(self, fid, memmap_specs, channels, i_start, i_stop, channel_indexes): + """ + Read a chunk from a segment that merges several data blocks. + + In the standard format, consecutive blocks that are not separated by a + significant gap are merged into one segment. Each block lives at its own file + offset, so a contiguous sample range can straddle several of them. Only the + blocks overlapping the requested range are mapped, sliced, and concatenated. + + Parameters + ---------- + fid : file-like + Open file object (must support .fileno()). + memmap_specs : list[dict] + One spec per block in the segment, in acquisition order. + channels : int + Number of channels (columns) in the int16 data matrix. + i_start : int or None + First sample index, relative to the segment. None means 0. + i_stop : int or None + Stop sample index (exclusive), relative to the segment. None means the end. + channel_indexes : slice or array-like + Which channels to return. + + Returns + ------- + np.ndarray + Signal data of shape (i_stop - i_start, len(channel_indexes)), dtype int16. + """ + total_samples = sum(spec["num_samples"] for spec in memmap_specs) + if i_start is None: + i_start = 0 + if i_stop is None: + i_stop = total_samples + + pieces = [] + block_start = 0 + for spec in memmap_specs: + block_stop = block_start + spec["num_samples"] + if block_stop > i_start and block_start < i_stop: + data = self._create_mmap_view( + fid=fid, + dtype=spec["dtype"], + offset=spec["offset"], + num_samples=spec["num_samples"], + num_channels=spec.get("num_channels", channels), + packet_size=spec.get("packet_size"), + item_offset=spec.get("item_offset", 0), + ) + local_start = max(0, i_start - block_start) + local_stop = min(spec["num_samples"], i_stop - block_start) + pieces.append(data[local_start:local_stop, channel_indexes]) + block_start = block_stop + + if len(pieces) == 1: + return pieces[0] + return np.concatenate(pieces, axis=0) + def _get_blackrock_timestamps(self, block_index, seg_index, i_start, i_stop, stream_index): """ Return timestamps in seconds for analog signal samples in the given range. @@ -1397,12 +1464,11 @@ def _segment_nsx_data(self, parsed_data_headers, nsx_nb): """ Segment NSX data based on timestamp gaps. - Takes the parsed data headers returned by _parse_nsx_data() and creates segments. - Segmentation logic depends on the file format: - - - Standard format (multiple blocks): Each block IS a segment - - PTP format (single block with timestamp array): Detect gaps in timestamps - - V2.1 format (no timestamps): Single segment + Takes the parsed data headers returned by _parse_nsx_data() and dispatches to + the segmenter for the file's format. The dispatch keys off the structure of the + parsed headers rather than the declared spec version, because + _parse_nsx_data_v30_ptp() falls back to standard parsing for files that declare + PTP but do not use it. Parameters ---------- @@ -1418,124 +1484,213 @@ def _segment_nsx_data(self, parsed_data_headers, nsx_nb): [ { "timestamp": scalar or None, + Raw timestamp ticks at the start of the segment. "nb_data_points": int, "header": int or None, "offset_to_data_block": None (deprecated but kept for compatibility), - "memmap_kwargs": dict, + "memmap_specs": list[dict], + One entry per data block making up the segment. Standard-format + segments hold several when consecutive blocks were merged; PTP + and v2.1 segments always hold exactly one. "timestamps_memmap_kwargs": dict (PTP only), }, ... ] """ - segments = [] + first_header = parsed_data_headers[0] + + if "ptp_timestamps_memmap_kwargs" in first_header: + return self._segment_nsx_ptp(parsed_data_headers, nsx_nb) + elif first_header["timestamp"] is None: + return self._segment_nsx_v21(parsed_data_headers) + else: + return self._segment_nsx_v22_v30(parsed_data_headers, nsx_nb) + + def _classify_gaps(self, time_differences, sampling_rate, nsx_nb, timestamps_in_seconds): + """ + Classify timestamp discontinuities into segment boundaries. + + Shared by the standard and PTP segmenters. Both hand over the same quantity: + the interval between the timestamps of two consecutive samples, which is one + sampling period when the data is contiguous. + + Two kinds of discontinuity are treated differently: + + - Forward gaps (a pause, or dropped samples) are detected with a fixed + threshold of 2 sampling periods, which sits above the jitter and rounding + noise floor of these files. Whether a detected forward gap actually becomes + a segment boundary is then up to gap_tolerance_ms. + - Backward jumps (the clock runs backwards between consecutive samples) are an + abnormal condition rather than missing data. They have no gap duration to + compare against a tolerance, so they always become segment boundaries. + + Parameters + ---------- + time_differences : np.ndarray + Intervals in seconds between the timestamps of consecutive samples. + sampling_rate : float + Sampling rate in Hz. + nsx_nb : int + NSX file number (for error reporting). + timestamps_in_seconds : np.ndarray + Timestamps in seconds, used to report where each gap falls. + + Returns + ------- + np.ndarray + Indices into time_differences at which a segment boundary falls; a boundary + at index i means a new segment starts at sample i + 1. + """ + # Detection threshold: use strict 2x sampling period to find ALL gaps + detection_threshold = 2.0 / sampling_rate + + forward_mask = time_differences > detection_threshold + backward_mask = time_differences < 0 + gap_indices = np.flatnonzero(forward_mask | backward_mask) + + if len(gap_indices) == 0: + return gap_indices + + gap_report = self._format_gap_report(gap_indices, timestamps_in_seconds, time_differences, nsx_nb) + + # Error by default - user must opt-in to segmentation + if self.gap_tolerance_ms is None: + raise ValueError( + f"Detected {len(gap_indices)} timestamp gaps in ns{nsx_nb} file.\n" + f"{gap_report}\n" + f"To load this data, provide gap_tolerance_ms parameter to automatically " + f"segment at gaps larger than the specified tolerance." + ) + + # User provided tolerance - keep forward gaps above it, and every backward jump + # (no warning - user opted in) + gap_tolerance_s = self.gap_tolerance_ms / 1000.0 + significant_forward_mask = forward_mask & (time_differences > gap_tolerance_s) + return np.flatnonzero(significant_forward_mask | backward_mask) + + def _segment_nsx_v21(self, parsed_data_headers): + """ + Segment v2.1 data. + + The format has a single continuous block and no timestamps, so there is nothing + to detect gaps in and the result is always one segment. + """ + block_info = parsed_data_headers[0] + return [ + { + "timestamp": block_info["timestamp"], + "nb_data_points": block_info["memmap_kwargs"]["num_samples"], + "header": None, + "offset_to_data_block": None, + "memmap_specs": [block_info["memmap_kwargs"]], + } + ] + + def _segment_nsx_v22_v30(self, parsed_data_headers, nsx_nb): + """ + Segment standard format data (v2.2, v2.3, v3.0) using block-level gap detection. + + The file already breaks data into a new block at every pause, so block + boundaries are the candidate segment boundaries and gap detection decides which + of them are real. Blocks separated by less than gap_tolerance_ms are merged into + one segment, which then keeps one memmap spec per contributing block because + those blocks sit at different offsets in the file. + """ + ts_res = float(self._nsx_basic_header[nsx_nb]["timestamp_resolution"]) + sampling_rate = self._nsx_sampling_frequency[nsx_nb] + + block_t_starts = np.array([float(header["timestamp"]) / ts_res for header in parsed_data_headers]) + block_sizes = np.array([header["memmap_kwargs"]["num_samples"] for header in parsed_data_headers]) - # Case 1: Multiple blocks (Standard format) - each block is a segment if len(parsed_data_headers) > 1: - for block_info in parsed_data_headers: - segments.append({ - "timestamp": block_info["timestamp"], - "nb_data_points": block_info["memmap_kwargs"]["num_samples"], - "header": 1, # Standard format has headers - "offset_to_data_block": None, # Not needed - "memmap_kwargs": block_info["memmap_kwargs"], - }) - - # Case 2: Single block - check if PTP (has ptp_timestamps_memmap_kwargs) or simple - elif len(parsed_data_headers) == 1: - block_info = parsed_data_headers[0] - - # PTP format: read timestamps on demand and detect gaps - if "ptp_timestamps_memmap_kwargs" in block_info: - ts_kw = block_info["ptp_timestamps_memmap_kwargs"] - samples_kw = block_info["memmap_kwargs"] - - # Read timestamps via strided mmap view for gap detection - fid = self._get_nsx_fid(nsx_nb) - timestamps = self._create_mmap_view( - fid=fid, - dtype=ts_kw["dtype"], - offset=ts_kw["offset"], - num_samples=ts_kw["num_samples"], - packet_size=ts_kw.get("packet_size"), - item_offset=ts_kw.get("item_offset", 0), - ) + # Interval between the last sample of a block and the first sample of the + # next one, so contiguous blocks give exactly one sampling period. This is + # the per-block equivalent of np.diff() over the PTP per-sample timestamps. + block_last_sample_times = block_t_starts + (block_sizes - 1) / sampling_rate + inter_block_intervals = block_t_starts[1:] - block_last_sample_times[:-1] + boundary_indices = self._classify_gaps(inter_block_intervals, sampling_rate, nsx_nb, block_t_starts) + else: + boundary_indices = np.array([], dtype=np.intp) - # Analyze timestamp gaps - sampling_rate = self._nsx_sampling_frequency[nsx_nb] + # Runs of blocks between boundaries become one segment each + segment_starts = np.hstack((0, boundary_indices + 1)).astype(np.intp) + segment_stops = np.hstack((boundary_indices + 1, len(parsed_data_headers))).astype(np.intp) - # Detection threshold: use strict 2x sampling period to find ALL gaps - detection_threshold = 2.0 / sampling_rate + segments = [] + for start, stop in zip(segment_starts, segment_stops): + block_infos = parsed_data_headers[start:stop] + segments.append({ + "timestamp": block_infos[0]["timestamp"], + "nb_data_points": int(sum(block["memmap_kwargs"]["num_samples"] for block in block_infos)), + "header": 1, # Standard format has headers + "offset_to_data_block": None, # Not needed + "memmap_specs": [block["memmap_kwargs"] for block in block_infos], + }) - timestamps_sampling_rate = self._nsx_basic_header[nsx_nb]["timestamp_resolution"] - timestamps_in_seconds = timestamps / timestamps_sampling_rate + return segments - time_differences = np.diff(timestamps_in_seconds) - gap_indices = np.argwhere(time_differences > detection_threshold).flatten() + def _segment_nsx_ptp(self, parsed_data_headers, nsx_nb): + """ + Segment PTP format data (v3.0) using sample-level gap detection. - # If gaps found, check user's tolerance - if len(gap_indices) > 0: - gap_report = self._format_gap_report(gap_indices, timestamps_in_seconds, time_differences, nsx_nb) + Every sample carries its own timestamp, so gaps are detected between samples and + each resulting segment is a slice of the single block the file stores. + """ + block_info = parsed_data_headers[0] + ts_kw = block_info["ptp_timestamps_memmap_kwargs"] + samples_kw = block_info["memmap_kwargs"] - # Error by default - user must opt-in to segmentation - if self.gap_tolerance_ms is None: - raise ValueError( - f"Detected {len(gap_indices)} timestamp gaps in ns{nsx_nb} file.\n" - f"{gap_report}\n" - f"To load this data, provide gap_tolerance_ms parameter to automatically " - f"segment at gaps larger than the specified tolerance." - ) + sampling_rate = self._nsx_sampling_frequency[nsx_nb] + ts_res = float(self._nsx_basic_header[nsx_nb]["timestamp_resolution"]) - # User provided tolerance - filter gaps and segment - gap_tolerance_s = self.gap_tolerance_ms / 1000.0 - significant_gap_mask = time_differences[gap_indices] > gap_tolerance_s - significant_gap_indices = gap_indices[significant_gap_mask] - - # Use significant gaps for segmentation (no warning - user opted in) - gap_indices = significant_gap_indices - - # Create segments based on gaps - num_total_samples = ts_kw["num_samples"] - segment_starts = np.hstack((0, gap_indices + 1)) - segment_boundaries = list(segment_starts) + [num_total_samples] - - packet_size = samples_kw["packet_size"] - base_samples_offset = samples_kw["offset"] - base_ts_offset = ts_kw["offset"] - - for seg_index, start in enumerate(segment_starts): - end = segment_boundaries[seg_index + 1] - seg_num_samples = end - start - - # Compute new file offset for this segment slice - seg_samples_offset = base_samples_offset + int(start) * packet_size - seg_ts_offset = base_ts_offset + int(start) * packet_size - - segments.append({ - "timestamp": None, # PTP timestamps read on demand - "nb_data_points": seg_num_samples, - "header": None, # PTP has no headers - "offset_to_data_block": None, - "memmap_kwargs": { - **samples_kw, - "offset": seg_samples_offset, - "num_samples": seg_num_samples, - }, - "timestamps_memmap_kwargs": { - **ts_kw, - "offset": seg_ts_offset, - "num_samples": seg_num_samples, - }, - }) - - # V2.1 or single block standard format: no segmentation needed - else: - segments.append({ - "timestamp": block_info["timestamp"], - "nb_data_points": block_info["memmap_kwargs"]["num_samples"], - "header": None, - "offset_to_data_block": None, - "memmap_kwargs": block_info["memmap_kwargs"], - }) + # Read timestamps via strided mmap view for gap detection + fid = self._get_nsx_fid(nsx_nb) + timestamps = self._create_mmap_view( + fid=fid, + dtype=ts_kw["dtype"], + offset=ts_kw["offset"], + num_samples=ts_kw["num_samples"], + packet_size=ts_kw.get("packet_size"), + item_offset=ts_kw.get("item_offset", 0), + ) + timestamps_in_seconds = timestamps.astype("float64") / ts_res + del timestamps # release the mmap buffer; timestamps_in_seconds is an independent copy + + time_differences = np.diff(timestamps_in_seconds) + boundary_indices = self._classify_gaps(time_differences, sampling_rate, nsx_nb, timestamps_in_seconds) + + num_total_samples = ts_kw["num_samples"] + segment_starts = np.hstack((0, boundary_indices + 1)).astype(np.intp) + segment_stops = np.hstack((boundary_indices + 1, num_total_samples)).astype(np.intp) + + packet_size = samples_kw["packet_size"] + base_samples_offset = samples_kw["offset"] + base_ts_offset = ts_kw["offset"] + + segments = [] + for start, stop in zip(segment_starts, segment_stops): + seg_num_samples = int(stop - start) + + # Compute new file offset for this segment slice + seg_samples_offset = base_samples_offset + int(start) * packet_size + seg_ts_offset = base_ts_offset + int(start) * packet_size + + segments.append({ + "timestamp": None, # PTP timestamps read on demand + "nb_data_points": seg_num_samples, + "header": None, # PTP has no headers + "offset_to_data_block": None, + "memmap_specs": [{ + **samples_kw, + "offset": seg_samples_offset, + "num_samples": seg_num_samples, + }], + "timestamps_memmap_kwargs": { + **ts_kw, + "offset": seg_ts_offset, + "num_samples": seg_num_samples, + }, + }) return segments diff --git a/neo/test/iotest/test_blackrockio.py b/neo/test/iotest/test_blackrockio.py index 3adccfd1d..74011d1a8 100644 --- a/neo/test/iotest/test_blackrockio.py +++ b/neo/test/iotest/test_blackrockio.py @@ -343,7 +343,9 @@ def test_segment_detection_pause(self): # And another one because there are spikes between segments with warnings.catch_warnings(record=True) as w: warnings.simplefilter("always") - reader = BlackrockIO(filename=filename, nsx_to_load=2, nev_override=filename_nev_outside_seg, gap_tolerance_ms=0) + reader = BlackrockIO( + filename=filename, nsx_to_load=2, nev_override=filename_nev_outside_seg, gap_tolerance_ms=0 + ) self.assertGreaterEqual(len(w), 2) # Check that warnings are correct diff --git a/neo/test/rawiotest/test_blackrockrawio.py b/neo/test/rawiotest/test_blackrockrawio.py index 03c3481ae..5485a031c 100644 --- a/neo/test/rawiotest/test_blackrockrawio.py +++ b/neo/test/rawiotest/test_blackrockrawio.py @@ -267,20 +267,30 @@ def test_get_blackrock_timestamps_ptp(self): np.testing.assert_array_equal(timestamps_ns6[:5], expected_ns6) def test_get_blackrock_timestamps_ptp_with_gaps(self): - """Test _get_blackrock_timestamps for PTP format with gaps and multiple segments.""" + """Test _get_blackrock_timestamps for PTP format with gaps and multiple segments. + + This file has three timestamp discontinuities, so it splits into four segments: + + - after sample 631: forward gap of ~0.97 ms (dropped samples) + - after sample 661: backward jump of ~1.9 ms (the clock runs backwards) + - after sample 689: forward gap of ~1.03 ms (dropped samples) + + Because of the backward jump, segment 2 starts earlier in time than segment 1 + ends, so the segments are not in ascending time order. + """ dirname = self.get_local_path("blackrock/blackrock_ptp_with_missing_samples/Hub1-NWBtestfile_neural_wspikes") gap_tolerance_ms = 0.5 reader = BlackrockRawIO(filename=dirname, nsx_to_load=6, gap_tolerance_ms=gap_tolerance_ms) reader.parse_header() n_segments = reader.segment_count(0) - self.assertEqual(n_segments, 3) + self.assertEqual(n_segments, 4) nanoseconds_per_second = 1_000_000_000.0 # This file has a single stream: ns6 (30 kHz, 1 channel) stream_index = 0 - expected_sizes = [632, 58, 310] + expected_sizes = [632, 30, 28, 310] # First 5 PTP clock values (nanoseconds since Unix epoch) per segment, read directly from file and hardcoded here for testing first_5_ptp_clock_ns_per_segment = [ @@ -304,6 +314,16 @@ def test_get_blackrock_timestamps_ptp_with_gaps(self): ], dtype="uint64", ), + np.array( + [ + 1752531864738809624, + 1752531864738842984, + 1752531864738876304, + 1752531864738909664, + 1752531864738942984, + ], + dtype="uint64", + ), np.array( [ 1752531864740742986, @@ -317,7 +337,7 @@ def test_get_blackrock_timestamps_ptp_with_gaps(self): ] # Last PTP clock value (nanoseconds since Unix epoch) of each segment, hardcoded here for testing last_ptp_clock_ns_per_segment = np.array( - [1752531864738776304, 1752531864739709625, 1752531864751042999], + [1752531864738776304, 1752531864740709666, 1752531864739709625, 1752531864751042999], dtype="uint64", ) @@ -333,15 +353,6 @@ def test_get_blackrock_timestamps_ptp_with_gaps(self): expected_last = last_ptp_clock_ns_per_segment[seg_index].astype("float64") / nanoseconds_per_second np.testing.assert_array_equal(timestamps[-1], expected_last) - # Verify the gaps between segments exceed the tolerance - for seg_index in range(n_segments - 1): - last_ts = last_ptp_clock_ns_per_segment[seg_index].astype("float64") / nanoseconds_per_second - first_ts_next = ( - first_5_ptp_clock_ns_per_segment[seg_index + 1][0].astype("float64") / nanoseconds_per_second - ) - gap_ms = (first_ts_next - last_ts) * 1000 - self.assertGreater(gap_ms, gap_tolerance_ms) - def test_gap_tolerance_ms_parameter(self): """ Test gap_tolerance_ms parameter for gap handling with files that have actual gaps. @@ -349,7 +360,8 @@ def test_gap_tolerance_ms_parameter(self): Tests the error-by-default behavior where files with timestamp gaps raise ValueError unless the user explicitly opts in with gap_tolerance_ms parameter. - See PR #1769 for the gap details on the example file used here. + See PR #1769 for the gap details on the example file used here. This file has two + forward gaps (~0.97 ms and ~1.03 ms) and one backward jump (~1.9 ms). """ # Use stubbed files with missing samples (timestamp gaps) from SimulatedSpikes data dirname = self.get_local_path("blackrock/blackrock_ptp_with_missing_samples/Hub1-NWBtestfile_neural_wspikes") @@ -360,19 +372,91 @@ def test_gap_tolerance_ms_parameter(self): reader = BlackrockRawIO(filename=dirname, nsx_to_load=6) reader.parse_header() - # Test 2: Explicit tolerance allows loading files with gaps - # User opts in by providing gap_tolerance_ms + # Test 2: A tolerance above every forward gap keeps the data either side of them + # together, but the backward jump splits regardless of the tolerance reader_with_tolerance = BlackrockRawIO(filename=dirname, nsx_to_load=6, gap_tolerance_ms=10.0) reader_with_tolerance.parse_header() segments_with_tolerance = reader_with_tolerance.segment_count(0) - self.assertEqual(1, segments_with_tolerance) # Gaps < 10ms are ignored + self.assertEqual(2, segments_with_tolerance) - # Test 3: Stricter tolerance creates 3 segments - # With strict tolerance (0.5ms), gaps > 0.5ms will create new segments + # Test 3: Stricter tolerance creates 4 segments + # With strict tolerance (0.5ms), both forward gaps split as well as the backward jump reader_strict = BlackrockRawIO(filename=dirname, nsx_to_load=6, gap_tolerance_ms=0.5) reader_strict.parse_header() segments_strict = reader_strict.segment_count(0) - self.assertEqual(segments_strict, 3) # + self.assertEqual(segments_strict, 4) + + def test_gap_tolerance_ms_standard_format(self): + """ + Test gap_tolerance_ms on a standard-format (v2.2/2.3/3.0) file. + + The pause_correct file has two data blocks separated by a ~27 s recording pause. + Unlike the PTP path, gaps here are detected between block timestamps rather than + between per-sample timestamps. + """ + dirname = self.get_local_path("blackrock/segment/PauseCorrect/pause_correct") + + # Default behavior (None) raises, same as for PTP files + with self.assertRaises(ValueError): + reader = BlackrockRawIO(filename=dirname, nsx_to_load=2) + reader.parse_header() + + # A tolerance below the pause keeps one segment per block + reader_strict = BlackrockRawIO(filename=dirname, nsx_to_load=2, gap_tolerance_ms=0) + reader_strict.parse_header() + self.assertEqual(reader_strict.segment_count(0), 2) + + # A tolerance above the pause merges both blocks into a single segment + reader_merged = BlackrockRawIO(filename=dirname, nsx_to_load=2, gap_tolerance_ms=100_000) + reader_merged.parse_header() + self.assertEqual(reader_merged.segment_count(0), 1) + + # The merged segment spans both blocks, which sit at different file offsets + stream_index = 0 + self.assertEqual(reader_merged.get_signal_size(0, 0, stream_index), 8000) # 4000 + 4000 + + def test_gap_tolerance_ms_standard_format_merged_chunk(self): + """ + Reading across blocks merged into one segment must stitch them together. + + A merged segment holds one memmap spec per block, so a requested sample range can + straddle several of them. Compare reads from the merged segment against the same + data read as two separate segments, which needs no stitching. + """ + dirname = self.get_local_path("blackrock/segment/PauseCorrect/pause_correct") + stream_index = 0 + + reader_split = BlackrockRawIO(filename=dirname, nsx_to_load=2, gap_tolerance_ms=0) + reader_split.parse_header() + block_0 = reader_split.get_analogsignal_chunk(seg_index=0, stream_index=stream_index) + block_1 = reader_split.get_analogsignal_chunk(seg_index=1, stream_index=stream_index) + expected = np.concatenate([block_0, block_1], axis=0) + + reader_merged = BlackrockRawIO(filename=dirname, nsx_to_load=2, gap_tolerance_ms=100_000) + reader_merged.parse_header() + + # Full read + full = reader_merged.get_analogsignal_chunk(seg_index=0, stream_index=stream_index) + self.assertEqual(full.shape[0], 8000) + np.testing.assert_array_equal(full, expected) + + # A range straddling the boundary between the two blocks + straddling = reader_merged.get_analogsignal_chunk( + seg_index=0, stream_index=stream_index, i_start=3990, i_stop=4010 + ) + np.testing.assert_array_equal(straddling, expected[3990:4010]) + + # A range wholly inside the second block + second_block_only = reader_merged.get_analogsignal_chunk( + seg_index=0, stream_index=stream_index, i_start=5000, i_stop=5100 + ) + np.testing.assert_array_equal(second_block_only, expected[5000:5100]) + + # A single channel selection across the boundary + one_channel = reader_merged.get_analogsignal_chunk( + seg_index=0, stream_index=stream_index, i_start=3990, i_stop=4010, channel_indexes=[0] + ) + np.testing.assert_array_equal(one_channel, expected[3990:4010, [0]]) if __name__ == "__main__": From f96b334d708f77a587660281ccd3e2a191904554 Mon Sep 17 00:00:00 2001 From: Heberto Mayorquin Date: Thu, 16 Jul 2026 17:54:00 -0600 Subject: [PATCH 5/9] Report gaps in sample terms, and make the report format-agnostic _format_gap_report derived a gap's position from a per-sample timestamp array, asking that array to be both the time of the sample before each gap and, via its first element, the start of the recording. Per-sample timestamps satisfy both, so the PTP path it was written for was correct. The standard path has no such array: its file stores one timestamp per block. It passed block timestamps, which got the origin right but not the position, so a gap was reported at the start of the block preceding it rather than at the sample it follows. In pause_correct that put a gap 4 s into the recording at 0.000000 s. Pass the resolved quantities instead, so the report renders what it is given and assumes nothing about how a format stores its timing. The standard path works out the sample before each gap from its block scalars, which is exact because within a block the timing is uniform by definition, and avoids synthesizing per-sample timestamps that would run to gigabytes on a long recording. Both paths now report in the same sample-level terms, which keeps the report reusable as the gap API reaches other readers (issue #1773). Also only build the report when about to raise, rather than on every load. --- neo/rawio/blackrockrawio.py | 95 ++++++++++++++++------- neo/test/rawiotest/test_blackrockrawio.py | 19 +++++ 2 files changed, 86 insertions(+), 28 deletions(-) diff --git a/neo/rawio/blackrockrawio.py b/neo/rawio/blackrockrawio.py index f4517de4e..8f9c6a0b2 100644 --- a/neo/rawio/blackrockrawio.py +++ b/neo/rawio/blackrockrawio.py @@ -1418,40 +1418,48 @@ def _parse_nsx_data_v30_ptp(self, nsx_nb): } ] - def _format_gap_report(self, gap_indices, timestamps_in_seconds, time_differences, nsx_nb): + def _format_gap_report( + self, gap_sample_indices, gap_positions_seconds, gap_durations_ms, stream_label, detection_threshold_ms + ): """ Format a detailed gap report showing where timestamp discontinuities occur. + This renders a report and nothing else: every quantity is passed in already + resolved, so it makes no assumption about how a format stores its timing. A + reader whose file carries one timestamp per sample and a reader whose file + carries one per block describe their gaps in the same sample-level terms, and + each works out those terms whichever way suits its own storage. That keeps the + report usable across readers as the gap API spreads to other formats (see + python-neo issue #1773). + Parameters ---------- - gap_indices : np.ndarray - Indices where gaps were detected - timestamps_in_seconds : np.ndarray - All timestamps converted to seconds - time_differences : np.ndarray - Time differences between consecutive timestamps - nsx_nb : int - NSX file number for the report + gap_sample_indices : np.ndarray + Index of the last sample before each gap. + gap_positions_seconds : np.ndarray + Time of that sample, in seconds since the start of the recording. + gap_durations_ms : np.ndarray + Length of each gap in milliseconds. Negative when the clock jumped backwards. + stream_label : str + Name of the stream the gaps were found in, e.g. "ns2". + detection_threshold_ms : float + Threshold that was used to detect the gaps, in milliseconds. Returns ------- str Formatted gap report with table """ - # Calculate gap details - gap_durations_seconds = time_differences[gap_indices] - gap_durations_ms = gap_durations_seconds * 1000 - gap_positions_seconds = timestamps_in_seconds[gap_indices] - timestamps_in_seconds[0] - # Build gap detail table gap_detail_lines = [ f"| {index:>15,} | {pos:>21.6f} | {dur:>21.3f} |\n" - for index, pos, dur in zip(gap_indices, gap_positions_seconds, gap_durations_ms) + for index, pos, dur in zip(gap_sample_indices, gap_positions_seconds, gap_durations_ms) ] return ( - f"Gap Report for ns{nsx_nb}:\n" - f"Found {len(gap_indices)} timestamp gaps (detection threshold: 2 x sampling period)\n\n" + f"Gap Report for {stream_label}:\n" + f"Found {len(gap_sample_indices)} timestamp gaps " + f"(detection threshold: {detection_threshold_ms:.6g} ms)\n\n" "Gap Details:\n" "+-----------------+-----------------------+-----------------------+\n" "| Sample Index | Sample at (Seconds) | Gap Size (ms) |\n" @@ -1506,13 +1514,17 @@ def _segment_nsx_data(self, parsed_data_headers, nsx_nb): else: return self._segment_nsx_v22_v30(parsed_data_headers, nsx_nb) - def _classify_gaps(self, time_differences, sampling_rate, nsx_nb, timestamps_in_seconds): + def _classify_gaps(self, time_differences, sampling_rate, nsx_nb, sample_indices, sample_positions_seconds): """ Classify timestamp discontinuities into segment boundaries. - Shared by the standard and PTP segmenters. Both hand over the same quantity: - the interval between the timestamps of two consecutive samples, which is one - sampling period when the data is contiguous. + Shared by the standard and PTP segmenters. Both hand over the interval between + the timestamps of two consecutive samples, which is one sampling period when the + data is contiguous, along with a sample-level description of where each interval + sits. The PTP path reads those descriptions from the per-sample timestamps in the + file; the standard path works them out from its block timestamps, since within a + block the timing is uniform by definition. Either way the gaps are detected, and + reported, in the same terms. Two kinds of discontinuity are treated differently: @@ -1532,14 +1544,17 @@ def _classify_gaps(self, time_differences, sampling_rate, nsx_nb, timestamps_in_ Sampling rate in Hz. nsx_nb : int NSX file number (for error reporting). - timestamps_in_seconds : np.ndarray - Timestamps in seconds, used to report where each gap falls. + sample_indices : np.ndarray + For each interval, the index of the sample that precedes it. + sample_positions_seconds : np.ndarray + For each interval, the time of that sample in seconds since the start of the + recording. Returns ------- np.ndarray Indices into time_differences at which a segment boundary falls; a boundary - at index i means a new segment starts at sample i + 1. + at index i means a new segment starts at i + 1. """ # Detection threshold: use strict 2x sampling period to find ALL gaps detection_threshold = 2.0 / sampling_rate @@ -1551,10 +1566,15 @@ def _classify_gaps(self, time_differences, sampling_rate, nsx_nb, timestamps_in_ if len(gap_indices) == 0: return gap_indices - gap_report = self._format_gap_report(gap_indices, timestamps_in_seconds, time_differences, nsx_nb) - # Error by default - user must opt-in to segmentation if self.gap_tolerance_ms is None: + gap_report = self._format_gap_report( + gap_sample_indices=sample_indices[gap_indices], + gap_positions_seconds=sample_positions_seconds[gap_indices], + gap_durations_ms=time_differences[gap_indices] * 1000, + stream_label=f"ns{nsx_nb}", + detection_threshold_ms=detection_threshold * 1000, + ) raise ValueError( f"Detected {len(gap_indices)} timestamp gaps in ns{nsx_nb} file.\n" f"{gap_report}\n" @@ -1608,7 +1628,18 @@ def _segment_nsx_v22_v30(self, parsed_data_headers, nsx_nb): # the per-block equivalent of np.diff() over the PTP per-sample timestamps. block_last_sample_times = block_t_starts + (block_sizes - 1) / sampling_rate inter_block_intervals = block_t_starts[1:] - block_last_sample_times[:-1] - boundary_indices = self._classify_gaps(inter_block_intervals, sampling_rate, nsx_nb, block_t_starts) + + # Describe each candidate boundary by the sample before it, the same way the + # PTP path does. Within a block the timing is uniform by definition, so the + # last sample of block i is sample cumsum(sizes)[i] - 1 and its time follows + # from the block's own timestamp. Working these out from the block scalars + # gives exactly what per-sample timestamps would say, without building them. + last_sample_indices = np.cumsum(block_sizes)[:-1] - 1 + last_sample_positions = block_last_sample_times[:-1] - block_t_starts[0] + + boundary_indices = self._classify_gaps( + inter_block_intervals, sampling_rate, nsx_nb, last_sample_indices, last_sample_positions + ) else: boundary_indices = np.array([], dtype=np.intp) @@ -1657,7 +1688,15 @@ def _segment_nsx_ptp(self, parsed_data_headers, nsx_nb): del timestamps # release the mmap buffer; timestamps_in_seconds is an independent copy time_differences = np.diff(timestamps_in_seconds) - boundary_indices = self._classify_gaps(time_differences, sampling_rate, nsx_nb, timestamps_in_seconds) + + # Every sample carries its own timestamp, so the sample before each interval and + # its time are read straight off the array + last_sample_indices = np.arange(len(time_differences)) + last_sample_positions = timestamps_in_seconds[:-1] - timestamps_in_seconds[0] + + boundary_indices = self._classify_gaps( + time_differences, sampling_rate, nsx_nb, last_sample_indices, last_sample_positions + ) num_total_samples = ts_kw["num_samples"] segment_starts = np.hstack((0, boundary_indices + 1)).astype(np.intp) diff --git a/neo/test/rawiotest/test_blackrockrawio.py b/neo/test/rawiotest/test_blackrockrawio.py index 5485a031c..b9c2c4609 100644 --- a/neo/test/rawiotest/test_blackrockrawio.py +++ b/neo/test/rawiotest/test_blackrockrawio.py @@ -415,6 +415,25 @@ def test_gap_tolerance_ms_standard_format(self): stream_index = 0 self.assertEqual(reader_merged.get_signal_size(0, 0, stream_index), 8000) # 4000 + 4000 + def test_gap_report_locates_gap_in_sample_terms(self): + """ + The gap report must locate gaps by sample, whatever the file stores. + + A standard-format file has one timestamp per block rather than one per sample, but + the report describes gaps the same way the PTP report does. In pause_correct the + first block holds 4000 samples at 1 kHz, so the gap follows sample 3999 at 3.999 s + rather than sitting at the start of the recording. + """ + dirname = self.get_local_path("blackrock/segment/PauseCorrect/pause_correct") + with self.assertRaises(ValueError) as context: + reader = BlackrockRawIO(filename=dirname, nsx_to_load=2) + reader.parse_header() + + report = str(context.exception) + self.assertIn("3,999", report) # last sample before the gap, not block index 0 + self.assertIn("3.999000", report) # when that sample was taken, not 0.000000 + self.assertIn("27009.700", report) # gap duration in ms + def test_gap_tolerance_ms_standard_format_merged_chunk(self): """ Reading across blocks merged into one segment must stitch them together. From 223688fc0344be384a568043ca27be755ddd6f3c Mon Sep 17 00:00:00 2001 From: Heberto Mayorquin Date: Thu, 16 Jul 2026 18:53:35 -0600 Subject: [PATCH 6/9] Fill one array in _read_multi_block_chunk instead of concatenating An empty range touched no block, leaving nothing to concatenate and raising ValueError, where an unmerged segment returns an empty array. Filling a preallocated output gives the empty range a correctly shaped result and keeps merged and unmerged segments behaving alike. The blocks are also no longer all held at once. With a channel list each block is a copy rather than a view, so concatenating held every one of them; now each is written into the output and released. A range inside a single block is still returned as a view, as before. --- neo/rawio/blackrockrawio.py | 25 ++++++++++++++++++----- neo/test/rawiotest/test_blackrockrawio.py | 10 +++++++++ 2 files changed, 30 insertions(+), 5 deletions(-) diff --git a/neo/rawio/blackrockrawio.py b/neo/rawio/blackrockrawio.py index 8f9c6a0b2..37218a274 100644 --- a/neo/rawio/blackrockrawio.py +++ b/neo/rawio/blackrockrawio.py @@ -791,6 +791,9 @@ def _read_multi_block_chunk(self, fid, memmap_specs, channels, i_start, i_stop, ------- np.ndarray Signal data of shape (i_stop - i_start, len(channel_indexes)), dtype int16. + A range falling inside a single block is returned as a view, as it would be + from a segment that was never merged; otherwise the blocks are copied into + one array. """ total_samples = sum(spec["num_samples"] for spec in memmap_specs) if i_start is None: @@ -798,7 +801,9 @@ def _read_multi_block_chunk(self, fid, memmap_specs, channels, i_start, i_stop, if i_stop is None: i_stop = total_samples - pieces = [] + # Map only the blocks the range touches. Slicing the rows keeps these as views, + # so nothing is read from disk until the data is actually used below. + block_views = [] block_start = 0 for spec in memmap_specs: block_stop = block_start + spec["num_samples"] @@ -814,12 +819,22 @@ def _read_multi_block_chunk(self, fid, memmap_specs, channels, i_start, i_stop, ) local_start = max(0, i_start - block_start) local_stop = min(spec["num_samples"], i_stop - block_start) - pieces.append(data[local_start:local_stop, channel_indexes]) + block_views.append(data[local_start:local_stop]) block_start = block_stop - if len(pieces) == 1: - return pieces[0] - return np.concatenate(pieces, axis=0) + if len(block_views) == 1: + return block_views[0][:, channel_indexes] + + # Fill one output array rather than concatenating, so the blocks are never all + # held at once, and a range touching no block still gives a correctly shaped + # empty result rather than failing + selected_channel_count = np.arange(channels)[channel_indexes].size + out = np.empty((i_stop - i_start, selected_channel_count), dtype=memmap_specs[0]["dtype"]) + position = 0 + for block_view in block_views: + out[position : position + len(block_view)] = block_view[:, channel_indexes] + position += len(block_view) + return out def _get_blackrock_timestamps(self, block_index, seg_index, i_start, i_stop, stream_index): """ diff --git a/neo/test/rawiotest/test_blackrockrawio.py b/neo/test/rawiotest/test_blackrockrawio.py index b9c2c4609..894c37dff 100644 --- a/neo/test/rawiotest/test_blackrockrawio.py +++ b/neo/test/rawiotest/test_blackrockrawio.py @@ -477,6 +477,16 @@ def test_gap_tolerance_ms_standard_format_merged_chunk(self): ) np.testing.assert_array_equal(one_channel, expected[3990:4010, [0]]) + # An empty range must behave the same as it does on an unmerged segment, which + # returns an empty array rather than failing + empty_merged = reader_merged.get_analogsignal_chunk( + seg_index=0, stream_index=stream_index, i_start=4000, i_stop=4000 + ) + empty_split = reader_split.get_analogsignal_chunk( + seg_index=0, stream_index=stream_index, i_start=4000, i_stop=4000 + ) + self.assertEqual(empty_merged.shape, empty_split.shape) + if __name__ == "__main__": unittest.main() From f5a480abfbac2550361f5923c970329d639263cf Mon Sep 17 00:00:00 2001 From: Heberto Mayorquin Date: Thu, 16 Jul 2026 20:29:44 -0600 Subject: [PATCH 7/9] Keep every merged block's timestamp so the accessor stays correct _get_blackrock_timestamps reconstructed non-PTP times as t_start + index/rate. That was correct only because segmentation gave every data block its own segment, and within a block the sampling is uniform by definition. Nothing stated that dependency. Merging blocks under gap_tolerance_ms removes it. A merged segment has one t_start and a pause somewhere inside it, so the arithmetic walks straight through the pause and reports every later sample early by its length. The accessor returned 4.0 s for a sample pause_correct timestamps at 31.0087 s, without a line of it changing. Keep each contributing block's timestamp on the segment and interpolate per block, which is where uniform spacing is a fact of the format rather than an artifact of the merge. gap_tolerance_ms decides how samples are grouped into segments; it does not decide what time a sample was recorded at. The segment's own t_start and t_stop were always right, being the first and last sample, and stay so. What a merged segment cannot express is the pause in between, which is exactly why these timestamps are worth reporting separately. --- neo/rawio/blackrockrawio.py | 43 ++++++++++++++++++++--- neo/test/rawiotest/test_blackrockrawio.py | 37 +++++++++++++++++++ 2 files changed, 76 insertions(+), 4 deletions(-) diff --git a/neo/rawio/blackrockrawio.py b/neo/rawio/blackrockrawio.py index 37218a274..8da92e5e5 100644 --- a/neo/rawio/blackrockrawio.py +++ b/neo/rawio/blackrockrawio.py @@ -854,11 +854,24 @@ def _get_blackrock_timestamps(self, block_index, seg_index, i_start, i_stop, str - **Standard formats (FileSpec 2.2, 2.3, 3.0 non-PTP):** Each data block has a single scalar timestamp for its first sample. All subsequent samples within the block are interpolated as - ``t_start + sample_index / sampling_rate``, assuming uniform spacing. + ``block_timestamp + sample_index_within_block / sampling_rate``, + assuming uniform spacing. Uniform spacing holds within a block + because that is what the format means by a block; it does not hold + across blocks, which are separated by a pause of recorded length. + A segment built with ``gap_tolerance_ms`` may merge several blocks, + so the interpolation is done per block and steps across each pause. - **FileSpec 2.1:** No timestamps are stored in the file. All samples are interpolated from ``t_start=0`` using the sampling rate. + The times returned here do not depend on ``gap_tolerance_ms``. That + parameter decides how samples are grouped into segments, not what time + a sample was recorded at. When a segment merges blocks across a pause, + its single ``t_start`` cannot express that pause, so the segment's own + ``t_start + sample_index / sampling_rate`` is wrong for every sample + after it while the times returned here stay correct. That gap is the + reason this method exists. + Parameters ---------- block_index : int @@ -901,10 +914,31 @@ def _get_blackrock_timestamps(self, block_index, seg_index, i_start, i_stop, str ts_res = float(self._nsx_basic_header[nsx_nb]["timestamp_resolution"]) return timestamps[i_start:i_stop].astype("float64") / ts_res else: - # Non-PTP: reconstruct from t_start + index / sampling_rate - t_start = self._sigs_t_starts[nsx_nb][seg_index] sr = self._nsx_sampling_frequency[nsx_nb] - return t_start + np.arange(i_start, i_stop, dtype="float64") / sr + block_timestamps = seg.get("block_timestamps") + + if block_timestamps is None or len(block_timestamps) == 1: + # One block, so the segment's own t_start is the block's recorded start + # and the whole stretch is uniform from it + t_start = self._sigs_t_starts[nsx_nb][seg_index] + return t_start + np.arange(i_start, i_stop, dtype="float64") / sr + + # Several blocks merged into this segment. The segment has a single t_start + # and cannot express the pauses between them, but each block recorded its own + # start, so use those: the times returned here step across each pause even + # though the signal's uniform time base does not. + ts_res = float(self._nsx_basic_header[nsx_nb]["timestamp_resolution"]) + timestamps = np.empty(i_stop - i_start, dtype="float64") + block_start = 0 + for block_timestamp, spec in zip(block_timestamps, seg["memmap_specs"]): + block_stop = block_start + spec["num_samples"] + first = max(i_start, block_start) + last = min(i_stop, block_stop) + if last > first: + within_block = np.arange(first - block_start, last - block_start, dtype="float64") + timestamps[first - i_start : last - i_start] = float(block_timestamp) / ts_res + within_block / sr + block_start = block_stop + return timestamps def _spike_count(self, block_index, seg_index, unit_index): channel_id, unit_id = self.internal_unit_ids[unit_index] @@ -1667,6 +1701,7 @@ def _segment_nsx_v22_v30(self, parsed_data_headers, nsx_nb): block_infos = parsed_data_headers[start:stop] segments.append({ "timestamp": block_infos[0]["timestamp"], + "block_timestamps": [block["timestamp"] for block in block_infos], "nb_data_points": int(sum(block["memmap_kwargs"]["num_samples"] for block in block_infos)), "header": 1, # Standard format has headers "offset_to_data_block": None, # Not needed diff --git a/neo/test/rawiotest/test_blackrockrawio.py b/neo/test/rawiotest/test_blackrockrawio.py index 894c37dff..dfc1f216d 100644 --- a/neo/test/rawiotest/test_blackrockrawio.py +++ b/neo/test/rawiotest/test_blackrockrawio.py @@ -415,6 +415,43 @@ def test_gap_tolerance_ms_standard_format(self): stream_index = 0 self.assertEqual(reader_merged.get_signal_size(0, 0, stream_index), 8000) # 4000 + 4000 + def test_blackrock_timestamps_do_not_depend_on_gap_tolerance(self): + """ + The timestamps accessor reports what the file recorded, whatever gap_tolerance_ms says. + + gap_tolerance_ms decides how samples are grouped into segments. It must not decide what + time a sample is reported at. Merging blocks into one segment gives that segment a single + t_start, which cannot express the pause between them, but every block recorded its own + start and the accessor still reports from those. + """ + dirname = self.get_local_path("blackrock/segment/PauseCorrect/pause_correct") + stream_index = 0 + + reader_split = BlackrockRawIO(filename=dirname, nsx_to_load=2, gap_tolerance_ms=0) + reader_split.parse_header() + expected = np.concatenate( + [ + reader_split._get_blackrock_timestamps(0, 0, None, None, stream_index), + reader_split._get_blackrock_timestamps(0, 1, None, None, stream_index), + ] + ) + + reader_merged = BlackrockRawIO(filename=dirname, nsx_to_load=2, gap_tolerance_ms=100_000) + reader_merged.parse_header() + merged = reader_merged._get_blackrock_timestamps(0, 0, None, None, stream_index) + + # Same samples, same recorded times, whether or not the blocks were merged + np.testing.assert_array_equal(merged, expected) + + # The pause is still visible in the timestamps even though the segment's uniform + # time base has no way to show it + self.assertAlmostEqual(merged[4000] - merged[3999], 27.0097, places=3) + + # Ranges straddling the block boundary agree with the full read + for i_start, i_stop in [(3990, 4010), (4000, 4001), (5000, 5100)]: + chunk = reader_merged._get_blackrock_timestamps(0, 0, i_start, i_stop, stream_index) + np.testing.assert_array_equal(chunk, expected[i_start:i_stop]) + def test_gap_report_locates_gap_in_sample_terms(self): """ The gap report must locate gaps by sample, whatever the file stores. From 69f66049183b9f7f9949ffa0968e99ecaa273d9d Mon Sep 17 00:00:00 2001 From: Heberto Mayorquin Date: Thu, 16 Jul 2026 21:20:09 -0600 Subject: [PATCH 8/9] Warn on backward jumps instead of raising, and keep them out of the gap report Raising by default earns its place when the reader cannot decide for the caller. A forward gap qualifies: merging across one distorts the segment's uniform time base, splitting gives more segments, and only the caller knows which is acceptable. A backward jump does not. It has no duration to compare against gap_tolerance_ms, so it always becomes a segment boundary and there is nothing to choose. Raising asked a question with one possible answer, and pointed the caller at a parameter that has no effect on it: whatever value they supplied, the answer was the same. Warn instead, naming what was found and what was done about it, and leave backward jumps out of the gap report, which is about gaps. The error now reports only what gap_tolerance_ms can act on, so its advice is true of everything in it. This narrows the breaking change to files with real forward gaps. A file whose only discontinuity is a clock reset loads as it did before, into the same segments, and test_segment_detection_reset goes back to master unchanged. The nev path already warns about resets on these files, so the signal path now says the same thing. --- neo/rawio/blackrockrawio.py | 48 ++++++++++++++++------- neo/test/iotest/test_blackrockio.py | 4 +- neo/test/rawiotest/test_blackrockrawio.py | 47 ++++++++++++++++++++++ 3 files changed, 83 insertions(+), 16 deletions(-) diff --git a/neo/rawio/blackrockrawio.py b/neo/rawio/blackrockrawio.py index 8da92e5e5..05edb8910 100644 --- a/neo/rawio/blackrockrawio.py +++ b/neo/rawio/blackrockrawio.py @@ -1575,15 +1575,19 @@ def _classify_gaps(self, time_differences, sampling_rate, nsx_nb, sample_indices block the timing is uniform by definition. Either way the gaps are detected, and reported, in the same terms. - Two kinds of discontinuity are treated differently: + Two kinds of discontinuity are treated differently, because only one of them + leaves the caller anything to decide: - Forward gaps (a pause, or dropped samples) are detected with a fixed threshold of 2 sampling periods, which sits above the jitter and rounding - noise floor of these files. Whether a detected forward gap actually becomes - a segment boundary is then up to gap_tolerance_ms. - - Backward jumps (the clock runs backwards between consecutive samples) are an - abnormal condition rather than missing data. They have no gap duration to - compare against a tolerance, so they always become segment boundaries. + noise floor of these files. Merging across one distorts the segment's uniform + time base while splitting gives more segments, and only the caller knows which + is acceptable, so with gap_tolerance_ms unset this raises and reports them. + - Backward jumps (the clock runs backwards between consecutive samples, which a + recording reset causes) have no duration to compare against a tolerance, so + there is no choice to offer: they always become segment boundaries. Raising + would pose a question with one possible answer, so they are warned about + instead, and left out of the gap report, which is about gaps. Parameters ---------- @@ -1610,22 +1614,38 @@ def _classify_gaps(self, time_differences, sampling_rate, nsx_nb, sample_indices forward_mask = time_differences > detection_threshold backward_mask = time_differences < 0 - gap_indices = np.flatnonzero(forward_mask | backward_mask) - if len(gap_indices) == 0: - return gap_indices + # Backward jumps always split, so there is nothing to ask about. Say what was + # found and what was done about it, and carry on. + backward_indices = np.flatnonzero(backward_mask) + if len(backward_indices) > 0: + locations = ", ".join( + f"sample {sample_indices[index]:,} at {sample_positions_seconds[index]:.6f} s " + f"(back {-time_differences[index] * 1000:.3f} ms)" + for index in backward_indices + ) + warnings.warn( + f"Detected {len(backward_indices)} backward timestamp jump(s) in ns{nsx_nb}: {locations}. " + f"The clock ran backwards, which a recording reset causes. Each one starts a new segment; " + f"gap_tolerance_ms does not apply to them because a backward jump has no gap duration." + ) + + forward_indices = np.flatnonzero(forward_mask) + if len(forward_indices) == 0: + return backward_indices - # Error by default - user must opt-in to segmentation + # Error by default on forward gaps - merging across one or splitting at it is a + # choice only the caller can make if self.gap_tolerance_ms is None: gap_report = self._format_gap_report( - gap_sample_indices=sample_indices[gap_indices], - gap_positions_seconds=sample_positions_seconds[gap_indices], - gap_durations_ms=time_differences[gap_indices] * 1000, + gap_sample_indices=sample_indices[forward_indices], + gap_positions_seconds=sample_positions_seconds[forward_indices], + gap_durations_ms=time_differences[forward_indices] * 1000, stream_label=f"ns{nsx_nb}", detection_threshold_ms=detection_threshold * 1000, ) raise ValueError( - f"Detected {len(gap_indices)} timestamp gaps in ns{nsx_nb} file.\n" + f"Detected {len(forward_indices)} timestamp gaps in ns{nsx_nb} file.\n" f"{gap_report}\n" f"To load this data, provide gap_tolerance_ms parameter to automatically " f"segment at gaps larger than the specified tolerance." diff --git a/neo/test/iotest/test_blackrockio.py b/neo/test/iotest/test_blackrockio.py index 74011d1a8..0579a187a 100644 --- a/neo/test/iotest/test_blackrockio.py +++ b/neo/test/iotest/test_blackrockio.py @@ -289,12 +289,12 @@ def test_segment_detection_reset(self): # This fails, because in the nev there is no way to separate two segments with self.assertRaises(NeoReadWriteError): - reader = BlackrockIO(filename=filename, nsx_to_load=2, nev_override=filename_nev_fail, gap_tolerance_ms=0) + reader = BlackrockIO(filename=filename, nsx_to_load=2, nev_override=filename_nev_fail) # The correct file will issue a warning because a reset has occurred # and could be detected, but was not explicitly documented in the file with warnings.catch_warnings(record=True) as w: - reader = BlackrockIO(filename=filename, nsx_to_load=2, gap_tolerance_ms=0) + reader = BlackrockIO(filename=filename, nsx_to_load=2) self.assertGreaterEqual(len(w), 1) messages = [str(warning.message) for warning in w if warning.category == UserWarning] self.assertIn("Detected 1 undocumented segments within nev data after " "timestamps [5451].", messages) diff --git a/neo/test/rawiotest/test_blackrockrawio.py b/neo/test/rawiotest/test_blackrockrawio.py index dfc1f216d..194beaaa4 100644 --- a/neo/test/rawiotest/test_blackrockrawio.py +++ b/neo/test/rawiotest/test_blackrockrawio.py @@ -3,6 +3,7 @@ """ import unittest +import warnings from neo.rawio.blackrockrawio import BlackrockRawIO from neo.test.rawiotest.common_rawio_test import BaseTestRawIO @@ -386,6 +387,52 @@ def test_gap_tolerance_ms_parameter(self): segments_strict = reader_strict.segment_count(0) self.assertEqual(segments_strict, 4) + def test_backward_jump_warns_and_splits_without_gap_tolerance_ms(self): + """ + A backward jump warns rather than raising, and splits whatever the tolerance says. + + Raising exists so the caller can choose between merging across a gap and splitting + at it. A backward jump offers no such choice: it has no duration to compare against + gap_tolerance_ms, so it always becomes a segment boundary. Raising would pose a + question with one possible answer, so the reader reports it and carries on. + + The reset file's clock restarts, so its two blocks are separated by a backward jump + rather than a gap. + """ + dirname = self.get_local_path("blackrock/segment/ResetCorrect/reset") + + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + reader = BlackrockRawIO(filename=dirname, nsx_to_load=2) + reader.parse_header() + + self.assertEqual(reader.segment_count(0), 2) + messages = [str(warning.message) for warning in caught] + self.assertTrue(any("backward timestamp jump" in message for message in messages)) + + # No tolerance merges the blocks either side of a reset + for gap_tolerance_ms in [0, 100.0, 10_000_000]: + reader = BlackrockRawIO(filename=dirname, nsx_to_load=2, gap_tolerance_ms=gap_tolerance_ms) + reader.parse_header() + self.assertEqual(reader.segment_count(0), 2) + + def test_gap_report_covers_gaps_only(self): + """ + The gap report is about gaps, so a backward jump is not counted among them. + + The PTP file holds two forward gaps and one backward jump. Only the gaps are + something gap_tolerance_ms can act on, so only they are reported and counted. + """ + dirname = self.get_local_path("blackrock/blackrock_ptp_with_missing_samples/Hub1-NWBtestfile_neural_wspikes") + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + with self.assertRaises(ValueError) as context: + reader = BlackrockRawIO(filename=dirname, nsx_to_load=6) + reader.parse_header() + + # Two forward gaps, not three discontinuities + self.assertIn("Detected 2 timestamp gaps", str(context.exception)) + def test_gap_tolerance_ms_standard_format(self): """ Test gap_tolerance_ms on a standard-format (v2.2/2.3/3.0) file. From 47a46fe06a654f5c39875699c64ced9f157d2380 Mon Sep 17 00:00:00 2001 From: Heberto Mayorquin Date: Thu, 16 Jul 2026 21:34:40 -0600 Subject: [PATCH 9/9] diff minimization --- neo/rawio/blackrockrawio.py | 135 +++++++--------------- neo/test/rawiotest/test_blackrockrawio.py | 66 +++-------- 2 files changed, 58 insertions(+), 143 deletions(-) diff --git a/neo/rawio/blackrockrawio.py b/neo/rawio/blackrockrawio.py index 05edb8910..c50b725ca 100644 --- a/neo/rawio/blackrockrawio.py +++ b/neo/rawio/blackrockrawio.py @@ -124,10 +124,6 @@ class BlackrockRawIO(BaseRawIO): - PTP format (v3.0-ptp): Gaps in per-sample timestamps - Standard format (v2.2/2.3/3.0): Gaps between data blocks - Backward jumps, where the clock runs backwards from one sample to the next, - are not gaps and have no duration to compare against this tolerance. They - always create a segment boundary, whatever the value given here. - Notes ----- * Note: This routine will handle files according to specification 2.1, 2.2, @@ -767,10 +763,8 @@ def _read_multi_block_chunk(self, fid, memmap_specs, channels, i_start, i_stop, """ Read a chunk from a segment that merges several data blocks. - In the standard format, consecutive blocks that are not separated by a - significant gap are merged into one segment. Each block lives at its own file - offset, so a contiguous sample range can straddle several of them. Only the - blocks overlapping the requested range are mapped, sliced, and concatenated. + Each block lives at its own file offset, so a contiguous sample range can + straddle several of them. Parameters ---------- @@ -791,9 +785,8 @@ def _read_multi_block_chunk(self, fid, memmap_specs, channels, i_start, i_stop, ------- np.ndarray Signal data of shape (i_stop - i_start, len(channel_indexes)), dtype int16. - A range falling inside a single block is returned as a view, as it would be - from a segment that was never merged; otherwise the blocks are copied into - one array. + A range inside a single block is a view; otherwise the blocks are copied + into one array. """ total_samples = sum(spec["num_samples"] for spec in memmap_specs) if i_start is None: @@ -801,8 +794,7 @@ def _read_multi_block_chunk(self, fid, memmap_specs, channels, i_start, i_stop, if i_stop is None: i_stop = total_samples - # Map only the blocks the range touches. Slicing the rows keeps these as views, - # so nothing is read from disk until the data is actually used below. + # Row slicing keeps these as views, so nothing is read until used below block_views = [] block_start = 0 for spec in memmap_specs: @@ -825,9 +817,8 @@ def _read_multi_block_chunk(self, fid, memmap_specs, channels, i_start, i_stop, if len(block_views) == 1: return block_views[0][:, channel_indexes] - # Fill one output array rather than concatenating, so the blocks are never all - # held at once, and a range touching no block still gives a correctly shaped - # empty result rather than failing + # Fill one array rather than concatenating: the blocks are never all held at + # once, and a range touching no block still gives a correctly shaped result selected_channel_count = np.arange(channels)[channel_indexes].size out = np.empty((i_stop - i_start, selected_channel_count), dtype=memmap_specs[0]["dtype"]) position = 0 @@ -855,22 +846,15 @@ def _get_blackrock_timestamps(self, block_index, seg_index, i_start, i_stop, str block has a single scalar timestamp for its first sample. All subsequent samples within the block are interpolated as ``block_timestamp + sample_index_within_block / sampling_rate``, - assuming uniform spacing. Uniform spacing holds within a block - because that is what the format means by a block; it does not hold - across blocks, which are separated by a pause of recorded length. - A segment built with ``gap_tolerance_ms`` may merge several blocks, - so the interpolation is done per block and steps across each pause. + assuming uniform spacing. A segment may merge several blocks, so the + interpolation is per block and steps across the pause between them. - **FileSpec 2.1:** No timestamps are stored in the file. All samples are interpolated from ``t_start=0`` using the sampling rate. - The times returned here do not depend on ``gap_tolerance_ms``. That - parameter decides how samples are grouped into segments, not what time - a sample was recorded at. When a segment merges blocks across a pause, - its single ``t_start`` cannot express that pause, so the segment's own - ``t_start + sample_index / sampling_rate`` is wrong for every sample - after it while the times returned here stay correct. That gap is the - reason this method exists. + The times returned do not depend on ``gap_tolerance_ms``, which decides + how samples are grouped into segments rather than when they were + recorded. Parameters ---------- @@ -918,15 +902,12 @@ def _get_blackrock_timestamps(self, block_index, seg_index, i_start, i_stop, str block_timestamps = seg.get("block_timestamps") if block_timestamps is None or len(block_timestamps) == 1: - # One block, so the segment's own t_start is the block's recorded start - # and the whole stretch is uniform from it + # One block, so the segment's t_start is that block's recorded start t_start = self._sigs_t_starts[nsx_nb][seg_index] return t_start + np.arange(i_start, i_stop, dtype="float64") / sr - # Several blocks merged into this segment. The segment has a single t_start - # and cannot express the pauses between them, but each block recorded its own - # start, so use those: the times returned here step across each pause even - # though the signal's uniform time base does not. + # Several blocks merged: the segment's single t_start cannot express the + # pauses between them, but each block recorded its own start ts_res = float(self._nsx_basic_header[nsx_nb]["timestamp_resolution"]) timestamps = np.empty(i_stop - i_start, dtype="float64") block_start = 0 @@ -1473,13 +1454,8 @@ def _format_gap_report( """ Format a detailed gap report showing where timestamp discontinuities occur. - This renders a report and nothing else: every quantity is passed in already - resolved, so it makes no assumption about how a format stores its timing. A - reader whose file carries one timestamp per sample and a reader whose file - carries one per block describe their gaps in the same sample-level terms, and - each works out those terms whichever way suits its own storage. That keeps the - report usable across readers as the gap API spreads to other formats (see - python-neo issue #1773). + Every quantity is passed in already resolved, so this makes no assumption about + how a format stores its timing and stays usable across readers (see #1773). Parameters ---------- @@ -1488,11 +1464,11 @@ def _format_gap_report( gap_positions_seconds : np.ndarray Time of that sample, in seconds since the start of the recording. gap_durations_ms : np.ndarray - Length of each gap in milliseconds. Negative when the clock jumped backwards. + Length of each gap in milliseconds. stream_label : str Name of the stream the gaps were found in, e.g. "ns2". detection_threshold_ms : float - Threshold that was used to detect the gaps, in milliseconds. + Threshold used to detect the gaps, in milliseconds. Returns ------- @@ -1522,10 +1498,9 @@ def _segment_nsx_data(self, parsed_data_headers, nsx_nb): Segment NSX data based on timestamp gaps. Takes the parsed data headers returned by _parse_nsx_data() and dispatches to - the segmenter for the file's format. The dispatch keys off the structure of the - parsed headers rather than the declared spec version, because - _parse_nsx_data_v30_ptp() falls back to standard parsing for files that declare - PTP but do not use it. + the segmenter for the file's format. Dispatch keys off the structure of the + headers rather than the declared spec, because _parse_nsx_data_v30_ptp() falls + back to standard parsing for files that declare PTP but do not use it. Parameters ---------- @@ -1541,14 +1516,11 @@ def _segment_nsx_data(self, parsed_data_headers, nsx_nb): [ { "timestamp": scalar or None, - Raw timestamp ticks at the start of the segment. + "block_timestamps": list (standard format only), "nb_data_points": int, "header": int or None, "offset_to_data_block": None (deprecated but kept for compatibility), - "memmap_specs": list[dict], - One entry per data block making up the segment. Standard-format - segments hold several when consecutive blocks were merged; PTP - and v2.1 segments always hold exactly one. + "memmap_specs": list[dict], one per data block in the segment, "timestamps_memmap_kwargs": dict (PTP only), }, ... @@ -1567,27 +1539,14 @@ def _classify_gaps(self, time_differences, sampling_rate, nsx_nb, sample_indices """ Classify timestamp discontinuities into segment boundaries. - Shared by the standard and PTP segmenters. Both hand over the interval between - the timestamps of two consecutive samples, which is one sampling period when the - data is contiguous, along with a sample-level description of where each interval - sits. The PTP path reads those descriptions from the per-sample timestamps in the - file; the standard path works them out from its block timestamps, since within a - block the timing is uniform by definition. Either way the gaps are detected, and - reported, in the same terms. - - Two kinds of discontinuity are treated differently, because only one of them - leaves the caller anything to decide: - - - Forward gaps (a pause, or dropped samples) are detected with a fixed - threshold of 2 sampling periods, which sits above the jitter and rounding - noise floor of these files. Merging across one distorts the segment's uniform - time base while splitting gives more segments, and only the caller knows which - is acceptable, so with gap_tolerance_ms unset this raises and reports them. - - Backward jumps (the clock runs backwards between consecutive samples, which a - recording reset causes) have no duration to compare against a tolerance, so - there is no choice to offer: they always become segment boundaries. Raising - would pose a question with one possible answer, so they are warned about - instead, and left out of the gap report, which is about gaps. + Shared by the standard and PTP segmenters, which both pass the interval between + consecutive samples, one sampling period when the data is contiguous. + + Forward gaps are detected at 2 sampling periods, above the jitter and rounding + noise floor of these files, and raise when gap_tolerance_ms is unset because + merging across one or splitting at it is a choice only the caller can make. + Backward jumps have no duration to compare against a tolerance, so they always + become boundaries and are warned about rather than raised on. Parameters ---------- @@ -1615,8 +1574,7 @@ def _classify_gaps(self, time_differences, sampling_rate, nsx_nb, sample_indices forward_mask = time_differences > detection_threshold backward_mask = time_differences < 0 - # Backward jumps always split, so there is nothing to ask about. Say what was - # found and what was done about it, and carry on. + # Backward jumps always split, so there is nothing to ask about backward_indices = np.flatnonzero(backward_mask) if len(backward_indices) > 0: locations = ", ".join( @@ -1634,8 +1592,7 @@ def _classify_gaps(self, time_differences, sampling_rate, nsx_nb, sample_indices if len(forward_indices) == 0: return backward_indices - # Error by default on forward gaps - merging across one or splitting at it is a - # choice only the caller can make + # Error by default on forward gaps - the caller must choose if self.gap_tolerance_ms is None: gap_report = self._format_gap_report( gap_sample_indices=sample_indices[forward_indices], @@ -1679,11 +1636,10 @@ def _segment_nsx_v22_v30(self, parsed_data_headers, nsx_nb): """ Segment standard format data (v2.2, v2.3, v3.0) using block-level gap detection. - The file already breaks data into a new block at every pause, so block - boundaries are the candidate segment boundaries and gap detection decides which - of them are real. Blocks separated by less than gap_tolerance_ms are merged into - one segment, which then keeps one memmap spec per contributing block because - those blocks sit at different offsets in the file. + The file starts a new block at every pause, so block boundaries are the + candidate segment boundaries. Blocks separated by less than gap_tolerance_ms + merge into one segment, which keeps one memmap spec per block because they sit + at different offsets in the file. """ ts_res = float(self._nsx_basic_header[nsx_nb]["timestamp_resolution"]) sampling_rate = self._nsx_sampling_frequency[nsx_nb] @@ -1692,17 +1648,13 @@ def _segment_nsx_v22_v30(self, parsed_data_headers, nsx_nb): block_sizes = np.array([header["memmap_kwargs"]["num_samples"] for header in parsed_data_headers]) if len(parsed_data_headers) > 1: - # Interval between the last sample of a block and the first sample of the - # next one, so contiguous blocks give exactly one sampling period. This is - # the per-block equivalent of np.diff() over the PTP per-sample timestamps. + # Interval between a block's last sample and the next block's first, so + # contiguous blocks give exactly one sampling period block_last_sample_times = block_t_starts + (block_sizes - 1) / sampling_rate inter_block_intervals = block_t_starts[1:] - block_last_sample_times[:-1] - # Describe each candidate boundary by the sample before it, the same way the - # PTP path does. Within a block the timing is uniform by definition, so the - # last sample of block i is sample cumsum(sizes)[i] - 1 and its time follows - # from the block's own timestamp. Working these out from the block scalars - # gives exactly what per-sample timestamps would say, without building them. + # Describe each boundary by the sample before it, as the PTP path does. + # Derived from the block scalars rather than built per sample. last_sample_indices = np.cumsum(block_sizes)[:-1] - 1 last_sample_positions = block_last_sample_times[:-1] - block_t_starts[0] @@ -1759,8 +1711,7 @@ def _segment_nsx_ptp(self, parsed_data_headers, nsx_nb): time_differences = np.diff(timestamps_in_seconds) - # Every sample carries its own timestamp, so the sample before each interval and - # its time are read straight off the array + # Every sample carries its own timestamp, so read these off the array last_sample_indices = np.arange(len(time_differences)) last_sample_positions = timestamps_in_seconds[:-1] - timestamps_in_seconds[0] diff --git a/neo/test/rawiotest/test_blackrockrawio.py b/neo/test/rawiotest/test_blackrockrawio.py index 194beaaa4..0ea5daa11 100644 --- a/neo/test/rawiotest/test_blackrockrawio.py +++ b/neo/test/rawiotest/test_blackrockrawio.py @@ -372,6 +372,8 @@ def test_gap_tolerance_ms_parameter(self): with self.assertRaises(ValueError) as context: reader = BlackrockRawIO(filename=dirname, nsx_to_load=6) reader.parse_header() + # The backward jump is not a gap, so it is not among them + self.assertIn("Detected 2 timestamp gaps", str(context.exception)) # Test 2: A tolerance above every forward gap keeps the data either side of them # together, but the backward jump splits regardless of the tolerance @@ -391,13 +393,9 @@ def test_backward_jump_warns_and_splits_without_gap_tolerance_ms(self): """ A backward jump warns rather than raising, and splits whatever the tolerance says. - Raising exists so the caller can choose between merging across a gap and splitting - at it. A backward jump offers no such choice: it has no duration to compare against - gap_tolerance_ms, so it always becomes a segment boundary. Raising would pose a - question with one possible answer, so the reader reports it and carries on. - - The reset file's clock restarts, so its two blocks are separated by a backward jump - rather than a gap. + It has no duration to compare against gap_tolerance_ms, so there is no choice to + offer the caller. The reset file's clock restarts, so its two blocks are separated + by a backward jump rather than a gap. """ dirname = self.get_local_path("blackrock/segment/ResetCorrect/reset") @@ -416,37 +414,23 @@ def test_backward_jump_warns_and_splits_without_gap_tolerance_ms(self): reader.parse_header() self.assertEqual(reader.segment_count(0), 2) - def test_gap_report_covers_gaps_only(self): - """ - The gap report is about gaps, so a backward jump is not counted among them. - - The PTP file holds two forward gaps and one backward jump. Only the gaps are - something gap_tolerance_ms can act on, so only they are reported and counted. - """ - dirname = self.get_local_path("blackrock/blackrock_ptp_with_missing_samples/Hub1-NWBtestfile_neural_wspikes") - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - with self.assertRaises(ValueError) as context: - reader = BlackrockRawIO(filename=dirname, nsx_to_load=6) - reader.parse_header() - - # Two forward gaps, not three discontinuities - self.assertIn("Detected 2 timestamp gaps", str(context.exception)) - def test_gap_tolerance_ms_standard_format(self): """ Test gap_tolerance_ms on a standard-format (v2.2/2.3/3.0) file. - The pause_correct file has two data blocks separated by a ~27 s recording pause. - Unlike the PTP path, gaps here are detected between block timestamps rather than - between per-sample timestamps. + The pause_correct file has two data blocks separated by a ~27 s recording pause, + detected between block timestamps rather than per-sample ones. """ dirname = self.get_local_path("blackrock/segment/PauseCorrect/pause_correct") - # Default behavior (None) raises, same as for PTP files - with self.assertRaises(ValueError): + # Default behavior (None) raises, same as for PTP files. The gap follows sample + # 3999 at 3.999 s, rather than sitting at the start of the recording + with self.assertRaises(ValueError) as context: reader = BlackrockRawIO(filename=dirname, nsx_to_load=2) reader.parse_header() + self.assertIn("3,999", str(context.exception)) + self.assertIn("3.999000", str(context.exception)) + self.assertIn("27009.700", str(context.exception)) # A tolerance below the pause keeps one segment per block reader_strict = BlackrockRawIO(filename=dirname, nsx_to_load=2, gap_tolerance_ms=0) @@ -499,32 +483,12 @@ def test_blackrock_timestamps_do_not_depend_on_gap_tolerance(self): chunk = reader_merged._get_blackrock_timestamps(0, 0, i_start, i_stop, stream_index) np.testing.assert_array_equal(chunk, expected[i_start:i_stop]) - def test_gap_report_locates_gap_in_sample_terms(self): - """ - The gap report must locate gaps by sample, whatever the file stores. - - A standard-format file has one timestamp per block rather than one per sample, but - the report describes gaps the same way the PTP report does. In pause_correct the - first block holds 4000 samples at 1 kHz, so the gap follows sample 3999 at 3.999 s - rather than sitting at the start of the recording. - """ - dirname = self.get_local_path("blackrock/segment/PauseCorrect/pause_correct") - with self.assertRaises(ValueError) as context: - reader = BlackrockRawIO(filename=dirname, nsx_to_load=2) - reader.parse_header() - - report = str(context.exception) - self.assertIn("3,999", report) # last sample before the gap, not block index 0 - self.assertIn("3.999000", report) # when that sample was taken, not 0.000000 - self.assertIn("27009.700", report) # gap duration in ms - def test_gap_tolerance_ms_standard_format_merged_chunk(self): """ Reading across blocks merged into one segment must stitch them together. - A merged segment holds one memmap spec per block, so a requested sample range can - straddle several of them. Compare reads from the merged segment against the same - data read as two separate segments, which needs no stitching. + Compared against the same data read as two separate segments, which needs no + stitching. """ dirname = self.get_local_path("blackrock/segment/PauseCorrect/pause_correct") stream_index = 0