Skip to content

Make source-to-site distance a configurable registry, and specify the magnitude calculation - #25

Merged
sgjholt merged 9 commits into
mainfrom
claude/distance-registry-and-magnitude-roadmap
Aug 8, 2026
Merged

Make source-to-site distance a configurable registry, and specify the magnitude calculation#25
sgjholt merged 9 commits into
mainfrom
claude/distance-registry-and-magnitude-roadmap

Conversation

@sgjholt

@sgjholt sgjholt commented Aug 8, 2026

Copy link
Copy Markdown
Owner

Follows PR #24. One piece of working code — the distance registry — and the magnitude roadmap it serves, sourced from the thesis rather than guessed.

The bug this started from

staged.InverseDistance, added in #24, hardcoded rhyp. Meanwhile [windows] distance_metric said repi and was read by nothing. So the project declared epicentral and computed hypocentral — the same defect class this refactor has now found seven times.

It is not bookkeeping at microseismic range. The nearest PNR station is 0.89 km epicentral against 2.30 km hypocentral, a factor of 2.57, while the farthest agree to 1.00. Since the weighting is by inverse distance, the two measures disagree hardest exactly where the weight is largest. Every number PR #24 quoted moves:

was (rhyp) now (repi)
event fc 12.751 Hz 11.585 Hz
weight carried by nearest 2 channels 22.8% 41.4%
dropping the nearest station +16%, 1.56× stress drop +34.7%, 2.44×

Epicentral is the right default beyond the config having said so. rhyp is built from the source depth and the station elevation, so it assumes every sensor sits at the surface; for a borehole deployment that is wrong by the burial depth and nothing announces it. The PNR inventory records channel depth as 123456.0 — a placeholder — so on that dataset rhyp is an assumption wearing a measurement's name.

specmod.distance

A registry, in the same shape as ESTIMATORS, NOISE_MODELS and WEIGHT_MODELS. It lives outside staged because geometric spreading needs distance too, and a second consumer restating the choice is how two things drift apart.

repi and rhyp are implemented. rrup and rjb are registered and raise, deliberately rather than being omitted:

rrup is not implemented. It needs a rupture surface — strike, dip, length
and width, and SpecMod carries a point source... For a point source rrup
degenerates to hypocentral; name that instead if it is what you mean,
rather than getting it by accident.

That degeneracy is the reason a silent fallback would be dangerous — plausible numbers, wrong for precisely the events large enough to warrant asking. A name that resolves to a clear failure is a better extension point than a name that does not resolve.

InverseDistance(measure=None) follows the configuration; explicit inverse_hypocentral_distance and inverse_epicentral_distance spellings remain for a study that wants to record which it used regardless.

[geometry] distance_measure

The setting moves out of [windows], which was the wrong home even when nothing read it — cutting a window does not depend on how distance is measured. Renamed as well as moved, since these are measures rather than metrics and the name now matches the registry.

studies/magna_2020_paper.toml carries the value across unchanged. The published run used hypocentral, and that is preserved as the record of what was done rather than edited toward the new default. Magna is a transcript, not a recommendation.

Caught while moving it: deleting the old field took p_velocity and s_velocity with it, since the three shared a comment block. Five test_config.py failures named the missing keys immediately — the value of a test that asserts the whole resolved configuration rather than the fields someone remembered to check.

Roadmap, sourced rather than guessed

Read from Holt (2019) — the Edwards et al. (2010) method this package implements. Chapters 1 and 2 have now been read in full. Two earlier claims in the plan were wrong and are corrected:

  • ΘλΦ = 0.55 is not a different average of the same quantity, it is a different quantity — the average radiation pattern of SH propagation, pairing with F = 2 for vertically incident SH. The textbook 0.63 is the RMS over total S. So it was never a departure from the standard value, and it couples directly to the rotation question: a formulation using the SH radiation pattern wants the transverse component.
  • R₀ = 1000 m is a reference source distance, not a kilometre-to-metre bridge as an earlier draft guessed.

