Caching channel interpolation matrix

[christianbrodbeck]

Interpolating bad channels on multiple evoked objects can be quite slow, and the bottleneck is computation of the interpolation matrix. When the evoked objects are from the same raw file they ought to have the same bad channels, so the interpolation matrix is redundantly computed multiple times. How about exposing the function to compute the interpolation matrix and adding a parameter to interpolate_channels() that allows providing a precomputed interpolation matrix?

[agramfort]

why not. @jasmainak had similar thoughts for autoreject work. One option I suggested was to use joblib.Memory but there is maybe a simpler way... do you have an API in mind?

[christianbrodbeck]

I was thinking of just adding a public compute_bad_channel_interpolation_matrix(info) and adding a parameter to the .interpolate_bads() to provide this matrix instead of computing it on the fly.

This would then make it possible to use different methods for caching depending on the situation, i.e., an autoreject function might locally cache them like you suggest, or a user could write them to disk. If there is a way to uniquely identify the matrix we could also implement module-wide caching like compute_morph_matrix(), but I think the approach would have to differ for MEG where the channels are constant within system, and EEG where they differ for each subject, so it might be too complex...

[agramfort]

sounds like a plan

[larsoner]

I was thinking of just adding a public compute_bad_channel_interpolation_matrix(info) and adding a parameter to the .interpolate_bads() to provide this matrix instead of computing it on the fly.

This would be a nice first step. Then when the parameter is used in interpolate_bads, it can check to make sure it's operating on a compatible input (see below)

If there is a way to uniquely identify the matrix

IIRC it will depend on info['dev_head_t'], the chosen expansion origin, and which channels are being interpolated from/to, although we might be able to relax this last assumption depending on the relative speed of the "self dots" and "cross dots" that are used (I think that's what we call them in the code somewhere).

At some point we might want to try using multipole moments (what MaxFilter does) instead of forward/inverse computations, it might be considerably faster.

I would say let's not bother with automating the process for now, I agree it adds a (potentially much) more difficult layer of complexity as you suggest.

[larsoner]

This is also potentially relevant for #4491

[jasmainak]

Great initiative! In autoreject, I ended up copying the code from MNE and creating a private function called _compute_dots() which is cached using joblibs. See the relevant portion of the code here. If someone has ideas on exposing this kind of API in MNE, it would be great. That way, it would remove some of the burden of copy-pasting code from MNE for me.

But, you don't necessarily need to use joblibs. I guess some kind of dictionary based memoizing might be good enough even for autoreject.

[christianbrodbeck]

when the parameter is used in interpolate_bads, it can check to make sure it's operating on a compatible input (see below)

Would it make sense to create a new object to return instead of an array? That way we could store a few properties (like the list of bad channels) that would make it possible to do some sanity checks before applying the matrix. But it would make saving more complicated.

@jasmainak yeah I would suggest for now to expose the function to compute and apply the matrix, and leave the caching mechanism to the specific application.

[larsoner]

But it would make saving more complicated.

If the object subclasses dict it remains almost trivial with read_hdf5. Same idea as in #3292 basically.