From 71a4e9c68c60c7ded332b6e2cdecd8da7c0b7a99 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 26 Sep 2026 18:21:35 +0000 Subject: [PATCH 1/3] feat(config): hold a configuration across an ambient read, and pin the goldens MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit REFACTOR_PLAN §6.6 records a claim whose mechanism was absent: "Regression tests pin an explicit config file, never the defaults, so changing a default cannot silently move a golden test." Nothing pinned anything. `load_config()` takes a `project_file`, but the pipeline does not call it that way: `spectrum_set_from_streams` and most of what it calls read settings from functions that accept no configuration argument, so there was nowhere to thread one through. The golden references were therefore captured under whatever the shipped defaults were, and moving a default moved the committed numbers — measured, bumping `smoothing.n_bins` from 151 to 158 broke 9 of 25, and changing the estimator and smoother defaults breaks 11. A deliberate change and an accidental regression were the same failure, wearing the same message. `config.using()` is the mechanism. A `ContextVar` rather than a module global, so a pin does not leak between threads or concurrent tasks and the reset token restores exactly the previous one however the block exits. Only an argument-free `load_config()` consults it: an explicit `project_file=` or override still wins, because a study config read inside a pinned block must still be the file it asked for. The local file and the environment are excluded inside a pin, since a pin a developer's shell could move would not be one. `load_config`'s `start` becomes positional-only, so a section named `start` in `**overrides` cannot collide with it. No caller passed it by keyword. `tests/golden/reference.toml` is the pinned configuration: every setting, not an overlay. An overlay would still inherit the shipped defaults for whatever it omitted, so a later default change to an unlisted key would move the numbers again — the same failure by a narrower path. It is generated from the stamp the references already carried, and round-trips to the same config hash, 7f602da31774. `tools/make_golden.py` captures inside the pin and the golden suite compares inside it. Verified by regenerating: every number reproduces to machine epsilon — 272 of 273 changed leaves at <= 1.4e-15, all BLAS-level — the exception being the known-unstable cwt band edge at `UR.AQ07.00.HHN`, which predates this change. The committed files are unchanged. `to_toml` no longer emits trailing whitespace, on a header's blank lines or after a null placeholder's `=`. It writes files that get committed, and the repository strips trailing whitespace on commit, which left a frozen file differing from what the function would write for it. Co-Authored-By: Claude Opus 5 --- src/specmod/config/__init__.py | 4 +- src/specmod/config/layers.py | 78 +++++++++++++++- src/specmod/config/serialize.py | 11 ++- tests/golden/reference.toml | 115 ++++++++++++++++++++++++ tests/test_config_pinning.py | 155 ++++++++++++++++++++++++++++++++ tests/test_golden_reference.py | 84 +++++++++-------- tools/make_golden.py | 49 ++++++---- 7 files changed, 438 insertions(+), 58 deletions(-) create mode 100644 tests/golden/reference.toml create mode 100644 tests/test_config_pinning.py diff --git a/src/specmod/config/__init__.py b/src/specmod/config/__init__.py index 51705ec..87e2465 100644 --- a/src/specmod/config/__init__.py +++ b/src/specmod/config/__init__.py @@ -10,7 +10,7 @@ from __future__ import annotations -from .layers import LAYER_NAMES, ResolvedConfig, load_config +from .layers import LAYER_NAMES, ResolvedConfig, load_config, pinned, using from .provenance import Provenance, config_hash from .sections import ( AcquireConfig, @@ -41,4 +41,6 @@ "WindowsConfig", "config_hash", "load_config", + "pinned", + "using", ] diff --git a/src/specmod/config/layers.py b/src/specmod/config/layers.py index e71c244..94ca3b0 100644 --- a/src/specmod/config/layers.py +++ b/src/specmod/config/layers.py @@ -12,6 +12,10 @@ irreproducible from the repository alone. That is resolved by recording the *resolved* configuration in every output rather than by forbidding local files; see :mod:`specmod.config.provenance`. + +:func:`using` holds one configuration across the ambient reads the pipeline +makes, which is how a regression test fixes the settings its numbers were +produced under. """ from __future__ import annotations @@ -19,6 +23,9 @@ import ast import os import tomllib +from collections.abc import Iterator +from contextlib import contextmanager +from contextvars import ContextVar from dataclasses import fields from pathlib import Path from typing import Any @@ -29,6 +36,7 @@ "LAYER_NAMES", "ResolvedConfig", "load_config", + "using", ] PROJECT_FILE = "specmod.toml" @@ -69,6 +77,52 @@ def explain(self) -> str: return "\n".join(lines) +#: The configuration :func:`using` is currently holding, if any. A +#: ``ContextVar`` rather than a module global so that a pin does not leak +#: between threads or between concurrent tasks — each gets its own view, and +#: the reset token restores exactly the previous one however the block exits. +_PINNED: ContextVar[ResolvedConfig | None] = ContextVar( + "specmod_pinned_config", default=None +) + + +@contextmanager +def using( + project_file: Path | str | None = None, **overrides: dict[str, Any] +) -> Iterator[ResolvedConfig]: + """Hold one configuration for every ambient read inside the block. + + The pipeline reads settings by calling :func:`load_config` with no + arguments, deep inside functions that take no configuration parameter. + Passing a config down to them is not possible, so a caller that needs to + fix the settings — a regression test, a study script — pins them here:: + + with config.using("tests/golden/reference.toml"): + result = spectrum_set_from_streams(signal, noise, estimator="fft") + + Everything the block calls sees that configuration, and only that: the + local file and the environment are excluded, since a pin whose value still + depended on the machine would not be one. + + Explicit arguments still win — ``load_config(project_file=...)`` inside the + block reads what it was asked for. The pin replaces the *ambient* answer, + not every answer. + """ + resolved = load_config( + project_file=project_file, use_local=False, use_env=False, **overrides + ) + token = _PINNED.set(resolved) + try: + yield resolved + finally: + _PINNED.reset(token) + + +def pinned() -> ResolvedConfig | None: + """The configuration :func:`using` is holding, or ``None``.""" + return _PINNED.get() + + def _read_toml(path: Path) -> dict[str, Any]: with path.open("rb") as fh: return tomllib.load(fh) @@ -116,6 +170,7 @@ def _merge( def load_config( start: Path | str | None = None, + /, *, project_file: Path | str | None = None, use_local: bool = True, @@ -129,7 +184,9 @@ def load_config( start Directory to search for config files. Defaults to the current directory. The search does not walk upwards — an implicit parent search - makes it unclear which file a run actually used. + makes it unclear which file a run actually used. Positional-only, so + that a section named ``start`` in ``**overrides`` cannot collide with + it. project_file Explicit path to a committed config, bypassing the search. This is how a study config (``studies/magna_2020_paper.toml``) is applied, and how @@ -139,7 +196,26 @@ def load_config( so a developer's machine cannot influence a result. **overrides Section-keyed dicts, e.g. ``snr={"tolerance": 4}``. Highest precedence. + + Notes + ----- + A call with no arguments is an *ambient* read, and returns whatever + :func:`using` is holding when it is holding something. Any argument at all + makes the call explicit and resolves the layers as documented above, so a + caller that asks for a particular file gets that file. """ + ambient = ( + start is None + and project_file is None + and use_local + and use_env + and not overrides + ) + if ambient: + held = _PINNED.get() + if held is not None: + return held + root = Path(start) if start is not None else Path.cwd() merged: dict[str, dict[str, Any]] = {} sources: dict[str, str] = {} diff --git a/src/specmod/config/serialize.py b/src/specmod/config/serialize.py index 40dea39..4665eae 100644 --- a/src/specmod/config/serialize.py +++ b/src/specmod/config/serialize.py @@ -40,7 +40,11 @@ def to_toml(config: Config, *, header: str | None = None) -> str: """ lines: list[str] = [] if header: - lines.extend(f"# {line}" for line in header.splitlines()) + # `"#"` rather than `"# "` on a blank line: a comment marker followed + # by a space is trailing whitespace, which the repository's pre-commit + # hook strips — leaving a frozen file that no longer matches what this + # function would write for it. + lines.extend(f"# {line}" if line else "#" for line in header.splitlines()) lines.append("") data = config.to_dict() @@ -49,7 +53,10 @@ def to_toml(config: Config, *, header: str | None = None) -> str: for key in sorted(data[section]): value = data[section][key] if value is None: - lines.append(f"# {key} = ") # TOML has no null + # No trailing space after the `=`: a frozen file is committed, + # and the repository strips trailing whitespace on commit, + # which would leave the file differing from what this writes. + lines.append(f"# {key} =") # TOML has no null elif isinstance(value, dict) and not value: lines.append(f"{key} = {{}}") else: diff --git a/tests/golden/reference.toml b/tests/golden/reference.toml new file mode 100644 index 0000000..848d52c --- /dev/null +++ b/tests/golden/reference.toml @@ -0,0 +1,115 @@ +# The configuration the golden references were captured under. +# +# Every setting, not an overlay: `tests/test_golden_reference.py` pins this file +# with `config.using()`, and a partial file would still inherit the shipped +# defaults for whatever it left out — so a later default change would move the +# committed numbers again, which is the whole failure this closes. +# +# Regenerated only alongside the references themselves, by +# `python tools/make_golden.py`. Changing a value here without regenerating is +# the same as changing the numbers, and the suite says so. + +[acquire] +channels = ["HH?", "BH?", "EN?"] +client = "IRIS" +# depth_km = +# event_id = +# latitude = +locations = ["*"] +# longitude = +# magnitude = +max_radius_km = 400.0 +networks = ["*"] +# origin_time = +remove_response = false +seconds_after = 300.0 +seconds_before = 60.0 +stations = ["*"] + +[fitting] +corner_frequency_min = 0.0 +event_parameter = "fc" +event_weighting = "inverse_distance" +exclude = [] +fit_bins = false +include = [] +initial_alpha = 1e-05 +initial_t_star = 0.01 +method = "powell" +require_pass = true +t_star_min = 0.0001 +weight_method = "none" + +[geometry] +distance_measure = "repi" + +[model] +frequency_dependent_attenuation = false +motion = "velocity" +source = "brune" + +[smoothing] +f_max = 200.0 +f_min = 0.001 +konno_ohmachi_bandwidth = 40.0 +method = "log_bins" +n_bins = 151 +octave_fraction = 0.3333333333333333 +savgol_points_per_decade = 200 +savgol_polyorder = 3 +savgol_window_length = 41 +statistic = "geometric" +window = "hann" + +[snr] +assert_bandwidths = false +bands = [[2.0, 4.0], [4.0, 6.0], [6.0, 8.0]] +bandwidth_method = "peak" +# fixed_band = +interpolate_noise = true +# max_nyquist_fraction = +min_points = 10 +resolution_floor = true +rotate_noise = true +rotation_increment = 0.05 +rotation_method = "boost" +rotation_space = [0.001, 1.001] +scale_parseval = true +tolerance = 3.0 + +[transform] +adaptive = true +dj = 0.125 +drop_dc = true +estimator = "multitaper" +mask_coi = true +# n_fft = +n_tapers = 5 +normalize_to_variance = false +omega0 = 6.0 +taper = "tukey" +taper_alpha = 0.05 +time_bandwidth = 3.0 +wavelet = "morlet" +# welch_segment_length = + +[viz] +plot_columns = 3 + +[windows] +emergency_ratio = 1.7 +noise_length = 1.0 +noise_shift = 0.2 +p_before = 0.0 +p_length = 0.8 +p_length_mode = "relative_time" +p_velocity = 5.9 +pad_seconds = 0.0 +pad_value = 0.0 +refine = true +refine_percentiles = [1.0, 99.0] +s_length = 20.0 +s_length_mode = "absolute_time" +s_start_ratio = 0.8 +s_velocity = 2.9 +station_shifts = {} diff --git a/tests/test_config_pinning.py b/tests/test_config_pinning.py new file mode 100644 index 0000000..073b7ea --- /dev/null +++ b/tests/test_config_pinning.py @@ -0,0 +1,155 @@ +"""Holding a configuration across the ambient reads the pipeline makes. + +`REFACTOR_PLAN` §6.6 records the claim this closes: + + "Regression tests pin an explicit config file, never the defaults, so + changing a default cannot silently move a golden test." + +The mechanism was absent. `load_config()` takes a `project_file`, but the +pipeline does not call it that way — `spectrum_set_from_streams` and the +functions under it read settings ambiently, from call sites that accept no +configuration argument. So a test could pin a config for itself and the code +under test would still read the shipped defaults, and bumping +`smoothing.n_bins` from 151 to 158 moved nine golden tests. +""" + +from __future__ import annotations + +import concurrent.futures +import textwrap +from pathlib import Path + +import pytest + +from specmod import config as cfg + + +@pytest.fixture +def pinned_file(tmp_path: Path) -> Path: + path = tmp_path / "pinned.toml" + path.write_text( + textwrap.dedent(""" + [transform] + estimator = "welch" + + [smoothing] + n_bins = 77 + """) + ) + return path + + +class TestThePinIsSeenByAnAmbientRead: + def test_a_read_with_no_arguments_returns_the_pinned_config( + self, pinned_file: Path + ) -> None: + """The property the whole thing exists for.""" + with cfg.using(pinned_file): + assert cfg.load_config().config.transform.estimator == "welch" + assert cfg.load_config().config.smoothing.n_bins == 77 + + def test_the_default_is_back_after_the_block(self, pinned_file: Path) -> None: + before = cfg.load_config().config.transform.estimator + with cfg.using(pinned_file): + pass + assert cfg.load_config().config.transform.estimator == before + + def test_it_is_restored_even_when_the_block_raises(self, pinned_file: Path) -> None: + before = cfg.load_config().config.transform.estimator + with pytest.raises(RuntimeError), cfg.using(pinned_file): + raise RuntimeError("boom") + assert cfg.load_config().config.transform.estimator == before + + def test_pins_nest(self, pinned_file: Path, tmp_path: Path) -> None: + inner = tmp_path / "inner.toml" + inner.write_text('[transform]\nestimator = "cwt"\n') + with cfg.using(pinned_file): + with cfg.using(inner): + assert cfg.load_config().config.transform.estimator == "cwt" + assert cfg.load_config().config.transform.estimator == "welch" + + +class TestWhatThePinDoesNotCapture: + def test_an_explicit_file_still_wins( + self, pinned_file: Path, tmp_path: Path + ) -> None: + """A caller asking for a particular file gets it, pin or no pin. + + The pin replaces the ambient answer. Anything else would make + `load_config(project_file=...)` lie inside a pinned block, which is + where a study config is read. + """ + other = tmp_path / "other.toml" + other.write_text('[transform]\nestimator = "prieto"\n') + with cfg.using(pinned_file): + resolved = cfg.load_config(project_file=other) + assert resolved.config.transform.estimator == "prieto" + + def test_explicit_overrides_still_win(self, pinned_file: Path) -> None: + with cfg.using(pinned_file): + resolved = cfg.load_config(transform={"estimator": "quadratic"}) + assert resolved.config.transform.estimator == "quadratic" + + def test_the_environment_cannot_reach_inside_a_pin( + self, pinned_file: Path, monkeypatch: pytest.MonkeyPatch + ) -> None: + """A pin a developer's shell could move would not be a pin.""" + monkeypatch.setenv("SPECMOD_TRANSFORM__ESTIMATOR", "cwt") + with cfg.using(pinned_file): + assert cfg.load_config().config.transform.estimator == "welch" + + def test_a_pin_is_not_visible_from_another_thread(self, pinned_file: Path) -> None: + """`ContextVar`, not a module global. + + A pin that leaked across threads would make a parallel test run depend + on which test happened to be inside a block. + """ + + def read() -> str: + return str(cfg.load_config().config.transform.estimator) + + with ( + cfg.using(pinned_file), + concurrent.futures.ThreadPoolExecutor(max_workers=1) as pool, + ): + elsewhere = pool.submit(read).result() + assert elsewhere != "welch" + + +class TestTheBlockReportsWhatItPinned: + def test_it_yields_the_resolved_configuration(self, pinned_file: Path) -> None: + with cfg.using(pinned_file) as resolved: + assert resolved.config.transform.estimator == "welch" + assert resolved is cfg.load_config() + + def test_pinned_says_whether_one_is_held(self, pinned_file: Path) -> None: + assert cfg.pinned() is None + with cfg.using(pinned_file) as resolved: + assert cfg.pinned() is resolved + assert cfg.pinned() is None + + def test_a_pin_with_no_file_is_the_defaults_plus_overrides(self) -> None: + """Pinning the defaults is a real request: it excludes the local file + and the environment, which an ambient read would include.""" + with cfg.using(transform={"estimator": "welch"}) as resolved: + assert resolved.config.transform.estimator == "welch" + assert cfg.load_config().config.transform.estimator == "welch" + + +def test_a_function_that_takes_no_configuration_sees_the_pin(tmp_path: Path) -> None: + """The end the claim is about, and the shape that made it false. + + `plot_columns()` takes no arguments and reads its setting ambiently. There + is no parameter to thread a configuration through, which is true of most + of the settings the pipeline consults; before the pin, such a reader could + only ever see the shipped defaults. + """ + from specmod.fitting.base import plot_columns # noqa: PLC0415 + + default = plot_columns() + path = tmp_path / "viz.toml" + path.write_text(f"[viz]\nplot_columns = {default + 4}\n") + + with cfg.using(path): + assert plot_columns() == default + 4 + assert plot_columns() == default diff --git a/tests/test_golden_reference.py b/tests/test_golden_reference.py index 46dc390..6a3b0f8 100644 --- a/tests/test_golden_reference.py +++ b/tests/test_golden_reference.py @@ -68,6 +68,7 @@ class from ``spectral.py`` into ``core/`` is checked here, on 28 real windows obspy = pytest.importorskip("obspy") import specmod.preprocess as pre # noqa: E402 +from specmod import config as cfg # noqa: E402 from specmod.datasets import PNR_2019 # noqa: E402 from specmod.pipeline import spectrum_set_from_streams # noqa: E402 @@ -75,6 +76,11 @@ class from ``spectral.py`` into ``core/`` is checked here, on 28 real windows _REFERENCE = _ROOT / "tests" / "golden" / "pipeline_reference.json" +#: The settings the references were captured under. Pinned rather than +#: inherited, so that changing a shipped default cannot move a committed +#: number — the property REFACTOR_PLAN §6.6 recorded as claimed but absent. +REFERENCE_CONFIG = _ROOT / "tests" / "golden" / "reference.toml" + # The event and its layout come from `specmod.datasets`, which is also what # `conftest.py` and `tools/make_golden.py` read. This module cannot import the # constants from `conftest.py` directly — a conftest is not an importable @@ -241,29 +247,23 @@ def _flatten(config: dict[str, Any]) -> dict[str, Any]: def test_the_settings_behind_the_reference_are_still_current() -> None: - """Which configuration produced these numbers, and is it still the one. - - The gap this closes is narrow and worth stating exactly. - ``docs/REFACTOR_PLAN.md`` §6.6 records that nothing pins a config here: the - pipeline reads ``load_config()`` ambiently, so this suite runs against - whatever the defaults currently are. The numbers are frozen by the - committed file, which is why moving a default fails the amplitude tests - loudly — 9 of 25 of them, measured, when ``smoothing.n_bins`` was bumped - from 151 to 158. - - What it did not do is say *why*. Nine numeric failures name stations and - ratios, not the setting behind them, and a **regenerated** reference - adopted the new default silently. So the generator now stamps the resolved - configuration into the file, and this test is the assertion over it: the - failure names the key that moved, and a regeneration under changed settings - shows up as a diff in ``_environment.config`` rather than as numbers alone. - - It is not the pin §6.6 asks for. Pinning needs a way to hold a - configuration across an ambient read, which does not exist yet; until then - the protection is "the numbers are frozen, and the settings behind them are - recorded", not "the settings are frozen". + """The pinned file and the stamp in the reference must describe one config. + + ``docs/REFACTOR_PLAN.md`` §6.6 recorded a claim whose mechanism was absent: + "regression tests pin an explicit config file, never the defaults, so + changing a default cannot silently move a golden test". Nothing pinned + anything — the pipeline reads ``load_config()`` ambiently — so bumping + ``smoothing.n_bins`` from 151 to 158 moved nine of these tests, and so did + any deliberate default change, indistinguishably. + + ``reference.toml`` is the pin, ``config.using()`` is what makes an ambient + read see it, and this is the assertion that the file still describes the + configuration the committed numbers were captured under. It compares the + two recorded artefacts, not the shipped defaults: **a default change is now + expected to leave every number here untouched**, and if it does not, the + pin has been bypassed rather than a setting having moved. """ - from specmod.config import config_hash, load_config # noqa: PLC0415 + from specmod.config import config_hash # noqa: PLC0415 stamp = _reference()["_environment"].get("config") assert stamp is not None, ( @@ -272,15 +272,21 @@ def test_the_settings_behind_the_reference_are_still_current() -> None: ) assert stamp["non_default"] == {}, ( f"the reference was generated with {sorted(stamp['non_default'])} set " - "by a layer above the defaults — a local file or the environment, " - "neither of which the repository can reproduce. Regenerate it without." + "by a layer above the pin — a local file or the environment, neither " + "of which the repository can reproduce. Regenerate it without." + ) + assert REFERENCE_CONFIG.is_file(), ( + f"{REFERENCE_CONFIG.name} is missing. It is the configuration these " + "numbers were captured under, and without it the suite falls back to " + "the shipped defaults — which is the failure it exists to prevent." ) - resolved = load_config() - if config_hash(resolved.config) == stamp["hash"]: - return + with cfg.using(REFERENCE_CONFIG) as resolved: + if config_hash(resolved.config) == stamp["hash"]: + return + now = _flatten(resolved.config.to_dict()) - was, now = _flatten(stamp["values"]), _flatten(resolved.config.to_dict()) + was = _flatten(stamp["values"]) moved = [ f" {key}: {was.get(key, '')!r} -> {now.get(key, '')!r}" for key in sorted(set(was) | set(now)) @@ -289,14 +295,15 @@ def test_the_settings_behind_the_reference_are_still_current() -> None: raise AssertionError( "\n".join( [ - f"{len(moved)} setting(s) differ from those the reference was " - "captured under:", + f"{len(moved)} setting(s) in {REFERENCE_CONFIG.name} differ " + "from those the reference was captured under:", *moved, "", - "If the change is intended, the numbers this suite compares " - "against were produced under the old ones. Regenerate with " - "`python tools/make_golden.py` and say in the commit message " - "which numbers moved and by how much.", + "These two are regenerated together. Editing the pinned config " + "without rerunning `python tools/make_golden.py` leaves the " + "committed numbers describing settings nothing uses; " + "regenerate, and say in the commit message which numbers moved " + "and by how much.", ] ) ) @@ -314,7 +321,14 @@ def _run(estimator: str) -> Any: estimators. """ signal, noise = _build_windows() - return spectrum_set_from_streams(signal.copy(), noise.copy(), estimator=estimator) + # Under the settings the reference was captured with, not the shipped + # defaults. Without this the committed numbers move whenever a default + # does, which makes a deliberate default change and an accidental + # regression the same failure — REFACTOR_PLAN §6.6. + with cfg.using(REFERENCE_CONFIG): + return spectrum_set_from_streams( + signal.copy(), noise.copy(), estimator=estimator + ) ESTIMATORS = ["fft", "welch", "multitaper", "quadratic", "cwt"] diff --git a/tools/make_golden.py b/tools/make_golden.py index 0a016f8..185c292 100644 --- a/tools/make_golden.py +++ b/tools/make_golden.py @@ -28,7 +28,8 @@ import scipy import specmod.preprocess as pre -from specmod.config import config_hash, load_config +from specmod import config as cfg +from specmod.config import config_hash from specmod.datasets import PNR_2019 from specmod.pipeline import spectrum_set_from_streams @@ -37,6 +38,12 @@ WINDOWS_OUT = ROOT / "tests" / "golden" / "window_reference.json" MOTION_OUT = ROOT / "tests" / "golden" / "motion_reference.json" +#: The settings every number below is captured under, pinned rather than read +#: from the shipped defaults. Capturing ambiently made the references hostage +#: to the defaults: changing one moved the committed numbers, so a deliberate +#: default change and an accidental regression were the same failure. +REFERENCE_CONFIG = ROOT / "tests" / "golden" / "reference.toml" + #: Read from `specmod.datasets` so that the script generating the reference and #: the test asserting against it cannot disagree about which event, or which #: hypocentre, produced it. @@ -209,27 +216,24 @@ def capture_windows() -> dict: def _configuration() -> dict: - """The settings these numbers were produced under.""" - # Nothing pins a config here — the pipeline reads `load_config()` ambiently - # — so a regenerated reference adopts whatever defaults were current. This - # does not close that gap (REFACTOR_PLAN §6.6); it makes a settings change - # land as a diff in this file and as a named key in the test failure, - # rather than as nine numeric failures with no cause attached. - resolved = load_config() - #: `sources` carries only what a layer above the defaults set, so anything - #: here came from a `specmod.toml`, a gitignored local file or the - #: environment — none of which the repository can reproduce. - non_default = {k: v for k, v in resolved.sources.items() if v != "default"} - if non_default: - print( - f"WARNING: {len(non_default)} setting(s) come from a layer above " - f"the defaults, e.g. {sorted(non_default)[:3]}. Regenerating under " - "an override bakes it into the reference — rerun without it." - ) + """The settings these numbers were produced under. + + Read back from inside the pin, so the stamp records what the capture + actually used rather than what it was meant to use. + """ + resolved = cfg.load_config() + #: Every value comes from `reference.toml`, which `using()` reads with the + #: local file and the environment excluded — so this is empty by + #: construction and stays in the file as the assertion that it is. + above_the_pin = { + k: v + for k, v in resolved.sources.items() + if v not in ("default", str(REFERENCE_CONFIG)) + } return { "hash": config_hash(resolved.config), "values": resolved.config.to_dict(), - "non_default": dict(sorted(non_default.items())), + "non_default": dict(sorted(above_the_pin.items())), } @@ -297,6 +301,13 @@ def _carry_forward_legacy(fresh: dict, section: str, key: str, path: Path) -> No def main() -> None: + # Everything below runs inside the pin, including the capture: the whole + # point is that these numbers do not depend on the shipped defaults. + with cfg.using(REFERENCE_CONFIG): + _capture_all() + + +def _capture_all() -> None: reference = { "_environment": _environment(), **{ From 108136484a55ea5cc6c3490ced0fc8cdaf0052b5 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 26 Sep 2026 18:21:52 +0000 Subject: [PATCH 2/3] feat(config)!: default to FFT with Konno-Ohmachi smoothing MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit REFACTOR_PLAN §8.1's open decision, resolved on the measurement it asked for. Against the shipped multitaper + `log_bins`, on the 28 PNR windows with each station fitted freely: fft + log_bins median dMw -0.039 per-station -0.85 to +1.90 fft + konno_ohmachi median dMw +0.027 per-station -0.52 to +0.51 The smoother earns the change, not the transform. `log_bins` reduces the frequency axis without reducing the variance — 151 bins over a few hundred Fourier samples is about one sample per bin — so an FFT read through it swings a station by up to 1.9 magnitude units. Pairing it with a real smoother halves the worst excursion and pulls the `fc` range in by a third. The median station barely moves either way, which is why the per-station spread is the argument. What this event cannot settle, and what the decision does not rest on: its event corner is unconstrained under all three (stage-one spread 895% to 1042% of the event value). On the ensemble, PNR's Mw moves -0.1615 under both constant sets — a pure Omega shift, the constants being an offset applied after. Against the catalogue's 2.9 the two sets move in opposite directions, the default constants from 0.194 to 0.032 away and §4.7's from 0.160 to 0.321, so they no longer both sit within 0.25. Recorded, not offered as evidence. `!` because every result moves for any configuration that did not pin both settings, while nothing is renamed and no signature changes: the break is in the numbers, not in the call. `docs/upgrading.md` says what to pin to keep the old ones. A coupling the measurement did not predict: `[fitting] fit_bins` fits the binned spectrum, and Konno-Ohmachi preserves the frequency axis, so under the new default the binned axis *is* the Fourier axis and the setting selects the same 98 points either way. It bites only under an axis-reducing smoother. The tests split the way the change does. Those asserting the pipeline's behaviour or a recorded measurement — the golden comparisons, the displacement reference, §4.7's magnitudes, the selection patterns, the `[fitting]` settings — now run under `tests/golden/reference.toml`, so they keep testing their subject rather than whichever pair ships. Those asserting what a user gets with no configuration moved, and `test_the_shipped_default_moves_omega_by_a_constant` is new: without it the default's own answer would be asserted nowhere. Co-Authored-By: Claude Opus 5 --- docs/REFACTOR_PLAN.md | 48 +++++++++++++++-------- docs/roadmap.md | 44 +++++++++++++++++++-- docs/upgrading.md | 63 ++++++++++++++++++++++++++++-- src/specmod/config/sections.py | 13 ++++++- src/specmod/fitting/spectrum.py | 9 +++-- tests/test_config.py | 26 ++++++++++--- tests/test_fitting_defaults.py | 21 +++++++++- tests/test_magnitude.py | 68 ++++++++++++++++++++++++++++++++- tests/test_pipeline.py | 17 ++++++++- tests/test_staged.py | 11 +++++- 10 files changed, 281 insertions(+), 39 deletions(-) diff --git a/docs/REFACTOR_PLAN.md b/docs/REFACTOR_PLAN.md index 16de3a9..627b890 100644 --- a/docs/REFACTOR_PLAN.md +++ b/docs/REFACTOR_PLAN.md @@ -3477,11 +3477,15 @@ bumping `smoothing.n_bins` from 151 to 158 and running the suite — and **9 of 25 golden tests fail**, loudly. The golden values live in a committed JSON file, so any behaviour change breaks them whatever moved it. -The residual gap is narrower than the claim suggests but real: the golden file -records no config, so a *regenerated* reference silently adopts whatever -defaults were current. Pinning a study file, as this plan proposes, is what -would close it. Until then the protection is "the numbers are frozen", not -"the settings are frozen". +**Closed in 0.6.** `specmod.config.using()` holds a configuration across the +ambient reads the pipeline makes — a `ContextVar`, so it does not leak between +threads, and explicit arguments still win — and `tests/golden/reference.toml` +carries every setting the committed numbers were captured under, not an +overlay that would inherit future default changes for whatever it omitted. +`tools/make_golden.py` captures inside the pin and the golden suite compares +inside it, so the property the claim asserted now holds: changing a shipped +default leaves every committed number untouched. Measured by doing it — the +0.6 change of default moved 11 golden tests before the pin and none after. **Claims that were checked and hold:** @@ -3717,10 +3721,8 @@ changes anything shipped. ### Still open -1. **Default estimator.** Multitaper (matching current behaviour) or FFT + - Konno–Ohmachi (more conventional in engineering seismology)? **No longer - blocked** — §6.6's dead `[smoothing] method` is wired up, so the second - option is selectable — and now measured, on the 28 PNR windows with each +1. ~~**Default estimator.**~~ **Resolved: FFT + Konno-Ohmachi**, from 0.6. The + measurement below is what decided it, on the 28 PNR windows with each station fitted freely: | Against the shipped multitaper + `log_bins` | Median ΔMw | Per-station range | `fc` ratio | @@ -3733,17 +3735,33 @@ changes anything shipped. because `log_bins` reduces the axis without reducing the variance — 151 log bins over a few hundred Fourier samples is roughly one sample per bin. Pairing it with a real smoother halves the worst excursion and pulls the - `fc` range in by a third. The median station barely moves either way. - - What this event cannot settle: the *event* corner is unconstrained under all - three (stage-one spread 895% to 1042% of the event value), so the ensemble - number is not evidence. Deciding on the ensemble needs an event where it is - constrained. + `fc` range in by a third. The median station barely moves either way, so the + per-station spread is the whole of the argument. + + What this event cannot settle, and what the decision therefore does not + rest on: the *event* corner is unconstrained under all three (stage-one + spread 895% to 1042% of the event value), so the ensemble number is not + evidence. On the ensemble, PNR's Mw moves −0.1615 under both constant sets + — a pure Omega shift, since the constants are an offset applied afterwards. + Against the catalogue's 2.9 the two sets move in opposite directions, the + default constants from 0.194 to 0.032 away and §4.7's from 0.160 to 0.321, + so they no longer both sit within 0.25. Recorded, not offered as evidence. Speed, measured rather than assumed and far too small to decide on: the transform is 0.209 s against 0.323 s for the 28 windows, and the smoothing adds 0.14 s (`log_bins`) to 0.83 s (`konno_ohmachi`) on top. + **What had to be built first.** Changing a default moved 11 golden tests, + because §6.6's config pin was a claim with no mechanism: the references were + captured under whatever the defaults were. `config.using()` and + `tests/golden/reference.toml` close that, so the references now hold their + own settings and a default change leaves every committed number untouched. + + **A coupling the measurement did not predict.** `[fitting] fit_bins` selects + the binned spectrum, which under an axis-preserving smoother is the Fourier + axis — so with Konno-Ohmachi the setting fits the same 98 points either way + and is inert. It only bites under `log_bins`. + 2. **Python floor.** 3.11 is proposed. Any users stuck on 3.9/3.10? 3. **History rewrite.** Deleting the 9.9 MB catalog and stripping the notebook (§5.1) shrinks the *working tree* but leaves both in history, so a fresh clone diff --git a/docs/roadmap.md b/docs/roadmap.md index 4a71944..1e67407 100644 --- a/docs/roadmap.md +++ b/docs/roadmap.md @@ -191,6 +191,44 @@ found by a live fetch rather than by the suite: as two independent records of one ground motion. `channel_priorities` and `location_priorities` rank them, first match wins per station. +## In progress — not yet released + +Merged on `main`, and it will name its version when a release goes out. + +### The default is FFT with Konno-Ohmachi smoothing + +The shipped pair was multitaper with log-binning. Measured on the 28 PNR +windows, with each station fitted freely, against that pair: + +| | Median ΔMw | Per-station range | `fc` ratio | +|---|---|---|---| +| `fft` + `log_bins` | −0.039 | −0.85 to +1.90 | 0.003 to 5.0 | +| `fft` + `konno_ohmachi` | +0.027 | −0.52 to +0.51 | 0.038 to 3.2 | + +The smoother is what earns the change. `log_bins` reduces the frequency axis +without reducing the variance — 151 bins over a few hundred Fourier samples is +about one sample per bin — so an FFT read through it swings a station by up to +1.9 magnitude units. A real smoother halves the worst excursion and pulls the +corner-frequency range in by a third. The median station barely moves either +way, which is why the per-station spread is the argument and the ensemble +number is not: this event's corner is unconstrained under all three. + +**Every result moves** unless a configuration already pinned both settings. +[Upgrading](upgrading.md#from-05-to-06) has what to pin to keep the old +numbers, and the one coupling this exposes: `[fitting] fit_bins` is inert +under a smoother that preserves the frequency axis. + +### A configuration can be held across an ambient read + +`config.using()` pins one configuration for everything inside the block, +including the many pipeline functions that read their settings from call sites +taking no configuration argument. That was the prerequisite for the change +above: the golden references were captured under whatever the defaults were, +so moving a default moved 11 committed numbers, and a deliberate change was +indistinguishable from a regression. They now carry their own settings in +`tests/golden/reference.toml`, and a default change leaves every one of them +untouched. + ## Planned ### 1.0 — the API stops moving @@ -258,9 +296,9 @@ guess written from documentation cannot substitute for: Each shipped entry names the version that carries it, so a reader can install that version and check the claim. Work that is merged but unreleased sits under *In progress* and gets its version when a release goes out — merged is not -shipped, and this page does not blur the two. **There is no such section right -now**, because nothing is merged and unreleased; that is the state the page -should be in between releases, not a section that went missing. +shipped, and this page does not blur the two. Between releases there is no such +section at all, which is a state the page should be in rather than a section +that went missing. The move at release time is the step this page got wrong once. The v0.3.0 work sat under *In progress — not yet released* for a week after v0.3.0 was diff --git a/docs/upgrading.md b/docs/upgrading.md index 6491ee9..181ed33 100644 --- a/docs/upgrading.md +++ b/docs/upgrading.md @@ -1,15 +1,70 @@ # Upgrading -Four moves are documented here. Most readers porting code want `0.2` to `0.3`, +Five moves are documented here. Most readers porting code want `0.2` to `0.3`, which renamed ten functions in `specmod.preprocess` and changed what they do -with the stream they are given. The two moves above it rename nothing: they are -settings that began doing what they always said, so the break is in the numbers -rather than in the call. +with the stream they are given. The three moves above it rename nothing: they +change what a setting does, or which one is the default, so the break is in the +numbers rather than in the call. `0.x` ships breaking changes in minor bumps without a deprecation cycle, because shimming an API still being worked out costs more than it protects. See [Releasing the software](releasing.md#what-a-version-number-means-while-this-is-0x). +## From 0.5 to 0.6 + +**The default estimator is `fft` and the default smoother is `konno_ohmachi`.** +They were `multitaper` and `log_bins`. Nothing is renamed and no signature +moves; **every result changes** unless your configuration already pinned both. + +On the 28 PNR windows, against the old pair: + +| | Median ΔMw | Per-station range | `fc` ratio | +|---|---|---|---| +| `fft` + `log_bins` | −0.039 | −0.85 to +1.90 | 0.003 to 5.0 | +| `fft` + `konno_ohmachi` | +0.027 | −0.52 to +0.51 | 0.038 to 3.2 | + +The smoother is what earns the change, not the transform. `log_bins` reduces +the frequency axis without reducing the variance — 151 bins over a few hundred +Fourier samples is about one sample per bin — so an FFT paired with it swings a +station by up to 1.9 magnitude units. A real smoother halves the worst +excursion. The median station barely moves either way. + +**To keep the old numbers**, pin both: + +```toml +[transform] +estimator = "multitaper" + +[smoothing] +method = "log_bins" +``` + +Two consequences worth knowing: + +- **`[fitting] fit_bins` is inert under the new default.** It fits the binned + spectrum, and Konno–Ohmachi preserves the frequency axis, so the "binned" + axis *is* the Fourier axis and the setting selects the same points either + way. It only bites under `log_bins`, or any other axis-reducing smoother. +- **Published results should pin the version anyway**, but this is the release + that shows why: a re-run under 0.6 of an analysis done under 0.5 will not + reproduce it, and nothing in the call changed to warn you. + +**New: `specmod.config.using()`**, a context manager that holds one +configuration across the ambient reads the pipeline makes. This is what lets a +regression test fix the settings its numbers were produced under — most of the +pipeline reads settings from functions that take no configuration argument, so +passing one down was not possible before: + +```python +from specmod import config + +with config.using("studies/my_study.toml"): + result = spectrum_set_from_streams(signal, noise) +``` + +The local file and the environment are excluded inside the block, and an +explicit `load_config(project_file=...)` still wins. + ## From 0.4 to 0.5 **`max_radius_km` and `min_radius_km` now mean kilometres.** They were converted diff --git a/src/specmod/config/sections.py b/src/specmod/config/sections.py index 9d501ae..260768d 100644 --- a/src/specmod/config/sections.py +++ b/src/specmod/config/sections.py @@ -99,9 +99,14 @@ class WindowsConfig: class TransformConfig: """Time-to-frequency conversion. Consumed by :mod:`specmod.transforms`.""" + #: FFT since 0.6, paired with the Konno-Ohmachi default in ``[smoothing]``. + #: The pairing is the point: an FFT read through a log-binner carries the + #: periodogram's variance into the fit, and it is the smoother that removes + #: it. ``docs/upgrading.md`` has the measured comparison, and + #: REFACTOR_PLAN §8.1 the decision. estimator: Literal[ "multitaper", "fft", "welch", "cwt", "prieto", "quadratic", "mtspec" - ] = "multitaper" + ] = "fft" #: Multitaper. ``time_bandwidth`` was previously the literal 3 passed #: positionally to mtspec, with no way to configure it. @@ -146,13 +151,17 @@ class SmoothingConfig: method keeps it. See :mod:`specmod.smoothing` for what each one does to a seismic spectrum and ``docs/choosing-a-transform.md`` for the measured comparison. + + That distinction reaches ``[fitting] fit_bins``, which fits the binned + spectrum: under an axis-preserving method there is no separate binned axis, + so the setting selects the same points either way. """ #: Resolved through `specmod.smoothing.SMOOTHERS`. Only the parameters #: belonging to the chosen method are read. method: Literal[ "log_bins", "konno_ohmachi", "log_window", "savitzky_golay", "none" - ] = "log_bins" + ] = "konno_ohmachi" #: Log bin edges, for `log_bins`. ``None`` derives them from the record: #: fmin from 1/T, fmax from Nyquist. The old code hardcoded 0.001-200 Hz diff --git a/src/specmod/fitting/spectrum.py b/src/specmod/fitting/spectrum.py index c89a085..808d42b 100644 --- a/src/specmod/fitting/spectrum.py +++ b/src/specmod/fitting/spectrum.py @@ -173,10 +173,11 @@ def __init_params(self, **params: Any) -> None: ``# self.set_bounds('fc', min=0)`` sat commented out here. It is not a poor fit but an unphysical one: a negative ``t*`` says the wave gained energy travelling, and a corner frequency below zero says nothing at - all. lmfit returns either if the misfit surface leans that way, and - with the shipped multitaper default it returned ``fc = -4.45 Hz`` on - one PNR station while ``pass_fitting`` reported success — because a - parameter with no bound cannot be *at* its bound. + all. lmfit returns either if the misfit surface leans that way: under + multitaper with log-binning, which was the default through 0.5, it + returned ``fc = -4.45 Hz`` on one PNR station while ``pass_fitting`` + reported success — because a parameter with no bound cannot be *at* + its bound. """ # `**params: Any` rather than `float`, because lmfit's `make_params` # takes a leading `verbose` argument: a model with a parameter of that diff --git a/tests/test_config.py b/tests/test_config.py index 339abbd..0fa7edf 100644 --- a/tests/test_config.py +++ b/tests/test_config.py @@ -305,13 +305,29 @@ def test_each_method_resolves( else: assert isinstance(smoother, SMOOTHERS[method]) - def test_the_defaults_still_bin(self) -> None: - """The shipped default is the path every committed number came from.""" + def test_the_shipped_default_smooths_rather_than_bins(self) -> None: + """What a user gets with no configuration at all. + + This used to assert the opposite, on the premise that the default was + "the path every committed number came from". That premise is gone: + the committed numbers now come from `tests/golden/reference.toml`, + pinned, so the default is free to be the one that measures better + rather than the one the references happen to hold. + """ from specmod.pipeline import _compare_settings # noqa: PLC0415 - settings = _compare_settings() - assert settings["smoother"] is None - assert settings["n_bins"] == load_config().config.smoothing.n_bins + assert isinstance(_compare_settings()["smoother"], SMOOTHERS["konno_ohmachi"]) + + def test_the_pinned_reference_config_still_bins(self) -> None: + """And the references keep the path they were captured under.""" + from specmod import config as cfg # noqa: PLC0415 + from specmod.pipeline import _compare_settings # noqa: PLC0415 + + reference = Path(__file__).parent / "golden" / "reference.toml" + with cfg.using(reference): + settings = _compare_settings() + assert settings["smoother"] is None + assert settings["n_bins"] == cfg.load_config().config.smoothing.n_bins def test_the_section_parameters_reach_the_smoother( self, isolated: Path, monkeypatch: pytest.MonkeyPatch diff --git a/tests/test_fitting_defaults.py b/tests/test_fitting_defaults.py index 4c83aa9..bd5f353 100644 --- a/tests/test_fitting_defaults.py +++ b/tests/test_fitting_defaults.py @@ -25,15 +25,24 @@ obspy = pytest.importorskip("obspy") +from pathlib import Path # noqa: E402 + +from specmod import config as cfg # noqa: E402 from specmod.config import load_config # noqa: E402 from specmod.fitting import FitSpectra, initial_guess # noqa: E402 from specmod.pipeline import spectrum_set_from_streams # noqa: E402 +#: The settings the committed references hold, which still bin. These tests +#: assert that a `[fitting]` value is read at all, so they need a spectrum +#: whose shape does not depend on what the shipped smoother happens to be. +REFERENCE_CONFIG = Path(__file__).parent / "golden" / "reference.toml" + @functools.cache def _spectra(windows: Any) -> Any: signal, noise = windows() - return spectrum_set_from_streams(signal, noise, estimator="fft") + with cfg.using(REFERENCE_CONFIG): + return spectrum_set_from_streams(signal, noise, estimator="fft") def _fit(spectra: Any, **kwargs: Any) -> Any: @@ -160,7 +169,15 @@ def test_a_parameter_resting_on_its_bound_fails_the_fit(self) -> None: assert load_config().config.fitting.corner_frequency_min == 0.0 def test_fit_bins_comes_from_configuration(self, pnr_windows: Any) -> None: - """Fitting the binned spectrum uses far fewer points than the raw one.""" + """Fitting the binned spectrum uses far fewer points than the raw one. + + Only where the smoother replaces the frequency axis, which is why this + runs under the pinned reference config rather than the shipped + defaults. Under an axis-preserving smoother — Konno-Ohmachi, which is + the default from 0.6 — the comparison axis *is* the Fourier axis, so + ``fit_bins`` selects the same 98 points either way and the setting is + inert. See ``docs/upgrading.md``. + """ spectra = _spectra(pnr_windows) raw = _fit(spectra, fit_bins=False) binned = _fit(spectra, fit_bins=True) diff --git a/tests/test_magnitude.py b/tests/test_magnitude.py index 68b3201..4d1d6d8 100644 --- a/tests/test_magnitude.py +++ b/tests/test_magnitude.py @@ -2,7 +2,9 @@ from __future__ import annotations +import contextlib import math +from pathlib import Path import numpy as np import pytest @@ -135,6 +137,10 @@ def test_plateau_and_distance_must_correspond(self) -> None: density=2500.0, velocity=2500.0, radiation_pattern=0.63, free_surface=2.0 ) +#: The settings §4.7's recorded magnitudes were measured under, and the ones +#: the golden references hold. +REFERENCE_CONFIG = Path(__file__).parent / "golden" / "reference.toml" + @pytest.mark.usefixtures("pnr_windows") class TestOnTheRealEvent: @@ -147,11 +153,22 @@ class TestOnTheRealEvent: @staticmethod def _staged(windows): # type: ignore[no-untyped-def] + """Under the settings §4.7 recorded its numbers with. + + Pinned, because what these tests pin is the magnitude computation and + the spreading exponent — not which estimator ships as the default. The + values below are quoted from the plan, and a default change that moved + them would make the quotation false without anything about magnitude + having changed. ``test_the_shipped_default_gives_a_different_number`` + is where the default's own answer is asserted. + """ + from specmod import config as cfg # noqa: PLC0415 from specmod.pipeline import spectrum_set_from_streams # noqa: PLC0415 from specmod.staged import fit_event # noqa: PLC0415 signal, noise = windows() - return fit_event(spectrum_set_from_streams(signal, noise)) + with cfg.using(REFERENCE_CONFIG): + return fit_event(spectrum_set_from_streams(signal, noise)) def test_the_event_magnitude_is_unchanged(self, pnr_windows) -> None: # type: ignore[no-untyped-def] event = event_magnitude(self._staged(pnr_windows)) @@ -261,6 +278,7 @@ def test_the_pnr_event_brackets_its_catalogue_magnitude(pnr_windows) -> None: # unlabelled number is what made an earlier version of this look 1.15 magnitude units out. """ + from specmod import config as cfg # noqa: PLC0415 from specmod.datasets import PNR_2019 # noqa: PLC0415 from specmod.pipeline import spectrum_set_from_streams # noqa: PLC0415 from specmod.staged import fit_event # noqa: PLC0415 @@ -270,9 +288,55 @@ def test_the_pnr_event_brackets_its_catalogue_magnitude(pnr_windows) -> None: # assert catalogue == 2.9 signal, noise = pnr_windows() - staged = fit_event(spectrum_set_from_streams(signal, noise)) + with cfg.using(REFERENCE_CONFIG): + staged = fit_event(spectrum_set_from_streams(signal, noise)) low = event_magnitude(staged, constants=PLAN_CONSTANTS).value high = event_magnitude(staged).value assert low < catalogue < high assert abs(low - catalogue) < 0.25 assert abs(high - catalogue) < 0.25 + + +def test_the_shipped_default_moves_omega_by_a_constant(pnr_windows) -> None: # type: ignore[no-untyped-def] + """What a user with no configuration gets, and how far it is from the pin. + + Everything else in this module runs under `tests/golden/reference.toml`, + so without this the default's own answer would be asserted nowhere and + the change of default would be invisible to the suite. + + The shift is **-0.1615 Mw, identical under both constant sets**, which is + the shape a change of estimator and smoother should have: it moves the + low-frequency plateau the fit reads Omega off, and the medium constants + are an offset applied afterwards. A number that moved differently under + the two would mean something had changed in the fitting instead. + + Against the catalogue's 2.9 the two sets move in opposite directions — + the default constants from 0.194 to 0.032 away, the plan's from 0.160 to + 0.321 — so the pair no longer both sit within 0.25. That is one event + whose corner is unconstrained (REFACTOR_PLAN §8.1), so it is recorded + rather than offered as evidence either way. + """ + from specmod import config as cfg # noqa: PLC0415 + from specmod.pipeline import spectrum_set_from_streams # noqa: PLC0415 + from specmod.staged import fit_event # noqa: PLC0415 + + signal, noise = pnr_windows() + + def staged_under(context): # type: ignore[no-untyped-def] + with context: + return fit_event(spectrum_set_from_streams(signal.copy(), noise.copy())) + + pinned = staged_under(cfg.using(REFERENCE_CONFIG)) + shipped = staged_under(contextlib.nullcontext()) + + assert event_magnitude(shipped).value == pytest.approx(2.932, abs=0.01) + assert event_magnitude(shipped, constants=PLAN_CONSTANTS).value == pytest.approx( + 2.579, abs=0.01 + ) + + shifts = [ + event_magnitude(shipped, **kw).value - event_magnitude(pinned, **kw).value + for kw in ({}, {"constants": PLAN_CONSTANTS}) + ] + assert shifts[0] == pytest.approx(-0.1615, abs=0.005) + assert shifts[0] == pytest.approx(shifts[1], abs=1e-9) diff --git a/tests/test_pipeline.py b/tests/test_pipeline.py index 6509623..9fa450e 100644 --- a/tests/test_pipeline.py +++ b/tests/test_pipeline.py @@ -23,6 +23,7 @@ obspy = pytest.importorskip("obspy") +from specmod import config as cfg # noqa: E402 from specmod.core import Spectrum, SpectrumPair, SpectrumSet # noqa: E402 from specmod.core.collection import log_bin # noqa: E402 from specmod.pipeline import ( # noqa: E402 @@ -54,10 +55,24 @@ def _summary(a: Any) -> dict[str, Any]: } +#: The settings the committed references were captured under; shared with +#: `test_golden_reference.py`. +REFERENCE_CONFIG = Path(__file__).parent / "golden" / "reference.toml" + + @functools.cache def _direct(estimator: str, windows: Any) -> SpectrumSet: + """Build a set under the configuration the committed references hold. + + Pinned rather than ambient because what these tests assert — the recorded + displacement, the Parseval contract, the band a noise floor selects — is + the pipeline's behaviour, not whichever estimator and smoother happen to + be the shipped defaults this month. `tests/golden/reference.toml` is the + same pin `test_golden_reference.py` uses; see REFACTOR_PLAN §6.6. + """ signal, noise = windows() - return spectrum_set_from_streams(signal, noise, estimator=estimator) + with cfg.using(REFERENCE_CONFIG): + return spectrum_set_from_streams(signal, noise, estimator=estimator) # ------------------------------------------------------------ the unit itself diff --git a/tests/test_staged.py b/tests/test_staged.py index ef6965f..96b02f2 100644 --- a/tests/test_staged.py +++ b/tests/test_staged.py @@ -19,6 +19,9 @@ obspy = pytest.importorskip("obspy") +from pathlib import Path # noqa: E402 + +from specmod import config as cfg # noqa: E402 from specmod.config import load_config # noqa: E402 from specmod.distance import ( # noqa: E402 Epicentral, @@ -37,11 +40,17 @@ get_weight_model, ) +#: These tests are about selection patterns and minimiser comparison, not +#: about which estimator ships as the default, so they pin the settings the +#: committed references hold and stay on one spectrum whatever the default +#: becomes. +REFERENCE_CONFIG = Path(__file__).parent / "golden" / "reference.toml" + @functools.cache def _spectra(windows: Any) -> Any: signal, noise = windows() - with contextlib.redirect_stdout(io.StringIO()): + with contextlib.redirect_stdout(io.StringIO()), cfg.using(REFERENCE_CONFIG): return spectrum_set_from_streams(signal, noise) From 33975ed1042e7374e7750c46aa879e666cf0296a Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 26 Sep 2026 18:53:31 +0000 Subject: [PATCH 3/3] docs(fitting): say what the fit reads, and document it as a modelling choice MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The smoothed spectra set the signal-to-noise ratio and the band. With `[fitting] fit_bins = false` — the default in 0.1.1 and since — the fit and the initial guess both read the *unsmoothed* amplitudes inside that band, so the smoother reaches `Omega` only through the band. 1081364's text said otherwise: that Konno-Ohmachi removes the FFT's variance from the fit, and that `fit_bins` is inert under it. Neither holds. Under Konno-Ohmachi `fit_bins` selects the same frequencies and switches the amplitudes, which differ by a median 0.10 dex on the first PNR channel. The comparison behind the default is re-measured on the path the pipeline and the tutorial actually use — velocity, from which `llpsp` is the displacement plateau — against multitaper + `log_bins`, stage one, 28 channels: fft + log_bins median dMw -0.079 per-station -0.71 to +0.23 band 1.55 dec fft + konno_ohmachi median dMw -0.068 per-station -0.25 to +0.19 band 1.68 dec ... with fit_bins=true median dMw -0.073 per-station -0.25 to +0.15 The two default-path FFT rows fit identical amplitudes, so their difference is the band alone. The earlier table was measured on displacement (-0.85/+0.70 and -0.56/+0.52 there): same ordering, wider spread. The tutorial gains "What the fit reads", which states the default and fits the event both ways, computed on every build: on PNR 93% of the unsmoothed points sit in the band's top decade, which is what `weight_method = "log"` exists to shift. It also stops saying the shipped default is multitaper, and says the comparison spectra are smoothed rather than binned. Edited in its builder, as `test_docs_are_current` requires. `fit_bins` now has a comment on the field itself, the corrected upgrading, roadmap and §8.1 text carries the velocity numbers, and a test pins what the setting does under Konno-Ohmachi: same frequencies, different amplitudes, the default reading the unsmoothed ones. Co-Authored-By: Claude Opus 5.5 --- docs/REFACTOR_PLAN.md | 48 +++++++---- docs/_builders/specmod_tutorial.py | 78 +++++++++++++++-- docs/roadmap.md | 24 +++--- docs/upgrading.md | 37 ++++---- src/specmod/config/sections.py | 21 +++-- tests/test_fitting_defaults.py | 33 +++++-- tests/test_magnitude.py | 10 +-- tutorial/specmod-tutorial.ipynb | 134 ++++++++++++++++++++++------- 8 files changed, 277 insertions(+), 108 deletions(-) diff --git a/docs/REFACTOR_PLAN.md b/docs/REFACTOR_PLAN.md index 627b890..640e510 100644 --- a/docs/REFACTOR_PLAN.md +++ b/docs/REFACTOR_PLAN.md @@ -3722,21 +3722,31 @@ changes anything shipped. ### Still open 1. ~~**Default estimator.**~~ **Resolved: FFT + Konno-Ohmachi**, from 0.6. The - measurement below is what decided it, on the 28 PNR windows with each - station fitted freely: - - | Against the shipped multitaper + `log_bins` | Median ΔMw | Per-station range | `fc` ratio | - |---|---|---|---| - | `fft` + `log_bins` | −0.039 | −0.85 to +1.90 | 0.003 to 5.0 | - | `fft` + `konno_ohmachi` | +0.027 | −0.52 to +0.51 | 0.038 to 3.2 | - - **The smoothing is what makes the switch defensible, not the estimator.** - FFT with the shipped binner moves a station by up to 1.9 magnitude units, - because `log_bins` reduces the axis without reducing the variance — 151 log - bins over a few hundred Fourier samples is roughly one sample per bin. - Pairing it with a real smoother halves the worst excursion and pulls the - `fc` range in by a third. The median station barely moves either way, so the - per-station spread is the whole of the argument. + measurement below is what decided it, on the 28 PNR channels, each fitted + independently (stage one) on velocity, which is the motion the pipeline fits + by default and from which `llpsp` is the displacement plateau: + + | Against the shipped multitaper + `log_bins` | Median ΔMw | Per-station ΔMw | `fc` ratio, median (range) | Median band (decades) | Median points fitted | + |---|---|---|---|---|---| + | `multitaper` + `log_bins` | — | — | — | 1.73 | 113 | + | `fft` + `log_bins` | −0.079 | −0.71 to +0.23 | 0.91 (0.18 to 3.9) | 1.55 | 107 | + | `fft` + `konno_ohmachi` | −0.068 | −0.25 to +0.19 | 0.82 (0.17 to 3.5) | 1.68 | 120 | + | `fft` + `konno_ohmachi`, `fit_bins = true` | −0.073 | −0.25 to +0.15 | 0.90 (0.18 to 2.2) | 1.68 | 120 | + + **What each setting reaches.** The estimator sets the spectrum that is + fitted. The smoother sets the band: the signal-to-noise ratio and the gate + are computed on smoothed spectra. With `[fitting] fit_bins = false`, the + default in 0.1.1 and since, the fit and the initial guess both read the + *unsmoothed* amplitudes inside that band, so the smoother reaches `Omega` + only through the band. The two default-path FFT rows therefore fit identical + amplitudes, and the whole difference between them — the per-station range + narrowing from −0.71/+0.23 to −0.25/+0.19 — comes from Konno–Ohmachi + selecting a wider band. Fitting the smoothed spectrum instead + (`fit_bins = true`) changes little further on this event. + + Measured on displacement instead, the same comparison gives −0.85/+0.70 for + `fft` + `log_bins` and −0.56/+0.52 for `fft` + `konno_ohmachi`: the same + ordering, a wider spread. What this event cannot settle, and what the decision therefore does not rest on: the *event* corner is unconstrained under all three (stage-one @@ -3757,10 +3767,10 @@ changes anything shipped. `tests/golden/reference.toml` close that, so the references now hold their own settings and a default change leaves every committed number untouched. - **A coupling the measurement did not predict.** `[fitting] fit_bins` selects - the binned spectrum, which under an axis-preserving smoother is the Fourier - axis — so with Konno-Ohmachi the setting fits the same 98 points either way - and is inert. It only bites under `log_bins`. + **The fitting choice is documented where it is made.** Whether the fit + reads the unsmoothed or the smoothed spectrum is a modelling decision a user + should make knowingly, so the tutorial's *What the fit reads* states the + default and fits the event both ways. 2. **Python floor.** 3.11 is proposed. Any users stuck on 3.9/3.10? 3. **History rewrite.** Deleting the 9.9 MB catalog and stripping the notebook diff --git a/docs/_builders/specmod_tutorial.py b/docs/_builders/specmod_tutorial.py index c89c774..869ec63 100644 --- a/docs/_builders/specmod_tutorial.py +++ b/docs/_builders/specmod_tutorial.py @@ -237,12 +237,14 @@ ``` `spectrum_set_from_streams` transforms both windows, puts the noise on the -signal's frequency axis, bins both, raises the noise to account for what sits -*under* the signal, takes the ratio and selects the band where it passes. +signal's frequency axis, smooths both, raises the noise to account for what +sits *under* the signal, takes the ratio and selects the band where it passes. +The smoothing is Konno–Ohmachi by default (`[smoothing] method`). It is used +for the ratio and the band, not for the fit — see *What the fit reads* below. Which estimator is used is configuration, not code: `fft`, `welch`, `multitaper`, `quadratic` and `cwt` all satisfy the same Parseval contract, so -they are interchangeable here. The shipped default is `multitaper`. +they are interchangeable here. The shipped default is `fft`. """) @@ -256,7 +258,7 @@ code(""" -# One station in detail: signal, noise, the binned spectra the ratio is +# One station in detail: signal, noise, the smoothed spectra the ratio is # actually computed on, the selected band (red) and the resolution floor (grey). from specmod.plotting import plot_pair, plot_set @@ -453,9 +455,10 @@ def band_report(spectra: SpectrumSet) -> "pd.DataFrame": print(f"velocity {spectra[station].signal.unit}") print(f"displacement {displacement[station].signal.unit}") -# The unbinned signal-to-noise ratio is invariant under this — both spectra are -# divided by the same 2*pi*f — but the *binned* ratio is not, because a bin -# holds a geometric mean. So a few bands do move. +# The unsmoothed signal-to-noise ratio is invariant under this — both spectra +# are divided by the same 2*pi*f — but the smoothed ratio is not: Konno-Ohmachi +# averages amplitudes across a window over which 2*pi*f changes, so some +# bands move. moved = [i for i in spectra.ids() if spectra[i].band != displacement[i].band] print(f"bands that moved: {len(moved)} of {len(spectra)}") """) @@ -473,9 +476,10 @@ def band_report(spectra: SpectrumSet) -> "pd.DataFrame": maximum falls. `FitSpectra(spectra)` then `fit_spectra()` is the whole thing. The minimiser -(Powell), the `t*` lower bound and whether to fit the binned or unbinned +(Powell), the `t*` lower bound and whether to fit the smoothed or unsmoothed spectrum all come from `[fitting]` in the configuration, so the defaults are -recorded rather than remembered. +recorded rather than remembered. The last of those is a modelling choice with +consequences; *What the fit reads* below compares the two. """) @@ -551,6 +555,62 @@ def band_report(spectra: SpectrumSet) -> "pd.DataFrame": """) +md(""" +### What the fit reads + +The smoothed spectra choose the band — the signal-to-noise ratio and the gate +are computed on them — but **the fit reads the unsmoothed spectrum inside that +band**. Smoothing decides *where* the model is fitted, not *what* it is fitted +to. That is `[fitting] fit_bins = false`, the default, and it is a modelling +choice: + +- **Unsmoothed** (the default). Every Fourier sample in the band, each with its + full scatter. The samples are evenly spaced in Hz, so most of them sit in the + top decade of the band and dominate an unweighted misfit; + `[fitting] weight_method = "log"` weights residuals by `1/f` to shift that. +- **Smoothed** (`fit_bins = true`). The spectrum the band was chosen on. Under + Konno–Ohmachi that is the same frequencies with the scatter averaged down, + but neighbouring points share samples, so they are not independent and the + misfit counts the same information more than once. Under `log_bins` it is + one value per occupied bin, evenly spaced in log frequency. + +The cell below fits the event both ways. +""") + + +code(""" +smoothed_fits = FitSpectra(spectra, fit_bins=True) +smoothed_fits.fit_spectra() + +first = spectra.ids()[0] +raw_fit, smoothed_fit = fits.models[first], smoothed_fits.models[first] +low, high = raw_fit.mod_freq[0], raw_fit.mod_freq[-1] +top_decade = np.mean(raw_fit.mod_freq >= high / 10) +print( + f"{first}: {raw_fit.mod_freq.size} points unsmoothed, " + f"{smoothed_fit.mod_freq.size} smoothed, over {low:.2f}-{high:.2f} Hz" +) +print(f" share of the points in the band's top decade: {top_decade:.0%}") +print( + " median |difference| in the log10 amplitudes the two fits see: " + f"{np.median(np.abs(raw_fit.mod_amp - smoothed_fit.mod_amp)):.3f}" +) + +raw_table = fits.table.set_index("id") +smoothed_table = smoothed_fits.table.set_index("id") +d_plateau = (smoothed_table["llpsp"] - raw_table["llpsp"]).abs() +fc_ratio = smoothed_table["fc"] / raw_table["fc"] +print(f"\\nacross {len(raw_table)} channels, smoothed against unsmoothed:") +print( + f" |d log10 plateau| median {d_plateau.median():.3f}, max {d_plateau.max():.3f}" +) +print( + f" fc ratio median {fc_ratio.median():.2f}, " + f"range {fc_ratio.min():.2f} to {fc_ratio.max():.2f}" +) +""") + + md(""" ### The fit is not unique, and that is not a detail diff --git a/docs/roadmap.md b/docs/roadmap.md index 1e67407..fe53d84 100644 --- a/docs/roadmap.md +++ b/docs/roadmap.md @@ -198,25 +198,23 @@ Merged on `main`, and it will name its version when a release goes out. ### The default is FFT with Konno-Ohmachi smoothing The shipped pair was multitaper with log-binning. Measured on the 28 PNR -windows, with each station fitted freely, against that pair: +channels, each fitted independently on velocity, against that pair: -| | Median ΔMw | Per-station range | `fc` ratio | +| | Median ΔMw | Per-station ΔMw | `fc` ratio, median (range) | |---|---|---|---| -| `fft` + `log_bins` | −0.039 | −0.85 to +1.90 | 0.003 to 5.0 | -| `fft` + `konno_ohmachi` | +0.027 | −0.52 to +0.51 | 0.038 to 3.2 | +| `fft` + `log_bins` | −0.079 | −0.71 to +0.23 | 0.91 (0.18 to 3.9) | +| `fft` + `konno_ohmachi` | −0.068 | −0.25 to +0.19 | 0.82 (0.17 to 3.5) | -The smoother is what earns the change. `log_bins` reduces the frequency axis -without reducing the variance — 151 bins over a few hundred Fourier samples is -about one sample per bin — so an FFT read through it swings a station by up to -1.9 magnitude units. A real smoother halves the worst excursion and pulls the -corner-frequency range in by a third. The median station barely moves either -way, which is why the per-station spread is the argument and the ensemble -number is not: this event's corner is unconstrained under all three. +By default the fit reads the unsmoothed spectrum inside the band the smoothed +signal-to-noise ratio selects, so the two FFT rows fit identical amplitudes +and differ only in the band. Konno–Ohmachi's is wider, and that is what +narrows the per-station spread. This event's corner is unconstrained under +every pair, so the event-level number is not evidence either way. **Every result moves** unless a configuration already pinned both settings. [Upgrading](upgrading.md#from-05-to-06) has what to pin to keep the old -numbers, and the one coupling this exposes: `[fitting] fit_bins` is inert -under a smoother that preserves the frequency axis. +numbers, and what `[fitting] fit_bins = true` changes: it fits the smoothed +spectrum instead of the unsmoothed one. ### A configuration can be held across an ambient read diff --git a/docs/upgrading.md b/docs/upgrading.md index 181ed33..f9a843c 100644 --- a/docs/upgrading.md +++ b/docs/upgrading.md @@ -16,18 +16,24 @@ See [Releasing the software](releasing.md#what-a-version-number-means-while-this They were `multitaper` and `log_bins`. Nothing is renamed and no signature moves; **every result changes** unless your configuration already pinned both. -On the 28 PNR windows, against the old pair: +The two do different jobs. The estimator sets the spectrum that is fitted. +The smoother sets the band the fit is confined to — the signal-to-noise ratio +and the gate are computed on the smoothed spectra — but **by default the fit +reads the unsmoothed spectrum inside that band** (`[fitting] fit_bins = false`), +so the smoother reaches the result only through the band. -| | Median ΔMw | Per-station range | `fc` ratio | -|---|---|---|---| -| `fft` + `log_bins` | −0.039 | −0.85 to +1.90 | 0.003 to 5.0 | -| `fft` + `konno_ohmachi` | +0.027 | −0.52 to +0.51 | 0.038 to 3.2 | +On the 28 PNR channels, each fitted independently on velocity, against the old +pair: -The smoother is what earns the change, not the transform. `log_bins` reduces -the frequency axis without reducing the variance — 151 bins over a few hundred -Fourier samples is about one sample per bin — so an FFT paired with it swings a -station by up to 1.9 magnitude units. A real smoother halves the worst -excursion. The median station barely moves either way. +| | Median ΔMw | Per-station ΔMw | `fc` ratio, median (range) | Median band (decades) | +|---|---|---|---|---| +| `multitaper` + `log_bins` (0.5) | — | — | — | 1.73 | +| `fft` + `log_bins` | −0.079 | −0.71 to +0.23 | 0.91 (0.18 to 3.9) | 1.55 | +| `fft` + `konno_ohmachi` | −0.068 | −0.25 to +0.19 | 0.82 (0.17 to 3.5) | 1.68 | + +The two FFT rows fit identical amplitudes; they differ only in the band. +Konno–Ohmachi gives a wider one — a median 120 points fitted per channel +against 107 — and that is what narrows the per-station spread. **To keep the old numbers**, pin both: @@ -39,12 +45,13 @@ estimator = "multitaper" method = "log_bins" ``` -Two consequences worth knowing: +Two things worth knowing: -- **`[fitting] fit_bins` is inert under the new default.** It fits the binned - spectrum, and Konno–Ohmachi preserves the frequency axis, so the "binned" - axis *is* the Fourier axis and the setting selects the same points either - way. It only bites under `log_bins`, or any other axis-reducing smoother. +- **`[fitting] fit_bins = true` fits the smoothed spectrum instead.** Under + Konno–Ohmachi that is the same frequencies with the scatter averaged down; + on PNR the per-station range becomes −0.25 to +0.15. Neighbouring smoothed + points share samples, so they are not independent. The tutorial's *What the + fit reads* compares the two. - **Published results should pin the version anyway**, but this is the release that shows why: a re-run under 0.6 of an analysis done under 0.5 will not reproduce it, and nothing in the call changed to warn you. diff --git a/src/specmod/config/sections.py b/src/specmod/config/sections.py index 260768d..a2fcbf2 100644 --- a/src/specmod/config/sections.py +++ b/src/specmod/config/sections.py @@ -99,11 +99,10 @@ class WindowsConfig: class TransformConfig: """Time-to-frequency conversion. Consumed by :mod:`specmod.transforms`.""" - #: FFT since 0.6, paired with the Konno-Ohmachi default in ``[smoothing]``. - #: The pairing is the point: an FFT read through a log-binner carries the - #: periodogram's variance into the fit, and it is the smoother that removes - #: it. ``docs/upgrading.md`` has the measured comparison, and - #: REFACTOR_PLAN §8.1 the decision. + #: FFT since 0.6, paired with Konno-Ohmachi in ``[smoothing]``. This is the + #: spectrum the fit reads; the smoother only chooses the band it is read + #: over, unless ``[fitting] fit_bins`` is on. ``docs/upgrading.md`` has the + #: measured comparison, and REFACTOR_PLAN §8.1 the decision. estimator: Literal[ "multitaper", "fft", "welch", "cwt", "prieto", "quadratic", "mtspec" ] = "fft" @@ -152,9 +151,11 @@ class SmoothingConfig: seismic spectrum and ``docs/choosing-a-transform.md`` for the measured comparison. - That distinction reaches ``[fitting] fit_bins``, which fits the binned - spectrum: under an axis-preserving method there is no separate binned axis, - so the setting selects the same points either way. + The smoothed spectra set the signal-to-noise ratio and the band. The fit + reads the unsmoothed spectrum inside that band unless ``[fitting] + fit_bins`` is on, in which case it reads the smoothed one: one value per + bin under ``log_bins``, the smoothed amplitude at every frequency under the + others. """ #: Resolved through `specmod.smoothing.SMOOTHERS`. Only the parameters @@ -264,6 +265,10 @@ class FittingConfig: """Minimisation.""" method: str = "powell" + #: Fit the smoothed spectrum rather than the unsmoothed one. Off by default: + #: the smoothed spectra choose the band, and the fit reads the unsmoothed + #: amplitudes inside it. A modelling choice, compared in the tutorial's + #: "What the fit reads". fit_bins: bool = False weight_method: Literal["none", "log"] = "none" diff --git a/tests/test_fitting_defaults.py b/tests/test_fitting_defaults.py index bd5f353..c23c406 100644 --- a/tests/test_fitting_defaults.py +++ b/tests/test_fitting_defaults.py @@ -171,12 +171,8 @@ def test_a_parameter_resting_on_its_bound_fails_the_fit(self) -> None: def test_fit_bins_comes_from_configuration(self, pnr_windows: Any) -> None: """Fitting the binned spectrum uses far fewer points than the raw one. - Only where the smoother replaces the frequency axis, which is why this - runs under the pinned reference config rather than the shipped - defaults. Under an axis-preserving smoother — Konno-Ohmachi, which is - the default from 0.6 — the comparison axis *is* the Fourier axis, so - ``fit_bins`` selects the same 98 points either way and the setting is - inert. See ``docs/upgrading.md``. + Under `log_bins`, which the pinned reference config uses: one value per + occupied bin against every Fourier sample in the band. """ spectra = _spectra(pnr_windows) raw = _fit(spectra, fit_bins=False) @@ -188,6 +184,31 @@ def test_fit_bins_comes_from_configuration(self, pnr_windows: Any) -> None: assert load_config().config.fitting.fit_bins is False assert _fit(spectra).models[id].mod_freq.size == raw.models[id].mod_freq.size + def test_under_konno_ohmachi_it_switches_the_amplitudes_not_the_points( + self, pnr_windows: Any + ) -> None: + """The shipped default smoother keeps the frequency axis. + + So ``fit_bins`` selects the same frequencies either way and changes the + amplitudes read at them: unsmoothed by default, smoothed when on. The + default reading the unsmoothed ones is the modelling choice the + tutorial's "What the fit reads" documents. + """ + signal, noise = pnr_windows() + with cfg.using(smoothing={"method": "konno_ohmachi"}): + spectra = spectrum_set_from_streams(signal, noise, estimator="fft") + unsmoothed = _fit(spectra, fit_bins=False) + smoothed = _fit(spectra, fit_bins=True) + + id = next(iter(unsmoothed.models)) + raw, smooth = unsmoothed.models[id], smoothed.models[id] + assert np.array_equal(raw.mod_freq, smooth.mod_freq) + assert not np.allclose(raw.mod_amp, smooth.mod_amp) + # And the amplitudes the default fit reads are the unsmoothed ones. + expected = np.log10(spectra[id].signal.amp) + inside = np.isin(spectra[id].signal.freq, raw.mod_freq) + np.testing.assert_array_equal(raw.mod_amp, expected[inside]) + class TestInitialGuess: def test_it_asks_the_model_which_parameters_it_takes( diff --git a/tests/test_magnitude.py b/tests/test_magnitude.py index 4d1d6d8..5ca7fc1 100644 --- a/tests/test_magnitude.py +++ b/tests/test_magnitude.py @@ -304,11 +304,11 @@ def test_the_shipped_default_moves_omega_by_a_constant(pnr_windows) -> None: # so without this the default's own answer would be asserted nowhere and the change of default would be invisible to the suite. - The shift is **-0.1615 Mw, identical under both constant sets**, which is - the shape a change of estimator and smoother should have: it moves the - low-frequency plateau the fit reads Omega off, and the medium constants - are an offset applied afterwards. A number that moved differently under - the two would mean something had changed in the fitting instead. + The shift is **-0.1615 Mw, identical under both constant sets**, as it + must be: the medium constants add the same offset to every station's + ``log10 M0``, so only the fitted plateaus move. They move because the + estimator changes the amplitudes the fit reads and the smoother changes + the band it reads them over. Against the catalogue's 2.9 the two sets move in opposite directions — the default constants from 0.194 to 0.032 away, the plan's from 0.160 to diff --git a/tutorial/specmod-tutorial.ipynb b/tutorial/specmod-tutorial.ipynb index 541254f..afbffab 100644 --- a/tutorial/specmod-tutorial.ipynb +++ b/tutorial/specmod-tutorial.ipynb @@ -261,12 +261,14 @@ "```\n", "\n", "`spectrum_set_from_streams` transforms both windows, puts the noise on the\n", - "signal's frequency axis, bins both, raises the noise to account for what sits\n", - "*under* the signal, takes the ratio and selects the band where it passes.\n", + "signal's frequency axis, smooths both, raises the noise to account for what\n", + "sits *under* the signal, takes the ratio and selects the band where it passes.\n", + "The smoothing is Konno–Ohmachi by default (`[smoothing] method`). It is used\n", + "for the ratio and the band, not for the fit — see *What the fit reads* below.\n", "\n", "Which estimator is used is configuration, not code: `fft`, `welch`,\n", "`multitaper`, `quadratic` and `cwt` all satisfy the same Parseval contract, so\n", - "they are interchangeable here. The shipped default is `multitaper`." + "they are interchangeable here. The shipped default is `fft`." ] }, { @@ -290,7 +292,7 @@ "metadata": {}, "outputs": [], "source": [ - "# One station in detail: signal, noise, the binned spectra the ratio is\n", + "# One station in detail: signal, noise, the smoothed spectra the ratio is\n", "# actually computed on, the selected band (red) and the resolution floor (grey).\n", "from specmod.plotting import plot_pair, plot_set\n", "\n", @@ -532,9 +534,10 @@ "print(f\"velocity {spectra[station].signal.unit}\")\n", "print(f\"displacement {displacement[station].signal.unit}\")\n", "\n", - "# The unbinned signal-to-noise ratio is invariant under this — both spectra are\n", - "# divided by the same 2*pi*f — but the *binned* ratio is not, because a bin\n", - "# holds a geometric mean. So a few bands do move.\n", + "# The unsmoothed signal-to-noise ratio is invariant under this — both spectra\n", + "# are divided by the same 2*pi*f — but the smoothed ratio is not: Konno-Ohmachi\n", + "# averages amplitudes across a window over which 2*pi*f changes, so some\n", + "# bands move.\n", "moved = [i for i in spectra.ids() if spectra[i].band != displacement[i].band]\n", "print(f\"bands that moved: {len(moved)} of {len(spectra)}\")" ] @@ -558,9 +561,10 @@ "maximum falls.\n", "\n", "`FitSpectra(spectra)` then `fit_spectra()` is the whole thing. The minimiser\n", - "(Powell), the `t*` lower bound and whether to fit the binned or unbinned\n", + "(Powell), the `t*` lower bound and whether to fit the smoothed or unsmoothed\n", "spectrum all come from `[fitting]` in the configuration, so the defaults are\n", - "recorded rather than remembered." + "recorded rather than remembered. The last of those is a modelling choice with\n", + "consequences; *What the fit reads* below compares the two." ] }, { @@ -672,6 +676,70 @@ "cell_type": "markdown", "id": "39", "metadata": {}, + "source": [ + "### What the fit reads\n", + "\n", + "The smoothed spectra choose the band — the signal-to-noise ratio and the gate\n", + "are computed on them — but **the fit reads the unsmoothed spectrum inside that\n", + "band**. Smoothing decides *where* the model is fitted, not *what* it is fitted\n", + "to. That is `[fitting] fit_bins = false`, the default, and it is a modelling\n", + "choice:\n", + "\n", + "- **Unsmoothed** (the default). Every Fourier sample in the band, each with its\n", + " full scatter. The samples are evenly spaced in Hz, so most of them sit in the\n", + " top decade of the band and dominate an unweighted misfit;\n", + " `[fitting] weight_method = \"log\"` weights residuals by `1/f` to shift that.\n", + "- **Smoothed** (`fit_bins = true`). The spectrum the band was chosen on. Under\n", + " Konno–Ohmachi that is the same frequencies with the scatter averaged down,\n", + " but neighbouring points share samples, so they are not independent and the\n", + " misfit counts the same information more than once. Under `log_bins` it is\n", + " one value per occupied bin, evenly spaced in log frequency.\n", + "\n", + "The cell below fits the event both ways." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "40", + "metadata": {}, + "outputs": [], + "source": [ + "smoothed_fits = FitSpectra(spectra, fit_bins=True)\n", + "smoothed_fits.fit_spectra()\n", + "\n", + "first = spectra.ids()[0]\n", + "raw_fit, smoothed_fit = fits.models[first], smoothed_fits.models[first]\n", + "low, high = raw_fit.mod_freq[0], raw_fit.mod_freq[-1]\n", + "top_decade = np.mean(raw_fit.mod_freq >= high / 10)\n", + "print(\n", + " f\"{first}: {raw_fit.mod_freq.size} points unsmoothed, \"\n", + " f\"{smoothed_fit.mod_freq.size} smoothed, over {low:.2f}-{high:.2f} Hz\"\n", + ")\n", + "print(f\" share of the points in the band's top decade: {top_decade:.0%}\")\n", + "print(\n", + " \" median |difference| in the log10 amplitudes the two fits see: \"\n", + " f\"{np.median(np.abs(raw_fit.mod_amp - smoothed_fit.mod_amp)):.3f}\"\n", + ")\n", + "\n", + "raw_table = fits.table.set_index(\"id\")\n", + "smoothed_table = smoothed_fits.table.set_index(\"id\")\n", + "d_plateau = (smoothed_table[\"llpsp\"] - raw_table[\"llpsp\"]).abs()\n", + "fc_ratio = smoothed_table[\"fc\"] / raw_table[\"fc\"]\n", + "print(f\"\\nacross {len(raw_table)} channels, smoothed against unsmoothed:\")\n", + "print(\n", + " f\" |d log10 plateau| median {d_plateau.median():.3f}, max {d_plateau.max():.3f}\"\n", + ")\n", + "print(\n", + " f\" fc ratio median {fc_ratio.median():.2f}, \"\n", + " f\"range {fc_ratio.min():.2f} to {fc_ratio.max():.2f}\"\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "41", + "metadata": {}, "source": [ "### The fit is not unique, and that is not a detail\n", "\n", @@ -694,7 +762,7 @@ { "cell_type": "code", "execution_count": null, - "id": "40", + "id": "42", "metadata": {}, "outputs": [], "source": [ @@ -722,7 +790,7 @@ }, { "cell_type": "markdown", - "id": "41", + "id": "43", "metadata": {}, "source": [ "Most stations agree to a fraction of a percent. A few do not, and the top of\n", @@ -738,7 +806,7 @@ { "cell_type": "code", "execution_count": null, - "id": "42", + "id": "44", "metadata": {}, "outputs": [], "source": [ @@ -758,7 +826,7 @@ { "cell_type": "code", "execution_count": null, - "id": "43", + "id": "45", "metadata": {}, "outputs": [], "source": [ @@ -774,7 +842,7 @@ }, { "cell_type": "markdown", - "id": "44", + "id": "46", "metadata": {}, "source": [ "Look at where they differ: high up the falling limb, where the source corner\n", @@ -799,7 +867,7 @@ { "cell_type": "code", "execution_count": null, - "id": "45", + "id": "47", "metadata": {}, "outputs": [], "source": [ @@ -813,7 +881,7 @@ }, { "cell_type": "markdown", - "id": "46", + "id": "48", "metadata": {}, "source": [ "### Why the answer comes from many stations, and from two stages\n", @@ -838,7 +906,7 @@ { "cell_type": "code", "execution_count": null, - "id": "47", + "id": "49", "metadata": {}, "outputs": [], "source": [ @@ -883,7 +951,7 @@ }, { "cell_type": "markdown", - "id": "48", + "id": "50", "metadata": {}, "source": [ "A factor of three in stress drop at the worst station becomes **under 2%**\n", @@ -898,7 +966,7 @@ { "cell_type": "code", "execution_count": null, - "id": "49", + "id": "51", "metadata": {}, "outputs": [], "source": [ @@ -925,7 +993,7 @@ }, { "cell_type": "markdown", - "id": "50", + "id": "52", "metadata": {}, "source": [ "0.3% in $t^*$ and about 0.002 in $\\log_{10}\\Omega$ — which is 0.003 magnitude\n", @@ -947,7 +1015,7 @@ }, { "cell_type": "markdown", - "id": "51", + "id": "53", "metadata": {}, "source": [ "### The same thing, as one call\n", @@ -966,7 +1034,7 @@ { "cell_type": "code", "execution_count": null, - "id": "52", + "id": "54", "metadata": {}, "outputs": [], "source": [ @@ -982,7 +1050,7 @@ }, { "cell_type": "markdown", - "id": "53", + "id": "55", "metadata": {}, "source": [ "`describe()` prints the spread as well as the mean, and that is deliberate. A\n", @@ -1011,7 +1079,7 @@ { "cell_type": "code", "execution_count": null, - "id": "54", + "id": "56", "metadata": {}, "outputs": [], "source": [ @@ -1033,7 +1101,7 @@ { "cell_type": "code", "execution_count": null, - "id": "55", + "id": "57", "metadata": {}, "outputs": [], "source": [ @@ -1061,7 +1129,7 @@ }, { "cell_type": "markdown", - "id": "56", + "id": "58", "metadata": {}, "source": [ "One quality-control decision, a factor of 1.5 in stress drop. That is not an\n", @@ -1082,7 +1150,7 @@ { "cell_type": "code", "execution_count": null, - "id": "57", + "id": "59", "metadata": {}, "outputs": [], "source": [ @@ -1104,7 +1172,7 @@ }, { "cell_type": "markdown", - "id": "58", + "id": "60", "metadata": {}, "source": [ "So changing the minimiser changes *which stations vote*, not only how each one\n", @@ -1120,7 +1188,7 @@ }, { "cell_type": "markdown", - "id": "59", + "id": "61", "metadata": {}, "source": [ "## 5. Save the results" @@ -1128,7 +1196,7 @@ }, { "cell_type": "markdown", - "id": "60", + "id": "62", "metadata": {}, "source": [ "Two formats, because the data is used two different ways.\n", @@ -1147,7 +1215,7 @@ { "cell_type": "code", "execution_count": null, - "id": "61", + "id": "63", "metadata": {}, "outputs": [], "source": [ @@ -1166,7 +1234,7 @@ { "cell_type": "code", "execution_count": null, - "id": "62", + "id": "64", "metadata": {}, "outputs": [], "source": [ @@ -1177,7 +1245,7 @@ }, { "cell_type": "markdown", - "id": "63", + "id": "65", "metadata": {}, "source": [ "---\n",