And that resolves the disagreement about the unit of R: there are two. R₀ in the moment expression is in metres alongside ρ in kg/m³ and β in m/s, giving M₀ in N·m. The geometric spreading model S(R) is a separate term and its R is in kilometres. Both answers were right about different terms.

Also recorded, each with the measurement behind it:

  • The spreading exponent is ~1, not 2. 1/r² gives median Mw +5.39 against a catalogue 1.6 — the distance term applied twice. The thesis's own Utah inversion gives 0.88 ± 0.02 over 1–40 km, independently.
  • Bilinear belongs in the design but one event cannot fit it. Regressing log10(Ω) on log10(R) over the 28 channels buys 5% of rms for two extra parameters across one decade of distance — a hinge finding scatter. The 0.265 log10 residual is site response and radiation pattern, inseparable from spreading when every station contributes one distance. Which is the argument for the non-parametric inversion.
  • Magna stays the calibration anchor, because its S(R) is non-parametric where Utah's is piecewise parametric — and S(R) is the term least checkable by other means. Consequence recorded: it must arrive as tabulated data, so the spreading registry has to accept a table and not only a functional form.
  • Components, phases and station/channel identity types (§4.7), including that back_azimuth is computed and read by nothing, that horizontals are fitted as independent measurements so each station is counted twice, and the N/E vs 1/2 trap.
  • Documentation equations do not render — MyST dollarmath syntax with no MyST build.

The equations cannot be verified from the PDFs, and that route is now closed

An earlier draft of the plan recorded the constants as unverified and suggested a PDF export of Chapters 1–2 under 10 MB would settle them. Those exports were obtained. They do not, and the reason rules out trying again.

It is not an extraction failure. The pages carrying Eq. 1.12, 1.13, 1.14, 2.6, 2.7 and 2.8 declare only Times New Roman fonts and no image, and rendering them produces the equation numbers against blank space. The Word-to-PDF export dropped the embedded equation objects, taking inline maths with it — which is why the prose reads "where, is the long period spectral displacement plateau at the source". Settling the algebra needs the original document, or Edwards et al. (2010).

The prose survived, and settles every constant in the table, each now read from the section that defines it. Two things follow:

  • The partition factor should not be folded in by default. Ch. 1 lists one in the generic C; Ch. 2, defining the Utah calculation this package follows, enumerates Eq. 2.7's constants without one. A study that wants it can set it.
  • Table 2.1 carries a preferred model the plan was missing — Holt et al. [R]: 0.90 / 2.57 / 0.44 / 1.54, breaking at 43/76/136 km. §2.7 recommends it over [O] on lower uncertainty across every slope and more events resolved. If a Utah spreading table ever ships as study data, [R] is the one.

Mw is a near miss: Ch. 1 confirms the SI (N·m) convention around it, but the numeric constant sits inside a dropped equation, so the 9.1 is still unverified from this source.

Four transform defaults now have a published source

The equations did not survive the export, but Ch. 2 §2.2's description of the processing did, and it is the more useful half here — four steps land on defaults §4.4 and §4.5 previously justified on internal measurement alone.

  • MultitaperEstimator defaults time_bandwidth=3.0, n_tapers=5 — exactly "five 3-π prolate tapers", NW=3 giving 2NW-1=5. §4.4 could only describe the 3 as the literal argument passed positionally to mtspec, with no way to configure it. It is now a considered choice with a citation.
  • "Multiply by the sampling period" is one_sided_fas's |spec| * dt, reached there from the requirement that padding change only the frequency sampling. Same normalisation, opposite direction.
  • Excluding zero-padding from the signal/noise length ratio is the same rule as that function's refusal to rescale for padding.
  • SNR ≥ 3 matches SnrConfig.tolerance, but applied to a different job: the thesis gates the fit band per spectrum where the Magna paper gates spectrum selection in fixed bands. That is already the distinction between bandwidth selection and assert_bandwidths, so both are expressible — said explicitly so the shared number is not misread.

