From 1a59aace94f5b5d19b04db4731594b6ae6ec1d8d Mon Sep 17 00:00:00 2001 From: Jung Dae Suh Date: Wed, 1 Jul 2026 13:31:10 -0400 Subject: [PATCH] fix(core): correct Crucible-confirmed QS + spectral physics defects MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - qs_bridge: fit-quality status used ok/warn/error, which the contract's Quality (good/warn/bad) and NodeView.QCOLOR can't color; route through contracts.quality_for_k (SSOT). - _slcontour/fit.py + core/quasistationary.py: per-coefficient SVD σ summed the wrong axis (axis=0 → 1/w_k² by singular position, not the covariance diagonal); axis=1 fixes the QS error bars and sizes fit_sigmas to n_cols so underdetermined fits no longer IndexError. - _slcontour/fit.py: 2D sinusoidal-integral basis divided by deg2rad(dx*dy)*n*m (π/180 once, no i²) → m≠0 columns off by -180/π; use the separable product. - core/quasistationary.py reconstruct_grid: used c instead of conj(c) vs its own forward model (and qs_bridge's sign) → toroidal mirror; use conj(c). - core/spectral.py: phase_sigma reported the inverse-variance (BLUE) variance while the intercept is amplitude/coherence-weighted; propagate σ through the same weights so the reported 1σ matches the returned estimator. Adds tests/test_crucible_regressions.py + a qs_bridge quality-vocab test, each failing against the pre-fix behavior. ruff + ty(src) clean; full suite green. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/magnetics/_slcontour/fit.py | 6 +- src/magnetics/core/qs_bridge.py | 16 ++-- src/magnetics/core/quasistationary.py | 11 ++- src/magnetics/core/spectral.py | 12 ++- tests/test_crucible_regressions.py | 110 ++++++++++++++++++++++++++ tests/test_qs_bridge.py | 10 +++ 6 files changed, 146 insertions(+), 19 deletions(-) create mode 100644 tests/test_crucible_regressions.py diff --git a/src/magnetics/_slcontour/fit.py b/src/magnetics/_slcontour/fit.py index 876bba6..f882873 100644 --- a/src/magnetics/_slcontour/fit.py +++ b/src/magnetics/_slcontour/fit.py @@ -88,7 +88,7 @@ def form_basis_function( return ( (np.exp(1j * m * np.deg2rad(y2)) - np.exp(1j * m * np.deg2rad(y1))) * (np.exp(1j * n * np.deg2rad(x2)) - np.exp(1j * n * np.deg2rad(x1))) - ) / (np.deg2rad(dx * dy) * n * m) + ) / (np.deg2rad(dx) * np.deg2rad(dy) * (1j * n) * (1j * m)) # ── gaussian-point ──────────────────────────────────────────────────────── if fit_basis == "gaussian-point": @@ -352,7 +352,9 @@ def _printv(*a): w_inv = 1.0 / w_a w_inv[~valid] = 0.0 w_inv = np.hstack((w_inv, [0.0] * max(Vh_a.shape[0] - w_a.shape[0], 0))) - fit_sigma2 = np.sum((Vh_a.T * w_inv) ** 2, axis=0) + # diag(cov) = diag(V W^-2 V^H) = Σ_k V[j,k]^2 w_inv[k]^2 — sum over the singular + # index k (axis=1 of Vh_a.T, shape [n_cols, k]), NOT over the coefficient index. + fit_sigma2 = np.sum((Vh_a.T * w_inv) ** 2, axis=1) fit_sigmas = np.sqrt(fit_sigma2) # (n_columns,) # ── reform per-mode complex (sinusoidal) or real (Gaussian) coefficients ── diff --git a/src/magnetics/core/qs_bridge.py b/src/magnetics/core/qs_bridge.py index 13358c0..b985e8c 100644 --- a/src/magnetics/core/qs_bridge.py +++ b/src/magnetics/core/qs_bridge.py @@ -55,14 +55,6 @@ def _mode_label(n: int | float, m: int | float) -> str: return f"n={n}" if m == 0 else f"m/n={m}/{n}" -def _quality_status(K: float) -> str: - if K > 20: - return "error" - if K > 10: - return "warn" - return "ok" - - def _reconstruct_grid( fit_ds, phi_grid: np.ndarray, theta_grid: np.ndarray, t_idx: int ) -> np.ndarray: @@ -257,8 +249,12 @@ def fit_to_fit_quality_node(fit_ds) -> dict: return contracts.metrics( title="fit quality", fields=[ - {"label": "K (raw)", "value": f"{K:.1f}", "status": _quality_status(K)}, - {"label": "K (eff)", "value": f"{eff_cn:.1f}", "status": _quality_status(eff_cn)}, + {"label": "K (raw)", "value": f"{K:.1f}", "status": contracts.quality_for_k(K)}, + { + "label": "K (eff)", + "value": f"{eff_cn:.1f}", + "status": contracts.quality_for_k(eff_cn), + }, {"label": "K cutoff", "value": f"{fit_cond:.0f}"}, {"label": "χ² (mean)", "value": f"{mean_chi2:.3f}"}, {"label": "channels", "value": n_ch}, diff --git a/src/magnetics/core/quasistationary.py b/src/magnetics/core/quasistationary.py index fb636c5..74e058e 100644 --- a/src/magnetics/core/quasistationary.py +++ b/src/magnetics/core/quasistationary.py @@ -170,7 +170,10 @@ def fit( # Per-coeff uncertainty from SVD pseudo-inverse w_inv = np.where(valid, 1.0 / np.where(w_a != 0, w_a, 1.0), 0.0) - fit_sigmas = np.sqrt(np.sum((Vh_a.T * w_inv) ** 2, axis=0)) # [n_cols] + # diag(cov) = diag(V W^-2 V^H) = Σ_k V[j,k]^2 w_inv[k]^2 — sum over the singular + # index k (axis=1 of Vh_a.T, shape [n_cols, k]). axis=1 also yields length n_cols + # for an underdetermined fit (k = n_ch < n_cols), so the reform loop can't overrun. + fit_sigmas = np.sqrt(np.sum((Vh_a.T * w_inv) ** 2, axis=1)) # [n_cols] # ── reform complex coefficients (one complex number per mode) ───────────── coeffs_c: list[np.ndarray] = [] @@ -250,7 +253,9 @@ def reconstruct_grid( z = np.zeros((len(theta_grid), len(phi_grid))) for i, (n, m) in enumerate(zip(result.ns, result.ms)): c = result.coeffs[i, t_idx] - # outer product: [n_theta, n_phi] + # The fit's forward model is A@x = x_r·cos ψ + x_i·sin ψ = Re(conj(c)·e^{iψ}) + # with c = x_r + i·x_i and ψ = mθ + nφ, so the reconstruction must use conj(c) + # (matching qs_bridge._reconstruct_grid); using c mirrors the map toroidally. basis = np.exp(1j * m * theta_rad)[:, None] * np.exp(1j * n * phi_rad)[None, :] - z += (c * basis).real + z += (np.conj(c) * basis).real return z diff --git a/src/magnetics/core/spectral.py b/src/magnetics/core/spectral.py index a2e5bb3..668c31c 100644 --- a/src/magnetics/core/spectral.py +++ b/src/magnetics/core/spectral.py @@ -800,9 +800,10 @@ def fit_toroidal_mode( best_n, best_R, best_c = candidates[best_j], float(rmag[best_j]), float(inter[best_j]) # Uncertainty (eigspec-style): propagate per-probe phase σ from the cross-spectral - # statistics into (a) the intercept's 1σ via inverse-variance combination and - # (b) a posterior over candidate n from each hypothesis's weighted χ². Falls back - # to unit σ when phase_error is absent so older callers still get sane numbers. + # statistics into (a) the intercept's 1σ by propagating σ through the SAME weights w + # that produced the intercept, and (b) a posterior over candidate n from each + # hypothesis's weighted χ². Falls back to unit σ when phase_error is absent so older + # callers still get sane numbers. sigma = np.asarray( mode_result.phase_error if mode_result.phase_error is not None @@ -825,7 +826,10 @@ def _wrap180(x): post = np.exp(-(chi2 - chi2.min()) / 2.0) post /= post.sum() n_confidence = float(post[best_j]) - phase_sigma = float(np.sqrt(1.0 / np.sum(1.0 / sigma**2))) + # 1σ of the ACTUALLY-returned intercept — the w-weighted circular mean, not the + # inverse-variance-optimal one: Var(c) = Σ(w_k·σ_k)² / (Σ w_k)². (Σw > 0: w is + # replaced by ones above if every weight was zero.) + phase_sigma = float(np.sqrt(np.sum((w * sigma) ** 2)) / w.sum()) return ToroidalFitResult( kind="toroidal_fit", diff --git a/tests/test_crucible_regressions.py b/tests/test_crucible_regressions.py new file mode 100644 index 0000000..94773fc --- /dev/null +++ b/tests/test_crucible_regressions.py @@ -0,0 +1,110 @@ +"""Regression tests for the 2026-07-01 Crucible review fixes. + +Each test names one confirmed defect and fails against the pre-fix behavior: + P3 SVD per-coefficient σ summed the wrong axis (wrong QS error bars) + P4 toroidal phase_sigma reported the BLUE variance, not the returned estimator's + P6 the 2D sinusoidal-integral basis was mis-normalized by -180/π for m≠0 modes + P7 reconstruct_grid used c instead of conj(c) (toroidal mirror); underdetermined + fits raised IndexError +""" + +from __future__ import annotations + +from types import SimpleNamespace + +import numpy as np +import pytest + +from magnetics._slcontour import fit as slfit +from magnetics.core import quasistationary +from magnetics.core.spectral import fit_toroidal_mode + + +# ── P6: 2D sinusoidal-integral basis normalization ────────────────────────── +def test_slcontour_2d_sinusoidal_integral_basis_is_separable_product(): + """The (n≠0, m≠0) integral basis must equal the product of the two 1D branches. + Pre-fix it divided by deg2rad(dx*dy)*n*m (π/180 once, no i²), off by -180/π.""" + x1 = np.array([10.0, 100.0]) + x2 = x1 + 20.0 + y1 = np.array([5.0, 60.0]) + y2 = y1 + 15.0 + n, m = 2, 3 + fmn_2d = slfit.form_basis_function(n, m, x1, x2, y1, y2, fit_basis="sinusoidal-integral") + fmn_n = slfit.form_basis_function(n, 0, x1, x2, y1, y2, fit_basis="sinusoidal-integral") + fmn_m = slfit.form_basis_function(0, m, x1, x2, y1, y2, fit_basis="sinusoidal-integral") + assert np.allclose(fmn_2d, fmn_n * fmn_m) + + +# ── P3: SVD per-coefficient uncertainty axis ──────────────────────────────── +def test_svd_covariance_diagonal_sums_the_singular_axis(): + """diag((AᵀA)⁻¹) = Σ_k V[j,k]² / w_k² sums over the singular index (axis=1 of + Vh.T). The old axis=0 collapsed to 1/w_k² by singular position — not the + covariance diagonal. Guards the reduction in fit.py:355 and quasistationary.py:173.""" + rng = np.random.default_rng(0) + a = rng.standard_normal((10, 4)) # non-orthogonal ⇒ condition number ≠ 1 + _, w, vh = np.linalg.svd(a, full_matrices=False) + w_inv = 1.0 / w + ref = np.diag(np.linalg.inv(a.T @ a)) + assert np.allclose(np.sum((vh.T * w_inv) ** 2, axis=1), ref) + assert not np.allclose(np.sum((vh.T * w_inv) ** 2, axis=0), ref) + + +# ── P7: underdetermined fit must not IndexError (resolved by the axis fix) ─── +def test_quasistationary_fit_underdetermined_sizes_sigmas_without_indexerror(): + phi = np.linspace(0.0, 300.0, 6) + theta = np.linspace(0.0, 150.0, 6) + signal = np.cos(np.deg2rad(phi))[:, None] + 0.3 * np.sin(np.deg2rad(theta))[:, None] + res = quasistationary.fit( + np.array([0.0]), + signal, + phi, + phi, + theta, + theta, + ns=(1, 2, 3), + ms=(0, 1, 2), # 9 modes → ~18 basis columns > 6 sensors + fit_basis="sinusoidal-point", + ) + assert res.sigmas.shape[0] == res.ns.shape[0] # one σ per mode; overran before the fix + assert np.all(np.isfinite(np.abs(res.sigmas))) + + +# ── P7: reconstruct_grid must not mirror a phase-shifted mode ──────────────── +def test_reconstruct_grid_roundtrips_phase_shifted_mode(): + """cos(φ − 40°) fit as n=1, m=0 then reconstructed must reproduce cos(φ − 40°). + Pre-fix (c instead of conj(c)) it produced cos(φ + 40°) — a toroidal mirror.""" + phi = np.linspace(0.0, 315.0, 8) + theta0 = np.zeros_like(phi) + signal = np.cos(np.deg2rad(phi - 40.0))[:, None] # [n_ch, 1] + res = quasistationary.fit( + np.array([0.0]), + signal, + phi, + phi, + theta0, + theta0, + ns=(1,), + ms=(0,), + fit_basis="sinusoidal-point", + ) + grid = np.linspace(0.0, 350.0, 36) + z = quasistationary.reconstruct_grid(res, grid, np.array([0.0]), 0) # [1, n_phi] + assert np.allclose(z[0], np.cos(np.deg2rad(grid - 40.0)), atol=1e-6) + + +# ── P4: toroidal phase_sigma reflects the actually-returned (weighted) estimator ─ +def test_toroidal_phase_sigma_matches_the_weighted_estimator_not_the_blue(): + mode = SimpleNamespace( + frequency=3000.0, + toroidal_angle=np.array([0.0, 120.0, 240.0]), + phase=np.array([0.0, 240.0, 120.0]), + amplitude=np.array([100.0, 1.0, 1.0]), # one loud, noisy probe dominates the weight + coherence=np.array([1.0, 1.0, 1.0]), + phase_error=np.array([30.0, 3.0, 3.0]), + ) + fit = fit_toroidal_mode(mode) # weights="amplitude" + w, sigma = mode.amplitude, mode.phase_error + expected = float(np.sqrt(np.sum((w * sigma) ** 2)) / w.sum()) + blue = float(np.sqrt(1.0 / np.sum(1.0 / sigma**2))) # the old, over-optimistic value + assert fit.phase_sigma == pytest.approx(expected) + assert fit.phase_sigma > 5 * blue diff --git a/tests/test_qs_bridge.py b/tests/test_qs_bridge.py index 2331958..1d1704f 100644 --- a/tests/test_qs_bridge.py +++ b/tests/test_qs_bridge.py @@ -45,6 +45,16 @@ def test_dropping_a_required_variable_raises_not_silently_misserves(fit_ds): qs_bridge.fit_to_qs_fit_node(broken) +def test_fit_quality_statuses_use_the_contract_vocabulary(fit_ds): + """The traffic-light status must be one of the GUI's Quality values + (good/warn/bad = contracts.quality_for_k), not the old 'ok'/'error' strings + which NodeView's QCOLOR can't color.""" + node = qs_bridge.fit_to_fit_quality_node(fit_ds) + statuses = [f["status"] for f in node["fields"] if "status" in f] + assert statuses # the K (raw)/K (eff) rows carry a status + assert all(s in {"good", "warn", "bad"} for s in statuses) + + def test_amplitude_sigma_is_finite(fit_ds): node = qs_bridge.fit_to_amplitude_node(fit_ds) sigma = np.asarray(node["meta"]["sigma"], dtype=float)