From 91de646fbdf75ceda3198eb17bfd0f74dc63cca8 Mon Sep 17 00:00:00 2001 From: Marine Denolle Date: Sun, 2 Aug 2026 10:10:44 +0200 Subject: [PATCH 1/4] Add return_cc option to run_pipeline for coherence-based error models run_pipeline(..., return_cc=True) returns (dvv, valid, cc) with the per-epoch stretching correlation coefficient, the input needed by uq_measurement.weaver_stretching_error. CC is collected for the fixed and moving references with the stretching estimator (the moving loop was already computing it and discarding it); NaN for other estimators and the inversion reference. Behavior-preserving: the default two-tuple return is unchanged, and CC-gating stays fixed-reference-only (the moving-reference CC is returned for error modelling but does not alter the valid mask). Verified against the golden expected-metrics suite (22 passed). Motivation: the noisepy-dvv-cloud pipeline needs the real per-epoch CC to replace a placeholder in its dvv_err_within (Weaver/Clarke) column. Co-Authored-By: Claude Fable 5 --- .claude/settings.json | 8 +++++ CHANGELOG.md | 11 +++++++ src/codameter/deviations.py | 58 +++++++++++++++++++++++++++++++------ tests/test_deviations.py | 40 +++++++++++++++++++++++-- 4 files changed, 106 insertions(+), 11 deletions(-) create mode 100644 .claude/settings.json diff --git a/.claude/settings.json b/.claude/settings.json new file mode 100644 index 0000000..cfaae4e --- /dev/null +++ b/.claude/settings.json @@ -0,0 +1,8 @@ +{ + "permissions": { + "allow": [ + "Skill(deep-research)", + "Skill(deep-research:*)" + ] + } +} diff --git a/CHANGELOG.md b/CHANGELOG.md index aaad0bc..673b40e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,17 @@ All notable changes to `codameter` will be documented in this file. The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project follows [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## Unreleased + +### Added + +- **`run_pipeline(..., return_cc=True)`** — optionally return the per-epoch + stretching correlation coefficient alongside `(dvv, valid)`, for + coherence-based error models (`uq_measurement.weaver_stretching_error`). + Works for fixed and moving references with the stretching estimator; NaN + otherwise. The default two-tuple return and all gating behavior are + unchanged (CC-gating remains fixed-reference-only). + ## 0.3.0 — 2026-07-27 ### Added diff --git a/src/codameter/deviations.py b/src/codameter/deviations.py index 7084197..9087cdf 100644 --- a/src/codameter/deviations.py +++ b/src/codameter/deviations.py @@ -31,6 +31,7 @@ import numpy as np from .synthetic_demo import ( + METHODS, YEAR_D, C, Synth, @@ -83,27 +84,51 @@ # --------------------------------------------------------------------------- # Run one pipeline configuration on a shared set of daily CCFs. # --------------------------------------------------------------------------- -def _moving_reference(name, ccfs, t, *, band, fs, window, ref_days=45, **kw): +def _moving_reference( + name, ccfs, t, *, band, fs, window, ref_days=45, collect_cc=False, **kw +): """Generic trailing-reference measurement for *any* estimator. A moving reference re-baselines each epoch against the previous ``ref_days`` — the deviation that erases slow trends (best_practices rule 7). + + With ``collect_cc=True``, also returns the per-epoch correlation + coefficient for estimators that produce one (stretching); NaN otherwise. """ ndays = ccfs.shape[0] out = np.full(ndays, np.nan) + cc_out = np.full(ndays, np.nan) for d in range(ref_days, ndays): ref = ccfs[d - ref_days : d].mean(axis=0) - val = measure(name, ccfs[d], ref, t, band=band, fs=fs, window=window, **kw) - out[d] = np.atleast_1d(val)[0] - return out + if collect_cc: + res = METHODS[name](ccfs[d], ref, t, band=band, fs=fs, window=window, **kw) + if isinstance(res, tuple): + out[d] = np.atleast_1d(res[0])[0] + cc_out[d] = np.atleast_1d(res[1])[0] + else: + out[d] = np.atleast_1d(res)[0] + else: + val = measure(name, ccfs[d], ref, t, band=band, fs=fs, window=window, **kw) + out[d] = np.atleast_1d(val)[0] + return (out, cc_out) if collect_cc else out -def run_pipeline(ccfs, t, fs, cfg, *, eps_max=0.05): +def run_pipeline(ccfs, t, fs, cfg, *, eps_max=0.05, return_cc=False): """Recover dv/v(t) under one processing configuration ``cfg``. Returns ``(dvv, valid)``: the per-day series and a boolean mask of epochs the pipeline actually produced (moving/inversion references have a warm-up gap; CC-gating drops low-coherence epochs). + + With ``return_cc=True``, returns ``(dvv, valid, cc)`` where ``cc`` is the + per-epoch stretching correlation coefficient — the input to coherence-based + error models such as :func:`codameter.uq_measurement.weaver_stretching_error`. + ``cc`` is NaN wherever the configuration does not produce one (non-stretching + estimators, the inversion reference, and warm-up epochs). + + CC-gating (``cfg["gate"]``) applies to the fixed reference only, as it + always has; the moving-reference CC is returned for error modelling but + does not change ``valid``. """ name = cfg["estimator"] band, window, k, ref = cfg["band"], cfg["window"], cfg["stack"], cfg["reference"] @@ -122,9 +147,21 @@ def run_pipeline(ccfs, t, fs, cfg, *, eps_max=0.05): name, stacked, reference, t, band=band, fs=fs, window=window, **extra ) elif ref == "moving": - dvv = _moving_reference( - name, stacked, t, band=band, fs=fs, window=window, **extra - ) + if name == "stretching (TS)": + dvv, cc = _moving_reference( + name, + stacked, + t, + band=band, + fs=fs, + window=window, + collect_cc=True, + **extra, + ) + else: + dvv = _moving_reference( + name, stacked, t, band=band, fs=fs, window=window, **extra + ) elif ref == "inversion": # Brenguier et al. (2014) joint inversion (stretching) dvv = measure_inversion(ccfs, t, band=band, fs=fs, window=window) else: @@ -132,9 +169,12 @@ def run_pipeline(ccfs, t, fs, cfg, *, eps_max=0.05): dvv = np.asarray(dvv, float) valid = np.isfinite(dvv) - if cfg.get("gate") and cc is not None: + if cfg.get("gate") and ref == "fixed" and cc is not None: keep = cc > 0.6 valid &= keep + if return_cc: + cc_arr = np.full(dvv.shape, np.nan) if cc is None else np.asarray(cc, float) + return dvv, valid, cc_arr return dvv, valid diff --git a/tests/test_deviations.py b/tests/test_deviations.py index f715dbb..f6577e1 100644 --- a/tests/test_deviations.py +++ b/tests/test_deviations.py @@ -4,10 +4,12 @@ import numpy as np import pytest - from codameter import deviations as D from codameter.synthetic_demo import ( - Synth, _days, daily_ccfs, volcano_truth, + Synth, + _days, + daily_ccfs, + volcano_truth, ) @@ -56,3 +58,37 @@ def test_multiverse_sobol_sums_sensible(): assert -1e-9 <= v <= 1.0 + 1e-9 # The pipeline spread is non-trivial (the whole point). assert np.nanstd(mv["rms"]) > 0 + + +class TestReturnCC: + def test_default_still_two_tuple(self, small_dataset): + s, days, truth, ccfs = small_dataset + out = D.run_pipeline(ccfs, s.t, s.fs, D.BASELINE) + assert len(out) == 2 + + def test_fixed_stretching_returns_cc(self, small_dataset): + s, days, truth, ccfs = small_dataset + dvv, valid, cc = D.run_pipeline(ccfs, s.t, s.fs, D.BASELINE, return_cc=True) + assert cc.shape == dvv.shape + # On a clean synthetic the coherence should be high wherever valid. + assert np.all(cc[valid] > 0.6) + # dvv/valid identical to the two-tuple call (return_cc is read-only). + dvv2, valid2 = D.run_pipeline(ccfs, s.t, s.fs, D.BASELINE) + np.testing.assert_array_equal(dvv, dvv2) + np.testing.assert_array_equal(valid, valid2) + + def test_moving_stretching_returns_cc_after_warmup(self, small_dataset): + s, days, truth, ccfs = small_dataset + cfg = dict(D.BASELINE, reference="moving") + dvv, valid, cc = D.run_pipeline(ccfs, s.t, s.fs, cfg, return_cc=True) + assert np.isnan(cc[:10]).all() # warm-up gap + assert np.isfinite(cc[valid]).all() + # Gating stays fixed-reference-only: valid must match the legacy call. + dvv2, valid2 = D.run_pipeline(ccfs, s.t, s.fs, cfg) + np.testing.assert_array_equal(valid, valid2) + + def test_non_stretching_cc_is_nan(self, small_dataset): + s, days, truth, ccfs = small_dataset + cfg = dict(D.BASELINE, estimator="MWCS", gate=False) + dvv, valid, cc = D.run_pipeline(ccfs, s.t, s.fs, cfg, return_cc=True) + assert np.isnan(cc).all() From b11c75cb1cf65487abf7195e4ed86a0d48fb88f7 Mon Sep 17 00:00:00 2001 From: Marine Denolle Date: Mon, 3 Aug 2026 07:57:46 +0200 Subject: [PATCH 2/4] Vectorize trailing stack, moving-reference stretching, and add filter-once ensembles Three performance fast paths, each verified to reproduce the per-day loop it replaces to ~1e-15 in dv/v (regression tests at rtol=0, atol=1e-12): - _trailing_stack: difference of float64 cumulative sums, O(ndays*nlag) independent of the stack length instead of O(ndays*k*nlag) (~2.3x at k=45). - measure_stretching_trailing: vectorized stretching against a trailing reference. The stretched sample positions t/(1+eps) are data-independent, so the linear-interpolation gather indices/weights are computed once per epsilon and applied to all days at once; trailing references come from a cumulative sum and the band-pass runs once over the whole matrix. deviations._moving_reference dispatches to it for "stretching (TS)" (measured 4.7x on the 3-year volcano synthetic, 12.2 -> 2.6 s), keeping the generic per-day loop for other estimators. collect_cc behavior from the return_cc branch is preserved. - run_pipeline(..., prefiltered=True): callers evaluating several stack/reference variants at the same band can band-pass the raw CCF matrix once; the estimators skip their internal band-pass. Exact because the band-pass is linear and commutes with linear stacking; restricted to the estimators whose band usage is that one linear filter (stretching, WCC, DTW, MWCS), ValueError otherwise. Combined, a 5-member same-band ensemble drops ~4x per pair-band. Co-Authored-By: Claude Fable 5 --- CHANGELOG.md | 23 +++++++ src/codameter/deviations.py | 47 +++++++++++++- src/codameter/synthetic_demo.py | 106 ++++++++++++++++++++++++++++---- tests/test_deviations.py | 76 +++++++++++++++++++++++ 4 files changed, 237 insertions(+), 15 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 673b40e..d339fc6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,29 @@ and this project follows [Semantic Versioning](https://semver.org/spec/v2.0.0.ht Works for fixed and moving references with the stretching estimator; NaN otherwise. The default two-tuple return and all gating behavior are unchanged (CC-gating remains fixed-reference-only). +- **`run_pipeline(..., prefiltered=True)`** — accept CCFs already band-passed + at `cfg["band"]` and skip the estimators' internal band-pass, so callers + evaluating several stack/reference variants at the same band filter the raw + matrix once. Exact to float rounding because the band-pass is linear and + commutes with linear stacking; only valid at an identical band and only for + the estimators whose band usage is that one linear filter (stretching, WCC, + DTW, MWCS — the wavelet estimators raise). +- **`measure_stretching_trailing`** — vectorized stretching against a trailing + (moving) reference. The stretched sample positions `t/(1+eps)` are + data-independent, so the interpolation gather indices/weights are computed + once per epsilon and applied to all days at once; trailing references come + from a cumulative sum and the band-pass runs once over the whole matrix. + `deviations._moving_reference` dispatches to it for the stretching + estimator (~4.9x on the 3-year volcano synthetic), keeping the generic + per-day loop for the other estimators. + +### Changed + +- **`_trailing_stack`** is now a difference of float64 cumulative sums — + O(ndays x nlag) independent of the stack length instead of + O(ndays x k x nlag) (~2.3x at k=45). All three fast paths reproduce the + replaced per-day loops to ~1e-15 in dv/v, enforced by regression tests at + atol=1e-12; combined, a 5-member same-band ensemble drops ~4x in runtime. ## 0.3.0 — 2026-07-27 diff --git a/src/codameter/deviations.py b/src/codameter/deviations.py index 9087cdf..817cd2f 100644 --- a/src/codameter/deviations.py +++ b/src/codameter/deviations.py @@ -41,6 +41,7 @@ measure, measure_inversion, measure_stretching, + measure_stretching_trailing, volcano_truth, ) @@ -94,7 +95,17 @@ def _moving_reference( With ``collect_cc=True``, also returns the per-epoch correlation coefficient for estimators that produce one (stretching); NaN otherwise. + + Stretching dispatches to the vectorized + :func:`codameter.synthetic_demo.measure_stretching_trailing` fast path + (identical to float rounding, ~5x faster); the generic per-day loop below + serves every other estimator. """ + if name == "stretching (TS)": + out, cc_out = measure_stretching_trailing( + ccfs, t, band=band, fs=fs, window=window, ref_days=ref_days, **kw + ) + return (out, cc_out) if collect_cc else out ndays = ccfs.shape[0] out = np.full(ndays, np.nan) cc_out = np.full(ndays, np.nan) @@ -113,7 +124,12 @@ def _moving_reference( return (out, cc_out) if collect_cc else out -def run_pipeline(ccfs, t, fs, cfg, *, eps_max=0.05, return_cc=False): +# Estimators whose only use of the band is one linear band-pass of the input +# waveforms, so a caller may apply that band-pass once and skip it here. +_PREFILTER_OK = {"stretching (TS)", "WCC", "DTW", "MWCS"} + + +def run_pipeline(ccfs, t, fs, cfg, *, eps_max=0.05, return_cc=False, prefiltered=False): """Recover dv/v(t) under one processing configuration ``cfg``. Returns ``(dvv, valid)``: the per-day series and a boolean mask of epochs the @@ -129,18 +145,41 @@ def run_pipeline(ccfs, t, fs, cfg, *, eps_max=0.05, return_cc=False): CC-gating (``cfg["gate"]``) applies to the fixed reference only, as it always has; the moving-reference CC is returned for error modelling but does not change ``valid``. + + With ``prefiltered=True``, ``ccfs`` is taken as already band-passed at + ``cfg["band"]`` and the estimator skips its internal band-pass. Callers + that evaluate several stack/reference variants at the *same* band can + band-pass the raw CCF matrix once and share it. This is exact (to float + rounding) because the band-pass is linear, so it commutes with the linear + stacking that builds trailing stacks and references — it is only valid at + an identical band and only for the estimators whose band usage is that one + linear filter (stretching, WCC, DTW, MWCS; the wavelet estimators apply no + such filter, so ``prefiltered`` raises for them). """ name = cfg["estimator"] band, window, k, ref = cfg["band"], cfg["window"], cfg["stack"], cfg["reference"] + if prefiltered and name not in _PREFILTER_OK: + raise ValueError( + f"prefiltered=True is only valid for {sorted(_PREFILTER_OK)}, not {name!r}" + ) stacked = _trailing_stack(ccfs, k) extra = {"eps_max": eps_max} if name in ("stretching (TS)", "WTS") else {} + if prefiltered: + extra["prefiltered"] = True cc = None if ref == "fixed": reference = ccfs[: int(0.6 * len(ccfs))].mean(axis=0) # long stable stack if name == "stretching (TS)": dvv, cc = measure_stretching( - stacked, reference, t, band=band, fs=fs, window=window, eps_max=eps_max + stacked, + reference, + t, + band=band, + fs=fs, + window=window, + eps_max=eps_max, + prefiltered=prefiltered, ) else: dvv = measure( @@ -163,7 +202,9 @@ def run_pipeline(ccfs, t, fs, cfg, *, eps_max=0.05, return_cc=False): name, stacked, t, band=band, fs=fs, window=window, **extra ) elif ref == "inversion": # Brenguier et al. (2014) joint inversion (stretching) - dvv = measure_inversion(ccfs, t, band=band, fs=fs, window=window) + dvv = measure_inversion( + ccfs, t, band=band, fs=fs, window=window, prefiltered=prefiltered + ) else: raise ValueError(ref) diff --git a/src/codameter/synthetic_demo.py b/src/codameter/synthetic_demo.py index bf307e7..554437f 100644 --- a/src/codameter/synthetic_demo.py +++ b/src/codameter/synthetic_demo.py @@ -345,6 +345,7 @@ def stretching_cc( branch: str = "both", eps_max: float = 0.06, n_eps: int = 161, + prefiltered: bool = False, ) -> tuple[np.ndarray, np.ndarray]: """The full correlation-coefficient image ``CC(epsilon, time)``. @@ -352,14 +353,19 @@ def stretching_cc( that aggregation workflows either reduce to a per-trace peak *before* averaging, or average *as images* before peak-picking (see :func:`peak_dvv` and the aggregation demo). + + With ``prefiltered=True``, ``cur_mat`` and ``ref`` are taken as already + band-passed at ``band`` and the internal band-pass is skipped (valid + because band-passing is linear and commutes with the linear stacking that + builds references — see :func:`codameter.deviations.run_pipeline`). """ cur_mat = np.atleast_2d(cur_mat) - reff = bandpass(ref, fs, *band) + reff = ref if prefiltered else bandpass(ref, fs, *band) es = np.linspace(-eps_max, eps_max, n_eps) sel = _window_mask(t, window, branch) trials = np.stack([np.interp(t / (1.0 + e), t, reff)[sel] for e in es]) trials = trials / (np.linalg.norm(trials, axis=1, keepdims=True) + 1e-12) - curf = bandpass(cur_mat, fs, *band)[:, sel] + curf = (cur_mat if prefiltered else bandpass(cur_mat, fs, *band))[:, sel] curf = curf / (np.linalg.norm(curf, axis=1, keepdims=True) + 1e-12) return es, curf @ trials.T # [ndays, n_eps] @@ -384,13 +390,15 @@ def measure_stretching( branch: str = "both", eps_max: float = 0.06, n_eps: int = 161, + prefiltered: bool = False, ) -> tuple[np.ndarray, np.ndarray]: """Stretching dv/v: grid-search the stretch maximizing windowed correlation. ``ref`` is a single reference vector (fixed-reference scheme). ``branch`` selects the causal, acausal, or both coda branches — measuring the two branches separately is the standard clock-error diagnostic. Returns the - per-day dv/v and the peak correlation coefficient. + per-day dv/v and the peak correlation coefficient. ``prefiltered`` is + forwarded to :func:`stretching_cc`. """ es, cc = stretching_cc( cur_mat, @@ -402,6 +410,7 @@ def measure_stretching( branch=branch, eps_max=eps_max, n_eps=n_eps, + prefiltered=prefiltered, ) return peak_dvv(es, cc) @@ -433,6 +442,68 @@ def measure_stretching_moving( return out +def measure_stretching_trailing( + cur_mat: np.ndarray, + t: np.ndarray, + *, + band: tuple[float, float], + fs: float, + window: tuple[float, float], + ref_days: int = 45, + branch: str = "both", + eps_max: float = 0.06, + n_eps: int = 161, + prefiltered: bool = False, +) -> tuple[np.ndarray, np.ndarray]: + """Vectorized stretching against a trailing reference (previous ``ref_days``). + + Numerically equivalent (to float rounding, ~1e-15 in dv/v) to calling + :func:`measure_stretching` day by day against + ``cur_mat[d - ref_days : d].mean(axis=0)``, but ~5x faster: the stretched + sample positions ``t / (1 + eps)`` are data-independent, so the + linear-interpolation gather indices and weights are computed once per + epsilon and applied to every day's band-passed trailing reference at once. + The trailing references are built as a difference of cumulative sums and + the band-pass runs over the whole matrix in one FFT. + + Returns ``(dvv, cc)`` over the full length of ``cur_mat``; the ``ref_days`` + warm-up epochs are NaN. + """ + cur_mat = np.atleast_2d(np.asarray(cur_mat, float)) + ndays, nlag = cur_mat.shape + dvv = np.full(ndays, np.nan) + cc_peak = np.full(ndays, np.nan) + if ndays <= ref_days: + return dvv, cc_peak + # Trailing reference for day d is the mean of rows [d - ref_days, d). + csum = np.cumsum(cur_mat, axis=0, dtype=np.float64) + head = csum[ref_days - 1 : ndays - 1] + tail = np.concatenate([np.zeros((1, nlag)), csum[: ndays - ref_days - 1]], axis=0) + refs = (head - tail) / float(ref_days) + reffs = refs if prefiltered else bandpass(refs, fs, *band) + + sel = _window_mask(t, window, branch) + tsel = t[sel] + curf = ( + cur_mat[ref_days:] if prefiltered else bandpass(cur_mat[ref_days:], fs, *band) + )[:, sel] + curf = curf / (np.linalg.norm(curf, axis=1, keepdims=True) + 1e-12) + + es = np.linspace(-eps_max, eps_max, n_eps) + cc_img = np.empty((ndays - ref_days, n_eps)) + for ei, e in enumerate(es): + # Gather indices/weights of np.interp(t / (1 + e), t, .) on the window, + # clamped at the grid ends exactly as np.interp clamps. + q = tsel / (1.0 + e) + j = np.clip(np.searchsorted(t, q, side="right") - 1, 0, t.size - 2) + w = np.clip((q - t[j]) / (t[j + 1] - t[j]), 0.0, 1.0) + trials = reffs[:, j] * (1.0 - w) + reffs[:, j + 1] * w + trials = trials / (np.linalg.norm(trials, axis=1, keepdims=True) + 1e-12) + cc_img[:, ei] = np.einsum("ij,ij->i", curf, trials) + dvv[ref_days:], cc_peak[ref_days:] = peak_dvv(es, cc_img) + return dvv, cc_peak + + def measure_mwcs( cur_mat: np.ndarray, ref: np.ndarray, @@ -443,6 +514,7 @@ def measure_mwcs( window: tuple[float, float], subwin_s: float = 6.0, step_s: float = 3.0, + prefiltered: bool = False, ) -> np.ndarray: """MWCS dv/v: cross-spectral phase delay per sub-window, slope of dt vs lapse. @@ -455,8 +527,8 @@ def measure_mwcs( dv/v (e.g. pre-failure landslides) where stretching stays robust. """ cur_mat = np.atleast_2d(cur_mat) - reff = bandpass(ref, fs, *band) - curf = bandpass(cur_mat, fs, *band) + reff = ref if prefiltered else bandpass(ref, fs, *band) + curf = cur_mat if prefiltered else bandpass(cur_mat, fs, *band) centers = np.arange(window[0] + subwin_s / 2, window[1] - subwin_s / 2, step_s) half = int(round(subwin_s / 2 * fs)) taper = np.hanning(2 * half) @@ -500,6 +572,7 @@ def measure_wcc( window: tuple[float, float], subwin_s: float = 6.0, step_s: float = 3.0, + prefiltered: bool = False, ) -> np.ndarray: """WCC dv/v: time-domain windowed cross-correlation delay, slope of dt vs lapse. @@ -510,8 +583,8 @@ def measure_wcc( the seven estimators in NoisePy's ``monitoring_methods`` (``wcc_dvv``). """ cur_mat = np.atleast_2d(cur_mat) - reff = bandpass(ref, fs, *band) - curf = bandpass(cur_mat, fs, *band) + reff = ref if prefiltered else bandpass(ref, fs, *band) + curf = cur_mat if prefiltered else bandpass(cur_mat, fs, *band) centers = np.arange(window[0] + subwin_s / 2, window[1] - subwin_s / 2, step_s) half = int(round(subwin_s / 2 * fs)) taper = np.hanning(2 * half) @@ -579,6 +652,7 @@ def measure_dtw( fs: float, window: tuple[float, float], max_lag_s: float = 0.8, + prefiltered: bool = False, ) -> np.ndarray: """DTW dv/v: warp the current trace onto the reference, slope of lag vs lapse. @@ -587,8 +661,8 @@ def measure_dtw( changes (Yuan et al. 2021). NoisePy ``dtw_dvv``. """ cur_mat = np.atleast_2d(cur_mat) - reff = bandpass(ref, fs, *band) - curf = bandpass(cur_mat, fs, *band) + reff = ref if prefiltered else bandpass(ref, fs, *band) + curf = cur_mat if prefiltered else bandpass(cur_mat, fs, *band) sel = (t >= window[0]) & (t <= window[1]) # causal branch only tt = t[sel] max_lag = int(round(max_lag_s * fs)) @@ -875,6 +949,7 @@ def measure_inversion( block_days: int = 7, max_lag_blocks: int = 10, smooth: float = 5.0, + prefiltered: bool = False, ) -> np.ndarray: """Brenguier et al. (2014)-style joint inversion for a continuous dv/v series. @@ -903,6 +978,7 @@ def measure_inversion( window=window, eps_max=0.03, n_eps=81, + prefiltered=prefiltered, ) for i in range(j + 1, min(m, j + max_lag_blocks + 1)): rows += [eq, eq] @@ -1059,9 +1135,15 @@ def _yrs(days: np.ndarray) -> np.ndarray: def _trailing_stack(ccfs: np.ndarray, k: int) -> np.ndarray: if k <= 1: return ccfs - out = np.empty_like(ccfs) - for d in range(ccfs.shape[0]): - out[d] = ccfs[max(0, d - k + 1) : d + 1].mean(axis=0) + # Trailing mean of the last k days (shorter at the start), as a difference + # of float64 cumulative sums: O(ndays * nlag) instead of O(ndays * k * nlag). + ndays = ccfs.shape[0] + csum = np.cumsum(ccfs, axis=0, dtype=np.float64) + out = np.empty_like(csum) + out[:k] = csum[:k] + np.subtract(csum[k:], csum[:-k], out=out[k:]) # window sum over [d-k+1, d] + counts = np.minimum(np.arange(1, ndays + 1), k).astype(np.float64) + out /= counts[:, None] return out diff --git a/tests/test_deviations.py b/tests/test_deviations.py index f6577e1..f5fbf65 100644 --- a/tests/test_deviations.py +++ b/tests/test_deviations.py @@ -8,7 +8,10 @@ from codameter.synthetic_demo import ( Synth, _days, + _trailing_stack, + bandpass, daily_ccfs, + measure_stretching, volcano_truth, ) @@ -92,3 +95,76 @@ def test_non_stretching_cc_is_nan(self, small_dataset): cfg = dict(D.BASELINE, estimator="MWCS", gate=False) dvv, valid, cc = D.run_pipeline(ccfs, s.t, s.fs, cfg, return_cc=True) assert np.isnan(cc).all() + + +class TestFastPathRegressions: + """The vectorized fast paths must reproduce the per-day loops they replace.""" + + def test_trailing_stack_matches_per_day_loop(self, small_dataset): + s, days, truth, ccfs = small_dataset + for k in (1, 2, 10, 45, ccfs.shape[0] + 5): + fast = _trailing_stack(ccfs, k) + slow = np.stack( + [ + ccfs[max(0, d - k + 1) : d + 1].mean(axis=0) + for d in range(ccfs.shape[0]) + ] + ) + np.testing.assert_allclose(fast, slow, rtol=0, atol=1e-12) + + def test_moving_reference_matches_generic_loop(self, small_dataset): + s, days, truth, ccfs = small_dataset + band, window = D.BASELINE["band"], D.BASELINE["window"] + stacked = _trailing_stack(ccfs, D.BASELINE["stack"]) + fast, fast_cc = D._moving_reference( + "stretching (TS)", + stacked, + s.t, + band=band, + fs=s.fs, + window=window, + collect_cc=True, + eps_max=0.05, + ) + ndays = stacked.shape[0] + slow = np.full(ndays, np.nan) + slow_cc = np.full(ndays, np.nan) + for d in range(45, ndays): + ref = stacked[d - 45 : d].mean(axis=0) + v, c = measure_stretching( + stacked[d], ref, s.t, band=band, fs=s.fs, window=window, eps_max=0.05 + ) + slow[d], slow_cc[d] = v[0], c[0] + np.testing.assert_allclose(fast, slow, rtol=0, atol=1e-12) + np.testing.assert_allclose(fast_cc, slow_cc, rtol=0, atol=1e-12) + + @pytest.mark.parametrize( + "cfg", + [ + D.BASELINE, + dict(D.BASELINE, reference="moving"), + dict(D.BASELINE, reference="inversion"), + dict(D.BASELINE, estimator="MWCS"), + ], + ids=["fixed", "moving", "inversion", "mwcs"], + ) + def test_prefiltered_matches_internal_bandpass(self, small_dataset, cfg): + # Band-passing is linear, so filtering the raw CCFs once outside must + # equal the estimator's internal band-pass of every stack/reference. + s, days, truth, ccfs = small_dataset + filt = bandpass(ccfs, s.fs, *cfg["band"]) + dvv_a, val_a, cc_a = D.run_pipeline(ccfs, s.t, s.fs, cfg, return_cc=True) + dvv_b, val_b, cc_b = D.run_pipeline( + filt, s.t, s.fs, cfg, return_cc=True, prefiltered=True + ) + np.testing.assert_allclose(dvv_b, dvv_a, rtol=0, atol=1e-12) + np.testing.assert_array_equal(val_b, val_a) + np.testing.assert_allclose(cc_b, cc_a, rtol=0, atol=1e-12) + + def test_prefiltered_rejects_estimators_without_linear_bandpass( + self, small_dataset + ): + s, days, truth, ccfs = small_dataset + cfg = dict(D.BASELINE, estimator="WTS") + with pytest.raises(ValueError, match="prefiltered"): + D.run_pipeline(ccfs, s.t, s.fs, cfg, prefiltered=True) From 0e2e9ee323f8e48da479e0cf8c0615cdc11b8533 Mon Sep 17 00:00:00 2001 From: Marine Denolle Date: Mon, 3 Aug 2026 11:55:18 +0200 Subject: [PATCH 3/4] Make equal_nan explicit in the prefiltered-CC regression test Copilot flagged assert_allclose(cc_b, cc_a, ...) as failing on the all-NaN cc arrays (inversion/mwcs cases, where no CC is collected) because equal_nan defaults to False. That's true of plain np.allclose, but numpy.testing.assert_allclose already defaults equal_nan=True -- verified all 4 parametrized cases (fixed/moving/inversion/mwcs) were already passing. Not a real bug, but making the default explicit (with a comment on why) so the next reader doesn't hit the same false alarm. --- tests/test_deviations.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/tests/test_deviations.py b/tests/test_deviations.py index f5fbf65..289bcfa 100644 --- a/tests/test_deviations.py +++ b/tests/test_deviations.py @@ -159,7 +159,10 @@ def test_prefiltered_matches_internal_bandpass(self, small_dataset, cfg): ) np.testing.assert_allclose(dvv_b, dvv_a, rtol=0, atol=1e-12) np.testing.assert_array_equal(val_b, val_a) - np.testing.assert_allclose(cc_b, cc_a, rtol=0, atol=1e-12) + # cc is all-NaN for "inversion"/"mwcs" (no CC collected for those); + # equal_nan=True (assert_allclose's default, unlike plain np.allclose) + # is what makes that comparison pass -- kept explicit here. + np.testing.assert_allclose(cc_b, cc_a, rtol=0, atol=1e-12, equal_nan=True) def test_prefiltered_rejects_estimators_without_linear_bandpass( self, small_dataset From 5fbb1e0715d71278102db81c60a313d1e5635541 Mon Sep 17 00:00:00 2001 From: Marine Denolle Date: Mon, 3 Aug 2026 12:05:39 +0200 Subject: [PATCH 4/4] Drop .claude/settings.json from this PR (out of scope) Same issue as codameter#32: this file grants Skill(deep-research) permissions, unrelated to the vectorized fast-paths change, isn't tracked on master, and doesn't belong bundled in here. --- .claude/settings.json | 8 -------- 1 file changed, 8 deletions(-) delete mode 100644 .claude/settings.json diff --git a/.claude/settings.json b/.claude/settings.json deleted file mode 100644 index cfaae4e..0000000 --- a/.claude/settings.json +++ /dev/null @@ -1,8 +0,0 @@ -{ - "permissions": { - "allow": [ - "Skill(deep-research)", - "Skill(deep-research:*)" - ] - } -}