One knowing divergence recorded: the recipe pads to 2^N where resolve_n_fft defaults to a 5-smooth length, safe for the same reason the recipe itself relies on, with pow2 available for exact reproduction. The remaining steps — the 20 s window at 80% of the Pg–Sg interval, the noise window at 75%, 2.5σ station rejection, 3-station minimum — are event-level aggregation this package does not do yet, written down so the magnitude module has a specification to build to instead of inventing thresholds.

Checks

  • 592 passed, 1 skipped
  • ruff check, ruff format --check, mypy clean
  • Tutorial notebook re-executed end to end; its numbers re-measured on epicentral distance

The last two commits are documentation only — no code changed since the checks above were run.

claude added 9 commits August 8, 2026 16:56
…y measurement

Records the design for a basic Mw and notes it is nearly free: `llpsp` from
the fit is already `log10(Omega)` on displacement — `sources.motion_scaling`
converts the displacement model to whatever motion was recorded, so the fitted
plateau needs no correction — and `rhyp` is on every trace.

One correction to the proposal, backed by measurement rather than assertion.
The spreading exponent should be **1, not 2**. Body-wave amplitude from a
point source decays as `1/r` in a whole space; `1/r**2` is how energy decays,
and the moment formula corrects an amplitude.

Measured on the 28 PNR windows through the two-stage fit, with round-number
constants (rho 2500, beta 2500, R 0.63, F 2):

    1/r      median M0 1.67e13 N m   median Mw +2.75
    1/r**2   median M0 1.55e17 N m   median Mw +5.39
    catalogue Mw 1.6                          +1.60

`1/r**2` is nearly four magnitude units out — the distance term applied twice.
`1/r` lands close enough that the remainder sits inside the constants: M0 goes
as rho * c**3, so beta from 2500 to 1600 m/s is a factor of 3.8 by itself,
about 0.4 magnitude units. The exponent is settled by this; the absolute
calibration is not, and wants the study's own velocity model rather than the
round numbers used here.

The exponent stays configurable with a default of 1.0 all the same, since a
study fitting an empirical spreading term is a legitimate thing to want and
should have to say so rather than being unable to.

Also notes the two preconditions from §4.7 that must hold before the number
means anything: Omega should be the combined horizontal rather than one
component, and the phase constants must match the phase measured. And the test
that makes the units discussion enforceable rather than advisory — asserting
the Mw of a known event, which is what stops a factor of 10^9 from a km/s
velocity going unnoticed.
…ation settles it

An earlier version of this section asserted that R must be in metres. That is
true of the all-SI formulation and not true in general, and the correction is
worth more than the original claim.

With rho in kg/m^3, c in m/s and Omega in m s, dimensional consistency does
force R into metres — kg/m^3 * m^3/s^3 * m * m s is kg m^2 s^-2, newton-metres.
But most of the classic literature writes the same formula in CGS-lineage
units, with rho in g/cm^3, c in km/s, R in km and Omega in cm s, giving
dyne-centimetres. Both are internally consistent. "R is in km" is a statement
about which formulation is in use, not about the physics.

Which is exactly why it matters: the two differ by 10^3 in R alone, two
magnitude units, and mixing them is the easiest way to be quietly wrong.
Measured on the same fits, R in metres gives median Mw +2.75 and R in
kilometres +0.75, with the catalogue Mw 1.6 sitting between them. So the unit
choice alone does not reconcile it, and neither does any single constant —
landing on 1.6 with R in metres and rho 2500 would need beta 667 m/s, which is
not a shale S velocity. Something else is also off, most likely that Omega
here is one horizontal component rather than the combined horizontal, and that
round numbers for the medium are not this event's.

The conclusion is the useful part. Two people can both be certain about the
unit of R and both be right about different formulations, and no amount of
arguing resolves it. Pin the study's own rho and c, state the unit of every
input, and assert the Mw of a known event in a test — the calibration is the
deliverable, not the derivation.
…oth R units

