hpos timepoints error

[drammock]

in trying to work around #964 (by passing --session a on the command line) I hit this unrelated error (happens for both subjects):

│14:42:54│ ❌ sub-215 ses-a run-noise A critical error occurred. The error message was: Head position time points must be greater than first sample offset, but found 43.3330 < 78.0000

Aborting pipeline run. The traceback is:

  File "/opt/conda/miniforge3/envs/badbaby/lib/python3.12/site-packages/mne_bids_pipeline/steps/preprocessing/_03_maxfilter.py", line 415, in run_maxwell_filter
    raw_sss = mne.preprocessing.maxwell_filter(raw, **mf_kws)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "<decorator-gen-423>", line 12, in maxwell_filter
  File "/opt/conda/miniforge3/envs/badbaby/lib/python3.12/site-packages/mne/preprocessing/maxwell.py", line 385, in maxwell_filter
    params = _prep_maxwell_filter(
             ^^^^^^^^^^^^^^^^^^^^^
  File "<decorator-gen-424>", line 12, in _prep_maxwell_filter
  File "/opt/conda/miniforge3/envs/badbaby/lib/python3.12/site-packages/mne/preprocessing/maxwell.py", line 472, in _prep_maxwell_filter
    head_pos = _check_pos(head_pos, head_frame, raw, st_fixed, raw.info["sfreq"])
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/conda/miniforge3/envs/badbaby/lib/python3.12/site-packages/mne/preprocessing/maxwell.py", line 1155, in _check_pos
    raise ValueError(

I don't know this dataset intimately yet (i.e. I haven't manually processed/plotted the files so I'm not sure what to expect; this pipeline run of only 2 subjs was meant as a way to do just that). Any suggestions what I should examine to see if it's a legitimate data problem vs a bug?

Here's the full pipeline output:

``` $ mne_bids_pipeline --config amtone.py --session a ┌────────┬ Welcome aboard MNE-BIDS-Pipeline! 👋 ────────────────────────────────────────────────────────────────────────────────────────────────── │14:35:21│ 📝 ses-a Using configuration: amtone.py │14:35:21│ ❌ Overriding config.sessions = ['a'] └────────┴ ┌────────┬ init/_01_init_derivatives_dir ───────────────────────────────────────────────────────────────────────────────────────────────────────── │14:35:21│ ⏳️ Initializing output directories. └────────┴ done (1s) ┌────────┬ init/_02_find_empty_room ────────────────────────────────────────────────────────────────────────────────────────────────────────────── └────────┴ done (1s) ┌────────┬ preprocessing/_01_data_quality ──────────────────────────────────────────────────────────────────────────────────────────────────────── │14:35:24│ ⏳️ sub-116 ses-a Reading experimental recording: sub-116_ses-a_task-AmplitudeModulatedTones │14:35:24│ ⏳️ sub-215 ses-a run-noise Reading empty-room recording: sub-emptyroom_ses-19130721_task-noise_meg.fif │14:35:24│ ⏳️ sub-215 ses-a Reading experimental recording: sub-215_ses-a_task-AmplitudeModulatedTones │14:35:24│ ⏳️ sub-116 ses-a run-noise Reading empty-room recording: sub-emptyroom_ses-19130718_task-noise_meg.fif │14:35:25│ ⏳️ sub-215 ses-a run-noise Finding flat channels and noisy channels using Maxwell filtering. │14:35:25│ ⏳️ sub-116 ses-a Finding flat channels and noisy channels using Maxwell filtering. │14:35:25│ ⏳️ sub-215 ses-a Finding flat channels and noisy channels using Maxwell filtering. │14:35:26│ ⏳️ sub-116 ses-a run-noise Finding flat channels and noisy channels using Maxwell filtering. │14:35:48│ ⏳️ sub-215 ses-a run-noise Found no flat channels. │14:35:48│ ⏳️ sub-215 ses-a run-noise Found 2 noisy channels: MEG2422, MEG2432 │14:35:48│ ⏳️ sub-215 ses-a run-noise Found 2 bad channels. │14:35:48│ ⏳️ sub-215 ses-a run-noise Initializing report HDF5 file │14:35:48│ ⏳️ sub-215 ses-a run-noise Adding original raw data to report │14:35:54│ ⏳️ sub-116 ses-a Found no flat channels. │14:35:54│ ⏳️ sub-116 ses-a Found 3 noisy channels: MEG1413, MEG2432, MEG2622 │14:35:54│ ⏳️ sub-116 ses-a Found 3 bad channels. │14:35:54│ ⏳️ sub-116 ses-a Initializing report HDF5 file │14:35:54│ ⏳️ sub-116 ses-a Adding original raw data to report │14:35:57│ ⏳️ sub-215 ses-a run-noise Adding noisy channel detection to report │14:36:01│ ⏳️ sub-215 ses-a run-noise Adding config and sys info to report [48/1526] │14:36:05│ ⏳️ sub-215 ses-a run-noise Saving report: /storage/badbaby-redux/bids-data/derivatives/mne-bids-pipeline/sub-215/ses-a/meg/sub-215_ses-a_task-AmplitudeModulatedTones_report.html │14:36:11│ ⏳️ sub-116 ses-a run-noise Found no flat channels. │14:36:11│ ⏳️ sub-116 ses-a run-noise Found 3 noisy channels: MEG1941, MEG2422, MEG2432 │14:36:11│ ⏳️ sub-116 ses-a run-noise Found 3 bad channels. │14:36:14│ ⏳️ sub-116 ses-a Adding noisy channel detection to report │14:36:14│ ⏳️ sub-215 ses-a Found no flat channels. │14:36:14│ ⏳️ sub-215 ses-a Found 5 noisy channels: MEG0633, MEG1042, MEG1041, MEG1112, MEG1123 │14:36:14│ ⏳️ sub-215 ses-a Found 5 bad channels. │14:36:14│ ⏳️ sub-215 ses-a Adding original raw data to report │14:36:17│ ⏳️ sub-116 ses-a Adding config and sys info to report │14:36:21│ ⏳️ sub-116 ses-a Saving report: /storage/badbaby-redux/bids-data/derivatives/mne-bids-pipeline/sub-116/ses-a/meg/sub-116_ses-a_task-AmplitudeModulatedTones_report.html │14:36:22│ ⏳️ sub-116 ses-a run-noise Adding original raw data to report │14:36:31│ ⏳️ sub-116 ses-a run-noise Adding noisy channel detection to report │14:36:34│ ⏳️ sub-215 ses-a Adding noisy channel detection to report │14:36:35│ ⏳️ sub-116 ses-a run-noise Adding config and sys info to report │14:36:38│ ⏳️ sub-215 ses-a Adding config and sys info to report │14:36:39│ ⏳️ sub-116 ses-a run-noise Saving report: /storage/badbaby-redux/bids-data/derivatives/mne-bids-pipeline/sub-116/ses-a/meg/sub-116_ses-a_task-AmplitudeModulatedTones_report.html │14:36:42│ ⏳️ sub-215 ses-a Saving report: /storage/badbaby-redux/bids-data/derivatives/mne-bids-pipeline/sub-215/ses-a/meg/sub-215_ses-a_task-AmplitudeModulatedTones_report.html └────────┴ done (1m 21s) ┌────────┬ preprocessing/_02_head_pos ──────────────────────────────────────────────────────────────────────────────────────────────────────────── │14:36:44│ ⏳️ sub-215 ses-a Marking 5 channels as bad. │14:36:44│ ⏳️ sub-215 ses-a Estimating cHPI amplitudes │14:36:44│ ⏳️ sub-116 ses-a Marking 3 channels as bad. │14:36:44│ ⏳️ sub-116 ses-a Estimating cHPI amplitudes │14:37:13│ ⏳️ sub-215 ses-a Estimating cHPI SNR │14:37:16│ ⏳️ sub-116 ses-a Estimating cHPI SNR │14:39:12│ ⏳️ sub-215 ses-a Estimating cHPI locations │14:39:18│ ⏳️ sub-116 ses-a Estimating cHPI locations │14:39:42│ ⏳️ sub-215 ses-a Estimating head positions │14:39:43│ ⏳️ sub-215 ses-a Adding cHPI SNR and head positions to report. │14:39:44│ ⏳️ sub-215 ses-a Saving report: /storage/badbaby-redux/bids-data/derivatives/mne-bids-pipeline/sub-215/ses-a/meg/sub-215_ses-a_task-AmplitudeModulatedTones_report.html │14:40:10│ ⏳️ sub-116 ses-a Estimating head positions │14:40:11│ ⏳️ sub-116 ses-a Adding cHPI SNR and head positions to report. │14:40:13│ ⏳️ sub-116 ses-a Saving report: /storage/badbaby-redux/bids-data/derivatives/mne-bids-pipeline/sub-116/ses-a/meg/sub-116_ses-a_task-AmplitudeModulatedTones_report.html └────────┴ done (3m 31s) ┌────────┬ preprocessing/_03_maxfilter ─────────────────────────────────────────────────────────────────────────────────────────────────────────── │14:40:13│ ⏳️ sub-116 ses-a Loading reference run: None. │14:40:13│ ⏳️ sub-215 ses-a Loading reference run: None. │14:40:13│ ⏳️ sub-116 ses-a Applying tSSS (10.0 sec, corr=0.98) with MC to experimental data │14:40:13│ ⏳️ sub-215 ses-a Applying tSSS (10.0 sec, corr=0.98) with MC to experimental data │14:40:14│ ⏳️ sub-215 ses-a Marking 5 channels as bad. │14:40:14│ ⏳️ sub-215 ses-a Destination is 0.0 mm and 0.0° from the original head position │14:40:14│ ⏳️ sub-116 ses-a Marking 3 channels as bad. │14:40:14│ ⏳️ sub-116 ses-a Destination is 0.0 mm and 0.0° from the original head position │14:41:41│ ⏳️ sub-215 ses-a Filtering cHPI │14:41:51│ ⏳️ sub-215 ses-a Writing sub-215/ses-a/meg/sub-215_ses-a_task-AmplitudeModulatedTones_proc-sss_raw.fif │14:41:54│ ⏳️ sub-215 ses-a Adding Maxwell filtered raw data to report. │14:42:07│ ⏳️ sub-215 ses-a Saving report: /storage/badbaby-redux/bids-data/derivatives/mne-bids-pipeline/sub-215/ses-a/meg/sub-215_ses-a_task-AmplitudeModulatedTones_report.html │14:42:25│ ⏳️ sub-116 ses-a Filtering cHPI │14:42:36│ ⏳️ sub-116 ses-a Writing sub-116/ses-a/meg/sub-116_ses-a_task-AmplitudeModulatedTones_proc-sss_raw.fif │14:42:39│ ⏳️ sub-116 ses-a Adding Maxwell filtered raw data to report. │14:42:51│ ⏳️ sub-116 ses-a Saving report: /storage/badbaby-redux/bids-data/derivatives/mne-bids-pipeline/sub-116/ses-a/meg/sub-116_ses-a_task-AmplitudeModulatedTones_report.html │14:42:52│ ⏳️ sub-215 ses-a run-noise Loading reference run: None. │14:42:52│ ⏳️ sub-116 ses-a run-noise Loading reference run: None. │14:42:52│ ⏳️ sub-215 ses-a run-noise Applying tSSS (10.0 sec, corr=0.98) with MC to empty-room data │14:42:52│ ⏳️ sub-116 ses-a run-noise Applying tSSS (10.0 sec, corr=0.98) with MC to empty-room data │14:42:54│ ⏳️ sub-215 ses-a run-noise Destination is 78.0 mm and 52.4° from the MEG device origin │14:42:54│ ❌ sub-215 ses-a run-noise A critical error occurred. The error message was: Head position time points must be greater than first sample offset, but found 43.3330 < 78.0000 Aborting pipeline run. The traceback is: File "/opt/conda/miniforge3/envs/badbaby/lib/python3.12/site-packages/mne_bids_pipeline/steps/preprocessing/_03_maxfilter.py", line 415, in run_maxwell_filter raw_sss = mne.preprocessing.maxwell_filter(raw, **mf_kws) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "", line 12, in maxwell_filter File "/opt/conda/miniforge3/envs/badbaby/lib/python3.12/site-packages/mne/preprocessing/maxwell.py", line 385, in maxwell_filter params = _prep_maxwell_filter( ^^^^^^^^^^^^^^^^^^^^^ File "", line 12, in _prep_maxwell_filter File "/opt/conda/miniforge3/envs/badbaby/lib/python3.12/site-packages/mne/preprocessing/maxwell.py", line 472, in _prep_maxwell_filter head_pos = _check_pos(head_pos, head_frame, raw, st_fixed, raw.info["sfreq"]) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/opt/conda/miniforge3/envs/badbaby/lib/python3.12/site-packages/mne/preprocessing/maxwell.py", line 1155, in _check_pos raise ValueError( │14:42:54│ ⏳️ sub-116 ses-a run-noise Destination is 96.3 mm and 38.4° from the MEG device origin │14:42:54│ ❌ sub-116 ses-a run-noise A critical error occurred. The error message was: Head position time points must be greater than first sample offset, but found 18.4130 < 33.0000 Aborting pipeline run. The traceback is: File "/opt/conda/miniforge3/envs/badbaby/lib/python3.12/site-packages/mne_bids_pipeline/steps/preprocessing/_03_maxfilter.py", line 415, in run_maxwell_filter raw_sss = mne.preprocessing.maxwell_filter(raw, **mf_kws) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "", line 12, in maxwell_filter File "/opt/conda/miniforge3/envs/badbaby/lib/python3.12/site-packages/mne/preprocessing/maxwell.py", line 385, in maxwell_filter params = _prep_maxwell_filter( ^^^^^^^^^^^^^^^^^^^^^ File "", line 12, in _prep_maxwell_filter File "/opt/conda/miniforge3/envs/badbaby/lib/python3.12/site-packages/mne/preprocessing/maxwell.py", line 472, in _prep_maxwell_filter head_pos = _check_pos(head_pos, head_frame, raw, st_fixed, raw.info["sfreq"]) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/opt/conda/miniforge3/envs/badbaby/lib/python3.12/site-packages/mne/preprocessing/maxwell.py", line 1155, in _check_pos raise ValueError( (badbaby) drmccloy@bieber:/storage/badbaby-redux/pipeline$ /opt/conda/miniforge3/envs/badbaby/lib/python3.12/multiprocessing/resource_tracker.py:254: UserWarning: resource_tracker: There appear to be 4 leaked semaphore objects to clean up at shutdown warnings.warn('resource_tracker: There appear to be %d ' ```

[larsoner]

First tip is to use on_error = 'debug' to get into a pdb prompt. In principle we should be using the empty room prep function to adjust the empty room file and/or head position file so that this conflict does not occur. The short version is that you want to process the empty room data the same way as the task data as much as possible, so you use the head positions from the task data and use movement compensation on the empty room data. I can't remember if the head position or the first_samp of the empty room data (or both) are supposed to be adjusted but it seems like this isn't happening properly, or we aren't passing an adjusted version or something.

[drammock]

First tip is to use on_error = 'debug' to get into a pdb prompt.

I've been using the --debug flag and it doesn't leave me in pdb prompt 😢 that was another thing I meant to mention / ask about, but forgot.

I'll see if I can get a dev install of MBP on the remote machine (should work, but it's kinda old and has limited disk space). Then I'll stick in some print statements and 1/0s to see what I can figure out

[larsoner]

I've been using the --debug flag and it doesn't leave me in pdb prompt 😢 that was another thing I meant to mention / ask about, but forgot.

On my Linux machine I added a 1/0 above line 415 raw_sss = mne.preprocessing.maxwell_filter(raw, **mf_kws) in _03_maxfilter.py then did:

$ mne_bids_pipeline mne_bids_pipeline/tests/configs/config_ds004229.py --debug
┌────────┬ Welcome aboard MNE-BIDS-Pipeline! 👋 ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
│14:02:53│ 📝 Using configuration: mne_bids_pipeline/tests/configs/config_ds004229.py
│14:02:53│ ❌ Overriding config.on_error = 'debug'
└────────┴ 
...
┌────────┬ preprocessing/_03_maxfilter ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
│14:04:13│ ⏳️ sub-102 run-noise Computing eSSS basis with 8 components
│14:04:14│ ⏳️ sub-102 run-noise Adding eSSS projectors to report.
│14:04:15│ ⏳️ sub-102 run-noise Saving report: /home/larsoner/mne_data/derivatives/mne-bids-pipeline/ds004229/sub-102/meg/sub-102_task-amnoise_report.html
│14:04:15│ ⏳️ sub-102 Loading reference run: None.
│14:04:15│ ⏳️ sub-102 Applying tSSS (10 sec, corr=0.98) with MC/eSSS to experimental data
│14:04:16│ ⏳️ sub-102 Marking 1 channel as bad.
│14:04:16│ ⏳️ sub-102 Destination is 17.1 mm and 19.2° from the original head position
│14:04:16│ 🐛 sub-102 A critical error occurred. The error message was: division by zero

Starting post-mortem debugger.
<traceback object at 0x71bd58e3f880>
> /home/larsoner/python/mne-bids-pipeline/mne_bids_pipeline/steps/preprocessing/_03_maxfilter.py(415)run_maxwell_filter()
-> 1/0
(Pdb) p raw
<Raw | sub-102_task-amnoise_meg.fif, 332 x 360001 (300.0 s), ~919.3 MB, data loaded>

and I tried modifying the config to have on_error = "debug" and ended up at the same prompt. Then I thought maybe having n_jobs=2 could break things but even there things work because it gets overridden:

┌────────┬ Welcome aboard MNE-BIDS-Pipeline! 👋 ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
│14:07:04│ 📝 Using configuration: mne_bids_pipeline/tests/configs/config_ds004229.py
│14:07:04│ ❌ Setting config.n_jobs=1 because config.on_error="debug"
└────────┴ 

So I'm not sure why it's not working for you :(

[drammock]

when running with an editable install of MBP, the --debug flag does indeed drop me into a postmortem debugger on error. Not sure why that wasn't working with a standard (conda-forge) install of 1.9.0

Summarizing the error: it happens when processing the emptyroom recording (which here has a raw._first_time of 33 seconds), and it loads derivatives/mne-bids-pipeline/sub-116/ses-a/meg/sub-116_ses-a_task-AmplitudeModulatedTones_headpos.txt (which is headpos from the task run) and finds that the first time of the headpos measurement is ~18 seconds. As @larsoner noted:

I can't remember if the head position or the first_samp of the empty room data (or both) are supposed to be adjusted but it seems like this isn't happening properly, or we aren't passing an adjusted version or something.

The problem turns out to be that raw was recorded at 1800 Hz, and empty room was recorded at 1000 Hz. So the adjustment of raw_er_prepared here:

https://github.com/mne-tools/mne-python/blob/43fb9d8c28193ad2c0e11c9282e65338e02f98b7/mne/preprocessing/maxwell.py#L184

does so without taking into account any possible difference in sfreq. So one solution would be to downsample all our raw files before starting the pipeline (since MNEBP does resampling after SSS). But it seems (to me at least, maybe I'm mistaken?) that it ought to be OK to have different sfreqs for the raw and the ER, and that we could make the prepare_emptyroom function account for such differences. WDYT @larsoner? Is having different sfreqs likely to cause other problems later?

cc @mdclarke

[larsoner]

I think having different sampling rates is okay in principle. Subsequent analyses should take them into account if necessary, but given that the most common use cases for empty-room data are computing a covariance or projectors -- neither of which should generally be affected by a sample rate difference -- I think it's okay to process it at a different sample rate.