PyPREP interpolates more channels than MATLAB PREP during re-referencing

[a-hurst]

Another difference I caught trying to unit-test my code for #77: during re-referencing, after the initial pass that detects and removes bad-by-NaN/flat/low-SNR channels from the initial average reference, MATLAB PREP's rereferencing loop doesn't interpolate channels that are bad by SNR or dropout (including any bad-by-SNR channels that were flagged on the first pass).

Basically, in PyPREP we always update the list of total bad channels after each re-reference iteration using NoisyChannels.get_bads(), which returns all bad channel types. In MatPREP, it's updated by the following function:

function ref = updateBadChannels(ref, noisy)
% Update the bad channel lists from ref based on bad channels in noisy
    ref.badChannelsFromNaNs = union(ref.badChannelsFromNaNs, ...
        noisy.badChannelsFromNaNs);
    ref.badChannelsFromNoData = union(ref.badChannelsFromNoData, ...
        noisy.badChannelsFromNoData);
    ref.badChannelsFromHFNoise = union(ref.badChannelsFromHFNoise, ...
        noisy.badChannelsFromHFNoise);
    ref.badChannelsFromCorrelation = union(ref.badChannelsFromCorrelation, ...
        noisy.badChannelsFromCorrelation);
    ref.badChannelsFromDeviation = union(ref.badChannelsFromDeviation, ...
        noisy.badChannelsFromDeviation);
    ref.badChannelsFromRansac = union(ref.badChannelsFromRansac, ...
        noisy.badChannelsFromRansac);
    ref.badChannelsFromDropOuts =  union(ref.badChannelsFromDropOuts, ...
        noisy.badChannelsFromDropOuts);
    ref.all = union(...
              union(ref.badChannelsFromNaNs, ref.badChannelsFromNoData), ...
              union( ...
                union(ref.badChannelsFromHFNoise, ref.badChannelsFromCorrelation), ...
              union( ...  
              union(ref.badChannelsFromDeviation,ref.badChannelsFromRansac), ...
                    ref.badChannelsFromRansac)));

As you can see,

1. Bad-by-SNR channels aren't updated with each iteration of re-referencing (make sense, since they're just the union of bad-by-HF and bad-by-correlation channels).
2. The set of all bads is the combination of all updated bad-by-NaN/flat/deviation/HF-noise/correlation/RANSAC channels instead of a union between ref.all and noisy.all. This means that any bad-by-SNR channels flagged during the first pass don't get carried over into the ref.all for future iterations if they aren't flagged bad by some other metric.
3. Dropout channels are updated with each iteration but aren't included in ref.all, possibly by mistake (see: the duplicate bad-by-RANSAC at the end).

Looking back at the MatPREP commit where that function was introduced, it looks like this function replaces simple updating of the full set of noisy channels so at least some of it seems to be deliberate?

Anyway, I'm not entirely sure how to best fix this because I can't easily tell whether PyPREP's logic or MatPREP's makes the most sense, and thus whether this should be fully or partly wrapped in matlab_strict. 3 certainly seems to be undesirable, since you probably wouldn't want a channel with intermittent flat regions being included in your estimated clean reference signal. What do you think?

[sappelhoff]

It seems like there is no bottom to this pit! Thanks for digging on :)

Dropout channels are updated with each iteration but aren't included in ref.all, possibly by mistake (see: the duplicate bad-by-RANSAC at the end).

😱 agreed, it looks like a mistake. I think you should open an issue about that on MatPREP.

I'm not entirely sure how to best fix this because I can't easily tell whether PyPREP's logic or MatPREP's makes the most sense, and thus whether this should be fully or partly wrapped in matlab_strict

I am not sure either to be honest.

[a-hurst]

It seems like there is no bottom to this pit! Thanks for digging on :)

Hopefully apart from CleanLine this is the bottom: with this issue addressed along with #92, the final post-interpolation reference signal matches the one from the MatPREP artifacts within a tolerance of 1e-4 (think there's maybe a minor rounding difference somewhere in the interpolation code), and all of the pre-interpolation and post-interpolation noisy channels match up exactly 🎉

I think you should open an issue about that on MatPREP.

I'll bring it up whenever Dr. Robbins gets back to me about #54 (I pinged her about it a week or so ago, she said she's currently deep in another project but will respond as soon as she has a chance to pull herself out of that headspace). In the meantime, I guess we could wrap it in matlab_strict until the bug is fixed in MatPREP itself?

I am not sure either to be honest.

I guess it hinges on whether a channel originally flagged as bad-by-SNR is a channel we'd want to exclude from the estimated robust reference no matter what: if a channel is initially bad by both HF noise and low inter-channel correlation, but after initial average referencing is no longer bad by either, I can't think of an a-priori reason to permanently exclude it from the average reference while still leaving in any channels bad by RANSAC/deviation/correlation/noise/dropout detected during that first pass (see #78). The only way that SNR channels are different from those is that they're excluded from the initial reference signal along with bad-by-NaN/flat channels, which are always excluded due to NaN and rank deficiency issues (respectively).

[sappelhoff]

In the meantime, I guess we could wrap it in matlab_strict until the bug is fixed in MatPREP itself?

agreed 👍 regarding opening an issue: your choice of course, but I recommend to keep as much as possible on public channels instead of Emails (which only very few people can profit from/help with) :)

I can't think of an a-priori reason to permanently exclude it from the average reference while still leaving in any channels bad by RANSAC/deviation/correlation/noise/dropout detected during that first pass

that sounds reasonable to me.