Replaces guesswork in §4.7 with what Holt (2019) actually says — Ch. 1 §1.4
and Ch. 2 Eq. 2.7, the Edwards et al. (2010) spectral method this package
implements. Two earlier claims in this document were wrong.

`Theta-lambda-Phi = 0.55` is not a different average of the same quantity, it
is a **different quantity**: the average radiation pattern of **SH**
propagation over the focal sphere, which pairs with `F = 2` for vertically
incident SH. The textbook 0.63 is the RMS over total S. So the coefficient is
not "not the textbook value" — it is the right value for the phase and
component being measured, and that couples directly to the rotation question:
a formulation using the SH radiation pattern wants the transverse component,
not an arbitrary horizontal.

`R_0 = 1000 m` is a **reference source distance** — where the source spectrum
is defined — not a kilometre-to-metre bridge, which an earlier draft guessed
it might be.

And that resolves the disagreement about the unit of R, because there are two
different distances in one workflow. `R_0` in the moment expression is in
metres, alongside rho in kg/m^3 and beta in m/s, giving M0 in newton-metres.
The geometric spreading model `S(R)` is a separate term in Eq. 1.14 and its R
is in **kilometres** — the thesis tabulates it piecewise that way. Both
answers were right about different terms.

The spreading table also supports the exponent argument independently: the
inverted near-field decay for Utah is 0.88 +/- 0.02 over 1 < R <= 40 km, close
to the theoretical body-wave 1 and nowhere near 2. It also shows spreading is
not a single power law but a piecewise empirical function inverted per region,
so `S(R)` should be a registry in its own right — theoretical 1/R as the
default, region-specific tables as what an operator supplies — in the same
shape as WEIGHT_MODELS and NOISE_MODELS.

Records the constants as a table with rho = 2600 kg/m^3 for Utah, and notes
the energy-partition factor (Boore, 2003) that none of this document mentioned
before.

Caveat recorded with the source: the equations are embedded objects that do
not survive text extraction, so the symbols are read from surrounding prose
and the equation images should be checked against the table before anything is
implemented.
I suggested pointing the Mw calibration at the thesis's Utah catalogue
instead, on the grounds that >200 published Mw with regional constants and
bilinear Mw-ML relationships is more validation material than Magna offers.
Declined, and the reason is better than the suggestion: Magna's geometric
spreading is the better-constrained half.

Utah's S(R) is the piecewise parametric table in Ch. 2 — segment boundaries
chosen, an exponent fitted per segment. Magna used the non-parametric G(R)
inversion already noted in §5.2.5, with no imposed functional form and no
chosen breakpoints.

That ordering matters more than catalogue size, because S(R) is the term this
calculation is least able to check by other means. Density, velocity and the
radiation pattern are each bounded by physics and by the literature; a
spreading function is bounded only by the inversion that produced it. So
validating against the dataset with the better-constrained S(R) tests the part
of the formula that most needs testing.

One consequence for scope, recorded: the non-parametric inversion is not in
this repository, so S(R) for the Magna comparison has to arrive as **data** —
a tabulated G(R) in studies/magna_2020.toml — rather than being recomputed.
The spreading registry therefore has to accept a table and not only a
functional form. That is a requirement, not a nicety.
…annot fit it

The point stands that microseismic distances are not the regional regime.
Every spreading model in the thesis breaks at 40-50 km, and the PNR data used
throughout this repository spans 2.3 to 22.9 km — entirely inside the first
segment of all of them. Regional tables therefore offer no guidance at these
distances, and a shape fitted over 1-400 km should not be assumed to extend
downward. A bilinear form with a break inside the microseismic range is the
right thing for the design to allow.

Measured before writing it in, because a single event is exactly the case
where this looks easy and is not. For one event `Omega_source` is a constant,
so regressing log10(Omega) on log10(R) over the 28 channels *is* the spreading:

    single segment    exponent 0.73                      rms 0.265 log10
    bilinear          0.41 then 1.89, break 13.5 km      rms 0.253 log10

