Add gap_tolerance_ms API to Neuralynx reader

neo/io/neuralynxio.py

@@ -36,6 +36,8 @@ def __init__(
         include_filenames=None,
         exclude_filenames=None,
         keep_original_times=False,
+        gap_tolerance_ms=None,
+        strict_gap_mode=None,
         filename=None,
         exclude_filename=None,
     ):
@@ -67,6 +69,13 @@ def __init__(
             Preserve original time stamps as in data files. By default datasets are
             shifted to begin at t_start = 0*pq.second.
             Default: False
+        gap_tolerance_ms : float | None, default: None
+            Controls how timestamp gaps in NCS files are handled.
+            If None (default), a ValueError is raised when gaps are detected, with a
+            detailed gap report. If a float value is provided, gaps smaller than this
+            threshold (in milliseconds) are ignored, and gaps larger create new segments.
+        strict_gap_mode : bool | None, default: None
+            Deprecated. Use gap_tolerance_ms instead.
         """
 
         if filename is not None:
@@ -83,6 +92,8 @@ def __init__(
             include_filenames=include_filenames,
             exclude_filenames=exclude_filenames,
             keep_original_times=keep_original_times,
+            gap_tolerance_ms=gap_tolerance_ms,
+            strict_gap_mode=strict_gap_mode,
             use_cache=use_cache,
             cache_path=cache_path,
         )

neo/rawio/neuralynxrawio/ncssections.py

@@ -35,6 +35,7 @@
 """
 
 import math
+import warnings
 import numpy as np
 
 from enum import IntEnum, auto
@@ -67,6 +68,7 @@ def __init__(self):
         self.sects = []
         self.sampFreqUsed = 0  # actual sampling frequency of samples
         self.microsPerSampUsed = 0  # microseconds per sample
+        self.detected_gaps = []  # list of (record_index, gap_size_us) for all detected gaps
 
     def __eq__(self, other):
         samp_eq = self.sampFreqUsed == other.sampFreqUsed
@@ -216,13 +218,25 @@ def _buildNcsSections(ncsMemMap, sampFreq, gapTolerance=0):
                 n_samples=n_samples,
             )
             ncsSects.sects.append(section0)
+            # No gaps detected in fast path
+            ncsSects.detected_gaps = []
 
         else:
             # need to parse all data block to detect gaps
             # check when the predicted timestamp is outside the tolerance
             delta = (ncsMemMap["timestamp"][1:] - ncsMemMap["timestamp"][:-1]).astype(np.int64)
             delta_prediction = ((ncsMemMap["nb_valid"][:-1] / sampFreq) * 1e6).astype(np.int64)
 
+            # Always detect all gaps using the strict threshold for reporting
+            strict_tolerance = round(NcsSectionsFactory._maxGapSampFrac * 1e6 / sampFreq)
+            all_gap_inds = np.flatnonzero(np.abs(delta - delta_prediction) > strict_tolerance)
+            gap_sizes_us = (delta - delta_prediction)[all_gap_inds]
+            ncsSects.detected_gaps = [
+                (int(record_index + 1), int(gap_size))
+                for record_index, gap_size in zip(all_gap_inds, gap_sizes_us)
+            ]
+
+            # Use user-provided tolerance for actual segmentation
             gap_inds = np.flatnonzero(np.abs(delta - delta_prediction) > gapTolerance)
             gap_inds += 1
 
@@ -245,7 +259,7 @@ def _buildNcsSections(ncsMemMap, sampFreq, gapTolerance=0):
         return ncsSects
 
     @staticmethod
-    def build_for_ncs_file(ncsMemMap, nlxHdr, gapTolerance=None, strict_gap_mode=True):
+    def build_for_ncs_file(ncsMemMap, nlxHdr, gap_tolerance_us=None, **kwargs):
         """
         Build an NcsSections object for an NcsFile, given as a memmap and NlxHeader,
         handling gap detection appropriately given the file type as specified by the header.
@@ -256,11 +270,37 @@ def build_for_ncs_file(ncsMemMap, nlxHdr, gapTolerance=None, strict_gap_mode=Tru
             memory map of file
         nlxHdr:
             NlxHeader from corresponding file.
+        gap_tolerance_us : float | None, default: None
+            Gap tolerance in microseconds for segmentation.
 
         Returns
         -------
         An NcsSections corresponding to the provided ncsMemMap and nlxHdr
         """
