search

mne-tools / mne-nirs

Process Near-Infrared Spectroscopy Data in MNE
https://mne.tools/mne-nirs/
BSD 3-Clause "New" or "Revised" License
79 stars 35 forks source link

Incorrect looping in scalp_coupling_index_windowed #530

Closed alexk101 closed 6 months ago

alexk101 commented 7 months ago

Describe the bug

This bug has to do with upstream changes in the MNE package that were never propagated to this project in the function.

https://github.com/mne-tools/mne-nirs/blob/3ba42cd47bb4f5e6bc6c15d4677bb2b2177cb6ae/mne_nirs/preprocessing/_scalp_coupling_segmented.py#L16

I believe that when this function was originally written ~3 years ago, that some assumptions were made about the order of channel names which no longer hold true, or perhaps never did, though I am not sure. Looking at a similar function which does this in the main MNE package at the overall channel level (ie not windowed), the method by which channels are looped over is quite different, with the non-windowed function looping this way

https://github.com/mne-tools/mne-python/blob/195a2cc8009160fd125e355b0280e903a941c874/mne/preprocessing/nirs/_scalp_coupling_index.py#L63

and the broken windowed function looping this way

https://github.com/mne-tools/mne-nirs/blob/3ba42cd47bb4f5e6bc6c15d4677bb2b2177cb6ae/mne_nirs/preprocessing/_scalp_coupling_segmented.py#L95

The reason I say that this is broken is because when using this function to calculate the windowed sci for my data, half of the matrix is empty. Looking at this initially, I thought that my channels may just be really bad, but looking at the signal by hand, I can see a very clear heartbeat, but a zero in every entry for that channel in sci. Looking even more closely, the channels which are empty are never visited during the previously mentioned loop.

I rewrote this function to be in line with the non-windowed sci function in the main mne repo, and it fully populates the matrix, though I am unsure of the correctness. I can open a pull request with that change if the maintainers would like.

Steps to reproduce

from mne_nirs.preprocessing import scalp_coupling_index_windowed
import mne

from mne.io import read_raw_snirf

target = '2024-01-19_007.snirf'
raw = read_raw_snirf(target, optode_frame='mri')
raw_od = mne.preprocessing.nirs.optical_density(raw)
_, scores, times = scalp_coupling_index_windowed(raw_od, time_window=60, h_freq=1.35, h_trans_bandwidth=0.1)
print(scores)

2024-01-19_007.snirf.zip

Expected results

The score matrix should be fully populated for all channels.

[[0.34803152 0.31609754 0.1840236  ... 0.27982484 0.27195905 0.29996399]
 [0.88531414 0.75531441 0.96895067 ... 0.80156894 0.96327637 0.63165717]
 [0.82527266 0.8338465  0.98065507 ... 0.71513703 0.95658578 0.82656547]
 ...
 [0.43396164 0.39283085 0.56580488 ... 0.45583541 0.37026547 0.30859768]
 [0.8878728  0.83802229 0.90493666 ... 0.80879638 0.82400417 0.83602326]
 [0.98034899 0.90249146 0.98654176 ... 0.96404636 0.97938728 0.97132779]]

Actual results

Half of the score matrix is 0, having never been visited during the traversal of channels in scalp_coupling_index_windowed

[[0.93074168 0.92672519 0.98359125 ... 0.87374969 0.98690225 0.87919764]
 [0.73833978 0.52241733 0.82550531 ... 0.61606484 0.76549395 0.54156114]
 [0.27467124 0.19207873 0.23785208 ... 0.31039999 0.41150754 0.37729682]
 ...
 [0.         0.         0.         ... 0.         0.         0.        ]
 [0.         0.         0.         ... 0.         0.         0.        ]
 [0.         0.         0.         ... 0.         0.         0.        ]]
alexk101 commented 7 months ago

Another note, this behavior is also present in the function

https://github.com/mne-tools/mne-nirs/blob/3ba42cd47bb4f5e6bc6c15d4677bb2b2177cb6ae/mne_nirs/preprocessing/_peak_power.py#L17

With half of the peak spectral power matrix being empty (0). The looping structure is the same as in the previously mentioned window sci function

https://github.com/mne-tools/mne-nirs/blob/3ba42cd47bb4f5e6bc6c15d4677bb2b2177cb6ae/mne_nirs/preprocessing/_peak_power.py#L93

larsoner commented 7 months ago

@alexk101 thanks for the thorough analysis, sounds like you're right. And also like our tests are not complete enough. A bunch of zeros appearing in the output should probably should have started failing!

@alexk101 are you up for making a PR to take a stab at fixing this?

alexk101 commented 7 months ago

@larsoner The changes have been made and are correct (afaik). CI is failing on some other unrelated things though.