Skip to content

Add gap_tolerance_ms API to Neuralynx reader #1822

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
h-mayorquin wants to merge 5 commits into
base: master
Choose a base branch
from
Open

Add gap_tolerance_ms API to Neuralynx reader #1822

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
0d3363d
Add gap segmentation API to neuralynx
h-mayorquin Feb 18, 2026
4484dd4
Merge branch 'master' into fix_oversegmentation_neuralynx
h-mayorquin Apr 22, 2026
a7a0f2d
Merge branch 'master' into fix_oversegmentation_neuralynx
zm711 Jul 10, 2026
e8724d9
zach review
h-mayorquin Jul 15, 2026
4ea24db
Merge branch 'master' into fix_oversegmentation_neuralynx
h-mayorquin Jul 15, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 11 additions & 0 deletions neo/io/neuralynxio.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
):
Expand Down Expand Up @@ -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:
Expand All @@ -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,
)
Expand Down
89 changes: 62 additions & 27 deletions neo/rawio/neuralynxrawio/ncssections.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@
"""

import math
import warnings
import numpy as np

from enum import IntEnum, auto
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand All @@ -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):

@zm711 zm711 Jul 9, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

above you start with milliseconds and now you switch to microseconds.

@h-mayorquin h-mayorquin Jul 15, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The microseconds isn't a choice I'm making here, it's the factory's existing convention. NcsSectionsFactory has always worked in microseconds because NCS record timestamps are natively integer microseconds, and its gap math was already in us before this PR (the old parameter was just an unlabeled gapTolerance). I'm respecting that internal unit and only renamed it to gap_tolerance_us to make it explicit. The single ms->us conversion happens once at the reader boundary, so the user never leaves milliseconds; only the internal layer that was always in us stays in us. This is basically to respect the other contributor way of thinking inside the internals.

"""
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.
Expand All @@ -256,27 +270,57 @@ 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.",

@zm711 zm711 Jul 10, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

again how will this get exposed if the top-level function is in milliseconds?

@h-mayorquin h-mayorquin Jul 15, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

gap_tolerance_us is the factory-level parameter, not the user-facing one. End users go through NeuralynxRawIO(gap_tolerance_ms=...); this deprecation message is only reachable by code that calls the factory's old gapTolerance directly, which is internal. Fair point on consistency though: I can make the message point end users at the reader's gap_tolerance_ms so there's no ambiguity about which unit belongs where.

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"]

# 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.
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

Expand All @@ -288,17 +332,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:
Expand All @@ -315,13 +354,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)

Expand Down
155 changes: 146 additions & 9 deletions neo/rawio/neuralynxrawio/neuralynxrawio.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
-----
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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"
)
Comment on lines +910 to +949

@zm711 zm711 Jul 10, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have a hard time reading this, but it looked fine in Intan so I trust you for this bit of code.


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()

@zm711 zm711 Jul 10, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How expensive is this?

@h-mayorquin h-mayorquin Jul 15, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Cheap. It copies only the timestamp column: one uint64 per 512-sample record, so roughly 1/512 of the data size (tens of KB for a typical file, not the sample array). The .copy() is deliberate so the returned array owns its memory and isn't a live view into the file memmap.


def scan_stream_ncs_files(self, ncs_filenames):
"""
Given a list of ncs files, read their basic structure.
Expand Down Expand Up @@ -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:

@zm711 zm711 Jul 10, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again I'm not sure if I like this strategy overall is one I like. Rather than decide between legacy vs modern why not just make "modern" handle the legacy case so that we can slowly delete all legacy based code. You have this as a private variable which is better than making it public, but I think the strategy I would prefer would be:

def (a, b, c)
 if c is not None:
    b= c if c is x else y

This way right at the top we handle the legacy case by converting the b variable based on the deprecated c. Then when the deprecation window is over we just delete c and the top little bit of code.

@h-mayorquin h-mayorquin Jul 15, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I consolidated the value translation into the factory (see above), but the reader-level _use_legacy_gap_mode flag can't be collapsed the same way, because it isn't only routing: it also guards a genuine behavior difference. In legacy mode (a user passed strict_gap_mode) the reader deliberately skips the new error-on-gaps default so deprecated code keeps its old behavior; in modern mode it applies error-by-default. Collapsing the flag would silently change what deprecated calls do. The reader also can't pre-translate strict_gap_mode into a numeric tolerance up front, because the value depends on per-file acquisition type and frequency that only the factory sees. So the reader forwards the deprecated parameter to the factory, where it's now handled in one place, and keeps the flag only to preserve legacy behavior.

# 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}

@zm711 zm711 Jul 10, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is what I'm trying to understand. Won't it be confusing for the user if they are requesting a gap in ms but then internal we do everything on us? If we give them an error in us they might be confused no?

@h-mayorquin h-mayorquin Jul 15, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The conversion here is ms -> us only to feed the factory. The user never sees us: the gap report is explicitly emitted in ms (gap size) and seconds (position in the recording), so an error is already in the unit they passed. I'll double-check every user-facing string stays in ms.

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()
Expand All @@ -919,9 +1029,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]
Expand Down
Loading
Loading