+        # Handle deprecated parameters
+        gapTolerance = kwargs.pop("gapTolerance", None)
+        strict_gap_mode = kwargs.pop("strict_gap_mode", None)
+        if kwargs:
+            raise TypeError(f"Unexpected keyword arguments: {list(kwargs.keys())}")
+
+        if gapTolerance is not None:
+            warnings.warn(
+                "The `gapTolerance` parameter is deprecated and will be removed in version 0.16. "
+                "Use `gap_tolerance_us` instead.",
+                DeprecationWarning,
+                stacklevel=2,
+            )
+            if gap_tolerance_us is None:
+                gap_tolerance_us = gapTolerance
+
+        if strict_gap_mode is not None:
+            warnings.warn(
+                "The `strict_gap_mode` parameter is deprecated and will be removed in version 0.16. "
+                "Use `gap_tolerance_us` instead.",
+                DeprecationWarning,
+                stacklevel=2,
+            )
+
         acqType = nlxHdr.type_of_recording()
         freq = nlxHdr["sampling_rate"]
 
@@ -269,14 +309,13 @@ def build_for_ncs_file(ncsMemMap, nlxHdr, gapTolerance=None, strict_gap_mode=Tru
             # restriction arose from the sampling being based on a master 1 MHz clock.
             microsPerSampUsed = math.floor(NcsSectionsFactory.get_micros_per_samp_for_freq(freq))
             sampFreqUsed = NcsSectionsFactory.get_freq_for_micros_per_samp(microsPerSampUsed)
-            if gapTolerance is None:
-                if strict_gap_mode:
-                    # this is the old behavior, maybe we could put 0.9 sample interval no ?
-                    gapTolerance = 0
+            if gap_tolerance_us is None:
+                if strict_gap_mode is not None and not strict_gap_mode:
+                    gap_tolerance_us = 0
                 else:
-                    gapTolerance = 0
+                    gap_tolerance_us = 0
 
-            ncsSects = NcsSectionsFactory._buildNcsSections(ncsMemMap, sampFreqUsed, gapTolerance=gapTolerance)
+            ncsSects = NcsSectionsFactory._buildNcsSections(ncsMemMap, sampFreqUsed, gapTolerance=gap_tolerance_us)
             ncsSects.sampFreqUsed = sampFreqUsed
             ncsSects.microsPerSampUsed = microsPerSampUsed
 
@@ -288,17 +327,16 @@ def build_for_ncs_file(ncsMemMap, nlxHdr, gapTolerance=None, strict_gap_mode=Tru
             AcqType.RAWDATAFILE,
         ]:
             # digital lynx style with fractional frequency and micros per samp determined from block times
-            if gapTolerance is None:
-                if strict_gap_mode:
-                    # this is the old behavior
-                    gapTolerance = round(NcsSectionsFactory._maxGapSampFrac * 1e6 / freq)
+            if gap_tolerance_us is None:
+                if strict_gap_mode is not None and not strict_gap_mode:
+                    # quarter of packet size is tolerated
+                    gap_tolerance_us = round(0.25 * NcsSection._RECORD_SIZE * 1e6 / freq)
                 else:
-                    # quarter of paquet size is tolerate
-                    gapTolerance = round(0.25 * NcsSection._RECORD_SIZE * 1e6 / freq)
-            ncsSects = NcsSectionsFactory._buildNcsSections(ncsMemMap, freq, gapTolerance=gapTolerance)
+                    # default: strict detection (0.2 of a sample interval)
+                    gap_tolerance_us = round(NcsSectionsFactory._maxGapSampFrac * 1e6 / freq)
+            ncsSects = NcsSectionsFactory._buildNcsSections(ncsMemMap, freq, gapTolerance=gap_tolerance_us)
 
-            # take longer data block to compute reaal sampling rate
-            # ind_max = np.argmax([section.n_samples for section in ncsSects.sects])
+            # take longer data block to compute real sampling rate
             ind_max = np.argmax([section.endRec - section.startRec for section in ncsSects.sects])
             section = ncsSects.sects[ind_max]
             if section.endRec != section.startRec:
@@ -315,13 +353,13 @@ def build_for_ncs_file(ncsMemMap, nlxHdr, gapTolerance=None, strict_gap_mode=Tru
 
         elif acqType == AcqType.BML or acqType == AcqType.ATLAS:
             # BML & ATLAS style with fractional frequency and micros per samp
-            if strict_gap_mode:
-                # this is the old behavior, maybe we could put 0.9 sample interval no ?
-                gapTolerance = 0
-            else:
-                # quarter of paquet size is tolerate
-                gapTolerance = round(0.25 * NcsSection._RECORD_SIZE * 1e6 / freq)
-            ncsSects = NcsSectionsFactory._buildNcsSections(ncsMemMap, freq, gapTolerance=gapTolerance)
+            if gap_tolerance_us is None:
+                if strict_gap_mode is not None and not strict_gap_mode:
+                    # quarter of packet size is tolerated
+                    gap_tolerance_us = round(0.25 * NcsSection._RECORD_SIZE * 1e6 / freq)
+                else:
+                    gap_tolerance_us = 0
+            ncsSects = NcsSectionsFactory._buildNcsSections(ncsMemMap, freq, gapTolerance=gap_tolerance_us)
             ncsSects.sampFreqUsed = freq
             ncsSects.microsPerSampUsed = NcsSectionsFactory.get_micros_per_samp_for_freq(freq)
 

neo/rawio/neuralynxrawio/neuralynxrawio.py

@@ -91,12 +91,19 @@ class NeuralynxRawIO(BaseRawIO):
     keep_original_times: bool, default: False
         If True, keep original start time as in files,
         Otherwise set 0 of time to first time in dataset
-    strict_gap_mode: bool, default: True
-        Detect gaps using strict mode or not.
-          * strict_gap_mode = True then a gap is consider when timstamp difference between two
-            consequtive data packet is more than one sample interval.
-          * strict_gap_mode = False then a gap has an increased tolerance. Some new system with different clock need this option
-            otherwise, too many gaps are detected
+    gap_tolerance_ms : float | None, default: None
+        Controls how timestamp gaps in NCS files are handled.
+        If None (default), a ValueError is raised when gaps are detected, with a
+        detailed gap report showing the number, size, and location of each gap.
+        If a float value is provided, gaps smaller than this threshold (in milliseconds)
+        are ignored, and gaps larger than this threshold create new segments.
+        Use gap_tolerance_ms=0.0 to segment on all detected gaps.
+    strict_gap_mode : bool | None, default: None
+        .. deprecated::
+            Use ``gap_tolerance_ms`` instead. Will be removed in version 0.16.
+        If explicitly set, uses legacy gap detection behavior:
+        strict_gap_mode=True uses tight tolerance (0.2 sample intervals),
+        strict_gap_mode=False uses loose tolerance (quarter of 512-sample packet).
 
     Notes
     -----
@@ -157,7 +164,8 @@ def __init__(
         include_filenames=None,
         exclude_filenames=None,
         keep_original_times=False,
-        strict_gap_mode=True,
+        gap_tolerance_ms=None,
+        strict_gap_mode=None,
         filename=None,
         exclude_filename=None,
         **kargs,
@@ -196,11 +204,32 @@ def __init__(
         else:
             self.rawmode = "one-dir"
 
+        # Handle gap_tolerance_ms and deprecated strict_gap_mode
+        if strict_gap_mode is not None:
+            warnings.warn(
+                "`strict_gap_mode` is deprecated and will be removed in version 0.16. "
+                "Use `gap_tolerance_ms` instead to control gap handling. "
+                "See issue #1773 for details.",
+                DeprecationWarning,
+                stacklevel=2,
+            )
+            if gap_tolerance_ms is not None:
+                warnings.warn(
+                    "Both `gap_tolerance_ms` and `strict_gap_mode` were provided. "
+                    "`gap_tolerance_ms` takes precedence.",
+                    UserWarning,
+                    stacklevel=2,
+                )
+            self._use_legacy_gap_mode = gap_tolerance_ms is None
+        else:
+            self._use_legacy_gap_mode = False
+
         self.dirname = dirname
         self.include_filenames = include_filenames
         self.exclude_filenames = exclude_filenames
         self.keep_original_times = keep_original_times
-        self.strict_gap_mode = strict_gap_mode
+        self.gap_tolerance_ms = gap_tolerance_ms
+        self.strict_gap_mode = strict_gap_mode if strict_gap_mode is not None else True
         BaseRawIO.__init__(self, **kargs)
 
     def _source_name(self):
@@ -878,6 +907,76 @@ def _rescale_event_timestamp(self, event_timestamps, dtype, event_channel_index)
         event_times -= self.global_t_start
         return event_times
 
+    def _format_gap_report(self, detected_gaps, sampling_frequency, filename):
+        """
+        Format a detailed gap report showing where timestamp discontinuities occur.
+
+        Parameters
+        ----------
+        detected_gaps : list of (int, int)
+            List of (record_index, gap_size_us) tuples from NcsSections.detected_gaps.
+        sampling_frequency : float
+            Sampling frequency in Hz used for the file.
+        filename : str
+            Path to the NCS file for the report header.
+
+        Returns
+        -------
+        str
+            Formatted gap report with table.
+        """
+        gap_durations_ms = [abs(gap_size) / 1000.0 for _, gap_size in detected_gaps]
+        gap_positions_seconds = [
+            record_index * NcsSection._RECORD_SIZE / sampling_frequency
+            for record_index, _ in detected_gaps
+        ]
+
+        gap_detail_lines = [
+            f"| {record_index:>15,} | {pos:>21.6f} | {dur:>21.3f} |\n"
+            for (record_index, _), pos, dur in zip(detected_gaps, gap_positions_seconds, gap_durations_ms)
+        ]
+
+        return (
+            f"Gap Report for {os.path.basename(filename)}:\n"
+            f"Found {len(detected_gaps)} timestamp gaps "
+            f"(detection threshold: {NcsSectionsFactory._maxGapSampFrac} sample intervals)\n\n"
+            "Gap Details:\n"
+            "+-----------------+-----------------------+-----------------------+\n"
+            "| Record Index    | Record at (Seconds)   | Gap Size (ms)         |\n"
+            "+-----------------+-----------------------+-----------------------+\n"
+            + "".join(gap_detail_lines)
+            + "+-----------------+-----------------------+-----------------------+\n"
+        )
+
+    def _get_neuralynx_timestamps(self, block_index, seg_index, stream_index):
+        """
+        Return original NCS record timestamps for the first channel in a stream/segment.
+
+        These are the raw hardware timestamps from the NCS file records, one timestamp
+        per 512-sample record, in microseconds.
+
+        Parameters
+        ----------
+        block_index : int
+            Block index (always 0 for Neuralynx).
+        seg_index : int
+            Segment index.
+        stream_index : int
+            Stream index.
+
+        Returns
+        -------
+        np.ndarray
+            Timestamps in microseconds from the NCS records, one per 512-sample record.
+        """
+        stream_id = self.header["signal_streams"][stream_index]["id"]
+        stream_mask = self.header["signal_channels"]["stream_id"] == stream_id
+        channel = self.header["signal_channels"][stream_mask][0]
+        chan_uid = (channel["name"], channel["id"])
+
+        data = self._sigs_memmaps[seg_index][chan_uid]
+        return data["timestamp"].copy()
+
     def scan_stream_ncs_files(self, ncs_filenames):
         """
         Given a list of ncs files, read their basic structure.
@@ -906,6 +1005,17 @@ def scan_stream_ncs_files(self, ncs_filenames):
         if len(ncs_filenames) == 0:
             return None, None, None
 
+        # Determine gap tolerance in microseconds for NcsSectionsFactory
+        if self._use_legacy_gap_mode:
+            # Legacy mode: let NcsSectionsFactory use strict_gap_mode defaults
+            factory_kwargs = {"strict_gap_mode": self.strict_gap_mode}
+        elif self.gap_tolerance_ms is not None:
+            # New API: convert ms to us
+            factory_kwargs = {"gap_tolerance_us": self.gap_tolerance_ms * 1000.0}
+        else:
+            # New default: detect all gaps strictly, we'll error later if gaps found
+            factory_kwargs = {}
+
         # Build dictionary of chan_uid to associated NcsSections, memmap and NlxHeaders. Only
         # construct new NcsSections when it is different from that for the preceding file.
         chanSectMap = dict()
@@ -918,9 +1028,36 @@ def scan_stream_ncs_files(self, ncs_filenames):
             verify_sec_struct = NcsSectionsFactory._verifySectionsStructure
             if not chanSectMap or (not verify_sec_struct(data, chan_ncs_sections)):
                 chan_ncs_sections = NcsSectionsFactory.build_for_ncs_file(
-                    data, nlxHeader, strict_gap_mode=self.strict_gap_mode
+                    data, nlxHeader, **factory_kwargs
                 )
 
+                # Check for gaps and handle according to gap_tolerance_ms
+                if not self._use_legacy_gap_mode and chan_ncs_sections.detected_gaps:
+                    if self.gap_tolerance_ms is None:
+                        # Default mode: only error on gaps >= 1 sample period.
+                        # Sub-sample deviations (from timestamp rounding, clock jitter, etc.)
+                        # are not real gaps and should not produce false positives.
+                        one_sample_us = 1e6 / chan_ncs_sections.sampFreqUsed
+                        significant_gaps = [
+                            (record_index, gap_us)
+                            for record_index, gap_us in chan_ncs_sections.detected_gaps
+                            if abs(gap_us) >= one_sample_us
+                        ]
+                        if significant_gaps:
+                            gap_report = self._format_gap_report(
+                                significant_gaps,
+                                chan_ncs_sections.sampFreqUsed,
+                                ncs_filename,
+                            )
+                            raise ValueError(
+                                f"Detected {len(significant_gaps)} timestamp gaps "
+                                f"in {os.path.basename(ncs_filename)}.\n"
+                                f"{gap_report}\n"
+                                f"To load this data, provide the gap_tolerance_ms parameter to "
+                                f"automatically segment at gaps larger than the specified tolerance.\n"
+                                f"Example: NeuralynxRawIO(dirname=..., gap_tolerance_ms=1.0)"
+                            )
+
             # register file section structure for all contained channels
             for chan_uid in zip(nlxHeader["channel_names"], np.asarray(nlxHeader["channel_ids"], dtype=str)):
                 chanSectMap[chan_uid] = [chan_ncs_sections, nlxHeader, ncs_filename]