From 0702ba719928f06108c5c5066df6e25aaeade091 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 8 Aug 2026 15:59:26 +0000 Subject: [PATCH 1/4] feat(staged): the two-stage event fit, with selection after quality control MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `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. --- docs/REFACTOR_PLAN.md | 39 ++- src/specmod/config/sections.py | 32 ++ src/specmod/staged.py | 528 +++++++++++++++++++++++++++++++++ tests/test_staged.py | 351 ++++++++++++++++++++++ 4 files changed, 949 insertions(+), 1 deletion(-) create mode 100644 src/specmod/staged.py create mode 100644 tests/test_staged.py diff --git a/docs/REFACTOR_PLAN.md b/docs/REFACTOR_PLAN.md index d782c02..ebf491a 100644 --- a/docs/REFACTOR_PLAN.md +++ b/docs/REFACTOR_PLAN.md @@ -1604,7 +1604,7 @@ SpecMod computes. The two-stage fit is, separately, a good candidate for the new API: it is the workflow the science actually uses, and it currently has to be rebuilt by hand -by every user. +by every user. **Built** — `specmod.staged.fit_event`; see below. The tutorial now does exactly that rebuilding-by-hand, deliberately, because it is the clearest available demonstration of *why* the second stage exists. On @@ -1624,6 +1624,43 @@ second stage rests on. Fifteen lines of notebook is the right amount of code for that to cost, and it is the argument for putting those fifteen lines behind an API rather than leaving each user to write them again. +#### 5.2.5b The two-stage fit, as built + +`specmod.staged.fit_event(spectra)` is the published workflow with every +argument defaulted from `[fitting]`. It reproduces the by-hand weighted mean +exactly, and both stages are returned rather than only the second — the +stage-1 spread is the evidence for how well constrained the event value is, +and returning one number without it would be reporting a measurement with no +error on it. + +**Selection is the part that could not be defaulted away.** Quality control is +a judgement, and a station that is confidently wrong — bad response, clipped +record, pick on the wrong phase — moves the event value for every other +station. `ChannelSelection` reads `include`/`exclude` globs from the study +file and matches them at whichever level they are written: `"AQ04"` is the +station, `"HHE"` the component, `"UR"` the network, `"UR.AQ04.00.HHE"` the one +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. + +**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 +`value +/- stderr` reaches one — and Powell, the shipped minimiser, estimates +no covariance matrix, 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 just how each is +fitted, and the naive cross-minimiser comparison gives 144% where the +like-for-like one gives 0.6%. The flag is doing the right thing when it fires; +the asymmetry is what needed writing down. + #### 5.2.6 The three-way comparison The published values were produced by code containing the bugs in §2.2 and §2.5, diff --git a/src/specmod/config/sections.py b/src/specmod/config/sections.py index f41c329..db79519 100644 --- a/src/specmod/config/sections.py +++ b/src/specmod/config/sections.py @@ -227,6 +227,38 @@ class FittingConfig: t_star_min: float = 1e-4 corner_frequency_min: float = 0.0 + #: The two-stage event fit; see :mod:`specmod.staged`. + #: + #: One spectrum cannot separate the source corner from the path + #: attenuation — they trade off on the falling limb — so the corner is + #: determined by the ensemble and then held fixed while each station + #: refits the rest. ``event_parameter`` is what the ensemble decides. + #: ``"fc"`` because that is the term belonging to the source; ``"ts"`` is + #: the meaningful alternative for a study with an independent handle on Q. + event_parameter: str = "fc" + #: 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" + + #: Which channels contribute to the event value, as shell globs matched + #: against the trace id and each of its SEED components — so ``"AQ07"`` + #: means the station, ``"HHE"`` means the component, ``"UR"`` means the + #: network. Empty ``include`` means "everything not excluded". + #: + #: These exist to be edited *after* looking at a first pass. Quality + #: control is a judgement — a clipped record, a bad response, a pick on + #: the wrong phase — and a station that is confidently wrong moves the + #: event value for every other station. Putting the decision in the study + #: file is what makes it part of the record rather than something done in + #: a notebook and forgotten. + include: tuple[str, ...] = () + exclude: tuple[str, ...] = () + #: Drop a station whose stage-1 fit ended with a parameter pinned against + #: one of its bounds. The value reported there is the bound rather than a + #: measurement, so averaging it in is averaging in a constant. + require_pass: bool = True + @dataclass(frozen=True, slots=True) class VizConfig: diff --git a/src/specmod/staged.py b/src/specmod/staged.py new file mode 100644 index 0000000..a25a7c7 --- /dev/null +++ b/src/specmod/staged.py @@ -0,0 +1,528 @@ +"""The two-stage event fit, and the channel selection that feeds it. + +Fitting a source model to one spectrum is not a unique inversion. The source +corner and the path attenuation trade off against each other on the falling +limb, and two minimisers can reach the same reduced chi-squared at corner +frequencies differing by tens of percent — a factor of several in stress drop, +which scales as ``fc**3``. On the 28 PNR windows, Powell and ``leastsq`` land +at 21.26 Hz and 14.75 Hz on one station at redchi 0.0259 against 0.0254. + +That is not resolvable from one spectrum by any minimiser, and the published +workflow does not try. ``f_c`` belongs to the **source**: every station sees +the same rupture, so there is one value of it for the event. ``t*`` belongs to +the **path**, and every station has a different one. So a station whose ``t*`` +came out too high returns a corner that is too high, and the next station's +error does not point the same way. Averaging over the ensemble is not +cosmetic smoothing — it uses the fact that the quantity being averaged is +common to all of them while the contaminating one is not. + +The two stages +-------------- +1. ``Omega``, ``f_c`` and ``t*`` free at every station independently. The + output is not the answer; it is N noisy estimates of one number plus N + estimates of N different numbers. +2. The event ``f_c`` — a weighted mean over the stations that survived + selection — is held fixed, and every station refits ``Omega`` and ``t*`` + against a corner it can no longer trade against. + +Measured on the same 28 windows, with the same 28 channels contributing to +both: the two minimisers differ by a factor 1.44 in ``f_c`` at the worst +station and by **0.6%** on the event value. After stage two they agree to +0.31% on ``t*`` and 2.3e-3 log10 units on ``Omega``. + +"With the same channels contributing" is load-bearing, and is the trap in this +module. See :attr:`ChannelSelection.require_pass` — comparing two minimisers +under the default selection compares two different ensembles, and gives 144% +rather than 0.6%. + +Be clear about what that last number is not. Fixing ``f_c`` removes the +parameter the minimisers were disagreeing about, so of course they then agree. +What it shows is that the residual two-parameter problem is well conditioned: +once the corner is pinned, ``Omega`` and ``t*`` are determined by the spectrum +rather than negotiable. The judgement is concentrated into one number for the +whole event, and that number came from the ensemble. + +Choosing which channels contribute +---------------------------------- +Selection is the part that cannot be automated away, because it is where +quality control enters. A station with a bad instrument response, a clipped +record or a pick on the wrong phase produces a corner frequency that is +confidently wrong, and averaging it in moves the event value for every other +station. + +So the ensemble is chosen by :class:`ChannelSelection`, which reads +``[fitting]`` and can be overridden per call. The order is fixed and each step +is recorded with a reason, so ``StagedFit.excluded`` says why any given +channel is not contributing: + +1. anything not matching ``include`` (when ``include`` is non-empty), +2. anything matching ``exclude``, +3. anything whose stage-1 fit failed its bounds, when ``require_pass``, +4. anything left with no fit at all. + +Patterns are shell globs, and they match at whichever level you write them. +A trace id is ``NET.STA.LOC.CHA``, and a pattern is tried against the whole +id, against ``NET.STA``, and against each component on its own — so all of +these do what they look like they do:: + + "AQ07" every channel of that station + "UR" every station of that network + "UR.AQ07" that station, spelled unambiguously + "HHE" every east component, at every station + "HH?" every high-gain broadband channel + "UR.AQ07.00.HHE" exactly that channel + "LV.L00[123]..HH?" a glob over the full id + +Writing the station code alone is the common case after quality control — +a clipped record or a bad response is a property of the instrument, not of one +component — and needing ``"UR.AQ07.*"`` for it is the kind of detail that gets +mistyped as ``"AQ07"`` and silently matches nothing. Station and channel codes +do not collide in practice; where a pattern could match at two levels it +matches, and :attr:`StagedFit.excluded` records which level it hit. +""" + +from __future__ import annotations + +import contextlib +import fnmatch +import io +from dataclasses import dataclass, field +from typing import TYPE_CHECKING, Any, Protocol, runtime_checkable + +import numpy as np + +from .config import load_config +from .fitting import FitSpectra + +if TYPE_CHECKING: # pragma: no cover + from collections.abc import Iterable, Mapping, Sequence + + from numpy.typing import NDArray + +__all__ = [ + "WEIGHT_MODELS", + "ChannelSelection", + "StagedFit", + "WeightModel", + "fit_event", + "get_weight_model", +] + + +@runtime_checkable +class WeightModel(Protocol): + """How much each station's stage-1 estimate counts toward the event value. + + Given the stage-1 table and the spectra it was fitted to, return one weight + per row. Weights need not sum to anything; they are normalised on use. + """ + + name: str + + def weights( + self, table: Any, spectra: Any, ids: Sequence[str] + ) -> NDArray[np.float64]: ... # pragma: no cover + + +@dataclass(frozen=True, slots=True) +class InverseDistance: + """Weight by ``1 / distance``, the published choice. + + The nearer station has less path between the source and the sensor, so + less of its high-frequency falloff can be attenuation and its corner is the + better constrained of the two. That is the argument for it; it is a + 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. + """ + + metric: str = "rhyp" + 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 + + +@dataclass(frozen=True, slots=True) +class Uniform: + """Every contributing station counts the same. + + The honest default when the geometry is unknown, and the right one when + the stations are at comparable distances — where inverse-distance weighting + adds variance without adding information. + """ + + name: str = "uniform" + + def weights( + self, table: Any, spectra: Any, ids: Sequence[str] + ) -> NDArray[np.float64]: + return np.ones(len(ids), dtype=np.float64) + + +@dataclass(frozen=True, slots=True) +class InverseVariance: + """Weight by ``1 / stderr**2`` on the aggregated parameter. + + Statistically the right answer when the uncertainties are real, and + unavailable under the shipped minimiser: Powell estimates no covariance + matrix, so lmfit reports ``stderr`` as ``None`` for every parameter. This + raises rather than falling back, because silently becoming uniform is the + kind of substitution that ends up in a paper. + """ + + name: str = "inverse_variance" + + def weights( + self, table: Any, spectra: Any, ids: Sequence[str] + ) -> NDArray[np.float64]: + indexed = table.set_index("id") + column = f"{_aggregated_parameter()}-stderr" + if column not in indexed: + raise ValueError(f"the fit table has no {column!r} column") + errors = indexed.loc[list(ids), column].to_numpy(dtype=np.float64) + if not np.isfinite(errors).all(): + raise ValueError( + "inverse-variance weighting needs an uncertainty on every " + "contributing station, and some are missing. The configured " + "minimiser is likely 'powell', which estimates no covariance " + "matrix — use 'leastsq', or weight another way." + ) + weights: NDArray[np.float64] = 1.0 / errors**2 + return 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"), + "uniform": Uniform, + "inverse_variance": InverseVariance, +} + + +def get_weight_model(name: str) -> WeightModel: + """Resolve a registered weighting by name, with its defaults.""" + try: + factory = WEIGHT_MODELS[name] + except KeyError: + raise ValueError( + f"Unknown weighting {name!r}. Available: {sorted(WEIGHT_MODELS)}." + ) from None + model: WeightModel = factory() + return model + + +def _aggregated_parameter() -> str: + return str(load_config().config.fitting.event_parameter) + + +@dataclass(frozen=True, slots=True) +class ChannelSelection: + """Which channels contribute to the event value. + + Defaults come from ``[fitting]``; pass one of these to override per call. + The point of it being a value rather than four arguments is that a + selection can be written down, compared and stored — a run that dropped + three stations should be able to say so afterwards. + """ + + include: tuple[str, ...] = () + exclude: tuple[str, ...] = () + #: Drop a station whose stage-1 fit ended with a parameter against a bound. + #: + #: **This interacts with the minimiser, and not symmetrically.** + #: ``pass_fitting`` asks whether ``value +/- stderr`` reaches a bound. + #: Powell — the shipped default — estimates no covariance matrix, so + #: ``stderr`` is ``None``, the spread is zero, and the test degenerates to + #: "is the value exactly on the bound", which essentially never fires. On + #: these 28 windows it drops **0** stations under Powell and **6** under + #: ``leastsq``. + #: + #: So changing the minimiser changes *which stations vote*, not just how + #: each one is fitted. Comparing two minimisers under the default + #: selection therefore compares two different ensembles: 12.751 Hz from 28 + #: channels against 5.228 Hz from 22, a 144% difference that says almost + #: nothing about the minimisers. Holding the ensemble fixed with + #: ``require_pass=False`` gives 12.751 against 12.675 — 0.6%, which is the + #: honest comparison and the number quoted in the module docstring. + #: + #: Left ``True`` because the flag is doing the right thing when it fires: + #: a pinned parameter reports its bound rather than a measurement, and + #: averaging that in averages in a constant. But a study comparing + #: minimisers, or reporting a corner frequency alongside one fitted + #: another way, has to set this ``False`` or account for it. + require_pass: bool = True + + @classmethod + def from_config(cls) -> ChannelSelection: + fitting = load_config().config.fitting + return cls( + include=tuple(fitting.include), + exclude=tuple(fitting.exclude), + require_pass=bool(fitting.require_pass), + ) + + def choose(self, fit: FitSpectra) -> tuple[list[str], dict[str, str]]: + """Split ``fit``'s stations into contributors and reasoned exclusions.""" + contributing: list[str] = [] + excluded: dict[str, str] = {} + + table = fit.table.set_index("id") if len(fit.table) else None + # `ids()` on a SpectrumSet, plain iteration on anything else — the same + # duck-typed contract the fitter itself accepts. + spectra = fit.spectra + ids = spectra.ids() if hasattr(spectra, "ids") else list(spectra) + for id in ids: + if id not in fit.models: + excluded[id] = "no fit: the station did not pass the band gate" + continue + if self.include and _match(id, self.include) is None: + excluded[id] = f"not matched by include={list(self.include)}" + continue + hit = _match(id, self.exclude) + if hit is not None: + level, pattern = hit + excluded[id] = f"matched exclude={pattern!r} at {level}" + continue + if table is None or id not in table.index: + excluded[id] = "no row in the stage-1 table: the fit did not run" + continue + if self.require_pass and not bool(table.loc[id, "pass_fitting"]): + excluded[id] = ( + "stage-1 fit failed: a parameter is pinned against a bound, " + "so the value reported is the bound rather than a measurement" + ) + continue + contributing.append(id) + + return contributing, excluded + + +#: The parts of ``NET.STA.LOC.CHA`` a pattern may be written against, in the +#: order a match is reported. Full id first so the most specific spelling wins +#: the explanation. +def _levels(id: str) -> list[tuple[str, str]]: + parts = id.split(".") + if len(parts) != 4: + # Not a SEED id — match it whole and say nothing more about it. + return [("id", id)] + net, sta, loc, cha = parts + return [ + ("id", id), + ("station", f"{net}.{sta}"), + ("network", net), + ("station", sta), + ("location", loc), + ("channel", cha), + ] + + +def _match(id: str, patterns: Iterable[str]) -> tuple[str, str] | None: + """The (level, pattern) a channel matched at, or ``None``. + + Returned rather than a bool so an exclusion can say *why*: "matched + exclude='AQ07' at station" is actionable where "excluded" is not. + """ + for pattern in patterns: + for level, value in _levels(id): + if fnmatch.fnmatch(value, pattern): + return level, pattern + return None + + +@dataclass(frozen=True, slots=True) +class StagedFit: + """The result of :func:`fit_event`: both stages and how they were joined.""" + + #: Every station fitted independently. Kept because the spread across it is + #: the evidence for how well constrained the event value is, and discarding + #: it would leave only a number with no error on it. + stage1: FitSpectra + #: Every station refitted with the event parameter held fixed. ``None`` when + #: no station survived selection, in which case there is nothing to fix it + #: to and the honest result is stage one alone. + stage2: FitSpectra | None + #: The aggregated value and what produced it. + parameter: str + value: float + weighting: str + contributing: tuple[str, ...] + excluded: Mapping[str, str] = field(default_factory=dict) + + @property + def table(self) -> Any: + """The stage-2 table where there is one, else stage one's.""" + return self.stage1.table if self.stage2 is None else self.stage2.table + + def spread(self) -> dict[str, float]: + """How much the contributing stations disagreed, before aggregation. + + The number worth reporting beside the event value. A 2% spread and a + 60% spread give the same weighted mean and mean very different things. + """ + if not self.contributing: + return {} + values = ( + self.stage1.table.set_index("id") + .loc[list(self.contributing), self.parameter] + .to_numpy(dtype=np.float64) + ) + return { + "n": float(values.size), + "min": float(values.min()), + "max": float(values.max()), + "median": float(np.median(values)), + "weighted_mean": self.value, + "relative_spread": float((values.max() - values.min()) / self.value) + if self.value + else float("nan"), + } + + def describe(self) -> str: + """One paragraph a caller can print or paste into a log.""" + lines = [ + f"{self.parameter} = {self.value:.4g} " + f"from {len(self.contributing)} channels, weighted by " + f"{self.weighting}", + ] + spread = self.spread() + if spread: + lines.append( + f" stage-1 range {spread['min']:.4g} to {spread['max']:.4g} " + f"({100 * spread['relative_spread']:.0f}% of the event value)" + ) + if self.excluded: + lines.append(f" {len(self.excluded)} excluded:") + lines += [f" {id}: {why}" for id, why in sorted(self.excluded.items())] + return "\n".join(lines) + + +def fit_event( + spectra: Any, + *, + model: Any = None, + guess: Mapping[str, Mapping[str, float]] | None = None, + fit_bins: bool | None = None, + parameter: str | None = None, + weighting: str | WeightModel | None = None, + selection: ChannelSelection | None = None, + quiet: bool = True, + **fit_kwargs: Any, +) -> StagedFit: + """Fit an event in two stages, with the ensemble deciding the source term. + + Every argument has a configured default, so ``fit_event(spectra)`` is the + published workflow and the rest is there for when the default is wrong. + + Parameters + ---------- + spectra + A :class:`~specmod.core.SpectrumSet`, or anything mapping trace ids to + paired spectra. + model, guess, fit_bins + Passed to :class:`~specmod.fitting.FitSpectra` for both stages, so the + two are fitting the same thing. + parameter + Which parameter the ensemble determines and stage two holds fixed. + ``[fitting] event_parameter``, default ``"fc"`` — the corner frequency + is the term that belongs to the source. ``"ts"`` is the meaningful + alternative, for a study with an independent handle on ``Q``. + weighting + A registered name or a :class:`WeightModel`. ``[fitting] + event_weighting``, default inverse hypocentral distance. + selection + Which channels contribute. Defaults to :meth:`ChannelSelection.from_config`. + quiet + Suppress the fitter's per-station chatter, which is two full passes of + it here. The failures are on :attr:`StagedFit.excluded` either way. + **fit_kwargs + Passed to :meth:`~specmod.fitting.FitSpectra.fit_spectra` for both + stages — ``method=`` most usefully. + + Notes + ----- + Stage two is skipped, with ``stage2=None``, when selection leaves nothing. + Fixing the parameter to a mean of no stations would be inventing the very + number the second stage exists to constrain, and a caller who asked for two + stages and got one should be able to see that rather than read a value. + """ + fitting = load_config().config.fitting + parameter = parameter or str(fitting.event_parameter) + if selection is None: + selection = ChannelSelection.from_config() + if weighting is None: + weighting = str(fitting.event_weighting) + model_ = get_weight_model(weighting) if isinstance(weighting, str) else weighting + + def run(constant: tuple[str, float] | None = None) -> FitSpectra: + with contextlib.redirect_stdout(io.StringIO() if quiet else None): + fit = FitSpectra(spectra, model=model, guess=guess, fit_bins=fit_bins) + if constant is not None: + fit.set_const(*constant) + fit.fit_spectra(**fit_kwargs) + return fit + + stage1 = run() + contributing, excluded = selection.choose(stage1) + + if not contributing: + return StagedFit( + stage1=stage1, + stage2=None, + parameter=parameter, + value=float("nan"), + weighting=getattr(model_, "name", str(weighting)), + contributing=(), + excluded=excluded, + ) + + table = stage1.table.set_index("id") + values = table.loc[contributing, parameter].to_numpy(dtype=np.float64) + weights = np.asarray( + model_.weights(stage1.table, spectra, contributing), dtype=np.float64 + ) + if weights.shape != values.shape: + raise ValueError( + f"{getattr(model_, 'name', weighting)} returned {weights.size} " + f"weights for {values.size} channels" + ) + total = float(weights.sum()) + if total <= 0: + raise ValueError( + f"{getattr(model_, 'name', weighting)} gave every contributing " + "channel zero weight, so there is no event value to compute" + ) + value = float((values * weights).sum() / total) + + return StagedFit( + stage1=stage1, + stage2=run(constant=(parameter, value)), + parameter=parameter, + value=value, + weighting=getattr(model_, "name", str(weighting)), + contributing=tuple(contributing), + excluded=excluded, + ) diff --git a/tests/test_staged.py b/tests/test_staged.py new file mode 100644 index 0000000..12ba95f --- /dev/null +++ b/tests/test_staged.py @@ -0,0 +1,351 @@ +"""The two-stage event fit, and the selection that decides who votes. + +The science this encodes is in :mod:`specmod.staged`. What these tests hold is +that the default path is the published workflow, that every departure from it +is deliberate and recorded, and that the aggregate cannot be computed from +nothing. +""" + +from __future__ import annotations + +import contextlib +import functools +import io +from dataclasses import dataclass +from typing import Any + +import numpy as np +import pytest + +obspy = pytest.importorskip("obspy") + +from specmod.config import load_config # noqa: E402 +from specmod.fitting import FitSpectra # noqa: E402 +from specmod.pipeline import spectrum_set_from_streams # noqa: E402 +from specmod.staged import ( # noqa: E402 + WEIGHT_MODELS, + ChannelSelection, + InverseDistance, + Uniform, + fit_event, + get_weight_model, +) + + +@functools.cache +def _spectra(windows: Any) -> Any: + signal, noise = windows() + with contextlib.redirect_stdout(io.StringIO()): + return spectrum_set_from_streams(signal, noise) + + +@functools.cache +def _event(windows: Any) -> Any: + return fit_event(_spectra(windows)) + + +class TestTheDefaultPath: + def test_it_reproduces_the_weighted_mean_computed_by_hand( + self, pnr_windows: Any + ) -> None: + """The whole point of the module is that nobody writes this again. + + So it has to give what writing it by hand gives — the inverse + hypocentral distance weighted mean of the stage-1 corner frequencies. + """ + spectra = _spectra(pnr_windows) + staged = _event(pnr_windows) + + with contextlib.redirect_stdout(io.StringIO()): + stage1 = FitSpectra(spectra) + stage1.fit_spectra() + table = stage1.table.set_index("id") + weights = np.array( + [1.0 / float(spectra[id].signal.meta["rhyp"]) for id in table.index] + ) + by_hand = float((table["fc"].to_numpy() * weights).sum() / weights.sum()) + + assert staged.value == pytest.approx(by_hand, rel=1e-9) + + def test_stage_two_holds_the_event_value_fixed_everywhere( + self, pnr_windows: Any + ) -> None: + """Fixed, not merely seeded. A parameter that is still free has not + been constrained by the ensemble, it has only been given a hint.""" + staged = _event(pnr_windows) + assert staged.stage2 is not None + for id, model in staged.stage2.models.items(): + assert model.params["fc"].vary is False, id + assert model.params["fc"].value == pytest.approx(staged.value), id + + def test_the_two_stages_are_both_kept(self, pnr_windows: Any) -> None: + """The stage-1 spread is the evidence for how well constrained the + event value is. Returning only stage two would leave a number with no + error on it.""" + staged = _event(pnr_windows) + assert len(staged.stage1.models) == 28 + assert staged.stage2 is not None + assert len(staged.stage2.models) == 28 + assert staged.table is staged.stage2.table + + 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.include == () + assert fitting.exclude == () + assert fitting.require_pass is True + + +class TestSelection: + """Patterns match at whichever level they are written.""" + + @pytest.mark.parametrize( + ("pattern", "expected"), + [ + ("AQ04", 2), # a station, by its bare code + ("UR.AQ04", 2), # the same station, spelled out + ("UR.AQ04.00.HHE", 1), # one channel + ("HHE", 14), # one component, everywhere + ("UR", 16), # a whole network (LV has the other 12) + ("HH?", 28), # a glob over the channel code + ("AQ0[45]", 4), # a glob over the station code + ], + ) + def test_a_pattern_excludes_what_it_names( + self, pattern: str, expected: int, pnr_windows: Any + ) -> None: + spectra = _spectra(pnr_windows) + staged = fit_event(spectra, selection=ChannelSelection(exclude=(pattern,))) + dropped = [ + id for id, why in staged.excluded.items() if "matched exclude" in why + ] + assert len(dropped) == expected, sorted(dropped) + + def test_include_keeps_only_what_it_names(self, pnr_windows: Any) -> None: + staged = fit_event( + _spectra(pnr_windows), selection=ChannelSelection(include=("HHE",)) + ) + assert len(staged.contributing) == 14 + assert all(id.endswith("HHE") for id in staged.contributing) + + def test_an_exclusion_says_why_and_at_which_level(self, pnr_windows: Any) -> None: + """ "Excluded" is not actionable; "matched exclude='AQ04' at station" is. + + It also disambiguates the case the docstring warns about — a pattern + that could match at two levels reports the one it hit. + """ + staged = fit_event( + _spectra(pnr_windows), selection=ChannelSelection(exclude=("AQ04",)) + ) + reasons = {id: why for id, why in staged.excluded.items() if "matched" in why} + assert reasons == { + "UR.AQ04.00.HHE": "matched exclude='AQ04' at station", + "UR.AQ04.00.HHN": "matched exclude='AQ04' at station", + } + + def test_selection_changes_the_event_value_it_is_meant_to_change( + self, pnr_windows: Any + ) -> None: + """Otherwise the whole feature is decoration. + + The nearest station carries the most weight by construction, so + dropping it is the single largest lever selection has — and this is + why the option has to exist rather than being an average over whatever + happened to be recorded. + """ + spectra = _spectra(pnr_windows) + everything = fit_event(spectra) + without = fit_event(spectra, selection=ChannelSelection(exclude=("AQ04",))) + + assert len(without.contributing) == len(everything.contributing) - 2 + moved = abs(without.value / everything.value - 1) + assert moved > 0.1, ( + f"dropping the nearest station moved the event value by only " + f"{100 * moved:.1f}%; if the weighting has genuinely become that " + "insensitive this test should say so instead" + ) + + def test_require_pass_is_nearly_inert_under_the_shipped_minimiser( + self, pnr_windows: Any + ) -> None: + """The trap, pinned so it cannot drift into being a surprise. + + `pass_fitting` asks whether ``value +/- stderr`` reaches a bound. + Powell estimates no covariance matrix, so the spread is zero and the + test degenerates to "is the value exactly on the bound" — which + essentially never fires. `leastsq` reports uncertainties and fails six. + + The consequence is that changing the minimiser changes *which stations + vote*, not just how each is fitted, and a naive comparison of two + minimisers compares two ensembles. That is worth an assertion rather + than a docstring, because the failure is a plausible-looking number. + """ + spectra = _spectra(pnr_windows) + assert len(fit_event(spectra, method="powell").contributing) == 28 + assert len(fit_event(spectra, method="leastsq").contributing) == 22 + + def test_holding_the_ensemble_fixed_is_what_makes_minimisers_comparable( + self, pnr_windows: Any + ) -> None: + """The module docstring's 0.6% is only true like-for-like. + + Under the default selection the same comparison gives 144%, and the + difference is entirely the six stations `leastsq` drops and Powell + keeps — not a property of either minimiser. + """ + spectra = _spectra(pnr_windows) + same = ChannelSelection(require_pass=False) + powell = fit_event(spectra, method="powell", selection=same) + leastsq = fit_event(spectra, method="leastsq", selection=same) + + assert len(powell.contributing) == len(leastsq.contributing) == 28 + assert abs(powell.value / leastsq.value - 1) < 0.01 + + default = abs( + fit_event(spectra, method="powell").value + / fit_event(spectra, method="leastsq").value + - 1 + ) + assert default > 1.0, ( + "the two ensembles no longer disagree wildly; if selection has " + "changed so that this is no longer a trap, say so here instead" + ) + + def test_require_pass_drops_a_fit_resting_on_its_bound( + self, pnr_windows: Any + ) -> None: + """A pinned parameter reports the bound, not a measurement, so + averaging it in averages in a constant. + + ``leastsq`` reports uncertainties and so fails several stations on the + bound check, which is what makes it the minimiser to test this with. + """ + spectra = _spectra(pnr_windows) + strict = fit_event(spectra, method="leastsq") + loose = fit_event( + spectra, method="leastsq", selection=ChannelSelection(require_pass=False) + ) + + assert len(strict.contributing) < len(loose.contributing) + assert any("stage-1 fit failed" in why for why in strict.excluded.values()) + + def test_a_station_with_no_band_is_excluded_with_that_reason(self) -> None: + """Not silently absent. A station that never got a fit and a station + that was deselected are different outcomes, and the caller needs to be + able to tell them apart without re-deriving either.""" + import pandas as pd # noqa: PLC0415 + + @dataclass + class Fit: + """Just enough of `FitSpectra` for selection to read.""" + + spectra: dict[str, object] + models: dict[str, object] + table: Any + + fit = Fit(spectra={"XX.TEST..HHZ": object()}, models={}, table=pd.DataFrame([])) + contributing, excluded = ChannelSelection().choose(fit) # type: ignore[arg-type] + assert contributing == [] + assert "no fit" in excluded["XX.TEST..HHZ"] + + +class TestWhenNobodyVotes: + def test_stage_two_is_skipped_rather_than_invented(self, pnr_windows: Any) -> None: + """Fixing the corner to a mean of no stations would be inventing the + number the second stage exists to constrain.""" + staged = fit_event( + _spectra(pnr_windows), selection=ChannelSelection(include=("NOTHING",)) + ) + assert staged.stage2 is None + assert np.isnan(staged.value) + assert staged.contributing == () + assert len(staged.excluded) == 28 + # Stage one is still there, so the caller has something to look at. + assert len(staged.stage1.models) == 28 + assert staged.table is staged.stage1.table + + def test_describe_still_says_something_useful(self, pnr_windows: Any) -> None: + staged = fit_event( + _spectra(pnr_windows), selection=ChannelSelection(include=("NOTHING",)) + ) + text = staged.describe() + assert "from 0 channels" in text + assert "28 excluded" in text + + +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_hypocentral_distance"), InverseDistance + ) + with pytest.raises(ValueError, match="Unknown weighting"): + get_weight_model("nope") + + def test_every_registered_weighting_is_usable(self, pnr_windows: Any) -> None: + """A registry entry that cannot be selected is decoration. + + ``inverse_variance`` is exercised separately: it needs uncertainties, + which the configured minimiser does not produce. + """ + spectra = _spectra(pnr_windows) + for name in WEIGHT_MODELS: + if name == "inverse_variance": + continue + assert fit_event(spectra, weighting=name).value > 0, name + + def test_the_choice_of_weighting_moves_the_answer(self, pnr_windows: Any) -> None: + """Which is why it is a choice and not a constant.""" + spectra = _spectra(pnr_windows) + by_distance = fit_event(spectra, weighting="inverse_hypocentral_distance").value + uniform = fit_event(spectra, weighting="uniform").value + assert abs(by_distance / uniform - 1) > 0.05 + + def test_inverse_variance_refuses_rather_than_falling_back( + self, pnr_windows: Any + ) -> None: + """Powell estimates no covariance matrix, so there is nothing to weight + by. Silently becoming uniform is the kind of substitution that ends up + in a paper.""" + with pytest.raises(ValueError, match="uncertainty on every"): + fit_event(_spectra(pnr_windows), weighting="inverse_variance") + + def test_inverse_variance_works_where_the_uncertainties_exist( + self, pnr_windows: Any + ) -> None: + value = fit_event( + _spectra(pnr_windows), weighting="inverse_variance", method="leastsq" + ).value + assert value > 0 + + def test_distance_weighting_names_the_missing_geometry(self) -> None: + """Rather than weighting it zero, which would quietly drop a station.""" + + @dataclass + class Signal: + meta: dict[str, Any] + + @dataclass + class Pair: + signal: Signal + + spectra = {"XX.A..HHZ": Pair(signal=Signal(meta={}))} + with pytest.raises(ValueError, match="set_stream_distance"): + InverseDistance().weights(None, spectra, ["XX.A..HHZ"]) + + +class TestReporting: + def test_spread_reports_what_the_mean_hides(self, pnr_windows: Any) -> None: + """A 2% spread and a 300% spread give the same weighted mean and mean + very different things.""" + spread = _event(pnr_windows).spread() + assert spread["n"] == 28 + assert spread["min"] < spread["median"] < spread["max"] + assert spread["relative_spread"] > 1.0 + + def test_describe_names_the_weighting_and_the_count(self, pnr_windows: Any) -> None: + text = _event(pnr_windows).describe() + assert "from 28 channels" in text + assert "inverse_distance" in text + assert "stage-1 range" in text From 49107ccf21dc9b80434b638ac4055e865b574dd3 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 8 Aug 2026 16:11:36 +0000 Subject: [PATCH 2/4] docs(tutorial): show the API and the fifteen lines it replaces, side by side MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- Tutorial/SpecModTutorial.ipynb | 580 ++++++++++++++++++++++++--------- docs/REFACTOR_PLAN.md | 181 +++++++++- 2 files changed, 594 insertions(+), 167 deletions(-) diff --git a/Tutorial/SpecModTutorial.ipynb b/Tutorial/SpecModTutorial.ipynb index c121055..38c4665 100644 --- a/Tutorial/SpecModTutorial.ipynb +++ b/Tutorial/SpecModTutorial.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "d4c61a00", + "id": "39eda767", "metadata": {}, "source": [ "# SpecMod tutorial\n", @@ -29,7 +29,7 @@ }, { "cell_type": "markdown", - "id": "dda32967", + "id": "59f58548", "metadata": {}, "source": [ "## 1. Read and prepare" @@ -38,13 +38,13 @@ { "cell_type": "code", "execution_count": 1, - "id": "e9b06ca7", + "id": "e3afbe8f", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:20.892893Z", - "iopub.status.busy": "2026-08-07T14:52:20.892669Z", - "iopub.status.idle": "2026-08-07T14:52:21.926810Z", - "shell.execute_reply": "2026-08-07T14:52:21.925246Z" + "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" } }, "outputs": [], @@ -65,7 +65,7 @@ }, { "cell_type": "markdown", - "id": "86ea830f", + "id": "4c94213f", "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": "8969b1ee", + "id": "f1bc897f", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:21.930254Z", - "iopub.status.busy": "2026-08-07T14:52:21.929791Z", - "iopub.status.idle": "2026-08-07T14:52:21.934333Z", - "shell.execute_reply": "2026-08-07T14:52:21.933091Z" + "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" } }, "outputs": [], @@ -100,13 +100,13 @@ { "cell_type": "code", "execution_count": 3, - "id": "35e12155", + "id": "732e9a38", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:21.936253Z", - "iopub.status.busy": "2026-08-07T14:52:21.936073Z", - "iopub.status.idle": "2026-08-07T14:52:22.051058Z", - "shell.execute_reply": "2026-08-07T14:52:22.048891Z" + "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" } }, "outputs": [ @@ -131,13 +131,13 @@ { "cell_type": "code", "execution_count": 4, - "id": "e543680a", + "id": "86bd514a", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:22.053721Z", - "iopub.status.busy": "2026-08-07T14:52:22.053186Z", - "iopub.status.idle": "2026-08-07T14:52:22.063980Z", - "shell.execute_reply": "2026-08-07T14:52:22.062654Z" + "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" } }, "outputs": [ @@ -176,7 +176,7 @@ }, { "cell_type": "markdown", - "id": "6973a219", + "id": "fac2ce5b", "metadata": {}, "source": [ "**Instrument correction happens here, in the notebook, not inside the\n", @@ -192,13 +192,13 @@ { "cell_type": "code", "execution_count": 5, - "id": "3d060faf", + "id": "80933dcc", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:22.066128Z", - "iopub.status.busy": "2026-08-07T14:52:22.065896Z", - "iopub.status.idle": "2026-08-07T14:52:22.663818Z", - "shell.execute_reply": "2026-08-07T14:52:22.662723Z" + "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" } }, "outputs": [ @@ -239,13 +239,13 @@ { "cell_type": "code", "execution_count": 6, - "id": "1dc3fa33", + "id": "a4af3926", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:22.666080Z", - "iopub.status.busy": "2026-08-07T14:52:22.665820Z", - "iopub.status.idle": "2026-08-07T14:52:24.636889Z", - "shell.execute_reply": "2026-08-07T14:52:24.635427Z" + "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" } }, "outputs": [], @@ -256,7 +256,7 @@ }, { "cell_type": "markdown", - "id": "2d07a061", + "id": "23e4c246", "metadata": {}, "source": [ "## 2. Cut the windows" @@ -264,7 +264,7 @@ }, { "cell_type": "markdown", - "id": "d566b6a6", + "id": "a9667906", "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": "cb92e730", + "id": "56e0fd1c", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:24.639248Z", - "iopub.status.busy": "2026-08-07T14:52:24.639045Z", - "iopub.status.idle": "2026-08-07T14:52:24.664232Z", - "shell.execute_reply": "2026-08-07T14:52:24.662571Z" + "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" } }, "outputs": [], @@ -300,13 +300,13 @@ { "cell_type": "code", "execution_count": 8, - "id": "44ec381a", + "id": "a4fe9b32", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:24.666266Z", - "iopub.status.busy": "2026-08-07T14:52:24.666072Z", - "iopub.status.idle": "2026-08-07T14:52:24.671009Z", - "shell.execute_reply": "2026-08-07T14:52:24.669953Z" + "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" } }, "outputs": [ @@ -338,13 +338,13 @@ { "cell_type": "code", "execution_count": 9, - "id": "adc23a0a", + "id": "e9929541", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:24.673064Z", - "iopub.status.busy": "2026-08-07T14:52:24.672887Z", - "iopub.status.idle": "2026-08-07T14:52:27.248790Z", - "shell.execute_reply": "2026-08-07T14:52:27.247611Z" + "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" } }, "outputs": [], @@ -355,7 +355,7 @@ }, { "cell_type": "markdown", - "id": "e939c79b", + "id": "a8dd38a4", "metadata": {}, "source": [ "## 3. Spectra and bandwidth" @@ -363,7 +363,7 @@ }, { "cell_type": "markdown", - "id": "fa5b9b49", + "id": "e3fee1df", "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": "75621d26", + "id": "87aebed5", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:27.251320Z", - "iopub.status.busy": "2026-08-07T14:52:27.251032Z", - "iopub.status.idle": "2026-08-07T14:52:27.507370Z", - "shell.execute_reply": "2026-08-07T14:52:27.506427Z" + "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" } }, "outputs": [ @@ -408,13 +408,13 @@ { "cell_type": "code", "execution_count": 11, - "id": "a1801a1e", + "id": "b93a32ef", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:27.510229Z", - "iopub.status.busy": "2026-08-07T14:52:27.510031Z", - "iopub.status.idle": "2026-08-07T14:52:27.536913Z", - "shell.execute_reply": "2026-08-07T14:52:27.535711Z" + "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" } }, "outputs": [], @@ -431,13 +431,13 @@ { "cell_type": "code", "execution_count": 12, - "id": "d5ccc5d0", + "id": "a93bce39", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:27.539193Z", - "iopub.status.busy": "2026-08-07T14:52:27.538998Z", - "iopub.status.idle": "2026-08-07T14:52:27.542873Z", - "shell.execute_reply": "2026-08-07T14:52:27.541858Z" + "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" } }, "outputs": [ @@ -460,7 +460,7 @@ }, { "cell_type": "markdown", - "id": "3feac95e", + "id": "dd6078f8", "metadata": {}, "source": [ "### Changing ground-motion domain\n", @@ -473,13 +473,13 @@ { "cell_type": "code", "execution_count": 13, - "id": "1b926764", + "id": "bc8fef2a", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:27.544678Z", - "iopub.status.busy": "2026-08-07T14:52:27.544492Z", - "iopub.status.idle": "2026-08-07T14:52:27.608587Z", - "shell.execute_reply": "2026-08-07T14:52:27.607428Z" + "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" } }, "outputs": [ @@ -489,7 +489,7 @@ "text": [ "velocity m/s*s\n", "displacement m*s\n", - "bands that moved: 8 of 28\n" + "bands that moved: 0 of 28\n" ] } ], @@ -507,7 +507,7 @@ }, { "cell_type": "markdown", - "id": "f90c7399", + "id": "1d8f7fea", "metadata": {}, "source": [ "## 4. Fit a source model" @@ -515,7 +515,7 @@ }, { "cell_type": "markdown", - "id": "3c568457", + "id": "3dfca7b2", "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": "f4581e61", + "id": "b5ac93c1", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:27.610682Z", - "iopub.status.busy": "2026-08-07T14:52:27.610397Z", - "iopub.status.idle": "2026-08-07T14:52:27.666850Z", - "shell.execute_reply": "2026-08-07T14:52:27.665625Z" + "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" } }, "outputs": [ @@ -563,13 +563,13 @@ { "cell_type": "code", "execution_count": 15, - "id": "2a98d798", + "id": "f1e68721", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:27.669495Z", - "iopub.status.busy": "2026-08-07T14:52:27.668895Z", - "iopub.status.idle": "2026-08-07T14:52:28.851006Z", - "shell.execute_reply": "2026-08-07T14:52:28.849488Z" + "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" } }, "outputs": [ @@ -589,7 +589,7 @@ }, { "cell_type": "markdown", - "id": "c39545e3", + "id": "6a7eb72f", "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": "8f9aa834", + "id": "03f89e94", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:28.853192Z", - "iopub.status.busy": "2026-08-07T14:52:28.852986Z", - "iopub.status.idle": "2026-08-07T14:52:28.873114Z", - "shell.execute_reply": "2026-08-07T14:52:28.871936Z" + "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" } }, "outputs": [ @@ -725,13 +725,13 @@ { "cell_type": "code", "execution_count": 17, - "id": "036b200b", + "id": "15af287a", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:28.875116Z", - "iopub.status.busy": "2026-08-07T14:52:28.874928Z", - "iopub.status.idle": "2026-08-07T14:52:28.896631Z", - "shell.execute_reply": "2026-08-07T14:52:28.895465Z" + "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" } }, "outputs": [], @@ -744,13 +744,13 @@ { "cell_type": "code", "execution_count": 18, - "id": "667c822c", + "id": "704b062a", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:28.898682Z", - "iopub.status.busy": "2026-08-07T14:52:28.898452Z", - "iopub.status.idle": "2026-08-07T14:52:34.187434Z", - "shell.execute_reply": "2026-08-07T14:52:34.186358Z" + "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" } }, "outputs": [], @@ -763,13 +763,13 @@ { "cell_type": "code", "execution_count": 19, - "id": "edaaccd9", + "id": "c8d5cdd1", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:34.191011Z", - "iopub.status.busy": "2026-08-07T14:52:34.190788Z", - "iopub.status.idle": "2026-08-07T14:52:34.204864Z", - "shell.execute_reply": "2026-08-07T14:52:34.203401Z" + "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" } }, "outputs": [ @@ -922,7 +922,7 @@ }, { "cell_type": "markdown", - "id": "550495fb", + "id": "4f76ba73", "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": "876a3b3d", + "id": "4c83be9e", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:34.207040Z", - "iopub.status.busy": "2026-08-07T14:52:34.206852Z", - "iopub.status.idle": "2026-08-07T14:52:35.792189Z", - "shell.execute_reply": "2026-08-07T14:52:35.790595Z" + "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" } }, "outputs": [ @@ -1086,7 +1086,7 @@ }, { "cell_type": "markdown", - "id": "798fde02", + "id": "0424928a", "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": "d9a40ef2", + "id": "372771f9", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:35.794291Z", - "iopub.status.busy": "2026-08-07T14:52:35.794006Z", - "iopub.status.idle": "2026-08-07T14:52:35.800106Z", - "shell.execute_reply": "2026-08-07T14:52:35.798886Z" + "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" } }, "outputs": [ @@ -1139,13 +1139,13 @@ { "cell_type": "code", "execution_count": 22, - "id": "7996a6d9", + "id": "18e64e8f", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:35.801970Z", - "iopub.status.busy": "2026-08-07T14:52:35.801793Z", - "iopub.status.idle": "2026-08-07T14:52:35.823777Z", - "shell.execute_reply": "2026-08-07T14:52:35.822399Z" + "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" } }, "outputs": [], @@ -1161,7 +1161,7 @@ }, { "cell_type": "markdown", - "id": "128fc32d", + "id": "f52d1081", "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": "6d673735", + "id": "ea8f5a9b", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:35.825818Z", - "iopub.status.busy": "2026-08-07T14:52:35.825634Z", - "iopub.status.idle": "2026-08-07T14:52:35.831154Z", - "shell.execute_reply": "2026-08-07T14:52:35.829964Z" + "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" } }, "outputs": [ @@ -1214,7 +1214,7 @@ }, { "cell_type": "markdown", - "id": "ded1305c", + "id": "f10f5659", "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": "b376739a", + "id": "82e8c326", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:35.833289Z", - "iopub.status.busy": "2026-08-07T14:52:35.833075Z", - "iopub.status.idle": "2026-08-07T14:52:35.847322Z", - "shell.execute_reply": "2026-08-07T14:52:35.846282Z" + "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" } }, "outputs": [ @@ -1287,7 +1287,7 @@ }, { "cell_type": "markdown", - "id": "4b9aaf12", + "id": "be9c2b9a", "metadata": {}, "source": [ "A factor of three in stress drop at the worst station becomes **under 2%**\n", @@ -1302,13 +1302,13 @@ { "cell_type": "code", "execution_count": 25, - "id": "d82c8ed2", + "id": "d0e5ba30", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:35.849663Z", - "iopub.status.busy": "2026-08-07T14:52:35.849460Z", - "iopub.status.idle": "2026-08-07T14:52:36.921262Z", - "shell.execute_reply": "2026-08-07T14:52:36.919736Z" + "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" } }, "outputs": [ @@ -1441,7 +1441,7 @@ }, { "cell_type": "markdown", - "id": "0c01e350", + "id": "7f80cf74", "metadata": {}, "source": [ "0.3% in $t^*$ and about 0.002 in $\\log_{10}\\Omega$ — which is 0.003 magnitude\n", @@ -1463,7 +1463,267 @@ }, { "cell_type": "markdown", - "id": "5be0ba96", + "id": "4ab068b5", + "metadata": {}, + "source": [ + "### The same thing, as one call\n", + "\n", + "Everything above is the workflow written out, and it is written out on purpose\n", + "— the two stages and the weighted mean between them are the method, not an\n", + "implementation detail, and a reader who has not seen them cannot judge a\n", + "result that came out of them.\n", + "\n", + "But nobody should have to retype it. `specmod.staged.fit_event` is those\n", + "fifteen lines with every choice defaulted from `[fitting]`, and it should\n", + "reproduce the number we just computed by hand exactly. If it does not, one of\n", + "the two is wrong." + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "id": "0082031d", + "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" + } + }, + "outputs": [ + { + "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", + "\n", + "by hand : 12.7512 Hz\n", + "API : 12.7512 Hz\n", + "agree : True\n" + ] + } + ], + "source": [ + "from specmod.staged import ChannelSelection, fit_event\n", + "\n", + "staged = fit_event(spectra) # the whole thing, configured defaults\n", + "print(staged.describe())\n", + "print()\n", + "print(f\"by hand : {event_fc['powell']:.4f} Hz\")\n", + "print(f\"API : {staged.value:.4f} Hz\")\n", + "print(f\"agree : {abs(staged.value / event_fc['powell'] - 1) < 1e-12}\")" + ] + }, + { + "cell_type": "markdown", + "id": "a1c9c1a4", + "metadata": {}, + "source": [ + "`describe()` prints the spread as well as the mean, and that is deliberate. A\n", + "2% spread and a 300% spread give the same weighted mean and mean completely\n", + "different things; a corner frequency reported without it is a number with no\n", + "error on it.\n", + "\n", + "### Which stations vote is a decision, and a big one\n", + "\n", + "Everything so far has averaged over all 28 channels. That is rarely what you\n", + "want after looking at the data. A clipped record, a bad instrument response or\n", + "a pick on the wrong phase gives a corner frequency that is *confidently*\n", + "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." + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "id": "23a48dfd", + "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" + } + }, + "outputs": [ + { + "name": "stdout", + "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", + "\n", + "the single nearest is UR.AQ04.00.HHN at 2.30 km\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "\n", + "ids = list(staged.contributing)\n", + "distance = np.array([spectra[i].signal.meta[\"rhyp\"] for i in ids])\n", + "w = 1 / distance\n", + "w = w / w.sum()\n", + "order = np.argsort(-w)\n", + "\n", + "print(\"weight carried by the nearest channels:\")\n", + "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\")" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "id": "4ad2371c", + "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" + } + }, + "outputs": [ + { + "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", + "\n", + " UR.AQ04.00.HHE: matched exclude='AQ04' at station\n", + " UR.AQ04.00.HHN: matched exclude='AQ04' at station\n" + ] + } + ], + "source": [ + "# Drop that station. A bare station code matches every channel it has —\n", + "# `\"HHE\"` would match a component, `\"UR\"` a network, `\"UR.AQ04.00.HHE\"` one\n", + "# channel.\n", + "# `nearest`, not `station` — that name is bound to a full trace id further up\n", + "# and is read again when the spectra are saved.\n", + "nearest = ids[order[0]].split(\".\")[1]\n", + "without = fit_event(spectra, selection=ChannelSelection(exclude=(nearest,)))\n", + "\n", + "print(f\"all channels fc = {staged.value:6.3f} Hz ({len(staged.contributing)} channels)\")\n", + "print(f\"without {nearest:7s} fc = {without.value:6.3f} Hz ({len(without.contributing)} channels)\")\n", + "ratio = without.value / staged.value\n", + "print(f\" change in fc {100 * (ratio - 1):+.1f}%\")\n", + "print(f\" change in stress drop {ratio**3:.2f}x\")\n", + "print()\n", + "for id, why in sorted(without.excluded.items()):\n", + " print(f\" {id}: {why}\")" + ] + }, + { + "cell_type": "markdown", + "id": "02dddca6", + "metadata": {}, + "source": [ + "One quality-control decision, a factor of 1.5 in stress drop. That is not an\n", + "argument against making the decision — it is an argument for making it\n", + "deliberately, writing it into the study file rather than a notebook cell, and\n", + "reporting it. `exclude` lives in `[fitting]` for exactly that reason, and every\n", + "exclusion comes back with the reason and the level it matched at.\n", + "\n", + "### One trap worth knowing about\n", + "\n", + "`require_pass` drops a station whose stage-1 fit ended with a parameter pinned\n", + "against one of its bounds — the minimiser saying \"further, if you would let\n", + "me\", with the bound reported instead of a measurement. Sensible. But the test\n", + "is whether $value \\pm \\sigma$ reaches the bound, and Powell estimates no\n", + "covariance matrix, so $\\sigma$ is missing and the test almost never fires." + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "id": "79a89ac6", + "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" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "powell 28 channels vote, event fc 12.751 Hz\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "leastsq 22 channels vote, event fc 5.228 Hz\n", + "\n", + "with the ensemble held fixed at all 28:\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "powell 28 channels vote, event fc 12.751 Hz\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "leastsq 28 channels vote, event fc 12.675 Hz\n" + ] + } + ], + "source": [ + "for method in (\"powell\", \"leastsq\"):\n", + " run = fit_event(spectra, method=method)\n", + " print(f\"{method:8s} {len(run.contributing):2d} channels vote, event fc {run.value:6.3f} Hz\")\n", + "\n", + "same = ChannelSelection(require_pass=False)\n", + "print()\n", + "print(\"with the ensemble held fixed at all 28:\")\n", + "for method in (\"powell\", \"leastsq\"):\n", + " run = fit_event(spectra, method=method, selection=same)\n", + " print(f\"{method:8s} {len(run.contributing):2d} channels vote, event fc {run.value:6.3f} Hz\")" + ] + }, + { + "cell_type": "markdown", + "id": "5fcd0afb", + "metadata": {}, + "source": [ + "So changing the minimiser changes *which stations vote*, not only how each one\n", + "is fitted. Compared naively the two look 144% apart; compared over the same\n", + "ensemble they agree to 0.6%. Almost all of that gap is the six stations\n", + "`leastsq` rejects and Powell cannot.\n", + "\n", + "Neither setting is wrong. `require_pass=True` is doing the right thing when it\n", + "fires. But a study comparing minimisers, or quoting a corner frequency\n", + "alongside one obtained another way, has to hold the ensemble fixed or say that\n", + "it did not." + ] + }, + { + "cell_type": "markdown", + "id": "a083b0a6", "metadata": {}, "source": [ "## 5. Save the results" @@ -1471,7 +1731,7 @@ }, { "cell_type": "markdown", - "id": "7e85258e", + "id": "a2aebd59", "metadata": {}, "source": [ "Two formats, because the data is used two different ways.\n", @@ -1489,14 +1749,14 @@ }, { "cell_type": "code", - "execution_count": 26, - "id": "b5f3de51", + "execution_count": 30, + "id": "76c8c038", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:36.923434Z", - "iopub.status.busy": "2026-08-07T14:52:36.923199Z", - "iopub.status.idle": "2026-08-07T14:52:37.098650Z", - "shell.execute_reply": "2026-08-07T14:52:37.096738Z" + "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" } }, "outputs": [ @@ -1522,14 +1782,14 @@ }, { "cell_type": "code", - "execution_count": 27, - "id": "035b5fb4", + "execution_count": 31, + "id": "935c8085", "metadata": { "execution": { - "iopub.execute_input": "2026-08-07T14:52:37.101731Z", - "iopub.status.busy": "2026-08-07T14:52:37.100636Z", - "iopub.status.idle": "2026-08-07T14:52:37.139932Z", - "shell.execute_reply": "2026-08-07T14:52:37.138347Z" + "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" } }, "outputs": [ @@ -1539,7 +1799,7 @@ "['2019-08-26T07:49:24.200000Z.csv', '2019-08-26T07:49:24.200000Z.parquet']" ] }, - "execution_count": 27, + "execution_count": 31, "metadata": {}, "output_type": "execute_result" } @@ -1552,7 +1812,7 @@ }, { "cell_type": "markdown", - "id": "1951b880", + "id": "d93d8ed8", "metadata": {}, "source": [ "---\n", diff --git a/docs/REFACTOR_PLAN.md b/docs/REFACTOR_PLAN.md index ebf491a..ddd82df 100644 --- a/docs/REFACTOR_PLAN.md +++ b/docs/REFACTOR_PLAN.md @@ -250,7 +250,7 @@ regenerated with the 1.0 code, or the discrepancy understood. ``` src/specmod/ __init__.py # public API + __version__ - config/ # one module per semantic group (§4.7) + config/ # one module per semantic group (§4.8) layers.py # defaults -> specmod.toml -> *.local.toml -> env -> kwargs provenance.py # resolved config + hash + version, stamped into outputs core/ @@ -951,7 +951,7 @@ queryable with DuckDB or polars **without loading it** — which matters at 11,226 rows and matters more at catalogue scale. CSV stays as an *export*, since journal supplements want it. -**Provenance → both.** The §4.7 record goes into HDF5 attributes *and* a JSON +**Provenance → both.** The §4.8 record goes into HDF5 attributes *and* a JSON sidecar, because the sidecar is greppable and diffable without opening the container. @@ -983,7 +983,7 @@ name time-dependent power spectral densities as an auxiliary-data use case. It is still the wrong **primary**, for two reasons. ASDF is waveform-centric — derived spectra live in the loose `auxiliary_data` bucket, so we would be fitting our data to a schema built for something else. And SEIS-PROV models -*processing* provenance, not the *configuration* provenance of §4.7, so it does +*processing* provenance, not the *configuration* provenance of §4.8, so it does not remove the need for our own record. So: HDF5 with our own schema as primary, and `specmod export --format asdf` @@ -1121,7 +1121,149 @@ preserved — a CSV column carrying one `None` comes back as object-dtype strings, so `pass_fitting` stops being boolean and a downstream `.sum()` counts the wrong thing. -### 4.7 Configuration: semantic groups, layered overrides, recorded provenance +### 4.7 Components, phases, and what the source model assumes about them + +**Not built. This is the largest remaining gap between what the code does and +what the published method does**, and unlike everything else in §4 it is not a +refactor — it is missing capability. + +#### What is there now + +`preprocess.set_stream_distance` computes and stores `azimuth` and +`back_azimuth` on every trace. Nothing reads either. That is the same shape as +the resolution floor before it was fixed: a value tracked from the start with +no reader, and it is the value component rotation needs. + +Everything downstream treats each channel as an independent measurement. On +the PNR event that gives **28 "channels" for 14 stations**, with `HHE` and +`HHN` entering the fit, the flat file and now the event ensemble as if they +were separate observations. They are not; they are two components of one +horizontal motion. + +That has a concrete consequence for §5.2.5's two-stage fit, which is built and +merged. Its ensemble weighting reports that "the nearest two channels carry +22.8% of the weight" — those two channels are one station's two components, so +each station is being counted twice, and a station whose components disagree +contributes that disagreement as if it were between-station scatter. The +weighting is doing what it says; what it is averaging over is wrong. + +#### Combining the horizontals + +The source model is fitted with an average radiation-pattern coefficient, so +what it expects to be given is the **full horizontal S-wave amplitude**, not +one arbitrary component of it. Fitting `HHE` alone underestimates it by an +azimuth-dependent factor, and fitting both independently produces two +different `Omega` for one station. + +The published work combines the two as a geometric mean, `sqrt(A_E * A_N)`, +which is worth noting is *exactly the operation the binner already performs* +in frequency — a bin holds the geometric mean of `log10(amp)` — so the same +convention would apply in two places and should be named once. + +One honest caveat to record now rather than discover later: the geometric mean +of two horizontals is **not** rotation-invariant. It depends on how the +instrument happened to be oriented, which is why ground-motion work moved to +rotation-invariant measures such as RotD50. Whether that matters at the +precision SpecMod reports is a question for measurement, not assertion, and +the answer belongs beside the implementation. + +#### Rotating to transverse + +The alternative, and the other thing to demonstrate: rotate `N`/`E` to +radial/transverse using the stored back-azimuth, and take the S-wave on the +**transverse** component, where the SH energy is. + +These two are not the same choice and should not be presented as +interchangeable. The transverse component isolates SH; the geometric mean +recovers total horizontal motion. Which one is right depends on which +radiation-pattern coefficient the source model was given, and that coupling is +the thing to make explicit — it is currently implicit in a constant. + +#### P-wave source parameters + +Also used in the published work, and needing a **different set of assumptions** +rather than the same model pointed at a different window. At minimum, and each +of these is currently either absent or a hardcoded S-wave value: + +| | S | P | +|---|---|---| +| phase velocity in `M0 = 4*pi*rho*c^3*R*Omega / (R_c * F)` | `beta` | `alpha` | +| average radiation pattern `R_c` | ~0.63 | ~0.52 | +| corner frequency to source radius | `k_S` | `k_P` | +| component the measurement is read on | transverse / horizontal | vertical, or L after LQT | + +`(alpha/beta)^3` is about 5.2 for a Poisson solid, so using the S velocity on +a P measurement is not a small error — it is a factor of five in seismic +moment, which is more than a magnitude unit. + +**The constants above are the standard textbook values and are written here to +show the shape of the problem, not to be adopted.** The Magna work pins +`F = 2` and `Theta-lambda-Phi = 0.55` (§5.2.5), which is already not the +textbook 0.63, so this project's own published choices are what the +implementation should take — read off the papers, recorded in +`studies/*.toml`, and never defaulted silently. A source model that reports a +moment without saying which phase and which constants produced it is the +defect class this whole document is about, applied to the number the package +exists to produce. + +#### Station and channel identity should be a type + +Everything above is hard to say clearly because the code has no notion of a +station. A trace id is a string, and every question about it is asked by +splitting on dots at the point of use. That is already happening in three +places, all written independently in the last week: + +- `staged._levels` splits `NET.STA.LOC.CHA` to let a selection pattern match + at the station or channel level, +- the tutorial does `ids[order[0]].split(".")[1]` to name a station, +- `preprocess.set_picks_from_pyrocko` does + `".".join([tr.stats.network, tr.stats.station])` to build a lookup key. + +Three spellings of one idea, none of which can be given a type. And the bug in +§4.7's opening — components counted as independent stations — is invisible +precisely because nothing in the type system distinguishes "a channel" from +"a station": both are `str`, so summing over channels when you meant stations +type-checks perfectly. + +ObsPy has some of this. `Trace.id` builds the string and `Stream.select` +filters on components, but there is no value type for an identifier, no way to +ask whether two ids are the same station, and no grouping by station. Which is +why every project writes the `split(".")` again. + +What would help, in the same spirit as the `Motion`/`AmplitudeKind` enums of +§4.2 — small, frozen, and doing one thing: + +- `ChannelId`, parsed once from a trace id, with `network`, `station`, + `location`, `channel` and a `station_id` that compares equal across an + instrument's components. Comparable, hashable, and printing back to the + SEED string so it can be a dict key and a column value unchanged. +- Grouping over a set of channels by station, which is what horizontal + combination and the two-stage ensemble both need and both currently lack. +- Pattern matching as a method on the type rather than a free function taking + strings, so `staged._levels` becomes one implementation instead of the + first of several. +- The band/instrument code (`HH`, `EH`, `BH`) and the component letter (`Z`, + `N`, `E`, `1`, `2`) separated, since the `1`/`2` spelling for + non-oriented horizontals is common in real inventories and is exactly the + case string matching on `"HHE"` silently misses. + +This is a prerequisite for §4.7 rather than a parallel nicety: "combine the +horizontals of each station" cannot be written honestly until "the horizontals +of each station" is something the code can express. + +#### Suggested shape + +- `preprocess.rotate_to_rt(st)` reading the stored back-azimuth, and a + horizontal-combination step producing one spectrum per station. +- `phase` as a first-class attribute alongside `motion` and `kind`, so a + spectrum knows whether it is P or S and the model can refuse a mismatch — + the same enum discipline §4.2 applied to units, applied to phases. +- Phase-dependent constants in `[model]`, with the study file supplying them. +- The tutorial demonstrating both horizontal treatments on the same event, in + the way it now demonstrates both minimisers: showing that the choice moves + the answer is what stops it being invisible. + +### 4.8 Configuration: semantic groups, layered overrides, recorded provenance Scientific parameters are currently scattered across three places with no coherent story: a `config.py` of three flat dicts, function defaults in @@ -1534,7 +1676,7 @@ tuned for studies after the paper, so this is drift rather than a defect — but is a trap for step 2 of §5.2.6, because running the current code with stock settings will **not** reproduce the paper. -The fix is §4.7: keep the current values as defaults, and pin the published run +The fix is §4.8: keep the current values as defaults, and pin the published run in `studies/magna_2020_paper.toml`. Each later study gets its own file alongside it, so "which settings produced this" stops being a question anyone has to reconstruct. @@ -1995,6 +2137,31 @@ existing only in the config file. - `intersphinx` to numpy, scipy, obspy, lmfit, matplotlib. - `sphinx.ext.doctest` — the units/normalisation examples in §4.2 and §4.4 are exactly the kind of thing that should be executable in the docs. +**The equations in `docs/` do not render today, and building this is the fix.** +`processing.md` and `choosing_a_transform.md` are written in LaTeX with +`$...$` and `$$...$$`, which is what MyST's `dollarmath` extension reads — and +that extension does not exist yet, because neither does the Sphinx build. The +only renderer these files currently meet is GitHub's, whose math support is +both newer and weaker, so the equations that state the Parseval contract and +the window refinement are being read as literal dollar signs and backslashes. + +Two things to do when this section is built rather than before, since neither +is verifiable without a renderer to check against: + +- Turn on `myst_enable_extensions = ["dollarmath", "amsmath"]`. Without + `dollarmath` MyST does not read `$...$` at all, so adding Sphinx without it + would change nothing. +- Fix the syntax that is wrong independently of the renderer. A scan finds one + display block in `processing.md` without a blank line before it and one + spanning multiple lines, both of which break under MyST as well as GitHub. + The 26 inline expressions containing underscores are the other risk: on + GitHub the emphasis parser can reach them before the math parser does. + +Worth stating the general point, because it is the same shape as §6.6. Prose +that has never been rendered is prose that has never been checked. These files +have been edited a dozen times in this refactor against a renderer nobody has +run. + - `sphinx-build -W` (warnings as errors) in CI: a broken cross-reference fails the build. **Not** an undocumented public symbol, which is what this line used to claim — `-W` promotes warnings that Sphinx already emits, and @@ -2242,7 +2409,7 @@ Each phase ends green on CI and is independently mergeable. |---|---|---|---| | **0. Safety net** | Freeze `master`, default branch → `main`, optional `v0.1.0` tag (§6.7); reproducible legacy env (`Dockerfile`: gfortran + ObsPy 1.2.0 / SciPy 1.4.1 / NumPy 1.18 / pandas 1.0.0 (§5.2.6)); write `datasets/magna_2020.toml` and a first cut of `specmod.acquire`, publish the artifact as a `data-v1` release asset (§5.2); capture golden outputs for PNR **and** Magna; reproduce Table S2 / Figure 2 with 0.1.1 (§5.2.6 step 2); convert any `.spec` files (§4.6) | — | 1.5–2 days | | **1. Make it installable** | `pyproject.toml` + hatch-vcs, `src/` layout, `__init__.py`; ruff config, one-shot `ruff format` + `.git-blame-ignore-revs`, module renames to snake_case; mypy skeleton; pre-commit; `test`/`build` CI; `.gitignore`, `CITATION.cff`; fix the three hard breakages (§1) and the four `F821` bugs ruff finds (§2.5); delete `Tests/Tutorial/`, strip notebook outputs, subset the inventory (§5.1) | 0 | 3–4 days | -| **2. De-globalise** | `config/` package per §4.7 — semantic groups, layer resolution, `config show`/`freeze`, provenance stamping; remove all module-level config reads (tracked by `PLW0603`); `Motion`/`AmplitudeKind` enums; `Spectrum` as a frozen dataclass with `duration`; mutable class attrs (`RUF012`); `isinstance` checks; `logging`. **Tag `v0.2.0`** | 1 | 3–4 days | +| **2. De-globalise** | `config/` package per §4.8 — semantic groups, layer resolution, `config show`/`freeze`, provenance stamping; remove all module-level config reads (tracked by `PLW0603`); `Motion`/`AmplitudeKind` enums; `Spectrum` as a frozen dataclass with `duration`; mutable class attrs (`RUF012`); `isinstance` checks; `logging`. **Tag `v0.2.0`** | 1 | 3–4 days | | **2b. Release plumbing** | Sphinx skeleton + `pydata-sphinx-theme` + autodoc/napoleon/intersphinx; `docs.yml` → GH Pages; release-please + `publish.yml` (PyPI Trusted Publishing); Zenodo webhook. Parallel with 2 | 1 | 1–2 days | | **3. Transform layer** | `SpectralEstimator` protocol; `FFTEstimator`, `WelchEstimator`, `MultitaperEstimator`; `smoothing/` incl. Konno–Ohmachi and `LogBinner`; mtspec demoted to optional legacy backend; Tier 1 + Tier 2 tests; theory docs page | 2 | 5–7 days | | **4. CWT** | `CWTEstimator` + `Scalogram`; COI handling; the Parseval/units calibration and its test; `time_average()`; `ScalogramQC` + the four QC checks; COI floor into `BandwidthSelector`; scalogram plotting; HDF5 scalogram storage | 3 | 6–8 days | @@ -2292,7 +2459,7 @@ end-to-end proves the pipeline while the stakes are zero. export (§4.6). - **Configuration** — semantic groups, layered overrides with local files gitignored by default, resolved config and version stamped into every output - (§4.7). Current behaviour stays the default. + (§4.8). Current behaviour stays the default. - **Tooling** — ruff for lint and format, mypy staged to strict, Sphinx for docs, automated versioning and publishing for both docs and package (§6). - **Branch layout** — `master` frozen as the pre-refactor record, `main` as the new trunk (§6.7). One of the three From d0d257bfa57a6b889b143bab2da85a9c255c01bf Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 8 Aug 2026 16:15:20 +0000 Subject: [PATCH 3/4] docs(plan): ground the station/channel type in FDSN/SEED, and record what checking it found MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- docs/REFACTOR_PLAN.md | 111 ++++++++++++++++++++++++++++++++---------- 1 file changed, 84 insertions(+), 27 deletions(-) diff --git a/docs/REFACTOR_PLAN.md b/docs/REFACTOR_PLAN.md index ddd82df..407e323 100644 --- a/docs/REFACTOR_PLAN.md +++ b/docs/REFACTOR_PLAN.md @@ -1206,7 +1206,7 @@ moment without saying which phase and which constants produced it is the defect class this whole document is about, applied to the number the package exists to produce. -#### Station and channel identity should be a type +#### Station and channel identity should be a type, following FDSN/SEED Everything above is hard to say clearly because the code has no notion of a station. A trace id is a string, and every question about it is asked by @@ -1221,35 +1221,92 @@ places, all written independently in the last week: Three spellings of one idea, none of which can be given a type. And the bug in §4.7's opening — components counted as independent stations — is invisible -precisely because nothing in the type system distinguishes "a channel" from -"a station": both are `str`, so summing over channels when you meant stations -type-checks perfectly. - -ObsPy has some of this. `Trace.id` builds the string and `Stream.select` -filters on components, but there is no value type for an identifier, no way to -ask whether two ids are the same station, and no grouping by station. Which is -why every project writes the `split(".")` again. - -What would help, in the same spirit as the `Motion`/`AmplitudeKind` enums of -§4.2 — small, frozen, and doing one thing: - -- `ChannelId`, parsed once from a trace id, with `network`, `station`, - `location`, `channel` and a `station_id` that compares equal across an - instrument's components. Comparable, hashable, and printing back to the - SEED string so it can be a dict key and a column value unchanged. -- Grouping over a set of channels by station, which is what horizontal - combination and the two-stage ensemble both need and both currently lack. -- Pattern matching as a method on the type rather than a free function taking - strings, so `staged._levels` becomes one implementation instead of the - first of several. -- The band/instrument code (`HH`, `EH`, `BH`) and the component letter (`Z`, - `N`, `E`, `1`, `2`) separated, since the `1`/`2` spelling for - non-oriented horizontals is common in real inventories and is exactly the - case string matching on `"HHE"` silently misses. +precisely because nothing distinguishes "a channel" from "a station": both are +`str`, so summing over channels when you meant stations type-checks perfectly. + +**The conventions already exist and should be followed rather than invented.** +FDSN/SEED defines the hierarchy and the semantics of every field, and each +level carries things SpecMod currently either ignores or re-derives: + +| Level | Code | Carries | +|---|---|---| +| Network | 1–2 chars, FDSN-assigned | operator, and an epoch — temporary network codes are reused, so network alone is not unique in time | +| Station | ≤5 chars | the **site**: latitude, longitude, elevation, and its own epochs | +| Location | 2 chars | which co-located sensor package — surface against borehole, or two instruments at one site | +| Channel | 3 chars | band, instrument, orientation — plus azimuth, dip, sample rate, sensor depth and **the response**, per epoch | + +Four points where that structure matters and the current code does not have it. + +**The channel code is three separate fields.** Band (`H` = high broadband, +`B` = broadband, `E`/`S` = short period, `L` = long period), instrument +(`H` = high-gain seismometer, `N` = accelerometer, `L` = low-gain), and +orientation. Matching `"HHE"` as an opaque string, which +`staged.ChannelSelection` currently does, cannot express "every high-gain +seismometer regardless of band" or "every accelerometer", both of which are +ordinary requests. + +**`N`/`E` versus `1`/`2` is a correctness trap, not a spelling.** The FDSN +convention is that `1`/`2` are used precisely when the horizontals are *not* +aligned to north and east. So code that globs for `HHE`/`HHN` — which this +repository's own tutorial does when reading waveforms — silently finds nothing +on such a network, and code that rotates to radial/transverse by assuming +`N` is 0° is simply wrong there. Rotation needs the channel's own azimuth and +dip from the inventory. + +Worth noting what that means for §4.7's rotation work on the committed data: +**the PNR inventory has `azimuth=None` and `dip=None` on every channel.** The +orientation is not recorded, so rotating this dataset means assuming +`HHN` = 0° and `HHE` = 90° — defensible for channels named `N`/`E`, but an +assumption the code should state rather than bury. Two other things that same +check turned up, both worth fixing whenever the inventory is next touched: +`depth` is `123456.0`, a placeholder rather than a measurement, and the +inventory mixes location codes `''` and `'00'` — which under SEED are +*different* locations, not synonyms. + +**Response and coordinates are epoch-scoped.** An instrument swap creates a +new channel epoch with a different response, and the correct one depends on +the time of the record being corrected. `Inventory.get_response(seed_id, +datetime)` takes the time for that reason. Any station type SpecMod grows has +to keep the epoch rather than flattening to "the response", or it will +silently apply the wrong one across a swap. + +**`azimuth` already means two different things.** `preprocess` writes +`tr.stats["azimuth"]` and `tr.stats["back_azimuth"]` as the *source-receiver* +geometry. `Channel.azimuth` in StationXML is the *component orientation*. +These are unrelated quantities one letter apart in the same namespace, and +rotation needs both at once — which is the moment the collision becomes a bug +rather than a confusion. + +**What ObsPy gives and does not.** `Inventory`/`Network`/`Station`/`Channel` +model the hierarchy well and `get_channel_metadata` resolves an epoch, so the +metadata side is largely solved and should be used rather than reimplemented. +What is missing is on the *identifier*: `Trace.id` builds the string, +`Stream.select` filters with wildcards, and that is all. There is no value +type, no equality at the station level, no grouping by station, and no +structured access to the band/instrument/orientation split. Which is why every +project writes `split(".")` again, including this one, three times. + +So the piece to add is small and specific, in the spirit of the +`Motion`/`AmplitudeKind` enums of §4.2: + +- `ChannelId`, parsed once from a trace id, exposing `network`, `station`, + `location`, `band`, `instrument`, `orientation`, and a `station_id` that + compares equal across an instrument's components. Hashable, and printing + back to the SEED string unchanged so it can be a dict key and a table value. +- Grouping a set of channels by station, which horizontal combination and the + two-stage ensemble both need and both currently lack. +- Matching as a method on the type, so `staged._levels` becomes one + implementation rather than the first of several. +- The FDSN Source Identifier spelling + (`FDSN:NET_STA_LOC_BAND_SOURCE_SUBSOURCE`) is where the standard is going + and separates those fields explicitly; worth parsing to, even while the + dotted SEED form stays what is printed. This is a prerequisite for §4.7 rather than a parallel nicety: "combine the horizontals of each station" cannot be written honestly until "the horizontals -of each station" is something the code can express. +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. #### Suggested shape From 1806ed971738ceff004e162a59084e5044632839 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 8 Aug 2026 16:18:14 +0000 Subject: [PATCH 4/4] docs(plan): separate the units hazard from the radiation-pattern coefficient MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- docs/REFACTOR_PLAN.md | 47 ++++++++++++++++++++++++++++++++++++------- 1 file changed, 40 insertions(+), 7 deletions(-) diff --git a/docs/REFACTOR_PLAN.md b/docs/REFACTOR_PLAN.md index 407e323..629cefc 100644 --- a/docs/REFACTOR_PLAN.md +++ b/docs/REFACTOR_PLAN.md @@ -1198,13 +1198,46 @@ moment, which is more than a magnitude unit. **The constants above are the standard textbook values and are written here to show the shape of the problem, not to be adopted.** The Magna work pins -`F = 2` and `Theta-lambda-Phi = 0.55` (§5.2.5), which is already not the -textbook 0.63, so this project's own published choices are what the -implementation should take — read off the papers, recorded in -`studies/*.toml`, and never defaulted silently. A source model that reports a -moment without saying which phase and which constants produced it is the -defect class this whole document is about, applied to the number the package -exists to produce. +`F = 2` and `Theta-lambda-Phi = 0.55` (§5.2.5), and this project's own +published choices are what the implementation should take — read off the +papers, recorded in `studies/*.toml`, and never defaulted silently. + +**On whether the difference from the textbook 0.63 is a unit conversion.** It +cannot be, at least not in that coefficient: the radiation pattern term is +dimensionless, and a change of units cannot move a dimensionless number. +0.63 is the RMS average of the S-wave radiation pattern over the focal sphere +and other averages of it are also in the literature, so the likelier +explanation is which average was taken — but that is a question for the paper, +not for this document to settle. + +Units are a real hazard here all the same, just in the other terms, and the +exposure is far larger than a factor of 0.87: + +| term | metric | CGS | factor | +|---|---|---|---| +| density | kg/m^3 | g/cm^3 | 10^3 | +| velocity, **cubed** in the moment formula | m/s | km/s | 10^9 | +| distance | m | km | 10^3 | +| the moment itself | N m | dyne cm | 10^7 | + +`r_0 = 1000 m` sitting among the recorded parameters is exactly the shape of a +constant that bridges kilometres and metres, which is worth checking rather +than assuming either way. + +There is a cheap arithmetic check that settles it without reading any code. +`M0 = 10 ** (1.5 * Mw + 9.1)` newton-metres, so the Preston New Road **Mw 1.6** +used throughout this repository should come out near **3e11 N m** — or 3e18 if +the pipeline is working in dyne-centimetres. Anything else is a unit error, +and a factor of 10^9 from a kilometre-per-second velocity is not subtle. + +The general point is that this is ambiguous *now* because the constants live +as bare floats in prose with no units attached — which is the same failure +this document keeps finding, applied to the number the package exists to +produce. §4.2 gave `Motion` and `AmplitudeKind` to spectra so a conversion +could not be applied twice or in the wrong direction. The moment calculation +needs the same discipline: units on the inputs, a declared output unit, and a +test asserting the Mw of a known event. Then the question answers itself +instead of being reasoned about. #### Station and channel identity should be a type, following FDSN/SEED