Equivalent sources specifically for magnetic data

[leouieda]

The EquivalentSources class implements a simple method that uses 1/distance as the kernel function. It's also based on the Verde model of fit/predict/grid. This means that it works for any potential field with few extra parameters. But it has some restrictions:

1. It can only fit and predict one data type (whatever field it fits, that's what it predicts).
2. It can't handle magnetic data specific operations like reduction to the pole because it's not really modelling magnetic fields.
3. It can't fit multiple data types, for example gravity + gravity gradients or magnetic total-field anomaly and anomalous B field components (from SWARM for example).
4. It can't predict fields it didn't fit. For example, predict anomalous field components from fitting total-field anomaly data.
5. The grid method is a bit strange because it needs the extra upward component. So our classes aren't really compatible with Verde. This style also doesn't work if we want to implement the other restrictions above (how do we specify inc/dec for total field anomaly or the field type?).

Me and @indiauppal have been playing with implementations for an upcoming paper that fits total-field anomaly and predicts the 3-component anomalous magnetic field (basically undo the projection onto the main field directions). While coding and using these magnetic-specific equivalent sources, we found that:

1. The fit/predict/grid model doesn't generalize. The main problem is .grid. It's so much simpler to just have fit and predict and let users call verde.grid_coordinates and verde.make_xarray_grid. This handles every possible use case without us having to guess what people want to do.
2. Passing data to fit be tricky. See below.
3. Predicting should be relatively easy. predict can take an extra field argument that can take a string or list of strings of fields to predict. For magnetics, we never need to predict total field anomaly since that can be calculated separately by function that takes the 3 field components and the main field direction. field can default to something or just be mandatory.
4. Basically ditch all pretense that equivalent sources are based on the Verde gridder model. This won't work anyway and has been discussed in the past (can't find a link right now).

Proposal for field-specific equivalent sources

1. Simplify the code classes for EquivalentSources and EquivalentSourcesGB. Make them not based on verde.base.BaseGridder. Inherit from scikit-learn directly.
2. Only implement fit and predict as public methods (maybe a few more the gradient-boosted version if needed). For the current classes, the arguments of these methods remain the same.
3. Start with magnetic equivalent sources since mag is more complicated than grav.
4. Arguments for fit could be a single data argument. It would be a dictionary with keys=fields and values = tuple of (coordinates, data, weights). For total-field anomaly it would need (coordinates, data, weights, main_field_direction).
5. predict takes coordinates, field only, both required.
6. To predict total-field anomaly, use field="whatever we use for all 3 components" and then we can call a separate function total_field_anomaly(field_components, main_field_direction). This one extra function call but makes everything so much simpler because we don't have to deal with the extra field direction argument that is only used for this (did I mention that magnetics suck?).
7. To reduce an anomaly to the pole correctly, we need to use the correct magnetization direction for the sources when fitting. Then, we need to predict while changing the direction of the sources to the pole and setting the main field direction also to the pole. This is a problem because the source direction is a parameter of the class. Using set_params would mean that we lose the original class so probably not the best option. A way to do this would be to add a method to_pole that checks if the class is fitted (issue a warning/error that you should fit first) and returns a copy of the class (keeping references to estimated attributes) with the source direction changed to the pole. It would be used as total_field_anomaly(eqs.fit(...).to_pole().predict(...), direction_at_pole).

Whatever works for magnetics will also work for gravity. So that's why we should do this one first. Tests of the API would be:

1. An example fitting and predicting total-field anomaly using the defaults.
2. An example reducing an anomaly to the pole.

Predicting the 3-field components correctly is tricky and @indiauppal is working on this as a research question.

[santisoler]

We've been discussing this issue in the Fatiando Meeting of 2024-01-04.

We think a nice interface for the new class would be something like this:

        class EquivalentMagSources:
            def __init__(self, ...):
                ...
                
            def fit(self, tfa=None, b_ee=None, ...):
                ...

            def predict(self, coordinates, field):
                ...
            
        eqs = hm.EquivalentMagSources(
            depth=..., 
            inc_sources: float =..., 
            dec_sources: float =...,
        )
        eqs.fit(
            tfa=(coordinates, data, weights, inclination, declination),
            b_ee=(coordinates, data, weights),
            b_nn=(coordinates, data, weights),
        )
        b_ee = eqs.predict(coordiantes, data, field="b_ee")
        (b_e, b_n, b_u) = eqs.predict(coordiantes, data, field="b")
        tfa = hm.total_field_anomaly((b_e, b_n, b_u), inc, dec)

Basically, the fit method should take keyword arguments (one for each field), with tuples containing the coordinates of the observation points, the observed data, their weights (optionally), and in case we are working if tfa (total field anomaly) we should also pass the inclination of declination of the background field.

The predict method will ask the coordinates and the field that should be computed. This "field" should be one of the ones available in the forward modelling functions: b, b_e, b_n, b_u (see #447 for more details). The tfa is not going to be a valid field: users should predict the three magnetic components and then compute the total field anomaly from them.