Add gap_tolerance_ms API to Neuralynx reader (#1822)

* Add gap segmentation API to neuralynx

* zach review

* zach review

---------

Co-authored-by: Zach McKenzie <92116279+zm711@users.noreply.github.com>

Author: h-mayorquin

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 and will be removed in version 0.16.0. 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,27 +270,60 @@ 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.0. "
+                "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.0. "
+                "Use `gap_tolerance_us` instead.",
+                DeprecationWarning,
+                stacklevel=2,
+            )
+
+        if gap_tolerance_us is not None and gap_tolerance_us < 0:
+            raise ValueError(f"`gap_tolerance_us` must be non-negative, got {gap_tolerance_us}")
+
         acqType = nlxHdr.type_of_recording()
         freq = nlxHdr["sampling_rate"]
 
+        # Deprecation shim for strict_gap_mode (the boolean predecessor of gap_tolerance_us).
+        # strict_gap_mode True/None already match the modern per-type defaults below; only
+        # strict_gap_mode=False differed, tolerating a quarter-packet gap (PRE4 was 0 either way).
+        # Translating that single case here keeps the per-type branches on gap_tolerance_us only.
+        # Remove this block when strict_gap_mode is dropped in v0.16.0.
+        if gap_tolerance_us is None and strict_gap_mode is not None and not strict_gap_mode and acqType != AcqType.PRE4:
+            gap_tolerance_us = round(0.25 * NcsSection._RECORD_SIZE * 1e6 / freq)
+
         if acqType == AcqType.PRE4:
             # Old Neuralynx style with truncated whole microseconds for actual sampling. This
             # 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
-                else:
-                    gapTolerance = 0
-
-            ncsSects = NcsSectionsFactory._buildNcsSections(ncsMemMap, sampFreqUsed, gapTolerance=gapTolerance)
+            if gap_tolerance_us is None:
+                gap_tolerance_us = 0
+
+            ncsSects = NcsSectionsFactory._buildNcsSections(ncsMemMap, sampFreqUsed, gapTolerance=gap_tolerance_us)
             ncsSects.sampFreqUsed = sampFreqUsed
             ncsSects.microsPerSampUsed = microsPerSampUsed
 
@@ -288,17 +335,12 @@ 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)
-                else:
-                    # quarter of paquet size is tolerate
-                    gapTolerance = round(0.25 * NcsSection._RECORD_SIZE * 1e6 / freq)
-            ncsSects = NcsSectionsFactory._buildNcsSections(ncsMemMap, freq, gapTolerance=gapTolerance)
-
-            # take longer data block to compute reaal sampling rate
-            # ind_max = np.argmax([section.n_samples for section in ncsSects.sects])
+            if gap_tolerance_us is None:
+                # 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 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 +357,9 @@ 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:
+                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)