BUG: handle temporal discontinuities in Neuralynx .ncs files

[KristijanArmeni]

Reference issue

Addresses #12247

What does this implement/fix?

Tries do deal with temporal gaps in .ncs files. Conceptually, read_raw_neuralynx() now checks for any temporal differences between consecutive neo.Block[0].segments start/stop times and if any of those are larger than the sampling period, call it a gap (i.e. missing > 1 sample). Then, infer the number of missing samples for each gap based on gap times, do the bookkeeping, generate arrays of 0's for each gap data segment and insert them between respective valid segments.

I can read-in a large dataset locally (and test_neuralyx.py passes), but now raw.plot() will choke (see more below).

Example behavior on a local dataset where largest discontinuity seems to be 30 msecs (sampling rate = 4kHz).

>>> raw = read_raw_neuralynx(dataset, exclude_fname_patterns=file_patterns, preload=False)
Checking files in /local/neuralynx/dataset/
Ignoring .ncs files:
[...]
N = 418 discontinuous Neo segments detected with delta > 0.00025 sec
(max = 0.030219078063964844 sec, min = 0.00031304359436035156)

Note, that this assumes that all segments returned by Neo are valid data. I tried asking about his in neo github, but haven't heard back.

Additional information: raw.plot() breaks

Essentially, calling .plot with some selection raw.plot(start=10, duration=2) will choke with the output below.

My hunch is the reloading of chunks for plotting somehow interacts with my implementation, where all_data array constructed internally from neo segments + times of inferred gaps, won't match the number of samples of the array generated as block = np.zeros(n_channels, stop-start).

