Open leouieda opened 1 year ago
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.
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: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:
.grid
. It's so much simpler to just havefit
andpredict
and let users callverde.grid_coordinates
andverde.make_xarray_grid
. This handles every possible use case without us having to guess what people want to do.fit
be tricky. See below.predict
can take an extrafield
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.Proposal for field-specific equivalent sources
EquivalentSources
andEquivalentSourcesGB
. Make them not based onverde.base.BaseGridder
. Inherit from scikit-learn directly.fit
andpredict
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.fit
could be a singledata
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).predict
takescoordinates, field
only, both required.field="whatever we use for all 3 components"
and then we can call a separate functiontotal_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?).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 methodto_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 astotal_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:
Predicting the 3-field components correctly is tricky and @indiauppal is working on this as a research question.