Five percent of rms for two extra parameters over one decade of distance. That
is a hinge finding scatter, not a break. And the scatter is the substance:
0.265 log10 units is a factor of 1.8, which for a single event is site
response and radiation pattern — neither separable from spreading when every
station contributes exactly one distance.

Which is the argument for the non-parametric inversion rather than against
bilinear, and reinforces the Magna choice recorded below. Spreading separates
from site only across a dataset where each station sees many distances and
each distance is sampled by many stations. The apparent 0.73 here is not a
measurement of anything and is recorded as such.

Design conclusion unchanged in shape, sharper in default: allow bilinear and
tabulated forms, default to theoretical 1/R, and resist fitting a spreading
exponent to a single event. The default is defensible precisely because it is
not fitted.
`[windows] distance_metric` has existed since the config package landed and
was read by nothing — the same defect this session has now found seven times.
Meanwhile `staged.InverseDistance`, written two commits ago, hardcoded `rhyp`.
So the project declared epicentral and computed hypocentral.

That is not bookkeeping at microseismic range. The nearest PNR station is
**0.89 km epicentral against 2.30 km hypocentral**, a factor of 2.57, while
the farthest agree to 1.00 — and since the weighting is by *inverse* distance,
the two measures disagree hardest exactly where the weight is largest. Every
number this branch has quoted moves:

    event fc          12.751 -> 11.585 Hz
    nearest 2 weight    22.8 ->   41.4 %
    drop AQ04           +16% ->  +34.7%  (1.56x -> 2.44x in stress drop)

Epicentral is the right default and not only because the config said so.
`rhyp` is built from the source depth and the station *elevation*, so it
assumes every sensor sits at the surface; for a borehole deployment that is
wrong by the burial depth and nothing announces it. The PNR inventory records
channel `depth` as `123456.0` — a placeholder — so on this dataset `rhyp` is
an assumption wearing a measurement's name.

New `specmod.distance` holds the registry, in the same shape as ESTIMATORS,
NOISE_MODELS and WEIGHT_MODELS. It lives outside `staged` because geometric
spreading needs distance too, and a second consumer restating the choice is
how the two drift apart.

`rrup` and `rjb` are registered and **raise**, deliberately rather than being
omitted. Both need a rupture surface, and for a point source they degenerate
exactly to hypocentral and epicentral — so a silent fallback would give
plausible numbers that are wrong for precisely the events big enough to
warrant asking. The error says what they would need, which puts the
requirement where whoever adds finite-fault support will read it. A name that
resolves to a clear failure is a better extension point than a name that does
not resolve.

`InverseDistance(measure=None)` follows the configuration; the explicit
`inverse_hypocentral_distance` and `inverse_epicentral_distance` spellings stay
for a study that wants to record which it used regardless. Notebook, tests and
plan re-measured throughout.
…tion

`distance_metric` sat in `[windows]`, which was the wrong home even when
nothing read it: cutting a window does not depend on how distance is measured.
It now has two readers — the ensemble weighting of the two-stage fit, and the
geometric spreading a moment calculation will correct for — so leaving it there
would have meant two stages reaching into a third's settings, which is how a
setting ends up copied and the copies disagree.

`[geometry] distance_measure`, resolved through
`specmod.distance.DISTANCE_MEASURES`. Renamed as well as moved: these are
measures rather than metrics, and the name now matches the registry.

`studies/magna_2020_paper.toml` carries the setting across unchanged. The
published run used **hypocentral** distance, and that is preserved as the
record of what was done rather than edited toward the new default — `repi`
ships as the default because it does not assume sensor depths that are often
unknown, but Magna is not a recommendation, it is a transcript.

Caught while doing it: removing the old field took `p_velocity` and
`s_velocity` with it, since the three shared a comment block. The five
`test_config.py` failures named the missing keys immediately, which is the
value of a test that asserts the whole resolved configuration rather than the
fields someone remembered to check.
… equations

Attempted to check the constants table against the equations, first through
text extraction of the 44 MB original and then against separate exports of
Chapters 1 and 2. **The equations cannot be verified from the PDFs at all,
and the reason is worth stating precisely, because it rules out the obvious
next attempt.**

It is not an extraction failure. The pages carrying Eq. 1.12, 1.13, 1.14,
2.6, 2.7 and 2.8 declare only Times New Roman fonts and no image XObject, and
rendering them produces the equation numbers against blank space. The
Word-to-PDF export dropped the embedded equation objects, taking inline maths
with it, which is why the prose reads "where, is the long period spectral
displacement plateau at the source". So the algebraic form of Eq. 1.12, 1.13,
2.7 and 2.8 remains unverified and no re-export will recover it. Settling it
needs the original document, or Edwards et al. (2010), which is the published
method this package implements anyway.

The prose, however, survived, and it settles every constant in the table --
each now read from the section that defines it rather than inferred. One
correction falls out: Ch. 2 enumerates the constants of Eq. 2.7 without a
partition factor, though Ch. 1 lists one in the generic C. The implementation
should not fold one in by default; a study that wants one can set it.

Table 2.1 survives as real text and confirms the exponents already quoted. It
also carries the refined model the thesis actually recommends, which was
missing here -- Holt et al. [R], lower uncertainty on every slope, more events
resolved. If a Utah spreading table is ever shipped as study data, [R] is the
one to ship.

Two things this section had inferred are now cited rather than deduced. The
spreading model's distance is in kilometres, stated outright by both Eq. 2.6's
description and Table 2.1's caption -- that is the term whose unit was in
question two commits ago. And the measurement is on the horizontal component
of the Sg phase, "referred to as SgH", which is exactly the pairing this
section infers from Theta-lambda-Phi being the SH average: phase, component
and radiation-pattern constant are one choice made once, which is the argument
for making them one setting rather than three.

Mw is a near miss. Ch. 1 confirms the SI convention around it, but the
numeric constant sits inside a dropped equation, so the 9.1 is still
unverified from this source.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The thesis equations did not survive its PDF export, but the prose
describing the *processing* did, and it is the more useful half for this
repository. Ch. 2 §2.2 specifies the workflow step by step, and four of
those steps land on defaults §4.4 and §4.5 currently justify on internal
measurement alone.

The multitaper pair is the one worth having. MultitaperEstimator defaults
to time_bandwidth=3.0 with n_tapers=5, which is exactly "five 3-pi prolate
tapers" -- NW=3 giving 2NW-1=5. §4.4 could previously only describe the 3
as the literal argument passed positionally to mtspec, with no way to
configure it. It is now a considered choice with a source, not an
inherited constant.

Two more agree from the other direction. "Multiply by the sampling period"
is one_sided_fas's |spec| * dt, arrived at there from the requirement that
padding change only the frequency sampling. And the recipe's care to
exclude zero-padding from the signal/noise length ratio is the same rule
as that function's refusal to rescale for padding.

The fourth, SNR >= 3, matches SnrConfig.tolerance but is applied to a
different job -- the thesis gates the fit band per spectrum where the
Magna paper gates spectrum selection in fixed bands. That distinction is
already the one between bandwidth selection and assert_bandwidths, so both
are expressible; the section says so rather than leaving the coincidence
of the number to be misread.

One knowing divergence is recorded: the recipe pads to 2^N where
resolve_n_fft defaults to a 5-smooth length, which is safe for the same
reason the recipe itself relies on, and pow2 stays available for exact
reproduction.

The remaining steps are event-level aggregation this package does not do
yet. They are written down so the magnitude module has a specification to
build to instead of inventing thresholds.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@sgjholt
sgjholt force-pushed the claude/distance-registry-and-magnitude-roadmap branch from cf8e1dc to 4236e42 Compare August 8, 2026 17:55
@sgjholt
sgjholt merged commit 1b4639b into main Aug 8, 2026
10 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants