diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml new file mode 100644 index 0000000..cfb57c7 --- /dev/null +++ b/.github/workflows/docs.yml @@ -0,0 +1,57 @@ +name: docs + +on: + push: + branches: [main] + pull_request: + workflow_dispatch: + +# Pages needs these; the default token is read-only for contents. +permissions: + contents: read + pages: write + id-token: write + +# A second push while a deploy is in flight should win, not queue behind it. +concurrency: + group: pages + cancel-in-progress: false + +jobs: + build: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + with: + fetch-depth: 0 # hatch-vcs derives the version from tags + - uses: astral-sh/setup-uv@v5 + - run: uv python install 3.11 + - run: uv venv # create a venv to install into + + # The package itself, not just the docs extra: autodoc imports every + # module it documents, so a docs build is also an import check. + - run: uv pip install -e ".[docs,io]" + + - name: Build + # No -W. Intersphinx resolves seven inventories over the network and + # warns when one is briefly unreachable, which would turn a third + # party's downtime into a red build. A genuinely broken build exits + # non-zero on its own. + run: uv run sphinx-build -b html docs docs/_build/html + + - uses: actions/upload-pages-artifact@v3 + with: + path: docs/_build/html + + deploy: + # Only main deploys. A pull request builds, which is what catches a broken + # reference before it lands, and stops there. + if: github.ref == 'refs/heads/main' + needs: build + runs-on: ubuntu-latest + environment: + name: github-pages + url: ${{ steps.deployment.outputs.page_url }} + steps: + - id: deployment + uses: actions/deploy-pages@v4 \ No newline at end of file diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 8fd4464..f6fc7ca 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,3 +1,9 @@ +# `pre-commit install` wires only the pre-commit stage by default, which left +# the commit-msg hook below inert in every clone that ran the documented setup +# — including the one that then committed three session links. Naming the +# stages here means one `pre-commit install` installs both. +default_install_hook_types: [pre-commit, commit-msg] + repos: # CI runs `uvx ruff` unpinned, so it always gets the newest release. Keep # this rev at that newest release or CI will fail on rules this does not have. diff --git a/.readthedocs.yaml b/.readthedocs.yaml new file mode 100644 index 0000000..332a3ab --- /dev/null +++ b/.readthedocs.yaml @@ -0,0 +1,47 @@ +# Read the Docs build configuration. +# +# Read the Docs publishes the documentation; GitHub Actions only checks that it +# builds. The reason for the split is versions: this package is alpha and users +# are told to pin an exact version, so the documentation for a release has to +# stay readable after the trunk has moved on. Read the Docs keeps one build per +# tag, a `stable` pointing at the newest, and a `latest` from `main`, with a +# version switcher between them — and it builds pull requests to their own +# throwaway URL. +# +# Schema: https://docs.readthedocs.com/platform/stable/config-file/v2.html + +version: 2 + +build: + os: ubuntu-24.04 + tools: + python: "3.11" # the version CI's lint, typecheck and docs jobs use + jobs: + post_checkout: + # Read the Docs clones shallow and without tags, to save time. hatch-vcs + # derives the version from `git describe`, and `docs/conf.py` reads it + # back through importlib.metadata — so without these the sidebar reads + # the fallback 0.0.0 on every build, including tagged ones. + # `|| true` because a repository that is already complete makes + # --unshallow exit non-zero, which would fail the build. + - git fetch --unshallow || true + - git fetch --tags || true + +sphinx: + configuration: docs/conf.py + # Deliberately off, for the same reason the CI job has no `-W`: intersphinx + # resolves seven inventories over the network and warns whenever one of them + # is briefly unreachable. That would make a third party's downtime a failed + # documentation build. A genuinely broken build still exits non-zero. + fail_on_warning: false + +python: + install: + - method: pip + path: . + extra_requirements: + # `io` alongside `docs` because autodoc imports every module it + # documents, and `specmod.io` imports h5py and pyarrow. Without it the + # API reference loses those pages to import errors. + - docs + - io diff --git a/.release-please-manifest.json b/.release-please-manifest.json new file mode 100644 index 0000000..a915e8c --- /dev/null +++ b/.release-please-manifest.json @@ -0,0 +1,3 @@ +{ + ".": "0.1.1" +} diff --git a/AGENTS.md b/AGENTS.md new file mode 100644 index 0000000..1d21c69 --- /dev/null +++ b/AGENTS.md @@ -0,0 +1,86 @@ +# Working on SpecMod with a coding agent + +Rules for Claude Code, Codex, and any other agent committing to this +repository. `CLAUDE.md` points here so there is one copy. + +The long version of everything below is +[`docs/development.md`](docs/development.md). This file is the part an agent +must not get wrong, and every entry is here because it has actually gone wrong. + +## Before the first commit + +```sh +uv venv && uv pip install -e ".[dev]" +pre-commit install # installs BOTH the pre-commit and commit-msg hooks +``` + +`pre-commit install` is not optional. In a fresh container it is easy to skip, +and the `commit-msg` hook is the only thing enforcing the rule below. + +## Never publish session links + +**No `Claude-Session:` trailers, session URLs, or agent-console links** in +commit messages, PR titles, PR bodies, code comments, or anything else that +lands in the repository. It is public; those links are private state. + +`Co-Authored-By:` is fine. If your harness appends a session trailer by +default, strip it — the repository's rule wins over the harness default. The +`commit-msg` hook rejects it, which is why installing the hooks comes first. + +## Commit messages are load-bearing + +[Conventional Commits](https://www.conventionalcommits.org). They are not a +style preference: `release-please` reads them to compute the version bump and +to generate `CHANGELOG.md`. See +[`docs/releasing.md`](docs/releasing.md). + +- `feat:` minor, `fix:` patch, `refactor:` / `docs:` / `build:` appear in the + changelog, `test:` / `ci:` / `chore:` are hidden. +- `!` or a `BREAKING CHANGE:` footer bumps the minor while the project is + `0.x`, not the major. +- Say *why*, with the measurement if there was one. The history is the record + of what was checked; a message that only restates the diff wastes it. + +## You cannot push workflow files + +A GitHub App token has no `workflows` permission, so any push touching +`.github/workflows/` is rejected outright. Write the intended file to +`ci/workflows/.yml` instead and say in the PR that it needs copying +across. `tools/check_ci_mirror.py` runs in the `lint` job and fails until the +copy is made — that failure is the reminder, not a fault. +See [`ci/README.md`](ci/README.md). + +## Verify before reporting + +Run these, and report what they actually printed: + +```sh +pytest -m "not dataset and not notebook" # the suite CI runs +pytest --without-optional-extras # what a default install sees +ruff check src/ tests/ tools/ && ruff format --check src/ tests/ tools/ +mypy +python tools/check_ci_mirror.py +sphinx-build -b html docs docs/_build/html # if docs/ changed +``` + +`--without-optional-extras` matters: a development environment with +`specmod[multitaper]` installed passes tests that CI fails. + +## Things that look like noise and are not + +- **Golden references.** `tests/golden/*.json` is a record of numbers this code + used to produce. Do not regenerate it to make a test pass. If a change moves + a number, that is the finding — say which number, by how much, and why, and + regenerate deliberately with `python tools/make_golden.py`. +- **Measured tables in the docs.** Numbers in `docs/*.md` are generated between + markers by `python tools/measure_docs.py`. Edit the tool, not the table. +- **Tolerances.** Several carry a comment explaining what was measured to + choose them. Widening one to get to green, without measuring, is the specific + failure `docs/REFACTOR_PLAN.md` §6.6 exists to catch. + +## Say what you did not check + +The plan's §6.6 is an audit of claims in this repository that turned out to +describe mechanisms nobody had built. Do not add to it. If something is +untested, unreproducible, or assumed, write that down next to the claim — a +bound with a number behind it beats a confident sentence. diff --git a/CLAUDE.md b/CLAUDE.md new file mode 100644 index 0000000..a3af6e2 --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1,5 @@ +# CLAUDE.md + +See [`AGENTS.md`](AGENTS.md) — one copy of the rules, for every agent. + +The long version is [`docs/development.md`](docs/development.md). diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 0000000..9add397 --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,90 @@ +# Contributing to SpecMod + +The working guide is [`docs/development.md`](docs/development.md) — setup, the +daily loop, the tooling, the CI jobs, and how development relates to releases. +If you are an AI coding agent, read [`AGENTS.md`](AGENTS.md) first. + +This file covers the one thing that is a promise rather than a practice. + +## The stability promise, and its exact scope + +SpecMod is `0.x` and says so loudly: names and signatures move between minor +releases without a deprecation cycle, because the API is still being worked +out and shipping shims for names that are about to move again costs more than +it protects. + +**`specmod.api` is the exception.** It exists so that downstream packages have +something that does not move while the internals do, and it is the only part of +SpecMod that carries a compatibility guarantee. + +| | Everything else | `specmod.api` | +|---|---|---| +| Rename or remove | any minor release, no notice | one minor release of `DeprecationWarning` first | +| Change a signature | any minor release, no notice | one minor release of `DeprecationWarning` first | +| Add | freely | freely, but see below | +| Behaviour change that moves a number | called out in the changelog | called out in the changelog **and** in the warning | + +"One minor release" means: if a removal is decided during `0.4.x`, the warning +ships in `0.5.0` and the removal is `0.6.0` at the earliest. A deprecation that +has not been through a release has not been announced. + +### Deprecating something in `specmod.api` + +```python +warnings.warn( + "specmod.api.old_name is deprecated and will be removed in 0.6.0; " + "use specmod.api.new_name, which takes the same arguments.", + DeprecationWarning, + stacklevel=2, +) +``` + +Three things that make the difference between a warning people act on and one +they filter out: + +- **Name the replacement**, or say plainly that there is none. +- **Name the release it goes in**, not "a future version". +- **`stacklevel=2`**, so the warning points at the caller's line rather than at + SpecMod's. + +Keep the old name working for the whole cycle. A `DeprecationWarning` on +something that already raises is not a deprecation, it is a breakage with a +note attached. + +### Adding to `specmod.api` + +Every export is a compatibility obligation, so the surface is deliberately +small and does not grow opportunistically. To add one: + +1. Have a caller that needs it. "Studio might want this" is not one. +2. Make it satisfy the five properties in the module docstring — path-free, + deterministic, non-mutating, quiet, typed errors. If the underlying + function does not, the wrapper is where that gets fixed, not the caller. +3. Add it to `EXPECTED_EXPORTS` in `tests/test_api_surface.py`. The list is + duplicated there on purpose, so an addition shows up as a diff in review + rather than as a passing test. +4. Give it a docstring with `Parameters`, `Returns` and `Raises`, and a type + annotation on everything. Both are tested. + +**Do not reach around the surface.** If a downstream package needs something +`specmod.api` does not export, extend `specmod.api` in its own pull request +with the reason stated. An import of `specmod.core` or `specmod.fitting` from +a downstream package is a bug in the boundary, not a shortcut. + +### What is *not* promised + +- **Internals.** `specmod.core`, `specmod.fitting`, `specmod.transforms`, + `specmod.picks`, `specmod.config` and everything else may change in any + release. They are documented because SpecMod's own users read them; that is + not a stability claim. +- **Numerical output.** The promise is about names and signatures. A bug fix + that moves a number is still a bug fix, and it will move it — that is what + the golden references and the changelog are for. If you depend on exact + values, pin an exact version and keep the config hash. +- **The objects the surface returns**, beyond the attributes its docstrings + name. `SpectrumPair` gaining a field is not a breaking change. + +## Everything else + +Conventional Commits, `pre-commit install` before your first commit, and the +rest of the mechanics are in [`docs/development.md`](docs/development.md). diff --git a/README.md b/README.md index 3f6b343..8f4af57 100644 --- a/README.md +++ b/README.md @@ -7,12 +7,17 @@ SpecMod estimates source parameters — long-period spectral level Ω, corner frequency `f_c`, and the attenuation operator `t*` — by fitting a Brune-type source model to direct-phase spectra. -> **Status: under active reconstruction.** -> The package is mid-refactor. The modern layers (`specmod.config`, -> `specmod.core`, `specmod.transforms`) are built and tested; the older -> pipeline modules still carry pre-refactor behaviour and are being replaced -> stage by stage. Expect breaking changes at every `0.x` release until the API -> settles at 1.0. See [`docs/REFACTOR_PLAN.md`](docs/REFACTOR_PLAN.md). +> **Status: alpha, and under active reconstruction.** +> The package is pre-1.0 and mid-refactor. The modern layers +> (`specmod.config`, `specmod.core`, `specmod.transforms`, `specmod.picks`, +> `specmod.fitting`) are built and tested; the older pipeline modules still +> carry pre-refactor behaviour and are being replaced stage by stage. Expect +> breaking changes at every `0.x` release until the API settles at 1.0 — they +> land in minor bumps by design, with no deprecation cycle. Pin an exact +> version for anything you intend to publish. +> [`docs/roadmap.md`](docs/roadmap.md) says which stages are done and what 1.0 +> will mean; [`docs/REFACTOR_PLAN.md`](docs/REFACTOR_PLAN.md) is the working +> document behind it. ## Installation @@ -141,10 +146,25 @@ Every output records the configuration that produced it, a hash of it, and the SpecMod version, so a locally-overridden run is still reproducible from its outputs. +## Documentation + +The full documentation — the pipeline with its equations, the estimator +comparison, pick formats, and an API reference — builds with Sphinx: + +```bash +uv pip install -e '.[docs]' +sphinx-build -b html docs docs/_build/html +``` + +`docs/REFACTOR_PLAN.md` is excluded from the built site on purpose: it is a +working document that records decisions and the measurements behind them, not +documentation for using the package. + ## Development ```sh uv venv && uv pip install -e ".[dev]" +pre-commit install # both hook types; not optional pytest # test suite pytest --without-optional-extras # as a default install and CI see it ruff check src/ tests/ tools/ # lint @@ -152,6 +172,11 @@ ruff format src/ tests/ tools/ mypy # strict on the rewritten modules ``` +[`docs/development.md`](docs/development.md) is the full guide — the repository +mapped, every tool and CI check, the branch and commit conventions, and where +development stops and releasing begins. [`AGENTS.md`](AGENTS.md) is the short +version that binds AI coding sessions. + Run `--without-optional-extras` before pushing. A development environment with `specmod[multitaper]` installed will pass tests that a default install fails, and CI installs only `[dev]`. @@ -178,6 +203,15 @@ published number fails the suite rather than quietly leaving the prose wrong. Measurements that read `tutorial/data/events/` are slower and opt-in via `--field`; refresh those by hand after changing an estimator. +### Releasing + +Commit messages follow [Conventional Commits](https://www.conventionalcommits.org), +which is what makes the changelog and the version automatic: `release-please` +opens a standing release pull request, and merging it creates the tag, the +GitHub Release, the PyPI upload and the Zenodo DOI. Nothing is released until +that merge. See [`docs/releasing.md`](docs/releasing.md), which also lists the +repository settings that have to be turned on once. + ## References Edwards, B., Allmann, B., Fäh, D., Clinton, J. (2010). Automatic computation of diff --git a/ci/README.md b/ci/README.md index 85c639a..66479c4 100644 --- a/ci/README.md +++ b/ci/README.md @@ -24,6 +24,31 @@ Copy the whole file over its counterpart. No merging, no partial application: | staged copy | destination | |---|---| | `ci/workflows/test.yml` | `.github/workflows/test.yml` | +| `ci/workflows/docs.yml` | `.github/workflows/docs.yml` | +| `ci/workflows/release.yml` | `.github/workflows/release.yml` | + +One of them needs a repository setting turning on as well, and it cannot be +done from a commit: + +- `release.yml` — **Settings → Actions → General → Allow GitHub Actions to + create and approve pull requests**. Without it release-please fails with + `GitHub Actions is not permitted to create or approve pull requests`. + +`docs.yml` needs nothing: it builds the site as a check and publishes nothing. +Read the Docs publishes, and its setup lives in +[`docs/documentation.md`](../docs/documentation.md). + +`release.yml` needs four more one-time steps before it can publish anything — +a `pypi` environment, a trusted publisher registered on PyPI, the Zenodo +webhook, and branch protection. They are listed in +[`docs/releasing.md`](../docs/releasing.md). Until they are done the workflow +opens a release pull request and stops there, which is inert rather than +wrong. + +**The PyPI trusted publisher names this file.** It is registered against the +workflow filename `release.yml`, matched from the OIDC token, so renaming the +workflow breaks authentication at upload time. Rename both in the same +sitting or not at all. A file here is the **intended** state, which is not necessarily the current one. `tools/check_ci_mirror.py` reports which of the two each pair is in, and diff --git a/ci/workflows/docs.yml b/ci/workflows/docs.yml new file mode 100644 index 0000000..65eac67 --- /dev/null +++ b/ci/workflows/docs.yml @@ -0,0 +1,47 @@ +name: docs + +on: + push: + branches: [main] + pull_request: + workflow_dispatch: + +# Read the Docs publishes the site; this job only checks that it builds. It is +# kept alongside Read the Docs' own pull request build because it is the fast, +# in-repository check that does not depend on a third-party service being up — +# and because it is an import check, which is worth having in CI regardless. +# See docs/documentation.md. + +concurrency: + group: docs-${{ github.ref }} + cancel-in-progress: true + +jobs: + build: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + with: + fetch-depth: 0 # hatch-vcs derives the version from tags + - uses: astral-sh/setup-uv@v5 + - run: uv python install 3.11 + - run: uv venv # create a venv to install into + + # The package itself, not just the docs extra: autodoc imports every + # module it documents, so a docs build is also an import check. + - run: uv pip install -e ".[docs,io]" + + - name: Build + # No -W. Intersphinx resolves seven inventories over the network and + # warns when one is briefly unreachable, which would turn a third + # party's downtime into a red build. A genuinely broken build exits + # non-zero on its own. `.readthedocs.yaml` sets fail_on_warning: false + # for the same reason — keep the two in step. + run: uv run sphinx-build -b html docs docs/_build/html + + # Not a Pages artefact: this is a plain zip anyone can download from the + # run and open locally, and it stays useful if Read the Docs is down. + - uses: actions/upload-artifact@v4 + with: + name: docs-html + path: docs/_build/html diff --git a/ci/workflows/release.yml b/ci/workflows/release.yml new file mode 100644 index 0000000..362dccb --- /dev/null +++ b/ci/workflows/release.yml @@ -0,0 +1,76 @@ +name: release + +on: + push: + branches: [main] + workflow_dispatch: + +# release-please pushes the release branch, opens the PR, and on merge creates +# the tag and the GitHub Release. The publish job below needs neither: PyPI +# Trusted Publishing authenticates with an OIDC token, so nothing long-lived +# lives in secrets. +permissions: + contents: write + pull-requests: write + +# Two pushes to main should not race to update the same release PR. +concurrency: + group: release + cancel-in-progress: false + +jobs: + release-please: + runs-on: ubuntu-latest + outputs: + release_created: ${{ steps.release.outputs.release_created }} + tag_name: ${{ steps.release.outputs.tag_name }} + steps: + - uses: googleapis/release-please-action@v4 + id: release + with: + config-file: release-please-config.json + manifest-file: .release-please-manifest.json + + # Deliberately in this workflow rather than a separate publish.yml keyed on + # `release: published`, which is what §6.5 of the plan used to describe. + # release-please creates the release with the default GITHUB_TOKEN, and + # per GitHub's docs "events triggered by the GITHUB_TOKEN will not create a + # new workflow run" — so that publish workflow would never fire. The + # alternative is a personal access token in secrets, which is the one thing + # Trusted Publishing exists to avoid. Gating on the action's own output + # keeps the token count at zero. + # + # Zenodo is unaffected: it listens to the release *webhook*, which is not + # subject to that restriction, so the DOI is still minted from the release. + publish: + needs: release-please + if: needs.release-please.outputs.release_created == 'true' + runs-on: ubuntu-latest + environment: + name: pypi + url: https://pypi.org/p/specmod + permissions: + id-token: write # mandatory for Trusted Publishing + contents: read + steps: + - uses: actions/checkout@v4 + with: + # The tag, not the branch. main has already moved on by one commit + # (the release PR merge), and hatch-vcs would derive a .postN.devN + # version from it. + ref: ${{ needs.release-please.outputs.tag_name }} + fetch-depth: 0 # hatch-vcs derives the version by describing tags + + - uses: astral-sh/setup-uv@v5 + - run: uv python install 3.11 + - run: uv build + + # The tag is the version, and pyproject's tag_regex decides whether + # hatch-vcs can read it. If those two ever disagree the wheel is built + # as 0.1.1.postN.devN and would be uploaded under that name — PyPI does + # not let it be taken back. tests/test_release_config.py checks the + # formats agree; this checks the artefact that is about to be published. + - name: The built version is the tag + run: python tools/check_built_version.py "${{ needs.release-please.outputs.tag_name }}" dist + + - uses: pypa/gh-action-pypi-publish@release/v1 diff --git a/docs/REFACTOR_PLAN.md b/docs/REFACTOR_PLAN.md index 74b3a9f..4c00618 100644 --- a/docs/REFACTOR_PLAN.md +++ b/docs/REFACTOR_PLAN.md @@ -628,6 +628,58 @@ failing checks to one, and the survivor was `cwt` alone, disagreeing by 1-2% on four stations. Small enough to look like tolerance, but it was another branch — the "already touching" guard inside the lift. +**The `cwt` residual, measured.** After those four, `cwt` was still held at +`rtol = 5e-2` on CI for a disagreement nobody could explain. Measured on a box +matching the reference environment exactly, every proposed mechanism fails: + +| Hypothesis | Measurement | Verdict | +|---|---|---| +| Linux disagrees with a Linux-generated reference | `cwt` reproduces it to **3.8e-16** here | it is that *runner*, not Linux | +| A remaining discontinuity in the `cwt` path | 1e-13 in moves the noise 8.7e-14 out, **no step, all 28 windows** | linear; so 1-2% out needs 1-2% *in*, which floating point cannot supply | +| A sample sitting on a bin edge | closest `cwt` sample is **5.7e-4 of a bin** from an interior edge; `fft` 1.4e-6 | not fragile, and `fft` is the fragile one | +| The window differs on CI | one sample fewer moves `fft` noise **8.5%** and `cwt` **3.6%** | ruled out — `fft` would disagree first, and it agrees exactly | +| Quantile fragility from short arrays (51 vs 109) | 1e-15 moves `cwt`'s worst quantile **2.4e-14**, `fft`'s 5.1e-13 | `cwt` is the *more* stable of the two | +| PyWavelets, or threading | the estimator is hand-written; the transform is a batched `numpy.fft.ifft` | neither applies | + +So the CWT is as stable as everything else, and `5e-2` was twelve orders of +magnitude looser than anything reproducible — wide enough to hide the real +difference rather than describe it. It was set to **`1e-3`** as an experiment: +if the residual was still there CI would report it with per-station detail, +and if it was not the entry could go. Either outcome is worth more than the +number that was there. + +**The experiment's answer: the residual is real, and it is the machine.** One +CI run, six test jobs, one commit: + +| Job | Result | +|---|---| +| ubuntu 3.11 | fails — 8 differences, 3 windows, worst 1.44e-2 | +| ubuntu 3.13 | fails — the *same* 8, the same windows, the same magnitudes to 3 s.f. | +| ubuntu 3.12 | passes | +| macOS 3.11, 3.12, 3.13 | pass | + +Two things follow that the earlier round could not establish. It is +deterministic rather than flaky, since two jobs on different Python versions +produce identical numbers. And it is not the OS, not the Python version, and +not a package version that tracks the Python version — a third ubuntu job in +the same run agrees with the reference exactly. What is left is which machine +the job landed on, which is what "it is that *runner*, not Linux" suspected +and this is the same-run control for. + +The sharpest datum is the one that localises it: on the machine that +disagrees, `fft`, `welch`, `multitaper` and `quadratic` all still reproduce +the reference byte for byte. It is in the CWT path, not in that machine's +arithmetic in general. + +The tolerance is now **`2e-2`** — the worst observed difference plus about 40% +of headroom, and 2.5 times tighter than the `5e-2` it replaces. It is a bound +on the residual, not an explanation of it, and it is labelled that way in the +test: while the mechanism is unidentified and reproduces on no box available +to work on, a bound with numbers behind it is the most that can honestly be +claimed. Taking it further needs the failing machine — the cwt path run on it +and on a passing one, diffing the three named windows' noise arrays *before* +binning rather than after. + **Measured after the fix**, end to end over all 28 windows: perturbing the input by 1e-15 moves the noise by 1.8e-11 and **no band edge at all**. The response is linear. The golden reference's exact noise comparison, which had @@ -2237,14 +2289,42 @@ numbers: record; the rest is the taper. That is **0.03 in Mw** — negligible, but one-directional, so it is asserted as a signed bias rather than absorbed into a symmetric tolerance. -- **The fit belongs on velocity, not displacement.** The model carries a motion - factor, so `llpsp` is the displacement plateau either way — but - `initial_guess` takes the spectral *peak* as the `f_c` guess, and a - displacement spectrum falls monotonically. Fitting the displacement set puts - the guess at the low band edge and the fit settles near it: `f_c` came back - **1.6 instead of 8.0**, and `Ω₀` 0.6 magnitude units low. Nothing warns. That - is a sharp edge on a public API, and a candidate for `initial_guess` to - either detect the motion or refuse it. +- **The fit is done in the sensor's natural units** — velocity for these + instruments — and that is a general rule rather than a quirk of this + pipeline. The model carries a motion factor, so `llpsp` is the displacement + plateau whichever domain is fitted; what changes is the record being fitted + and the noise on it. + + Converting first is not a neutral change of view. Integrating to + displacement divides by `2πf`, which implicitly low-passes: it damps + high-frequency noise, but it also destroys the feature the guess is read + from. Differentiating to acceleration multiplies by `2πf` and blows the + high-frequency noise up instead. + + Velocity is also the convenient domain, because **it peaks at `f_c`** — so + `initial_guess` taking the spectral peak is exact rather than approximate. + Maximising `f·[1 + (f/f_c)^{γn}]^{-1/γ}` puts the stationary point at + `(n − 1)·(f/f_c)^{γn} = 1`, which is `f = f_c` **whenever `n = 2`, for any + `γ`** — the corner's sharpness does not enter. So it holds for the whole + omega-squared family, not just Brune: measured on a 400001-point grid, both + registered sources peak at 8.0000 Hz for a true 8.0, `brune` (`γ=1`) and + `boatwright` (`γ=2`) alike. + + That is worth knowing before a source model with `n ≠ 2` is registered, + because the guess silently stops being exact at that point — the peak moves + to `f_c·(n−1)^{−1/(γn)}`. + + In displacement and acceleration the spectrum is monotonic over the band, so + the "peak" is just whichever edge it was handed — 0.1 Hz and 100 Hz on the + same grid. + + So `FitSpectra(spectra.to_motion("displacement"))` got a guess at the bottom + of the band and settled near it: `f_c` came back **1.6 instead of 8.0**, `Ω₀` + 0.6 magnitude units low, with nothing warning. The design was right and the + silence was not, so `initial_guess` now warns when the spectrum's motion is + not one the peak means anything in. `FittableView` grew a `motion` property + to make that possible — it reads through to the pair like the rest, and a + fitter has to know the domain for exactly this reason. **Tier 3 — golden/regression.** Run the *current* code on the tutorial event and on Magna (§5.2.4), and snapshot `freq`, `amp`, `bsnr`, `ubfreqs` and the fit @@ -3085,17 +3165,49 @@ existing only in the config file. ### 6.3 Documentation (Sphinx) +**Built**, in `docs/conf.py`, `.readthedocs.yaml` and `ci/workflows/docs.yml`. +What follows is the plan as written, annotated with what building it changed. + - **Sphinx** with `pydata-sphinx-theme` (the NumPy/SciPy/ObsPy house style — familiar to this audience, good API-reference layout). - `myst-parser` so prose pages can stay Markdown; `myst-nb` to execute and render the tutorial notebook as a docs page, which makes the tutorial a **tested** - artefact rather than a snapshot that silently rots. -- `autodoc` + `napoleon` (numpydoc style) + `sphinx-autodoc-typehints`, so - signatures come from the annotations rather than being hand-maintained. -- `intersphinx` to numpy, scipy, obspy, lmfit, matplotlib. + artefact rather than a snapshot that silently rots. `myst-nb` waits for + Phase 6; until then `docs/notebooks/` is excluded, because copying an + unexecuted notebook into the site is worse than leaving it out. +- `autodoc` + `napoleon` (numpydoc style), signatures coming from the + annotations rather than being hand-maintained. **Not** + `sphinx-autodoc-typehints`, which this line named until it was measured: + with `autodoc_typehints = "description"` the built-in extension produced the + same 367 documented objects, while the third-party one calls an API Sphinx + 10 removes and emits a deprecation warning per module. Dropped, after + checking rather than assuming — removing a dependency that was quietly doing + work would have been a bad trade. +- `intersphinx` to python, numpy, scipy, pandas, matplotlib, obspy and lmfit. - `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.** + exactly the kind of thing that should be executable in the docs. Not turned + on yet; the examples have to be written as doctests first. + +Four things the first build found, each fixed rather than tolerated: + +- **Ambiguous cross-references.** Documenting a package *and* its submodules + gave every re-exported name two targets, so `PickSet`, `SensorID`, + `Resolution` and `FitSpectra` were all ambiguous. Packages are now documented + at the path you import from. +- **A link into `REFACTOR_PLAN.md`**, which is excluded on purpose — it is a + working document, not documentation. `processing.md` now links to it on + GitHub instead. +- **`notes/` was excluded on that same reasoning, and that was wrong.** + `choosing_a_transform.md` links to `notes/window_position.md` for a per-trace + table, which makes it documentation. Now built, and reachable through a + hidden toctree on the page that cites it. +- **`HOLT_2019_UTAH` broke autodoc.** It is a callable dataclass instance + documented as module data, and autodoc's signature formatter raises on it + where `inspect.signature` handles it fine. Excluded, with the value written + out in prose. + +**The equations in `docs/` did not render before this, and building it was 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 @@ -3103,35 +3215,73 @@ 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 +Two things were listed here to do at build time, since neither was verifiable +without a renderer to check against. Measured against the built site: + +- **`dollarmath` was necessary, `amsmath` was not.** `myst_enable_extensions` + turns on `dollarmath` (plus `colon_fence`, `deflist` and `substitution`). + `amsmath` covers bare `\begin{align}` outside `$` delimiters, and there are + none — the one `\begin{cases}` in `processing.md` sits inside `$$` and + renders without it. +- **The syntax predicted to break does not.** The prediction was one display + block without a preceding blank line, one spanning two lines, and 26 inline + expressions containing underscores. Built: `processing.html` and + `choosing_a_transform.html` contain **no** literal `$` at all and 120 math + nodes between them, the two-line block and the `cases` block included. The + risk was real on GitHub's renderer; MyST's parses all of it. + +The general point stands anyway, 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 - autodoc emits none for a symbol it was never asked to document. Catching - that needs `sphinx.ext.coverage` with `coverage_show_missing_items`, or - `nitpicky` for unresolved references. Worth having; it is a different - setting, and writing it as a property of `-W` would have meant discovering - the gap only after trusting it. +were edited a dozen times in this refactor against a renderer nobody had run, +and running it found four separate faults (above) even though the equations +turned out fine. + +- **No `sphinx-build -W`**, though this line asked for it twice. `-W` promotes + warnings Sphinx emits — and intersphinx emits one whenever any of the seven + inventories is briefly unreachable, which makes a third party's downtime a + red build. That is flakiness, not a check; a genuinely broken build exits + non-zero on its own. (The line before that claimed `-W` fails on an + undocumented public symbol, which it does not: autodoc emits no warning for + a symbol it was never asked to document. That needs `sphinx.ext.coverage` + with `coverage_show_missing_items`, or `nitpicky` for unresolved references. + Still worth having, and still a different setting.) - Structure: Getting started → User guide (preprocessing, transforms, SNR, fitting) → **Theory** (the normalisation conventions, one page, with the Parseval contract stated explicitly) → Tutorial → API reference → Migration - guide from 0.x → Changelog. + guide from 0.x → Changelog. Built so far: an index, the four existing prose + pages plus `notes/`, `releasing.md`, `roadmap.md`, `development.md`, + `documentation.md`, and the API reference. Theory, tutorial and migration + pages are Phase 6. + +**Published by Read the Docs, not GitHub Pages.** The plan said Pages, and +Pages was built first; it was replaced before anything was deployed, which is +the cheapest moment to change a decision like this. The reason is versions. +`actions/deploy-pages` publishes an artefact that *becomes* the whole site, so +every deploy replaces everything — one site, always showing `main`. That is +incompatible with telling alpha users to pin an exact version, because the +documentation they need is then never the documentation they get. Keeping +history on Pages means either accumulating directories on a `gh-pages` branch +or rebuilding every tag on every deploy; Read the Docs does versions per tag, a +`stable` alias, a switcher, cross-version search and per-pull-request previews +without any of that machinery. + +What stays in CI is the `docs` job, reduced to a build check. It is worth +keeping for two reasons that survive the move: it does not depend on a third +party being up, and autodoc imports every module it documents, so it is an +import check as much as a docs check. + +One trap, recorded because it is silent rather than loud: Read the Docs clones +shallow and without tags, and `hatch-vcs` derives the version from +`git describe`. Without the `post_checkout` unshallow in `.readthedocs.yaml` +every build reports the `0.0.0` fallback in the sidebar — including tagged +ones, which is exactly where it would be believed. + +`roadmap.md` is this section's §7 restated for a reader rather than for +whoever is doing the work: stages, no durations, and no phase numbers. The +durations here are estimates that have already been wrong; publishing them on +the site would turn an estimate into a promise. It says explicitly that the +stages become milestones against released versions once there are releases to +anchor them to. The theory page matters more than usual here. The units question that prompted this refactor is not obvious from the code, and if the conventions are only @@ -3149,8 +3299,8 @@ version string is ever committed, so there is nothing to forget to bump and no **Version *decision* — `release-please` (GitHub Action).** It parses [Conventional Commits](https://www.conventionalcommits.org) since the last release, works out the SemVer bump, and opens a standing "release PR" carrying -the generated `CHANGELOG.md`. Merging that PR creates the tag; the tag triggers -publication. Nothing is released until a human merges. +the generated `CHANGELOG.md`. Merging that PR creates the tag and the GitHub +Release. Nothing is released until a human merges. That human gate is the reason to prefer `release-please` over `python-semantic-release` (which tags on every qualifying push to `main`) @@ -3160,23 +3310,71 @@ mint a citable version of the software. Use `python-semantic-release` only if yo would rather have zero-touch releases and accept that. This does impose Conventional Commits (`feat:`, `fix:`, `refactor:`, `docs:`, -`feat!:` for breaking) on commit messages, enforced by a `commitlint` pre-commit -hook. It is a small discipline and it is what makes the changelog automatic. +`feat!:` for breaking) on commit messages. The convention is followed by hand: +there is no `commitlint` hook, and this paragraph claimed one until §6.6 went +looking. It is a small discipline and it is what makes the changelog automatic. + +**Built, in `release-please-config.json` and `ci/workflows/release.yml`.** +Three settings there are load-bearing, and each was chosen against a measured +consequence rather than a default: + +- **`include-component-in-tag: false`.** Left on, release-please tags + `specmod-v0.2.0`, which `pyproject.toml`'s `--match v[0-9]*` does not + describe and its `tag_regex` does not parse. The wheel would then be built + as `0.1.1.postN.devN` and uploaded under that name, which PyPI does not let + you take back. `tests/test_release_config.py` asserts the two formats agree; + `tools/check_built_version.py` re-checks the artefact between the build and + the upload, because a configuration test cannot see a build. +- **`bump-minor-pre-major: true`.** The history holds two breaking commits + (`30b4e89`, `628d36d`). Without this, either one proposes `1.0.0` — a + version that says the API has stopped moving, with a DOI attached. +- **Explicit `changelog-sections`.** The default preset hides `refactor`, + `docs`, `build`, `test`, `ci`, `style` and `chore`. Over this repository's + 146 conventional commits that prints 73 and drops 73, so a release that is + mostly refactoring would ship an almost empty changelog. `refactor`, `docs` + and `build` are shown. + +The `simple` release type is what suits a project with no version string to +update: its only other updater targets `version.txt` with `createIfMissing: +false`, so with no such file it writes the changelog and nothing else. Checked +in release-please's source rather than assumed, since a strategy that created +a second source of truth for the version would defeat `hatch-vcs`. + +The manifest starts at `0.1.1` — the version the Magna paper cites, and the +last one this code had. There are no tags in the repository at all, so that is +a statement about history rather than something derived; the first release PR +therefore proposes `0.2.0` and carries the whole refactor as its changelog, +which is what the delta from 0.1.1 actually is. ### 6.5 CI (GitHub Actions) | Workflow | Trigger | Does | |---|---|---| | `test.yml` | PR, push | `lint` (ruff check + format), `typecheck` (mypy), `test` matrix 3.11/3.12/3.13 × ubuntu/macos (pytest + coverage → Codecov), `floors` (`--resolution lowest-direct`, see §6.6). Five job names, not one matrix — the row used to describe them as a single matrix step | -| `docs.yml` | PR, push | `sphinx-build -W`; on `main`, deploy to GitHub Pages. Builds on PRs too, so doc breakage is caught before merge | +| `docs.yml` | PR, push | `sphinx-build` (no `-W`, see §6.3), as a check only — Read the Docs publishes. Also an import check, since autodoc imports every module it documents | | `build.yml` | PR, push | sdist + wheel, `twine check`, install-from-wheel smoke test in a clean env (catches missing package data) | -| `release-please.yml` | push to `main` | maintains the release PR; creates tag + GitHub Release on merge | -| `publish.yml` | GitHub Release published | PyPI via Trusted Publishing (OIDC — no long-lived token in secrets) | - -Zenodo is wired to the GitHub Release webhook, so the DOI is minted from the same +| `release.yml` | push to `main` | `release-please` maintains the release PR and creates tag + GitHub Release on merge; a gated `publish` job then builds from the tag and uploads to PyPI via Trusted Publishing (OIDC — no long-lived token in secrets) | + +**That last row was two workflows until it was built.** The plan had +`publish.yml` triggered by `release: published`, which never fires: +release-please creates the release with the default `GITHUB_TOKEN`, and per +GitHub's documentation "events triggered by the `GITHUB_TOKEN` will not create +a new workflow run". The fix is either a personal access token in secrets — +the one thing Trusted Publishing exists to avoid — or putting the publish job +in the same workflow, gated on release-please's `release_created` output. +The second, so the token count stays at zero. + +Zenodo is unaffected by that restriction, because it listens to the release +*webhook* rather than running a workflow: the DOI is still minted from the same event as the PyPI upload. Branch protection on `main`: require `test`, `docs` and `build` green. +Six repository settings have to be turned on by hand before a release can +happen — Actions-may-open-PRs, the `pypi` environment, the PyPI trusted +publisher, the Zenodo webhook, and branch protection. None of them is +expressible in a commit, so they are written down in `docs/releasing.md` +instead of being remembered. + **Versioning policy.** SemVer. Stay on `0.x` for the whole refactor — breaking changes are expected and permitted in minor bumps there, so no deprecation cycle is needed (§3.1). Tag `v0.2.0` at the end of Phase 2 to prove the pipeline; @@ -3213,6 +3411,8 @@ been audited once is worth more than one that reads confidently throughout. | "Conventional Commits, **`commitlint`-enforced**" (§7) | No `commitlint` hook. The convention is followed by hand. §7 corrected. | | "`sphinx-build -W` … an **undocumented public symbol** fails the build" | `-W` promotes warnings Sphinx emits; autodoc emits none for a symbol it was never asked to document. Needs `sphinx.ext.coverage` or `nitpicky`. §6.3 corrected. | | `test.yml` "matrix … → ruff, mypy, pytest" | Five jobs, not one matrix: `lint`, `typecheck`, `test`, `floors`, plus `build.yml`. §6.5 corrected. | +| "`commitlint`-enforced" again, in §6.4 | The §7 instance was corrected and this one was left standing, in the same document. §6.4 corrected. | +| `publish.yml` "GitHub Release published → PyPI" | It would never have run. release-please creates the release with the default `GITHUB_TOKEN`, and events triggered by that token do not start a workflow. The publish job moved into `release.yml`, gated on `release_created`. §6.5 corrected. | **A claim whose mechanism is absent but whose property holds — for a different reason, which matters:** @@ -3368,10 +3568,10 @@ 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.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 | +| **2b. Release plumbing** ✅ | ~~Sphinx skeleton + `pydata-sphinx-theme` + autodoc/napoleon/intersphinx~~ ✅; ~~`docs.yml` → a published site~~ ✅ — Read the Docs rather than GH Pages, for versions (§6.3); ~~release-please + PyPI Trusted Publishing~~ ✅ — one `release.yml`, not two workflows (§6.5); ~~Zenodo webhook~~ ✅ documented. All three workflows are staged in `ci/` and need copying across, and six repository settings have to be turned on by hand: `docs/releasing.md` lists them. 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 | -| **5. Decompose** | Split `Spectral.py` (655 lines) into `core/` + `snr/`; `Fitting.py` → `fitting/`; models as objects; `io/`; `viz/`; non-mutating operations; mypy override list → empty; `picks/` per §4.9 — resolution first, then the registry, the ObsPy delegate and the plugin entry point | 3 | 4–6 days | +| **5. Decompose** | Split `Spectral.py` (655 lines) into `core/` + `snr/` ✅; ~~`Fitting.py` → `fitting/`~~ ✅ — `base`/`guess`/`spectrum`/`event`, 745 lines to four modules of 30–300, public names unchanged; models as objects ✅; `io/`; `viz/`; non-mutating operations; mypy override list → empty ✅; `picks/` per §4.9 ✅ | 3 | 4–6 days | | **6. Ship** | Full docs content, tutorial rewritten as an executed `myst-nb` page with no `os.chdir`, 0.1→1.0 "what changed" page, **1.0 release** | 4, 5 | 2–3 days | **Executing the tutorial was pulled forward out of Phase 6.** The `myst-nb` diff --git a/docs/api.md b/docs/api.md new file mode 100644 index 0000000..bf05f16 --- /dev/null +++ b/docs/api.md @@ -0,0 +1,121 @@ +# API reference + +Grouped by what a run does in order: get the data, cut it, transform it, fit +it, write it out. + +**If you are writing a package on top of SpecMod, start with +[`specmod.api`](#the-stable-surface) instead.** It is a small, frozen subset of +what follows, and the only part that carries a compatibility promise — the rest +of this page documents internals that move between `0.x` releases. See +[`CONTRIBUTING.md`](https://github.com/sgjholt/SpecMod/blob/main/CONTRIBUTING.md) +for the exact scope of that promise. + +## The stable surface + +```{eval-rst} +.. automodule:: specmod.api + :exclude-members: AmplitudeKind, Config, InternalError, InvalidInputError, + MissingBackendError, Motion, ResolvedConfig, SpecModError, + Spectrum, SpectrumPair, config_hash, load_config, + make_window, window_correction +.. automodule:: specmod.exceptions +``` + +The names excluded above are re-exports, documented at the path they are +defined — `Spectrum` and `SpectrumPair` under [Spectra](#spectra), `Config` +and `load_config` under [Configuration](#configuration), `make_window` and +`window_correction` under [Transforms](#transforms). Documenting them twice +gives every cross-reference to them two targets and makes all of them +ambiguous, which is the same trap package-level `automodule` set earlier on +this page. `specmod.api.__all__` is the authoritative list, and +`tests/test_api_surface.py` asserts it. + +Packages are documented at the path you import from — `specmod.picks.PickSet`, +not `specmod.picks.base.PickSet`. Documenting both the package and its +submodules gave every re-exported name two targets and made every +cross-reference to it ambiguous. + +## Getting data + +```{eval-rst} +.. automodule:: specmod.datasets +.. automodule:: specmod.acquire +``` + +## Picks + +```{eval-rst} +.. automodule:: specmod.picks + :imported-members: +``` + +## Preparing waveforms + +```{eval-rst} +.. automodule:: specmod.preprocess +``` + +## Spectra + +```{eval-rst} +.. automodule:: specmod.pipeline +.. automodule:: specmod.core.spectrum +.. automodule:: specmod.core.collection +.. automodule:: specmod.core.units +.. automodule:: specmod.core.scalogram +``` + +## Transforms and smoothing + +```{eval-rst} +.. automodule:: specmod.transforms + :imported-members: +.. automodule:: specmod.transforms.base +.. automodule:: specmod.smoothing + :imported-members: +``` + +## Source models + +```{eval-rst} +.. automodule:: specmod.sources + :imported-members: +``` + +## Fitting + +```{eval-rst} +.. automodule:: specmod.fitting + :imported-members: +.. automodule:: specmod.staged +``` + +## Magnitude + +```{eval-rst} +.. automodule:: specmod.magnitude +.. automodule:: specmod.spreading + :exclude-members: HOLT_2019_UTAH +``` + +`HOLT_2019_UTAH` is the piecewise model fitted in Holt (2019), pre-built: +`Piecewise(segments=((0.90, 43.0), (2.57, 76.0), (0.44, 136.0), (1.54, 400.0)))`. +It is excluded above because autodoc cannot format a signature for a callable +dataclass *instance* documented as module data — `inspect.signature` handles it +fine, autodoc's own formatter raises. Excluding one constant is cheaper than +working around that. + +## Output + +```{eval-rst} +.. automodule:: specmod.io +.. automodule:: specmod.tables +.. automodule:: specmod.plotting +``` + +## Configuration + +```{eval-rst} +.. automodule:: specmod.config +.. automodule:: specmod.utils +``` diff --git a/docs/choosing_a_transform.md b/docs/choosing_a_transform.md index d9f0684..62f965b 100644 --- a/docs/choosing_a_transform.md +++ b/docs/choosing_a_transform.md @@ -600,3 +600,9 @@ normalize_to_variance = true # mtspec's convention The first of those is unverified — see `docs/REFACTOR_PLAN.md` §5.2.5. The 0.1.1 re-run tests it directly. + +```{toctree} +:hidden: + +notes/window_position +``` diff --git a/docs/conf.py b/docs/conf.py new file mode 100644 index 0000000..519888a --- /dev/null +++ b/docs/conf.py @@ -0,0 +1,110 @@ +"""Sphinx configuration for the SpecMod documentation. + +Build with:: + + uv pip install -e '.[docs]' + sphinx-build -b html docs docs/_build/html + +``-W`` is deliberately **not** used. Intersphinx resolves seven inventories +over the network, and a warning is emitted whenever one of them is briefly +unreachable — turning a third party's downtime into a red build. The docs job +fails on a non-zero exit instead, which is what a genuinely broken build +gives. +""" + +from __future__ import annotations + +from importlib.metadata import PackageNotFoundError, version + +project = "SpecMod" +author = "James Holt" +#: Sphinx substitutes `%Y` with the build year, so this does not need editing +#: — and it honours `SOURCE_DATE_EPOCH`, so a reproducible build stamps the +#: source date rather than the day it happened to run. 2020 is the year in +#: `LICENSE` and the year the history starts. +copyright = "2020-%Y, James Holt" + +try: + release = version("specmod") +except PackageNotFoundError: # pragma: no cover - docs built without an install + release = "0.0.0" +#: The short X.Y form, which is what the sidebar shows. +version = ".".join(release.split(".")[:2]) + +extensions = [ + "myst_parser", + "sphinx.ext.autodoc", + "sphinx.ext.napoleon", + "sphinx.ext.intersphinx", + "sphinx.ext.viewcode", +] +#: `sphinx_autodoc_typehints` is deliberately not used. `autodoc_typehints` +#: below is built into `sphinx.ext.autodoc` and does the same job here, and the +#: extension calls an API Sphinx 10 removes — it emits a deprecation warning +#: per module on Sphinx 9. One less dependency for no loss. + +#: `REFACTOR_PLAN.md` is a working document, not documentation — it is written +#: for whoever is doing the refactor and records decisions and their evidence. +#: `notebooks/` is built by Phase 6 with myst-nb; until then the `.ipynb` files +#: would be copied in without being executed, which is worse than leaving them +#: out. `notes/` *is* included: `choosing_a_transform.md` links to it for a +#: per-trace table, so excluding it broke that link. +exclude_patterns = [ + "_build", + "REFACTOR_PLAN.md", + "notebooks/*", + "Thumbs.db", + ".DS_Store", +] + +#: Markdown only. Every page in `docs/` is already written that way, and +#: allowing both means two syntaxes for the same job. +source_suffix = {".md": "markdown"} + +#: `linkify` is deliberately absent: it needs `linkify-it-py` and every link in +#: these pages is already explicit. +myst_enable_extensions = [ + "colon_fence", # ::: fences, so a directive can hold a code block + "deflist", + "dollarmath", + "substitution", +] +#: Heading anchors down to h3, so `docs/*.md` can link to each other's sections. +myst_heading_anchors = 3 + +# ------------------------------------------------------------------ autodoc + +autodoc_typehints = "description" +autodoc_member_order = "bysource" +autodoc_default_options = { + "members": True, + "show-inheritance": True, + "member-order": "bysource", +} + +napoleon_google_docstring = False +napoleon_numpy_docstring = True + +intersphinx_mapping = { + "python": ("https://docs.python.org/3", None), + "numpy": ("https://numpy.org/doc/stable", None), + "scipy": ("https://docs.scipy.org/doc/scipy", None), + "pandas": ("https://pandas.pydata.org/docs", None), + "matplotlib": ("https://matplotlib.org/stable", None), + "obspy": ("https://docs.obspy.org", None), + "lmfit": ("https://lmfit.github.io/lmfit-py", None), +} +#: Resolving these needs the network. A build without it still produces a +#: site, with the cross-references left as plain text. +intersphinx_disabled_reftypes = ["*"] + +# --------------------------------------------------------------------- html + +html_theme = "pydata_sphinx_theme" +html_title = f"{project} {version}" +html_theme_options = { + "github_url": "https://github.com/sgjholt/SpecMod", + "show_prev_next": True, + "navigation_with_keys": False, +} +html_static_path: list[str] = [] diff --git a/docs/development.md b/docs/development.md new file mode 100644 index 0000000..719d9f3 --- /dev/null +++ b/docs/development.md @@ -0,0 +1,437 @@ +# Developer guide + +Everything needed to work on SpecMod: the tools, the loop, the conventions, +and — the part that is easiest to get wrong — where development stops and +releasing begins. + +This page is the hub. The three companions are +[Documentation workflow](documentation.md) (previewing, publishing, versions), +[Releasing the software](releasing.md) (tags, PyPI, DOI) and +[Publishing a dataset](releasing-data.md) (`data-v*` artefacts). + +**Where to look for what** + +| I want to… | Go to | +|---|---| +| Get a working checkout | [Quick start](#quick-start) | +| Know what lives where | [The repository, mapped](#the-repository-mapped) | +| Make a change | [The daily loop](#the-daily-loop) | +| Understand a tool or a check | [Tooling reference](#tooling-reference) | +| Write or fix a test | [Testing](#testing) | +| Read a CI failure | [What CI runs](#what-ci-runs) | +| Change a workflow file | [The `ci/` mirror](#the-ci-mirror) | +| Understand versions and releases | [Development versus release](#development-versus-release) | +| Preview or publish docs | [Documentation workflow](documentation.md) | +| Use Claude or Codex on this repo | [Working with agents](#working-with-agents) | +| Build a package on top of SpecMod | [The stable surface](#the-stable-surface) | + +## Quick start + +Requires Python 3.11+ and [uv](https://docs.astral.sh/uv/). + +```sh +git clone https://github.com/sgjholt/SpecMod.git +cd SpecMod +uv venv && uv pip install -e ".[dev]" +pre-commit install # both hook types; see below +pytest -m "not dataset and not notebook" # ~3 minutes, should be all green +``` + +Four things about that, in order of how often they bite: + +1. **`pre-commit install` is a required step, not a nicety.** It installs the + `pre-commit` *and* `commit-msg` hooks — the config names both stages under + `default_install_hook_types`, because a plain `pre-commit install` used to + wire only the first and left the commit-msg check inert. +2. **Install with `[dev]`, not bare.** The I/O suite needs `h5py` and + `pyarrow`, which `[dev]` pulls in; CI does the same. +3. **`-m "not dataset and not notebook"`** is what CI runs. `dataset` tests + need a network download and `notebook` executes the tutorial with a Jupyter + kernel (~40 s). +4. **The editable install must come from a git checkout with history.** The + version is derived by `hatch-vcs` from `git describe`; a tarball without + `.git` reports the fallback `0.0.0`. + +Optional extras, if you are working on those code paths: + +| Extra | Adds | Needed for | +|---|---|---| +| `multitaper` | [`multitaper`](https://github.com/gaprieto/multitaper) | Prieto's estimator, jackknife CIs | +| `wavelet` | [PyWavelets](https://pywavelets.readthedocs.io) | wavelet families beyond the built-in Morlet | +| `io` | h5py, pyarrow | HDF5 and Parquet persistence (in `[dev]` already) | +| `docs` | Sphinx and friends | building this site | +| `tutorial` | ipykernel, nbclient | executing the tutorial notebook | +| `mcmc` | [emcee](https://emcee.readthedocs.io) | sampling-based fits | + +## The repository, mapped + +``` +src/specmod/ the package + config/ layered settings, provenance stamping + core/ Spectrum, collections, noise, bandwidth, scalogram, units + transforms/ FFT, Welch, multitaper, Prieto, quadratic, CWT + smoothing/ Konno–Ohmachi, log binning + sources/ source models, attenuation, motion factors + fitting/ the fitter: base, guess, spectrum, event + picks/ pick readers, sensor resolution, the registry + acquire.py the only module that touches the network + datasets.py hash-pinned published datasets, via pooch + cli.py the `specmod` command +tests/ the suite, plus tests/golden/ reference numbers +tools/ repository scripts, each with a CI job or test behind it +ci/workflows/ staged copies of .github/workflows (see below) +docs/ this site + REFACTOR_PLAN.md the working document — not part of the built site +datasets/ dataset definitions for `specmod fetch` +tutorial/ the tutorial notebook and its data +stubs/ hand-written ObsPy type stubs +``` + +Two files that are not what they look like: + +- **`docs/REFACTOR_PLAN.md`** is a working document, deliberately excluded from + the built site. It records decisions, the measurements behind them, and an + audit (§6.6) of claims in it that turned out to be false. When something here + says "why", that is usually where the long answer is. +- **`ci/workflows/`** holds complete copies of the live GitHub Actions + workflows. See [The `ci/` mirror](#the-ci-mirror). + +## The daily loop + +```sh +git switch -c my-change # branch off main +# ... edit ... +pytest -m "not dataset and not notebook" # the suite CI runs +pytest --without-optional-extras # what a default install sees +mypy # strict, on the whole package +git commit # Conventional Commits; hooks run +git push -u origin my-change # open a PR against main +``` + +`ruff` runs automatically on commit via pre-commit; run it by hand with +`ruff check src/ tests/ tools/` and `ruff format src/ tests/ tools/`. + +**Run `--without-optional-extras` before pushing.** A development environment +with `specmod[multitaper]` installed passes tests that CI, which installs only +`[dev]`, fails. It has happened twice. + +### Branches + +- **`main`** is the trunk. Everything lands here, and it is the default branch. +- **`master`** is frozen: the permanent record of the pre-refactor code, doing + the job a `v0.1.0` tag would have done. Never commit to it. +- Feature branches are short-lived and merge into `main` via pull request. + +Every pull request must target **`sgjholt/SpecMod`**. See §6.7 of the plan for +why that is worth checking rather than assuming. + +### Commit messages + +[Conventional Commits](https://www.conventionalcommits.org), because +`release-please` parses them to decide the version and write the changelog — +see [Development versus release](#development-versus-release). The type +controls where the commit lands: + +| Type | Effect on the release | +|---|---| +| `feat:` | minor bump; **Features** | +| `fix:` | patch bump; **Bug Fixes** | +| `perf:` | patch bump; **Performance** | +| `refactor:`, `docs:`, `build:` | no bump; shown in the changelog | +| `test:`, `ci:`, `style:`, `chore:` | no bump; hidden | +| `feat!:`, or a `BREAKING CHANGE:` footer | **minor** bump while below 1.0 | + +There is no `commitlint` hook — the convention is followed by hand, and the +plan's §6.6 records that as a claim it once made falsely. + +**No session URLs from AI coding tools** in commit messages or anywhere else +published. The repository is public and those links are private state; a +`commit-msg` hook rejects them, which is the reason `pre-commit install` +appears in the quick start rather than further down. + +## Tooling reference + +| Tool | Runs | Configured in | +|---|---|---| +| [uv](https://docs.astral.sh/uv/) | environments, installs, builds | `pyproject.toml` | +| [ruff](https://docs.astral.sh/ruff/) | lint + format, pre-commit and CI | `[tool.ruff]` | +| [mypy](https://mypy.readthedocs.io) | strict types over `src/specmod` | `[tool.mypy]` | +| [pytest](https://docs.pytest.org) | the suite | `[tool.pytest.ini_options]` | +| [hypothesis](https://hypothesis.readthedocs.io) | property tests | in-test | +| [pre-commit](https://pre-commit.com) | hooks on commit | `.pre-commit-config.yaml` | +| [hatch-vcs](https://github.com/ofek/hatch-vcs) | version from git tags | `[tool.hatch.version]` | +| [Sphinx](https://www.sphinx-doc.org) + [MyST](https://myst-parser.readthedocs.io) | this site | `docs/conf.py` | +| [release-please](https://github.com/googleapis/release-please) | changelog and version decision | `release-please-config.json` | + +### The `tools/` scripts + +Each one exists because something was claimed and not enforced. All are +standard-library-only unless noted. + +| Script | Does | Enforced by | +|---|---|---| +| `check_ci_mirror.py` | staged workflows match the live ones | `lint` job | +| `check_floors.py` | the installed versions really are the declared minimums | `floors` job | +| `check_built_version.py` | the built wheel's version is the tag | `publish` job | +| `make_golden.py` | regenerates `tests/golden/*.json` | run by hand, deliberately | +| `measure_docs.py` | regenerates the measured tables in `docs/` | `tests/test_docs_are_current.py` | + +`measure_docs.py` is worth knowing about before editing a table by hand: + +```sh +python tools/measure_docs.py show # print the tables +python tools/measure_docs.py write # refresh the docs in place +python tools/measure_docs.py check # fail if any table is stale +``` + +Numbers that came from a measurement live between markers and are generated. +`tests/test_docs_are_current.py` runs `check`, so a change that moves a +published number fails the suite instead of leaving the prose quietly wrong. +The `--field` measurements read `tutorial/data/events/` and are slower; +refresh those by hand after changing an estimator. + +### Configuration and provenance + +Settings live in `src/specmod/config/` as semantic sections with layered +overrides, not module-level constants read at import time. Two commands: + +```sh +specmod config show # the resolved configuration, with its layers +specmod config freeze # write it out, pinned +``` + +Every output records the configuration that produced it, a hash of it, and the +SpecMod version. That is what makes a locally-overridden run reproducible from +its own outputs, and it is the mechanism that lets the package stay alpha +without published results becoming unrepeatable. + +## Testing + +```sh +pytest -m "not dataset and not notebook" # what CI runs +pytest --without-optional-extras # as a default install sees it +pytest -m notebook # executes the tutorial (~40 s) +pytest -m dataset # needs a network download +pytest tests/test_transforms.py -q # one module +``` + +Markers are declared in `pyproject.toml` and `--strict-markers` is on, so a +typo in a marker name is an error rather than a silently-skipped filter. + +**The tiers**, as the plan lays them out: + +1. **Property tests** (`hypothesis`) encoding the physics — Parseval, scaling, + units — which hold for any input rather than one recorded case. +2. **Synthetic end-to-end**: generate a spectrum with known source parameters, + run the whole pipeline, recover them. +3. **Golden/regression**: run the current code on the tutorial event and on + Magna, and compare against committed summaries in `tests/golden/`. +4. **Unit tests** for specific bugs, each written as a failing test first. + +### Golden references, and the one rule about them + +`tests/golden/*.json` records what this code produced at a known-good point. +It is compared as a **distributional summary** — median, quantile profile, +length — with a relative tolerance, not as a byte digest: an earlier version +hashed the raw float64 bytes and failed on every CI runner, because a +different numpy or BLAS build produces last-bit differences on identical +input. A reference that only holds on the machine that generated it is not a +reference. + +**Do not regenerate it to make a test pass.** If a change moves a number, that +is the finding: say which number, by how much, and why. Then regenerate +deliberately: + +```sh +python tools/make_golden.py # and commit the result, with the reason +``` + +Tolerances carry comments explaining what was measured to choose them. The +`cwt` entry is the worked example — it records a per-runner residual that is +bounded rather than explained, and says so in as many words. + +## What CI runs + +Five jobs in `test.yml`, plus two more workflows. All of them run on every pull +request. + +| Job | Workflow | Does | +|---|---|---| +| `lint` | `test.yml` | `ruff check`, `ruff format --check`, and the `ci/` mirror check | +| `typecheck` | `test.yml` | `mypy`, strict, over the whole package | +| `test` | `test.yml` | pytest on 3.11/3.12/3.13 × ubuntu/macOS, coverage to Codecov from one cell | +| `floors` | `test.yml` | installs the *declared minimum* versions and runs the suite | +| `notebook` | `test.yml` | executes the tutorial notebook | +| `build` | `build.yml` | sdist + wheel, `twine check`, install-from-wheel smoke test | +| `docs` | `docs.yml` | builds the site as a check; publishing is Read the Docs' job | +| `release` | `release.yml` | the release PR, and publishing — see below | + +Two of those are worth understanding before you read a failure from them: + +**`floors`** installs `--resolution lowest-direct`, exercising the oldest +dependency set the project claims to support. It caught two floors that could +never have worked: `lmfit>=1.2` with `numpy>=2.0` (lmfit below 1.3 calls +`np.asfarray`, removed in NumPy 2), and `scipy>=1.13` silently breaking the +quadratic multitaper. If you raise or add a dependency, this is the job that +tells you whether the floor you wrote is real. + +**`docs`** deliberately does **not** use `-W`. Intersphinx resolves seven +inventories over the network and warns whenever one is briefly unreachable; +turning a third party's downtime into a red build is flakiness, not a check. + +## The `ci/` mirror + +`ci/workflows/*.yml` holds complete, ready-to-paste copies of +`.github/workflows/*.yml`. The reason is narrow: the GitHub App token used by +AI coding sessions has no `workflows` permission, so a push touching +`.github/workflows/` is rejected. Rather than describe an edit in a comment and +hope it is applied correctly, the intended file is committed in full. + +To apply one, copy the whole file over its counterpart — no merging, no partial +application. `tools/check_ci_mirror.py` runs in the `lint` job and fails while +a pair differs, so a staged change that has not been copied across shows up as +a red build rather than being forgotten. **That failure is the reminder.** It +clears the moment the files match. + +The mirror is the *intended* state, which is not always the current one — in +either direction. If you fix a workflow through the GitHub web editor, copy it +back into `ci/` so the next person staging a change starts from the working +version. + +Full detail in [`ci/README.md`](https://github.com/sgjholt/SpecMod/blob/main/ci/README.md). + +## Development versus release + +The thing to hold onto: **merging to `main` releases nothing.** `main` is +continuously integrated and continuously *documented*, but it is not +continuously published. Three separate clocks: + +| | What moves it | What it produces | Who decides | +|---|---|---|---| +| **Development** | any merge to `main` | updated `main`, updated docs site | whoever merges the PR | +| **Software release** | merging the release PR | a `v*` tag, a GitHub Release, a PyPI upload, a DOI | a human, deliberately | +| **Data release** | creating a `data-v*` tag by hand | a dataset artefact pinned by hash | a human, deliberately | + +### How a version comes to exist + +No version string is committed anywhere. `hatch-vcs` derives it from +`git describe`, so **the tag is the version**: + +- On `main` between releases you get `.postN.devN` — a version that + claims nothing. +- On a `v*` tag you get exactly that tag without its `v`. + +`pyproject.toml` constrains which tags count, with both a `--match v[0-9]*` on +the describe command and a `tag_regex` on the parse. Both are needed: a +`data-v1` tag was measured to take the package version from +`0.1.0.post1.dev173` to `1`, because setuptools-scm's default pattern strips +the `data-` prefix and reads what is left. + +That is why the two release channels use different tag prefixes, and why +`tests/test_versioning.py` pins the parse. + +### What a release actually does + +`release-please` watches `main` and keeps a standing pull request titled +`chore(main): release `, carrying the generated `CHANGELOG.md`. +Merging it creates the tag and the GitHub Release; a gated job then builds from +the tag, checks the built version against it, and uploads to PyPI via Trusted +Publishing; Zenodo mints a DOI from the release webhook. + +Nothing publishes until that merge. The gate exists because **a DOI cannot be +retracted** — fully automatic tagging plus Zenodo means a typo fix can mint a +citable version of the software. + +The step-by-step, including the six repository settings that have to be turned +on once and cannot be expressed in a commit, is in +[Releasing the software](releasing.md). + +### What this means day to day + +- **Land work whenever it is ready.** The changelog accrues; you are not + choosing a version when you merge. +- **Write the commit message for the changelog**, because that is where it ends + up verbatim. +- **Breaking changes are allowed** and land in minor bumps while the project is + `0.x`. There is no deprecation cycle, by design — see the + [roadmap](roadmap.md) for what 1.0 will change about that. +- **Datasets are versioned by registry name**, not by the package version: + `magna_2020_v1` and `magna_2020_v2` are separate entries, so a published + result pinned to v1 keeps fetching v1 forever. Nothing revises an entry in + place. + +## The stable surface + +`specmod.api` is a small re-export module, and the only part of SpecMod that +carries a compatibility promise: one minor release of `DeprecationWarning` +before anything on it is removed or changes signature, even while the package +is `0.x`. Everything else may move in any release. + +It exists for downstream packages. SpecMod's internals are still being +refactored; a package that imports `specmod.core` or `specmod.fitting` directly +takes on that churn, and one that imports `specmod.api` does not. + +Five properties hold across the surface, and they are enforced by +`tests/test_api_surface.py` rather than promised in prose: **path-free** +(nothing on it opens a file — a consumer that stores its data on S3 has to be +able to hand in arrays), **deterministic**, **non-mutating**, **quiet** (no +`print`), and **typed errors** rooted at `SpecModError`, distinguishing a +caller's bad input from a missing optional backend from an internal bug. + +Adding an export is a compatibility obligation, and the procedure is in +[`CONTRIBUTING.md`](https://github.com/sgjholt/SpecMod/blob/main/CONTRIBUTING.md). + +The audit that established what could go on it — path coupling, hidden state, +determinism, and what one multitaper estimate actually costs — is in +[Audit: what `specmod.api` found in core](notes/api_audit.md). Two of its +findings were defects in core rather than in the surface, both since fixed and +both now guarded package-wide by `tests/test_ambient_state.py`: **no module +reads configuration at import time**, and **no module prints**. Those two +properties are worth knowing before adding code — a config read at module +level freezes the working directory the process started in, and a `print` is +invisible to a caller capturing logs. + +```{toctree} +:hidden: + +notes/api_audit +``` + +## Working with agents + +Claude Code, Codex and similar tools are used on this repository. The durable +rules live in [`AGENTS.md`](https://github.com/sgjholt/SpecMod/blob/main/AGENTS.md) +at the repository root, with `CLAUDE.md` pointing at it so there is one copy; +agents read those automatically. What follows is the context for a human +supervising one. + +**The failure modes are environmental, not intellectual.** Every one of these +has happened here: + +- **A fresh container has no git hooks.** `pre-commit install` has to be run in + the session, or the `commit-msg` check that rejects session links is simply + not there. Three commits went out with session trailers before this was + noticed; the config now installs both hook types from one command, and the + first thing `AGENTS.md` says is to run it. +- **An agent's token cannot push `.github/workflows/`.** This is what the + `ci/` mirror exists for. An agent that does not know about it will either + fail the push or, worse, quietly drop the change. +- **A development container often lacks the optional extras**, so an agent's + green run can be greener than CI's. `--without-optional-extras` is the check. +- **A harness may append its own commit trailers.** The repository's rules take + precedence over a tool's defaults, and this one is a publishing rule rather + than a style preference. + +**What to ask for in review.** The habit this repository is built around is +saying what was checked and what was not. An agent that reports "fixed" should +be able to show the command and its output; one that widens a tolerance or +regenerates a golden file to reach green has moved the goalposts rather than +found the problem. §6.6 of the plan is an audit of exactly that failure — three +claims stated as fact with no mechanism behind them — and it is worth reading +once before delegating anything that touches a check. + +**What agents are good at here.** The mechanical, checkable work: splitting +modules while keeping public names, writing the test that pins a behaviour +before changing it, running the same verification five ways, and the tedious +correctness of the docs — which is the same skill as the tests, since numbers +in prose go stale silently. diff --git a/docs/documentation.md b/docs/documentation.md new file mode 100644 index 0000000..7027ef3 --- /dev/null +++ b/docs/documentation.md @@ -0,0 +1,211 @@ +# Documentation workflow + +How this site is written, previewed, checked and published — and what the +difference is between a build you made to look at, the build CI makes on a +pull request, and the versions people actually read. + +Companion to the [Developer guide](development.md). + +## Where each build goes + +| Build | Made by | Lives | Official? | +|---|---|---|---| +| **Local** | you, `sphinx-build` | `docs/_build/html/`, gitignored | no — nobody else can see it | +| **Pull request check** | the `docs` job in GitHub Actions | a downloadable artefact on the run | no — it checks, it does not publish | +| **Pull request preview** | Read the Docs | its own temporary URL, linked from the PR | no — it disappears when the PR closes | +| **`latest`** | Read the Docs, from `main` | `/en/latest/` | yes, but it is the trunk | +| **`stable`** | Read the Docs, from the newest `v*` tag | `/en/stable/`, and the default | yes — this is the site | +| **`v0.2.0`, `v0.3.0`, …** | Read the Docs, from each tag | `/en/v0.2.0/` | yes, and frozen | + +**Read the Docs publishes; GitHub Actions only checks.** That split is +deliberate, and the reason is versions. This package is alpha and the docs tell +people to pin an exact version, so the documentation for a release has to stay +readable after the trunk moves on. One site that always shows `main` — which is +what deploying to GitHub Pages gives you — means someone pinned to `0.2.0` +reads about code they do not have. + +Merging to `main` updates `latest` within a few minutes and changes `stable` +not at all. `stable` moves when a release is tagged. See +[Development versus release](development.md#development-versus-release) for why +those are separate clocks. + +## Building it locally + +```sh +uv pip install -e '.[docs]' +sphinx-build -b html docs docs/_build/html +python -m http.server -d docs/_build/html 8000 # then open localhost:8000 +``` + +Rebuild after editing; Sphinx is incremental, so only changed pages are redone. +To force a clean build — worth doing before trusting a warning count — +`rm -rf docs/_build` first. + +**Expect warnings about intersphinx.** The build resolves seven inventories +(python, numpy, scipy, pandas, matplotlib, obspy, lmfit) over the network, and +each unreachable one is a warning. Offline you will see seven of them and a +site that builds correctly, with cross-references to other projects left as +plain text. That is also why neither the CI job nor `.readthedocs.yaml` fails +on warnings: turning a third party's downtime into a red build is flakiness +rather than a check. + +Anything *other* than those seven is a real warning. A clean build today looks +like: + +``` +build succeeded, 7 warnings. # all of them intersphinx, offline +``` + +## The two automatic builds on a pull request + +**Read the Docs builds a preview** at its own URL and posts the link as a +status check on the pull request. That is the one to click when you want to +*look* at a change. It is torn down when the pull request closes. + +**The `docs` job in GitHub Actions builds the same site and keeps the HTML** +as an artefact. It is kept alongside the Read the Docs preview for two reasons: +it does not depend on a third-party service being up, and because autodoc +imports every module it documents, it doubles as an import check. A module that +no longer imports fails here. + +To open what Actions built: the pull request → **Checks** → **docs** → the run +summary → download the **docs-html** artefact, then + +```sh +unzip docs-html.zip -d preview +python -m http.server -d preview 8000 +``` + +## Read the Docs setup + +One-time, at the point the first version is published. None of it can be done +from a commit; all of it is on . + +1. **Import the project.** Log in with GitHub, *Import a Project*, pick + `sgjholt/SpecMod`. The slug chosen here becomes the URL — + `https://.readthedocs.io` — so it is worth getting right first time. +2. **Check it found the config.** `.readthedocs.yaml` is in the repository + root and is read automatically; the first build's log says which + configuration file it used. If it did not find one, Read the Docs falls back + to defaults that will not install this package. +3. **Turn on pull request builds.** *Admin → Settings → Advanced settings → + Build pull requests for this project*. Without it there are no previews. +4. **Add an automation rule for tags.** *Admin → Automation Rules → Add rule*, + version type **Tag**, action **Activate version**, matching `^v.*`. New + releases then publish themselves; without a rule each tag has to be + activated by hand before anyone can read it. +5. **Set the default version to `stable`** once a tag exists. *Admin → + Settings → Default version*. Until then it is `latest`, which is correct + while there are no releases. + +The build configuration itself is +[`.readthedocs.yaml`](https://github.com/sgjholt/SpecMod/blob/main/.readthedocs.yaml). +Two things in it are load-bearing and easy to lose: + +- **`post_checkout` unshallows the clone and fetches tags.** Read the Docs + clones shallow to save time, and `hatch-vcs` derives the version from + `git describe`. Without those two lines every build — including a tagged one + — reports the fallback `0.0.0` in the sidebar. +- **`extra_requirements: [docs, io]`.** The `io` extra is there because autodoc + imports `specmod.io`, which imports h5py and pyarrow. The CI job installs the + same pair for the same reason. + +## Versions, and what to link to + +- **Link to `/en/stable/`** in papers, READMEs and anywhere durable. It follows + releases without becoming stale. +- **Link to `/en/v0.2.0/`** when the reference is to a specific version's + behaviour — for a published result, this is the honest link. +- **`/en/latest/` is the trunk**, ahead of any release during alpha. It is what + contributors should read and what nobody should cite. + +A version's build is frozen at what it said when that tag was cut, which is the +point: rebuilding 0.2's documentation later would not make it more true. + +## Writing for the site + +### What belongs here, and what does not + +The site is for someone *using or developing* SpecMod. `REFACTOR_PLAN.md` is +excluded from it on purpose: it is a working document, written for whoever is +doing the refactor, recording decisions and the evidence behind them. Link to +it on GitHub rather than adding it to the build. + +`docs/notes/` **is** included, which was a correction: `choosing_a_transform.md` +links to `notes/window_position.md` for a per-trace table, and a page another +page depends on is documentation whatever its folder is called. + +### Adding a page + +1. Write `docs/.md` in MyST Markdown. +2. Add it to the `toctree` at the bottom of `docs/index.md`, and usually to the + "Where to start" list above it. +3. Build. A page in no toctree builds but warns, and is reachable only by a + direct link. + +If a page is a supporting note rather than a top-level one, put it in a hidden +toctree on the page that cites it — that is how `notes/window_position.md` is +attached to `choosing_a_transform.md`: + +````markdown +```{toctree} +:hidden: + +notes/window_position +``` +```` + +### Markdown, math and directives + +[MyST](https://myst-parser.readthedocs.io) with `colon_fence`, `deflist`, +`dollarmath` and `substitution` enabled, and heading anchors down to `h3` so +pages can link to each other's sections. + +- Inline math is `$...$`, display math `$$...$$`. `dollarmath` is what makes + those work; without it MyST does not read `$` at all. +- `amsmath` is **not** enabled — nothing here uses a bare `\begin{align}` + outside `$` delimiters. A `\begin{cases}` inside `$$` renders without it. +- Markdown only. `source_suffix` maps `.md`; there is no reStructuredText in + `docs/`, so there is one syntax rather than two. +- For a Sphinx directive that has no MyST spelling, drop into rST: + +````markdown +```{eval-rst} +.. automodule:: specmod.picks + :imported-members: +``` +```` + +### The API reference + +`docs/api.md` is hand-ordered — grouped by what a run does in order, rather +than alphabetically — and each entry is an `automodule`. Two rules learned from +the first build: + +- **Document a package at the path you import from.** `specmod.picks.PickSet`, + not `specmod.picks.base.PickSet`. Documenting both the package and its + submodules gives every re-exported name two targets and makes every + cross-reference to it ambiguous. +- **Some objects break autodoc's signature formatter.** `HOLT_2019_UTAH` is a + callable dataclass instance documented as module data; `inspect.signature` + handles it, autodoc's own formatter raises. It is excluded with + `:exclude-members:` and written out in prose instead. + +Type hints come from the annotations via `autodoc_typehints = "description"`. +`sphinx-autodoc-typehints` is deliberately **not** used: measured, it produced +the same 367 documented objects while calling an API Sphinx 10 removes. + +### Numbers in prose + +Any table that came from a measurement is generated, not typed. Edit +`tools/measure_docs.py` and run `python tools/measure_docs.py write`; +`tests/test_docs_are_current.py` fails if a table drifts from what the code +does. See the [Developer guide](development.md#the-tools-scripts). + +## Changing the docs workflow + +`docs.yml` cannot be pushed by an AI coding session — no `workflows` token +permission — so the intended file lives at `ci/workflows/docs.yml` and is +copied across by hand. The `lint` job fails while the two differ, in **either** +direction: if you fix the live workflow in the web editor, copy it back into +`ci/`. See [The `ci/` mirror](development.md#the-ci-mirror). diff --git a/docs/index.md b/docs/index.md new file mode 100644 index 0000000..261b790 --- /dev/null +++ b/docs/index.md @@ -0,0 +1,109 @@ +# SpecMod + +A Python toolbox for processing and modelling seismic spectra: cut a window, +estimate its spectrum, decide which part of it is above the noise, and fit a +source model to that band. + +:::{warning} +**Alpha. Pre-1.0, and mid-refactor.** SpecMod is being rebuilt in the open, so +treat everything here as provisional until the API settles at 1.0: + +- **Names and signatures move between `0.x` releases**, without a deprecation + cycle. Breaking changes land in minor bumps by design — that is what `0.x` + is for, and the alternative is a deprecation shim on an API that is still + being worked out. +- **Some numbers still move too.** The modern layers (`specmod.config`, + `specmod.core`, `specmod.transforms`, `specmod.picks`, `specmod.fitting`) + are built and tested against golden references; the modules that have not + been reached yet carry pre-refactor behaviour, and a fix there can change a + result. Changes that move a published number are called out in the + changelog. +- **Pin an exact version for anything you intend to publish**, and keep the + configuration stamp that every output carries. Together they are what make a + run reproducible while the package underneath is still moving. + +What will *not* change silently: the units conventions and the Parseval +contract are pinned by tests, and the golden references fail loudly rather +than drifting. The [roadmap](roadmap.md) says which stages are done and what +1.0 will mean. +::: + +```python +import specmod.preprocess as pre +from specmod.pipeline import spectrum_set_from_streams +from specmod.fitting import FitSpectra + +pre.set_picks(stream, "event.xml") +signal = pre.get_signal(stream, pre.cut_s, rafp=0.8, tafs=20) +noise = pre.get_noise_p(stream, signal) + +spectra = spectrum_set_from_streams(signal, noise) +fits = FitSpectra(spectra) +fits.fit_spectra() +print(fits.table[["id", "llpsp", "fc", "ts"]]) +``` + +## Where to start + +[Processing](processing.md) +: Every step of the pipeline with the equation it implements — what a window + is, how the noise is compared against it, and what the bandwidth selector + does. + +[Choosing a transform](choosing_a_transform.md) +: What each estimator does to your data, measured. The choices here change + recovered amplitude by factors of three on real windows, which is about 0.3 + magnitude units. + +[Reading picks](pick-formats.md) +: What arrival formats are read out of the box, how to add one, and how a pick + is matched to a trace. + +[Publishing a dataset](releasing-data.md) +: Taking an event from an FDSN archive to a hash-pinned entry in the registry. + +[Releasing the software](releasing.md) +: How a merged commit becomes a tag, a PyPI release and a DOI, and the six + settings that have to be turned on once. + +[Roadmap](roadmap.md) +: What is built, what is being worked on, and what 1.0 will mean. Stages, not + dates. + +## Working on SpecMod + +[Developer guide](development.md) +: Quick start, the repository mapped, the daily loop, every tool and check, and + where development stops and releasing begins. + +[Documentation workflow](documentation.md) +: Building the site locally, previewing a pull request's build, what the + official site is, and how versions would work. + +## Two things worth knowing early + +**Units are typed.** A spectrum carries its ground-motion domain and amplitude +convention as attributes. Converting between them is a method that returns a +new spectrum, so a moment computed from the wrong domain is a type error rather +than a wrong number. + +**Fit in the units the sensor recorded.** Integrating to displacement +implicitly low-passes and differentiating to acceleration amplifies +high-frequency noise, so the record to fit is the one that was measured. The +model carries a motion factor, so the plateau it reports is the displacement +one either way. + +```{toctree} +:maxdepth: 2 +:hidden: + +processing +choosing_a_transform +pick-formats +development +documentation +releasing-data +releasing +roadmap +api +``` diff --git a/docs/notes/api_audit.md b/docs/notes/api_audit.md new file mode 100644 index 0000000..35c3460 --- /dev/null +++ b/docs/notes/api_audit.md @@ -0,0 +1,201 @@ +# Audit: what `specmod.api` found in core + +Answers to the five questions that had to be settled before writing +`specmod.api`, each measured rather than reasoned about. The surface itself is +small; this was the work. + +It is kept because the properties it audits — path coupling, hidden state, +determinism — are the ones core's own provenance claims rest on, whether or not +anything downstream ever consumes them. + +## 1. Joint per-event inversion, or per-spectrum fitting? + +**Per-spectrum only.** `FitSpectra` is a loop, not a joint solver: + +```python +for name, mod in self.models.items(): + mod.fit_mod(**kwargs) +``` + +Each station gets its own `FitSpectrum` with its own parameters. Nothing is +shared between them — no common `t*`, no common Ω₀, no event-level term. The +"event fit" is an aggregation of independent fits into one table. + +A joint inversion is therefore not in core and `specmod.api` exposes the +per-spectrum primitive, `fit_spectrum`, plus everything a joint solver needs to +build its own problem: the per-bin SNR, the selected band, the model, and the +covariance of each single-station fit. + +## 2. Per-bin SNR, or a scalar bandwidth? + +**Per-bin, already.** `SpectrumPair.snr` is an array aligned with +`binned_signal.freq`, computed as the element-wise ratio of the two binned +spectra: + +```python +snr = binned_signal.amp / binned_noise.amp +band = find_bandwidth(binned_signal.freq, snr, threshold, method=bandwidth) +``` + +The scalar `band` is *derived* from the curve, and both are kept. So the thing +that cannot be un-collapsed later was never collapsed: a consumer can apply its +own threshold, or admit data bin by bin instead of over a contiguous interval, +without a reprocess. + +Measured on one synthetic station: 75 bins of `snr` against a `band` of +(6.47, 9.15) Hz. `compare_spectra` returns the pair whole — the curve, the +binned noise spectrum, and the resolution floor — rather than a summary. + +**What core does not have** is a `valid_mask`. Bins are excluded by falling +outside the band or below the resolution floor, and there is no per-bin defect +flag. That is a consumer-side concept and building it does not need core. + +## 3. Which functions take paths with no in-memory form? + +Listed rather than changed. None of them is on `specmod.api`, and the +estimation and fitting paths were already path-free — `SpectralEstimator.estimate` +takes `(data, dt)`. + +| Function | Takes | In-memory form today | +|---|---|---| +| `preprocess.read_picks` | path | none — wraps `picks.read` | +| `preprocess.set_picks` | path | none, though `picks.read` accepts a `Catalog` | +| `preprocess.rstfl` | paths | none | +| `picks.read` | path or `Catalog` | **yes**, `Catalog` | +| `picks.detect_reader`, every `PickReader.read` | path | none — they sniff the file | +| `tables.read_table` / `write_table` | path | none | +| `io.*` | path | none | +| `datasets.*`, `acquire.*` | paths, network | not applicable — that *is* their job | + +The one worth fixing first is `set_picks`: `picks.read` already accepts an +in-memory `Catalog`, so the path-free form exists one layer down and is not +plumbed through. A caller holding a `Catalog` has to reach past `preprocess` to +use it. + +## 4. Hidden global or module-level state + +Four instances. Three are benign; the fourth was not, and is fixed. + +**A module-level config read, at import time — fixed.** `fitting/base.py` +had: + +```python +PLOT_COLUMNS = cfg.load_config().config.viz.plot_columns +``` + +`load_config()` with no arguments resolves against the *current working +directory* and the environment. At module level that answer is frozen for the +life of the interpreter, so two jobs in one worker with different project +directories both get the first one's value. Reproduced before fixing: +importing from a project whose `specmod.toml` said `plot_columns = 5`, then +moving to one that resolves to 3, left the constant at 5. + +It is now `fitting.plot_columns()`, resolved per call; the old name still +imports and emits a `DeprecationWarning` naming the replacement and the +release it goes in. `tests/test_ambient_state.py` parses **every** module in +the package and fails on a module-level `load_config()` anywhere, because this +is the kind of defect that comes back one file at a time. + +That test found nothing else — but only after its own first version was +wrong. It used `ast.walk`, which descends into method bodies, and flagged +three modules that read configuration perfectly properly at call time. The +walk now stops at a function body while still checking default arguments and +decorators, which *do* run at import and are where this would hide next. + +**Implicit config reads at call time**, in twelve places including +`fitting/event.py`, `fitting/guess.py`, `fitting/spectrum.py`, `pipeline.py` +and `sources/composite.py`. Not import-time, so not frozen, but still +working-directory-dependent. `specmod.api` closes this where it can — every +estimation and comparison argument is explicit — and documents it where it +cannot: `fit_spectrum` reads `[fitting]` for the minimiser and the initial +guess, and its docstring says so. + +**The pick-reader registry.** `PICK_READERS` is a module-level dict mutated by +`register_reader` and by entry-point discovery, guarded by a `_plugins_loaded` +flag. Which readers exist is a property of what is installed, and it is not +recorded in provenance. It affects reading, never numbers, and nothing on +`specmod.api` touches it. + +**No randomness anywhere.** Nothing in `src/specmod` imports `random` or +`numpy.random`, and nothing seeds. The `emcee` extra is declared but unused, so +the first sampler added is the moment an explicit `seed` argument has to be +required rather than recommended. + +### Two more things the audit turned up + +Not asked for, but found while looking, and both affect a consumer: + +**Nine `print()` calls — fixed** in `fitting/event.py` and `utils.py`, on +paths a caller reaches: an unrecognised weight method, a station that failed +to fit, a missing id. A service capturing logs per job got nothing from them, +and a CLI writing to a pipe got its output corrupted. + +Each became a `warnings.warn` or a module-logger call, and which one is not a +matter of taste. **`warnings` deduplicates by code location**, so the +per-station failure inside `fit_spectra`'s loop would report the first station +and silently drop the rest — and a station that could not be fitted is missing +from the results, so that is precisely the line that must not collapse. That +site logs; the caller-actionable ones warn. Nothing calls `logging.basicConfig` +anywhere in the package: a library that configures logging decides formatting +and destination for its host process. + +The same package-wide test asserts no module calls `print` at all. The count +is now zero across `src/specmod`, `cli.py` included — it uses `click.echo`. + +**Uncertainty depends on the minimiser, and the default provides none.** +`[fitting] method` ships as `powell`, which estimates no covariance matrix, so +every parameter's `stderr` is `None` and there is no correlation to report. +Measured on one synthetic station: + +| `method` | `fc` | `fc` error | `fc`–`t*` correlation | +|---|---|---|---| +| `powell` (default) | 7.925 | — | — | +| `nelder` | 7.925 | — | — | +| `leastsq` | 7.925 | 0.129 | 0.837 | +| `least_squares` | 7.925 | 0.129 | 0.837 | + +All four agree on the point estimate to three decimals against a true 8.0 Hz. +Only the least-squares family answers "how well". `SpectrumFit` reports the +absence as an absence — empty `stderr`, `covariance=None`, +`correlation()` returning `None` rather than zero — because a zero error is a +claim, and the wrong one. The 0.84 correlation is the reason neither `fc` nor +`t*` should be quoted alone. + +## 5. What does one multitaper estimate cost? + +Measured on the 28 real S windows of the Preston New Road event, through +`MultitaperEstimator` (DPSS, `scipy.signal.windows.dpss`). + +**Re-estimation on a window change** — the full path: estimate the signal, +estimate the noise, then the Parseval rescale, interpolation onto a common +axis, log binning, the per-bin SNR and the band search. + +| Operation | Median | p95 | +|---|---|---| +| `multitaper` estimate, one window | 3.05 ms | 3.31 ms | +| `fft` estimate, one window | 0.12 ms | 0.20 ms | +| **multitaper ×2 + full compare** | **7.02 ms** | 7.55 ms | +| `fft` ×2 + full compare | 2.32 ms | 2.71 ms | + +Scaling is close to linear in window length, not quadratic: + +| Samples | Duration | Median | +|---|---|---| +| 512 | 2.6 s | 1.92 ms | +| 1024 | 5.1 s | 3.63 ms | +| 2048 | 10.2 s | 6.31 ms | +| 4096 | 20.5 s | 11.54 ms | +| 8192 | 41.0 s | 20.85 ms | + +So a 20-second window re-estimates in roughly 25 ms end to end, and a 3.7-second +one in 7 ms. Whatever budget a live window editor has, this is not what spends +it. + +**Configuration.** One machine: x86_64, 4 cores, Python 3.11.15, numpy 2.4.6, +scipy 1.17.1, obspy 1.5.0, local disk, records read from the repository's own +test data. Record under test: `UR.AQ06.00.HHN`, 737 samples at 200 Hz. + +**This is one configuration, and a shared cloud container at that.** Anything +that has to hold on a minimum-spec target or against remote object storage +needs measuring there; the numbers above answer "is this milliseconds or +seconds", which is the question that was blocking, and nothing more. diff --git a/docs/processing.md b/docs/processing.md index 8e78eea..5960bea 100644 --- a/docs/processing.md +++ b/docs/processing.md @@ -312,7 +312,7 @@ fixed increment and stopped at the first trial past the touch, with exactly the irreproducibility described above. Two departures from the legacy implementation are recorded in -[REFACTOR_PLAN §4.5.3](REFACTOR_PLAN.md#453-rotate-ported-the-noise-registry-is-complete): +[REFACTOR_PLAN §4.5.3](https://github.com/sgjholt/SpecMod/blob/main/docs/REFACTOR_PLAN.md): the solved angle, and taking the low/high split from the signal rather than the noise. Neither can move a published number — `ROT_METHOD = 1` was commented out on `master` and has never produced one. diff --git a/docs/releasing.md b/docs/releasing.md new file mode 100644 index 0000000..ec428b9 --- /dev/null +++ b/docs/releasing.md @@ -0,0 +1,160 @@ +# Releasing the software + +How a merged commit becomes a tag, a PyPI release and a DOI. The companion to +[Publishing a dataset](releasing-data.md), which covers the `data-v*` tags +instead — the two are deliberately separate, and `pyproject.toml` enforces the +separation. + +Nothing here is triggered by hand. What a human does is **merge the release +pull request**, and everything after that is automatic. That gate exists +because of Zenodo: every GitHub Release mints a DOI, and a DOI cannot be +retracted. A typo fix should not be able to mint a citable version of the +software without someone deciding it should. + +## What the pieces are + +| Piece | Decides | +|---|---| +| Conventional Commit messages | which section of the changelog a commit lands in, and whether it bumps | +| `release-please-config.json` | the bump rules, the changelog sections, the tag format | +| `.release-please-manifest.json` | the version we are on now | +| `.github/workflows/release.yml` | opens and maintains the release PR; publishes when it merges | +| `pyproject.toml` (`tag_regex`) | which tags `hatch-vcs` will read as a version | + +No version string is committed anywhere. `hatch-vcs` derives it from +`git describe`, so the tag *is* the version and there is nothing to forget to +bump. + +## One-time setup + +Six things, none of which can be done from a commit. Until they are done the +release workflow is inert rather than wrong — it opens a release PR and stops. + +1. **Copy the workflow into place.** `ci/workflows/release.yml` → + `.github/workflows/release.yml`. See [`ci/README.md`](https://github.com/sgjholt/SpecMod/blob/main/ci/README.md) + for why the file is staged rather than pushed. + +2. **Let Actions open pull requests.** Settings → Actions → General → + Workflow permissions → tick *Allow GitHub Actions to create and approve pull + requests*. Without it release-please fails with `GitHub Actions is not + permitted to create or approve pull requests`. + +3. **Create the `pypi` environment.** Settings → Environments → New + environment → `pypi`. Add required reviewers here if you want a second + human gate on the upload itself. + +4. **Register the trusted publisher on PyPI.** On the project page (or, for + the first ever upload, under *Publishing* → *Add a new pending publisher*): + + | Field | Value | + |---|---| + | Owner | `sgjholt` | + | Repository | `SpecMod` | + | Workflow | `release.yml` | + | Environment | `pypi` | + + The workflow name is matched against the OIDC token, so it must be the file + that actually runs. If the workflow is ever renamed, this has to be changed + in the same sitting or the upload fails authentication. + +5. **Turn on the Zenodo webhook.** Log into Zenodo with GitHub, find `SpecMod` + in the repository list, and flip the switch. It takes effect for releases + made *after* it is on, so do it before the first release rather than after. + `CITATION.cff` supplies the metadata. + +6. **Require the checks on `main`.** Branch protection → require the status + checks from `test.yml`, `docs.yml` and `build.yml`. The release pull request + is an ordinary pull request and goes through them like any other. + +## What a version number means while this is 0.x + +SemVer, with the pre-1.0 rule taken literally: **a breaking change bumps the +minor**, and there is no deprecation cycle. The package is alpha, the API is +still being worked out, and shipping shims for names that are about to move +again would cost more than it protects. + +So `0.2.0 → 0.3.0` may rename things, and `0.2.0 → 0.2.1` will not. Anyone +depending on this for published work should pin an exact version; the +configuration stamp on every output covers the rest. + +`1.0.0` is the release that says the API has stopped moving, which is why +`bump-minor-pre-major` exists — it must be a decision, not something a `feat!:` +does on its way past. + +## The cycle + +1. Merge work into `main` with Conventional Commit messages. Nothing is + released. +2. `release.yml` opens (or updates) a pull request titled + `chore(main): release `, carrying the generated `CHANGELOG.md` and + the new version in `.release-please-manifest.json`. It accrues every commit + since the last release. +3. Read the changelog. This is the review, and the only one — after the merge + nothing else is asked. +4. Merge it. release-please creates the tag and the GitHub Release; the + `publish` job builds the sdist and wheel from **the tag**, checks their + version against it, and uploads to PyPI; Zenodo mints the DOI from the + release webhook. + +To skip a release, do not merge the PR. It stays open and keeps accruing. + +## Things that are the way they are for a reason + +**The publish job is in `release.yml`, not a separate `publish.yml`.** The +obvious design — a workflow keyed on `release: published` — never fires. +release-please creates the release with the default `GITHUB_TOKEN`, and +"events triggered by the `GITHUB_TOKEN` will not create a new workflow run". +The workaround is a personal access token in secrets, which is the one thing +Trusted Publishing exists to avoid, so the publish job gates on +release-please's own `release_created` output instead. Zenodo is unaffected: +it listens to the release *webhook*, which is not subject to that restriction. + +**`include-component-in-tag` is `false`.** Left on, the tag is +`specmod-v0.2.0`, which `pyproject.toml`'s `--match v[0-9]*` does not describe +and its `tag_regex` does not parse — so the build would fall back to +`0.1.1.postN.devN` and upload under that name, permanently. +`tests/test_release_config.py` asserts the two formats agree, and +`tools/check_built_version.py` checks the actual artefact between the build +and the upload. + +**`bump-minor-pre-major` is `true`.** The history already contains two +breaking commits (`feat!: make the package installable and importable`, +`fix!: default multitaper adaptive weighting off`). Without this setting +either one reads as a `1.0.0` — a number that says the API has stopped +moving, minted with a DOI that cannot be withdrawn. Below 1.0 a breaking +change bumps the minor instead, which is what §3.1 of the refactor plan +assumes throughout. + +**The changelog sections are listed explicitly.** The default preset hides +`refactor`, `docs`, `build`, `test`, `ci`, `style` and `chore`. Measured over +this repository's 146 conventional commits, that default would print 73 of +them and drop the other 73 — a release that is mostly a refactor would ship an +almost empty changelog. `refactor`, `docs` and `build` are shown; `test`, +`ci`, `style` and `chore` stay hidden. + +**The first release will be enormous, and that is correct.** There are no tags +in this repository, so the manifest declares `0.1.1` — the version the Magna +paper cites, which is the pre-refactor code. Everything since then genuinely +is the delta, so the first changelog covers the whole refactor. With the two +breaking commits and `bump-minor-pre-major`, the version it proposes is +`0.2.0`, which is what §7 of the plan expects at the end of Phase 2. + +## Checking a release went out + +- **PyPI**: the version appears at , and + `pip install specmod==` in a clean environment works. +- **The version is real**: `python -c "import specmod; print(specmod.__version__)"` + from that install prints the tag without its `v`, not a `.postN.devN`. +- **Zenodo**: a new DOI under the concept DOI, with `CITATION.cff`'s metadata. +- **Docs**: Read the Docs builds the new tag as its own version and moves + `stable` onto it. Two things to check the first few times: that the tag was + activated (an automation rule does it, see + [Documentation workflow](documentation.md#read-the-docs-setup)), and that the + new version's sidebar shows the release number rather than `0.0.0` — the + latter means the shallow-clone fix in `.readthedocs.yaml` did not take. + `latest` moved earlier, when the release PR merged. + +If the `publish` job fails after the release exists — a transient PyPI error, +say — re-run that job from the Actions UI. Do not re-run `release-please` +expecting a second attempt: `release_created` is only true on the run that +created the release. diff --git a/docs/roadmap.md b/docs/roadmap.md new file mode 100644 index 0000000..b4d6606 --- /dev/null +++ b/docs/roadmap.md @@ -0,0 +1,109 @@ +# Roadmap + +SpecMod is being rebuilt in stages, each one independently usable. This page +says what is done, what is being worked on, and what 1.0 will mean. + +**No dates.** The stages are ordered by dependency, not by calendar, and the +order after the current one can still change. A stage is listed as done when +it is merged and tested, not when it is designed — +[`REFACTOR_PLAN.md`](https://github.com/sgjholt/SpecMod/blob/main/docs/REFACTOR_PLAN.md) +is the working document behind this one and carries the reasoning, the +measurements and the open questions. + +## Where it is now + +Alpha, pre-1.0. Everything below marked done is merged, tested against golden +references, and usable — but names and signatures still move between `0.x` +releases without a deprecation cycle. See +[Releasing the software](releasing.md) for what a version number means here. + +## The stages + +### 1. An installable package ✅ + +A real `pyproject.toml`, `src/` layout, snake_case modules, linting, type +checking and a test suite on CI. Before this, the package could not be +installed or imported without editing paths by hand. + +### 2. Configuration without globals ✅ + +Settings live in a `config/` package with semantic sections and layered +overrides, instead of module-level constants read at import time. Every output +records the configuration that produced it, so a locally-overridden run is +reproducible from its own outputs. + +### 3. Publishing: docs, PyPI, DOI ✅ + +This site, built from the repository and deployed on merge; automated +changelog and version derivation from tags; PyPI upload through Trusted +Publishing; a Zenodo DOI per release. Deliberately built early, while the +package is small enough that debugging the pipeline is cheap. + +The repository side is complete. Publishing also needs a handful of one-time +repository settings, which are listed in +[Releasing the software](releasing.md). + +### 4. The transform layer ✅ + +One `SpectralEstimator` protocol with interchangeable implementations — +`FFTEstimator`, `WelchEstimator`, `MultitaperEstimator` and Prieto's — plus +Konno–Ohmachi smoothing and log-binning as separate, composable steps. The +`mtspec` Fortran dependency, which no longer builds on current toolchains, is +demoted to an optional legacy backend rather than being the only path. + +What each estimator does to real data is measured in +[Choosing a transform](choosing_a_transform.md); the differences are large +enough to matter to a magnitude. + +### 5. Wavelets ✅ + +A continuous wavelet transform alongside the Fourier estimators: +`CWTEstimator`, a `Scalogram` that tracks its cone of influence, quality +checks over it, and time-averaging back to a spectrum with ground-motion units +preserved. + +### 6. Decomposition and typed I/O 🚧 + +Breaking the remaining large modules into packages with narrow +responsibilities, and replacing pickle-based persistence with typed, portable +formats. + +Done so far: the spectral core, the fitting layer and the pick readers are +packages; the type-checking backlog is empty. Still to come: the I/O and +plotting layers, and making the remaining operations return new objects +instead of mutating in place. + +### 7. 1.0 — the API stops moving + +The theory page that states the normalisation and units conventions +explicitly, the tutorial rebuilt as an executed notebook so it cannot rot +silently, and a "what changed since 0.1" guide for anyone upgrading. Then the +1.0 release, which is the promise that names and signatures stop moving +without a deprecation cycle. + +## After 1.0 + +Designed but deliberately not on the path to 1.0, because each is blocked on +an input rather than on effort — mostly a real file from a real tool, which a +guess written from documentation cannot substitute for: + +- **Presets for picker output** (PhaseNet, EQTransformer, SeisBench), so those + users write no column mapping. The configurable reader already exists; a + preset is a mapping and a test against one real file. +- **Confirming two more pick formats** (NonLinLoc `.hyp`, IMS/GSE bulletins) + that cannot be round-tripped through ObsPy and so are listed as unconfirmed + rather than claimed. +- **SeisComP picks**, which ObsPy's SCML reader currently discards. A test + fails when that starts working upstream. +- **Recording pick-resolution policy in the configuration**, alongside the + provenance stamp. +- **A fuller `acquire --verify`** that re-fetches and diffs against the + manifest, rather than only detecting local tampering. + +## How this page will change + +Stages are the honest unit while nothing has been released: there is no +version to point at, so "done" can only mean merged. Once releases exist each +completed stage gets the version it shipped in, and this page becomes a list +of milestones against those versions rather than a list of stages — the same +content, anchored to something a reader can install and check. diff --git a/pyproject.toml b/pyproject.toml index 43e5173..f004a7d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -67,12 +67,16 @@ dev = [ "mypy>=1.11", "pre-commit>=3.7", ] +# `sphinx-autodoc-typehints` was here and is not any more: measured on this +# package it produced the same 367 documented objects as `sphinx.ext.autodoc` +# with `autodoc_typehints = "description"`, while calling an API Sphinx 10 +# removes. `myst-nb` stays for the Phase 6 notebook pages, which are not built +# yet — see `docs/conf.py`. docs = [ "sphinx>=7.3", "pydata-sphinx-theme>=0.15", "myst-parser>=3.0", "myst-nb>=1.1", - "sphinx-autodoc-typehints>=2.2", ] tutorial = [ # for running the tutorial notebook, which is not part of the package "ipykernel>=7.0", @@ -232,6 +236,16 @@ filterwarnings = [ "error::DeprecationWarning:specmod.*", # obspy 1.5 on py3.11 trips importlib.metadata's SelectableGroups warning. "ignore:SelectableGroups dict interface is deprecated:DeprecationWarning", + # matplotlib's mathtext calling pyparsing's camelCase API, which pyparsing + # 3.3 deprecates. Only the `floors` job sees it, and it drowns that job: + # 439 warnings against 4 elsewhere. Measured with pyparsing 3.3.2 + # throughout — matplotlib 3.9.0, 3.9.4 and 3.10.0 each produce 328 of them + # in `test_io_and_plotting.py` alone, and 3.11.1 produces none. So it is + # fixed above the declared `matplotlib>=3.9` floor, and raising the floor + # to silence a warning that says nothing about this package would drop + # support for two matplotlib minors for cosmetics. Delete this line when + # the floor passes 3.11. + "ignore::pyparsing.warnings.PyparsingDeprecationWarning", ] [tool.coverage.run] diff --git a/release-please-config.json b/release-please-config.json new file mode 100644 index 0000000..e8c8a89 --- /dev/null +++ b/release-please-config.json @@ -0,0 +1,26 @@ +{ + "$schema": "https://raw.githubusercontent.com/googleapis/release-please/main/schemas/config.json", + "packages": { + ".": { + "release-type": "simple", + "package-name": "specmod", + "include-component-in-tag": false, + "bump-minor-pre-major": true, + "bump-patch-for-minor-pre-major": false, + "changelog-path": "CHANGELOG.md", + "changelog-sections": [ + { "type": "feat", "section": "Features" }, + { "type": "fix", "section": "Bug Fixes" }, + { "type": "perf", "section": "Performance" }, + { "type": "refactor", "section": "Code Refactoring" }, + { "type": "docs", "section": "Documentation" }, + { "type": "build", "section": "Build System" }, + { "type": "revert", "section": "Reverts" }, + { "type": "test", "section": "Tests", "hidden": true }, + { "type": "ci", "section": "Continuous Integration", "hidden": true }, + { "type": "style", "section": "Styles", "hidden": true }, + { "type": "chore", "section": "Miscellaneous", "hidden": true } + ] + } + } +} diff --git a/src/specmod/api.py b/src/specmod/api.py new file mode 100644 index 0000000..20b90dc --- /dev/null +++ b/src/specmod/api.py @@ -0,0 +1,480 @@ +"""Stable public surface for downstream packages. + +Anything not exported here is internal and may change without notice. Anything +exported here follows the deprecation policy in ``CONTRIBUTING.md``: one minor +version of ``DeprecationWarning`` before a removal or a signature change, even +while SpecMod is ``0.x``. + +The point of the module is containment. SpecMod's internals are still being +refactored and its own documentation warns of breaking changes at every ``0.x`` +release; downstream packages import *this* and nothing else, so an internal +rename costs a line here instead of a release there. + +Five properties hold for everything below, and they are what make the surface +usable from a service that owns its own IO and has to be able to replay a job: + +1. **Path-free.** Every function takes in-memory data — arrays, or ObsPy + objects. None of them opens a file. Convenience wrappers that take paths + live elsewhere in the package. +2. **Deterministic.** The same inputs and the same explicit arguments produce + the same outputs. Nothing here reads the working directory or the + environment, and nothing draws random numbers. See the caveat on + :func:`fit_spectrum`. +3. **Non-mutating.** Inputs are left as they were found; results are new + objects. +4. **Quiet.** Nothing prints. Diagnostics go through :mod:`logging` and + :mod:`warnings`. +5. **Typed errors.** Failures are :class:`~specmod.exceptions.SpecModError` + subclasses — see :mod:`specmod.exceptions` for which of the three, and why + the distinction is the useful part. + +Examples +-------- +>>> import numpy as np +>>> from specmod import api +>>> rng = np.random.default_rng(0) +>>> signal = api.estimate_spectrum(rng.normal(size=2048), 0.01, +... estimator="multitaper") +>>> noise = api.estimate_spectrum(rng.normal(size=1024), 0.01, +... estimator="multitaper") +>>> pair = api.compare_spectra(signal, noise) +>>> pair.snr.shape == pair.binned_signal.freq.shape +True +""" + +from __future__ import annotations + +import importlib.util +from collections.abc import Iterator, Mapping +from contextlib import contextmanager +from dataclasses import dataclass +from typing import Any + +import numpy as np +from numpy.typing import ArrayLike, NDArray + +from . import __version__ +from .config import Config, ResolvedConfig, config_hash, load_config +from .config.serialize import to_toml as _to_toml +from .core.collection import SpectrumPair +from .core.spectrum import Spectrum +from .core.units import AmplitudeKind, Motion +from .exceptions import ( + InternalError, + InvalidInputError, + MissingBackendError, + SpecModError, +) +from .fitting import FitSpectrum, fittable_signal, initial_guess +from .transforms import ESTIMATORS, get_estimator +from .transforms.base import make_window, window_correction + +__all__ = [ + "AmplitudeKind", + "Config", + "InternalError", + "InvalidInputError", + "MissingBackendError", + "Motion", + "ResolvedConfig", + "SpecModError", + "Spectrum", + "SpectrumFit", + "SpectrumPair", + "__version__", + "available_estimators", + "compare_spectra", + "config_hash", + "config_to_toml", + "estimate_spectrum", + "fit_spectrum", + "load_config", + "make_window", + "window_correction", +] + +#: Which distribution each estimator needs beyond a default install. +#: +#: Measured rather than inferred, by constructing every registered estimator +#: and running it in an environment with none of the extras present: `cwt` and +#: `quadratic` are implemented against numpy and scipy and work without their +#: nominal extras, and only `prieto` actually requires one. Guessing from the +#: extras table in ``pyproject.toml`` would have marked three unavailable. +_ESTIMATOR_REQUIRES: Mapping[str, str | None] = { + "fft": None, + "welch": None, + "multitaper": None, + "quadratic": None, + "cwt": None, + "prieto": "multitaper", +} + + +@contextmanager +def _typed_errors() -> Iterator[None]: + """Translate the builtins internals raise into the documented hierarchy. + + At the boundary rather than inside, because the internals are still moving + and this module is the thing that is supposed to stay still. + """ + try: + yield + except SpecModError: + raise + except ImportError as error: + raise MissingBackendError(str(error)) from error + except (ValueError, TypeError, KeyError) as error: + raise InvalidInputError(str(error)) from error + + +def available_estimators() -> tuple[str, ...]: + """The estimators that can actually run in this environment, sorted. + + SpecMod installs without its optional backends, so the registry is not the + same question as what will work. Ask this before offering a choice to a + user, rather than discovering the answer as a failed job. + + Returns + ------- + tuple of str + Names accepted by ``estimator=`` on :func:`estimate_spectrum`. + + Examples + -------- + >>> "fft" in available_estimators() + True + """ + available = [] + for name, requires in sorted(_ESTIMATOR_REQUIRES.items()): + if name not in ESTIMATORS: # pragma: no cover - registry drift + continue + if requires is None: + available.append(name) + continue + try: + found = importlib.util.find_spec(requires) is not None + except (ImportError, ValueError): + # A blocked or broken module. `--without-optional-extras` installs + # a finder that raises ModuleNotFoundError from `find_spec`, which + # is exactly the "not available" answer. + found = False + if found: + available.append(name) + return tuple(available) + + +def estimate_spectrum( + data: ArrayLike, + dt: float, + *, + estimator: str, + motion: Motion | str = Motion.VELOCITY, + meta: Mapping[str, Any] | None = None, + **options: Any, +) -> Spectrum: + """Estimate the amplitude spectrum of one in-memory record. + + Parameters + ---------- + data + The record, as a 1-D array of samples. Not modified. + dt + Sample interval in seconds. + estimator + Which backend, from :func:`available_estimators`. Required rather than + defaulted: the configured default is a property of a study, and a + service that resolves it silently cannot replay a job it did not + record. + motion + The ground-motion domain the record is in. Carried on the result, and + what makes converting between domains a typed operation later. + meta + Extra metadata to attach to the spectrum. Copied, not held. + **options + Passed to the estimator's constructor — ``n_tapers``, + ``time_bandwidth`` and so on. Backend-specific. + + Returns + ------- + Spectrum + Frequency axis, amplitude, and the metadata needed to interpret both. + + Raises + ------ + InvalidInputError + The record is empty, not 1-D, contains non-finite values, or the + estimator name is not known. + MissingBackendError + The estimator needs an optional extra that is not installed. + """ + with _typed_errors(): + backend = get_estimator(estimator, **options) + return backend.estimate( + np.asarray(data), dt, motion=motion, meta=dict(meta) if meta else None + ) + + +def compare_spectra( + signal: Spectrum, + noise: Spectrum, + **settings: Any, +) -> SpectrumPair: + """Judge a signal spectrum against its noise window. + + Returns the pair, **including the per-bin signal-to-noise curve** rather + than only the band derived from it: ``pair.snr`` is an array aligned with + ``pair.binned_signal.freq``, and ``pair.band`` is one summary of it. A + consumer that needs a different threshold, or that admits data bin by bin + rather than over a contiguous interval, needs the curve — and a curve + cannot be recovered from a stored interval. + + Parameters + ---------- + signal, noise + Spectra from :func:`estimate_spectrum`. Neither is modified. + **settings + ``threshold``, ``f_min``, ``f_max``, ``n_bins``, ``noise_model``, + ``bandwidth`` and the rest of + :meth:`specmod.core.collection.SpectrumPair.compare`. All have explicit + defaults; none is read from configuration. + + Returns + ------- + SpectrumPair + With ``binned_signal``, ``binned_noise``, ``snr``, ``band`` and + ``resolution_floor``. + + Raises + ------ + InvalidInputError + The two spectra do not describe the same record geometry — a frequency + axis above its own Nyquist, most often from pairing windows that came + from different sampling rates. + """ + with _typed_errors(): + return SpectrumPair.compare(signal, noise, **settings) + + +@dataclass(frozen=True, slots=True) +class SpectrumFit: + """The result of fitting a source model to one spectrum. + + Frozen, and holding plain numbers rather than the fitter's own objects, so + it can be serialised and compared without depending on lmfit's API. + + Attributes + ---------- + params + Fitted values, keyed by name — ``llpsp`` (the long-period spectral + level Ω₀, as its base-10 logarithm), ``fc``, ``ts`` (t\\*). + stderr + One standard error per parameter, where the fitter could estimate one. + **Empty under some minimisers** — see the note on + :func:`fit_spectrum`. Absent means not measured, and is left absent + rather than filled with a zero that would read as "certain". + covariance + The covariance matrix, with ``names`` giving its row and column order, + or ``None`` when the minimiser produced none. The ``fc``-``t*`` + correlation lives here, and reporting either parameter without it + overstates both. + chisqr, redchi + Misfit, and misfit per degree of freedom. + n_points + How many spectral samples the fit actually used. + success + Whether the minimiser reported convergence. + """ + + params: Mapping[str, float] + stderr: Mapping[str, float] + covariance: NDArray[np.float64] | None + names: tuple[str, ...] + chisqr: float + redchi: float + n_points: int + success: bool + + def correlation(self, a: str, b: str) -> float | None: + """Correlation between two fitted parameters, or ``None``. + + ``None`` when there is no covariance matrix, or when either parameter + has no variance to correlate — not zero, which would read as + "independent" rather than "not measured". + """ + if self.covariance is None or a not in self.names or b not in self.names: + return None + i, j = self.names.index(a), self.names.index(b) + denominator = np.sqrt(self.covariance[i, i] * self.covariance[j, j]) + if not np.isfinite(denominator) or denominator == 0: + return None + return float(self.covariance[i, j] / denominator) + + +def fit_spectrum( + pair: SpectrumPair, + *, + id: str = "", + model: Any = None, + guess: Mapping[str, float] | None = None, + fit_bins: bool = False, + method: str | None = None, + weight_method: str | None = None, + **fit_options: Any, +) -> SpectrumFit: + """Fit a source model to one spectrum, with its uncertainty. + + This is the **per-spectrum** fit. SpecMod does not do a joint per-event + inversion: :class:`specmod.fitting.FitSpectra` loops over stations and fits + each independently, sharing no parameters between them. A joint solver + belongs to whoever needs one, on top of this. + + Parameters + ---------- + pair + From :func:`compare_spectra`. Its selected band is what gets fitted. + Not modified. + id + Label carried into the result's metadata. + model + A model object, or ``None`` for the configured default. + guess + Starting values for the fitted parameters. ``None`` derives them from + the spectrum with :func:`specmod.fitting.initial_guess`, which is what + :class:`specmod.fitting.FitSpectra` does. Do not skip it: without a + starting corner frequency the minimiser walks ``fc`` to zero and the + model evaluates to NaN, so an unguessed fit does not merely fit worse, + it raises. + fit_bins + Fit the log-binned spectrum rather than the full-resolution one. + method + Minimiser name, passed to lmfit. ``None`` takes ``[fitting] method`` + from configuration, which is what :class:`specmod.fitting.FitSpectra` + does — so a single-spectrum fit here matches the same station's fit in + an event run. Naming it explicitly is what makes the call reproducible + somewhere else, and the default matters: on the 28 PNR windows lmfit's + own default returns a negative corner frequency on one station where + the configured ``powell`` does not. + weight_method + ``"log"`` weights residuals by ``1/f``; ``"none"`` does not. ``None`` + takes ``[fitting] weight_method`` from configuration. + **fit_options + Anything else lmfit's ``fit`` accepts. + + Returns + ------- + SpectrumFit + Point estimates *and* their errors and covariance. Frozen. + + Raises + ------ + InvalidInputError + The pair has no usable band, or the spectrum is missing an attribute + the model needs. + + Notes + ----- + **Uncertainty depends on the minimiser, and the configured default does + not provide it.** Only the least-squares family produces a covariance + matrix. Measured on one synthetic station, all four agree on the corner + frequency and only two report an error for it: + + ========== ======== ============ =================== + ``method`` ``fc`` ``fc`` error ``fc``-``t*`` corr. + ========== ======== ============ =================== + powell 7.925 -- -- + nelder 7.925 -- -- + leastsq 7.925 0.129 0.837 + ========== ======== ============ =================== + + ``[fitting] method`` ships as ``powell``, so a default fit returns point + estimates with an empty ``stderr`` and no covariance. Pass + ``method="leastsq"`` when the uncertainty is the point. That correlation is + not incidental: 0.84 between ``fc`` and ``t*`` is why neither should be + quoted alone. + + **One determinism caveat, and it is the only one on this surface.** The + initial guess and the default minimiser are read from configuration by + internals, through :func:`specmod.config.load_config`, which resolves + against the current working directory and the environment. Two runs in the + same process with the same working directory agree exactly; two runs in + different directories may not, if a ``specmod.toml`` differs between them. + + Pass ``model`` and the minimiser options explicitly to close that gap, and + record :func:`config_hash` alongside any result you intend to replay. + """ + with _typed_errors(): + if method is None or weight_method is None: + fitting = load_config().config.fitting + method = fitting.method if method is None else method + weight_method = ( + fitting.weight_method if weight_method is None else weight_method + ) + if weight_method not in ("log", "none"): + raise InvalidInputError( + f"Unknown weight_method {weight_method!r}; expected 'log' or 'none'." + ) + + signal = fittable_signal(pair, id) + if signal is None: + raise InvalidInputError( + f"{id or 'this pair'} has no usable band, so there is nothing " + "to fit. `SpectrumPair.passes` reports that before you get here." + ) + if guess is None: + guess = initial_guess({id: pair}, model).get(id, {}) + + fitter = FitSpectrum(signal, model, fit_bins, **guess) + if weight_method == "log": + fit_options["weights"] = 1 / fitter.mod_freq + fitter.fit_mod(method=method, **fit_options) + return _as_fit(fitter) + + +def _as_fit(fitter: FitSpectrum) -> SpectrumFit: + """Extract the numbers from lmfit's result object. + + Nothing is computed here that the fitter did not already produce; this is + a projection, so that the public type does not have lmfit in it. + """ + result = fitter.result + if result is None: # pragma: no cover - fit_mod always assigns + raise InternalError("fit produced no result object") + + names = tuple(result.params) + params = {name: float(result.params[name].value) for name in names} + + stderr = {} + for name in names: + error = result.params[name].stderr + if error is not None: + stderr[name] = float(error) + + covariance = None + if result.covar is not None: + covariance = np.asarray(result.covar, dtype=np.float64) + # lmfit's covariance covers only the parameters it varied, in their + # order — not every parameter in `params`, which may include fixed + # ones. Using `names` for its axes would mislabel the matrix. + names = tuple(name for name in names if result.params[name].vary) + + return SpectrumFit( + params=params, + stderr=stderr, + covariance=covariance, + names=names, + chisqr=float(result.chisqr), + redchi=float(result.redchi), + n_points=int(result.ndata), + success=bool(result.success), + ) + + +def config_to_toml(config: Config, *, header: str | None = None) -> str: + """Serialise a configuration to TOML, as ``specmod config freeze`` does. + + Returns the text rather than writing it, so the caller decides where it + goes — which for anything but a local filesystem is the only workable + arrangement. + """ + with _typed_errors(): + return _to_toml(config, header=header) diff --git a/src/specmod/core/collection.py b/src/specmod/core/collection.py index ff3afed..d72d43e 100644 --- a/src/specmod/core/collection.py +++ b/src/specmod/core/collection.py @@ -492,6 +492,16 @@ def band(self) -> tuple[float, float] | None: def passes(self) -> bool: return self.pair.passes + @property + def motion(self) -> Motion: + """The ground-motion domain the arrays are in. + + Read through like the rest, because a fitter has to know it: the + initial guess for ``fc`` is the frequency of the spectral peak, and + that is the corner only in velocity. + """ + return self.pair.signal.motion + @dataclass(frozen=True) class SpectrumSet: diff --git a/src/specmod/exceptions.py b/src/specmod/exceptions.py new file mode 100644 index 0000000..b9a2ebe --- /dev/null +++ b/src/specmod/exceptions.py @@ -0,0 +1,58 @@ +"""The exception hierarchy :mod:`specmod.api` raises. + +Three kinds, because a caller does three different things with them: + +- :class:`InvalidInputError` — the caller's data or arguments are wrong. A + program shows a form error; a person fixes the input. +- :class:`MissingBackendError` — the code is fine and the environment is not: + an optional extra is not installed. Fixed by installing something, and + avoidable up front with :func:`specmod.api.available_estimators`. +- :class:`InternalError` — an invariant inside SpecMod is broken. Nothing the + caller can do; it is a bug report. + +Each also inherits the builtin exception the corresponding internal code +raises today, so ``except ValueError`` keeps working and internals can migrate +one at a time without a flag day. + +**Internals still raise the builtins.** :mod:`specmod.api` translates at its +own boundary, so the guarantee is specific: functions reached *through +`specmod.api`* raise this hierarchy. Reaching around it gets the builtins. +""" + +from __future__ import annotations + +__all__ = [ + "InternalError", + "InvalidInputError", + "MissingBackendError", + "SpecModError", +] + + +class SpecModError(Exception): + """Base class for every error SpecMod raises deliberately.""" + + +class InvalidInputError(SpecModError, ValueError): + """The data or arguments given to SpecMod are not usable. + + A record containing NaN, a frequency axis that does not belong to its + record, an unknown estimator name, a band with no samples in it. + """ + + +class MissingBackendError(SpecModError, ImportError): + """An optional backend is not installed. + + Raised at call time rather than import time, so a default install stays + importable. :func:`specmod.api.available_estimators` answers the same + question without provoking the error. + """ + + +class InternalError(SpecModError, RuntimeError): + """An invariant inside SpecMod does not hold. + + Not caused by the caller and not fixable by them. If one of these reaches + you, it is a bug in SpecMod. + """ diff --git a/src/specmod/fitting.py b/src/specmod/fitting.py deleted file mode 100644 index 4111511..0000000 --- a/src/specmod/fitting.py +++ /dev/null @@ -1,712 +0,0 @@ -"""Fitting a source model to one spectrum, and to a whole event. - -The lmfit surface used here is declared in ``stubs/lmfit``; see -``stubs/README.md``. lmfit ships no annotations, so without those a -``ModelResult`` is `Any` and nothing checks that ``result.redchi`` exists or -that ``Parameter.stderr`` can be `None` — which it is under every minimiser -that estimates no covariance matrix, including the shipped default. -""" - -from __future__ import annotations - -import inspect -from copy import deepcopy -from typing import TYPE_CHECKING, Any - -import lmfit as lm -import matplotlib.pyplot as plt -import numpy as np -import pandas as pd -from matplotlib.ticker import NullFormatter, StrMethodFormatter - -from . import config as cfg -from . import sources -from .tables import read_table, write_table - -if TYPE_CHECKING: # pragma: no cover - from collections.abc import Mapping - from pathlib import Path - - from matplotlib.axes import Axes - from numpy.typing import NDArray - -#: What a spectrum-like object handed to :class:`FitSpectrum` looks like. -#: Structural rather than nominal on purpose — see :meth:`FitSpectrum.__check_input`. -Spectrumish = Any -#: A container mapping trace ids to paired spectra; see -#: :meth:`FitSpectra.__check_spectra`. -SpectraLike = Any - -# global variables -# One home for this: it used to be defined in *both* the SPECTRAL and FITTING -# dicts, and the two copies could disagree. -PLOT_COLUMNS = cfg.load_config().config.viz.plot_columns - -#: What :class:`FitSpectrum` reads off whatever it is given. Kept as data so -#: the requirement is stated once and can be asserted against. -REQUIRED_SPECTRUM_ATTRIBUTES = ("id", "meta", "freq", "amp", "bfreq", "bamp") - - -def fittable_signal(pair: Any, id: str = "") -> Spectrumish | None: - """The signal to fit from a paired spectrum, or ``None`` to skip it. - - Skipping is a decision the container should not have to spell out at every - call site: a pair is unfittable when the signal-to-noise gate rejected it. - - What comes back for a :class:`~specmod.core.SpectrumPair` is its - :class:`~specmod.core.collection.FittableView`, not its ``signal``. The - pair keeps the unbinned and binned spectra as separate objects, which is - right for the comparison and wrong for a fitter that wants ``freq``, - ``amp``, ``bfreq`` and ``bamp`` side by side; the view is what puts them - there. ``id`` names the station on it, since a frozen pair does not carry - one of its own. - - The ``getattr`` fallback below is what a spectrum-like object that is not - a pair takes — a bare view, or anything else presenting the same - attributes. It is not a legacy shim; it is what lets the fitter be given - something constructed by hand. - """ - view = getattr(pair, "for_fitting", None) - if view is not None: - return pair.for_fitting(id) if pair.passes else None - - signal = getattr(pair, "signal", pair) - passes = getattr(pair, "passes", None) - if passes is None: - passes = getattr(signal, "pass_snr", True) - return signal if passes else None - - -def initial_guess( - spectra: SpectraLike, model: Any = None -) -> dict[str, dict[str, float]]: - """Starting parameters for every fittable spectrum in ``spectra``. - - Replaces ``model_guess.create_simple_guess`` and its ``_fdep`` twin, which - were two near-identical functions differing only in whether they added an - ``a`` for frequency-dependent Q — so adding a third model meant writing a - third guess function, and picking the wrong one gave lmfit a parameter the - model did not take. - - **Which parameters are needed is asked of the model, not assumed.** The - fitted callable declares them in its signature, so a model gets exactly the - guesses it takes and nothing else. Values that cannot be read off the - spectrum come from ``[fitting]`` in the configuration. - - The two that *are* read off the spectrum: - - ``llpsp`` - ``log10`` of the largest amplitude inside the selected band — the - long-period plateau, which is what ``Omega`` is. - ``fc`` - the frequency at which that maximum falls. - - Both assume a **velocity** spectrum, where the peak sits near the corner. - On a displacement spectrum the peak is at the low-frequency end and ``fc`` - would start at the bottom of the band; that is the pre-existing assumption, - made explicit here rather than left in a function name. - - Stations with no band are omitted rather than given ``None`` guesses. The - old version emitted ``{"llpsp": None, "fc": None, "ts": None}`` on - ``IndexError``, which lmfit cannot use — the failure simply moved to the - fit call. - """ - if model is None: - model = sources.from_config() - callable_ = ( - model.as_callable() if isinstance(model, sources.SpectralModel) else model - ) - wanted = set(inspect.signature(callable_).parameters) - {"f"} - - fitting = cfg.load_config().config.fitting - #: Parameters no spectrum can suggest a value for. - defaults = { - "ts": fitting.initial_t_star, - "a": fitting.initial_alpha, - } - - guesses: dict[str, dict[str, float]] = {} - for id in spectra: - signal = fittable_signal(spectra[id], id) - if signal is None: - continue - band = selected_band(signal) - if band is None: - continue - inside = (signal.freq >= band[0]) & (signal.freq <= band[1]) - if not inside.any(): - continue - - amp, freq = signal.amp[inside], signal.freq[inside] - peak = int(amp.argmax()) - available = { - "llpsp": float(np.log10(amp[peak])), - "fc": float(freq[peak]), - **defaults, - } - missing = wanted - set(available) - if missing: - raise ValueError( - f"no initial guess is defined for {sorted(missing)}, which " - f"{getattr(model, 'describe', lambda: callable_.__name__)()} " - f"takes. Add it to specmod.config.FittingConfig and to " - f"`initial_guess`, or pass explicit guesses." - ) - guesses[id] = {k: v for k, v in available.items() if k in wanted} - - return guesses - - -def selected_band(spectrum: Any) -> tuple[float, float] | None: - """The band to fit over, or ``None`` to fit everything available. - - ``None`` rather than an empty array, because "no band survived" and "a band - from 0 to 0" are different claims and the legacy spelling — an empty - ``ubfreqs`` — could be read as either. - """ - band = getattr(spectrum, "band", None) - if band is None: - return None - return (float(band[0]), float(band[1])) - - -class FitSpectrum: - """Fit a source model to one spectrum with lmfit. - - Takes anything carrying :data:`REQUIRED_SPECTRUM_ATTRIBUTES` — in practice - a :class:`~specmod.core.collection.FittableView` from - :func:`fittable_signal`. - """ - - #: Declarations, not defaults. These were class attributes carrying `None` - #: and `{}`, which meant two things at once: every read had to cope with a - #: `None` that `__init__` had in fact replaced, and `meta = {}` was one - #: dictionary shared by every instance ever constructed. `__init__` assigns - #: all of them, so the type is what it is after construction — and the - #: shared-mutable-default hazard is gone rather than merely unreached. - sig: Spectrumish - mod: lm.Model - params: lm.Parameters - #: `None` until :meth:`fit_mod` runs. This one really is optional, and - #: callers test it — see :func:`specmod.plotting.plot_pair`. - result: lm.ModelResult | None - mod_freq: NDArray[np.float64] - mod_amp: NDArray[np.float64] - pass_fitting: bool - fit_bins: bool - meta: dict[str, Any] - #: The :class:`specmod.sources.SpectralModel` behind the fit, when there is - #: one. ``None`` if a bare callable was supplied. - spectral_model: sources.SpectralModel | None - - def __init__( - self, - signal: Spectrumish, - model: Any = None, - fit_bins: bool = False, - **params: float, - ) -> None: - self.result = None - self.pass_fitting = True - self.meta = {} - self.spectral_model = None - self.mod_freq = np.array([]) - self.mod_amp = np.array([]) - self.fit_bins = fit_bins - self.set_signal(signal) - self.set_model(model, **params) - - def fit_mod(self, **kwargs: Any) -> None: - """Fit, judge the result, then record it — in that order. - - The judgement used to be made *after* the recording, so the - ``pass_fitting`` column of every flat file held the value from before - the fit ran — ``True``, the class default, on a fresh `FitSpectrum`. - The attribute and the table disagreed, and the table is what gets - written out and regressed on. - """ - self.result = self.mod.fit(self.mod_amp, self.params, f=self.mod_freq, **kwargs) - self.__determine_pass_or_fail() - self.__set_results_to_meta() - - def set_signal(self, signal: Spectrumish) -> None: - if self.__check_input(signal): - self.sig = signal - self.__set_meta(signal.meta) - # if setting a new signal - assess and adjust the freq bounds - self.__set_mod_amp_freq() - - def set_model(self, model: Any = None, **params: float) -> None: - """Set the model to fit. - - Accepts a :class:`specmod.sources.SpectralModel`, a bare callable, or - ``None`` — in which case the model is whatever ``[model]`` in the - configuration asks for. That default is the point: before it existed, - ``config.model.source`` was read by nothing and the caller had to pass - the right function by hand, so a study file saying - ``source = "boatwright"`` silently got Brune. - - A bare callable still works, because fitting an ad-hoc shape is a - legitimate thing to want. It simply carries no provenance: - :attr:`spectral_model` is ``None`` and nothing can report what was fitted. - """ - if model is None: - model = sources.from_config() - - if isinstance(model, sources.SpectralModel): - self.spectral_model = model - model = model.as_callable() - else: - self.spectral_model = None - - self.mod = lm.Model(model) - # whenever a model is set the inital params must be set also - self.__init_params(**params) - - @property - def fitted(self) -> lm.ModelResult: - """The fit result, or a message saying it has not been fitted. - - Every private reader below went straight through ``self.result``, - which is ``None`` until :meth:`fit_mod` runs — so calling - :meth:`quick_vis` on an unfitted spectrum raised ``AttributeError: - 'NoneType' object has no attribute 'best_fit'``, from a line that - names neither the station nor the missing step. - """ - if self.result is None: - raise RuntimeError( - f"{getattr(self.sig, 'id', 'this spectrum')} has not been " - "fitted yet; call fit_mod() first" - ) - return self.result - - def describe_model(self) -> str | None: - """What is being fitted, or ``None`` for a bare callable.""" - return None if self.spectral_model is None else self.spectral_model.describe() - - def set_const(self, pname: str, value: float) -> None: - self.params[pname].value = value - self.params[pname].vary = False - - def set_bounds( - self, pname: str, min: float | None = None, max: float | None = None - ) -> None: - if min is not None: - self.params[pname].min = min - if max is not None: - self.params[pname].max = max - - def __set_meta(self, meta: Mapping[str, Any]) -> None: - self.meta = deepcopy(dict(meta)) - - def __init_params(self, **params: Any) -> None: - """Seed the parameters, and floor ``t*`` where the configuration says. - - ``fitting.t_star_min`` existed and was read by nothing. The tutorial - did ``fits.set_bounds("ts", min=0.0001)`` by hand and the config value - is 1e-4 — the same number — so the setting was a written-down record of - something every caller had to remember. Applied here, forgetting it is - no longer possible. - - The same applies to ``fc``, and the legacy code knew it — the line - ``# 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. - """ - # `**params: Any` rather than `float`, because lmfit's `make_params` - # takes a leading `verbose` argument: a model with a parameter of that - # name would have its seed swallowed as a flag. Not a hazard for any - # source model here, and not one this package can fix. - self.params = self.mod.make_params(**params) - fitting = cfg.load_config().config.fitting - for name, floor in ( - ("ts", fitting.t_star_min), - ("fc", fitting.corner_frequency_min), - ): - if name in self.params and floor is not None: - self.set_bounds(name, min=floor) - - def reset(self) -> None: - for par in self.params.values(): - par.vary = True - par.min = -np.inf - par.max = np.inf - - def __check_input(self, signal: Spectrumish) -> bool: - """Accept anything carrying what the fit reads, not one named class. - - This used to be ``isinstance(signal, spectral.Signal)``, which is the - coupling that kept the container holding the legacy pair — nothing - could be handed to the fitter unless it *was* that class. What the fit - actually needs is the six attributes below, so that is what is checked. - - Named explicitly rather than left to fail at first use: a missing - ``bamp`` should say so here, not surface as an AttributeError from - inside a band selection three calls later. - """ - missing = [ - name for name in REQUIRED_SPECTRUM_ATTRIBUTES if not hasattr(signal, name) - ] - if missing: - raise ValueError( - f"{type(signal).__name__} cannot be fitted: missing " - f"{', '.join(missing)}. A fittable spectrum needs " - f"{', '.join(REQUIRED_SPECTRUM_ATTRIBUTES)}." - ) - return True - - def __set_mod_amp_freq(self) -> None: - """ - Only fit between signal limits if they are specified. - """ - - if self.fit_bins: - freq = self.sig.bfreq - amp = self.sig.bamp - else: - freq = self.sig.freq - amp = self.sig.amp - - band = selected_band(self.sig) - if band is not None: - inds = np.where((freq >= band[0]) & (freq <= band[1])) - self.mod_freq = freq[inds] - self.mod_amp = amp[inds] - else: - self.mod_freq = freq - self.mod_amp = amp - - self.mod_amp = np.log10(self.mod_amp) - - def __param_string(self) -> str: - """``name: value+/-2sigma`` per parameter, or ``name: value`` alone. - - The old version computed ``2 * k.stderr`` unconditionally inside a - bare ``except Exception``. ``stderr`` is ``None`` whenever the - minimiser estimated no covariance matrix — which Powell, the shipped - default, never does — so this raised ``TypeError`` on every fit made - with the default configuration, swallowed it, and titled the plot - ``NaN``. A missing uncertainty is a property of the method, not a - failed fit, so the value is still worth printing. - """ - parts = [] - for k in self.fitted.params.values(): - if k.stderr is None: - parts.append(f"{k.name}: {k.value:.3f}") - else: - parts.append(f"{k.name}: {k.value:.3f}+/-{2 * k.stderr:.3f}") - return ", ".join(parts) - - def quick_vis(self, ax: Axes | None = None) -> Axes: - if ax is None: - _fig, ax = plt.subplots(1, 1) - - ax.loglog(self.mod_freq, 10**self.mod_amp, color="grey", label=self.sig.id) - ax.loglog(self.mod_freq, 10**self.fitted.best_fit, "k--", label="model") - ax.xaxis.set_major_formatter(StrMethodFormatter("{x:.2f}")) - ax.xaxis.set_minor_formatter(NullFormatter()) - ax.set_title(self.__param_string()) - ax.set_xlabel("freq [Hz]") - ax.set_ylabel("spectral amp") - ax.legend() - return ax - - def __get_pars(self) -> dict[str, Any]: - p: dict[str, Any] = {} - for k in self.fitted.params.values(): - p.update({k.name: k.value}) - p.update({k.name + "-stderr": k.stderr}) - return p - - def __get_fit_stats(self) -> dict[str, float]: - res = self.fitted - s: dict[str, float] = {} - s.update({"aic": res.aic}) - s.update({"bic": res.bic}) - s.update({"chisqr": res.chisqr}) - s.update({"redchi": res.redchi}) - return s - - def __get_test_results(self) -> dict[str, bool]: - t: dict[str, bool] = {} - t.update({"pass_fitting": self.pass_fitting}) - return t - - def __set_results_to_meta(self) -> None: - self.meta.update(self.__get_pars()) - self.meta.update(self.__get_fit_stats()) - self.meta.update(self.__get_test_results()) - - def __determine_pass_or_fail(self) -> None: - """A fit fails when a parameter is pinned against one of its bounds. - - Which is the useful question: a corner frequency resting on its floor - is the minimiser saying "lower, if you would let me", and the value it - reports is the bound rather than a measurement. - - Reset first. ``pass_fitting`` starts as a class attribute and was only - ever set *False*, so a `FitSpectrum` that failed once could never pass - again however many times it was refitted. - - **Where there is no uncertainty, the value itself is compared.** The - old version treated a missing ``stderr`` as a failure, which would mark - every fit failed under the shipped configuration: Powell does not - estimate a covariance matrix, so lmfit has no uncertainties to report. - That is a property of the minimiser, not a fault in the fit. Asking - whether the value sits on the bound is the same question with the - error bar removed. - """ - self.pass_fitting = True - for _par, vals in self.fitted.params.items(): - if not vals.vary: - continue - spread = vals.stderr if vals.stderr is not None else 0.0 - if (vals.value - spread <= vals.min) or (vals.value + spread >= vals.max): - self.pass_fitting = False - - -class FitSpectra: - """Fit every passing station in an event.""" - - #: Declarations, as on :class:`FitSpectrum`. `models = {}` at class level - #: was one dictionary shared by every `FitSpectra` ever built; `__init__` - #: rebinds it, so nothing reached the shared copy, but nothing prevented it - #: either. `guess = {}` was never assigned anywhere at all — a class - #: attribute recording a constructor argument that is not kept. - spectra: SpectraLike - models: dict[str, FitSpectrum] - table: pd.DataFrame - - def __init__( - self, - spectra: SpectraLike, - model: Any = None, - guess: Mapping[str, Mapping[str, float]] | None = None, - fit_bins: bool | None = None, - ) -> None: - """``guess=None`` derives one, rather than fitting nothing. - - It used to skip `init_fitting` entirely, so `FitSpectra(spectra)` built - an object with no models and `fit_spectra()` silently did nothing and - produced an empty table. There is a sensible guess available — see - :func:`initial_guess` — so that is now the default and an explicit - ``guess={}`` is how you say "none". - """ - self.models = {} - self.table = pd.DataFrame([]) - self.set_spectra(spectra) - if fit_bins is None: - fit_bins = cfg.load_config().config.fitting.fit_bins - if guess is None: - guess = initial_guess(spectra, model) - self.init_fitting(model, guess, fit_bins) - - def __len__(self) -> int: - return len(self.models) - - def set_spectra(self, spectra: SpectraLike) -> None: - if self.__check_spectra(spectra): - self.spectra = spectra - - def get_spectra(self) -> SpectraLike: - return self.spectra - - def get_fit(self, id: str) -> FitSpectrum | None: - if id.upper() in self.models: - return self.models[id.upper()] - print(f"WARNING: {id.upper()} not in group of available fits.") - return None - - def fit_spectra(self, weight_method: str | None = None, **kwargs: Any) -> None: - """Fit every station, with the configured minimiser unless told otherwise. - - ``method`` and ``weight_method`` both come from ``[fitting]`` when not - given. Neither used to: `fit_spectra()` fell through to lmfit's default - minimiser, so a study file saying ``method = "powell"`` was ignored and - the caller had to remember ``fit_spectra(method="powell")`` — which the - tutorial does and nothing enforced. - - It matters. On the 28 PNR windows lmfit's default returns a **negative - corner frequency** on one station where Powell does not; a corner - frequency below zero is not a degraded measurement but a meaningless - one, and nothing downstream rejects it. - """ - fitting = cfg.load_config().config.fitting - if weight_method is None: - weight_method = fitting.weight_method - kwargs.setdefault("method", fitting.method) - wm = self.__check_wm(weight_method) - for name, mod in self.models.items(): - try: - if wm == "log": - mod.fit_mod(weights=1 / mod.mod_freq, **kwargs) - else: - mod.fit_mod(**kwargs) - except ValueError as msg: - print(msg) - print("-" * 40) - print(f"Skipping {name}") - - self.__set_fit_models_to_spectrum() - self.__generate_group_fit_table() - - def init_fitting( - self, - model: Any, - guess: Mapping[str, Mapping[str, float]], - fit_bins: bool, - ) -> None: - """Build a fit per passing station. - - ``model=None`` resolves through the configuration once per station, - which is cheap and keeps every fit in a run agreeing on what it is - fitting. - """ - # Iterate the container rather than reaching into `.group`. `Spectra` - # and `core.SpectrumSet` both present this interface, which is what - # lets the container be swapped underneath without touching the fitter. - # A station is fitted when it passed the gate *and* has a guess. - # Indexing `guess[id]` unconditionally made a partial guess dict a - # `KeyError` naming a station, rather than a way to fit a subset — - # and made `guess={}` a crash instead of "fit nothing". - tmp: dict[str, FitSpectrum] = {} - for id in self.spectra: - signal = fittable_signal(self.spectra[id], id) - if signal is None or id not in guess: - continue - tmp[id] = FitSpectrum(signal, model, **guess[id], fit_bins=fit_bins) - self.models = tmp - - def set_const(self, pname: str, value: float, id: str | None = None) -> None: - if id is None: - for mod in self.models.values(): - mod.set_const(pname, value) - elif id in self.models: - self.models[id].set_const(pname, value) - - def set_bounds( - self, pname: str, min: float | None = None, max: float | None = None - ) -> None: - for mod in self.models.values(): - mod.set_bounds(pname, min, max) - - def reset(self, name: str = "all") -> None: - """Unbind every parameter, on one station or all of them. - - The lookup tested ``name.upper()`` for membership and then indexed with - ``name``, so any id not already upper-case passed the check and raised - ``KeyError`` on the next line. Station ids are upper-case in practice, - which is why it never fired. - """ - if name.upper() == "ALL": - for mod in self.models.values(): - mod.reset() - return - - id = name.upper() - if id in self.models: - self.models[id].reset() - else: - print(f"WARNING: {id} not in available channels.") - - def quick_vis(self, save: str | None = None) -> None: - rows = self.__num_rows() - fig, axes = plt.subplots(rows, PLOT_COLUMNS, figsize=(17, int(rows * 5))) - # `strict=False`: the grid is rounded up to whole rows, so there are - # more axes than models by construction. - for ax, mod in zip(axes.flatten(), self.models.values(), strict=False): - if mod.result is None or not mod.pass_fitting: - ax.set_title(f"Fitting Failed for {mod.sig.id}") - else: - mod.quick_vis(ax) - - if save is not None: - if type(save) is str: - fig.savefig(save) - else: - raise ValueError("Must provide valid path as str.") - - @staticmethod - def write_flatfile(path: str | Path, fits: FitSpectra) -> Path: - """Write the group fit table, in the format ``path``'s suffix names. - - ``.parquet`` is typed, compressed and queryable without loading; - ``.csv`` is what journal supplements want. See :mod:`specmod.tables`. - - The previous implementation was ``os.makedirs(os.path.join( - *path.split("/")[:-1]))``, which raised ``TypeError: join() missing 1 - required positional argument`` for any path without a directory - component — ``write_flatfile("out.csv", fits)`` could not work. It also - split on ``/`` literally, so it did nothing useful on Windows. - """ - return write_table(path, fits.table) - - @staticmethod - def read_flatfile(path: str | Path) -> pd.DataFrame: - """Read a fit table back. Format follows the suffix.""" - return read_table(path) - - def __check_wm(self, wm: str) -> str: - if wm not in ["log", "none"]: - print(f"WARNING: did not recognise weight method {wm}.") - print("Setting to none...") - wm = "none" - return wm - - def __generate_group_fit_table(self) -> None: - ds = [m.meta for m in self.models.values()] - df1 = pd.DataFrame([]) - for i, d in enumerate(ds): - df1 = pd.concat( - [df1, pd.DataFrame(d, index=[i])], ignore_index=True, sort=False - ) - self.table = df1 - - def __set_fit_models_to_spectrum(self) -> None: - """Hand each fit back to the spectrum it came from, where that is possible. - - The legacy `Signal` carries its own fit so that plotting and - serialisation can reach it from the spectrum. `core.SpectrumPair` is - frozen and cannot, by design — a result writing itself back into its - own input is how a container stops being trustworthy. - - Nothing is lost by skipping it: `self.models` is the source of truth - either way, and the write-back was only ever a convenience. So this - writes where the container accepts it and moves on where it does not, - rather than requiring every container to be mutable. - """ - for id, mod in self.models.items(): - spectrum = self.spectra[id] - signal = getattr(spectrum, "signal", spectrum) - setter = getattr(signal, "set_model", None) - if setter is not None: - setter(mod) - - def __check_spectra(self, spectra: SpectraLike) -> bool: - """Accept anything that maps trace ids to paired spectra. - - Was ``isinstance(spectra, spectral.Spectra)``, which is why the fitter - could not be handed a :class:`~specmod.core.SpectrumSet` even though it - only ever iterates and indexes. Requiring one concrete class was the - last thing tying the fitter to the legacy module. - """ - required = ("__iter__", "__getitem__", "__len__") - missing = [name for name in required if not hasattr(spectra, name)] - if missing: - raise ValueError( - f"{type(spectra).__name__} cannot be fitted: it must map trace " - f"ids to paired spectra, and is missing {', '.join(missing)}. " - f"Use specmod.pipeline.spectrum_set_from_streams." - ) - return True - - def __num_rows(self) -> int: - count = len(self) - cols = PLOT_COLUMNS - if count % cols > 0: - return int((cols * (int(count / cols) + 1)) / cols) - return int(count / cols) diff --git a/src/specmod/fitting/__init__.py b/src/specmod/fitting/__init__.py new file mode 100644 index 0000000..977ff9f --- /dev/null +++ b/src/specmod/fitting/__init__.py @@ -0,0 +1,58 @@ +"""Fitting a source model to one spectrum, and to a whole event. + +:func:`fittable_signal` decides what to fit, :func:`initial_guess` where to +start, :class:`FitSpectrum` fits one station and :class:`FitSpectra` an event. +""" + +from __future__ import annotations + +import warnings + +from .base import ( + REQUIRED_SPECTRUM_ATTRIBUTES, + SpectraLike, + Spectrumish, + plot_columns, +) +from .event import FitSpectra +from .guess import fittable_signal, initial_guess, selected_band +from .spectrum import FitSpectrum + +__all__ = [ + "PLOT_COLUMNS", + "REQUIRED_SPECTRUM_ATTRIBUTES", + "FitSpectra", + "FitSpectrum", + "SpectraLike", + "Spectrumish", + "fittable_signal", + "initial_guess", + "plot_columns", + "selected_band", +] + + +def __getattr__(name: str) -> object: + """Keep ``PLOT_COLUMNS`` importable, resolved at each access. + + It was a module-level constant evaluated at import time, which froze the + configuration of whatever directory the process started in. Reading it now + resolves configuration per access, so the value is at least correct — but + a name that looks like a constant and performs a lookup is a poor bargain + either way, hence the warning and :func:`plot_columns`. + + Deliberately *not* ``from .base import ...``: a from-import would run this + once at import time and rebind the result, which is the frozen behaviour + this replaced, reintroduced one level up. + """ + if name == "PLOT_COLUMNS": + warnings.warn( + "specmod.fitting.PLOT_COLUMNS is deprecated and will be removed " + "in 0.4.0; call specmod.fitting.plot_columns() instead, which " + "resolves the current configuration rather than the one that " + "happened to be in effect at import.", + DeprecationWarning, + stacklevel=2, + ) + return plot_columns() + raise AttributeError(f"module {__name__!r} has no attribute {name!r}") diff --git a/src/specmod/fitting/base.py b/src/specmod/fitting/base.py new file mode 100644 index 0000000..f8b9281 --- /dev/null +++ b/src/specmod/fitting/base.py @@ -0,0 +1,50 @@ +"""What the fitting layer expects of whatever it is handed. + +Structural rather than nominal: :class:`~specmod.fitting.FitSpectrum` reads +attributes off its input and does not care what class provides them, which is +what lets it take a :class:`~specmod.core.collection.FittableView`, a bare +spectrum, or something assembled by hand in a notebook. +""" + +from __future__ import annotations + +from typing import Any + +from .. import config as cfg + +__all__ = [ + "REQUIRED_SPECTRUM_ATTRIBUTES", + "SpectraLike", + "Spectrumish", + "plot_columns", +] + +#: What a spectrum-like object handed to :class:`FitSpectrum` looks like. +#: Structural rather than nominal on purpose — see :meth:`FitSpectrum.__check_input`. +Spectrumish = Any +#: A container mapping trace ids to paired spectra; see +#: :meth:`FitSpectra.__check_spectra`. +SpectraLike = Any + + +def plot_columns() -> int: + """How many columns a multi-panel figure uses, from ``[viz]``. + + A function rather than a constant, and that is the whole point. This was + ``PLOT_COLUMNS = cfg.load_config().config.viz.plot_columns`` evaluated at + import time, so importing :mod:`specmod.fitting` resolved configuration + against whatever directory the process happened to start in and froze the + answer for the life of the interpreter. Measured: importing from a project + whose ``specmod.toml`` says 5, then moving to one that resolves to 3, left + the constant at 5 — a worker serving two projects would use the first + one's layout for both. + + One home for the setting either way: it used to be defined in *both* the + SPECTRAL and FITTING dicts, and the two copies could disagree. + """ + return int(cfg.load_config().config.viz.plot_columns) + + +#: What :class:`FitSpectrum` reads off whatever it is given. Kept as data so +#: the requirement is stated once and can be asserted against. +REQUIRED_SPECTRUM_ATTRIBUTES = ("id", "meta", "freq", "amp", "bfreq", "bamp") diff --git a/src/specmod/fitting/event.py b/src/specmod/fitting/event.py new file mode 100644 index 0000000..d186116 --- /dev/null +++ b/src/specmod/fitting/event.py @@ -0,0 +1,284 @@ +"""Fitting every passing station in an event, and the table that comes out.""" + +from __future__ import annotations + +import logging +import warnings +from typing import TYPE_CHECKING, Any + +import matplotlib.pyplot as plt +import pandas as pd + +from .. import config as cfg +from ..tables import read_table, write_table +from .base import SpectraLike, plot_columns +from .guess import fittable_signal, initial_guess +from .spectrum import FitSpectrum + +if TYPE_CHECKING: # pragma: no cover + from collections.abc import Mapping + from pathlib import Path + +__all__ = ["FitSpectra"] + + +#: Per-station progress goes here rather than to stdout. A library must not +#: configure logging for its host, so there is no `basicConfig` anywhere in +#: this package: a caller that wants to see these calls `logging.basicConfig` +#: itself, and one that does not is not written to by surprise. +logger = logging.getLogger(__name__) + + +class FitSpectra: + """Fit every passing station in an event.""" + + #: Declarations, as on :class:`FitSpectrum`. `models = {}` at class level + #: was one dictionary shared by every `FitSpectra` ever built; `__init__` + #: rebinds it, so nothing reached the shared copy, but nothing prevented it + #: either. `guess = {}` was never assigned anywhere at all — a class + #: attribute recording a constructor argument that is not kept. + spectra: SpectraLike + models: dict[str, FitSpectrum] + table: pd.DataFrame + + def __init__( + self, + spectra: SpectraLike, + model: Any = None, + guess: Mapping[str, Mapping[str, float]] | None = None, + fit_bins: bool | None = None, + ) -> None: + """``guess=None`` derives one, rather than fitting nothing. + + It used to skip `init_fitting` entirely, so `FitSpectra(spectra)` built + an object with no models and `fit_spectra()` silently did nothing and + produced an empty table. There is a sensible guess available — see + :func:`initial_guess` — so that is now the default and an explicit + ``guess={}`` is how you say "none". + """ + self.models = {} + self.table = pd.DataFrame([]) + self.set_spectra(spectra) + if fit_bins is None: + fit_bins = cfg.load_config().config.fitting.fit_bins + if guess is None: + guess = initial_guess(spectra, model) + self.init_fitting(model, guess, fit_bins) + + def __len__(self) -> int: + return len(self.models) + + def set_spectra(self, spectra: SpectraLike) -> None: + if self.__check_spectra(spectra): + self.spectra = spectra + + def get_spectra(self) -> SpectraLike: + return self.spectra + + def get_fit(self, id: str) -> FitSpectrum | None: + if id.upper() in self.models: + return self.models[id.upper()] + warnings.warn( + f"{id.upper()} is not among the fitted stations " + f"({', '.join(sorted(self.models)) or 'none'}); returning None.", + stacklevel=2, + ) + return None + + def fit_spectra(self, weight_method: str | None = None, **kwargs: Any) -> None: + """Fit every station, with the configured minimiser unless told otherwise. + + ``method`` and ``weight_method`` both come from ``[fitting]`` when not + given. Neither used to: `fit_spectra()` fell through to lmfit's default + minimiser, so a study file saying ``method = "powell"`` was ignored and + the caller had to remember ``fit_spectra(method="powell")`` — which the + tutorial does and nothing enforced. + + It matters. On the 28 PNR windows lmfit's default returns a **negative + corner frequency** on one station where Powell does not; a corner + frequency below zero is not a degraded measurement but a meaningless + one, and nothing downstream rejects it. + """ + fitting = cfg.load_config().config.fitting + if weight_method is None: + weight_method = fitting.weight_method + kwargs.setdefault("method", fitting.method) + wm = self.__check_wm(weight_method) + for name, mod in self.models.items(): + try: + if wm == "log": + mod.fit_mod(weights=1 / mod.mod_freq, **kwargs) + else: + mod.fit_mod(**kwargs) + except ValueError as error: + # `logging`, not `warnings`, and the difference matters here: + # warnings are deduplicated per code location by default, so a + # run that skipped twenty stations would report one. Each skip + # is a station missing from the results and has to be visible. + logger.warning("skipping %s: %s", name, error) + + self.__set_fit_models_to_spectrum() + self.__generate_group_fit_table() + + def init_fitting( + self, + model: Any, + guess: Mapping[str, Mapping[str, float]], + fit_bins: bool, + ) -> None: + """Build a fit per passing station. + + ``model=None`` resolves through the configuration once per station, + which is cheap and keeps every fit in a run agreeing on what it is + fitting. + """ + # Iterate the container rather than reaching into `.group`. `Spectra` + # and `core.SpectrumSet` both present this interface, which is what + # lets the container be swapped underneath without touching the fitter. + # A station is fitted when it passed the gate *and* has a guess. + # Indexing `guess[id]` unconditionally made a partial guess dict a + # `KeyError` naming a station, rather than a way to fit a subset — + # and made `guess={}` a crash instead of "fit nothing". + tmp: dict[str, FitSpectrum] = {} + for id in self.spectra: + signal = fittable_signal(self.spectra[id], id) + if signal is None or id not in guess: + continue + tmp[id] = FitSpectrum(signal, model, **guess[id], fit_bins=fit_bins) + self.models = tmp + + def set_const(self, pname: str, value: float, id: str | None = None) -> None: + if id is None: + for mod in self.models.values(): + mod.set_const(pname, value) + elif id in self.models: + self.models[id].set_const(pname, value) + + def set_bounds( + self, pname: str, min: float | None = None, max: float | None = None + ) -> None: + for mod in self.models.values(): + mod.set_bounds(pname, min, max) + + def reset(self, name: str = "all") -> None: + """Unbind every parameter, on one station or all of them. + + The lookup tested ``name.upper()`` for membership and then indexed with + ``name``, so any id not already upper-case passed the check and raised + ``KeyError`` on the next line. Station ids are upper-case in practice, + which is why it never fired. + """ + if name.upper() == "ALL": + for mod in self.models.values(): + mod.reset() + return + + id = name.upper() + if id in self.models: + self.models[id].reset() + else: + warnings.warn( + f"{id} is not among the fitted stations; nothing was reset.", + stacklevel=2, + ) + + def quick_vis(self, save: str | None = None) -> None: + rows = self.__num_rows() + fig, axes = plt.subplots(rows, plot_columns(), figsize=(17, int(rows * 5))) + # `strict=False`: the grid is rounded up to whole rows, so there are + # more axes than models by construction. + for ax, mod in zip(axes.flatten(), self.models.values(), strict=False): + if mod.result is None or not mod.pass_fitting: + ax.set_title(f"Fitting Failed for {mod.sig.id}") + else: + mod.quick_vis(ax) + + if save is not None: + if type(save) is str: + fig.savefig(save) + else: + raise ValueError("Must provide valid path as str.") + + @staticmethod + def write_flatfile(path: str | Path, fits: FitSpectra) -> Path: + """Write the group fit table, in the format ``path``'s suffix names. + + ``.parquet`` is typed, compressed and queryable without loading; + ``.csv`` is what journal supplements want. See :mod:`specmod.tables`. + + The previous implementation was ``os.makedirs(os.path.join( + *path.split("/")[:-1]))``, which raised ``TypeError: join() missing 1 + required positional argument`` for any path without a directory + component — ``write_flatfile("out.csv", fits)`` could not work. It also + split on ``/`` literally, so it did nothing useful on Windows. + """ + return write_table(path, fits.table) + + @staticmethod + def read_flatfile(path: str | Path) -> pd.DataFrame: + """Read a fit table back. Format follows the suffix.""" + return read_table(path) + + def __check_wm(self, wm: str) -> str: + if wm not in ["log", "none"]: + warnings.warn( + f"Unknown weight method {wm!r}; expected 'log' or 'none'. " + "Falling back to 'none'.", + stacklevel=3, + ) + wm = "none" + return wm + + def __generate_group_fit_table(self) -> None: + ds = [m.meta for m in self.models.values()] + df1 = pd.DataFrame([]) + for i, d in enumerate(ds): + df1 = pd.concat( + [df1, pd.DataFrame(d, index=[i])], ignore_index=True, sort=False + ) + self.table = df1 + + def __set_fit_models_to_spectrum(self) -> None: + """Hand each fit back to the spectrum it came from, where that is possible. + + The legacy `Signal` carries its own fit so that plotting and + serialisation can reach it from the spectrum. `core.SpectrumPair` is + frozen and cannot, by design — a result writing itself back into its + own input is how a container stops being trustworthy. + + Nothing is lost by skipping it: `self.models` is the source of truth + either way, and the write-back was only ever a convenience. So this + writes where the container accepts it and moves on where it does not, + rather than requiring every container to be mutable. + """ + for id, mod in self.models.items(): + spectrum = self.spectra[id] + signal = getattr(spectrum, "signal", spectrum) + setter = getattr(signal, "set_model", None) + if setter is not None: + setter(mod) + + def __check_spectra(self, spectra: SpectraLike) -> bool: + """Accept anything that maps trace ids to paired spectra. + + Was ``isinstance(spectra, spectral.Spectra)``, which is why the fitter + could not be handed a :class:`~specmod.core.SpectrumSet` even though it + only ever iterates and indexes. Requiring one concrete class was the + last thing tying the fitter to the legacy module. + """ + required = ("__iter__", "__getitem__", "__len__") + missing = [name for name in required if not hasattr(spectra, name)] + if missing: + raise ValueError( + f"{type(spectra).__name__} cannot be fitted: it must map trace " + f"ids to paired spectra, and is missing {', '.join(missing)}. " + f"Use specmod.pipeline.spectrum_set_from_streams." + ) + return True + + def __num_rows(self) -> int: + count = len(self) + cols = plot_columns() + if count % cols > 0: + return int((cols * (int(count / cols) + 1)) / cols) + return int(count / cols) diff --git a/src/specmod/fitting/guess.py b/src/specmod/fitting/guess.py new file mode 100644 index 0000000..a160fda --- /dev/null +++ b/src/specmod/fitting/guess.py @@ -0,0 +1,170 @@ +"""Choosing what to fit, and where the fit starts from.""" + +from __future__ import annotations + +import inspect +import warnings +from typing import Any + +import numpy as np + +from .. import config as cfg +from .. import sources +from ..core.units import Motion +from .base import SpectraLike, Spectrumish + +__all__ = ["fittable_signal", "initial_guess", "selected_band"] + + +def fittable_signal(pair: Any, id: str = "") -> Spectrumish | None: + """The signal to fit from a paired spectrum, or ``None`` to skip it. + + Skipping is a decision the container should not have to spell out at every + call site: a pair is unfittable when the signal-to-noise gate rejected it. + + What comes back for a :class:`~specmod.core.SpectrumPair` is its + :class:`~specmod.core.collection.FittableView`, not its ``signal``. The + pair keeps the unbinned and binned spectra as separate objects, which is + right for the comparison and wrong for a fitter that wants ``freq``, + ``amp``, ``bfreq`` and ``bamp`` side by side; the view is what puts them + there. ``id`` names the station on it, since a frozen pair does not carry + one of its own. + + The ``getattr`` fallback below is what a spectrum-like object that is not + a pair takes — a bare view, or anything else presenting the same + attributes. It is not a legacy shim; it is what lets the fitter be given + something constructed by hand. + """ + view = getattr(pair, "for_fitting", None) + if view is not None: + return pair.for_fitting(id) if pair.passes else None + + signal = getattr(pair, "signal", pair) + passes = getattr(pair, "passes", None) + if passes is None: + passes = getattr(signal, "pass_snr", True) + return signal if passes else None + + +def _warn_if_peak_is_meaningless(signal: Spectrumish, id: str) -> None: + """Warn when the peak-as-``fc`` guess is being read off the wrong domain. + + A no-op for velocity, and for a spectrum that does not say what motion it + carries — something assembled by hand is the caller's business. + """ + motion = getattr(signal, "motion", None) + if motion is None or Motion(motion) is Motion.VELOCITY: + return + warnings.warn( + f"{id or 'this spectrum'} is in {Motion(motion).value}, and the " + f"initial guess for fc is the frequency of the spectral peak — which " + f"is the corner only in velocity. A {Motion(motion).value} spectrum " + f"falls monotonically across the band, so the guess will be a band " + f"edge and the fit will settle near it. Fit the velocity spectrum; " + f"`llpsp` is the displacement plateau either way.", + stacklevel=3, + ) + + +def initial_guess( + spectra: SpectraLike, model: Any = None +) -> dict[str, dict[str, float]]: + """Starting parameters for every fittable spectrum in ``spectra``. + + Replaces ``model_guess.create_simple_guess`` and its ``_fdep`` twin, which + were two near-identical functions differing only in whether they added an + ``a`` for frequency-dependent Q — so adding a third model meant writing a + third guess function, and picking the wrong one gave lmfit a parameter the + model did not take. + + **Which parameters are needed is asked of the model, not assumed.** The + fitted callable declares them in its signature, so a model gets exactly the + guesses it takes and nothing else. Values that cannot be read off the + spectrum come from ``[fitting]`` in the configuration. + + The two that *are* read off the spectrum: + + ``llpsp`` + ``log10`` of the largest amplitude inside the selected band — the + long-period plateau, which is what ``Omega`` is. + ``fc`` + the frequency at which that maximum falls. + + Both assume a **velocity** spectrum, which is where a fit belongs anyway: + the model carries a motion factor, so ``llpsp`` is the displacement plateau + whichever domain is fitted, but converting first is not a neutral change of + view — integrating implicitly low-passes and differentiating amplifies + high-frequency noise, so the record to fit is the one the sensor recorded. + + In velocity the peak is not merely near the corner, it *is* the corner, for + any omega-squared source: the stationary point of + ``f * [1 + (f/fc)**(gamma*n)]**(-1/gamma)`` sits at ``f = fc`` whenever + ``n == 2``, whatever the corner sharpness. In displacement and acceleration + the spectrum is monotonic across the band, so the peak is whichever band + edge it was handed and the guess is meaningless. Handed one of those, this + warns rather than proceeding quietly. + + Stations with no band are omitted rather than given ``None`` guesses. The + old version emitted ``{"llpsp": None, "fc": None, "ts": None}`` on + ``IndexError``, which lmfit cannot use — the failure simply moved to the + fit call. + """ + if model is None: + model = sources.from_config() + callable_ = ( + model.as_callable() if isinstance(model, sources.SpectralModel) else model + ) + wanted = set(inspect.signature(callable_).parameters) - {"f"} + + fitting = cfg.load_config().config.fitting + #: Parameters no spectrum can suggest a value for. + defaults = { + "ts": fitting.initial_t_star, + "a": fitting.initial_alpha, + } + + guesses: dict[str, dict[str, float]] = {} + for id in spectra: + signal = fittable_signal(spectra[id], id) + if signal is None: + continue + band = selected_band(signal) + if band is None: + continue + inside = (signal.freq >= band[0]) & (signal.freq <= band[1]) + if not inside.any(): + continue + + _warn_if_peak_is_meaningless(signal, id) + + amp, freq = signal.amp[inside], signal.freq[inside] + peak = int(amp.argmax()) + available = { + "llpsp": float(np.log10(amp[peak])), + "fc": float(freq[peak]), + **defaults, + } + missing = wanted - set(available) + if missing: + raise ValueError( + f"no initial guess is defined for {sorted(missing)}, which " + f"{getattr(model, 'describe', lambda: callable_.__name__)()} " + f"takes. Add it to specmod.config.FittingConfig and to " + f"`initial_guess`, or pass explicit guesses." + ) + guesses[id] = {k: v for k, v in available.items() if k in wanted} + + return guesses + + +def selected_band(spectrum: Any) -> tuple[float, float] | None: + """The band to fit over, or ``None`` to fit everything available. + + ``None`` rather than an empty array, because "no band survived" and "a band + from 0 to 0" are different claims and the legacy spelling — an empty + ``ubfreqs`` — could be read as either. + """ + band = getattr(spectrum, "band", None) + if band is None: + return None + return (float(band[0]), float(band[1])) diff --git a/src/specmod/fitting/spectrum.py b/src/specmod/fitting/spectrum.py new file mode 100644 index 0000000..c89a085 --- /dev/null +++ b/src/specmod/fitting/spectrum.py @@ -0,0 +1,330 @@ +"""Fitting a source model to a single spectrum. + +The lmfit surface used here is declared in ``stubs/lmfit``; see +``stubs/README.md``. lmfit ships no annotations, so without those a +``ModelResult`` is `Any` and nothing checks that ``result.redchi`` exists or +that ``Parameter.stderr`` can be `None` — which it is under every minimiser +that estimates no covariance matrix, including the shipped default. +""" + +from __future__ import annotations + +from copy import deepcopy +from typing import TYPE_CHECKING, Any + +import lmfit as lm +import matplotlib.pyplot as plt +import numpy as np +from matplotlib.ticker import NullFormatter, StrMethodFormatter + +from .. import config as cfg +from .. import sources +from .base import REQUIRED_SPECTRUM_ATTRIBUTES, Spectrumish +from .guess import selected_band + +if TYPE_CHECKING: # pragma: no cover + from collections.abc import Mapping + + from matplotlib.axes import Axes + from numpy.typing import NDArray + +__all__ = ["FitSpectrum"] + + +class FitSpectrum: + """Fit a source model to one spectrum with lmfit. + + Takes anything carrying :data:`REQUIRED_SPECTRUM_ATTRIBUTES` — in practice + a :class:`~specmod.core.collection.FittableView` from + :func:`fittable_signal`. + """ + + #: Declarations, not defaults. These were class attributes carrying `None` + #: and `{}`, which meant two things at once: every read had to cope with a + #: `None` that `__init__` had in fact replaced, and `meta = {}` was one + #: dictionary shared by every instance ever constructed. `__init__` assigns + #: all of them, so the type is what it is after construction — and the + #: shared-mutable-default hazard is gone rather than merely unreached. + sig: Spectrumish + mod: lm.Model + params: lm.Parameters + #: `None` until :meth:`fit_mod` runs. This one really is optional, and + #: callers test it — see :func:`specmod.plotting.plot_pair`. + result: lm.ModelResult | None + mod_freq: NDArray[np.float64] + mod_amp: NDArray[np.float64] + pass_fitting: bool + fit_bins: bool + meta: dict[str, Any] + #: The :class:`specmod.sources.SpectralModel` behind the fit, when there is + #: one. ``None`` if a bare callable was supplied. + spectral_model: sources.SpectralModel | None + + def __init__( + self, + signal: Spectrumish, + model: Any = None, + fit_bins: bool = False, + **params: float, + ) -> None: + self.result = None + self.pass_fitting = True + self.meta = {} + self.spectral_model = None + self.mod_freq = np.array([]) + self.mod_amp = np.array([]) + self.fit_bins = fit_bins + self.set_signal(signal) + self.set_model(model, **params) + + def fit_mod(self, **kwargs: Any) -> None: + """Fit, judge the result, then record it — in that order. + + The judgement used to be made *after* the recording, so the + ``pass_fitting`` column of every flat file held the value from before + the fit ran — ``True``, the class default, on a fresh `FitSpectrum`. + The attribute and the table disagreed, and the table is what gets + written out and regressed on. + """ + self.result = self.mod.fit(self.mod_amp, self.params, f=self.mod_freq, **kwargs) + self.__determine_pass_or_fail() + self.__set_results_to_meta() + + def set_signal(self, signal: Spectrumish) -> None: + if self.__check_input(signal): + self.sig = signal + self.__set_meta(signal.meta) + # if setting a new signal - assess and adjust the freq bounds + self.__set_mod_amp_freq() + + def set_model(self, model: Any = None, **params: float) -> None: + """Set the model to fit. + + Accepts a :class:`specmod.sources.SpectralModel`, a bare callable, or + ``None`` — in which case the model is whatever ``[model]`` in the + configuration asks for. That default is the point: before it existed, + ``config.model.source`` was read by nothing and the caller had to pass + the right function by hand, so a study file saying + ``source = "boatwright"`` silently got Brune. + + A bare callable still works, because fitting an ad-hoc shape is a + legitimate thing to want. It simply carries no provenance: + :attr:`spectral_model` is ``None`` and nothing can report what was fitted. + """ + if model is None: + model = sources.from_config() + + if isinstance(model, sources.SpectralModel): + self.spectral_model = model + model = model.as_callable() + else: + self.spectral_model = None + + self.mod = lm.Model(model) + # whenever a model is set the inital params must be set also + self.__init_params(**params) + + @property + def fitted(self) -> lm.ModelResult: + """The fit result, or a message saying it has not been fitted. + + Every private reader below went straight through ``self.result``, + which is ``None`` until :meth:`fit_mod` runs — so calling + :meth:`quick_vis` on an unfitted spectrum raised ``AttributeError: + 'NoneType' object has no attribute 'best_fit'``, from a line that + names neither the station nor the missing step. + """ + if self.result is None: + raise RuntimeError( + f"{getattr(self.sig, 'id', 'this spectrum')} has not been " + "fitted yet; call fit_mod() first" + ) + return self.result + + def describe_model(self) -> str | None: + """What is being fitted, or ``None`` for a bare callable.""" + return None if self.spectral_model is None else self.spectral_model.describe() + + def set_const(self, pname: str, value: float) -> None: + self.params[pname].value = value + self.params[pname].vary = False + + def set_bounds( + self, pname: str, min: float | None = None, max: float | None = None + ) -> None: + if min is not None: + self.params[pname].min = min + if max is not None: + self.params[pname].max = max + + def __set_meta(self, meta: Mapping[str, Any]) -> None: + self.meta = deepcopy(dict(meta)) + + def __init_params(self, **params: Any) -> None: + """Seed the parameters, and floor ``t*`` where the configuration says. + + ``fitting.t_star_min`` existed and was read by nothing. The tutorial + did ``fits.set_bounds("ts", min=0.0001)`` by hand and the config value + is 1e-4 — the same number — so the setting was a written-down record of + something every caller had to remember. Applied here, forgetting it is + no longer possible. + + The same applies to ``fc``, and the legacy code knew it — the line + ``# 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. + """ + # `**params: Any` rather than `float`, because lmfit's `make_params` + # takes a leading `verbose` argument: a model with a parameter of that + # name would have its seed swallowed as a flag. Not a hazard for any + # source model here, and not one this package can fix. + self.params = self.mod.make_params(**params) + fitting = cfg.load_config().config.fitting + for name, floor in ( + ("ts", fitting.t_star_min), + ("fc", fitting.corner_frequency_min), + ): + if name in self.params and floor is not None: + self.set_bounds(name, min=floor) + + def reset(self) -> None: + for par in self.params.values(): + par.vary = True + par.min = -np.inf + par.max = np.inf + + def __check_input(self, signal: Spectrumish) -> bool: + """Accept anything carrying what the fit reads, not one named class. + + This used to be ``isinstance(signal, spectral.Signal)``, which is the + coupling that kept the container holding the legacy pair — nothing + could be handed to the fitter unless it *was* that class. What the fit + actually needs is the six attributes below, so that is what is checked. + + Named explicitly rather than left to fail at first use: a missing + ``bamp`` should say so here, not surface as an AttributeError from + inside a band selection three calls later. + """ + missing = [ + name for name in REQUIRED_SPECTRUM_ATTRIBUTES if not hasattr(signal, name) + ] + if missing: + raise ValueError( + f"{type(signal).__name__} cannot be fitted: missing " + f"{', '.join(missing)}. A fittable spectrum needs " + f"{', '.join(REQUIRED_SPECTRUM_ATTRIBUTES)}." + ) + return True + + def __set_mod_amp_freq(self) -> None: + """ + Only fit between signal limits if they are specified. + """ + + if self.fit_bins: + freq = self.sig.bfreq + amp = self.sig.bamp + else: + freq = self.sig.freq + amp = self.sig.amp + + band = selected_band(self.sig) + if band is not None: + inds = np.where((freq >= band[0]) & (freq <= band[1])) + self.mod_freq = freq[inds] + self.mod_amp = amp[inds] + else: + self.mod_freq = freq + self.mod_amp = amp + + self.mod_amp = np.log10(self.mod_amp) + + def __param_string(self) -> str: + """``name: value+/-2sigma`` per parameter, or ``name: value`` alone. + + The old version computed ``2 * k.stderr`` unconditionally inside a + bare ``except Exception``. ``stderr`` is ``None`` whenever the + minimiser estimated no covariance matrix — which Powell, the shipped + default, never does — so this raised ``TypeError`` on every fit made + with the default configuration, swallowed it, and titled the plot + ``NaN``. A missing uncertainty is a property of the method, not a + failed fit, so the value is still worth printing. + """ + parts = [] + for k in self.fitted.params.values(): + if k.stderr is None: + parts.append(f"{k.name}: {k.value:.3f}") + else: + parts.append(f"{k.name}: {k.value:.3f}+/-{2 * k.stderr:.3f}") + return ", ".join(parts) + + def quick_vis(self, ax: Axes | None = None) -> Axes: + if ax is None: + _fig, ax = plt.subplots(1, 1) + + ax.loglog(self.mod_freq, 10**self.mod_amp, color="grey", label=self.sig.id) + ax.loglog(self.mod_freq, 10**self.fitted.best_fit, "k--", label="model") + ax.xaxis.set_major_formatter(StrMethodFormatter("{x:.2f}")) + ax.xaxis.set_minor_formatter(NullFormatter()) + ax.set_title(self.__param_string()) + ax.set_xlabel("freq [Hz]") + ax.set_ylabel("spectral amp") + ax.legend() + return ax + + def __get_pars(self) -> dict[str, Any]: + p: dict[str, Any] = {} + for k in self.fitted.params.values(): + p.update({k.name: k.value}) + p.update({k.name + "-stderr": k.stderr}) + return p + + def __get_fit_stats(self) -> dict[str, float]: + res = self.fitted + s: dict[str, float] = {} + s.update({"aic": res.aic}) + s.update({"bic": res.bic}) + s.update({"chisqr": res.chisqr}) + s.update({"redchi": res.redchi}) + return s + + def __get_test_results(self) -> dict[str, bool]: + t: dict[str, bool] = {} + t.update({"pass_fitting": self.pass_fitting}) + return t + + def __set_results_to_meta(self) -> None: + self.meta.update(self.__get_pars()) + self.meta.update(self.__get_fit_stats()) + self.meta.update(self.__get_test_results()) + + def __determine_pass_or_fail(self) -> None: + """A fit fails when a parameter is pinned against one of its bounds. + + Which is the useful question: a corner frequency resting on its floor + is the minimiser saying "lower, if you would let me", and the value it + reports is the bound rather than a measurement. + + Reset first. ``pass_fitting`` starts as a class attribute and was only + ever set *False*, so a `FitSpectrum` that failed once could never pass + again however many times it was refitted. + + **Where there is no uncertainty, the value itself is compared.** The + old version treated a missing ``stderr`` as a failure, which would mark + every fit failed under the shipped configuration: Powell does not + estimate a covariance matrix, so lmfit has no uncertainties to report. + That is a property of the minimiser, not a fault in the fit. Asking + whether the value sits on the bound is the same question with the + error bar removed. + """ + self.pass_fitting = True + for _par, vals in self.fitted.params.items(): + if not vals.vary: + continue + spread = vals.stderr if vals.stderr is not None else 0.0 + if (vals.value - spread <= vals.min) or (vals.value + spread >= vals.max): + self.pass_fitting = False diff --git a/src/specmod/utils.py b/src/specmod/utils.py index be4b2da..1b94848 100644 --- a/src/specmod/utils.py +++ b/src/specmod/utils.py @@ -13,6 +13,8 @@ from __future__ import annotations import contextlib +import logging +import warnings from typing import TYPE_CHECKING, Any import matplotlib @@ -26,6 +28,10 @@ from specmod.picks import PickSet, SnufflerReader, from_catalog, select_event +#: See the note in `specmod.fitting.event`: diagnostics go here, and nothing in +#: this package configures logging on a caller's behalf. +logger = logging.getLogger(__name__) + if TYPE_CHECKING: # pragma: no cover from collections.abc import Sequence from datetime import datetime @@ -247,7 +253,7 @@ def plot_traces( fig.savefig(save) fig.clear() plt.close(fig) - print("deleted td fig") + logger.debug("closed the time-domain figure after saving to %s", save) def stream_distance_sort(st: Stream, dist_met: str = "repi") -> Stream: @@ -260,7 +266,11 @@ def stream_distance_sort(st: Stream, dist_met: str = "repi") -> Stream: try: st = obspy.Stream(sorted(st, key=lambda x: x.stats[dist_met])) except KeyError: - print("WARNING: No distance info, stream not sorted by distance.") + warnings.warn( + f"No {dist_met!r} on these traces, so the stream is returned " + "unsorted. Set distances with `preprocess.set_stream_distance`.", + stacklevel=2, + ) return st.copy() diff --git a/stubs/lmfit/model.pyi b/stubs/lmfit/model.pyi index b6855e1..c508a9f 100644 --- a/stubs/lmfit/model.pyi +++ b/stubs/lmfit/model.pyi @@ -21,6 +21,13 @@ class ModelResult: #: the shipped default, never does. `stderr` on every parameter is then #: `None`, and that is a property of the method rather than a failed fit. errorbars: bool + #: The covariance matrix over the *varied* parameters, in their order, or + #: `None` under a minimiser that estimates none — the same condition + #: `errorbars` reports. Its axes are not `params`, which may also hold + #: fixed ones. + covar: NDArray[np.float64] | None + #: Number of data points the fit actually used. + ndata: int def fit_report(self, **kwargs: Any) -> str: ... def eval(self, params: Parameters | None = ..., **kwargs: Any) -> Any: ... diff --git a/tests/test_ambient_state.py b/tests/test_ambient_state.py new file mode 100644 index 0000000..98f3aa6 --- /dev/null +++ b/tests/test_ambient_state.py @@ -0,0 +1,221 @@ +"""Two properties the package has to hold everywhere, not just where fixed. + +Both were found by auditing for what a long-lived process would do with this +library — see ``docs/notes/api_audit.md`` — and both are the kind of defect +that reappears one file at a time. So each is checked by walking the package +rather than by testing the one site that had it. + +1. **No configuration is read at import time.** ``load_config()`` resolves + against the current working directory and the environment. At module level + that answer is frozen for the life of the interpreter, so a worker serving + two projects uses the first one's settings for both. +2. **Nothing prints.** A service capturing logs per job gets nothing from + ``print``, and a CLI writing to a pipe gets its output corrupted. + Diagnostics go through :mod:`warnings` or :mod:`logging`. +""" + +from __future__ import annotations + +import ast +import logging +import warnings +from pathlib import Path +from typing import Any + +import pytest + +from specmod import fitting, utils + +SOURCE = Path(fitting.__file__).resolve().parent.parent +#: Every module in the installed package. +MODULES = sorted(path for path in SOURCE.rglob("*.py") if "_vendor" not in path.parts) + + +def _module_level_calls(node: ast.AST) -> list[ast.Call]: + """Calls evaluated when the module is imported. + + Recursion stops at a function body, which runs when something calls it — + that is the whole point of moving a config read into one. It does *not* + stop at a class body, which executes at import like any other statement, + nor at a function's default arguments and decorators, which are evaluated + at definition time and are a favourite hiding place for exactly this. + + Written as an explicit walk rather than `ast.walk`, because `ast.walk` + descends into method bodies and would flag every call-time read in the + package. It did, on the first run of this test. + """ + found: list[ast.Call] = [] + for child in ast.iter_child_nodes(node): + if isinstance(child, ast.FunctionDef | ast.AsyncFunctionDef | ast.Lambda): + for default in [*child.args.defaults, *child.args.kw_defaults]: + if default is not None: + found.extend(_calls_in(default)) + for decorator in getattr(child, "decorator_list", []): + found.extend(_calls_in(decorator)) + continue + if isinstance(child, ast.Call): + found.append(child) + found.extend(_module_level_calls(child)) + return found + + +def _calls_in(node: ast.AST) -> list[ast.Call]: + """Every call inside an expression that runs at import.""" + return [n for n in ast.walk(node) if isinstance(n, ast.Call)] + + +def _callee(call: ast.Call) -> str: + """The dotted name being called, as written.""" + node: Any = call.func + parts = [] + while isinstance(node, ast.Attribute): + parts.append(node.attr) + node = node.value + if isinstance(node, ast.Name): + parts.append(node.id) + return ".".join(reversed(parts)) + + +@pytest.mark.parametrize("path", MODULES, ids=lambda p: p.name) +class TestTheModuleIsInert: + def test_it_reads_no_configuration_at_import(self, path: Path) -> None: + """`fitting/base.py` did this, and froze `[viz] plot_columns` at + whatever the importing directory said. Measured before the fix: + importing from a project whose `specmod.toml` said 5, then moving to + one resolving to 3, left it at 5.""" + tree = ast.parse(path.read_text(), filename=str(path)) + offenders = [ + _callee(call) + for call in _module_level_calls(tree) + if _callee(call).endswith("load_config") + ] + assert not offenders, ( + f"{path.name} resolves configuration at import time " + f"({', '.join(offenders)}). Read it inside the function that needs " + "it, so the value follows the caller rather than the importer." + ) + + def test_it_does_not_print(self, path: Path) -> None: + tree = ast.parse(path.read_text(), filename=str(path)) + printed = [ + node.lineno + for node in ast.walk(tree) + if isinstance(node, ast.Call) and _callee(node) == "print" + ] + assert not printed, ( + f"{path.name} calls print() at line(s) " + f"{', '.join(map(str, printed))}. Use `warnings.warn` for something " + "the caller should act on, or the module logger for progress." + ) + + +class TestPlotColumnsFollowsTheConfiguration: + def test_it_reads_the_current_directory( + self, monkeypatch: pytest.MonkeyPatch, tmp_path: Path + ) -> None: + (tmp_path / "specmod.toml").write_text("[viz]\nplot_columns = 7\n") + monkeypatch.chdir(tmp_path) + assert fitting.plot_columns() == 7 + + def test_it_changes_when_the_directory_does( + self, monkeypatch: pytest.MonkeyPatch, tmp_path: Path + ) -> None: + """The defect itself: the value used to be whatever the *first* + resolution said, forever.""" + first, second = tmp_path / "a", tmp_path / "b" + for directory, columns in ((first, 5), (second, 2)): + directory.mkdir() + (directory / "specmod.toml").write_text( + f"[viz]\nplot_columns = {columns}\n" + ) + + monkeypatch.chdir(first) + assert fitting.plot_columns() == 5 + monkeypatch.chdir(second) + assert fitting.plot_columns() == 2 + + def test_the_old_name_still_works_and_says_it_is_going(self) -> None: + with pytest.warns(DeprecationWarning, match="plot_columns"): + value = fitting.PLOT_COLUMNS + assert value == fitting.plot_columns() + + def test_an_unknown_attribute_is_still_an_attribute_error(self) -> None: + """The `__getattr__` must not swallow typos into something else.""" + with pytest.raises(AttributeError): + fitting.NOT_A_REAL_NAME # noqa: B018 + + +class TestDiagnosticsAreAudible: + """Each path that used to print. What replaced it is chosen per site: + `warnings` for something the caller should act on, `logging` for per-item + progress — because warnings are deduplicated by code location, so a run + skipping twenty stations would report one.""" + + def test_an_unknown_weight_method_warns_and_falls_back(self) -> None: + from specmod.fitting.event import FitSpectra # noqa: PLC0415 + + checker = FitSpectra.__dict__["_FitSpectra__check_wm"] + with pytest.warns(UserWarning, match="Unknown weight method"): + assert checker(None, "sideways") == "none" + + def test_a_missing_distance_warns_rather_than_printing( + self, capsys: pytest.CaptureFixture[str] + ) -> None: + import obspy # noqa: PLC0415 + + stream = obspy.Stream([obspy.Trace()]) + with pytest.warns(UserWarning, match="unsorted"): + utils.stream_distance_sort(stream) + assert capsys.readouterr().out == "" + + def test_the_module_loggers_exist_and_are_not_configured(self) -> None: + """A library must not call `basicConfig` for its host: that decides + formatting and destination for the whole process.""" + from specmod.fitting import event # noqa: PLC0415 + + for module in (event, utils): + assert isinstance(module.logger, logging.Logger) + assert module.logger.name.startswith("specmod") + assert not module.logger.handlers + + def test_every_failed_station_is_reported_not_just_the_first( + self, caplog: pytest.LogCaptureFixture, monkeypatch: pytest.MonkeyPatch + ) -> None: + """The reason this site is logging rather than a warning. + + Driven through `fit_spectra` itself with two stations that both fail: + a station that could not be fitted is missing from the results, and + one line covering both would be a report that hides the second. + """ + from specmod.fitting.event import FitSpectra # noqa: PLC0415 + + class _Fails: + def fit_mod(self, **kwargs: Any) -> None: + raise ValueError("no usable band") + + fitter = object.__new__(FitSpectra) + fitter.models = {"XX.A..HHZ": _Fails(), "XX.B..HHZ": _Fails()} + # The two bookkeeping steps after the loop need real models. + monkeypatch.setattr( + FitSpectra, "_FitSpectra__set_fit_models_to_spectrum", lambda self: None + ) + monkeypatch.setattr( + FitSpectra, "_FitSpectra__generate_group_fit_table", lambda self: None + ) + + with caplog.at_level(logging.WARNING, logger="specmod.fitting.event"): + fitter.fit_spectra(weight_method="none") + + reported = [record.getMessage() for record in caplog.records] + assert len(reported) == 2, reported + assert any("XX.A..HHZ" in line for line in reported) + assert any("XX.B..HHZ" in line for line in reported) + assert all("no usable band" in line for line in reported) + + def test_warnings_really_would_have_collapsed_that(self) -> None: + """Not an assumption about `warnings`; the behaviour it is avoiding.""" + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("default") + for _ in range(3): + warnings.warn("same site, same message", stacklevel=1) + assert len(caught) == 1 diff --git a/tests/test_api_surface.py b/tests/test_api_surface.py new file mode 100644 index 0000000..de02c30 --- /dev/null +++ b/tests/test_api_surface.py @@ -0,0 +1,318 @@ +"""The promises :mod:`specmod.api` makes to downstream packages. + +The surface exists so that SpecMod's internals can keep moving while +`specmod-studio`, `specmod-model` and `specmod-git` import one thing that does +not. That is only true if the promises are checked, so each one here is a test +rather than a sentence in a docstring: + +- the export list is frozen, so adding or removing one shows up in review; +- every export is documented and annotated; +- the same input twice gives the same output; +- nothing touches the filesystem; +- :func:`~specmod.api.available_estimators` answers for the environment it is + actually in, including under ``--without-optional-extras``. +""" + +from __future__ import annotations + +import builtins +import inspect +import io +from pathlib import Path +from typing import Any + +import numpy as np +import pytest + +from specmod import api +from specmod.exceptions import InvalidInputError, MissingBackendError, SpecModError + +from .test_end_to_end import WINDOW_S, _trace + +#: Every name the surface promises. Adding to this list is a compatibility +#: obligation under the policy in CONTRIBUTING.md; removing one needs a +#: deprecation cycle first. The literal is duplicated here on purpose — +#: comparing `api.__all__` against itself would assert nothing. +EXPECTED_EXPORTS = ( + "AmplitudeKind", + "Config", + "InternalError", + "InvalidInputError", + "MissingBackendError", + "Motion", + "ResolvedConfig", + "SpecModError", + "Spectrum", + "SpectrumFit", + "SpectrumPair", + "__version__", + "available_estimators", + "compare_spectra", + "config_hash", + "config_to_toml", + "estimate_spectrum", + "fit_spectrum", + "load_config", + "make_window", + "window_correction", +) + + +@pytest.fixture(scope="module") +def pair() -> Any: + """A real signal/noise pair, through the public surface only.""" + import specmod.preprocess as pre # noqa: PLC0415 + + import obspy # noqa: PLC0415, isort: skip + + stream = obspy.Stream([_trace("S00", seed=11)]) + signal = pre.get_signal( + stream, pre.cut_s, rafp=0.0, tafs=WINDOW_S, time_after="absolute_time" + ) + noise = pre.get_noise_p(stream, signal) + sig, noi = signal[0], noise[0] + return api.compare_spectra( + api.estimate_spectrum(sig.data, float(sig.stats.delta), estimator="multitaper"), + api.estimate_spectrum(noi.data, float(noi.stats.delta), estimator="multitaper"), + ) + + +class TestTheExportList: + def test_it_is_exactly_what_was_agreed(self) -> None: + assert tuple(sorted(api.__all__)) == EXPECTED_EXPORTS + + @pytest.mark.parametrize("name", EXPECTED_EXPORTS) + def test_every_export_exists(self, name: str) -> None: + assert hasattr(api, name), f"{name} is promised and missing" + + @pytest.mark.parametrize("name", EXPECTED_EXPORTS) + def test_every_export_is_documented(self, name: str) -> None: + obj = getattr(api, name) + if name == "__version__": + pytest.skip("a string, not a documentable object") + assert obj.__doc__, f"{name} has no docstring" + + @pytest.mark.parametrize("name", EXPECTED_EXPORTS) + def test_every_exported_function_is_annotated(self, name: str) -> None: + obj = getattr(api, name) + if not inspect.isfunction(obj): + pytest.skip("not a function") + signature = inspect.signature(obj) + assert signature.return_annotation is not inspect.Signature.empty, name + for parameter in signature.parameters.values(): + assert parameter.annotation is not inspect.Signature.empty, ( + f"{name}({parameter.name}) is not annotated" + ) + + +class TestCapabilities: + def test_it_reports_what_actually_runs(self) -> None: + """The point of the function: no name it returns may fail to run. + + This is the invariant under both installs. With the extras present + `prieto` is in the list and works; under `--without-optional-extras` + it is absent, and the ones that remain still work. + """ + data = np.random.default_rng(0).normal(size=512) + for name in api.available_estimators(): + spectrum = api.estimate_spectrum(data, 0.01, estimator=name) + assert spectrum.freq.size > 0, name + + def test_the_backends_needing_nothing_are_always_there(self) -> None: + assert {"fft", "welch", "multitaper"} <= set(api.available_estimators()) + + def test_it_is_sorted_and_hashable(self) -> None: + names = api.available_estimators() + assert isinstance(names, tuple) + assert list(names) == sorted(names) + + def test_an_unavailable_backend_says_so_in_the_type(self) -> None: + """`prieto` needs an extra. Absent, it must raise the typed error and + not a bare ImportError, so a caller can tell 'install something' from + 'your input is wrong'.""" + if "prieto" in api.available_estimators(): + pytest.skip("the multitaper extra is installed here") + with pytest.raises(MissingBackendError): + api.estimate_spectrum(np.zeros(512) + 1.0, 0.01, estimator="prieto") + + +class TestDeterminism: + def test_estimation_repeats_exactly(self) -> None: + data = np.random.default_rng(7).normal(size=1024) + first = api.estimate_spectrum(data, 0.01, estimator="multitaper") + second = api.estimate_spectrum(data, 0.01, estimator="multitaper") + assert np.array_equal(first.amp, second.amp) + assert np.array_equal(first.freq, second.freq) + + def test_comparison_repeats_exactly(self, pair: Any) -> None: + again = api.compare_spectra(pair.signal, pair.noise) + assert np.array_equal(again.snr, pair.snr) + assert again.band == pair.band + + def test_fitting_repeats_exactly(self, pair: Any) -> None: + first = api.fit_spectrum(pair, id="XX.S00..HHN") + second = api.fit_spectrum(pair, id="XX.S00..HHN") + assert first.params == second.params + assert first.chisqr == second.chisqr + + def test_the_config_hash_is_stable(self) -> None: + config = api.load_config(use_local=False, use_env=False).config + assert api.config_hash(config) == api.config_hash(config) + + +class TestItDoesNotTouchTheFilesystem: + """Studio owns its IO, so that projects can live on S3, Azure or GCS. + + A core function that opens a path itself defeats that, and the failure is + silent on a workstation — it only shows up in a deployment where the path + does not exist. So the check is mechanical: make opening a file an error + and call the surface. + """ + + @pytest.fixture + def no_open(self, monkeypatch: pytest.MonkeyPatch, tmp_path: Path) -> None: + # An empty directory, so `load_config` finds no config to legitimately + # read, and warm caches before the ban so a first-use import inside + # numpy or lmfit is not blamed on the surface. + monkeypatch.chdir(tmp_path) + api.estimate_spectrum(np.arange(64.0), 0.01, estimator="fft") + + def refuse(*args: Any, **kwargs: Any) -> Any: + raise AssertionError(f"opened a file: {args[0]!r}") + + monkeypatch.setattr(builtins, "open", refuse) + monkeypatch.setattr(io, "open", refuse) + + def test_estimation(self, no_open: None) -> None: + data = np.random.default_rng(1).normal(size=512) + for name in api.available_estimators(): + api.estimate_spectrum(data, 0.01, estimator=name) + + def test_comparison_and_fitting(self, no_open: None, pair: Any) -> None: + again = api.compare_spectra(pair.signal, pair.noise) + fit = api.fit_spectrum(again, id="XX.S00..HHN") + assert fit.n_points > 0 + + def test_config_and_capabilities(self, no_open: None) -> None: + config = api.load_config(use_local=False, use_env=False).config + api.config_hash(config) + api.config_to_toml(config) + api.available_estimators() + + def test_tapers(self, no_open: None) -> None: + assert api.make_window("tukey", 128).size == 128 + assert api.window_correction(api.make_window("tukey", 128), "energy") > 0 + + +class TestItDoesNotMutateItsInputs: + def test_estimation_leaves_the_record_alone(self) -> None: + data = np.random.default_rng(3).normal(size=512) + before = data.copy() + api.estimate_spectrum(data, 0.01, estimator="multitaper") + assert np.array_equal(data, before) + + def test_comparison_leaves_the_spectra_alone(self, pair: Any) -> None: + signal_amp = pair.signal.amp.copy() + noise_amp = pair.noise.amp.copy() + api.compare_spectra(pair.signal, pair.noise) + assert np.array_equal(pair.signal.amp, signal_amp) + assert np.array_equal(pair.noise.amp, noise_amp) + + +class TestTheFitCarriesItsUncertainty: + """A point estimate without a bound is not a measurement, and the f_c-t* + correlation is the one downstream consumers are told to display. + + Which is available depends on the minimiser, and the configured default is + `powell`, which produces no covariance. Both halves are pinned here: asking + for it gets it, and not asking gets an honest absence rather than a + fabricated zero. + """ + + def test_the_fit_finds_the_parameters(self, pair: Any) -> None: + fit = api.fit_spectrum(pair, id="XX.S00..HHN") + assert fit.success + assert set(fit.params) >= {"llpsp", "fc", "ts"} + assert fit.n_points > 0 + + def test_least_squares_reports_errors_and_a_covariance(self, pair: Any) -> None: + fit = api.fit_spectrum(pair, id="XX.S00..HHN", method="leastsq") + assert set(fit.stderr) >= {"llpsp", "fc", "ts"} + assert fit.covariance is not None + assert fit.covariance.shape == (len(fit.names), len(fit.names)) + + def test_the_correlation_is_reachable(self, pair: Any) -> None: + fit = api.fit_spectrum(pair, id="XX.S00..HHN", method="leastsq") + correlation = fit.correlation("fc", "ts") + assert correlation is not None + assert -1.0 <= correlation <= 1.0 + + def test_a_minimiser_without_a_covariance_says_so(self, pair: Any) -> None: + """`powell` gives no covariance. The result must show that, not + invent one — a zero error is a claim, and the wrong one.""" + fit = api.fit_spectrum(pair, id="XX.S00..HHN", method="powell") + assert fit.stderr == {} + assert fit.covariance is None + assert fit.correlation("fc", "ts") is None + + def test_an_unmeasured_correlation_is_none_not_zero(self, pair: Any) -> None: + """Zero would read as 'independent' rather than 'not measured'.""" + fit = api.fit_spectrum(pair, id="XX.S00..HHN", method="leastsq") + assert fit.correlation("fc", "not_a_parameter") is None + + def test_the_result_is_frozen(self, pair: Any) -> None: + fit = api.fit_spectrum(pair, id="XX.S00..HHN") + with pytest.raises((AttributeError, TypeError)): + fit.chisqr = 0.0 # type: ignore[misc] + + +class TestTheSnrIsPerBin: + """§3.5 of the Studio design: store the curve, derive the intervals. A + scalar band cannot be un-collapsed later.""" + + def test_the_curve_is_aligned_with_the_frequency_axis(self, pair: Any) -> None: + assert pair.snr.shape == pair.binned_signal.freq.shape + assert pair.snr.size > 1 + + def test_the_noise_spectrum_comes_back_too(self, pair: Any) -> None: + assert pair.binned_noise.amp.shape == pair.binned_signal.amp.shape + + def test_a_band_is_derivable_at_any_threshold(self, pair: Any) -> None: + """What a consumer does with the curve, and cannot do with a band.""" + for threshold in (2.0, 3.0, 5.0): + admitted = pair.snr >= threshold + assert admitted.dtype == bool + + +class TestErrorsAreTyped: + @pytest.mark.parametrize( + ("data", "reason"), + [ + (np.array([np.nan, 1.0, 2.0]), "non-finite"), + (np.array([1.0]), "too short"), + (np.zeros((4, 4)), "not 1-D"), + ], + ) + def test_a_bad_record_is_invalid_input(self, data: np.ndarray, reason: str) -> None: + with pytest.raises(InvalidInputError): + api.estimate_spectrum(data, 0.01, estimator="fft") + + def test_an_unknown_estimator_is_invalid_input(self) -> None: + with pytest.raises(InvalidInputError): + api.estimate_spectrum(np.arange(64.0), 0.01, estimator="nope") + + def test_an_unknown_weight_method_is_invalid_input(self, pair: Any) -> None: + with pytest.raises(InvalidInputError): + api.fit_spectrum(pair, weight_method="sideways") + + def test_every_typed_error_is_a_specmod_error(self) -> None: + """One `except` clause has to be able to catch all of them.""" + assert issubclass(InvalidInputError, SpecModError) + assert issubclass(MissingBackendError, SpecModError) + assert issubclass(api.InternalError, SpecModError) + + def test_they_still_look_like_the_builtins_they_replace(self) -> None: + """Existing callers catch ValueError and ImportError. They keep working + rather than being broken by the introduction of the hierarchy.""" + assert issubclass(InvalidInputError, ValueError) + assert issubclass(MissingBackendError, ImportError) diff --git a/tests/test_end_to_end.py b/tests/test_end_to_end.py index fad1017..b187cd4 100644 --- a/tests/test_end_to_end.py +++ b/tests/test_end_to_end.py @@ -12,6 +12,7 @@ from __future__ import annotations +import warnings from typing import Any import numpy as np @@ -141,12 +142,19 @@ def measured() -> Any: def recovered(measured: Any) -> Any: """The fit table. - Fitted on **velocity**, which is the pipeline's convention: the model - carries a motion factor, so ``llpsp`` is the displacement plateau either - way. Fitting the displacement set instead makes `initial_guess` useless — - it takes the spectral peak as the ``fc`` guess, and a displacement spectrum - falls monotonically, so the guess lands at the low band edge and the fit - settles near it. + Fitted in the motion the sensor recorded — velocity. The model carries a + motion factor, so ``llpsp`` is the displacement plateau whichever domain is + fitted, but converting first is not a neutral change of view: integrating + to displacement implicitly low-passes, and differentiating to acceleration + amplifies high-frequency noise. + + Velocity is also the convenient domain: it peaks at ``fc``, so + ``initial_guess`` taking the spectral peak is exact rather than + approximate. That holds for any omega-squared source — both registered + models — since the stationary point of ``f * [1 + (f/fc)**(g*n)]**(-1/g)`` + sits at ``f = fc`` whenever ``n == 2``, whatever the corner sharpness + ``g``. In displacement the spectrum is monotonic, so the peak is whichever + band edge it was handed. """ from specmod.fitting import FitSpectra # noqa: PLC0415 @@ -259,3 +267,67 @@ def test_the_target_falls_as_f_squared_above_the_corner(self) -> None: freq = np.array([40.0, 80.0]) amp = _brune_displacement_fas(freq) * np.exp(np.pi * freq * TRUE_TSTAR) assert amp[0] / amp[1] == pytest.approx(4.0, rel=0.05) + + +class TestFittingTheWrongMotionWarns: + """The one silent failure this exercise turned up, now audible. + + `FitSpectra(spectra.to_motion("displacement"))` is a natural thing to + write and used to return `fc` 1.6 against a true 8.0 with nothing said. + """ + + def test_a_displacement_spectrum_warns(self, measured: Any) -> None: + from specmod.fitting import initial_guess # noqa: PLC0415 + + with pytest.warns(UserWarning, match="the corner only in velocity"): + initial_guess(measured.to_motion("displacement")) + + def test_the_warning_names_the_station_and_the_motion(self, measured: Any) -> None: + from specmod.fitting import initial_guess # noqa: PLC0415 + + # One warning per station, not one for the set — otherwise a run over + # a mixed collection names whichever spectrum happened to be first. + # Asserting on all three also stops the two that `match=` does not + # select being re-emitted into pytest's warning summary. + with pytest.warns(UserWarning, match="is in displacement") as records: + initial_guess(measured.to_motion("displacement")) + + messages = [str(record.message) for record in records] + assert len(messages) == 3 + for station in ("XX.S00..HHN", "XX.S01..HHN", "XX.S02..HHN"): + assert any( + f"{station} is in displacement" in message for message in messages + ), f"nothing warned about {station}" + + def test_velocity_is_silent(self, measured: Any) -> None: + from specmod.fitting import initial_guess # noqa: PLC0415 + + with warnings.catch_warnings(): + warnings.simplefilter("error") + guesses = initial_guess(measured) + assert len(guesses) == 3 + + def test_the_guess_it_warns_about_is_the_one_that_ruins_the_fit( + self, measured: Any + ) -> None: + """Why it is worth a warning: the guess lands on the low band edge. + + The contrast is the point, not the precision. Peak-equals-corner is + exact for the noiseless model; on a measured spectrum the argmax + wanders, so velocity gives a guess in the right neighbourhood — a + starting point, which is all it has to be. Displacement gives one an + order of magnitude out, and the fit does not recover from it. + """ + from specmod.fitting import initial_guess # noqa: PLC0415 + + with pytest.warns(UserWarning, match="the corner only in velocity"): + displacement = initial_guess(measured.to_motion("displacement")) + velocity = initial_guess(measured) + + for id, guess in displacement.items(): + # Measured on these three stations: displacement guesses 0.71, + # 2.12 and 1.02 against velocity's 6.47, 7.83 and 5.72, for a true + # 8.0. Bounding each side is the contrast; a ratio between them + # would only be a brittle way of saying the same. + assert guess["fc"] < TRUE_FC / 3.0 + assert velocity[id]["fc"] == pytest.approx(TRUE_FC, rel=0.35) diff --git a/tests/test_golden_reference.py b/tests/test_golden_reference.py index 391a9d7..c0a414a 100644 --- a/tests/test_golden_reference.py +++ b/tests/test_golden_reference.py @@ -104,28 +104,66 @@ class from ``spectral.py`` into ``core/`` is checked here, on 28 real windows #: If a runner ever fails on tolerance alone, that bound is what to revisit. RTOL = 1e-6 -#: ``cwt`` alone is not exactly reproducible across machines, and this is an -#: open question rather than a tolerance that was tuned until green. +#: ``cwt`` was held at 5e-2 for a residual disagreement on CI that nothing +#: could explain. A round of measurement (below) failed to reproduce it and +#: ruled out every mechanism proposed for it, so the tolerance is being +#: tightened by four orders of magnitude to find out whether it still exists. #: -#: What is known. Four discontinuities were found and fixed (see -#: ``docs/REFACTOR_PLAN.md`` §4.5.2); they were worth 41-82% and this residual -#: is 1-2% on the quantile profile of 4 of 28 stations, with sums agreeing to -#: 1e-4 and array lengths identical. The other four estimators are exact -#: everywhere, and ``cwt``'s *signal* amplitudes and frequency axis are exact -#: too — it is only the post-rotation noise that moves. +#: **What was measured**, on a Linux box matching the reference environment +#: exactly — same system, arch, Python 3.11, numpy 2.4.6, scipy 1.17.1: #: -#: What is ruled out. Not a library-version effect: the failing runner matches -#: the reference machine exactly, down to Python 3.11.15, numpy 2.4.6, scipy -#: 1.17.1, obspy 1.5.0 and x86_64. Not a remaining branch in the lift, as far -#: as sweeping across the "already touching" boundary and across the centroid -#: split can show. Not local amplification: perturbing the input by 1e-15 -#: moves ``cwt`` noise by 2e-13, *better* than multitaper. +#: - ``cwt`` reproduces the committed reference to **3.8e-16**, machine +#: epsilon. The disagreement is not "Linux versus macOS"; it is that runner +#: versus this one. +#: - The ``cwt`` noise path is **linear** in the input: a 1e-13 perturbation +#: moves the noise by 8.7e-14, with no step, on all 28 windows. So a 1-2% +#: output difference needs a 1-2% *input* difference — which is not +#: something last-bit floating point can produce. +#: - Not bin-edge fragility: the closest ``cwt`` sample sits 5.7e-4 of a bin +#: from an interior edge, against 1.4e-6 for ``fft``. +#: - Not a differing window. One sample fewer moves ``fft``'s noise by 8.5% +#: and ``cwt``'s by 3.6%, so a window that differed on CI would show up in +#: ``fft`` *first* — and ``fft`` agrees exactly. +#: - Not quantile fragility from ``cwt``'s shorter arrays (51 samples against +#: 109): a 1e-15 perturbation moves its worst quantile by 2.4e-14, better +#: than ``fft``'s 5.1e-13. +#: - Not PyWavelets, which the estimator does not use, and not threading: +#: the transform is a batched ``numpy.fft.ifft``. #: -#: What is not known. Why macOS agrees with the reference and Linux does not, -#: when the reference was generated on Linux. That is the thread to pull next. -#: Until then this is loose enough to pass and tight enough that the 41-82% -#: class of defect cannot come back unnoticed. -RTOL_BY_ESTIMATOR = {"cwt": 5e-2} +#: So the CWT is as numerically stable here as every other estimator, and +#: 5e-2 was 12 orders of magnitude looser than anything measurable. Holding it +#: there hides whatever the real difference was rather than describing it. +#: +#: **1e-3 was an experiment**, not a calibration, and it has now returned its +#: answer: the residual is real. What one CI run showed, across six test jobs +#: on the same commit: +#: +#: - It fails on **ubuntu 3.11 and 3.13** with the *same* 8 differences, the +#: same three windows, and the same magnitudes to three significant figures. +#: Deterministic, not flaky. +#: - It **passes on ubuntu 3.12** and on macOS 3.11, 3.12 and 3.13, in that +#: same run. So it is not the OS, not the Python version, and not a package +#: version that tracks the Python version — it is which machine the job +#: landed on, which is what "that runner, not Linux" above suspected and +#: this is the same-run control for. +#: - On the machine that disagrees, ``fft``, ``welch``, ``multitaper`` and +#: ``quadratic`` all still reproduce the reference exactly. Whatever it is, +#: it is in the CWT path and not in that machine's arithmetic generally. +#: - Worst observed: 1.44e-2, on ``UR.AQ10.00.HHN bsnr``. Three of 28 windows +#: move at all (``LV.L007..HHN``, ``UR.AQ01.00.HHE``, ``UR.AQ10.00.HHN``). +#: +#: **2e-2 is a bound on that, not an explanation of it.** It is the worst +#: observed difference with about 40% of headroom, and 2.5 times tighter than +#: the 5e-2 it replaces — which is the most that can honestly be claimed while +#: the mechanism is still unidentified and unreproducible on any box available +#: to work on. A real regression in the CWT noise path smaller than 2% would +#: pass here; the four estimators held at 1e-6 are what covers the pipeline +#: those windows share. +#: +#: To take it further, the next measurement needs the failing machine: run the +#: cwt path on it against a passing one, on the three named windows, and diff +#: the noise arrays before binning rather than after. +RTOL_BY_ESTIMATOR = {"cwt": 2e-2} QUANTILES = np.linspace(0.0, 1.0, 33) diff --git a/tests/test_pick_readers.py b/tests/test_pick_readers.py index 4de205a..5a65933 100644 --- a/tests/test_pick_readers.py +++ b/tests/test_pick_readers.py @@ -162,6 +162,12 @@ def test_a_catalog_needs_no_reader(self) -> None: assert len(pk.read(_catalog())[0]) == 4 +@pytest.mark.filterwarnings( + # ObsPy's NORDIC writer, on a pick that carries no evaluation mode. It is + # third-party, it says nothing about what these tests assert, and it is + # matched by message so a different NORDIC warning would still surface. + "ignore:Evaluation mode None is not mappable:UserWarning" +) class TestRoster: """What §4.9.6 claims, measured rather than assumed. diff --git a/tests/test_release_config.py b/tests/test_release_config.py new file mode 100644 index 0000000..38bc4ba --- /dev/null +++ b/tests/test_release_config.py @@ -0,0 +1,167 @@ +"""The release configuration must agree with how the version is derived. + +Two files decide a release between them, and nothing makes them talk: +``release-please-config.json`` decides what the tag is *called*, and +``pyproject.toml`` decides which tags ``hatch-vcs`` will read (see +``test_versioning.py``). If they disagree the failure is quiet and expensive — +the tag is created, the release is published, and the wheel built from it +carries the fallback version ``0.1.1.postN.devN`` instead of ``0.2.0``. PyPI +will accept that upload and will not let the filename be reused. + +So the coupling is asserted here rather than left to be discovered once. +``tools/check_built_version.py`` is the same check made again at publish time, +on the artefact instead of the configuration, and its logic is tested below. +""" + +from __future__ import annotations + +import importlib.util +import json +import re +import sys +import tomllib +from pathlib import Path +from types import ModuleType + +import pytest + +ROOT = Path(__file__).resolve().parent.parent +CONFIG = ROOT / "release-please-config.json" +MANIFEST = ROOT / ".release-please-manifest.json" +WORKFLOW = ROOT / "ci" / "workflows" / "release.yml" + + +@pytest.fixture(scope="module") +def package_config() -> dict: + """The root package's entry — this is a single-package manifest.""" + packages = json.loads(CONFIG.read_text())["packages"] + assert list(packages) == ["."], "one package, at the repository root" + return packages["."] + + +@pytest.fixture(scope="module") +def tag_regex() -> re.Pattern[str]: + pyproject = tomllib.loads((ROOT / "pyproject.toml").read_text()) + return re.compile(pyproject["tool"]["hatch"]["version"]["raw-options"]["tag_regex"]) + + +def _load_check_built_version() -> ModuleType: + """Import ``tools/check_built_version.py``, which is a script rather than + a package member and so is not importable by name.""" + path = ROOT / "tools" / "check_built_version.py" + spec = importlib.util.spec_from_file_location("check_built_version", path) + assert spec is not None + assert spec.loader is not None + module = importlib.util.module_from_spec(spec) + sys.modules["check_built_version"] = module + spec.loader.exec_module(module) + return module + + +@pytest.fixture(scope="module") +def checker() -> ModuleType: + return _load_check_built_version() + + +class TestTheTagFormat: + def test_the_component_stays_out_of_the_tag(self, package_config: dict) -> None: + """With it left on, release-please tags ``specmod-v0.2.0``, which + pyproject's ``--match v[0-9]*`` does not describe and its ``tag_regex`` + does not parse. The default differs between release-please's own + modes, so it is set explicitly rather than relied on.""" + assert package_config["include-component-in-tag"] is False + + @pytest.mark.parametrize("version", ["0.2.0", "0.10.0", "1.0.0", "0.2.0-rc1"]) + def test_the_tag_release_please_will_create_is_one_hatch_vcs_reads( + self, tag_regex: re.Pattern[str], version: str + ) -> None: + # `include-component-in-tag: false` means the tag is exactly this. + tag = f"v{version}" + match = tag_regex.match(tag) + assert match is not None, f"hatch-vcs would not parse the tag {tag}" + assert match.group("version") == version + + +class TestTheVersionItWillPropose: + def test_the_manifest_names_the_same_single_package(self) -> None: + assert list(json.loads(MANIFEST.read_text())) == ["."] + + def test_the_recorded_version_is_a_release_version(self) -> None: + """release-please computes the next version from this one, so it has + to be a version and not a description of one.""" + version = json.loads(MANIFEST.read_text())["."] + assert re.fullmatch(r"\d+\.\d+\.\d+", version), version + + def test_a_breaking_change_cannot_mint_1_0_0_by_itself( + self, package_config: dict + ) -> None: + """The plan's §6.5 keeps the project on 0.x for the whole refactor. + There are already two `!` commits in the history, and the default + behaviour would read either as a 1.0.0 — a version that says the API + has stopped moving, minted with a DOI that cannot be retracted. + + This only matters below 1.0, which is exactly where the manifest is. + """ + major = int(json.loads(MANIFEST.read_text())["."].split(".")[0]) + if major >= 1: + pytest.skip("past 1.0; a breaking change should bump the major") + assert package_config["bump-minor-pre-major"] is True + + +class TestTheWorkflowUsesTheseFiles: + """A rename here is silent: release-please falls back to its defaults and + releases with none of the settings above.""" + + @pytest.mark.parametrize("path", [CONFIG, MANIFEST]) + def test_the_workflow_names_the_config_files(self, path: Path) -> None: + assert path.name in WORKFLOW.read_text() + + def test_the_publish_step_runs_the_built_version_check(self) -> None: + assert "tools/check_built_version.py" in WORKFLOW.read_text() + + +class TestTheBuiltVersionCheck: + """``tools/check_built_version.py`` runs between `uv build` and the + upload, where a wrong answer is permanent.""" + + def _dist(self, tmp_path: Path, *names: str) -> Path: + dist = tmp_path / "dist" + dist.mkdir() + for name in names: + (dist / name).touch() + return dist + + def test_a_matching_pair_passes(self, checker: ModuleType, tmp_path: Path) -> None: + dist = self._dist( + tmp_path, "specmod-0.2.0-py3-none-any.whl", "specmod-0.2.0.tar.gz" + ) + assert checker.check("v0.2.0", dist) == [] + + def test_the_fallback_version_is_caught( + self, checker: ModuleType, tmp_path: Path + ) -> None: + """What a tag `tag_regex` cannot parse actually produces.""" + dist = self._dist(tmp_path, "specmod-0.1.1.post1.dev4-py3-none-any.whl") + problems = checker.check("v0.2.0", dist) + assert len(problems) == 1 + assert "0.1.1.post1.dev4" in problems[0] + + def test_one_bad_artefact_among_good_ones_is_caught( + self, checker: ModuleType, tmp_path: Path + ) -> None: + dist = self._dist( + tmp_path, "specmod-0.2.0-py3-none-any.whl", "specmod-0.1.1.tar.gz" + ) + assert len(checker.check("v0.2.0", dist)) == 1 + + def test_an_empty_dist_is_a_failure_not_a_pass( + self, checker: ModuleType, tmp_path: Path + ) -> None: + """Nothing to compare must not read as nothing wrong.""" + assert checker.check("v0.2.0", self._dist(tmp_path)) != [] + + def test_a_tag_without_the_v_is_refused( + self, checker: ModuleType, tmp_path: Path + ) -> None: + dist = self._dist(tmp_path, "specmod-0.2.0-py3-none-any.whl") + assert checker.check("0.2.0", dist) != [] diff --git a/tests/test_utils.py b/tests/test_utils.py index 237bfcb..eef0336 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -211,13 +211,18 @@ def test_it_returns_a_copy_not_the_input(self) -> None: assert [tr.stats["repi"] for tr in st] == [30.0, 10.0], "input was mutated" def test_a_stream_without_distances_comes_back_unsorted(self) -> None: + """And says so. It used to `print` the fact, which is invisible to a + caller capturing logs and corrupts a CLI writing to a pipe — the + stream coming back in input order is otherwise indistinguishable from + a stream that was already in distance order.""" st = obspy.Stream( [ obspy.Trace(np.zeros(10), header={"station": "B"}), obspy.Trace(np.zeros(10), header={"station": "A"}), ] ) - got = ut.stream_distance_sort(st) + with pytest.warns(UserWarning, match="returned unsorted"): + got = ut.stream_distance_sort(st) assert [tr.stats.station for tr in got] == ["B", "A"] diff --git a/tools/check_built_version.py b/tools/check_built_version.py new file mode 100644 index 0000000..3f464c7 --- /dev/null +++ b/tools/check_built_version.py @@ -0,0 +1,84 @@ +#!/usr/bin/env python +"""Check that the distributions in ``dist/`` carry the version of the tag. + +Run as:: + + python tools/check_built_version.py v0.2.0 dist + +The version is never written down: ``hatch-vcs`` derives it from +``git describe``, filtered by the ``tag_regex`` in ``pyproject.toml``. That +makes one silent failure possible — a tag the regex does not match leaves the +build falling back to ``.postN.devN``, which is a perfectly valid +version string and would be uploaded under that name. PyPI does not allow a +filename to be reused, so a wrong upload is permanent. + +This runs between ``uv build`` and the upload, on the artefacts themselves +rather than on the configuration that produced them. Standard library only: +the publish job has no virtualenv. +""" + +from __future__ import annotations + +import re +import sys +from pathlib import Path + +#: ``specmod-0.2.0-py3-none-any.whl`` and ``specmod-0.2.0.tar.gz``. Both put +#: the version in the second ``-``/``.``-delimited field of the stem, and both +#: are normalised by the build backend, so the comparison is exact. +_WHEEL = re.compile(r"^(?P[^-]+)-(?P[^-]+)-.*\.whl$") +_SDIST = re.compile(r"^(?P[^-]+)-(?P.+)\.tar\.gz$") + + +def versions_in(dist: Path) -> dict[str, str]: + """Map each distribution filename in ``dist`` to the version it declares.""" + found: dict[str, str] = {} + for path in sorted(dist.iterdir()): + for pattern in (_WHEEL, _SDIST): + match = pattern.match(path.name) + if match: + found[path.name] = match.group("version") + break + return found + + +def check(tag: str, dist: Path) -> list[str]: + """Return the problems found; an empty list means the artefacts are good.""" + if not tag.startswith("v"): + return [f"tag {tag!r} does not start with 'v', which pyproject requires"] + expected = tag[1:] + + found = versions_in(dist) + if not found: + return [f"no wheel or sdist in {dist}"] + + return [ + f"{filename} is version {version}, expected {expected} from tag {tag}" + for filename, version in found.items() + if version != expected + ] + + +def main(argv: list[str]) -> int: + if len(argv) != 3: + print(__doc__) + return 2 + + tag, dist = argv[1], Path(argv[2]) + problems = check(tag, dist) + if problems: + for problem in problems: + print(f"ERROR {problem}") + print( + "\nThe usual cause is a tag pyproject.toml's tag_regex does not " + "match, which leaves hatch-vcs on its fallback version." + ) + return 1 + + for filename, version in versions_in(dist).items(): + print(f"ok {filename} is {version}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main(sys.argv))