feat(staged): the two-stage event fit, with selection after quality control - #24
Merged
Conversation
…ontrol `fit_event(spectra)` is the published workflow: fit every station free, take the inverse-hypocentral-distance weighted mean of the corner frequencies, fix it, refit. It reproduces the by-hand calculation exactly, which is the point — it has been rebuilt by hand by every user and by this repository's own tutorial in fifteen lines. Both stages are returned. The stage-1 spread is the evidence for how well constrained the event value is; returning only stage two would be a number with no error on it. `spread()` reports the range beside the mean, because a 2% spread and a 300% spread give the same weighted mean and mean very different things — on these windows it is 307%. Selection is the part that could not be defaulted away, and is what this adds beyond the notebook version. Quality control is a judgement, and a station that is confidently wrong moves the event value for every other station. `ChannelSelection` reads include/exclude globs from `[fitting]` and matches them at whichever level they are written — `"AQ04"` the station, `"HHE"` the component, `"UR"` the network, `"UR.AQ04.00.HHE"` the one channel. Writing the station code alone is the common case after QC, since a clipped record or a bad response is a property of the instrument rather than of one component, and requiring `"UR.AQ04.*"` for it is the kind of detail that gets mistyped as `"AQ04"` and silently matches nothing. Every exclusion carries a reason naming the level it matched at. "Excluded" is not actionable; "matched exclude='AQ04' at station" is. Weighting is a registry — inverse hypocentral (the published choice), inverse epicentral, uniform, inverse variance — because it is a modelling choice. That it matters is measured: the nearest two channels carry 22.8% of the weight and the nearest four 36.1%, so dropping the single nearest station moves the event corner 12.751 -> 14.774 Hz, 16%, which is 1.56x in stress drop. Weighting choice moves it comparably: 12.751 / 11.585 / 11.048. `inverse_variance` raises rather than falling back to uniform when the uncertainties are missing. Powell estimates no covariance matrix, so under the shipped default they always are, and silently becoming uniform is the kind of substitution that ends up in a paper. Stage two is skipped, with `stage2=None`, when selection leaves nothing. Fixing the corner to a mean of no stations would invent the number the second stage exists to constrain. One trap, found by testing the module against its own docstring and now pinned. `require_pass` drops a station whose stage-1 fit ended against a bound, and `pass_fitting` asks whether `value +/- stderr` reaches one. Powell reports no `stderr`, so the spread is zero and the test almost never fires: it drops 0 stations under Powell and 6 under `leastsq`. Changing the minimiser therefore changes *which stations vote*, not only how each is fitted — the naive comparison gives 144% where the like-for-like one gives 0.6%. The docstring's numbers now carry that condition explicitly, and two tests hold both halves of it.
…by side
The notebook keeps the hand-written two-stage fit and adds `fit_event`
underneath it. The long version stays because the two stages and the weighted
mean between them are the method rather than an implementation detail — a
reader who has not seen them cannot judge a result that came out of them — and
the short version is there because nobody should retype it. The notebook
asserts the two agree to 1e-12 rather than claiming it, so the demonstration
doubles as a check that the API does what the worked example does.
Then the part the notebook could not show before, because the capability did
not exist: which channels vote. Inverse-distance weighting turns out to be
concentrated — the nearest two channels carry 22.8% of the weight and the
nearest eight carry 56.3% — so dropping the single nearest station moves the
event corner 12.751 -> 14.774 Hz, which is 1.56x in stress drop from one
quality-control decision. Shown with a bare station code, since that is the
form a reader will reach for.
And the `require_pass` trap, demonstrated rather than described: 28 channels
vote under Powell and 22 under `leastsq`, so the two look 144% apart until the
ensemble is held fixed, at which point they agree to 0.6%.
Fixed while writing it: the new section bound `station`, which a cell twenty
lines later reads as a full trace id when checking the HDF5 round trip. The
notebook is executed end to end in CI-adjacent form here, which is the only
reason a name collision three sections apart was caught at all.
Also records three roadmap items raised while reviewing this, none of them
built:
- **Documentation equations do not render.** `processing.md` uses `$...$` and
`$$...$$`, which is MyST `dollarmath` syntax, and neither the extension nor
the Sphinx build exists — so the only renderer these files have ever met is
GitHub's. A scan also finds one display block without a preceding blank line
and one spanning multiple lines, which break under MyST too. Prose that has
never been rendered has never been checked.
- **Components and phases (new §4.7).** `back_azimuth` is computed on every
trace and read by nothing, which is the value component rotation needs.
Horizontal channels are fitted as independent measurements — 28 "channels"
for 14 stations — so each station is counted twice in the ensemble the
previous commit added, and between-component disagreement enters as if it
were between-station scatter. Records both horizontal treatments the
published work uses, notes that the geometric mean is not rotation-invariant
before someone assumes it is, and sets out what P-wave parameters need that
S-wave ones do not: `(alpha/beta)^3` is about 5.2, so the phase velocity
alone is a factor of five in seismic moment. The constants are written down
as the shape of the problem and explicitly not as values to adopt — the
Magna work already pins `F = 2` and `Theta-lambda-Phi = 0.55`, which is not
the textbook number.
- **Station and channel identity should be a type.** Three independent
`split(".")` spellings already exist, written in the last week. More to the
point, "components counted as independent stations" is invisible because a
channel and a station are both `str` — summing over one when you meant the
other type-checks perfectly. ObsPy builds the id string and filters on
components but has no value type, which is why every project rewrites this.
Made a prerequisite for §4.7 rather than a parallel nicety: "the horizontals
of each station" has to be expressible before it can be implemented.
Configuration moves to §4.8; all six cross-references repointed.
…what checking it found
Rewrites §4.7's identity subsection against the actual conventions rather than
against the dotted string. The hierarchy is Network / Station / Location /
Channel and each level carries something the code currently ignores or
re-derives: the network code is reused across epochs, the station is the site
and holds the coordinates, the location code distinguishes co-located sensor
packages, and the channel holds azimuth, dip, sample rate and the response —
per epoch.
Four consequences worth having written down before any of this is built.
The channel code is three fields, not an opaque string. Band, instrument and
orientation. `staged.ChannelSelection` matches `"HHE"` whole, so it cannot
express "every high-gain seismometer regardless of band" or "every
accelerometer", both ordinary requests.
`N`/`E` against `1`/`2` is a correctness trap. FDSN uses `1`/`2` precisely
when the horizontals are *not* aligned to north and east — so globbing for
`HHE`/`HHN`, which this repository's own tutorial does when reading waveforms,
finds nothing on such a network, and rotating by assuming `N` is 0 degrees is
wrong there. Rotation needs the channel's own azimuth and dip.
Response and coordinates are epoch-scoped. An instrument swap makes a new
channel epoch, and the right response depends on the time of the record.
Flattening to "the response" silently applies the wrong one across a swap.
`azimuth` already means two things. `preprocess` writes `tr.stats["azimuth"]`
as source-receiver geometry; `Channel.azimuth` is component orientation.
Unrelated quantities, one letter apart, and rotation needs both at once.
Checking the committed inventory to write this turned up three things, all
recorded: every channel has `azimuth=None` and `dip=None`, so rotating the
PNR data means assuming the `N`/`E` orientation rather than reading it;
`depth` is `123456.0`, a placeholder rather than a measurement; and the file
mixes location codes `''` and `'00'`, which under SEED are different
locations rather than synonyms.
Also states what ObsPy does and does not give, since that was the open
question. The metadata side is solved — `Inventory`/`Station`/`Channel` model
the hierarchy and `get_channel_metadata` resolves an epoch — and should be
used rather than reimplemented. What is missing is the identifier: `Trace.id`
builds a string, `Stream.select` filters it, and there is no value type, no
station-level equality, no grouping, and no access to the band/instrument/
orientation split. Hence three `split(".")` implementations in this codebase
alone.
…ficient Corrects a framing in §4.7. It noted that the Magna work uses `Theta-lambda-Phi = 0.55` against the textbook 0.63 and left that hanging as though the discrepancy were explained. It is not, and a unit conversion cannot explain it: the radiation pattern term is dimensionless, and changing units cannot move a dimensionless number. 0.63 is the RMS average over the focal sphere and other averages are also in the literature, so which average was taken is the likelier answer — a question for the paper rather than for this document to settle. Units are a real hazard in the same formula, though, and much larger than the 0.87 that coefficient would account for. Density carries 10^3 between kg/m^3 and g/cm^3, distance another 10^3, the moment itself 10^7 between newton-metres and dyne-centimetres, and velocity — which is **cubed** — carries 10^9 between m/s and km/s. `r_0 = 1000 m` sitting among the recorded parameters is exactly the shape of a constant bridging kilometres and metres, so it is worth checking rather than assuming in either direction. Records the arithmetic check that settles it without reading any code: `M0 = 10 ** (1.5 * Mw + 9.1)` newton-metres, so the Mw 1.6 event this repository uses throughout should come out near 3e11 N m, or 3e18 in dyne-centimetres. A factor of 10^9 from a km/s velocity is not subtle. SpecMod computes neither M0 nor Mw today — the constants appear only in this plan, describing analysis code that was never part of the package — so nothing is wrong in the repository. The point is what to build: the ambiguity exists because the constants live as bare floats in prose with no units attached, which is this document's recurring failure applied to the number the package exists to produce. §4.2 gave spectra `Motion` and `AmplitudeKind` so a conversion could not be applied twice or backwards. The moment calculation needs the same: units on the inputs, a declared output unit, and a test asserting the Mw of a known event, so the question answers itself rather than being reasoned about.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
fit_event(spectra)is the published workflow: fit every station free, takethe inverse-hypocentral-distance weighted mean of the corner frequencies, fix
it, refit. It reproduces the by-hand calculation exactly, which is the point —
it has been rebuilt by hand by every user and by this repository's own
tutorial in fifteen lines.
Both stages are returned. The stage-1 spread is the evidence for how well
constrained the event value is; returning only stage two would be a number
with no error on it.
spread()reports the range beside the mean, because a2% spread and a 300% spread give the same weighted mean and mean very
different things — on these windows it is 307%.
Selection is the part that could not be defaulted away, and is what this adds
beyond the notebook version. Quality control is a judgement, and a station
that is confidently wrong moves the event value for every other station.
ChannelSelectionreads include/exclude globs from[fitting]and matchesthem at whichever level they are written —
"AQ04"the station,"HHE"thecomponent,
"UR"the network,"UR.AQ04.00.HHE"the one channel. Writing thestation code alone is the common case after QC, since a clipped record or a
bad response is a property of the instrument rather than of one component, and
requiring
"UR.AQ04.*"for it is the kind of detail that gets mistyped as"AQ04"and silently matches nothing.Every exclusion carries a reason naming the level it matched at. "Excluded" is
not actionable; "matched exclude='AQ04' at station" is.
Weighting is a registry — inverse hypocentral (the published choice), inverse
epicentral, uniform, inverse variance — because it is a modelling choice. That
it matters is measured: the nearest two channels carry 22.8% of the weight and
the nearest four 36.1%, so dropping the single nearest station moves the event
corner 12.751 -> 14.774 Hz, 16%, which is 1.56x in stress drop. Weighting
choice moves it comparably: 12.751 / 11.585 / 11.048.
inverse_varianceraises rather than falling back to uniform when theuncertainties are missing. Powell estimates no covariance matrix, so under the
shipped default they always are, and silently becoming uniform is the kind of
substitution that ends up in a paper.
Stage two is skipped, with
stage2=None, when selection leaves nothing.Fixing the corner to a mean of no stations would invent the number the second
stage exists to constrain.
One trap, found by testing the module against its own docstring and now
pinned.
require_passdrops a station whose stage-1 fit ended against abound, and
pass_fittingasks whethervalue +/- stderrreaches one. Powellreports no
stderr, so the spread is zero and the test almost never fires: itdrops 0 stations under Powell and 6 under
leastsq. Changing the minimisertherefore changes which stations vote, not only how each is fitted — the
naive comparison gives 144% where the like-for-like one gives 0.6%. The
docstring's numbers now carry that condition explicitly, and two tests hold
both halves of it.