Not sure how to best troubleshoot this (I can't share the data I use locally to reproduce and the testing dataset can't be used for this since it's continuous), but hopefully there's something obvious I'm missing?

Traceback (most recent call last):
  File "/home/kriarm/miniconda3/envs/mnedev/lib/python3.11/site-packages/mne_qt_browser/_pg_figure.py", line 2958, in run
    data_chunk, times_chunk = browser._load_data(start, stop)
                              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/kriarm/git/mne-python/mne/viz/_figure.py", line 321, in _load_data
    return self.mne.inst[:, start : stop + 2]
           ~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^
  File "/home/kriarm/git/mne-python/mne/io/base.py", line 853, in __getitem__
    return self._getitem(item)
           ^^^^^^^^^^^^^^^^^^^
  File "/home/kriarm/git/mne-python/mne/io/base.py", line 860, in _getitem
    data = self._read_segment(start=start, stop=stop, sel=sel)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "<decorator-gen-165>", line 12, in _read_segment
  File "/home/kriarm/git/mne-python/mne/io/base.py", line 478, in _read_segment
    _ReadSegmentFileProtector(self)._read_segment_file(
  File "/home/kriarm/git/mne-python/mne/io/base.py", line 2525, in _read_segment_file
    return self.__raw.__class__._read_segment_file(
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/kriarm/git/mne-python/mne/io/neuralynx/neuralynx.py", line 349, in _read_segment_file
    block[idx] = all_data  # shape = (n_channels, n_samples)
    ~~~~~^^^^^
ValueError: could not broadcast input array from shape (64,990915) into shape (64,990957)

[larsoner]

Not sure how to best troubleshoot this (I can't share the data I use locally to reproduce and the testing dataset can't be used for this since it's continuous), but hopefully there's something obvious I'm missing?

First, to cut raw.plot out of the equation, doing raw[:, 2 * int(round(raw.info["sfreq"]))] or so should also replicate the problem. Or maybe even just a raw.load_data() would do it.

Not sure how to best troubleshoot this (I can't share the data I use locally to reproduce and the testing dataset can't be used for this since it's continuous), but hopefully there's something obvious I'm missing?

I would make sure first that raw.load_data() works. This will force you to have a model for the complete set of data, and when something is wrong, it will be clear because you can think of things from sample 0 to sample N-1 of the entire recording. Once that works, you need to map from start, stop into that segment properly. The problem could be with the initial construction of the whole signal, or with the mapping of start, stop into a particular contiguous segment of the whole signal.

[KristijanArmeni]

Thanks! Indeed, raw.load_data() gives the same error (I guess I was testing with raw.get_data(), but could be I mixed up something). Anyway, it's more tractable to debug now. Will report back once there's some progress.

[KristijanArmeni]

@larsoner Fixed the bug on my end and added commits to detect gaps in Neo segments, infer missing samples, add them as 0's, and annotate as BAD_ACQ_SKIP.

Example gaps

For context, this is a ~40-min recording of story listening.

I cannot comment on whether these look like a typical Neuralynx ring buffer error or not since I have never eye-balled these before. But it seems they can be burst-like, interspersed with a very short (1 or 2 sample) valid samples (this is a 4Khz dataset, so timescale is stretched out to microsecond precision in plots). And these are very short, mostly microsecond-level gaps in the case of my datasets.

Creating annotations

Right now, I add the all information for annotations in _raw_extras and annotations can be created like so:

gap_annot = raw._raw_extras[0]["gap_annotations"]
annot = mne.Annotations(
    onset=gap_annot["onset"],
    duration=gap_annot["duration"],
    description=gap_annot["descriptions"],
)
raw.set_annotations(annot)

This requires the user to set annotations. Not sure, if it is preferred that annotations should be added under the hood.

[larsoner]

Yes I think the gaps should be automatically annotated. Then the user can decide what to do with them (including removing them if they want)

[KristijanArmeni]

@larsoner read_raw_neuralynx now automatically sets annotations if gaps are detected. See in a32ff97 how its done and if in good order.

Example call and feedback on a local dataset with gaps detected:

raw = read_raw_neuralynx("/path/to/local/dataset",
   ...:                  exclude_fname_patterns=glob_patterns,
   ...:                  preload=False)
Checking files in /path/to/local/dataset
Ignoring .ncs files:
CSC2.ncs
CSC3.ncs
... # prints all files
xBLT1.ncs
Checking for temporal discontinuities in Neo data segments.
N = 418 discontinuous Neo segments detected with delta > 0.00025 sec.
(max = 0.030219078063964844 sec, min = 0.00031304359436035156)
Annotating gaps as BAD_ACQ_SKIP.

In [2]: raw.annotations
Out[2]: <Annotations | 418 segments: BAD_ACQ_SKIP (418)>

[larsoner]

Looks good to me so far! One idea for a potential refactoring and a couple small cleanup ideas

[KristijanArmeni]

@larsoner updated datasets/config.py and added test in e3d4d.

It seems CIs are failing due to some neo import issues (via _soft_import or so?). Looking at the error messages, can't really figure why -- is it anything related to the test script/read_raw_neuralynx? (it passes locally)

Otherwise, a recap: Created 3 gaps (130 invalid samples) in 2 .ncs files, then read them in with read_raw_neuralynx, check that it inferred the correct number of missing samples and test the loaded data against the "ground truth" files (.mat loaded in via scipy).

Also added a line to devel.rst.

[larsoner]

This is one problem

Suggested change

	# ignore files with no gaps
	pytest.importorskip("neo")
	# ignore files with no gaps

but a worse problem is that it seems like neo is not being installed on any of those CIs. I'll push a commit

[larsoner]

Feel free to take out of draft mode if you're done then ping me @KristijanArmeni and I'll merge assuming CIs are green (other than the Windows pip-pre one, which is having unrelated issues)

[KristijanArmeni]

Feel free to take out of draft mode if you're done then ping me @KristijanArmeni

Thanks and done! @larsoner

and I'll merge assuming CIs are green (other than the Windows pip-pre one, which is having unrelated issues)

I pulled this latest branch locally and test_neuralynx.py passes fine. This is good for me as far as testing and implementation is concerned -- as far as I can understand the failing CIs are unrelated directly to this PR (but have sth to do with neo install failing?).

[larsoner]

Some are, for example:

  └─ neo   does not exist (perhaps a typo or a missing channel).

is because I changed environment.yml in this PR to add neo, but used the wrong name. After I push a commit to fix this, all CIs should come back green other than windows pip-pre.

[KristijanArmeni]

is because I changed environment.yml in this PR to add neo, but used the wrong name. After I push a commit to fix this, all CIs should come back green other than windows pip-pre.

Yeah, looks like the CIs are green, now @larsoner. Thanks! (This branch still needs a merge from main, but other than that it seems good)

[larsoner]

Updated to latest main and marking for merge when green, thanks in advance @KristijanArmeni !