ENH: Add support for Harmonic Field correction

[georgeoneill]

Fixes #10780.

What does this implement/fix?

• Add support for Homogenous Field Correction and extending the harmonic basis set to higher orders.
  - Creates an FieldCorrector() class which can be extend to correct forward models (important for deeper sources).
• Generates a Projector() for each column of the basis set.

To do

1. Add unit test(s)
2. update OPM preprocessing tutorial to incorporate HFC (results will look way better I think).
3. ?

Sorry this has taken this long to get around to this, I've been spinning a lot of plates recently, but this initial PR should keep this in the forefront of my mind for a while!

[larsoner]

Could you run HFC (of different orders, e.g., maybe just 1, 2, and 3) using your MATLAB code on the mne-testing-data dataset we have for read_raw_fil? Then we can check the Python implementation to numerical precision :)

[georgeoneill]

$\ell=1$

Max absolute difference: 1.56248141e-17
Max relative difference: 5.95234384e-08

$\ell=2$

Max absolute difference: 1.56249363e-17
Max relative difference: 5.95987598e-08

$\ell=3$

Max absolute difference: 4.2325848e-15
Max relative difference: 551.81531484

Interesting that $\ell=3$ has such a large difference. The lower orders all pass the default np.testing.assert_allclose(). This requires some investigation.

[larsoner]

The lower orders all pass the default np.testing.assert_allclose(). This requires some investigation.

For this sort of thing I often boil it down to a single use case to test (sounds like you're already there!) and then

1. save ... in MATLAB at some suitable place in the execution
2. In Python insert a scipy.io.loadmat(...)['whatever'] and use the MATLAB array in place of the Python one
3. Repeat 1-2 starting very late in the process (e.g., here, using the pseudo-inverted matrix) where you get equivalence, iteratively save .../loadmating earlier in the process until you find the culprit.

You can do the same thing by manually inspecting variables but that's not always easier...

[georgeoneill]

Below are the results for the three orders when for varying amount of field correction computation done in MNE.

Top table is max absolute error, bottom is max relative error. The numbers correspoding to Basis & Projector are smaller than the values posted last week due to a way I save out the data involving fewer casting steps.

MAX ABS DIFFERENCE

Order	Basis & Projector	Projector Only	None
$\ell=1$	1.28277618e-24	1.04366929e-24	5.62224323e-25
$\ell=2$	9.919705e-25	8.1748709e-25	6.33310156e-25
$\ell=3$	4.23163525e-15	1.04690046e-24	5.68686671e-25

MAX REL DIFFERENCE

Order	Basis & Projector	Projector Only	None
$\ell=1$	1.14841825e-08	9.00876775e-09	3.62863866e-09
$\ell=2$	2.57761766e-08	2.47442535e-08	5.75180087e-09
$\ell=3$	551.81529938	1.06555737e-08	8.69080621e-09

Currently pointing to an unexpected divergence between _sss_basis and spm_opm_vslm.

[georgeoneill]

I'm going to see what happens if we tell MNE to assume these are point magnetometers. Will report back.

[larsoner]

Wait the max relative difference is now 1e-9? That seems fine to me! I mean ideally for double float we'd end up at something closer to 1e-12 but 1e-9 should be good enough (especially since we routinely save to disk in single precision which will lead to ~1e-6 relative errors anyway)

[larsoner]

I'm going to see what happens if we tell MNE to assume these are point magnetometers. Will report back.

This could definitely be it, so feel free. But as above I think 1e-9 relative error is already enough. We routinely check for equivalence with MNE-C inverses for example at around 1e-7 relative precision or so

[georgeoneill]

Wait the max relative difference is now 1e-9? That seems fine to me! I mean ideally for double float we'd end up at something closer to 1e-12 but 1e-9 should be good enough (especially since we routinely save to disk in single precision which will lead to ~1e-6 relative errors anyway)

For $\ell=3$ the max relative difference 1e-9 when you use the basis set from SPM, but as soon as you calculate the basis set in MNE the relative error becomes 1e3, which equates to femtoTesla level differences.

[larsoner]

but as soon as you calculate the basis set in MNE the relative error becomes 1e3, which equates to femtoTesla level differences.

Ahh okay great that you isolated it!

BTW there has always been a slight mismatch when regularizing the maxwell filtering (SSS) basis set between what MNE chooses and what MaxFilter by MEGIN chooses. I wonder if we have some slight lingering scale factor or sign bug that changes things slightly... let me know what you find.

[larsoner]

(Also in maxwell_filter we column normalize the bases so scale factor differences might not matter there, whereas they could for this function)

[georgeoneill]

Its the coil definition!

Accurate max rel error: 551.81529938

Normal max rel error: 576.00884121

Point Magnetometer max rel error: 1.88025786e-07

[larsoner]

@georgeoneill awesome!

A bit scary that the coil def matters this much. But in theory at least MNE has a better definition I think (?), so maybe we're good here.

If you want, in tests you could set the coil type manually to be a point magnetometer and test against some hard-coded MATLAB values that you either write in ASCII (ouch) or add to MNE-testing-data. Then you can have a nice satisfying rtol=1e-8 or so in tests.

[georgeoneill]

I'll have a think about how we proceed with this. I'm not entirely sure how well definied the accurate coil definitions for QuSpin sensors actually are. There's the further complication for the triaxial sensors that there is a second laser used to generate the third axis so the 'accurate' coil definition may not entirely reflect the volume we should be integrating over.

I'll probably add a custom point magnetometer function to go with the accurate coil definition found in _prep_mf_coils (that we we dont need to touch mne.forward) with a flag which defaults to point magnetometer behaviour... for now.

[georgeoneill]

Still have the tests to make but wanted to at least get a near-final HFC script into the cloud!

[larsoner]

Awesome, looks very close now! Feel free to implement my (hopefully) final suggestions or, if you're sick of iterating, I can quickly push 😄

[larsoner]

Great! This might be what you actually end up wanting to do in practice with your data, but for now it allows us to decide the right course of action. If it is something you'll do a lot, we can do something like mne.preprocessing.mark_invalid_position_bad(inst) that looks through all data channels (which would include MEG) and looks at ch['loc'] and marks ones as bad that don't have valid positions 👍

[georgeoneill]

Hi @larsoner I think I agree with/understand most of the motivation. Happy for you to push those changes as I think you'll be able to get them to stick with a bit more precision than I can.

[larsoner]

@georgeoneill sorry for the slow response here, I plan to tackle this tomorrow

[georgeoneill]

No worries, no rush? Unless we need to get this in for the 1.4 deadline and I can have a go tonight?

[larsoner]

It would be nice to have for 1.4 but we're going to push that release date/deadline at least a week I think so. I'll just adjust it on GitHub now -- we have a lot of PRs/issues marked that we need to sort through.

[georgeoneill]

Okay, if you need me to step back in to make the revisions, let me know.

[larsoner]

Okay @georgeoneill I pushed some simplifying commits. Can you take a look at the diff now and see if you're happy with my changes? Hopefully you can see most (all?) places I changed simplified some stuff

[georgeoneill]

Hi @larsoner, this looks excellent. I've just renamed one of the test functions. I think its ready to go!

[larsoner]

Okay marked for auto-merge when green, thanks for working through this one @georgeoneill !

[georgeoneill]

Thank you for the assistance on this one.