Conversation
Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## main #1779 +/- ##
==========================================
- Coverage 81.75% 81.74% -0.02%
==========================================
Files 422 422
Lines 55756 56056 +300
Branches 5294 5345 +51
==========================================
+ Hits 45585 45824 +239
- Misses 8758 8797 +39
- Partials 1413 1435 +22 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
santisoler
left a comment
santisoler
left a comment
There was a problem hiding this comment.
Thanks for putting this together, @dccowan!
I just pushed some changes to the branch, mainly improving the docstring of the new receivers, and removing some uneeded noqa comments.
I left a few comments below. Let me know what do you think and if I could help with anything. I'm leaving some more general comments here too.
Breaking change in ApparentConductivity?
I was wondering if you checked that the ApparentConductivity in this PR (with the default arguments) generates the same data as in main.
Since now we are introducing an option to change how the data is computed, I'm wondering if the cross_product_amplitude method returns the same values as before.
I would strongly suggest to avoid changing those values since that would introduce a breaking change: users' code using ApparentConductivity would generate different values after updating.
Repeated mathematical operations
I noticed that we are repeating very similar mathematical operations, like computing the determinant of some matrices. I think we could create some private functions to handle those and avoid potential typos every time we repeat them. We can also add unit tests for those functions to check that individually they work as intended.
Structure
As we talked in last meeting, I'm not a big fan of the current structure of the classes. The _BaseOrientationInvariant is handling too many things, while the other receiver classes just work as dispatchers. This makes the _BaseOrientationInvariant too complex, and moves the implementation of each one of the receivers to a parent class.
I understand the spirit of the idea: since these methods are going to be used by each one of the new receiver classes (RootGramDeterminant, CrossProductAmplitude, and HorizontalDeterminant) and also by the ApparentConductivity class, it makes sense not to repeat things.
But what if we have private functions to compute each one of these quantities, and those functions can get reused by each one of those classes.
For example, we could have a function like:
def _get_root_gram_determinant(receiver, src, mesh, f):
if mesh.dim < 3:
raise NotImplementedError(
"'RootGramDeterminant' transfer function only for 3D simulation."
)
h = f[src, "h"]
hx = receiver.getP(mesh, "Fx", 0) @ h
hy = receiver.getP(mesh, "Fy", 0) @ h
hz = receiver.getP(mesh, "Fz", 0) @ h
if receiver.base_type == "magnetic":
bx = receiver.getP(mesh, "Fx", 1) @ h
by = receiver.getP(mesh, "Fy", 1) @ h
else:
e = f[src, "e"]
bx = receiver.getP(mesh, "Ex", 1) @ e
by = receiver.getP(mesh, "Ey", 1) @ e
# abs(det(H H*))
top = (
(np.abs(hx[:, 0] ** 2) + np.abs(hy[:, 0] ** 2) + np.abs(hz[:, 0] ** 2))
* (np.abs(hx[:, 1] ** 2) + np.abs(hy[:, 1] ** 2) + np.abs(hz[:, 1] ** 2))
) - np.abs(
hx[:, 0] * hx[:, 1].conjugate()
+ hy[:, 0] * hy[:, 1].conjugate()
+ hz[:, 0] * hz[:, 1].conjugate()
) ** 2
# abs(det(B B*)) = abs(det(B))**2
bot = np.abs(bx[:, 0] * by[:, 1] - bx[:, 1] * by[:, 0]) ** 2
return np.sqrt(top / bot)That gets used by RootGramDeterminant:
class RootGramDeterminant(...):
def eval(self, src, mesh, f):
return _get_root_gram_determinant(self, src, mesh, f)And also by ApparentConductivity:
class ApparentConductivity(...):
def eval(self, src, mesh, f): # noqa: A003 D102
if self._component == "root_gram_determinant":
root_gram_determinant = self._get_root_gram_determinant(self, src, mesh, f)
return _alpha(src) ** -1 * root_gram_determinant
...We could also keep those private functions simpler: we can remove those if statements that perform action if the receiver is an ApparentConductivity, and let it handle those extra steps.
Nonetheless, these changes would be applied to the internal structure of how these classes are implemented and wouldn't affect the public interface. So we can leave them here as ideas to implement later if we are happy with the public API of these new receivers, and the changes to the ApparentConductivity one.
|
Thanks so much @dccowan and @santisoler for your work on this! Just adding a quick comment on the apparent conductivity computation - I would be fine with us introducing a breaking change here. The previous version computed |H|/|E|, which is a quantity that cannot be computed from timeseries data in MT. So it was an incorrect value. The number of users who have used this so far is quite small because we haven't had any examples online yet, so we can put a big warning in the release notes, and if we wish, put a warning in receiver call, but I would consider the previous implementation to have been a bug in our understanding of the problem |
|
@santisoler regarding breaking change to ApparentConductivity: The way that the apparent conductivity datum is being generated in main is currently incorrect. While mathematically it works, it is not possible to generate the datum in this way in practice. The discrepancy between the apparent conductivity datum in main and the new ones is actually very small so I doubt anyone would notice. But since the way the apparent conductivity datum computed on the main branch is incorrect, I want to remove it as soon as possible. |
|
Just to confirm then general reorganization of functionality:
|
|
Thanks for clarifying the breaking change. Sounds good then! I would add an admonition to the .. admonition:: Breaking change
Since SimPEG v0.XX.0 the :class:`~simpeg.electromagnetics.natural_source.receivers.ApparentConductivity` receiver computes the apparent conductivity as [...]We could use SimPEG v0.26.0 since we are planning to merge this before the next release. I would also explicit such breaking change and the bug it fixed in the description of the PR, so it'll be included in the commit message, and we can include it in the changelog once we make the release.
If you are down with refactoring the code in this PR, please go ahead! Yes, my suggestion was that you could get rid of the simpeg/simpeg/electromagnetics/natural_source/receivers.py Lines 1258 to 1259 in b15d6fe I'm not 100% sure if this is possible since I haven't dug deep in the math, but if it's possible to apply those changes afterwards, I think it would make the code cleaner.
I wouldn't do that, since not necessarily every natural source receiver has a base station. For example the I wouldn't mind having the |
3 participants
Summary
This PR is intended to add the complete suite of rotation invariant transfer transfer functions to our receiver classes. These are primarily used for airborne NSEM systems but that can be used for general survey configurations.
PR Checklist
expect style.
to a Pull Request
@simpeg/simpeg-developerswhen ready for review.Supplementary math
Consider an AirMT system that measures three-component airborne magnetic fields. For a magnetic base station, we can generate the fundamental set of transfer functions$\mathbf{T}$ , where $\mathbf{T}$ is a 2x3 matrix. Likewise, for an electric base station, we can generate the fundamental admittance tensor $\mathbf{Y}$ . A multitude of orientation invariant transfer functions can be derived from the fundamental set. These are helpful when receiver orientation cannot be precisely resolved.
Root Gram Determinant:$\widehat{\mathbf{T}} = \sqrt{\mathbf{TT^\dagger}}$ where $\dagger$ represents the Hermitian. Similarly, the root Gram determinant can be generated for an electric base station $\widehat{\mathbf{Y}} = \sqrt{\mathbf{YY^\dagger}}$
Cross Product Amplitude: Orientation invariant transfer functions$|T|$ and $|Y|$ can be generated by taking the amplitude of the cross-product of the two rows of $\mathbf{T}$ and $\mathbf{Y}$ , respectively. See Mackie and Soyer (2026).
Horizontal Determinant: Taking the determinant of the horizontal transfer functions provides something invariant to z-axis rotation$det(T_H) = T_{xx} T_{yy} - T_{xy} T_{yx}$ and $det(Y_H) = Y_{xx} Y_{yy} - Y_{xy} Y_{yx}$ .
Apparent Conductivity: For an electric base station, all of these approaches can be used to generate apparent conductivity data.
This PR adds all the aforementioned transfer functions.
What is being added
A
_BaseOrientationInvariantwill host the underlying functionality for taking the root Gram determinant, cross product amplitude and horizontal determinant. There will be a property that sets whether an electric or magnetic base station is used.From this, I created child classes
RootGramDeterminant,CrossProductAmplitudeandHorizontalDeterminantthat users can interact with.There is also an
ApparentConductivityclass. Here, the component property is used to set the approach used to generate the apparent conductivity datum