-
Notifications
You must be signed in to change notification settings - Fork 281
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.
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
base: master
Are you sure you want to change the base?
Changes from all commits
0d3363d
4484dd4
a7a0f2d
e8724d9
4ea24db
3e661eb
09b7dc3
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -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.", | ||
|
Contributor
There was a problem hiding this comment. Choose a reason for hiding this commentThe 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?
Contributor
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
|
||
| 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) | ||
|
Contributor
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Do we need to control the round precision or just let it do what it wants with float?
Contributor
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. That line is verbatim from master's strict_gap_mode=False branch, so keeping round() reproduces the deprecated behavior exactly. Precision isn't really at stake: at 30 kHz it's 4266.67 against 4267 us, and the comparison is on int64 microseconds, so rounding can only flip the verdict for that one integer. |
||
|
|
||
| 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: | ||
|
Contributor
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Should we have some protection to ensure gap_tolerance_us is positive? What if someone puts in -10 will that error or pass silently with nonsense?
Contributor
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I think this makes sense and I added the assertion and a test. |
||
| 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) | ||
|
|
||
|
|
||
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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.
NcsSectionsFactoryhas 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 unlabeledgapTolerance). I'm respecting that internal unit and only renamed it togap_tolerance_usto 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.