Ensure times at which chpi are fitted are valid samples

[wmvanvliet]

I encountered this while maxfiltering 40 subjects. Some subjects failed, claiming Head position time points must be greater than first sample offset, but found 500.0005 < 500. A bit of a misleading error, as it was not the first time point that caused the problem, but the final time point.

When the continuous head position is fitted through compute_chpi_amplitudes, the first column of the result contains the times at which each fit was performed. However, these times were not necessarily aligned on the times of the original raw object, meaning some times fell in between two samples. This caused problems later on in maxwell_filter as a check is done whether the times are all valid, as the final timestamp can be slightly larger than the final timestamp of the raw object.

In this PR, compute_chpi_amplitudes is changed to always produce times that are aligned with those in the original raw object.

[larsoner]

the final timestamp can be slightly larger than the final timestamp of the raw object... In this PR, compute_chpi_amplitudes is changed to always produce times that are aligned with those in the original raw object.

I actually don't mind the fractional samples too much as it accurately reflects the center of the window (except at the edges, but this is by design I think). Instead the problem is maybe that fit_idxs can be chosen such that the center of a window is outside the bounds of the object. In other words, the fix should probably be in the calculation of fit_idxs rather than in rounding later.

WDYT? If you agree I feel free to give it a shot or I can push a commit to try it

[wmvanvliet]

I simplified a bit in my explanation.

The errors occurs in mne.utils.numerics._time_mask, which does not get the raw object passed into it, but just sfreq. Since the timestamps are not aligned on samples, in some cases, the include_tmax=True behavior can fail. This behavior is programmed as:

sfreq = 1000
tmax = times[-1]  # 810.6805
tmax = int(round(tmax * sfreq)) / sfreq  # 810.68
tmax += (0.5 if include_tmax else -0.5) / sfreq  # 810.6804999999999999, uhoh! 

maybe round should become ceil?

[larsoner]

The errors occurs in mne.utils.numerics._time_mask

To me this is just a symptom of the problem. We shouldn't have any time in quats[:, 0] that is greater than raw.times[-1]. It means that we centered a window at a point that is beyond the limit of the raw data. So I think this still points to the problem being that we had one too many points in fit_idxs, probably because we keep mixing float and int arithmetic all over the place in that function. I'd rather see if it's possible to stick with int arithmetic based on n_window and len(raw.times) until the very end and compute times directly from from fit_idxs (maybe with a +0.5 if we have an even number of window samples -- but the inte-based fit_idxs computation could account for adjusting stop in an arange call if that's the case).

[larsoner]

... or we just use floor arithmetic for the window length, then for even-length windows instead of getting 500.5 you get 500 and instead of going beyond len(raw.times) - 1 by 0.5 you land right at len(raw.times) - 1. That to me seems like the easiest solution. Then the fit_idxs that includes the last time point means you've used exactly half the window in your fit (which is what we do at the beginning, too, I think).

[larsoner]

Clearly I am a bit conflicted on the right approach here as I am contradicting myself in these discussions...

I think I'm happy with merging this as is, and at some point I might try to simplify the calculations so that there are fewer rounds and casts going on. If we can stick with ints longer I think the behavior and code will be clearer.

[wmvanvliet]

I agree and will see if I can fix it higher up.

[wmvanvliet]

This can actually not be fixed higher up. This is due to rounding errors. Perhaps we should just merge this because it works?

[larsoner]

Okay in it goes, someday maybe I'll try to simplify the logic but I don't think it should matter if we have to change by a sample or two. Thanks @wmvanvliet !