Add P-amplitude and S/P-amplitude-ratio inversion support#348
Open
MHkhosravi wants to merge 3 commits into
Open
Add P-amplitude and S/P-amplitude-ratio inversion support#348MHkhosravi wants to merge 3 commits into
MHkhosravi wants to merge 3 commits into
Conversation
- Add PolarityPampSPampRatioMisfit for joint polarity, P-amplitude, and S/P-amplitude-ratio objectives - Add multicomponent grid search with per-component weighting, optional normalization, and component result output - Add example for waveform, polarity, P-amplitude, and S/P-ratio inversion - Add beachball/data plotting support for polarity-amplitude results - Include mtuq subpackages in wheel packaging - Add focused tests for polarity-amplitude helpers
There was a problem hiding this comment.
Pull request overview
This PR extends MTUQ’s inversion workflow to optionally include P-wave amplitude and S/P amplitude-ratio information alongside existing waveform and polarity inversions, adding new misfit/grid-search utilities and plotting helpers to support joint analyses.
Changes:
- Added
PolarityPampSPampRatioMisfit(polarity + P-amplitude + log10(S/P)) with dictionary-to-observation helpers and unit tests. - Added
grid_search_multicomponentplus batching/progress-bar improvements to support multi-term misfit evaluation and combination. - Expanded plotting utilities (beachballs and lune plots) and added a new end-to-end example script.
Reviewed changes
Copilot reviewed 15 out of 17 changed files in this pull request and generated 12 comments.
Show a summary per file
| File | Description |
|---|---|
mtuq/misfit/polarity_amplitudes.py |
New joint polarity/P-amp/S/P misfit implementation and helpers |
tests/test_polarity_amplitudes.py |
Data-free unit tests for new misfit and helpers |
mtuq/grid_search.py |
Adds multi-component grid search + optional batching |
mtuq/grid/moment_tensor.py |
Adds orientation-grid helpers and changes random-grid defaults |
mtuq/util/math.py |
Adds MT decomposition/lune-coordinate utilities used by new plots |
mtuq/util/__init__.py |
Upgrades ProgressCallback; refactors DataArray min/max helpers |
mtuq/misfit/waveform/level2.py |
Optional chunking when using the new progress bar |
mtuq/graphics/beachball.py |
Accepts array-like MTs; adds best-DC overlay and used-data plotting |
mtuq/graphics/uq/_matplotlib.py |
Adds lune reference points/labels; tweaks lune axes/figure size |
mtuq/graphics/uq/_gmt/plot_lune |
Adds lune reference points/labels in GMT backend |
mtuq/misfit/__init__.py |
Re-exports new misfit class and dict helper |
mtuq/__init__.py |
Exposes new misfit from top-level API |
mtuq/graphics/__init__.py |
Exports new beachball plotting functions |
examples/Waveforms+Polarities+Pamp+SPratio.py |
New end-to-end joint inversion example |
pyproject.toml |
Fixes setuptools package discovery for mtuq.* |
.gitignore |
Adjusts ignored benchmark/example data + adds pytest cache |
Comments suppressed due to low confidence (2)
mtuq/grid/moment_tensor.py:24
- The default npts for FullMomentTensorGridRandom changed from 1,000,000 to 50,000. This is a significant behavioral change for users relying on default grid density and can affect inversion quality. Please confirm this is intentional and consider keeping the old default or documenting the change prominently (e.g., release notes).
def FullMomentTensorGridRandom(magnitudes=[1.], npts=50000):
""" Grid with randomly-drawn full moment tensors
Given input parameters ``magnitudes`` (`list`) and ``npts`` (`int`),
returns an ``UnstructuredGrid`` of size `npts*len(magnitudes)`.
mtuq/grid/moment_tensor.py:100
- The default npts for DeviatoricGridRandom changed from 100,000 to 50,000. Like the full MT grid, this reduces default search resolution and may change results for existing scripts. Please confirm this is intentional and consider documenting it or retaining prior defaults for backward compatibility.
def DeviatoricGridRandom(magnitudes=[1.], npts=50000):
""" Grid with randomly-drawn deviatoric moment tensors
Given input parameters ``magnitudes`` (`list`) and ``npts`` (`int`),
returns an ``UnstructuredGrid`` of size `npts*len(magnitudes)`.
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Comment on lines
+63
to
+67
| - ``'waveform'`` | ||
| Determine predicted polarity and amplitudes directly from full-waveform | ||
| synthetics | ||
|
|
||
|
|
Comment on lines
+649
to
+651
| if isinstance(data, tuple): | ||
| return True | ||
|
|
Comment on lines
+689
to
+696
| if norm == 'L2': | ||
| denominator = np.sum(weighted_obs**2) | ||
| if denominator <= 0.0: | ||
| return np.zeros(n_mt, dtype=np.float64) | ||
| numerator = np.sum( | ||
| (weights[None, :] * residual)**2, | ||
| axis=1) | ||
| return numerator / denominator |
| if r_i <= 0.0 or rho_i <= 0.0 or vp_i <= 0.0 or vs_i <= 0.0: | ||
| p_amplitude[:, _i] = Rp | ||
| sp_amplitude_ratio[:, _i] = np.log10( | ||
| np.sqrt(Rsv**2 + Rsh**2) / np.maximum(np.abs(Rp), eps) |
Comment on lines
+945
to
+956
| return | ||
|
|
||
| #model = _model_type(greens) | ||
| #method = _method_type(method) | ||
|
|
||
| #if model != method: | ||
| # print() | ||
| # print(' Possible inconsistency?') | ||
| # print(' Green''s functions are from: %s' % model) | ||
| # print(' Polarities are from: %s' % method) | ||
| # print() | ||
|
|
Comment on lines
+639
to
+646
| def _decompose_moment_tensor_ned(mt, vector=True): | ||
| mt_array = _full_mt_array_ned(mt) | ||
| eigenvalues, eigenvectors = np.linalg.eig(mt_array) | ||
| max_value = eigenvalues[np.argsort(np.abs(eigenvalues))[2]] | ||
| vol = np.sum(eigenvalues) / 3.0 | ||
| m_star = eigenvalues - vol | ||
| sorted_index = np.argsort(np.abs(m_star)) | ||
| max_value_star = m_star[sorted_index[2]] |
Comment on lines
+673
to
+678
| def _mom2other_ned(mt): | ||
| mt_array = _full_mt_array_ned(mt) | ||
| eigenvalues, eigenvectors = np.linalg.eig(mt_array) | ||
| sorted_index = np.argsort(eigenvalues) | ||
| p_axis = _v2trpl(eigenvectors[:, sorted_index[0]].copy()) | ||
| t_axis = _v2trpl(eigenvectors[:, sorted_index[2]].copy()) |
Comment on lines
+601
to
+603
| scale = np.percentile(values, percentile) | ||
| if scale <= 0.0: | ||
| return values |
Comment on lines
+288
to
+299
| results_polarity_amplitude = grid_search_multicomponent( | ||
| polarity_amp_data, greens, polarity_amplitude_misfit, origin, grid | ||
| ) | ||
|
|
||
| if comm.rank == 0: | ||
| results = results_bw + results_sw | ||
|
|
||
| # `grid` index corresponding to minimum misfit | ||
| idx = results.source_idxmin() | ||
|
|
||
| best_mt = grid.get(idx) | ||
| lune_dict = grid.get_dict(idx) |
| """Returns coordinates for the first finite global min/max of a DataArray.""" | ||
| values = np.asarray(da.values) | ||
| if values.size == 0: | ||
| raise ValueError("Cannot index an empty DataArray") |
Author
|
Addressed the Copilot review feedback in commit 54490c9. Main fixes:
|
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
PR: Add P-amplitude and S/P-amplitude-ratio inversions to MTUQ
Branch:
add-ps-amplitude-inversionFork: https://github.com/MHkhosravi/mtuq
Target:
mtuqorg/mtuqmasterSummary
This PR adds optional P-amplitude and S/P-amplitude-ratio misfits to MTUQ's
existing waveform and polarity inversion workflows.
The main addition is
PolarityPampSPampRatioMisfit, a new misfit class that cancombine first-motion polarity, P-amplitude, and S/P-amplitude-ratio information
with MTUQ's grid-search workflow. All three terms are independently weighted and
can be used in any combination, including polarity-only (backward-compatible with
the existing
PolarityMisfituse pattern) and all three jointly.Main changes
New files
mtuq/misfit/polarity_amplitudes.pyPolarityPampSPampRatioMisfit, radiation-pattern calculator_pol_pamp_spamp_mt, and helper utilitiesexamples/Waveforms+Polarities+Pamp+SPratio.pytests/test_polarity_amplitudes.pyModified files
mtuq/__init__.pyPolarityPampSPampRatioMisfitandpolarities_pamp_spamp_ratio_from_dictfrom the top-level APImtuq/misfit/__init__.pyPolarityMisfitandWaveformMisfitmtuq/grid_search.pygrid_search_multicomponentfor multi-term weighted misfit grid searchesmtuq/grid/moment_tensor.pymtuq/util/math.pyto_delta,to_gamma,lune_*coordinate utilities used by new plotsmtuq/util/__init__.pymtuq/util/beachball.pymtuq/graphics/beachball.pyplot_beachball_withbestDC,plot_beachballs_useddata; relaxplot_beachballtype constraint to accept array-like moment tensorsmtuq/graphics/__init__.pymtuq/graphics/uq/_matplotlib.py_add_marker,_generate_lunehelpers used by new scatter plotsmtuq/graphics/uq/_gmt/plot_lunemtuq/misfit/waveform/level2.pypyproject.toml[tool.setuptools]→[tool.setuptools.packages.find]so that all sub-packages (mtuq.*) are found by setuptoolsMotivation
MTUQ already supports waveform-based and polarity-based moment tensor workflows.
This PR extends those capabilities by allowing users to include P-amplitude and
S/P-amplitude-ratio observations as additional misfits.
P-wave first-motion amplitude and S/P amplitude ratio can provide complementary
source information, particularly for smaller events where surface-wave signal is
limited or waveform fits are non-unique. These additional observations can be
read from station picks (e.g., CAP format files) and incorporated alongside the
waveform grid search with user-controlled weights, so the contribution of each
data type can be tuned.
Implementation notes
Radiation pattern
Predicted P polarity, P amplitude, and S/P amplitude ratio are computed from
far-field radiation patterns in NED coordinates using the Aki & Richards
formulation. USE-convention MTUQ moment tensors are converted to NED internally.
Takeoff angles and azimuths come from TauP (default) or FK/CPS metadata.
Optional geometric spreading (
rho,vp,vs,distance) converts radiationpatterns to physical amplitude units. When
method='taup', these parameters areextracted automatically from the TauP velocity model at the source depth.
Zero S-radiation robustness (bugfix in this PR)
Isotropic or S-nodal sources produce exactly zero S radiation at some or all
stations. Without a guard,
log10(0)propagates as-infthrough the S/P ratioand poisons the misfit surface. A
np.maximum(s_mag, eps)floor on the Snumerator (both geometric-spreading and plain-radiation branches) keeps the ratio
finite. A dedicated unit test (
test_isotropic_sp_ratio_is_finite) covers bothbranches.
Misfit normalisation
Each term (polarity, P-amplitude, S/P-ratio) is independently L1- or L2-normalised
by the sum of weighted absolute observed values. This makes each term dimensionless
and comparable when combined via
lambda_pol,lambda_pamp,lambda_sprat.Stations with zero weight, zero polarity (unpicked), missing amplitude, or
non-positive S/P ratio are automatically excluded from the relevant term.
Backward compatibility
PolarityMisfitis preserved unchanged.PolarityPampSPampRatioMisfitaccepts polarity-only data input fordrop-in compatibility with existing polarity-based workflows.
plot_beachballnow accepts bothMomentTensorobjects and six-elementUSE-convention arrays; previous code only accepted
MomentTensor.are unaffected.
Verification
Compilation and import
Unit tests (data-free, fast)
Waveform misfit regression
Example (requires MPI, data download, and ~4 minutes runtime)
Full end-to-end validation (polarity+amplitude grid search over 210k sources)
A lightweight validation script was run and produced:
Open questions / flags for reviewer
The following items are scientifically or technically ambiguous and should be
reviewed by domain experts before merging:
P-amplitude units — When
lambda_pamp > 0, predicted P amplitudes (fromradiation pattern × geometric spreading) must have the same units as the
observed values supplied in
pamp_dict. The code does not enforce or documentthis requirement beyond a docstring note. Users who supply raw waveform
peak amplitudes without normalising by the same 1/(4πρv³r) factor will get
a misfit that is numerically large but physically meaningless. Consider adding
a warning or example guidance.
sp_obs_is_log10flag in example — The example setssp_obs_is_log10=Trueand supplies
spamp_dictvalues in the range 0.46–2.07. A comment says"This CAP field already stores log10(S/P)", which would mean the actual S/P
ratios are 10^0.46 ≈ 2.9 to 10^2.07 ≈ 117. Whether the original CAP data
stores log10 or linear ratios should be confirmed against the data source.
FK/CPS methods declared but not implemented —
method='FK_metadata'andmethod='CPS_metadata'are accepted by__init__but raiseNotImplementedErrorimmediately when any calculation is attempted. Considerraising in
__init__instead to give an earlier, clearer error.Default grid size reduction —
FullMomentTensorGridRandomdefaultnptswas reduced from 1,000,000 to 50,000 and
DeviatoricGridRandomfrom 100,000to 50,000. Confirm this is intentional and not a debugging leftover.
Files NOT included in this PR
data/tests/benchmark_cap/,data/tests/benchmark_cps/(downloaded test data,now gitignored)
data/examples/(gitignored)build/,*.egg-info/,__pycache__/,.pytest_cache/(gitignored)*.png,*.nc,*.jsonoutput filesmtuq_test/validation outputsRebase note
At time of preparation this branch (
add-ps-amplitude-inversion) is 7 commitsbehind
upstream/master. The upstream commits touchmtuq/__init__.py,mtuq/graphics/beachball.py, andpyproject.toml— all in different lines fromour changes.
git merge-treeconfirms a clean three-way merge.Rebase before opening the PR: