From 32c39a8307dffc7858712195d28888c60f430f18 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 8 Aug 2026 16:26:07 +0000 Subject: [PATCH 1/9] docs(plan): specify the basic Mw, and settle the spreading exponent by measurement MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- docs/REFACTOR_PLAN.md | 46 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/docs/REFACTOR_PLAN.md b/docs/REFACTOR_PLAN.md index 629cefc..97220be 100644 --- a/docs/REFACTOR_PLAN.md +++ b/docs/REFACTOR_PLAN.md @@ -1341,6 +1341,52 @@ of each station" is something the code can express, and "rotate to transverse" cannot be written correctly until the code can ask a channel which way it points. +#### A basic Mw, and the spreading exponent + +Worth building, and nearly free once the constants have units: everything it +needs is already produced. `llpsp` from the fit is `log10(Omega)` on +**displacement** — `sources.motion_scaling` converts the displacement model to +whatever motion was recorded, so the fitted plateau is `Omega` directly rather +than a velocity plateau needing correction — and `rhyp` is on every trace. + + M0 = 4 * pi * rho * c**3 * r**n * Omega / (R_c * F) + Mw = (2/3) * (log10(M0) - 9.1) + +**The exponent `n` is 1, not 2, and this is worth stating because it is easy +to say the other one.** Body-wave *amplitude* from a point source decays as +`1/r` in a homogeneous whole space; `1/r**2` is how *energy* decays. The +moment formula corrects an amplitude, so it multiplies by `r`. + +Measured on the 28 PNR windows through the two-stage fit, with `rho = 2500` +kg/m^3, `beta = 2500` m/s, `R_c = 0.63`, `F = 2`: + +| spreading | median M0 | median Mw | +|---|---|---| +| `1/r` | 1.67e13 N m | **+2.75** | +| `1/r**2` | 1.55e17 N m | **+5.39** | +| catalogue Mw 1.6 | 3.16e11 N m | +1.60 | + +`1/r**2` is nearly four magnitude units out, which is the distance term +applied twice; `1/r` lands close enough that the remainder sits inside the +constants. Note what the remainder is worth: `M0` goes as `rho * c**3`, so +taking `beta` from 2500 to 1600 m/s — entirely plausible for shallow Bowland +shale — is a factor of 3.8 on its own, about 0.4 magnitude units. **The +exponent is settled by this; the absolute calibration is not, and needs the +study's own velocity model** rather than the round numbers used here. + +Two other things that must be true before the number means anything, both +from §4.7 above: `Omega` should be the combined horizontal rather than one +component, and the phase constants must match the phase actually measured. + +So the shape is a `magnitude` module taking a `StagedFit` and a small typed +set of medium constants — density, velocity, radiation pattern, free-surface +factor, spreading exponent — with the exponent **configurable and defaulting +to 1.0**, since a study fitting an empirical spreading term is a legitimate +thing to want and should have to say so. `[model]` holds them, `studies/*.toml` +pins them, and the output carries its unit. A test asserting the Mw of the PNR +event to within a stated tolerance is what stops a factor of 10^9 from a km/s +velocity going unnoticed — the whole point of §4.7's units discussion. + #### Suggested shape - `preprocess.rotate_to_rt(st)` reading the stored back-azimuth, and a From fa3a4badde12920cee60263a1722c2252c7866a9 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 8 Aug 2026 16:29:06 +0000 Subject: [PATCH 2/9] docs(plan): the unit of R depends on the formulation, and only calibration settles it MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- docs/REFACTOR_PLAN.md | 42 ++++++++++++++++++++++++++++++++++++------ 1 file changed, 36 insertions(+), 6 deletions(-) diff --git a/docs/REFACTOR_PLAN.md b/docs/REFACTOR_PLAN.md index 97220be..f4f9c74 100644 --- a/docs/REFACTOR_PLAN.md +++ b/docs/REFACTOR_PLAN.md @@ -1367,12 +1367,42 @@ kg/m^3, `beta = 2500` m/s, `R_c = 0.63`, `F = 2`: | catalogue Mw 1.6 | 3.16e11 N m | +1.60 | `1/r**2` is nearly four magnitude units out, which is the distance term -applied twice; `1/r` lands close enough that the remainder sits inside the -constants. Note what the remainder is worth: `M0` goes as `rho * c**3`, so -taking `beta` from 2500 to 1600 m/s — entirely plausible for shallow Bowland -shale — is a factor of 3.8 on its own, about 0.4 magnitude units. **The -exponent is settled by this; the absolute calibration is not, and needs the -study's own velocity model** rather than the round numbers used here. +applied twice; `1/r` lands within reach of the constants. **The exponent is +settled by this; the absolute calibration is not.** + +**And the unit of `R` cannot be settled by argument, which is the useful +lesson.** With `rho` in kg/m^3, `c` in m/s and `Omega` in m s, dimensional +consistency forces `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, `rho` in g/cm^3, `c` in km/s, `R` in **km** +and `Omega` in cm s, giving dyne-centimetres. Both are right; they are +different formulations, and "R is in km" is a statement about which one is in +use rather than about the physics. + +That matters because the two differ by 10^3 in `R` alone — two magnitude units +— and mixing them is the single easiest way to be quietly wrong. Measured, with +the same fits: + +| rho kg/m^3 | beta m/s | R in | median M0 (N m) | median Mw | +|---|---|---|---|---| +| 2500 | 2500 | m | 1.67e13 | +2.75 | +| 2500 | 2500 | km | 1.67e10 | +0.75 | +| 2500 | 1600 | m | 4.37e12 | +2.36 | +| 2700 | 3500 | km | 4.94e10 | +1.06 | + +The catalogue Mw 1.6 sits *between* the metre and kilometre answers, 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 — the likeliest candidates being that `Omega` here is one horizontal +component rather than the combined horizontal (§4.7), and that the round +numbers for the medium are not this event's. + +So the calibration is the deliverable, not the derivation. 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. Two people can both be certain about the unit of `R` and both +be right about different formulations; only the calibration tells you which +one the code is in. Two other things that must be true before the number means anything, both from §4.7 above: `Omega` should be the combined horizontal rather than one From 391ccc0ad0602c5065a0573a4c60f9dbc17b42d5 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 8 Aug 2026 16:34:35 +0000 Subject: [PATCH 3/9] docs(plan): source the moment constants from the thesis, and settle both R units MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- docs/REFACTOR_PLAN.md | 97 ++++++++++++++++++++++++++++--------------- 1 file changed, 64 insertions(+), 33 deletions(-) diff --git a/docs/REFACTOR_PLAN.md b/docs/REFACTOR_PLAN.md index f4f9c74..d499aff 100644 --- a/docs/REFACTOR_PLAN.md +++ b/docs/REFACTOR_PLAN.md @@ -1370,39 +1370,70 @@ kg/m^3, `beta = 2500` m/s, `R_c = 0.63`, `F = 2`: applied twice; `1/r` lands within reach of the constants. **The exponent is settled by this; the absolute calibration is not.** -**And the unit of `R` cannot be settled by argument, which is the useful -lesson.** With `rho` in kg/m^3, `c` in m/s and `Omega` in m s, dimensional -consistency forces `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, `rho` in g/cm^3, `c` in km/s, `R` in **km** -and `Omega` in cm s, giving dyne-centimetres. Both are right; they are -different formulations, and "R is in km" is a statement about which one is in -use rather than about the physics. - -That matters because the two differ by 10^3 in `R` alone — two magnitude units -— and mixing them is the single easiest way to be quietly wrong. Measured, with -the same fits: - -| rho kg/m^3 | beta m/s | R in | median M0 (N m) | median Mw | -|---|---|---|---|---| -| 2500 | 2500 | m | 1.67e13 | +2.75 | -| 2500 | 2500 | km | 1.67e10 | +0.75 | -| 2500 | 1600 | m | 4.37e12 | +2.36 | -| 2700 | 3500 | km | 4.94e10 | +1.06 | - -The catalogue Mw 1.6 sits *between* the metre and kilometre answers, 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 — the likeliest candidates being that `Omega` here is one horizontal -component rather than the combined horizontal (§4.7), and that the round -numbers for the medium are not this event's. - -So the calibration is the deliverable, not the derivation. 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. Two people can both be certain about the unit of `R` and both -be right about different formulations; only the calibration tells you which -one the code is in. +**The constants and the units are settled — from Holt (2019), the author's +doctoral thesis, Chapter 1 §1.4 and Chapter 2 Eq. 2.7**, which is the +Edwards et al. (2010) spectral method this package implements. Read directly +rather than inferred: + +| symbol | value | what it is | +|---|---|---| +| `rho` | 2600 kg/m^3 (Utah); 2800 generic | density **at the source** | +| `beta` | 3500 m/s generic | S velocity **at the source** | +| `R_0` | **1000 m** | *reference source distance* — the distance at which the source spectrum is defined | +| `F` | 2 | free surface, for **vertically incident SH** | +| `Theta-lambda-Phi` | **0.55** | average radiation pattern **of SH propagation** over the focal sphere (Boatwright, 1978) | +| partition factor | (Boore, 2003) | splits energy between vertical and horizontal ground motion | + +Two things this corrects in earlier drafts of this section. + +**`0.55` is not a different average of the same quantity — it is a different +quantity.** It is the **SH** radiation pattern, 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 actually being measured — and that couples straight back +to the rotation question above: a formulation using the SH radiation pattern +wants the transverse component, not an arbitrary horizontal. + +**`R_0 = 1000 m` is a reference distance, not a unit bridge.** An earlier +draft of this document guessed it might be bridging kilometres and metres. It +is not; it is where the source spectrum is defined, and the spreading model +carries the observation from there to the site. + +**And that is why both answers about the unit of `R` are right.** There are two +distances. `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, e.g. Holt et al. +[O] for Utah: + + 0.88 +/- 0.02 1 < R <= 40 km + 2.93 +/- 0.28 40 < R <= 63 km + 0.50 +/- 0.33 63 < R <= 100 km + 1.36 +/- 0.07 100 < R <= 400 km + +Worth noting the first segment independently supports the exponent argument +above: the inverted near-field decay is **0.88**, close to the theoretical +body-wave 1, and nowhere near 2. + +It also shows the spreading is not a single power law but a piecewise +empirical function, inverted per region. So `S(R)` should be a registered +model in its own right — theoretical `1/R` as the default a study can start +from, with region-specific piecewise tables as the thing an operator supplies, +in the same shape as `WEIGHT_MODELS` and `NOISE_MODELS`. + +The calibration is still the deliverable rather than the derivation. 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. This section is the case for it: two people can each +be certain about "the unit of R" and both be right, because there are two +different `R`s in the same workflow, and only a test that runs the numbers +end to end distinguishes them. + +*Sourced from Holt (2019), "Addressing Uncertainty in Earthquake Magnitudes +Commonly Used in Modern Seismic Hazard Assessment", University of Liverpool — +Ch. 1 §1.4 and Ch. 2 Eq. 2.7-2.8 and Table 2.1. The equations themselves are +embedded objects that do not survive text extraction, so the symbols above are +read from the surrounding prose; the equation images should be checked against +this table before any of it is implemented.* Two other things that must be true before the number means anything, both from §4.7 above: `Omega` should be the combined horizontal rather than one From 73f630a5d101bf7f9260a9b5fdd2fda4fddebb6d Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 8 Aug 2026 16:40:37 +0000 Subject: [PATCH 4/9] docs(plan): keep Magna as the calibration anchor, and say why MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- docs/REFACTOR_PLAN.md | 27 +++++++++++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) diff --git a/docs/REFACTOR_PLAN.md b/docs/REFACTOR_PLAN.md index d499aff..ff57fff 100644 --- a/docs/REFACTOR_PLAN.md +++ b/docs/REFACTOR_PLAN.md @@ -1418,8 +1418,31 @@ body-wave 1, and nowhere near 2. It also shows the spreading is not a single power law but a piecewise empirical function, inverted per region. So `S(R)` should be a registered model in its own right — theoretical `1/R` as the default a study can start -from, with region-specific piecewise tables as the thing an operator supplies, -in the same shape as `WEIGHT_MODELS` and `NOISE_MODELS`. +from, with region-specific spreading as the thing an operator supplies, in the +same shape as `WEIGHT_MODELS` and `NOISE_MODELS`. + +**The calibration anchor stays Magna, not the thesis's Utah catalogue.** Utah +is tempting — it is a published catalogue of >200 `Mw` with regional constants +and bilinear `Mw`-`ML` relationships attached, which is more validation +material than Magna offers. It is not the right anchor anyway, because the +spreading estimate is the weaker half of it: Utah's `S(R)` is the piecewise +parametric table above, with segment boundaries chosen and an exponent fitted +per segment, while Magna used the **non-parametric G(R) inversion** already +noted in §5.2.5 — no imposed functional form and no chosen breakpoints. + +That ordering matters more than the size of the catalogue, because `S(R)` is +the term this section is least able to check by other means. `rho`, `beta` and +the radiation pattern are all bounded by physics and by the literature; a +spreading function is bounded by the inversion that produced it. Validating +against the dataset with the better-constrained `S(R)` tests the part of the +formula that most needs testing. + +The consequence for scope is worth stating plainly: the non-parametric +inversion is not in this repository (§5.2.5), so `S(R)` for the Magna +comparison has to come in as **data** — a tabulated `G(R)` in +`studies/magna_2020.toml` — rather than being recomputed. A spreading registry +that can be fed a table, not only a functional form, is therefore a +requirement rather than a nicety. The calibration is still the deliverable rather than the derivation. Pin the study's own `rho` and `c`, state the unit of every input, and assert the Mw of From 9ce99b78f225aaf9046536ba370f78ba017d0542 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 8 Aug 2026 16:43:54 +0000 Subject: [PATCH 5/9] docs(plan): bilinear spreading belongs in the design, but one event cannot fit it MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- docs/REFACTOR_PLAN.md | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/docs/REFACTOR_PLAN.md b/docs/REFACTOR_PLAN.md index ff57fff..7570f59 100644 --- a/docs/REFACTOR_PLAN.md +++ b/docs/REFACTOR_PLAN.md @@ -1421,6 +1421,41 @@ model in its own right — theoretical `1/R` as the default a study can start from, with region-specific spreading as the thing an operator supplies, in the same shape as `WEIGHT_MODELS` and `NOISE_MODELS`. +**Bilinear, and why the microseismic regime is not the regional one.** Every +spreading model in the thesis breaks at 40-50 km. The PNR data used throughout +this repository spans **2.3 to 22.9 km** — entirely inside the first segment of +all of them — so the regional tables offer no guidance at these distances, and +a shape fitted for 1-400 km should not be assumed to extend downward. A +bilinear form with a break inside the microseismic range is the right thing to +allow. + +It is not something a single event can determine, though, and it is worth +recording why rather than fitting one and believing it. Measured on the 28 PNR +channels, regressing `log10(Omega)` on `log10(R)` — which for one event is +exactly the spreading, since `Omega_source` is a constant: + +| model | exponents | rms residual | +|---|---|---| +| single segment | 0.73 | 0.265 log10 | +| bilinear, break searched | 0.41 then 1.89, break 13.5 km | 0.253 log10 | + +The bilinear fit buys **5%** in rms for two extra parameters, across **one +decade** of distance. That is not evidence of a break; it is a hinge finding +scatter. And the scatter is the point: 0.265 log10 units is a factor of 1.8, +which for a single event is site response and radiation pattern, neither of +which is separable from spreading when every station contributes exactly one +distance. + +That is the argument for the non-parametric inversion rather than against +bilinear. Spreading becomes separable from site only across a dataset where +each station sees many distances and each distance is sampled by many +stations. One event cannot do it, and the apparent 0.73 here should not be +read as a measurement of anything. + +So: 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. + **The calibration anchor stays Magna, not the thesis's Utah catalogue.** Utah is tempting — it is a published catalogue of >200 `Mw` with regional constants and bilinear `Mw`-`ML` relationships attached, which is more validation From 782760cbf44800b1f71711d65dcce81f3da2303d Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 8 Aug 2026 16:52:00 +0000 Subject: [PATCH 6/9] feat(distance): make the source-to-site distance a registry, and read it MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `[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. --- Tutorial/SpecModTutorial.ipynb | 495 +++++++++++++++++---------------- docs/REFACTOR_PLAN.md | 31 ++- src/specmod/config/sections.py | 2 +- src/specmod/distance.py | 189 +++++++++++++ src/specmod/staged.py | 43 ++- tests/test_staged.py | 60 +++- 6 files changed, 545 insertions(+), 275 deletions(-) create mode 100644 src/specmod/distance.py diff --git a/Tutorial/SpecModTutorial.ipynb b/Tutorial/SpecModTutorial.ipynb index 38c4665..f86d11d 100644 --- a/Tutorial/SpecModTutorial.ipynb +++ b/Tutorial/SpecModTutorial.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "39eda767", + "id": "80661db4", "metadata": {}, "source": [ "# SpecMod tutorial\n", @@ -29,7 +29,7 @@ }, { "cell_type": "markdown", - "id": "59f58548", + "id": "01a3ea27", "metadata": {}, "source": [ "## 1. Read and prepare" @@ -38,13 +38,13 @@ { "cell_type": "code", "execution_count": 1, - "id": "e3afbe8f", + "id": "b06e3ef2", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:20.194120Z", - "iopub.status.busy": "2026-08-08T16:05:20.193836Z", - "iopub.status.idle": "2026-08-08T16:05:21.133218Z", - "shell.execute_reply": "2026-08-08T16:05:21.131949Z" + "iopub.execute_input": "2026-08-08T16:49:51.779439Z", + "iopub.status.busy": "2026-08-08T16:49:51.779238Z", + "iopub.status.idle": "2026-08-08T16:49:53.063613Z", + "shell.execute_reply": "2026-08-08T16:49:53.062068Z" } }, "outputs": [], @@ -65,7 +65,7 @@ }, { "cell_type": "markdown", - "id": "4c94213f", + "id": "6606fad9", "metadata": {}, "source": [ "The origin is what distances and theoretical arrivals are measured from.\n", @@ -82,13 +82,13 @@ { "cell_type": "code", "execution_count": 2, - "id": "f1bc897f", + "id": "864322f1", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:21.135876Z", - "iopub.status.busy": "2026-08-08T16:05:21.135432Z", - "iopub.status.idle": "2026-08-08T16:05:21.139895Z", - "shell.execute_reply": "2026-08-08T16:05:21.138517Z" + "iopub.execute_input": "2026-08-08T16:49:53.066533Z", + "iopub.status.busy": "2026-08-08T16:49:53.066048Z", + "iopub.status.idle": "2026-08-08T16:49:53.071135Z", + "shell.execute_reply": "2026-08-08T16:49:53.069373Z" } }, "outputs": [], @@ -100,13 +100,13 @@ { "cell_type": "code", "execution_count": 3, - "id": "732e9a38", + "id": "b4e74062", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:21.141881Z", - "iopub.status.busy": "2026-08-08T16:05:21.141655Z", - "iopub.status.idle": "2026-08-08T16:05:21.243229Z", - "shell.execute_reply": "2026-08-08T16:05:21.242244Z" + "iopub.execute_input": "2026-08-08T16:49:53.073283Z", + "iopub.status.busy": "2026-08-08T16:49:53.073067Z", + "iopub.status.idle": "2026-08-08T16:49:53.178931Z", + "shell.execute_reply": "2026-08-08T16:49:53.177544Z" } }, "outputs": [ @@ -131,13 +131,13 @@ { "cell_type": "code", "execution_count": 4, - "id": "86bd514a", + "id": "7f37be38", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:21.245322Z", - "iopub.status.busy": "2026-08-08T16:05:21.244903Z", - "iopub.status.idle": "2026-08-08T16:05:21.254042Z", - "shell.execute_reply": "2026-08-08T16:05:21.252921Z" + "iopub.execute_input": "2026-08-08T16:49:53.181679Z", + "iopub.status.busy": "2026-08-08T16:49:53.181107Z", + "iopub.status.idle": "2026-08-08T16:49:53.191799Z", + "shell.execute_reply": "2026-08-08T16:49:53.190749Z" } }, "outputs": [ @@ -176,7 +176,7 @@ }, { "cell_type": "markdown", - "id": "fac2ce5b", + "id": "fbef660f", "metadata": {}, "source": [ "**Instrument correction happens here, in the notebook, not inside the\n", @@ -192,13 +192,13 @@ { "cell_type": "code", "execution_count": 5, - "id": "80933dcc", + "id": "df25562f", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:21.256143Z", - "iopub.status.busy": "2026-08-08T16:05:21.255921Z", - "iopub.status.idle": "2026-08-08T16:05:21.765679Z", - "shell.execute_reply": "2026-08-08T16:05:21.764595Z" + "iopub.execute_input": "2026-08-08T16:49:53.194187Z", + "iopub.status.busy": "2026-08-08T16:49:53.193855Z", + "iopub.status.idle": "2026-08-08T16:49:53.782359Z", + "shell.execute_reply": "2026-08-08T16:49:53.781148Z" } }, "outputs": [ @@ -239,13 +239,13 @@ { "cell_type": "code", "execution_count": 6, - "id": "a4af3926", + "id": "396193c1", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:21.767678Z", - "iopub.status.busy": "2026-08-08T16:05:21.767471Z", - "iopub.status.idle": "2026-08-08T16:05:23.302059Z", - "shell.execute_reply": "2026-08-08T16:05:23.300791Z" + "iopub.execute_input": "2026-08-08T16:49:53.784573Z", + "iopub.status.busy": "2026-08-08T16:49:53.784327Z", + "iopub.status.idle": "2026-08-08T16:49:55.334948Z", + "shell.execute_reply": "2026-08-08T16:49:55.333489Z" } }, "outputs": [], @@ -256,7 +256,7 @@ }, { "cell_type": "markdown", - "id": "23e4c246", + "id": "1286883e", "metadata": {}, "source": [ "## 2. Cut the windows" @@ -264,7 +264,7 @@ }, { "cell_type": "markdown", - "id": "a9667906", + "id": "aeb08889", "metadata": {}, "source": [ "The S-window opens at a fixed fraction of the elapsed P–S time after the P\n", @@ -280,13 +280,13 @@ { "cell_type": "code", "execution_count": 7, - "id": "56e0fd1c", + "id": "45304586", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:23.305017Z", - "iopub.status.busy": "2026-08-08T16:05:23.304814Z", - "iopub.status.idle": "2026-08-08T16:05:23.322224Z", - "shell.execute_reply": "2026-08-08T16:05:23.321012Z" + "iopub.execute_input": "2026-08-08T16:49:55.337432Z", + "iopub.status.busy": "2026-08-08T16:49:55.337141Z", + "iopub.status.idle": "2026-08-08T16:49:55.358539Z", + "shell.execute_reply": "2026-08-08T16:49:55.357428Z" } }, "outputs": [], @@ -300,13 +300,13 @@ { "cell_type": "code", "execution_count": 8, - "id": "a4fe9b32", + "id": "720530dd", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:23.324105Z", - "iopub.status.busy": "2026-08-08T16:05:23.323905Z", - "iopub.status.idle": "2026-08-08T16:05:23.328110Z", - "shell.execute_reply": "2026-08-08T16:05:23.327014Z" + "iopub.execute_input": "2026-08-08T16:49:55.361393Z", + "iopub.status.busy": "2026-08-08T16:49:55.361161Z", + "iopub.status.idle": "2026-08-08T16:49:55.366743Z", + "shell.execute_reply": "2026-08-08T16:49:55.365013Z" } }, "outputs": [ @@ -338,13 +338,13 @@ { "cell_type": "code", "execution_count": 9, - "id": "e9929541", + "id": "9381e25c", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:23.329901Z", - "iopub.status.busy": "2026-08-08T16:05:23.329699Z", - "iopub.status.idle": "2026-08-08T16:05:25.155951Z", - "shell.execute_reply": "2026-08-08T16:05:25.154634Z" + "iopub.execute_input": "2026-08-08T16:49:55.368922Z", + "iopub.status.busy": "2026-08-08T16:49:55.368732Z", + "iopub.status.idle": "2026-08-08T16:49:56.997523Z", + "shell.execute_reply": "2026-08-08T16:49:56.996091Z" } }, "outputs": [], @@ -355,7 +355,7 @@ }, { "cell_type": "markdown", - "id": "a8dd38a4", + "id": "57733e86", "metadata": {}, "source": [ "## 3. Spectra and bandwidth" @@ -363,7 +363,7 @@ }, { "cell_type": "markdown", - "id": "e3fee1df", + "id": "08db6076", "metadata": {}, "source": [ "`spectrum_set_from_streams` transforms both windows, puts the noise on the\n", @@ -378,13 +378,13 @@ { "cell_type": "code", "execution_count": 10, - "id": "87aebed5", + "id": "0c1fa24f", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:25.158367Z", - "iopub.status.busy": "2026-08-08T16:05:25.158112Z", - "iopub.status.idle": "2026-08-08T16:05:25.355124Z", - "shell.execute_reply": "2026-08-08T16:05:25.353938Z" + "iopub.execute_input": "2026-08-08T16:49:57.000245Z", + "iopub.status.busy": "2026-08-08T16:49:56.999975Z", + "iopub.status.idle": "2026-08-08T16:49:57.201338Z", + "shell.execute_reply": "2026-08-08T16:49:57.199880Z" } }, "outputs": [ @@ -408,13 +408,13 @@ { "cell_type": "code", "execution_count": 11, - "id": "b93a32ef", + "id": "1e547105", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:25.357445Z", - "iopub.status.busy": "2026-08-08T16:05:25.357245Z", - "iopub.status.idle": "2026-08-08T16:05:25.377656Z", - "shell.execute_reply": "2026-08-08T16:05:25.376409Z" + "iopub.execute_input": "2026-08-08T16:49:57.204138Z", + "iopub.status.busy": "2026-08-08T16:49:57.203843Z", + "iopub.status.idle": "2026-08-08T16:49:57.231532Z", + "shell.execute_reply": "2026-08-08T16:49:57.230099Z" } }, "outputs": [], @@ -431,13 +431,13 @@ { "cell_type": "code", "execution_count": 12, - "id": "a93bce39", + "id": "5af0fd70", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:25.379784Z", - "iopub.status.busy": "2026-08-08T16:05:25.379584Z", - "iopub.status.idle": "2026-08-08T16:05:25.383672Z", - "shell.execute_reply": "2026-08-08T16:05:25.382549Z" + "iopub.execute_input": "2026-08-08T16:49:57.234367Z", + "iopub.status.busy": "2026-08-08T16:49:57.234148Z", + "iopub.status.idle": "2026-08-08T16:49:57.239171Z", + "shell.execute_reply": "2026-08-08T16:49:57.237476Z" } }, "outputs": [ @@ -460,7 +460,7 @@ }, { "cell_type": "markdown", - "id": "dd6078f8", + "id": "cb81de32", "metadata": {}, "source": [ "### Changing ground-motion domain\n", @@ -473,13 +473,13 @@ { "cell_type": "code", "execution_count": 13, - "id": "bc8fef2a", + "id": "508d3f73", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:25.385510Z", - "iopub.status.busy": "2026-08-08T16:05:25.385339Z", - "iopub.status.idle": "2026-08-08T16:05:25.424633Z", - "shell.execute_reply": "2026-08-08T16:05:25.423478Z" + "iopub.execute_input": "2026-08-08T16:49:57.241350Z", + "iopub.status.busy": "2026-08-08T16:49:57.241166Z", + "iopub.status.idle": "2026-08-08T16:49:57.280015Z", + "shell.execute_reply": "2026-08-08T16:49:57.278838Z" } }, "outputs": [ @@ -507,7 +507,7 @@ }, { "cell_type": "markdown", - "id": "1d8f7fea", + "id": "75de79d6", "metadata": {}, "source": [ "## 4. Fit a source model" @@ -515,7 +515,7 @@ }, { "cell_type": "markdown", - "id": "3dfca7b2", + "id": "b92f4418", "metadata": {}, "source": [ "The model comes from configuration — a Brune source with constant Q by\n", @@ -532,13 +532,13 @@ { "cell_type": "code", "execution_count": 14, - "id": "b5ac93c1", + "id": "000db0ef", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:25.426632Z", - "iopub.status.busy": "2026-08-08T16:05:25.426437Z", - "iopub.status.idle": "2026-08-08T16:05:25.477264Z", - "shell.execute_reply": "2026-08-08T16:05:25.476242Z" + "iopub.execute_input": "2026-08-08T16:49:57.282340Z", + "iopub.status.busy": "2026-08-08T16:49:57.282151Z", + "iopub.status.idle": "2026-08-08T16:49:57.318987Z", + "shell.execute_reply": "2026-08-08T16:49:57.317867Z" } }, "outputs": [ @@ -563,13 +563,13 @@ { "cell_type": "code", "execution_count": 15, - "id": "f1e68721", + "id": "6abdc313", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:25.479686Z", - "iopub.status.busy": "2026-08-08T16:05:25.479095Z", - "iopub.status.idle": "2026-08-08T16:05:26.383789Z", - "shell.execute_reply": "2026-08-08T16:05:26.382633Z" + "iopub.execute_input": "2026-08-08T16:49:57.322514Z", + "iopub.status.busy": "2026-08-08T16:49:57.321509Z", + "iopub.status.idle": "2026-08-08T16:49:58.166284Z", + "shell.execute_reply": "2026-08-08T16:49:58.164994Z" } }, "outputs": [ @@ -589,7 +589,7 @@ }, { "cell_type": "markdown", - "id": "6a7eb72f", + "id": "d19e4f9b", "metadata": {}, "source": [ "The guess is only a starting point, and a crude one: it takes the largest\n", @@ -602,13 +602,13 @@ { "cell_type": "code", "execution_count": 16, - "id": "03f89e94", + "id": "f33f5eda", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:26.385827Z", - "iopub.status.busy": "2026-08-08T16:05:26.385633Z", - "iopub.status.idle": "2026-08-08T16:05:26.401731Z", - "shell.execute_reply": "2026-08-08T16:05:26.400676Z" + "iopub.execute_input": "2026-08-08T16:49:58.168880Z", + "iopub.status.busy": "2026-08-08T16:49:58.168643Z", + "iopub.status.idle": "2026-08-08T16:49:58.188051Z", + "shell.execute_reply": "2026-08-08T16:49:58.186792Z" } }, "outputs": [ @@ -725,13 +725,13 @@ { "cell_type": "code", "execution_count": 17, - "id": "15af287a", + "id": "095d9120", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:26.403508Z", - "iopub.status.busy": "2026-08-08T16:05:26.403329Z", - "iopub.status.idle": "2026-08-08T16:05:26.421423Z", - "shell.execute_reply": "2026-08-08T16:05:26.420372Z" + "iopub.execute_input": "2026-08-08T16:49:58.190521Z", + "iopub.status.busy": "2026-08-08T16:49:58.190301Z", + "iopub.status.idle": "2026-08-08T16:49:58.210382Z", + "shell.execute_reply": "2026-08-08T16:49:58.208850Z" } }, "outputs": [], @@ -744,13 +744,13 @@ { "cell_type": "code", "execution_count": 18, - "id": "704b062a", + "id": "514a66be", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:26.423709Z", - "iopub.status.busy": "2026-08-08T16:05:26.423496Z", - "iopub.status.idle": "2026-08-08T16:05:30.421284Z", - "shell.execute_reply": "2026-08-08T16:05:30.419781Z" + "iopub.execute_input": "2026-08-08T16:49:58.212770Z", + "iopub.status.busy": "2026-08-08T16:49:58.212534Z", + "iopub.status.idle": "2026-08-08T16:50:01.491642Z", + "shell.execute_reply": "2026-08-08T16:50:01.490434Z" } }, "outputs": [], @@ -763,13 +763,13 @@ { "cell_type": "code", "execution_count": 19, - "id": "c8d5cdd1", + "id": "54bf2daa", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:30.424820Z", - "iopub.status.busy": "2026-08-08T16:05:30.424540Z", - "iopub.status.idle": "2026-08-08T16:05:30.439861Z", - "shell.execute_reply": "2026-08-08T16:05:30.438685Z" + "iopub.execute_input": "2026-08-08T16:50:01.493730Z", + "iopub.status.busy": "2026-08-08T16:50:01.493496Z", + "iopub.status.idle": "2026-08-08T16:50:01.505023Z", + "shell.execute_reply": "2026-08-08T16:50:01.503985Z" } }, "outputs": [ @@ -922,7 +922,7 @@ }, { "cell_type": "markdown", - "id": "4f76ba73", + "id": "1d7cdb71", "metadata": {}, "source": [ "### The fit is not unique, and that is not a detail\n", @@ -946,13 +946,13 @@ { "cell_type": "code", "execution_count": 20, - "id": "4c83be9e", + "id": "9b806a45", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:30.442104Z", - "iopub.status.busy": "2026-08-08T16:05:30.441839Z", - "iopub.status.idle": "2026-08-08T16:05:31.463489Z", - "shell.execute_reply": "2026-08-08T16:05:31.462097Z" + "iopub.execute_input": "2026-08-08T16:50:01.507171Z", + "iopub.status.busy": "2026-08-08T16:50:01.506985Z", + "iopub.status.idle": "2026-08-08T16:50:02.470569Z", + "shell.execute_reply": "2026-08-08T16:50:02.469397Z" } }, "outputs": [ @@ -1086,7 +1086,7 @@ }, { "cell_type": "markdown", - "id": "0424928a", + "id": "c2c148f3", "metadata": {}, "source": [ "Most stations agree to a fraction of a percent. A few do not, and the top of\n", @@ -1102,13 +1102,13 @@ { "cell_type": "code", "execution_count": 21, - "id": "372771f9", + "id": "493d0dea", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:31.466011Z", - "iopub.status.busy": "2026-08-08T16:05:31.465712Z", - "iopub.status.idle": "2026-08-08T16:05:31.471262Z", - "shell.execute_reply": "2026-08-08T16:05:31.470165Z" + "iopub.execute_input": "2026-08-08T16:50:02.472588Z", + "iopub.status.busy": "2026-08-08T16:50:02.472382Z", + "iopub.status.idle": "2026-08-08T16:50:02.477907Z", + "shell.execute_reply": "2026-08-08T16:50:02.476853Z" } }, "outputs": [ @@ -1139,13 +1139,13 @@ { "cell_type": "code", "execution_count": 22, - "id": "18e64e8f", + "id": "056e59ac", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:31.473285Z", - "iopub.status.busy": "2026-08-08T16:05:31.473029Z", - "iopub.status.idle": "2026-08-08T16:05:31.491080Z", - "shell.execute_reply": "2026-08-08T16:05:31.489606Z" + "iopub.execute_input": "2026-08-08T16:50:02.480019Z", + "iopub.status.busy": "2026-08-08T16:50:02.479814Z", + "iopub.status.idle": "2026-08-08T16:50:02.500292Z", + "shell.execute_reply": "2026-08-08T16:50:02.498830Z" } }, "outputs": [], @@ -1161,7 +1161,7 @@ }, { "cell_type": "markdown", - "id": "f52d1081", + "id": "171b6750", "metadata": {}, "source": [ "Look at where they differ: high up the falling limb, where the source corner\n", @@ -1186,13 +1186,13 @@ { "cell_type": "code", "execution_count": 23, - "id": "ea8f5a9b", + "id": "f4f10578", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:31.493795Z", - "iopub.status.busy": "2026-08-08T16:05:31.493595Z", - "iopub.status.idle": "2026-08-08T16:05:31.498573Z", - "shell.execute_reply": "2026-08-08T16:05:31.497367Z" + "iopub.execute_input": "2026-08-08T16:50:02.503081Z", + "iopub.status.busy": "2026-08-08T16:50:02.502866Z", + "iopub.status.idle": "2026-08-08T16:50:02.508180Z", + "shell.execute_reply": "2026-08-08T16:50:02.507142Z" } }, "outputs": [ @@ -1214,7 +1214,7 @@ }, { "cell_type": "markdown", - "id": "f10f5659", + "id": "4ce09cca", "metadata": {}, "source": [ "### Why the answer comes from many stations, and from two stages\n", @@ -1239,13 +1239,13 @@ { "cell_type": "code", "execution_count": 24, - "id": "82e8c326", + "id": "6b32068b", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:31.500823Z", - "iopub.status.busy": "2026-08-08T16:05:31.500622Z", - "iopub.status.idle": "2026-08-08T16:05:31.512549Z", - "shell.execute_reply": "2026-08-08T16:05:31.511462Z" + "iopub.execute_input": "2026-08-08T16:50:02.510699Z", + "iopub.status.busy": "2026-08-08T16:50:02.510493Z", + "iopub.status.idle": "2026-08-08T16:50:02.522429Z", + "shell.execute_reply": "2026-08-08T16:50:02.521541Z" } }, "outputs": [ @@ -1253,20 +1253,28 @@ "name": "stdout", "output_type": "stream", "text": [ - "powell event fc = 12.751 Hz\n", - "leastsq event fc = 12.675 Hz\n", + "powell event fc = 11.585 Hz\n", + "leastsq event fc = 11.535 Hz\n", "\n", "worst single station: fc ratio 1.441 -> stress drop 2.99x\n", - "event ensemble: fc ratio 1.006 -> stress drop 1.018x\n" + "event ensemble: fc ratio 1.004 -> stress drop 1.013x\n" ] } ], "source": [ - "# Weighted by inverse hypocentral distance: the nearer station has less path\n", - "# between the source and the sensor, so less of its high-frequency falloff can\n", - "# be attenuation, and its corner is the better constrained of the two.\n", + "# Weighted by inverse distance: the nearer station has less path between the\n", + "# source and the sensor, so less of its high-frequency falloff can be\n", + "# attenuation, and its corner is the better constrained of the two.\n", + "#\n", + "# `repi` — epicentral — because that is what `[windows] distance_metric` says.\n", + "# Which distance you use is not a detail at short range: here the nearest\n", + "# station is 0.89 km epicentral against 2.30 km hypocentral. `rhyp` is built\n", + "# from the source depth and the station *elevation*, so it assumes every\n", + "# sensor is at the surface; where sensor depths are unknown, as they are here,\n", + "# epicentral is the honest choice.\n", + "metric = \"repi\"\n", "weight = pd.Series(\n", - " {id: 1.0 / spectra[id].signal.meta[\"rhyp\"] for id in fits.models}\n", + " {id: 1.0 / spectra[id].signal.meta[metric] for id in fits.models}\n", ")\n", "\n", "event_fc = {}\n", @@ -1287,7 +1295,7 @@ }, { "cell_type": "markdown", - "id": "be9c2b9a", + "id": "c308fef1", "metadata": {}, "source": [ "A factor of three in stress drop at the worst station becomes **under 2%**\n", @@ -1302,13 +1310,13 @@ { "cell_type": "code", "execution_count": 25, - "id": "d0e5ba30", + "id": "d67f1cef", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:31.514662Z", - "iopub.status.busy": "2026-08-08T16:05:31.514484Z", - "iopub.status.idle": "2026-08-08T16:05:32.299622Z", - "shell.execute_reply": "2026-08-08T16:05:32.298368Z" + "iopub.execute_input": "2026-08-08T16:50:02.524750Z", + "iopub.status.busy": "2026-08-08T16:50:02.524533Z", + "iopub.status.idle": "2026-08-08T16:50:03.295928Z", + "shell.execute_reply": "2026-08-08T16:50:03.294940Z" } }, "outputs": [ @@ -1316,8 +1324,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "t* worst disagreement between minimisers: 0.31%\n", - "omega worst disagreement between minimisers: 2.3e-03 log10 units\n" + "t* worst disagreement between minimisers: 0.23%\n", + "omega worst disagreement between minimisers: 1.8e-03 log10 units\n" ] }, { @@ -1357,45 +1365,45 @@ " \n", " \n", " LV.L001..HHE\n", - " 0.0338\n", - " 0.0337\n", - " -4.9827\n", - " -4.9818\n", + " 0.0326\n", + " 0.0326\n", + " -4.9657\n", + " -4.9649\n", " \n", " \n", " LV.L002..HHE\n", - " 0.0547\n", - " 0.0546\n", - " -5.3057\n", - " -5.3053\n", + " 0.0529\n", + " 0.0529\n", + " -5.2984\n", + " -5.2979\n", " \n", " \n", " LV.L006..HHE\n", - " 0.0834\n", - " 0.0833\n", - " -5.1199\n", - " -5.1193\n", + " 0.0820\n", + " 0.0819\n", + " -5.1084\n", + " -5.1078\n", " \n", " \n", " LV.L007..HHE\n", - " 0.0774\n", - " 0.0772\n", - " -5.6150\n", - " -5.6150\n", + " 0.0752\n", + " 0.0751\n", + " -5.6139\n", + " -5.6138\n", " \n", " \n", " LV.L008..HHE\n", - " 0.0559\n", - " 0.0557\n", - " -5.8243\n", - " -5.8246\n", + " 0.0533\n", + " 0.0532\n", + " -5.8285\n", + " -5.8286\n", " \n", " \n", " LV.L009..HHE\n", - " 0.0439\n", - " 0.0439\n", - " -5.4172\n", - " -5.4162\n", + " 0.0427\n", + " 0.0427\n", + " -5.4009\n", + " -5.4001\n", " \n", " \n", "\n", @@ -1404,12 +1412,12 @@ "text/plain": [ " t* powell t* leastsq log10 omega powell log10 omega leastsq\n", "id \n", - "LV.L001..HHE 0.0338 0.0337 -4.9827 -4.9818\n", - "LV.L002..HHE 0.0547 0.0546 -5.3057 -5.3053\n", - "LV.L006..HHE 0.0834 0.0833 -5.1199 -5.1193\n", - "LV.L007..HHE 0.0774 0.0772 -5.6150 -5.6150\n", - "LV.L008..HHE 0.0559 0.0557 -5.8243 -5.8246\n", - "LV.L009..HHE 0.0439 0.0439 -5.4172 -5.4162" + "LV.L001..HHE 0.0326 0.0326 -4.9657 -4.9649\n", + "LV.L002..HHE 0.0529 0.0529 -5.2984 -5.2979\n", + "LV.L006..HHE 0.0820 0.0819 -5.1084 -5.1078\n", + "LV.L007..HHE 0.0752 0.0751 -5.6139 -5.6138\n", + "LV.L008..HHE 0.0533 0.0532 -5.8285 -5.8286\n", + "LV.L009..HHE 0.0427 0.0427 -5.4009 -5.4001" ] }, "execution_count": 25, @@ -1441,7 +1449,7 @@ }, { "cell_type": "markdown", - "id": "7f80cf74", + "id": "7efdbc23", "metadata": {}, "source": [ "0.3% in $t^*$ and about 0.002 in $\\log_{10}\\Omega$ — which is 0.003 magnitude\n", @@ -1463,7 +1471,7 @@ }, { "cell_type": "markdown", - "id": "4ab068b5", + "id": "9bb2ff9c", "metadata": {}, "source": [ "### The same thing, as one call\n", @@ -1482,13 +1490,13 @@ { "cell_type": "code", "execution_count": 26, - "id": "0082031d", + "id": "34091db5", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:32.301740Z", - "iopub.status.busy": "2026-08-08T16:05:32.301520Z", - "iopub.status.idle": "2026-08-08T16:05:33.379213Z", - "shell.execute_reply": "2026-08-08T16:05:33.377893Z" + "iopub.execute_input": "2026-08-08T16:50:03.298051Z", + "iopub.status.busy": "2026-08-08T16:50:03.297784Z", + "iopub.status.idle": "2026-08-08T16:50:04.380268Z", + "shell.execute_reply": "2026-08-08T16:50:04.378693Z" } }, "outputs": [ @@ -1496,11 +1504,11 @@ "name": "stdout", "output_type": "stream", "text": [ - "fc = 12.75 from 28 channels, weighted by inverse_distance\n", - " stage-1 range 2.835 to 59.64 (445% of the event value)\n", + "fc = 11.59 from 28 channels, weighted by inverse_distance\n", + " stage-1 range 2.835 to 59.64 (490% of the event value)\n", "\n", - "by hand : 12.7512 Hz\n", - "API : 12.7512 Hz\n", + "by hand : 11.5851 Hz\n", + "API : 11.5851 Hz\n", "agree : True\n" ] } @@ -1518,7 +1526,7 @@ }, { "cell_type": "markdown", - "id": "a1c9c1a4", + "id": "a8ef817f", "metadata": {}, "source": [ "`describe()` prints the spread as well as the mean, and that is deliberate. A\n", @@ -1534,21 +1542,26 @@ "wrong, and averaging it in moves the event value for every other station.\n", "\n", "Two things make that lever bigger than it first looks. Inverse-distance\n", - "weighting is concentrated — the nearest two channels carry about a fifth of\n", - "the total weight — and stress drop goes as $f_c^3$, so a modest shift in the\n", - "corner is a large shift in the thing you are reporting." + "weighting is concentrated — on epicentral distance here the nearest two\n", + "channels carry over 40% of the total weight — and stress drop goes as $f_c^3$,\n", + "so a modest shift in the corner is a large shift in the thing you report.\n", + "\n", + "Which distance measure you choose feeds straight into this, and\n", + "`specmod.distance` makes it a registry for that reason: `repi` and `rhyp` are\n", + "implemented, and `rrup`/`rjb` are registered but raise, since a point source\n", + "has no rupture surface to measure from." ] }, { "cell_type": "code", "execution_count": 27, - "id": "23a48dfd", + "id": "4d563fb6", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:33.381166Z", - "iopub.status.busy": "2026-08-08T16:05:33.380995Z", - "iopub.status.idle": "2026-08-08T16:05:33.386019Z", - "shell.execute_reply": "2026-08-08T16:05:33.384855Z" + "iopub.execute_input": "2026-08-08T16:50:04.382465Z", + "iopub.status.busy": "2026-08-08T16:50:04.382271Z", + "iopub.status.idle": "2026-08-08T16:50:04.387143Z", + "shell.execute_reply": "2026-08-08T16:50:04.385991Z" } }, "outputs": [ @@ -1557,12 +1570,12 @@ "output_type": "stream", "text": [ "weight carried by the nearest channels:\n", - " nearest 1: 11.4%\n", - " nearest 2: 22.8%\n", - " nearest 4: 36.1%\n", - " nearest 8: 56.3%\n", + " nearest 1: 20.7%\n", + " nearest 2: 41.4%\n", + " nearest 4: 52.5%\n", + " nearest 8: 68.2%\n", "\n", - "the single nearest is UR.AQ04.00.HHN at 2.30 km\n" + "the single nearest is UR.AQ04.00.HHN at 0.89 km (repi)\n" ] } ], @@ -1570,7 +1583,7 @@ "import numpy as np\n", "\n", "ids = list(staged.contributing)\n", - "distance = np.array([spectra[i].signal.meta[\"rhyp\"] for i in ids])\n", + "distance = np.array([spectra[i].signal.meta[metric] for i in ids])\n", "w = 1 / distance\n", "w = w / w.sum()\n", "order = np.argsort(-w)\n", @@ -1579,19 +1592,19 @@ "for k in (1, 2, 4, 8):\n", " print(f\" nearest {k:2d}: {100 * w[order[:k]].sum():5.1f}%\")\n", "print()\n", - "print(f\"the single nearest is {ids[order[0]]} at {distance[order[0]]:.2f} km\")" + "print(f\"the single nearest is {ids[order[0]]} at {distance[order[0]]:.2f} km ({metric})\")" ] }, { "cell_type": "code", "execution_count": 28, - "id": "4ad2371c", + "id": "add4bb15", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:33.387853Z", - "iopub.status.busy": "2026-08-08T16:05:33.387680Z", - "iopub.status.idle": "2026-08-08T16:05:34.417373Z", - "shell.execute_reply": "2026-08-08T16:05:34.416061Z" + "iopub.execute_input": "2026-08-08T16:50:04.389289Z", + "iopub.status.busy": "2026-08-08T16:50:04.389072Z", + "iopub.status.idle": "2026-08-08T16:50:05.666198Z", + "shell.execute_reply": "2026-08-08T16:50:05.664533Z" } }, "outputs": [ @@ -1599,10 +1612,10 @@ "name": "stdout", "output_type": "stream", "text": [ - "all channels fc = 12.751 Hz (28 channels)\n", - "without AQ04 fc = 14.774 Hz (26 channels)\n", - " change in fc +15.9%\n", - " change in stress drop 1.56x\n", + "all channels fc = 11.585 Hz (28 channels)\n", + "without AQ04 fc = 15.602 Hz (26 channels)\n", + " change in fc +34.7%\n", + " change in stress drop 2.44x\n", "\n", " UR.AQ04.00.HHE: matched exclude='AQ04' at station\n", " UR.AQ04.00.HHN: matched exclude='AQ04' at station\n" @@ -1630,7 +1643,7 @@ }, { "cell_type": "markdown", - "id": "02dddca6", + "id": "3f9b319d", "metadata": {}, "source": [ "One quality-control decision, a factor of 1.5 in stress drop. That is not an\n", @@ -1651,13 +1664,13 @@ { "cell_type": "code", "execution_count": 29, - "id": "79a89ac6", + "id": "ba244a42", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:34.419665Z", - "iopub.status.busy": "2026-08-08T16:05:34.419462Z", - "iopub.status.idle": "2026-08-08T16:05:37.185366Z", - "shell.execute_reply": "2026-08-08T16:05:37.184045Z" + "iopub.execute_input": "2026-08-08T16:50:05.668903Z", + "iopub.status.busy": "2026-08-08T16:50:05.668683Z", + "iopub.status.idle": "2026-08-08T16:50:08.758298Z", + "shell.execute_reply": "2026-08-08T16:50:08.756603Z" } }, "outputs": [ @@ -1665,14 +1678,14 @@ "name": "stdout", "output_type": "stream", "text": [ - "powell 28 channels vote, event fc 12.751 Hz\n" + "powell 28 channels vote, event fc 11.585 Hz\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ - "leastsq 22 channels vote, event fc 5.228 Hz\n", + "leastsq 22 channels vote, event fc 5.407 Hz\n", "\n", "with the ensemble held fixed at all 28:\n" ] @@ -1681,14 +1694,14 @@ "name": "stdout", "output_type": "stream", "text": [ - "powell 28 channels vote, event fc 12.751 Hz\n" + "powell 28 channels vote, event fc 11.585 Hz\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ - "leastsq 28 channels vote, event fc 12.675 Hz\n" + "leastsq 28 channels vote, event fc 11.535 Hz\n" ] } ], @@ -1707,7 +1720,7 @@ }, { "cell_type": "markdown", - "id": "5fcd0afb", + "id": "71acbbd4", "metadata": {}, "source": [ "So changing the minimiser changes *which stations vote*, not only how each one\n", @@ -1723,7 +1736,7 @@ }, { "cell_type": "markdown", - "id": "a083b0a6", + "id": "54c400bd", "metadata": {}, "source": [ "## 5. Save the results" @@ -1731,7 +1744,7 @@ }, { "cell_type": "markdown", - "id": "a2aebd59", + "id": "bd03f6be", "metadata": {}, "source": [ "Two formats, because the data is used two different ways.\n", @@ -1750,13 +1763,13 @@ { "cell_type": "code", "execution_count": 30, - "id": "76c8c038", + "id": "9071b6ce", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:37.189039Z", - "iopub.status.busy": "2026-08-08T16:05:37.188751Z", - "iopub.status.idle": "2026-08-08T16:05:37.383084Z", - "shell.execute_reply": "2026-08-08T16:05:37.381766Z" + "iopub.execute_input": "2026-08-08T16:50:08.762125Z", + "iopub.status.busy": "2026-08-08T16:50:08.761787Z", + "iopub.status.idle": "2026-08-08T16:50:08.943258Z", + "shell.execute_reply": "2026-08-08T16:50:08.942038Z" } }, "outputs": [ @@ -1783,13 +1796,13 @@ { "cell_type": "code", "execution_count": 31, - "id": "935c8085", + "id": "9093db2d", "metadata": { "execution": { - "iopub.execute_input": "2026-08-08T16:05:37.386461Z", - "iopub.status.busy": "2026-08-08T16:05:37.385575Z", - "iopub.status.idle": "2026-08-08T16:05:37.425304Z", - "shell.execute_reply": "2026-08-08T16:05:37.424051Z" + "iopub.execute_input": "2026-08-08T16:50:08.946449Z", + "iopub.status.busy": "2026-08-08T16:50:08.945512Z", + "iopub.status.idle": "2026-08-08T16:50:08.983333Z", + "shell.execute_reply": "2026-08-08T16:50:08.982016Z" } }, "outputs": [ @@ -1812,7 +1825,7 @@ }, { "cell_type": "markdown", - "id": "d93d8ed8", + "id": "38f5c621", "metadata": {}, "source": [ "---\n", diff --git a/docs/REFACTOR_PLAN.md b/docs/REFACTOR_PLAN.md index 7570f59..40a06a3 100644 --- a/docs/REFACTOR_PLAN.md +++ b/docs/REFACTOR_PLAN.md @@ -2040,13 +2040,30 @@ channel. Every exclusion is recorded with a reason naming the level it matched at, because "excluded" is not actionable and "matched exclude='AQ04' at station" is. -That this matters is measured, not assumed. Under inverse hypocentral distance -weighting the nearest two channels carry 22.8% of the weight and the nearest -four carry 36.1%, so dropping the single nearest station moves the event -corner from 12.751 Hz to 14.774 Hz — 16%, which is 1.56x in stress drop. The -choice of weighting moves it too: 12.751 (inverse hypocentral), 11.585 -(inverse epicentral), 11.048 (uniform). Both are therefore registry choices -with the published one as the default, not constants. +That this matters is measured, not assumed. Under the configured weighting — +inverse *epicentral* distance — the nearest two channels carry **41.4%** of the +weight and the nearest four 52.5%, so dropping the single nearest station moves +the event corner from 11.585 Hz to 15.602 Hz: **35%, which is 2.44x in stress +drop**. The choice of weighting moves it too: 11.585 (epicentral), 12.751 +(hypocentral), 11.048 (uniform). + +**Which distance is itself configured, and `specmod.distance` is the registry +for it.** `[windows] distance_metric` existed from the start and was read by +nothing; it is read now. That is not bookkeeping at these ranges: 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 — so the two measures disagree most +exactly where the inverse-distance weight is largest. Epicentral is the honest +default here because `rhyp` is built from the source depth and the station +*elevation*, assuming every sensor sits at the surface; where sensor depths are +unknown, as they are for this deployment, that assumption is unverifiable and +the inventory's placeholder channel `depth` of `123456.0` is the tell. + +`rrup` and `rjb` are registered and **raise**. Both need a rupture surface, +and for a point source they degenerate exactly to hypocentral and epicentral — +so a silent fallback would produce plausible numbers that are wrong for any +event large enough to justify asking for them. The error names what they would +need, which puts the requirement where whoever adds finite-fault support will +read it. **One trap, found while testing and now pinned.** `require_pass` drops a station whose stage-1 fit ended against a bound. `pass_fitting` asks whether diff --git a/src/specmod/config/sections.py b/src/specmod/config/sections.py index db79519..7534103 100644 --- a/src/specmod/config/sections.py +++ b/src/specmod/config/sections.py @@ -239,7 +239,7 @@ class FittingConfig: #: How stations are weighted into the event value. The published choice is #: inverse hypocentral distance: the nearer station has less path, so less #: of its falloff can be attenuation. See ``specmod.staged.WEIGHT_MODELS``. - event_weighting: str = "inverse_hypocentral_distance" + event_weighting: str = "inverse_distance" #: Which channels contribute to the event value, as shell globs matched #: against the trace id and each of its SEED components — so ``"AQ07"`` diff --git a/src/specmod/distance.py b/src/specmod/distance.py new file mode 100644 index 0000000..2bbcce4 --- /dev/null +++ b/src/specmod/distance.py @@ -0,0 +1,189 @@ +"""Source-to-site distance, as a registry rather than a stat name. + +Which distance you mean is a modelling choice, and at short range it is not a +small one. On the PNR data the nearest station is **0.89 km epicentral against +2.30 km hypocentral** — a factor of 2.57 — while the farthest agree to 1.00. +Anything weighted by inverse distance, or corrected for geometric spreading, +therefore depends on the choice most strongly at exactly the station that +matters most. + +Two are implemented here because they are the two a point source supports. +Both read a value :func:`specmod.preprocess.set_stream_distance` has already +computed: + +``repi`` + Epicentral. Horizontal distance from the epicentre. +``rhyp`` + Hypocentral. Slant distance from the hypocentre. + +**Epicentral is the honest choice when sensor depths are unknown**, and that +is more often than it sounds. ``rhyp`` is built from the source depth and the +station *elevation*, which silently assumes every sensor sits at the surface. +For a borehole deployment that is wrong by the burial depth, and nothing in +the metadata 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. + +Finite-fault measures +--------------------- +``Rrup`` (closest distance to the rupture surface) and ``Rjb`` (Joyner-Boore, +closest horizontal distance to the surface projection of the rupture) are the +measures ground-motion work generally wants, and they are **not implemented** +— deliberately, rather than by omission. + +Both need a rupture *surface*: strike, dip, length, width and a hypocentre +position on it. SpecMod carries a point source, so there is nothing to compute +them from, and a version that quietly degenerated to ``rhyp`` and ``repi`` +would be worse than an error — those are exactly what `Rrup` and `Rjb` reduce +to for a point source, so the substitution would be invisible in the output +and wrong for any event large enough to warrant asking. + +They are registered all the same, raising with what they would need. A name +that resolves to a clear failure is a better extension point than a name that +does not resolve at all, and it puts the requirement where someone adding +finite-fault support will read it. + +The registry is the same shape as :data:`specmod.transforms.ESTIMATORS`, +:data:`specmod.core.noise.NOISE_MODELS` and +:data:`specmod.staged.WEIGHT_MODELS`, so a study names a distance the way it +names anything else and the choice travels with the resolved configuration. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import TYPE_CHECKING, Any, Protocol, runtime_checkable + +import numpy as np + +from .config import load_config + +if TYPE_CHECKING: # pragma: no cover + from collections.abc import Sequence + + from numpy.typing import NDArray + +__all__ = [ + "DISTANCE_MEASURES", + "DistanceMeasure", + "Epicentral", + "FiniteFaultDistance", + "Hypocentral", + "get_distance_measure", + "resolve_distance_measure", +] + + +@runtime_checkable +class DistanceMeasure(Protocol): + """One distance per channel, in kilometres.""" + + name: str + + def distances( + self, spectra: Any, ids: Sequence[str] + ) -> NDArray[np.float64]: ... # pragma: no cover + + +@dataclass(frozen=True, slots=True) +class _FromMeta: + """A distance already computed onto the trace metadata.""" + + key: str + name: str + + def distances(self, spectra: Any, ids: Sequence[str]) -> NDArray[np.float64]: + out = np.empty(len(ids), dtype=np.float64) + for i, id in enumerate(ids): + meta = spectra[id].signal.meta + if self.key not in meta: + raise ValueError( + f"{id} carries no {self.key!r}, so its {self.name} distance " + f"is unknown. Set the geometry with " + f"specmod.preprocess.set_stream_distance." + ) + value = float(meta[self.key]) + if value <= 0: + raise ValueError( + f"{id} has {self.key}={value}, which is not a distance" + ) + out[i] = value + return out + + +@dataclass(frozen=True, slots=True) +class Epicentral(_FromMeta): + key: str = "repi" + name: str = "epicentral" + + +@dataclass(frozen=True, slots=True) +class Hypocentral(_FromMeta): + key: str = "rhyp" + name: str = "hypocentral" + + +@dataclass(frozen=True, slots=True) +class FiniteFaultDistance: + """``Rrup`` and ``Rjb``: registered, and not implemented. + + Raising here rather than omitting the name is the point. For a point source + these degenerate exactly to hypocentral and epicentral, so an + implementation that silently fell back would produce plausible numbers that + are wrong for any event big enough to justify asking for them. + """ + + name: str + needs: str + + def distances(self, spectra: Any, ids: Sequence[str]) -> NDArray[np.float64]: + raise NotImplementedError( + f"{self.name} is not implemented. It needs {self.needs}, and " + f"SpecMod carries a point source — there is no rupture surface to " + f"measure from. For a point source {self.name} degenerates to " + f"{'hypocentral' if self.name == 'rrup' else 'epicentral'}; name " + f"that instead if it is what you mean, rather than getting it by " + f"accident." + ) + + +#: Registered distance measures, resolved by name from configuration. +DISTANCE_MEASURES: dict[str, Any] = { + "repi": Epicentral, + "rhyp": Hypocentral, + "rrup": lambda: FiniteFaultDistance( + name="rrup", needs="a rupture surface — strike, dip, length and width" + ), + "rjb": lambda: FiniteFaultDistance( + name="rjb", + needs="the surface projection of a rupture — strike, dip, length and width", + ), +} + + +def get_distance_measure(name: str) -> DistanceMeasure: + """Resolve a registered measure by name.""" + try: + factory = DISTANCE_MEASURES[name] + except KeyError: + raise ValueError( + f"Unknown distance measure {name!r}. " + f"Available: {sorted(DISTANCE_MEASURES)}." + ) from None + measure: DistanceMeasure = factory() + return measure + + +def resolve_distance_measure( + measure: str | DistanceMeasure | None = None, +) -> DistanceMeasure: + """A measure from a name, an instance, or the configuration. + + ``None`` takes ``[windows] distance_metric``, which is the project-wide + choice and had no reader at all before this. + """ + if measure is None: + measure = str(load_config().config.windows.distance_metric) + if isinstance(measure, str): + return get_distance_measure(measure) + return measure diff --git a/src/specmod/staged.py b/src/specmod/staged.py index a25a7c7..1de6421 100644 --- a/src/specmod/staged.py +++ b/src/specmod/staged.py @@ -92,6 +92,7 @@ import numpy as np from .config import load_config +from .distance import DistanceMeasure, resolve_distance_measure from .fitting import FitSpectra if TYPE_CHECKING: # pragma: no cover @@ -134,35 +135,24 @@ class InverseDistance: modelling choice rather than a derivation, which is why this is a registry and not a hardcoded expression. - ``metric`` names a trace stat — ``rhyp`` for hypocentral, ``repi`` for - epicentral. A station missing it is a hard error rather than a silent - weight of zero: a geometry that was never set is a broken run, not a - station that should quietly stop contributing. + **Which distance is itself a choice**, and at short range not a small one: + see :mod:`specmod.distance`. ``measure=None`` takes the project-wide + setting, so a study that has decided on epicentral does not have to say so + again here. """ - metric: str = "rhyp" + #: ``None`` means "whatever the configuration says". Distance is needed by + #: geometric spreading as well as by weighting, so the choice belongs in + #: one place rather than being restated per consumer. + measure: str | DistanceMeasure | None = None name: str = "inverse_distance" def weights( self, table: Any, spectra: Any, ids: Sequence[str] ) -> NDArray[np.float64]: - out = np.empty(len(ids), dtype=np.float64) - for i, id in enumerate(ids): - meta = spectra[id].signal.meta - if self.metric not in meta: - raise ValueError( - f"{id} has no {self.metric!r}, so it cannot be weighted by " - f"distance. Set the geometry with " - f"specmod.preprocess.set_stream_distance, or choose a " - f"weighting that does not need it." - ) - distance = float(meta[self.metric]) - if distance <= 0: - raise ValueError( - f"{id} has {self.metric}={distance}, which is not a distance" - ) - out[i] = 1.0 / distance - return out + distances = resolve_distance_measure(self.measure).distances(spectra, ids) + weights: NDArray[np.float64] = 1.0 / distances + return weights @dataclass(frozen=True, slots=True) @@ -216,8 +206,13 @@ def weights( #: Registered weightings, resolved by name from ``[fitting] event_weighting``. WEIGHT_MODELS: dict[str, Any] = { - "inverse_hypocentral_distance": lambda: InverseDistance(metric="rhyp"), - "inverse_epicentral_distance": lambda: InverseDistance(metric="repi"), + # Follows the configured distance measure, so a project-wide choice is + # honoured in one place. The shipped default. + "inverse_distance": InverseDistance, + # And explicit spellings, for a study that wants to say which it used + # regardless of what the rest of the configuration says. + "inverse_hypocentral_distance": lambda: InverseDistance(measure="rhyp"), + "inverse_epicentral_distance": lambda: InverseDistance(measure="repi"), "uniform": Uniform, "inverse_variance": InverseVariance, } diff --git a/tests/test_staged.py b/tests/test_staged.py index 12ba95f..c789265 100644 --- a/tests/test_staged.py +++ b/tests/test_staged.py @@ -20,6 +20,12 @@ obspy = pytest.importorskip("obspy") from specmod.config import load_config # noqa: E402 +from specmod.distance import ( # noqa: E402 + Epicentral, + Hypocentral, + get_distance_measure, + resolve_distance_measure, +) from specmod.fitting import FitSpectra # noqa: E402 from specmod.pipeline import spectrum_set_from_streams # noqa: E402 from specmod.staged import ( # noqa: E402 @@ -60,8 +66,10 @@ def test_it_reproduces_the_weighted_mean_computed_by_hand( stage1 = FitSpectra(spectra) stage1.fit_spectra() table = stage1.table.set_index("id") + # `repi`, because that is what `[windows] distance_metric` says and the + # weighting now reads it. It used to hardcode `rhyp`. weights = np.array( - [1.0 / float(spectra[id].signal.meta["rhyp"]) for id in table.index] + [1.0 / float(spectra[id].signal.meta["repi"]) for id in table.index] ) by_hand = float((table["fc"].to_numpy() * weights).sum() / weights.sum()) @@ -91,7 +99,7 @@ def test_the_two_stages_are_both_kept(self, pnr_windows: Any) -> None: def test_it_reads_the_configured_parameter_and_weighting(self) -> None: fitting = load_config().config.fitting assert fitting.event_parameter == "fc" - assert fitting.event_weighting == "inverse_hypocentral_distance" + assert fitting.event_weighting == "inverse_distance" assert fitting.include == () assert fitting.exclude == () assert fitting.require_pass is True @@ -277,6 +285,7 @@ def test_describe_still_says_something_useful(self, pnr_windows: Any) -> None: class TestWeighting: def test_the_registry_resolves_and_rejects_by_name(self) -> None: assert isinstance(get_weight_model("uniform"), Uniform) + assert isinstance(get_weight_model("inverse_distance"), InverseDistance) assert isinstance( get_weight_model("inverse_hypocentral_distance"), InverseDistance ) @@ -334,6 +343,53 @@ class Pair: with pytest.raises(ValueError, match="set_stream_distance"): InverseDistance().weights(None, spectra, ["XX.A..HHZ"]) + def test_the_configured_distance_measure_is_what_gets_used( + self, pnr_windows: Any + ) -> None: + """`[windows] distance_metric` had no reader before this. + + It matters at short range. Hypocentral and epicentral converge far from + the source and diverge near it — on these windows the nearest station + is 0.89 km epicentral against 2.30 km hypocentral — and since this + weighting is by *inverse* distance the disagreement lands hardest on + the station carrying the most weight. + """ + spectra = _spectra(pnr_windows) + assert load_config().config.windows.distance_metric == "repi" + + configured = fit_event(spectra).value + epicentral = fit_event(spectra, weighting="inverse_epicentral_distance").value + hypocentral = fit_event(spectra, weighting="inverse_hypocentral_distance").value + + assert configured == pytest.approx(epicentral, rel=1e-12) + assert abs(hypocentral / epicentral - 1) > 0.05, ( + "the two distance measures now agree; if the geometry has changed " + "so that this no longer matters, say so here instead" + ) + + +class TestDistanceMeasures: + def test_the_registry_resolves_and_rejects(self) -> None: + assert isinstance(get_distance_measure("repi"), Epicentral) + assert isinstance(get_distance_measure("rhyp"), Hypocentral) + with pytest.raises(ValueError, match="Unknown distance measure"): + get_distance_measure("nope") + + @pytest.mark.parametrize("name", ["rrup", "rjb"]) + def test_finite_fault_measures_refuse_rather_than_degenerate( + self, name: str + ) -> None: + """For a point source these reduce exactly to hypocentral and + epicentral, so a silent fallback would give plausible numbers that are + wrong for any event large enough to justify asking for them.""" + with pytest.raises(NotImplementedError, match="rupture surface"): + get_distance_measure(name).distances(None, ["XX.A..HHZ"]) + + def test_resolve_takes_a_name_an_instance_or_the_configuration(self) -> None: + assert resolve_distance_measure("rhyp").name == "hypocentral" + assert resolve_distance_measure(Epicentral()).name == "epicentral" + assert resolve_distance_measure().name == "epicentral" # configured + class TestReporting: def test_spread_reports_what_the_mean_hides(self, pnr_windows: Any) -> None: From 5e9a3661b4ca0d5cbf77d74d3c3397dbefb146ba Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 8 Aug 2026 17:02:59 +0000 Subject: [PATCH 7/9] refactor(config): move the distance measure to its own [geometry] section MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `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. --- docs/REFACTOR_PLAN.md | 5 +++-- src/specmod/config/__init__.py | 2 ++ src/specmod/config/sections.py | 34 +++++++++++++++++++++++++++++++++- src/specmod/distance.py | 7 ++++--- studies/magna_2020_paper.toml | 7 ++++++- tests/test_staged.py | 8 ++++---- 6 files changed, 52 insertions(+), 11 deletions(-) diff --git a/docs/REFACTOR_PLAN.md b/docs/REFACTOR_PLAN.md index 40a06a3..d41d4e5 100644 --- a/docs/REFACTOR_PLAN.md +++ b/docs/REFACTOR_PLAN.md @@ -2048,8 +2048,9 @@ drop**. The choice of weighting moves it too: 11.585 (epicentral), 12.751 (hypocentral), 11.048 (uniform). **Which distance is itself configured, and `specmod.distance` is the registry -for it.** `[windows] distance_metric` existed from the start and was read by -nothing; it is read now. That is not bookkeeping at these ranges: the nearest +for it.** `[geometry] distance_measure` is read now; it existed from the start as +`[windows] distance_metric` and was read by nothing, in the wrong section — +cutting a window does not depend on how distance is measured. That is not bookkeeping at these ranges: 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 — so the two measures disagree most exactly where the inverse-distance weight is largest. Epicentral is the honest diff --git a/src/specmod/config/__init__.py b/src/specmod/config/__init__.py index a1f7b01..51705ec 100644 --- a/src/specmod/config/__init__.py +++ b/src/specmod/config/__init__.py @@ -16,6 +16,7 @@ AcquireConfig, Config, FittingConfig, + GeometryConfig, ModelConfig, SmoothingConfig, SnrConfig, @@ -29,6 +30,7 @@ "AcquireConfig", "Config", "FittingConfig", + "GeometryConfig", "ModelConfig", "Provenance", "ResolvedConfig", diff --git a/src/specmod/config/sections.py b/src/specmod/config/sections.py index 7534103..4da429c 100644 --- a/src/specmod/config/sections.py +++ b/src/specmod/config/sections.py @@ -66,7 +66,6 @@ class WindowsConfig: #: run used s=3.4; 2.9 is the shipped default and is kept as such. p_velocity: float = 5.9 s_velocity: float = 2.9 - distance_metric: Literal["repi", "rhyp"] = "repi" #: Used when an S pick is missing: s_time = p_time + emergency_ratio * (p - o). emergency_ratio: float = 1.7 @@ -271,12 +270,45 @@ class VizConfig: plot_columns: int = 3 +@dataclass(frozen=True, slots=True) +class GeometryConfig: + """Source-to-site geometry. + + Its own section because more than one stage needs it. Distance feeds the + ensemble weighting of the two-stage fit (:mod:`specmod.staged`) and the + geometric spreading a moment calculation corrects for, and a setting two + consumers each keep their own copy of is how the two come to disagree. + + It lived in ``[windows]`` until there was a second reader, which was the + wrong home even then: cutting a window does not depend on how distance is + measured. + """ + + #: Which distance, resolved through :data:`specmod.distance.DISTANCE_MEASURES`. + #: + #: ``repi`` is the default and is the honest one wherever sensor depths are + #: not known. ``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 in + #: the metadata says so. + #: + #: The choice is not a detail at short range: on the PNR data the nearest + #: station is 0.89 km epicentral against 2.30 km hypocentral, a factor of + #: 2.57, while the farthest agree to 1.00 — so anything weighted by inverse + #: distance is most sensitive to it exactly where it matters most. + #: + #: ``rrup`` and ``rjb`` are registered and raise: both need a rupture + #: surface, and for a point source they degenerate to ``rhyp`` and ``repi``. + distance_measure: str = "repi" + + @dataclass(frozen=True, slots=True) class Config: """The whole resolved configuration.""" acquire: AcquireConfig = field(default_factory=AcquireConfig) windows: WindowsConfig = field(default_factory=WindowsConfig) + geometry: GeometryConfig = field(default_factory=GeometryConfig) transform: TransformConfig = field(default_factory=TransformConfig) smoothing: SmoothingConfig = field(default_factory=SmoothingConfig) snr: SnrConfig = field(default_factory=SnrConfig) diff --git a/src/specmod/distance.py b/src/specmod/distance.py index 2bbcce4..99afd2c 100644 --- a/src/specmod/distance.py +++ b/src/specmod/distance.py @@ -179,11 +179,12 @@ def resolve_distance_measure( ) -> DistanceMeasure: """A measure from a name, an instance, or the configuration. - ``None`` takes ``[windows] distance_metric``, which is the project-wide - choice and had no reader at all before this. + ``None`` takes ``[geometry] distance_measure``, which is the project-wide + choice. It lived in ``[windows]`` and had no reader at all until this + module; cutting a window does not depend on how distance is measured. """ if measure is None: - measure = str(load_config().config.windows.distance_metric) + measure = str(load_config().config.geometry.distance_measure) if isinstance(measure, str): return get_distance_measure(measure) return measure diff --git a/studies/magna_2020_paper.toml b/studies/magna_2020_paper.toml index 754f11b..0763da7 100644 --- a/studies/magna_2020_paper.toml +++ b/studies/magna_2020_paper.toml @@ -47,7 +47,6 @@ remove_response = false # Group velocities from Pechmann et al. (2007). p_velocity = 5.9 s_velocity = 3.4 -distance_metric = "rhyp" # 20 s window opening at 80% of the elapsed Pg-Sg time, on the transverse. s_start_ratio = 0.8 s_length = 20.0 @@ -129,3 +128,9 @@ motion = "velocity" [fitting] # "We use Powell's minimization technique (Powell, 1964; Press et al., 1997)." method = "powell" + +[geometry] +# The published run used hypocentral distance. Kept as it was: this is the +# record of what was done, not a recommendation. `repi` is the shipped default +# because it does not assume sensor depths that are often unknown. +distance_measure = "rhyp" diff --git a/tests/test_staged.py b/tests/test_staged.py index c789265..c3aa673 100644 --- a/tests/test_staged.py +++ b/tests/test_staged.py @@ -66,8 +66,8 @@ def test_it_reproduces_the_weighted_mean_computed_by_hand( stage1 = FitSpectra(spectra) stage1.fit_spectra() table = stage1.table.set_index("id") - # `repi`, because that is what `[windows] distance_metric` says and the - # weighting now reads it. It used to hardcode `rhyp`. + # `repi`, because that is what `[geometry] distance_measure` says and + # the weighting now reads it. It used to hardcode `rhyp`. weights = np.array( [1.0 / float(spectra[id].signal.meta["repi"]) for id in table.index] ) @@ -346,7 +346,7 @@ class Pair: def test_the_configured_distance_measure_is_what_gets_used( self, pnr_windows: Any ) -> None: - """`[windows] distance_metric` had no reader before this. + """`[geometry] distance_measure` had no reader before this. It matters at short range. Hypocentral and epicentral converge far from the source and diverge near it — on these windows the nearest station @@ -355,7 +355,7 @@ def test_the_configured_distance_measure_is_what_gets_used( the station carrying the most weight. """ spectra = _spectra(pnr_windows) - assert load_config().config.windows.distance_metric == "repi" + assert load_config().config.geometry.distance_measure == "repi" configured = fit_event(spectra).value epicentral = fit_event(spectra, weighting="inverse_epicentral_distance").value From b1090cf47ce9c44e77c4db130d6afd108d3c617d Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 8 Aug 2026 17:07:21 +0000 Subject: [PATCH 8/9] docs(plan): record what could and could not be verified in the thesis 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 --- docs/REFACTOR_PLAN.md | 110 +++++++++++++++++++++++++++++++++++------- 1 file changed, 92 insertions(+), 18 deletions(-) diff --git a/docs/REFACTOR_PLAN.md b/docs/REFACTOR_PLAN.md index d41d4e5..f36db90 100644 --- a/docs/REFACTOR_PLAN.md +++ b/docs/REFACTOR_PLAN.md @@ -1375,14 +1375,31 @@ doctoral thesis, Chapter 1 §1.4 and Chapter 2 Eq. 2.7**, which is the Edwards et al. (2010) spectral method this package implements. Read directly rather than inferred: -| symbol | value | what it is | -|---|---|---| -| `rho` | 2600 kg/m^3 (Utah); 2800 generic | density **at the source** | -| `beta` | 3500 m/s generic | S velocity **at the source** | -| `R_0` | **1000 m** | *reference source distance* — the distance at which the source spectrum is defined | -| `F` | 2 | free surface, for **vertically incident SH** | -| `Theta-lambda-Phi` | **0.55** | average radiation pattern **of SH propagation** over the focal sphere (Boatwright, 1978) | -| partition factor | (Boore, 2003) | splits energy between vertical and horizontal ground motion | +| symbol | value | what it is | stated in | +|---|---|---|---| +| `rho` | 2600 kg/m^3 (Utah); 2800 generic | density **at the source** | Ch. 2 §2.2; Ch. 1 §1.4 | +| `beta` | 3500 m/s generic; 3520 m/s for the Utah normalisation | S velocity **at the source** | Ch. 1 §1.4; Ch. 2 §2.3 | +| `R_0` | **1000 m** | *reference source distance* — the distance at which the source spectrum is defined | Ch. 2 §2.2 | +| `F` | 2 | free surface, for **vertically incident SH** | Ch. 2 §2.2; the SH qualifier is Ch. 1 §1.4 | +| `Theta-lambda-Phi` | **0.55** | average radiation pattern **of SH propagation** over the focal sphere | Ch. 2 §2.2; Boatwright (1978) is cited in Ch. 1 §1.4 | +| partition factor | (Boore, 2003) | splits energy between vertical and horizontal ground motion | **Ch. 1 §1.4 only — not among the constants Ch. 2 enumerates for Eq. 2.7** | + +Every row above was read from the chapters directly, and the last one is the +correction. Ch. 1 describes the generic constant `C` and lists a partition +factor among its parts; Ch. 2, defining the Utah calculation this package +follows, enumerates the constants of Eq. 2.7 as `Omega_0`, `beta` at the +source, `rho`, `R_0`, `F` and the radiation pattern — and **no partition +factor**. So the implementation should not carry one by default. If a study +wants one it is a `[model]` entry it has to set, not a constant folded into +the formula. + +`beta = 3520 m/s` is worth recording alongside the generic 3500 because it is +where the Utah numbers actually come from: Ch. 2 §2.3 normalises the observed +plateaus to a theoretical Mw 3.5 source at `VS = 3520 m/s`, the value nearest +3500 in the Herrmann et al. (2011) WUS model at source depth. A cube of a +velocity is in the moment expression, so a 0.6% difference in `beta` is a 1.7% +difference in `M0` — negligible against everything else here, and worth +knowing is negligible rather than assuming so. Two things this corrects in earlier drafts of this section. @@ -1403,17 +1420,33 @@ carries the observation from there to the site. distances. `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, e.g. Holt et al. -[O] for Utah: +**kilometres** — Table 2.1's caption says so outright, "R in all cells is +hypocentral distance in kilometres", and tabulates the exponent piecewise. +Both of the thesis's own models, read from the table: - 0.88 +/- 0.02 1 < R <= 40 km - 2.93 +/- 0.28 40 < R <= 63 km + Holt et al. [O] — original + 0.88 +/- 0.02 1 < R <= 40 km + 2.93 +/- 0.28 40 < R <= 63 km 0.50 +/- 0.33 63 < R <= 100 km 1.36 +/- 0.07 100 < R <= 400 km + Holt et al. [R] — refined, and the one the thesis recommends + 0.90 +/- 0.01 1 < R <= 43 km + 2.57 +/- 0.07 43 < R <= 76 km + 0.44 +/- 0.08 76 < R <= 136 km + 1.54 +/- 0.04 136 < R <= 400 km + +The column heading is `-alpha`, so these are decay exponents: amplitude goes +as `R**-alpha`. **`[R]` is the preferred model**, and §2.7 gives the grounds — +lower uncertainty on every slope, more events resolved (218 against 201), and +the `Mw`-`Mc` relation moving closer to `Mw`-`ML`, which it should. An earlier +draft of this section quoted `[O]` alone; if a Utah spreading table is ever +shipped as study data, `[R]` is the one to ship. + Worth noting the first segment independently supports the exponent argument -above: the inverted near-field decay is **0.88**, close to the theoretical -body-wave 1, and nowhere near 2. +above, and slightly more strongly in the preferred model: the inverted +near-field decay is **0.88** and **0.90**, close to the theoretical body-wave +1, and nowhere near 2. It also shows the spreading is not a single power law but a piecewise empirical function, inverted per region. So `S(R)` should be a registered @@ -1488,10 +1521,51 @@ end to end distinguishes them. *Sourced from Holt (2019), "Addressing Uncertainty in Earthquake Magnitudes Commonly Used in Modern Seismic Hazard Assessment", University of Liverpool — -Ch. 1 §1.4 and Ch. 2 Eq. 2.7-2.8 and Table 2.1. The equations themselves are -embedded objects that do not survive text extraction, so the symbols above are -read from the surrounding prose; the equation images should be checked against -this table before any of it is implemented.* +Ch. 1 §1.4, Ch. 2 §2.2 (method and constants), §2.3 and Table 2.1 (spreading), +and §2.7 (which spreading model is preferred).* + +**Verification status, checked rather than assumed.** Chapters 1 and 2 have +since been read directly, as separate PDF exports (2.0 MB and 5.5 MB, both +well under the limit that blocked the 44 MB original). An earlier draft of +this section expected that to settle the algebraic forms. **It does not, and +the reason has changed: the equations are not in the file at all.** + +This was checked rather than inferred from a failed extraction. 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** — no math font, no vector drawing, nothing +an extractor could be failing to decode. Rendering those pages produces the +equation *numbers* against blank space. The Word-to-PDF export dropped the +embedded equation objects, and inline math went with them, 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 of the PDF will fix it** — that closes off the +route this section previously suggested. Settling it needs the original +document, or Edwards et al. (2010), which is the published method anyway. + +What the chapters *do* settle is every constant in the table above, each now +read from the prose of the section that defines it rather than inferred, plus +Table 2.1 in full — it survives as real text and matches the exponents quoted +above exactly. The one substantive correction is the partition factor, noted +under the table. + +`Mw` itself is a near miss. Ch. 1's footnote 10 confirms the thesis quotes +Hanks and Kanamori in **SI**, "the equivalent relation in SI units of Newton +meters (N·m) where 1 N·m = 1x10^7 dyne·cm" — so `M0` in N·m is the right +input and the formula is the SI one. The numeric constant is inside the +dropped equation, so **the 9.1 in `Mw = (2/3)(log10(M0) - 9.1)` is still +unverified from this source**, even though the unit convention around it now +is. + +Two claims in this section *are* independently confirmed by that prose, and +were inferences before: + +- **The spreading model's distance is in kilometres.** Eq. 2.6's description + states it outright — "is a geometrical spreading model and is distance + (km)" — rather than it being read off the units of Table 2.1. +- **The measurement is on the horizontal component of the Sg phase**, "referred + to as SgH". That is the pairing this section infers from `Theta-lambda-Phi` + being the SH average: the phase, the component and the radiation-pattern + constant are one choice, made once, in the published method. Two other things that must be true before the number means anything, both from §4.7 above: `Omega` should be the combined horizontal rather than one From 4236e4271a2f25fa02d57a1613252070cf32c2df Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 8 Aug 2026 17:55:34 +0000 Subject: [PATCH 9/9] docs(plan): source four transform defaults from the published method MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 --- docs/REFACTOR_PLAN.md | 66 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) diff --git a/docs/REFACTOR_PLAN.md b/docs/REFACTOR_PLAN.md index f36db90..a95dc23 100644 --- a/docs/REFACTOR_PLAN.md +++ b/docs/REFACTOR_PLAN.md @@ -1592,6 +1592,72 @@ velocity going unnoticed — the whole point of §4.7's units discussion. the way it now demonstrates both minimisers: showing that the choice moves the answer is what stops it being invisible. +#### The published processing recipe, and the four defaults it sources + +The equations did not survive the export, but Ch. 2 §2.2 describes the +*processing* in prose, and that survived intact. It is the most directly +useful thing the chapters give this repository, because it is a published +specification for choices §4.4 and §4.5 currently justify on internal +measurement alone. In order, for each record: + +1. Broadband stations (`HH` and `BH`) within 400 km; instrument response + deconvolved; horizontals rotated to the **transverse** component to isolate + SgH. +2. Pg and Sg velocities of 5.9 and 3.3 km/s predict the arrivals; a **20 s + window** starts at **80% of the predicted Pg-Sg interval**, the 20% margin + being deliberate slack for pick uncertainty. +3. Transform: **zero-pad to 2^N**, de-mean, **five 3-pi prolate tapers** + (Lees and Park, 1995), FFT, **multiply by the sampling period**. +4. Noise window: record start to **75% of the interval between start and the + predicted Pg**; same transform; normalised by the ratio of signal and noise + window lengths, **excluding the zero-padding**. +5. Fit only where spectral **SNR >= 3**, by Powell minimisation, with a common + source corner frequency taken from the station of smallest misfit. +6. Event `Mw` is the mean of station estimates after rejecting beyond + **2.5 sigma**, and requires **at least 3 stations**. + +Four of those land on defaults this repository already has, which is worth +recording because in each case the justification here was internal: + +- **`MultitaperEstimator` defaults to `time_bandwidth=3.0, n_tapers=5`** — + exactly "five 3-pi prolate tapers", `NW = 3` giving `2NW - 1 = 5`. §4.4 + could previously only say the 3 was "the literal 3 passed positionally to + `mtspec`, with no way to configure it". It now has a published source, and + the pair is a considered choice rather than an inherited constant. +- **"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. Same normalisation, reached from the other end. +- **The zero-padding exclusion agrees with that too.** The recipe is careful + that the signal/noise length ratio uses unpadded lengths — padding must not + move an amplitude. That is the same rule as `one_sided_fas`'s refusal to + rescale for padding, and `SnrConfig.scale_parseval`'s + `sqrt(len(signal)/len(noise))` is its amplitude-domain form. Worth checking + when the noise path is next touched that the ratio there is taken on the + cut lengths and not on `n_fft`. +- **SNR >= 3 matches `SnrConfig.tolerance = 3.0`** — though *applied + differently*, and the difference is exactly the distinction the config + already draws. The thesis gates the **fit band** per spectrum: fit where the + ratio holds. The Magna paper gates **spectrum selection**: keep spectra + above 3 in three fixed bands, which is `assert_bandwidths` with `bands`. + Same threshold, two different jobs, and the code can express both. + +One divergence, and it is a knowing one. The recipe pads to `2^N` where +`resolve_n_fft` defaults to `"fast"`, the next 5-smooth length. That is a +speed choice measured in §4.4 — `pow2` overshoots a 65537-sample record to +131072 where 65610 will do — and it is safe precisely because of the point +above: with normalisation keyed to `dt`, padding changes the frequency grid +and nothing else, so `"fast"` and `"pow2"` differ in sampling, not amplitude. +**`n_fft="pow2"` remains available and is what reproduces the published +workflow exactly**, which is the setting a study file should pin when +reproduction is the goal. + +The remaining steps — the 20 s window at 80% of the Pg-Sg interval, the noise +window at 75%, the 2.5-sigma station rejection and the 3-station minimum — +are event-level aggregation this package does not do yet, and they belong with +the `magnitude` module above rather than with the estimators. They are +recorded here so the module has a specification to build to instead of +inventing thresholds. + ### 4.8 Configuration: semantic groups, layered overrides, recorded provenance Scientific parameters are currently scattered across three places with no