From 3750aacb40ec2def3ba9b0ff1d8088cf846c5e22 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 10 Jul 2026 00:35:19 -0400 Subject: [PATCH 01/16] spectral: wrap delta_phi and pick an alias-safe 2-point pair round(phase/delta_phi) assumed the unwrapped n*delta_phi stayed inside +/-180 deg, but the pair picker chose the widest raw separation (320 deg on fetched DIII-D shots), aliasing every rotating mode to n = 0. Wrap delta_phi to its signed principal value in the core estimators and pick the widest pair that still resolves |n| <= 5 unaliased. Co-Authored-By: Claude Fable 5 --- src/magnetics/core/spectral.py | 25 +++++++++++--- src/magnetics/service/nodes.py | 40 +++++++++++++++++----- tests/test_spectral.py | 61 ++++++++++++++++++++++++++++++++++ 3 files changed, 114 insertions(+), 12 deletions(-) diff --git a/src/magnetics/core/spectral.py b/src/magnetics/core/spectral.py index 04a28f6..ff1e7ff 100644 --- a/src/magnetics/core/spectral.py +++ b/src/magnetics/core/spectral.py @@ -165,6 +165,17 @@ def _cross_spectral_errors( return phase_err, amp_err +def wrap_angle_deg(angle: float) -> float: + """Reduce an angle (deg) to the signed principal value in (-180, 180]. + + A mode's cross-phase n·Δφ is only defined modulo 360°, so the 2-point estimator + n = round(phase / Δφ) needs the *short-way, signed* probe separation: probes at + 20° and 340° are 320° apart measured one way, but the wrapped phase sees -40°. + Feeding the long-way separation aliases every mode toward n = 0.""" + wrapped = (float(angle) + 180.0) % 360.0 - 180.0 + return 180.0 if wrapped == -180.0 else wrapped + + def cross_spectrum( sig1: NDArray[np.floating], sig2: NDArray[np.floating], @@ -180,7 +191,10 @@ def cross_spectrum( sig1 (ndarray): first probe time series. sig2 (ndarray): second probe time series. sample_rate (float): sampling rate (Hz). - delta_phi (float | None): toroidal separation (deg); enables mode-number output. + delta_phi (float | None): toroidal separation (deg), phi2 - phi1; enables + mode-number output. Wrapped internally to the principal value (-180, 180], + so either the short- or long-way separation may be passed. Modes are only + resolved unaliased for |n| <= 180 / |wrapped delta_phi|. nperseg (int | None): Welch segment length for csd/coherence; None uses scipy's default (256). Set below the signal length on short slices to keep multiple averaging segments — coherence is meaningless from a single segment. @@ -209,8 +223,9 @@ def cross_spectrum( n_segments=k, ) + delta_phi = wrap_angle_deg(delta_phi) if delta_phi == 0: - raise ValueError("delta_phi must be non-zero to compute mode numbers") + raise ValueError("delta_phi must be non-zero (mod 360°) to compute mode numbers") mode = np.rint(phase / delta_phi).astype(np.intp) @@ -331,7 +346,8 @@ def compute_spectrogram( time (ndarray): sample times (s); assumed uniformly sampled. sig1 (ndarray): first probe time series. sig2 (ndarray): second probe time series. - delta_phi (float): toroidal separation (deg); must be non-zero. + delta_phi (float): toroidal separation (deg), phi2 - phi1; must be non-zero + mod 360°. Wrapped internally to (-180, 180] like ``cross_spectrum``. slice_duration (float): FFT window width (s) — sets the frequency resolution. window (str): scipy window name for the taper. max_columns (int): cap on spectrogram time columns (decimation lever). @@ -340,8 +356,9 @@ def compute_spectrogram( result (SpectrogramResult): power/coherence/mode_number as (n_times, n_freqs) arrays, with time/frequency/rms_by_mode/mode_indices. """ + delta_phi = wrap_angle_deg(delta_phi) if delta_phi == 0: - raise ValueError("delta_phi must be non-zero to compute mode numbers") + raise ValueError("delta_phi must be non-zero (mod 360°) to compute mode numbers") time = np.asarray(time) s1 = np.ascontiguousarray(sig1, dtype=np.float32) diff --git a/src/magnetics/service/nodes.py b/src/magnetics/service/nodes.py index e27028e..8cd011c 100644 --- a/src/magnetics/service/nodes.py +++ b/src/magnetics/service/nodes.py @@ -339,21 +339,45 @@ def _named_sets(dg, substr: str) -> list[str]: # ── spectrogram: real 2-point MODESPEC cross-spectrogram ───────────────────── +_PAIR_N_MAX = 5 # the pair must resolve |n| <= this without aliasing (matches the n-map) + + +def _best_pair(arr) -> tuple[tuple[str, float], tuple[str, float]] | None: + """Best 2-point probe pair from a (name, phi)-sorted array: the widest *wrapped* + separation that still resolves |n| <= _PAIR_N_MAX unaliased. round(phase/Δφ) is + exact for all |n| <= N iff (N + 0.5)·|Δφ| <= 180, and within that limit a wider + Δφ dilutes phase noise (σ_n = σ_φ/Δφ) — so: widest under the limit; if the array + has no such pair, the smallest non-zero separation (best available n range).""" + limit = 180.0 / (_PAIR_N_MAX + 0.5) + best = smallest = None + for i in range(len(arr)): + for j in range(i + 1, len(arr)): + sep = abs(spectral.wrap_angle_deg(arr[j][1] - arr[i][1])) + if sep == 0: + continue + if sep <= limit and (best is None or sep > best[0]): + best = (sep, arr[i], arr[j]) + if smallest is None or sep < smallest[0]: + smallest = (sep, arr[i], arr[j]) + pick = best or smallest + return None if pick is None else (pick[1], pick[2]) + + def _pick_pair(shot) -> tuple[tuple[str, float], tuple[str, float]]: """Two toroidally-separated probes for the 2-point cross-spectrogram. Prefer the fast Mirnov dB/dt array (DIII-D MPI_BDOT / an NSTX toroidal set / a KSTAR declared toroidal array), then integrated Bp. Returns ((name1, phi1), - (name2, phi2)) with the widest non-zero separation.""" + (name2, phi2)) chosen by ``_best_pair`` (aliasing-safe, noise-optimal).""" dev, arrays = _arrays(shot) if arrays: # device declares an explicit toroidal set (KSTAR) → use it precisely - arr = _set_channels(dev, arrays["toroidal"], shot) - if len(arr) >= 2 and arr[0][1] != arr[-1][1]: - return arr[0], arr[-1] + pair = _best_pair(_set_channels(dev, arrays["toroidal"], shot)) + if pair is not None: + return pair raise ValueError("need two toroidally-separated probes for a spectrogram") for families in _pick_pair_prefs(shot): - arr = _array_channels(shot, families) # (name, phi), sorted by phi - if len(arr) >= 2 and arr[0][1] != arr[-1][1]: - return arr[0], arr[-1] + pair = _best_pair(_array_channels(shot, families)) # (name, phi), sorted by phi + if pair is not None: + return pair raise ValueError("need two toroidally-separated probes for a spectrogram") @@ -394,7 +418,7 @@ def _spec_result(shot: str, slice_duration: float, coherence_smooth: int, max_co coherence_smooth=coherence_smooth, max_columns=max_columns, ) - return res, (n1, n2), round(float(phi2 - phi1), 1) + return res, (n1, n2), round(spectral.wrap_angle_deg(phi2 - phi1), 1) def _prep_spec(shot, params): diff --git a/tests/test_spectral.py b/tests/test_spectral.py index a2b810b..aced544 100644 --- a/tests/test_spectral.py +++ b/tests/test_spectral.py @@ -688,3 +688,64 @@ def test_smooth_mode_spectrogram(self): assert sm.quality.shape == ms.quality.shape assert np.all((sm.quality >= 0.0) & (sm.quality <= 1.0 + 1e-6)) np.testing.assert_array_equal(sm.mode_number, ms.mode_number) + + +# ----------------------------------------------------------------------- +# wide-pair aliasing regression (probes 320° apart the long way = -40° wrapped) +# ----------------------------------------------------------------------- + + +class TestWidePairAliasing: + """Regression for the production DIII-D pair MPI66M020D/340D (φ = 19.52°/339.7°): + Δφ = φ2 − φ1 = 320.2° the long way, −39.8° wrapped. Unwrapped, round(phase/Δφ) + reported n = 0 for every real mode. n is SIGNED here — sig ~ sin(ωt − nφ) with + Δφ = φ2 − φ1 (the service convention) must recover +n.""" + + PHI1, PHI2 = 19.52, 339.7 # real MPI66M020D / MPI66M340D toroidal angles + + @staticmethod + def _pair(n, phi1, phi2, fs=50_000, f_mode=3_000.0, duration=0.1): + t = np.linspace(0, duration, int(fs * duration), endpoint=False) + w = 2 * np.pi * f_mode + sig1 = np.sin(w * t - np.deg2rad(n * phi1)) + sig2 = np.sin(w * t - np.deg2rad(n * phi2)) + return t, sig1, sig2, float(fs) + + @pytest.mark.parametrize("n_true", [1, 2, 3, 4, -2]) + def test_cross_spectrum_recovers_signed_n(self, n_true): + _, sig1, sig2, fs = self._pair(n_true, self.PHI1, self.PHI2) + result = cross_spectrum(sig1, sig2, fs, delta_phi=self.PHI2 - self.PHI1) + peak = np.argmax(result.power) + assert result.mode_number[peak] == n_true + + @pytest.mark.parametrize("n_true", [1, 3]) + def test_spectrogram_recovers_signed_n(self, n_true): + t, sig1, sig2, fs = self._pair(n_true, self.PHI1, self.PHI2) + res = compute_spectrogram(t, sig1, sig2, delta_phi=self.PHI2 - self.PHI1) + it, fi = np.unravel_index(np.argmax(res.power), res.power.shape) + assert res.mode_number[it, fi] == n_true + + def test_long_and_short_way_separations_agree(self): + _, sig1, sig2, fs = self._pair(2, self.PHI1, self.PHI2) + long_way = cross_spectrum(sig1, sig2, fs, delta_phi=self.PHI2 - self.PHI1) + short_way = cross_spectrum(sig1, sig2, fs, delta_phi=self.PHI2 - self.PHI1 - 360.0) + np.testing.assert_array_equal(long_way.mode_number, short_way.mode_number) + np.testing.assert_array_equal(long_way.mode_indices, short_way.mode_indices) + + def test_full_turn_separation_raises(self): + _, sig1, sig2, fs = self._pair(1, 0.0, 360.0) + with pytest.raises(ValueError, match="mod 360"): + cross_spectrum(sig1, sig2, fs, delta_phi=360.0) + + +class TestWrapAngleDeg: + def test_values(self): + from magnetics.core.spectral import wrap_angle_deg + + assert wrap_angle_deg(320.18) == pytest.approx(-39.82) + assert wrap_angle_deg(-200.0) == pytest.approx(160.0) + assert wrap_angle_deg(33.0) == pytest.approx(33.0) + assert wrap_angle_deg(180.0) == 180.0 + assert wrap_angle_deg(-180.0) == 180.0 + assert wrap_angle_deg(540.0) == 180.0 + assert wrap_angle_deg(720.0) == 0.0 From 140a28dee689b69be060b0ac0226e2b4baa28e80 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 10 Jul 2026 00:50:27 -0400 Subject: [PATCH 02/16] spectral: pin one signed-n convention across all estimators A cos(n*phi - omega*t) mode measured through conj(sig)*ref cross-phases ramps as +n*phi, but fit_toroidal_mode fitted c - n*phi, array_mode_spectrogram projected onto an unconjugated template, and mac_n_spectrum matched e^{-in*phi} shapes - so the phase-fit family reported -n while the 2-point family reported +n, despite docstrings claiming they agree. Flip the fit family to +n, fix the service ramp and poloidal detrend, use the production delta_phi = phi2 - phi1 convention in fixtures, and pin all five estimators to the same signed n with a cross-family test. Co-Authored-By: Claude Fable 5 --- src/magnetics/core/mode_shape.py | 10 ++-- src/magnetics/core/spectral.py | 31 +++++++----- src/magnetics/service/nodes.py | 15 +++--- tests/conftest.py | 6 ++- tests/test_mode_shape.py | 7 +-- tests/test_spectral.py | 82 ++++++++++++++++++++++++++++---- 6 files changed, 115 insertions(+), 36 deletions(-) diff --git a/src/magnetics/core/mode_shape.py b/src/magnetics/core/mode_shape.py index 2a8d70f..8ddcadb 100644 --- a/src/magnetics/core/mode_shape.py +++ b/src/magnetics/core/mode_shape.py @@ -297,16 +297,20 @@ def mac_n_spectrum( *, n_range: tuple[int, int] = (-6, 6), ) -> tuple[NDArray[np.integer], NDArray[np.floating], int]: - """MAC of a measured shape vector against ideal pure-mode templates e^{−inφ}. + """MAC of a measured shape vector against ideal pure-mode templates e^{+inφ}. Gives a *shape-based*, geometry-aware toroidal mode-number identification (how closely the measured array pattern resembles each pure rotating mode), a useful - cross-check on the cross-phase fit. Returns (n_values, mac_values, best_n). + cross-check on the cross-phase fit. The measured shape comes from the + conj(sig)·ref cross-phases (``shape_vector`` of ``extract_mode_at_frequency`` / + ``mode_from_spectrum`` output), where a ``cos(nφ - ωt)`` mode appears as a + ``+nφ`` phase ramp, i.e. ``z ∝ e^{+inφ}`` — so matching that template reports + the same signed n as the cross-phase fit. Returns (n_values, mac_values, best_n). """ phi = np.deg2rad(np.asarray(angle_deg, dtype=np.float64)) z = np.asarray(complex_shape, dtype=np.complex128) ns = np.arange(n_range[0], n_range[1] + 1) - macs = np.array([mac(z, np.exp(-1j * n * phi)) for n in ns]) + macs = np.array([mac(z, np.exp(1j * n * phi)) for n in ns]) return ns, macs, int(ns[int(np.argmax(macs))]) diff --git a/src/magnetics/core/spectral.py b/src/magnetics/core/spectral.py index ff1e7ff..679070b 100644 --- a/src/magnetics/core/spectral.py +++ b/src/magnetics/core/spectral.py @@ -622,7 +622,7 @@ class ArrayModeSpectrogram: time: NDArray[np.floating] # (n_times,) window-center times (s) freq_band: NDArray[np.floating] # (n_band,) Hz kept mode_number: NDArray[np.integer] # (n_times, n_band) best-fit toroidal n - amplitude: NDArray[np.floating] # (n_times, n_band) |Σ_p Z_p e^{-inφ}| + amplitude: NDArray[np.floating] # (n_times, n_band) |⟨Z, e^{-inφ}⟩| = |Σ_p Z_p e^{+inφ}| quality: NDArray[np.floating] # (n_times, n_band) harmonic energy fraction ∈ [1/M,1] @@ -642,11 +642,13 @@ def array_mode_spectrogram( The per-cell mode-coherence is the *energy fraction* the single best-fit n captures of the summed toroidal-harmonic power, ``|R_{n*}|² / Σ_n |R_n|² ∈ [1/M, 1]`` with - ``R_n = Σ_p Z_p e^{-inφ_p}`` and ``M = 2·n_max + 1`` candidates: 1 = a pure single-n + ``R_n = Σ_p Z_p e^{+inφ_p}`` and ``M = 2·n_max + 1`` candidates: 1 = a pure single-n pattern, 1/M = white across harmonics (incoherent noise). This is a spectral- concentration ratio (more noise/signal contrast than the resultant length ``|R_{n*}| / - Σ_p|Z_p|``, whose noise floor sits higher at ~1/√P). The ``exp(-i n φ)`` sign matches - ``mode_from_spectrum`` and the phase fit, so the reported n agrees with them. + Σ_p|Z_p|``, whose noise floor sits higher at ~1/√P). ``R_n`` is the inner product + ``⟨Z, e^{-inφ}⟩`` (template conjugated): a ``cos(nφ - ωt)`` mode's positive-frequency + STFT pattern is ``Z_p ∝ e^{-inφ_p}``, so the projection peaks at the same signed +n + that ``cross_spectrum`` and the phase fit report. Inputs: spectrum (ArrayShapeSpectrum): the per-probe complex STFT band. @@ -658,8 +660,8 @@ def array_mode_spectrogram( z = np.asarray(spectrum.spec) # (P, T, F) complex phi = np.deg2rad(np.asarray(angle_deg, dtype=np.float64)) ns = np.arange(-int(n_max), int(n_max) + 1) - basis = np.exp(-1j * ns[:, None] * phi[None, :]) # (M, P) - proj = np.einsum("mp,ptf->mtf", basis, z) # (M, T, F) = R_n + basis = np.exp(1j * ns[:, None] * phi[None, :]) # (M, P) = conj(e^{-inφ_p}) templates + proj = np.einsum("mp,ptf->mtf", basis, z) # (M, T, F) = R_n = ⟨Z, e^{-inφ}⟩ amp = np.abs(proj) k = np.argmax(amp, axis=0) # (T, F) peak = np.take_along_axis(amp, k[None], axis=0)[0] # |R_{n*}| (T, F) @@ -884,11 +886,14 @@ def fit_toroidal_mode( """Fit a toroidal mode number to per-probe phase-vs-angle data. A rotating mode of toroidal number n imprints a linear phase ramp across the - toroidal array: ``phase(phi) = c - n * phi`` (deg). Rather than unwrap the - (mod-360) phases — ill-posed for sparse arrays — this scans integer candidates - and picks the n whose residual ``phase + n * phi`` clusters most tightly on the - circle (largest amplitude-weighted resultant). The intercept c is the angle of - that resultant, giving a wrap-free fit line for the GUI to draw. + toroidal array. With the conj(sig)·ref cross-phase convention used by + ``extract_mode_at_frequency`` / ``mode_from_spectrum``, a ``cos(nφ - ωt)`` mode + measures as ``phase(phi) = c + n * phi`` (deg) — the same signed n the 2-point + ``cross_spectrum`` reports. Rather than unwrap the (mod-360) phases — ill-posed + for sparse arrays — this scans integer candidates and picks the n whose residual + ``phase - n * phi`` clusters most tightly on the circle (largest + amplitude-weighted resultant). The intercept c is the angle of that resultant, + giving a wrap-free fit line for the GUI to draw. Inputs: mode_result (ModeAtFrequencyResult): per-probe phase/amplitude/angle, e.g. @@ -919,7 +924,7 @@ def fit_toroidal_mode( inter = np.empty(len(candidates)) # per-candidate intercept c_n (deg) rmag = np.empty(len(candidates)) # per-candidate resultant length for j, n in enumerate(candidates): - resultant = np.sum(w * np.exp(1j * np.deg2rad(phase + n * phi))) / w.sum() + resultant = np.sum(w * np.exp(1j * np.deg2rad(phase - n * phi))) / w.sum() rmag[j] = np.abs(resultant) inter[j] = np.rad2deg(np.angle(resultant)) best_j = int(np.argmax(rmag)) @@ -944,7 +949,7 @@ def _wrap180(x): chi2 = np.array( [ - np.sum((_wrap180(phase + n * phi - inter[j]) / sigma) ** 2) + np.sum((_wrap180(phase - n * phi - inter[j]) / sigma) ** 2) for j, n in enumerate(candidates) ] ) diff --git a/src/magnetics/service/nodes.py b/src/magnetics/service/nodes.py index 8cd011c..f7f11ae 100644 --- a/src/magnetics/service/nodes.py +++ b/src/magnetics/service/nodes.py @@ -884,7 +884,7 @@ def _array_mode_spec(shot, names, phis, slice_duration): # ── toroidal mode at one frequency/cursor (shared by phase_fit & mode_shape) ── def _toroidal_arr(shot): """The toroidal (midplane) array for the n-fit: ONE consistent probe type, all at - θ≈0 so the phase is a clean −nφ ramp. Prefer the fast-Mirnov dB/dt array; a "both" + θ≈0 so the phase is a clean +nφ ramp. Prefer the fast-Mirnov dB/dt array; a "both" pull also brings the integrated-Bp (MPID) family and the off-midplane *poloidal* probes — mixing those in (different units, 90° dB/dt-vs-B offset, m·θ dependence) scrambles the fit, so they're excluded.""" @@ -933,12 +933,13 @@ def _toroidal_mode(shot, params): # ── phase_fit: phase-vs-φ at one frequency, at the GUI time cursor ──────────── def _wrapped_ramp(intercept_deg, slope_n) -> dict: - """Fitted line phase(a) = (c − n·a) mod 360 over a∈[0,360], WRAPPED so it traces - the same |n| sawteeth as the (wrapped) data instead of one line shooting off-axis. - A null break is inserted at each 0/360 wrap so the polyline doesn't draw a vertical + """Fitted line phase(a) = (c + n·a) mod 360 over a∈[0,360] (the conj(sig)·ref + cross-phase ramp ``fit_toroidal_mode`` fits), WRAPPED so it traces the same |n| + sawteeth as the (wrapped) data instead of one line shooting off-axis. A null + break is inserted at each 0/360 wrap so the polyline doesn't draw a vertical jump across the panel.""" a = np.linspace(0.0, 360.0, 361) - y = (intercept_deg - slope_n * a) % 360.0 + y = (intercept_deg + slope_n * a) % 360.0 fx: list = [] fy: list = [] prev = None @@ -1105,7 +1106,7 @@ def _toroidal_n(shot, t0_s, f_khz): def _poloidal_mode(shot, params): """Per-probe phase/amplitude across the poloidal array vs θ. The probes span φ, so - the toroidal nφ ramp is removed (phase += n·φ → −m·θ + const) using the toroidal n + the toroidal +nφ ramp is removed (phase −= n·φ → m·θ + const) using the toroidal n at the same (t0, f). The probe angles are mapped to the elongation-corrected θ* (using the EFIT κ at the cursor) when available, so an `m` mode is a clean sinusoid rather than a κ-distorted one. Returns (arr, mode, f_khz, t0_ms, kappa).""" @@ -1121,7 +1122,7 @@ def _poloidal_mode(shot, params): spec = _array_spectrum(str(shot), tuple(n for n, _ in arr)) t0_s = (t0_ms * 1e-3) if t0_ms is not None else float(spec.time[spec.time.size // 2]) mode = spectral.mode_from_spectrum(spec, thetas, t0_s, f_khz * 1e3) - detrended = (mode.phase + _toroidal_n(str(shot), t0_s, f_khz) * pphis) % 360.0 + detrended = (mode.phase - _toroidal_n(str(shot), t0_s, f_khz) * pphis) % 360.0 mode = dataclasses.replace(mode, phase=detrended) return arr, mode, f_khz, t0_ms, kappa diff --git a/tests/conftest.py b/tests/conftest.py index a945cd6..55607a5 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -127,14 +127,16 @@ def synthetic_n2(): t = np.linspace(0, duration, int(fs * duration), endpoint=False) f_mode = 3_000.0 n = 2 - phi1, phi2 = 30.0, 63.0 # Δφ = -33°, mimics the real 307/340 pair separation + phi1, phi2 = 30.0, 63.0 # Δφ = +33°, mimics the real 307/340 pair separation sig1 = np.sin(2 * np.pi * f_mode * t - np.deg2rad(n * phi1)) sig2 = np.sin(2 * np.pi * f_mode * t - np.deg2rad(n * phi2)) return { "time": t, "sig1": sig1, "sig2": sig2, - "delta_phi": phi1 - phi2, + # φ2 − φ1: the service convention (nodes._spec_result), so estimators + # recover the SIGNED n of a cos(nφ − ωt) mode as +n. + "delta_phi": phi2 - phi1, "fs": float(fs), "f_mode": f_mode, "n_true": n, diff --git a/tests/test_mode_shape.py b/tests/test_mode_shape.py index 9d73086..a671476 100644 --- a/tests/test_mode_shape.py +++ b/tests/test_mode_shape.py @@ -152,9 +152,10 @@ class TestMacNSpectrum: def test_peaks_at_true_mode(self): phi = np.array([0.0, 33.0, 66.0, 120.0, 200.0, 300.0]) n_true = 3 - z = shape_vector(-n_true * phi, np.ones_like(phi)) # phase = -n·φ + # phase = +n·φ: the conj(sig)·ref cross-phase ramp a cos(nφ − ωt) mode measures + z = shape_vector(n_true * phi, np.ones_like(phi)) ns, macs, best = mac_n_spectrum(phi, z) - assert abs(best) == n_true + assert best == n_true # signed assert macs.max() > 0.99 assert ns.shape == macs.shape @@ -259,7 +260,7 @@ def test_recovers_mode_number_during_mode(self): sigs, phi, t = self._array() tr = track_mode_shape(sigs, phi, t, frequency=8000.0, n_slices=20) - assert abs(int(np.median(tr.n_by_time[tr.t_ms < 25.0]))) == 2 + assert int(np.median(tr.n_by_time[tr.t_ms < 25.0])) == 2 # signed class TestGPScaleInvariance: diff --git a/tests/test_spectral.py b/tests/test_spectral.py index aced544..a10adfd 100644 --- a/tests/test_spectral.py +++ b/tests/test_spectral.py @@ -101,7 +101,7 @@ def test_recovers_synthetic_mode_number(self, synthetic_n2): d = synthetic_n2 result = cross_spectrum(d["sig1"], d["sig2"], d["fs"], delta_phi=d["delta_phi"]) peak_idx = np.argmax(result.power) - assert abs(result.mode_number[peak_idx]) == d["n_true"] + assert result.mode_number[peak_idx] == d["n_true"] # signed def test_peak_frequency_near_expected(self, synthetic_n2): d = synthetic_n2 @@ -113,7 +113,7 @@ def test_rms_by_mode_peaks_at_true_n(self, synthetic_n2): d = synthetic_n2 result = cross_spectrum(d["sig1"], d["sig2"], d["fs"], delta_phi=d["delta_phi"]) peak_mode_idx = np.argmax(result.rms_by_mode) - assert abs(result.mode_indices[peak_mode_idx]) == d["n_true"] + assert result.mode_indices[peak_mode_idx] == d["n_true"] # signed def test_coherence_high_for_correlated_signals(self, synthetic_n2): d = synthetic_n2 @@ -168,11 +168,11 @@ def test_recovers_mode_across_all_slices(self, synthetic_n2): ) for i in range(len(result.time)): peak_idx = np.argmax(result.power[i]) - assert abs(result.mode_number[i, peak_idx]) == d["n_true"] + assert result.mode_number[i, peak_idx] == d["n_true"] # signed def test_on_real_data(self, shot_174446): d = shot_174446 - delta_phi = d["phi_307"] - d["phi_340"] + delta_phi = d["phi_340"] - d["phi_307"] # φ2 − φ1, service convention result = compute_spectrogram( d["time_s"], d["sig_307"], @@ -290,7 +290,7 @@ def test_preserves_dominant_mode_peak(self, synthetic_n2): dn = denoise_spectrogram(spec, coherence_min=0.5, power_floor_k=None) peak = np.unravel_index(np.argmax(spec.power), spec.power.shape) assert dn.power[peak] > 0 - assert abs(dn.mode_number[peak]) == d["n_true"] + assert dn.mode_number[peak] == d["n_true"] # signed def test_power_floor_removes_more_than_coherence_alone(self, synthetic_n2): spec = self._spec(synthetic_n2) @@ -300,7 +300,7 @@ def test_power_floor_removes_more_than_coherence_alone(self, synthetic_n2): def test_on_real_data(self, shot_174446): d = shot_174446 - delta_phi = d["phi_307"] - d["phi_340"] + delta_phi = d["phi_340"] - d["phi_307"] # φ2 − φ1, service convention spec = compute_spectrogram( d["time_s"], d["sig_307"], d["sig_340"], delta_phi, slice_duration=0.004 ) @@ -513,7 +513,7 @@ def test_reference_probe_flagged_nan(self): def test_fit_reports_confidence_and_sigma(self): mode, n_true = self._modes() fit = fit_toroidal_mode(mode) - assert abs(fit.n) == n_true + assert fit.n == n_true # signed assert fit.phase_sigma is not None and fit.phase_sigma > 0 assert fit.n_confidence is not None assert 0.0 < fit.n_confidence <= 1.0 @@ -656,7 +656,7 @@ def test_recovers_high_mode_number_at_peak(self, n_true): fi = int(np.argmin(np.abs(ms.freq_band - 8000.0))) # the injected mode bin it = ms.time.size // 2 n_fit = fit_toroidal_mode(mode_from_spectrum(spec, phi, t[it], 8000.0)).n - assert abs(ms.mode_number[it, fi]) == n_true # right |n| (not aliased) + assert ms.mode_number[it, fi] == n_true # right SIGNED n (not aliased) assert ms.mode_number[it, fi] == n_fit # agrees with the phase fit assert ms.quality[it, fi] > 0.9 # clean single-n → high q @@ -749,3 +749,69 @@ def test_values(self): assert wrap_angle_deg(-180.0) == 180.0 assert wrap_angle_deg(540.0) == 180.0 assert wrap_angle_deg(720.0) == 0.0 + + +# ----------------------------------------------------------------------- +# cross-family signed-n agreement — every estimator, one physical mode +# ----------------------------------------------------------------------- + + +class TestSignedNConventionAgreement: + """One physical mode b ∝ cos(nφ − ωt) (n = +2, rotating toward +φ) must be + reported with the SAME signed n by every estimator family: the 2-point + cross-spectrum/spectrogram, the phase-ramp fit, the array projection, and the + shape MAC. Regression for the families disagreeing in sign (docstrings claimed + agreement while the fit family returned −n).""" + + N_TRUE = 2 + F_MODE = 8_000.0 + + @pytest.fixture() + def array_mode(self): + fs = 200_000 + t = np.linspace(0, 0.05, int(fs * 0.05), endpoint=False) + phi = np.linspace(0.0, 330.0, 12) + rng = np.random.default_rng(11) + sigs = np.vstack( + [ + np.sin(2 * np.pi * self.F_MODE * t - np.deg2rad(self.N_TRUE * p)) + + 0.02 * rng.standard_normal(t.size) + for p in phi + ] + ) + return sigs, phi, t, float(fs) + + def test_all_families_agree_on_signed_n(self, array_mode): + from magnetics.core.mode_shape import mac_n_spectrum, shape_vector + from magnetics.core.spectral import ( + array_mode_spectrogram, + array_shape_spectrum, + ) + + sigs, phi, t, fs = array_mode + + # family A1: 2-point cross-spectrum (adjacent probes, Δφ = φ2 − φ1 = 30°) + two_pt = cross_spectrum(sigs[0], sigs[1], fs, delta_phi=phi[1] - phi[0]) + n_2pt = int(two_pt.mode_number[np.argmax(two_pt.power)]) + + # family A2: 2-point spectrogram + spec2 = compute_spectrogram(t, sigs[0], sigs[1], delta_phi=phi[1] - phi[0]) + it, fi = np.unravel_index(np.argmax(spec2.power), spec2.power.shape) + n_specgram = int(spec2.mode_number[it, fi]) + + # family B1: phase-ramp fit on the full array + mode = extract_mode_at_frequency(sigs, phi, t, frequency=self.F_MODE) + n_fit = fit_toroidal_mode(mode).n + + # family B2: array harmonic projection + aspec = array_shape_spectrum(sigs, t) + ms = array_mode_spectrogram(aspec, phi) + fb = int(np.argmin(np.abs(ms.freq_band - self.F_MODE))) + n_array = int(ms.mode_number[ms.time.size // 2, fb]) + + # family B3: shape MAC + _, _, n_mac = mac_n_spectrum(phi, shape_vector(mode.phase, mode.amplitude)) + + assert n_2pt == n_specgram == n_fit == n_array == n_mac == self.N_TRUE, ( + f"2pt={n_2pt} specgram={n_specgram} fit={n_fit} array={n_array} mac={n_mac}" + ) From 04613ad01cb89154801894654fce106937f56a49 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 10 Jul 2026 00:54:40 -0400 Subject: [PATCH 03/16] fetch: merge remote pulls into an existing shot file instead of clobbering The remote backend dispatched before the incremental-merge bookkeeping and rsynced the cluster result directly over shot_.h5 - since the cluster file holds only the requested selection, a one-signal custom pull (the GUI's QS custom-signal panel on the default DIII-D remote backend) silently destroyed every previously fetched channel. Dispatch the remote branch after the merge/skip resolution, stage the cluster file beside the local one when a compatible file exists, and fold it in through _write_h5's merge path (time-base dedup, replace-on-conflict, fetched/missing bookkeeping identical to a local incremental fetch). Co-Authored-By: Claude Fable 5 --- src/magnetics/data/fetch/toksearch.py | 121 ++++++++++++++++++++------ tests/test_toksearch_writer.py | 109 +++++++++++++++++++++++ 2 files changed, 205 insertions(+), 25 deletions(-) diff --git a/src/magnetics/data/fetch/toksearch.py b/src/magnetics/data/fetch/toksearch.py index 54d6ab4..cfbc697 100644 --- a/src/magnetics/data/fetch/toksearch.py +++ b/src/magnetics/data/fetch/toksearch.py @@ -1086,6 +1086,60 @@ def _write_h5( return got, missing +def _merge_h5_files(src, dst, *, compression="lzf"): + """Fold every channel of shot file ``src`` into the existing shot file ``dst`` + (same shot/window/decimation — the caller verifies via ``_existing_channels``). + + Channels are re-written through ``_write_h5``'s merge path so time-base dedup, + replace-on-conflict, and the fetched/missing bookkeeping behave exactly like a + local incremental fetch. Used by the remote backend: the cluster always ships a + fresh file holding only the requested selection, which must never replace a + local file holding earlier pulls. + """ + import h5py + + channels: list[Channel] = [] + query_names: dict[str, str] = {} + geometry: dict[str, dict] = {} + with h5py.File(src, "r") as h5: + analysis = h5.attrs.get("analysis", "") + analysis = analysis.decode() if isinstance(analysis, bytes) else str(analysis) + tmin = h5.attrs.get("tmin") + tmax = h5.attrs.get("tmax") + stride = int(h5.attrs.get("decimate", 1)) + tmin = None if isinstance(tmin, (bytes, str)) else tmin + tmax = None if isinstance(tmax, (bytes, str)) else tmax + for name, g in h5.items(): + if name == "_timebases" or not isinstance(g, h5py.Group) or "data" not in g: + continue + channels.append(Channel(name, time=g["time"][()], data=g["data"][()], ok=True)) + attrs = dict(g.attrs) + pn = attrs.pop("pointname", None) + if pn is not None: + query_names[name] = pn.decode() if isinstance(pn, bytes) else str(pn) + geometry[name] = { + k: v + for k, v in attrs.items() + if k != "time_units" and isinstance(v, (int, float, np.integer, np.floating)) + } + for m in sorted(_attr_names(h5, "channels_missing")): + channels.append(Channel(m, ok=False, error="missing in remote pull")) + return _write_h5( + dst, + None, # shot/device attrs already present in dst (merge never rewrites them) + analysis, + "remote", + channels, + compression=compression, + tmin=tmin, + tmax=tmax, + stride=stride, + query_names=query_names, + channel_geometry=geometry, + merge=True, + ) + + def _existing_channels(path, tmin, tmax, stride): """Names already attempted in an existing shot file IF its reduction matches. @@ -1203,31 +1257,10 @@ def fetch_shot( progress(0.0, f"{backend}→mdsthin (tree device)") backend = "mdsthin" - if backend == "remote" and not _tree_transport: - # Orchestrate a pull on the cluster from here; remote side runs this same - # script with --backend toksearch and writes the file we copy back. - from . import remote as remote_run - - kw = {} if remote_python is None else {"python": remote_python} - return remote_run.run_remote( - shot, - analysis, - host=remote_host, - jump=ssh_jump, - username=username, - password=password, - duo=duo, - remote_dir=remote_dir, - tmin=tmin, - tmax=tmax, - decimate=decimate, - device=device, - sensor_set=sensor_set, - raw_pointnames=raw_pointnames, - local_out_dir=(str(Path(out).parent) if out else None), - progress=progress, - **kw, - ) + # NOTE: the remote backend is dispatched BELOW, after the output path and the + # incremental-merge bookkeeping are resolved — the cluster ships a fresh file + # holding only the requested selection, so it must be merged into (never copied + # over) an existing local shot file. # Device config is the source of truth for mdsip addresses; an explicit # gateway/server (CLI or caller) overrides the device file. (`dev`, `device_name`, @@ -1367,6 +1400,44 @@ def fetch_shot( sys.stderr.write(f"All {n_skipped} requested signals already in {out}; nothing to fetch.\n") return out + if backend == "remote" and not _tree_transport: + # Orchestrate a pull on the cluster from here; remote side runs this same + # script with --backend toksearch and writes the file we copy back. When a + # compatible local file already exists (merge), the cluster result is staged + # beside it and folded in channel-by-channel — copying it over the local file + # would silently destroy every previously pulled channel set. + from . import remote as remote_run + + kw = {} if remote_python is None else {"python": remote_python} + stage_dir = out_path.parent / "_remote_stage" if merge else out_path.parent + fetched = remote_run.run_remote( + shot, + analysis, + host=remote_host, + jump=ssh_jump, + username=username, + password=password, + duo=duo, + remote_dir=remote_dir, + tmin=tmin, + tmax=tmax, + decimate=decimate, + device=device, + sensor_set=sensor_set, + raw_pointnames=raw_pointnames, + local_out_dir=str(stage_dir), + progress=progress, + **kw, + ) + if merge: + progress(0.95, f"merging into {out_path.name}") + _merge_h5_files(fetched, out) + Path(fetched).unlink(missing_ok=True) + with contextlib.suppress(OSError): + stage_dir.rmdir() + return out + return fetched + t0 = time.perf_counter() if _tree_transport: # Tree device reached via its own transport (KSTAR: KFE VPN + nkstar tunnel). diff --git a/tests/test_toksearch_writer.py b/tests/test_toksearch_writer.py index 461948b..b716d55 100644 --- a/tests/test_toksearch_writer.py +++ b/tests/test_toksearch_writer.py @@ -182,3 +182,112 @@ def compute_serial(self): assert by_name["PT"].ok is True assert by_name["kappa"].ok is False assert by_name["kappa"].error == "no tree server configured" + + +# --------------------------------------------------------------------------- +# _merge_h5_files — the remote-backend result must fold into an existing file +# --------------------------------------------------------------------------- + + +def test_merge_h5_files_preserves_existing_channels(tmp_path): + """Regression: a remote pull used to rsync the cluster file OVER the local shot + file, so a one-signal custom pull (e.g. 'betan') destroyed every previously + fetched channel. The staged file must merge in, replacing on name conflicts and + keeping everything else.""" + import h5py + + from magnetics.data.fetch.toksearch import _merge_h5_files + + t = np.linspace(0.0, 10.0, 100) + a = Channel("MPI66M307D", t, np.ones(100, np.float32), ok=True) + b = Channel("MPI66M340D", t, np.full(100, 2.0, np.float32), ok=True) + dst, _ = _write(tmp_path, [a, b]) + + stage = tmp_path / "stage" + stage.mkdir() + src = str(stage / "shot.h5") + t2 = np.linspace(0.0, 10.0, 50) + new = Channel("betan", t2, np.full(50, 1.5, np.float32), ok=True) + refetch = Channel("MPI66M340D", t, np.full(100, 9.0, np.float32), ok=True) + gone = Channel("ip", ok=False, error="no data") + _write_h5( + src, + 123, + "custom", + "toksearch", + [new, refetch, gone], + compression="lzf", + tmin=None, + tmax=None, + stride=1, + ) + + _merge_h5_files(src, dst) + + with h5py.File(dst, "r") as h5: + # pre-existing channel untouched, re-fetched channel replaced, new one added + assert set(h5) >= {"MPI66M307D", "MPI66M340D", "betan"} + np.testing.assert_array_equal(h5["MPI66M307D/data"][()], 1.0) + np.testing.assert_array_equal(h5["MPI66M340D/data"][()], 9.0) + np.testing.assert_array_equal(h5["betan/data"][()], 1.5) + np.testing.assert_allclose(h5["betan/time"][()], t2) + fetched = {x.decode() for x in h5.attrs["channels_fetched"]} + missing = {x.decode() for x in h5.attrs["channels_missing"]} + assert {"MPI66M307D", "MPI66M340D", "betan"} <= fetched + assert missing == {"ip"} + # both selection labels recorded + assert "custom" in str(h5.attrs["analysis"]) + + +def test_fetch_shot_remote_merges_instead_of_clobbering(tmp_path, monkeypatch): + """End-to-end dispatch regression: fetch_shot(backend='remote') with an existing + compatible shot file must stage the cluster result and merge it — the local file + keeps its earlier channels.""" + import h5py + + # existing local file with one 'earlier pull' channel, default window/stride + out = tmp_path / "shot_184927.h5" + t = np.linspace(0.0, 10.0, 100) + _write_h5( + str(out), + 184927, + "both", + "toksearch", + [Channel("MPI66M307D", t, np.ones(100, np.float32), ok=True)], + compression="lzf", + tmin=None, + tmax=None, + stride=1, + ) + + def fake_run_remote(shot, analysis="both", *, local_out_dir=None, progress=None, **kw): + # the cluster always writes a FRESH file with only the requested selection + path = str(f"{local_out_dir}/shot_{shot}.h5") + import os + + os.makedirs(local_out_dir, exist_ok=True) + t2 = np.linspace(0.0, 10.0, 50) + _write_h5( + path, + int(shot), + "custom", + "toksearch", + [Channel("betan", t2, np.full(50, 1.5, np.float32), ok=True)], + compression="lzf", + tmin=None, + tmax=None, + stride=1, + ) + return path + + from magnetics.data.fetch import remote as remote_mod + + monkeypatch.setattr(remote_mod, "run_remote", fake_run_remote) + + result = toksearch.fetch_shot(184927, backend="remote", raw_pointnames=["betan"], out=str(out)) + + assert result == str(out) + with h5py.File(out, "r") as h5: + assert "MPI66M307D" in h5, "remote pull clobbered the existing channel" + assert "betan" in h5 + assert not (tmp_path / "_remote_stage" / "shot_184927.h5").exists() From 2187d46dc8420f78b00c556b7a760fe585bc84a8 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 10 Jul 2026 00:58:18 -0400 Subject: [PATCH 04/16] qs_fit: de-normalize fit_signal by the same sigma that built the RHS The design matrix and RHS are normalized by the effective sigma (which honors sigma_override and NaN fills), but fit_signal was de-normalized by the dataset's raw signal_sigma - so any GUI edit of the "uncertainty (sigma)" box scaled fit_signal by signal_sigma/override, corrupting residual, chi_sq, and red_chi_sq while the coefficients stayed correct. Co-Authored-By: Claude Fable 5 --- src/magnetics/core/qs_fit.py | 15 ++++++++++----- tests/test_qs_bridge.py | 25 +++++++++++++++++++++++++ 2 files changed, 35 insertions(+), 5 deletions(-) diff --git a/src/magnetics/core/qs_fit.py b/src/magnetics/core/qs_fit.py index 14113cf..17d4c88 100644 --- a/src/magnetics/core/qs_fit.py +++ b/src/magnetics/core/qs_fit.py @@ -449,9 +449,15 @@ def _printv(*a): # ── assemble output Dataset ─────────────────────────────────────────────── fit_b = np.dot(A, fit_coeffs).real.reshape(ds["channel"].shape[0], -1) - ds["fit_signal"] = xr.DataArray(fit_b, coords=ds["signal"].coords, dims=ds["signal"].dims) * ds[ - "signal_sigma" - ].fillna(float(sigma.mean())) + # De-normalize by the SAME effective sigma that normalized the RHS `b` above + # (it honors sigma_override and NaN fills). Using the dataset's raw + # signal_sigma here would scale fit_signal by signal_sigma/override whenever + # they differ, corrupting residual and chi_sq while the coefficients stay + # correct. + sig_eff = xr.DataArray(sigma, coords={"channel": ds["channel"]}, dims=("channel",)) + ds["fit_signal"] = ( + xr.DataArray(fit_b, coords=ds["signal"].coords, dims=ds["signal"].dims) * sig_eff + ) ds["fit_ns"] = xr.DataArray(nms_arr[:, 0], coords={"mode": np.arange(len(nms))}, dims=("mode",)) ds["fit_ms"] = xr.DataArray(nms_arr[:, 1], coords={"mode": np.arange(len(nms))}, dims=("mode",)) ds["fit_coeffs"] = xr.DataArray( @@ -462,8 +468,7 @@ def _printv(*a): ) ds["residual"] = ds["signal"] - ds["fit_signal"] - sig_for_chi = xr.DataArray(sigma, coords={"channel": ds["channel"]}, dims=("channel",)) - ds["chi_sq"] = ((ds["residual"] / sig_for_chi) ** 2).sum("channel") + ds["chi_sq"] = ((ds["residual"] / sig_eff) ** 2).sum("channel") nu = max(b.shape[0] - rank_fit, 1) ds["red_chi_sq"] = ds["chi_sq"] / nu diff --git a/tests/test_qs_bridge.py b/tests/test_qs_bridge.py index f66fb19..6e87daf 100644 --- a/tests/test_qs_bridge.py +++ b/tests/test_qs_bridge.py @@ -130,3 +130,28 @@ def test_reconstruction_uses_minus_i_sign_convention(): f"reconstructed peak at φ={peak_phi}°, expected ~{delta}° (−i convention); " f"a peak near {360 - delta}° means the exp(+i…) sign bug is back" ) + + +def test_sigma_override_does_not_corrupt_fit_signal(synthetic_shot): + """Regression: fit_signal was de-normalized by the dataset signal_sigma while the + design matrix and RHS were normalized by the override sigma, scaling fit_signal + by signal_sigma/override and corrupting residual and chi_sq (the coefficients + stayed correct). With a uniform sigma the weighted LS solution is + override-invariant: fit_signal must be identical, and chi_sq must scale by + exactly (sigma_default / sigma_override)**2.""" + default = nodes._prep_qs_ds(synthetic_shot, {}).fit + override = 1e-4 + overridden = nodes._prep_qs_ds(synthetic_shot, {"sigma": str(override)}).fit + + np.testing.assert_allclose( + overridden["fit_coeffs"].values, default["fit_coeffs"].values, rtol=1e-9 + ) + np.testing.assert_allclose( + overridden["fit_signal"].values, default["fit_signal"].values, rtol=1e-9 + ) + sig0 = float(np.nanmean(default["signal_sigma"].values)) + np.testing.assert_allclose( + overridden["chi_sq"].values, + default["chi_sq"].values * (sig0 / override) ** 2, + rtol=1e-6, + ) From d9f58c46321983133a6faf00184a27ba4c655650 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 10 Jul 2026 01:01:49 -0400 Subject: [PATCH 05/16] qs_fit: restore OMFIT's 1e3 inversion cutoff as the fit_cond default fit_cond is 1/rcond for the lstsq inversion, but its default (10) was the GUI's "warn when K > 10" trust threshold, not OMFIT SLCONTOUR's 1e3 reference cutoff. Any fit with K in (10, 1e3) silently zeroed basis directions the reference fit keeps, and the reported K(eff) stayed <= 10 by construction - the quality panel looked healthiest exactly when the regularization was distorting the fit. Align the core, service, and GUI defaults on 1e3 and leave the trust threshold to contracts.quality_for_k. Co-Authored-By: Claude Fable 5 --- .../src/components/tabs/QuasiStationaryTab.tsx | 2 +- src/magnetics/core/qs_bridge.py | 4 ++-- src/magnetics/core/qs_fit.py | 7 ++++++- src/magnetics/service/nodes.py | 4 +++- tests/test_qs_bridge.py | 17 +++++++++++++++++ 5 files changed, 29 insertions(+), 5 deletions(-) diff --git a/gui/web/src/components/tabs/QuasiStationaryTab.tsx b/gui/web/src/components/tabs/QuasiStationaryTab.tsx index 15b186f..cc92814 100644 --- a/gui/web/src/components/tabs/QuasiStationaryTab.tsx +++ b/gui/web/src/components/tabs/QuasiStationaryTab.tsx @@ -171,7 +171,7 @@ export default function QuasiStationaryTab({ machine }: { machine: string }) { const [uncertainty, setUncertainty] = useState("2e-5"); const [energyFraction, setEnergyFraction] = useState("0.98"); const [fitBasis, setFitBasis] = useState("sinusoidal-integral"); - const [fitCond, setFitCond] = useState("10.0"); + const [fitCond, setFitCond] = useState("1000"); // OMFIT SLCONTOUR inversion cutoff (1/rcond), not the K>10 trust threshold const [cutoffLo, setCutoffLo] = useState("5.0"); const [cutoffHi, setCutoffHi] = useState("250.0"); diff --git a/src/magnetics/core/qs_bridge.py b/src/magnetics/core/qs_bridge.py index 3c42c88..301690b 100644 --- a/src/magnetics/core/qs_bridge.py +++ b/src/magnetics/core/qs_bridge.py @@ -248,7 +248,7 @@ def fit_to_fit_quality_node(fit_ds) -> dict: eff_cn = float(fit_ds.attrs.get("eff_cn", K)) mean_chi2 = float(np.nanmean(fit_ds["red_chi_sq"].values)) n_ch = int(fit_ds.sizes.get("channel", fit_ds["fit_sigmas"].shape[0])) - fit_cond = float(fit_ds.attrs.get("fit_condition", 10.0)) + fit_cond = float(fit_ds.attrs.get("fit_condition", 1e3)) n_modes = int(fit_ds.sizes["mode"]) return contracts.metrics( @@ -378,7 +378,7 @@ def fit_to_svd_condition_node(fit_ds) -> dict: w_a = np.linalg.svd(A, compute_uv=False) cond = np.abs(w_a[0] / w_a) idx = np.arange(1, len(cond) + 1) - fit_cond = float(fit_ds.attrs.get("fit_condition", 10.0)) + fit_cond = float(fit_ds.attrs.get("fit_condition", 1e3)) series = { "name": "condition number", diff --git a/src/magnetics/core/qs_fit.py b/src/magnetics/core/qs_fit.py index 17d4c88..cecc589 100644 --- a/src/magnetics/core/qs_fit.py +++ b/src/magnetics/core/qs_fit.py @@ -171,7 +171,7 @@ def fit( fit_exclude=(), fit_basis="sinusoidal-integral", fit_geometry="cylindrical", - fit_cond=10.0, + fit_cond=1e3, sigma_override=None, ncenters=6, mcenters=1, @@ -191,6 +191,11 @@ def fit( ``'gaussian-point'``, or ``'gaussian-integral'``. :param fit_geometry: ``'cylindrical'`` (phi, theta) or ``'vertical'`` (phi, z). :param fit_cond: condition-number cutoff for the lstsq inversion (= 1/rcond). + Basis directions with condition > fit_cond are zeroed. The default matches + OMFIT SLCONTOUR's 1e3 — this is the *inversion* cutoff, NOT the "warn when + K > 10" trust threshold (``contracts.quality_for_k``); setting it that low + silently truncates directions the reference fit would keep, while making + the reported K(eff) <= 10 by construction. :param sigma_override: when given, use this uniform measurement uncertainty for every channel instead of the per-channel ``signal_sigma`` baked into the dataset at load time. diff --git a/src/magnetics/service/nodes.py b/src/magnetics/service/nodes.py index f7f11ae..b294125 100644 --- a/src/magnetics/service/nodes.py +++ b/src/magnetics/service/nodes.py @@ -1446,7 +1446,9 @@ def _prep_qs_ds(shot, params): cutoff_hi = float(params.get("cutoff_hi", 250.0)) if params else 250.0 energy = float(params.get("energy", 0.98)) if params else 0.98 fit_basis = params.get("fit_basis", "sinusoidal-integral") if params else "sinusoidal-integral" - fit_cond = float(params.get("fit_cond", 10.0)) if params else 10.0 + # 1e3 = OMFIT SLCONTOUR's inversion cutoff (1/rcond); the "warn when K > 10" + # trust threshold lives in contracts.quality_for_k, not here. + fit_cond = float(params.get("fit_cond", 1e3)) if params else 1e3 sigma_str = params.get("sigma") if params else None sigma = float(sigma_str) if sigma_str is not None else None diff --git a/tests/test_qs_bridge.py b/tests/test_qs_bridge.py index 6e87daf..3475670 100644 --- a/tests/test_qs_bridge.py +++ b/tests/test_qs_bridge.py @@ -155,3 +155,20 @@ def test_sigma_override_does_not_corrupt_fit_signal(synthetic_shot): default["chi_sq"].values * (sig0 / override) ** 2, rtol=1e-6, ) + + +def test_default_fit_cond_is_inversion_cutoff_not_trust_threshold(synthetic_shot): + """The default cutoff must be OMFIT SLCONTOUR's 1e3 inversion cutoff (1/rcond) + at every layer. When it was 10 (the GUI's K-trust threshold), any fit with K in + (10, 1e3) silently zeroed basis directions the reference fit keeps — and K(eff) + read <= 10 by construction, so the quality panel looked healthiest exactly when + the regularization was distorting the fit.""" + import inspect + + from magnetics.core import qs_fit + + # core default (what a direct fit() call regularizes with) + assert inspect.signature(qs_fit.fit).parameters["fit_cond"].default == 1e3 + # service default (what a GUI request without fit_cond regularizes with) + run = nodes._prep_qs_ds(synthetic_shot, {}) + assert run.fit.attrs["fit_condition"] == 1e3 From 3fabb68dbc3d8e589ae440b90e0db42346a97bfd Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 10 Jul 2026 01:09:05 -0400 Subject: [PATCH 06/16] fetch: fetch KSTAR plasma tree_signals; guard pointname alias collisions Two silent-loss bugs on the tree paths: - The KSTAR _tree_transport branch never passed tree_signals to _fetch_mds_tree, so Ip/B_T/kappa were dropped from every KSTAR pull (absent from fetched AND missing lists) and helicity fell back to its default. Fetch them on the same connection after the sensors, first usable candidate wins, window-sliced (not resample()d onto the sensors' microsecond grid) - mirroring the NSTX path. - _resolve_pointnames mapped queried-pointname -> canonical id in a plain dict, so a per-era alias equal to another sensor's pointname (kstar.json: \MC1T10 -> \MC1P03 since shot 17377) collided: both fetched channels were relabeled to one name and the HDF5 writer silently dropped one. Query the node once, let the sensor literally named by the pointname keep it, and report the alias as skipped. KSTAR live behavior needs on-device verification (new plasma-node requests over the VPN); the changes are otherwise covered by offline stub-connection tests. Co-Authored-By: Claude Fable 5 --- src/magnetics/data/fetch/toksearch.py | 76 ++++++++++++++++-- tests/test_kstar_routing.py | 108 ++++++++++++++++++++++++++ 2 files changed, 178 insertions(+), 6 deletions(-) diff --git a/src/magnetics/data/fetch/toksearch.py b/src/magnetics/data/fetch/toksearch.py index cfbc697..46a169b 100644 --- a/src/magnetics/data/fetch/toksearch.py +++ b/src/magnetics/data/fetch/toksearch.py @@ -99,6 +99,26 @@ def _resolve_pointnames(dev, pointnames, shot): if pt is None: # out of range / NotAvailable skipped.append(cid) continue + prev = canonical_of.get(pt) + if prev is not None and prev != cid: + # A per-era alias collides with another sensor's queried pointname + # (e.g. KSTAR \MC1T10 -> \MC1P03 since shot 17377, while \MC1P03 is + # itself in the pull). Fetching one node under two canonical ids + # previously corrupted the relabel step — both channels got the same + # name and the writer silently dropped one. Query the node once, + # prefer the sensor literally named `pt` (identity mapping), and + # report the losing id as skipped instead of losing it silently. + if pt == cid: + skipped.append(prev) + canonical_of[pt] = cid + else: + skipped.append(cid) + sys.stderr.write( + f"pointname collision at shot {shot_i}: {skipped[-1]!r} aliases " + f"{pt!r}, which is fetched as {canonical_of[pt]!r}; skipping the " + "alias (check the device file's per-era pointname overrides).\n" + ) + continue query.append(pt) canonical_of[pt] = cid return query, canonical_of, skipped @@ -744,7 +764,9 @@ def _node_tree_map(dev): return node_tree, default_tree -def _fetch_mds_tree(shot, items, *, connect, tmin, tmax, progress, resample_dt=1e-6): +def _fetch_mds_tree( + shot, items, *, connect, tmin, tmax, progress, resample_dt=1e-6, tree_signals=None +): """Fetch tree nodes one-by-one over a single mdsthin connection. ``items`` is a list of ``(canonical_id, node, tree, gain)``: open each node's @@ -759,20 +781,28 @@ def _fetch_mds_tree(shot, items, *, connect, tmin, tmax, progress, resample_dt=1 ``resample(node, tmin, tmax, dt)`` when both bounds are given: it bounds the wire transfer (a raw 2 MHz channel is ~10^7 pts) at a cadence far above rotating-mode frequencies (kHz). ``resample_dt`` seconds sets that cadence (default 1 µs = - 1 MHz). With an open/absent window it falls back to a raw dimension slice.""" + 1 MHz). With an open/absent window it falls back to a raw dimension slice. + + ``tree_signals`` — plasma/equilibrium quantities ``{name: [(tree, node), …]}`` + (Ip, B_T, kappa from the device's "plasma pointnames"): fetched on the same + connection after the sensors, first usable candidate wins, mirroring + ``_fetch_mdsthin_tree``. Windowed with a raw dimension slice, NOT resample() — + these are slow (ms-cadence) traces and resampling them onto the sensors' µs + grid would interpolate them ~1000x.""" - def _slice(expr): + def _slice(expr, resample=True): if tmin is None and tmax is None: return expr lo = "*" if tmin is None else repr(float(tmin)) hi = "*" if tmax is None else repr(float(tmax)) - if tmin is not None and tmax is not None and resample_dt: + if resample and tmin is not None and tmax is not None and resample_dt: return f"resample({expr}, {lo}, {hi}, {float(resample_dt)!r})" return f"({expr})[{lo} : {hi}]" conn = connect() out: list[Channel] = [] - n = len(items) or 1 + tree_signals = tree_signals or {} + n = (len(items) + len(tree_signals)) or 1 current_tree = None for i, (canon, node, tree, gain) in enumerate(items, 1): try: @@ -805,6 +835,34 @@ def _slice(expr): current_tree = None # a failed openTree/get can leave the conn unsure out.append(Channel(canon, ok=False, error=str(exc))) progress(i / n, canon) + + for j, (name, candidates) in enumerate(tree_signals.items(), len(items) + 1): + chp = Channel(name, ok=False, error="no data") + for tree, node in candidates: + try: + if tree != current_tree: + conn.openTree(tree, int(shot)) + current_tree = tree + sliced = _slice(node, resample=False) + y = np.atleast_1d(np.asarray(conn.get(sliced).data(), dtype=np.float32)) + t = np.atleast_1d(np.asarray(conn.get(f"DIM_OF({sliced})").data())) + # A slow trace can legitimately hold few samples in a narrow window + # (helicity only needs the median sign), so allow n >= 1 here. + if y.size < 1 or t.size != y.size or not np.all(np.isfinite(t)): + chp.error = f"{tree}:{node}: degenerate result (n={y.size})" + continue + chp = Channel( + name, + (t * 1000.0).astype(np.float64), # tree time is s; h5 stores ms + y.astype(np.float32, copy=False), + ok=True, + ) + break + except Exception as exc: # noqa: BLE001 — plasma extras must not sink the pull + current_tree = None + chp.error = f"{tree}:{node}: {exc}" + out.append(chp) + progress(j / n, f"tree:{name} ({'ok' if chp.ok else 'missing'})") return out @@ -1478,7 +1536,13 @@ def _connect(): tmin_s = None if tmin is None else float(tmin) / 1000.0 tmax_s = None if tmax is None else float(tmax) / 1000.0 channels = _fetch_mds_tree( - shot_i, items, connect=_connect, tmin=tmin_s, tmax=tmax_s, progress=progress + shot_i, + items, + connect=_connect, + tmin=tmin_s, + tmax=tmax_s, + progress=progress, + tree_signals=tree_signals, ) elif tree_access: # NSTX/NSTX-U: sensors are fastmag tree nodes fetched with a value-window diff --git a/tests/test_kstar_routing.py b/tests/test_kstar_routing.py index 7e864bf..d20c142 100644 --- a/tests/test_kstar_routing.py +++ b/tests/test_kstar_routing.py @@ -69,3 +69,111 @@ def fake_remote(*a, **k): force=True, ) assert hit["remote"] is True, "DIII-D remote pull no longer reaches run_remote" + + +# --------------------------------------------------------------------------- +# tree_signals (Ip / B_T / kappa) must be fetched on the KSTAR tree path, and a +# per-era pointname alias that collides with another sensor must not silently +# drop a channel. Both run fully offline against stub connections. +# --------------------------------------------------------------------------- + +import numpy as np + +from magnetics.data.fetch.toksearch import _fetch_mds_tree, _resolve_pointnames + + +class _DefaultConn: + """Stub mdsthin connection: every get() returns a valid (data, DIM_OF) pair + unless the expression is in `fail` (then it raises). Records expressions.""" + + def __init__(self, fail=()): + self.fail = set(fail) + self.exprs = [] + self.opened = [] + + def openTree(self, tree, shot): # noqa: N802 (MDSplus API name) + self.opened.append(tree) + + def get(self, expr): + self.exprs.append(expr) + if any(f in expr for f in self.fail): + raise RuntimeError(f"TreeNNF: {expr}") + + class _R: + @staticmethod + def data(): + if expr.startswith("DIM_OF("): + return np.linspace(0.25, 0.35, 64) + return np.ones(64, np.float32) + + return _R() + + +def test_tree_signals_are_fetched_on_the_kstar_path(): + """Regression: the _tree_transport branch never passed tree_signals to + _fetch_mds_tree, so Ip/B_T/kappa were silently dropped for every KSTAR shot — + absent from channels_fetched AND channels_missing — and helicity fell back to + its default.""" + conn = _DefaultConn() + channels = _fetch_mds_tree( + 12345, + [("\\MC1T01", "\\MC1T01", "kstar", 1.0)], + connect=lambda: conn, + tmin=0.25, + tmax=0.35, + progress=lambda *_: None, + tree_signals={ + "Ip": [("PCS_KSTAR", "\\PC:PCRC03*(-1)*(1e-6)")], + "kappa": [("efitrt1", "\\kappa"), ("efitrt1", "\\top.results.aeqdsk:kappa")], + }, + ) + by_name = {c.name: c for c in channels} + assert by_name["Ip"].ok and by_name["kappa"].ok + # tree time is seconds; stored time must be ms + np.testing.assert_allclose(by_name["Ip"].time[0], 250.0) + # plasma traces are window-SLICED, never resample()d onto the sensor µs grid + ip_exprs = [e for e in conn.exprs if "PCRC03" in e and not e.startswith("DIM_OF")] + assert ip_exprs == ["(\\PC:PCRC03*(-1)*(1e-6))[0.25 : 0.35]"] + # the sensor still goes through resample() + assert any(e.startswith("resample(\\MC1T01") for e in conn.exprs) + + +def test_tree_signal_candidate_fallback_and_missing_are_reported(): + conn = _DefaultConn(fail=("\\kappa",)) # first kappa candidate fails + channels = _fetch_mds_tree( + 12345, + [], + connect=lambda: conn, + tmin=None, + tmax=None, + progress=lambda *_: None, + tree_signals={ + "kappa": [("efitrt1", "\\kappa"), ("efitrt1", "\\top.results.aeqdsk:kappa")], + "B_T": [("PCS_KSTAR", "\\kappa_not_there")], # also matches fail + }, + ) + by_name = {c.name: c for c in channels} + assert by_name["kappa"].ok # second candidate won + assert not by_name["B_T"].ok # recorded as missing, not lost + assert "B_T" in {c.name for c in channels if not c.ok} + + +def test_alias_collision_is_skipped_and_identity_wins(): + """Regression: kstar.json maps \\MC1T10 -> pointname \\MC1P03 (since shot 17377) + while \\MC1P03 is itself in the pull; both fetched channels got relabeled to the + same name and the HDF5 writer silently dropped one. The node must be queried + once, the sensor literally named \\MC1P03 keeps it, and the alias is reported + as skipped — in either request order.""" + dev = { + "sensors": { + "\\MC1P03": {"segments": [{"since_shot": 0, "phi": 10.0}]}, + "\\MC1T10": { + "segments": [{"since_shot": 17377, "phi": 160.5, "pointname": "\\MC1P03"}] + }, + } + } + for order in (["\\MC1T10", "\\MC1P03"], ["\\MC1P03", "\\MC1T10"]): + query, canon, skipped = _resolve_pointnames(dev, order, 20000) + assert query == ["\\MC1P03"], order + assert canon == {"\\MC1P03": "\\MC1P03"}, order + assert skipped == ["\\MC1T10"], order From 43de267e7dea11c590da00e3116cf67835d9a696 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 10 Jul 2026 01:11:41 -0400 Subject: [PATCH 07/16] fetch: stop caching transport failures as missing channels (#63) A dropped tunnel/ControlMaster mid-pull marked every remaining channel ok=False, and _write_h5 recorded them in channels_missing - so the next same-window pull treated them as already-attempted and silently skipped them, leaving the shot file permanently incomplete without --force. Classify failures at every fetch site (transport vs no_data, by exception type with message-word fallback) and exclude transport failures from the cached missing set; they are retried on the next pull while genuine no-data channels stay cached. Closes #63. Co-Authored-By: Claude Fable 5 --- src/magnetics/data/fetch/toksearch.py | 55 +++++++++++++++++++++++---- tests/test_toksearch_writer.py | 41 ++++++++++++++++++++ 2 files changed, 88 insertions(+), 8 deletions(-) diff --git a/src/magnetics/data/fetch/toksearch.py b/src/magnetics/data/fetch/toksearch.py index 46a169b..584291f 100644 --- a/src/magnetics/data/fetch/toksearch.py +++ b/src/magnetics/data/fetch/toksearch.py @@ -233,16 +233,44 @@ def _default_progress(frac: float, msg: str) -> None: # --- channel record ----------------------------------------------------------- class Channel: - """One fetched pointname: time (ms, float64) + data (float32).""" + """One fetched pointname: time (ms, float64) + data (float32). - __slots__ = ("name", "time", "data", "ok", "error") + ``error_kind`` classifies a failure (``ok=False``): ``"no_data"`` — the server + answered and the channel genuinely has nothing (cached as missing so the next + incremental pull skips it); ``"transport"`` — the connection/tunnel died (NOT + cached, so the next pull retries without --force; see issue #63).""" - def __init__(self, name, time=None, data=None, ok=False, error=""): + __slots__ = ("name", "time", "data", "ok", "error", "error_kind") + + def __init__(self, name, time=None, data=None, ok=False, error="", error_kind="no_data"): self.name = name self.time = time self.data = data self.ok = ok self.error = error + self.error_kind = error_kind + + +def _error_kind(exc) -> str: + """Classify a fetch exception: ``"transport"`` for connection-level failures + (socket/tunnel/timeout — worth retrying next run), ``"no_data"`` for a genuine + MDSplus no-data / no-node answer (cacheable as missing). Message matching is + the fallback for errors wrapped in strings by intermediate layers.""" + if isinstance(exc, (OSError, EOFError, TimeoutError, ConnectionError)): + return "transport" + msg = str(exc).lower() + transport_words = ( + "connection", + "timed out", + "timeout", + "broken pipe", + "reset by peer", + "eof", + "socket", + "tunnel", + "ssh", + ) + return "transport" if any(w in msg for w in transport_words) else "no_data" def _reduce(t: np.ndarray, y: np.ndarray, tmin, tmax, stride: int): @@ -402,6 +430,7 @@ def _mdsthin_tree_channels(connect, shot, tree_signals, tmin, tmax, progress): conn.openTree(tree, shot) except Exception as exc: ch.error = f"open {tree}: {exc}" + ch.error_kind = _error_kind(exc) continue try: y = np.atleast_1d(conn.get(node).data()) @@ -418,6 +447,7 @@ def _mdsthin_tree_channels(connect, shot, tree_signals, tmin, tmax, progress): ch.error = f"{tree}{node}: degenerate (n={y.size}, t={t.size})" except Exception as exc: ch.error = f"{tree}{node}: {exc}" + ch.error_kind = _error_kind(exc) out.append(ch) progress(1.0, f"tree:{name} ({'ok' if ch.ok else 'missing'})") return out @@ -548,7 +578,9 @@ def fetch_batch(conn, batch, use_many): t = _arr(gm.get(f"t{i}")) out.append(Channel(pt, t, y.astype(np.float32, copy=False), ok=True)) except Exception as exc: # per-channel error in the batch - out.append(Channel(pt, ok=False, error=str(exc))) + out.append( + Channel(pt, ok=False, error=str(exc), error_kind=_error_kind(exc)) + ) tick(batch[-1], k=len(batch)) return out, True except Exception as exc: # whole-batch failure -> fall back @@ -563,7 +595,7 @@ def fetch_batch(conn, batch, use_many): else: out.append(Channel(pt, ok=False, error="no data")) except Exception as exc: - out.append(Channel(pt, ok=False, error=str(exc))) + out.append(Channel(pt, ok=False, error=str(exc), error_kind=_error_kind(exc))) tick(pt) return out, False @@ -711,6 +743,7 @@ def fetch_all(conn): ch.error = f"degenerate (n={y.size}, t={t.size})" except Exception as exc: ch.error = str(exc) + ch.error_kind = _error_kind(exc) out.append(ch) tick(node) # plasma / equilibrium tree signals (efit01, ...): first candidate that opens @@ -734,6 +767,7 @@ def fetch_all(conn): chp.error = f"{tree}{node}: degenerate" except Exception as exc: chp.error = f"{tree}{node}: {exc}" + chp.error_kind = _error_kind(exc) out.append(chp) progress(1.0, f"tree:{name} ({'ok' if chp.ok else 'missing'})") return out @@ -833,7 +867,7 @@ def _slice(expr, resample=True): ) except Exception as exc: # noqa: BLE001 — one bad node shouldn't sink the pull current_tree = None # a failed openTree/get can leave the conn unsure - out.append(Channel(canon, ok=False, error=str(exc))) + out.append(Channel(canon, ok=False, error=str(exc), error_kind=_error_kind(exc))) progress(i / n, canon) for j, (name, candidates) in enumerate(tree_signals.items(), len(items) + 1): @@ -861,6 +895,7 @@ def _slice(expr, resample=True): except Exception as exc: # noqa: BLE001 — plasma extras must not sink the pull current_tree = None chp.error = f"{tree}:{node}: {exc}" + chp.error_kind = _error_kind(exc) out.append(chp) progress(j / n, f"tree:{name} ({'ok' if chp.ok else 'missing'})") return out @@ -1136,9 +1171,13 @@ def _write_h5( g.attrs[k] = float(v) # Union new channels with whatever the file already recorded; a name that - # is now fetched is removed from the missing list. + # is now fetched is removed from the missing list. Transport failures + # (dropped tunnel mid-pull, issue #63) are deliberately NOT recorded as + # missing — caching them would make every later incremental pull skip + # them, leaving the file permanently incomplete without --force. + no_data = {c.name for c in missing if getattr(c, "error_kind", "no_data") != "transport"} fetched = _attr_names(h5, "channels_fetched") | {c.name for c in got} - missing_names = (_attr_names(h5, "channels_missing") | {c.name for c in missing}) - fetched + missing_names = (_attr_names(h5, "channels_missing") | no_data) - fetched h5.attrs["channels_fetched"] = np.array(sorted(fetched), dtype="S") h5.attrs["channels_missing"] = np.array(sorted(missing_names), dtype="S") return got, missing diff --git a/tests/test_toksearch_writer.py b/tests/test_toksearch_writer.py index b716d55..668a4c9 100644 --- a/tests/test_toksearch_writer.py +++ b/tests/test_toksearch_writer.py @@ -291,3 +291,44 @@ def fake_run_remote(shot, analysis="both", *, local_out_dir=None, progress=None, assert "MPI66M307D" in h5, "remote pull clobbered the existing channel" assert "betan" in h5 assert not (tmp_path / "_remote_stage" / "shot_184927.h5").exists() + + +# --------------------------------------------------------------------------- +# issue #63 — transport failures must not poison the incremental cache +# --------------------------------------------------------------------------- + + +def test_transport_failures_are_not_cached_as_missing(tmp_path): + """A channel that failed because the tunnel died must be retried on the next + incremental pull: it may appear in the run's report, but never in the file's + channels_missing attr. Genuine no-data channels ARE cached.""" + import h5py + + from magnetics.data.fetch.toksearch import _existing_channels + + t = np.linspace(0.0, 10.0, 100) + got = Channel("GOOD", t, np.ones(100, np.float32), ok=True) + dead = Channel("DROPPED", ok=False, error="Connection reset by peer", error_kind="transport") + empty = Channel("EMPTY", ok=False, error="no data") # default kind: no_data + + out, (ok_ch, missing) = _write(tmp_path, [got, dead, empty]) + + # the run report still surfaces both failures... + assert {c.name for c in missing} == {"DROPPED", "EMPTY"} + with h5py.File(out, "r") as h5: + cached = {x.decode() for x in h5.attrs["channels_missing"]} + # ...but only the genuine no-data channel is cached as missing + assert cached == {"EMPTY"} + # so the next same-window pull retries DROPPED and skips GOOD + EMPTY + assert _existing_channels(out, None, None, 1) == {"GOOD", "EMPTY"} + + +def test_error_kind_classifier(): + from magnetics.data.fetch.toksearch import _error_kind + + assert _error_kind(ConnectionResetError("reset")) == "transport" + assert _error_kind(TimeoutError()) == "transport" + assert _error_kind(OSError(32, "Broken pipe")) == "transport" + assert _error_kind(RuntimeError("ssh tunnel exited")) == "transport" + assert _error_kind(RuntimeError("%TREE-W-NODATA, no data available")) == "no_data" + assert _error_kind(RuntimeError("TreeNNF: node not found")) == "no_data" From e349f4ea9151a6970d0f4dc8568c70defe668884 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 10 Jul 2026 01:14:08 -0400 Subject: [PATCH 08/16] spectral: scale spectrogram cross-power to csd's one-sided density compute_spectrogram used raw |conj(S2)*S1| with no window-power or fs scaling, yet rms_by_mode applied the same sqrt(sum(power)*df) formula as cross_spectrum (whose csd is a proper density) - so the n-spectrum axis labeled "rms amplitude" was off by a data-dependent factor of thousands. Apply csd's density scaling (1/(fs*sum(win^2)), doubled off DC/Nyquist); a unit sine now reports 0.707 RMS through both estimators. Log-power views only shift by a constant; coherence and phase are ratios/angles and are unaffected. Co-Authored-By: Claude Fable 5 --- src/magnetics/core/spectral.py | 19 +++++++++++++---- tests/test_spectral.py | 38 ++++++++++++++++++++++++++++++++++ 2 files changed, 53 insertions(+), 4 deletions(-) diff --git a/src/magnetics/core/spectral.py b/src/magnetics/core/spectral.py index 679070b..a41b5cf 100644 --- a/src/magnetics/core/spectral.py +++ b/src/magnetics/core/spectral.py @@ -338,9 +338,10 @@ def compute_spectrogram( Engine: one batched, single-precision short-time FFT per probe (no per-window loop), with the column count decimated to ``max_columns`` so cost scales with the display, - not the record length. Physics: cross-power = |conj(S1)·S2|, frequency-smoothed - coherence, n = round(phase / delta_phi), and per-mode RMS — identical to the 2-point - definitions in ``cross_spectrum``. + not the record length. Physics: cross-power = |conj(S1)·S2| scaled to a one-sided + cross-spectral density (csd's ``scaling="density"``), frequency-smoothed coherence, + n = round(phase / delta_phi), and per-mode RMS — identical to the 2-point + definitions in ``cross_spectrum``, so amplitudes are comparable across the two. Inputs: time (ndarray): sample times (s); assumed uniformly sampled. @@ -383,7 +384,17 @@ def compute_spectrogram( # order matches cross_spectrum's scipy.signal.csd(sig2, sig1) convention so the # spectrogram and the single-window 2-point analysis report the same *signed* n. cross = np.conj(spec2) * spec1 - power = np.abs(cross) + # Scale |cross| to a one-sided cross-spectral DENSITY exactly like scipy's + # csd(scaling="density"): 1/(fs·Σwin²), doubled off the DC/Nyquist bins. Without + # this, rms_by_mode's sqrt(Σ power·df) is off by a data-dependent factor of + # thousands vs the 2-point cross_spectrum RMS it claims to match. (Coherence and + # phase are ratios/angles — unaffected.) + scale = 1.0 / (sample_rate * float(np.sum(win.astype(np.float64) ** 2))) + power = np.abs(cross) * scale + if n_fft % 2 == 0: + power[:, 1:-1] *= 2.0 + else: + power[:, 1:] *= 2.0 phase = np.rad2deg(np.angle(cross)) # Coherence needs averaging; smooth the auto/cross spectra over frequency bins. diff --git a/tests/test_spectral.py b/tests/test_spectral.py index a10adfd..7b1535b 100644 --- a/tests/test_spectral.py +++ b/tests/test_spectral.py @@ -815,3 +815,41 @@ def test_all_families_agree_on_signed_n(self, array_mode): assert n_2pt == n_specgram == n_fit == n_array == n_mac == self.N_TRUE, ( f"2pt={n_2pt} specgram={n_specgram} fit={n_fit} array={n_array} mac={n_mac}" ) + + +# ----------------------------------------------------------------------- +# spectrogram power normalization — must match csd's density scaling +# ----------------------------------------------------------------------- + + +class TestSpectrogramPowerNormalization: + """Regression: compute_spectrogram used raw |conj(S2)·S1| with no window-power + or fs scaling yet applied cross_spectrum's density-RMS formula, so its + 'rms amplitude' was off by a data-dependent factor of thousands.""" + + def _pair(self, amp=1.0, fs=50_000, f_mode=3_000.0, duration=0.2): + t = np.linspace(0, duration, int(fs * duration), endpoint=False) + w = 2 * np.pi * f_mode + n, phi1, phi2 = 2, 30.0, 63.0 + sig1 = amp * np.sin(w * t - np.deg2rad(n * phi1)) + sig2 = amp * np.sin(w * t - np.deg2rad(n * phi2)) + return t, sig1, sig2, float(fs), phi2 - phi1 + + def test_rms_by_mode_recovers_physical_rms(self): + # a unit-amplitude sine has RMS 1/sqrt(2); both estimators must report it + t, sig1, sig2, fs, dphi = self._pair(amp=1.0) + two_pt = cross_spectrum(sig1, sig2, fs, delta_phi=dphi, nperseg=1024) + spec = compute_spectrogram(t, sig1, sig2, dphi, slice_duration=0.02) + rms_2pt = float(np.max(two_pt.rms_by_mode)) + rms_spec = float(np.max(spec.rms_by_mode)) + assert rms_2pt == pytest.approx(2**-0.5, rel=0.05) + assert rms_spec == pytest.approx(2**-0.5, rel=0.05) + + def test_rms_scales_linearly_with_amplitude(self): + _, s1a, s2a, fs, dphi = self._pair(amp=1.0) + t, s1b, s2b, _, _ = self._pair(amp=3.0) + a = compute_spectrogram(t, s1a, s2a, dphi, slice_duration=0.02) + b = compute_spectrogram(t, s1b, s2b, dphi, slice_duration=0.02) + assert float(np.max(b.rms_by_mode)) == pytest.approx( + 3.0 * float(np.max(a.rms_by_mode)), rel=1e-3 + ) From 9dad3f258674dc634c3d966906118ead9bbba349 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 10 Jul 2026 01:21:22 -0400 Subject: [PATCH 09/16] service: refresh clears theta cache; QS bands via contract; NaN -> null Four service-seam fixes from the review: - refresh() cleared every lru_cache except _real_theta, so a re-pull with more probes kept serving the stale channel->theta map and the poloidal nodes 422'd until server restart. - The QS amplitude/phase_t 1-sigma bands rode in untyped meta.sigma, which export.py drops - the downloaded HDF5 lacked the band shown on screen. Move them to the contract's per-series lower/upper (already exported, same as the rotating mode-shape nodes) and read those in the GUI. - channel_usage marked kappa/ip/bt and the QS midplane array as "unused / droppable" although the theta* correction, QS helicity, and the QS fit consume them. - NaN/Inf anywhere in a node payload was an opaque 500 (starlette serializes with allow_nan=False). Raw-signal nodes now emit non-finite samples as JSON null (a Plotly gap) via contracts.json_finite; zrange/chi^2/K scalars are guarded; and quality_for_k now matches contract.ts on +/-inf. Co-Authored-By: Claude Fable 5 --- .../components/tabs/QuasiStationaryTab.tsx | 10 +-- src/magnetics/core/contracts.py | 22 +++++- src/magnetics/core/qs_bridge.py | 30 +++++--- src/magnetics/service/nodes.py | 54 ++++++++++++-- tests/test_contract_meta.py | 9 ++- tests/test_nodes.py | 74 +++++++++++++++++++ tests/test_qs_bridge.py | 18 +++-- 7 files changed, 184 insertions(+), 33 deletions(-) diff --git a/gui/web/src/components/tabs/QuasiStationaryTab.tsx b/gui/web/src/components/tabs/QuasiStationaryTab.tsx index cc92814..e41e84c 100644 --- a/gui/web/src/components/tabs/QuasiStationaryTab.tsx +++ b/gui/web/src/components/tabs/QuasiStationaryTab.tsx @@ -58,26 +58,26 @@ function lineTraces( node: LineNode, opts?: { visible?: boolean[]; opacity?: number[]; palette?: string[] }, ): Partial[] { - const sigma = node.meta?.sigma as number[][] | undefined; const pal = opts?.palette ?? LINE_PALETTE; const traces: Partial[] = []; node.series.forEach((s, i) => { const color = pal[i % pal.length]; - const sig = sigma?.[i]; const vis = opts?.visible?.[i] !== false; const opacity = opts?.opacity?.[i] ?? 1; - if (sig && vis) { + // ±1σ band from the contract's typed lower/upper fields (also what the + // HDF5 export writes, so the band on screen matches the downloaded data). + if (s.lower && s.upper && vis) { traces.push({ type: "scatter", mode: "lines", x: s.x, - y: s.y.map((v, j) => v + sig[j]), + y: s.upper, line: { width: 0, color }, showlegend: false, hoverinfo: "skip", opacity, } as Partial); traces.push({ type: "scatter", mode: "lines", x: s.x, - y: s.y.map((v, j) => v - sig[j]), + y: s.lower, fill: "tonexty", fillcolor: hexToRgba(color, 0.45 * opacity), line: { width: 0, color }, showlegend: false, hoverinfo: "skip", opacity, diff --git a/src/magnetics/core/contracts.py b/src/magnetics/core/contracts.py index d5c80b0..3bf1f8c 100644 --- a/src/magnetics/core/contracts.py +++ b/src/magnetics/core/contracts.py @@ -8,8 +8,11 @@ from __future__ import annotations +import math from typing import Any +import numpy as np + def _clean(d: dict[str, Any]) -> dict[str, Any]: """Drop None-valued optional keys so the JSON matches the TS optionals.""" @@ -98,11 +101,24 @@ def metrics(title, fields, *, meta=None) -> dict: def quality_for_k(k: float) -> str: """SLCONTOUR condition-number thresholds (warn > 10, error > 20). - Mirrors `qualityForK` in contract.ts so the GUI's traffic-light coloring and - the backend agree. + Mirrors `qualityForK` in contract.ts (which tests `!Number.isFinite(K)`) so + the GUI's traffic-light coloring and the backend agree — including ±inf. """ - if not (k == k) or k > 20: # NaN or too ill-conditioned + if not math.isfinite(k) or k > 20: # NaN/±inf or too ill-conditioned return "bad" if k > 10: return "warn" return "good" + + +def json_finite(a): + """Array → nested lists with every non-finite value as None (JSON null). + + The HTTP layer serializes with ``allow_nan=False``, so a single NaN sample in + a raw fetched channel otherwise turns a whole node response into an opaque + 500. The GUI's Plotly renderers treat null as a gap, which is the honest + picture of missing samples.""" + a = np.asarray(a, dtype=float) + out = a.astype(object) + out[~np.isfinite(a)] = None + return out.tolist() diff --git a/src/magnetics/core/qs_bridge.py b/src/magnetics/core/qs_bridge.py index 301690b..277d4ff 100644 --- a/src/magnetics/core/qs_bridge.py +++ b/src/magnetics/core/qs_bridge.py @@ -117,7 +117,8 @@ def fit_to_qs_fit_node( theta_grid = np.linspace(0, 360, 49) z = _reconstruct_grid(fit_ds, phi_grid, theta_grid, t_idx) * _T_TO_G - zmax = float(np.nanmax(np.abs(z))) or 1.0 + finite = np.abs(z[np.isfinite(z)]) + zmax = float(finite.max()) if finite.size and finite.max() > 0 else 1.0 overlay = None if sensor_phis is not None and sensor_thetas is not None: @@ -147,7 +148,7 @@ def fit_to_qs_fit_node( "n": int(round(ns[dominant_idx])), "m": int(round(ms[dominant_idx])), "condition_number": round(K, 2), - "red_chi_sq": round(chi2_t, 3), + "red_chi_sq": round(chi2_t, 3) if np.isfinite(chi2_t) else None, "shot": str(fit_ds.attrs.get("shot", "")), }, ) @@ -156,8 +157,11 @@ def fit_to_qs_fit_node( def fit_to_amplitude_node(fit_ds) -> dict: """Mode amplitude ± 1σ vs time for each fitted mode → LineNode. - meta.sigma[i] is the per-series error band consumed by lineTraces() in the GUI. - meta.legend_title is "n" when all poloidal m=0, else "m/n" (mirrors plot_fit_modes). + The ±1σ band rides in each series' typed ``lower``/``upper`` fields (the + contract's error-band channel, same as the rotating mode-shape nodes) so the + GUI band and the HDF5 export are the same data — not an untyped meta side + channel the exporter drops. meta.legend_title is "n" when all poloidal m=0, + else "m/n" (mirrors plot_fit_modes). """ from ..core import contracts @@ -168,17 +172,19 @@ def fit_to_amplitude_node(fit_ds) -> dict: sigmas = fit_ds["fit_sigmas"].values # [mode, time] complex (constant along time) series = [] - sigma_bands = [] for i, (n, m) in enumerate(zip(ns, ms_vals)): amp, amp_err, _, _ = _amp_phase(coeffs[i], sigmas[i]) + y = amp * _T_TO_G + err = amp_err * _T_TO_G series.append( { "name": _mode_label(n, m), "x": np.round(t_ms, 2).tolist(), - "y": np.round(amp * _T_TO_G, 4).tolist(), + "y": np.round(y, 4).tolist(), + "lower": np.round(y - err, 4).tolist(), + "upper": np.round(y + err, 4).tolist(), } ) - sigma_bands.append(np.round(amp_err * _T_TO_G, 4).tolist()) legend_title = "n" if np.all(ms_vals == 0) else "m/n" K = float(fit_ds.attrs.get("condition_number", fit_ds.attrs.get("raw_cn", 0.0))) @@ -186,7 +192,6 @@ def fit_to_amplitude_node(fit_ds) -> dict: series, {"x": "time (ms)", "y": "amplitude (G)"}, meta={ - "sigma": sigma_bands, "legend_title": legend_title, "condition_number": round(K, 2), "shot": str(fit_ds.attrs.get("shot", "")), @@ -215,7 +220,6 @@ def fit_to_phase_t_node(fit_ds) -> dict: phase_visible = [bool(a > 0.1 * max_amp) for a in p90_amps] series = [] - sigma_bands = [] for i, (n, m) in enumerate(zip(ns, ms_vals)): _, _, phase, phase_err = _amp_phase(coeffs[i], sigmas[i]) series.append( @@ -223,16 +227,16 @@ def fit_to_phase_t_node(fit_ds) -> dict: "name": _mode_label(n, m), "x": np.round(t_ms, 2).tolist(), "y": np.round(phase, 3).tolist(), + "lower": np.round(phase - phase_err, 3).tolist(), + "upper": np.round(phase + phase_err, 3).tolist(), } ) - sigma_bands.append(np.round(phase_err, 3).tolist()) K = float(fit_ds.attrs.get("condition_number", fit_ds.attrs.get("raw_cn", 0.0))) return contracts.line( series, {"x": "time (ms)", "y": "phase (deg)"}, meta={ - "sigma": sigma_bands, "phase_visible": phase_visible, "condition_number": round(K, 2), "shot": str(fit_ds.attrs.get("shot", "")), @@ -297,7 +301,9 @@ def fit_to_phi_t_node(fit_ds, theta_fixed_deg: float = 0.0, n_phi: int = 73) -> z += (coeffs[i][None, :] * phase_phi * phase_theta).real z *= _T_TO_G - zmax = float(np.nanpercentile(np.abs(z), 99)) or 1.0 + finite = np.abs(z[np.isfinite(z)]) + zmax = float(np.percentile(finite, 99)) if finite.size else 1.0 + zmax = zmax or 1.0 K = float(fit_ds.attrs.get("condition_number", fit_ds.attrs.get("raw_cn", 0.0))) # z shape: [n_phi, n_time] but ContourNode z is [n_y][n_x] diff --git a/src/magnetics/service/nodes.py b/src/magnetics/service/nodes.py index b294125..97c2fbc 100644 --- a/src/magnetics/service/nodes.py +++ b/src/magnetics/service/nodes.py @@ -111,6 +111,39 @@ def tag(names, role): tag([n for n, _ in _poloidal_arr(str(shot))], "poloidal array") except ValueError: pass + # The quasi-stationary fit's default sensor set (same resolution as _prep_qs_ds). + try: + from ..core import qs_device + + dev, _ = _device_key(str(shot)) + default_cf = "Bp_LFS_midplane" + if dev: + default_cf = dev.get("qs_default_set") or ( + "quasi_stationary" + if "quasi_stationary" in dev.get("sensor_sets", {}) + else default_cf + ) + pats = np.atleast_1d( + qs_device.resolve_channel_filter( + default_cf, h5source.meta(shot).get("device", "DIII-D") + ) + ) + tag( + [nm for nm in all_names if any(re.match(str(p), nm) for p in pats)], + "QS fit array", + ) + except Exception: # noqa: BLE001 — QS not applicable to this shot/device + pass + # Plasma context channels: not part of any probe array, but consumed by the + # analyses — dropping them silently degrades the physics (θ* falls back to + # geometric θ; QS helicity falls back to its default) with no error anywhere. + for nm, role in ( + ("kappa", "θ* elongation correction (poloidal nodes)"), + ("ip", "QS helicity (sign of Ip·Bt)"), + ("bt", "QS helicity (sign of Ip·Bt)"), + ): + if nm in all_names: + tag([nm], role) used = [{"name": nm, "roles": roles[nm]} for nm in all_names if nm in roles] unused = [nm for nm in all_names if nm not in roles] @@ -126,7 +159,15 @@ def tag(names, role): def refresh() -> None: """Forget cached state (call after a new fetch writes a file).""" h5source.refresh() - for fn in (_spec_result, _stack_cached, _array_spectrum, _array_mode_spec, _qs_run, _dev_geom): + for fn in ( + _spec_result, + _stack_cached, + _array_spectrum, + _array_mode_spec, + _qs_run, + _dev_geom, + _real_theta, # channel→θ map derives from channel_names: stale after a re-pull + ): fn.cache_clear() @@ -643,12 +684,13 @@ def _toroidal_grid(shot): # ── contour: raw δBp(φ, t) — x=φ, y=time (QS contour hero plot) ────────────── def _contour(shot, params=None) -> dict: t_sub, phi_grid, z, phis, n_ch = _toroidal_grid(shot) - zmax = float(np.nanmax(np.abs(z))) or 1.0 + finite = np.abs(z[np.isfinite(z)]) + zmax = float(finite.max()) if finite.size and finite.max() > 0 else 1.0 overlay = {"points": [{"x": float(p), "y": float(t_sub[0])} for p in phis], "symbol": "square"} return contracts.contour( phi_grid.tolist(), t_sub.tolist(), - z.tolist(), + contracts.json_finite(z), {"x": "φ (deg)", "y": "time (ms)", "z": "δBp (G)"}, zrange=[-zmax, zmax], overlay=overlay, @@ -698,7 +740,7 @@ def _toroidal_stripes(shot, params=None) -> dict: return contracts.heatmap( t_sub.tolist(), ang.tolist(), - z.tolist(), + contracts.json_finite(z), {"x": "time (ms)", "y": "φ (deg)", "z": "δBp (a.u.)"}, discrete=False, meta={ @@ -722,7 +764,7 @@ def _poloidal_stripes(shot, params=None) -> dict: return contracts.heatmap( t_sub.tolist(), ang.tolist(), - z.tolist(), + contracts.json_finite(z), {"x": "time (ms)", "y": "θ (deg)", "z": "δBp (a.u.)"}, discrete=False, meta={ @@ -754,7 +796,7 @@ def _raw_trace(shot, params=None) -> dict: sel = np.arange(max(0, c - 2000), min(t_ms.size, c + 2000)) if sel.size > 2000: # keep the line light sel = sel[np.linspace(0, sel.size - 1, 2000).astype(int)] - series = [{"name": name, "x": t_ms[sel].tolist(), "y": d[sel].tolist()}] + series = [{"name": name, "x": t_ms[sel].tolist(), "y": contracts.json_finite(d[sel])}] return contracts.line( series, {"x": "time (ms)", "y": "dB/dt (a.u.)"}, diff --git a/tests/test_contract_meta.py b/tests/test_contract_meta.py index 35a03e6..a1bb8f7 100644 --- a/tests/test_contract_meta.py +++ b/tests/test_contract_meta.py @@ -33,12 +33,17 @@ def test_signal_conditioning_meta_pairs(synthetic_shot): def test_amplitude_meta_fields(synthetic_shot): meta = _meta(synthetic_shot, "amplitude") - assert "sigma" in meta and "legend_title" in meta + assert "legend_title" in meta + # the ±1σ band rides in the typed per-series lower/upper, not meta.sigma + node = nodes.build_node(synthetic_shot, "amplitude") + assert all("lower" in s and "upper" in s for s in node["series"]) def test_phase_t_meta_fields(synthetic_shot): meta = _meta(synthetic_shot, "phase_t") - assert "sigma" in meta and "phase_visible" in meta + assert "phase_visible" in meta + node = nodes.build_node(synthetic_shot, "phase_t") + assert all("lower" in s and "upper" in s for s in node["series"]) def test_phase_fit_reports_n_estimate(synthetic_shot): diff --git a/tests/test_nodes.py b/tests/test_nodes.py index ece9701..07b78cd 100644 --- a/tests/test_nodes.py +++ b/tests/test_nodes.py @@ -313,3 +313,77 @@ def test_quality_for_k_thresholds(): assert contracts.quality_for_k(15) == "warn" assert contracts.quality_for_k(25) == "bad" assert contracts.quality_for_k(float("nan")) == "bad" + + +# --------------------------------------------------------------------------- +# NaN samples in fetched channels must serialize as JSON null, not 500 +# --------------------------------------------------------------------------- + + +def test_raw_trace_with_nan_samples_serializes(synthetic_shot, monkeypatch): + """starlette serializes with allow_nan=False, so a NaN gap in a raw channel + used to turn the whole node response into an opaque 500. Non-finite samples + must come out as JSON null (a Plotly gap).""" + import json + + from magnetics.data import h5source + + orig = h5source.load_channel + + def nan_load(shot, name): + t, d = orig(shot, name) + d = np.asarray(d, dtype=float).copy() + d[::7] = np.nan + return t, d + + monkeypatch.setattr(nodes.h5source, "load_channel", nan_load) + try: + node = nodes.build_node(synthetic_shot, "raw_trace") + json.dumps(node, allow_nan=False) # must not raise + assert any(v is None for v in node["series"][0]["y"]) + assert any(v is not None for v in node["series"][0]["y"]) + finally: + nodes.refresh() # don't leak NaN-poisoned caches into other tests + + +def test_json_finite_helper(): + from magnetics.core.contracts import json_finite + + assert json_finite([1.0, np.nan, np.inf, -np.inf, 2.5]) == [1.0, None, None, None, 2.5] + assert json_finite([[1.0, np.nan], [3.0, 4.0]]) == [[1.0, None], [3.0, 4.0]] + + +def test_quality_for_k_matches_ts_on_nonfinite(): + """contract.ts uses !Number.isFinite(K) -> "bad"; the Python twin must agree + for NaN and BOTH infinities (K = -inf used to return "good").""" + assert contracts.quality_for_k(float("nan")) == "bad" + assert contracts.quality_for_k(float("inf")) == "bad" + assert contracts.quality_for_k(float("-inf")) == "bad" + assert contracts.quality_for_k(5.0) == "good" + assert contracts.quality_for_k(15.0) == "warn" + assert contracts.quality_for_k(25.0) == "bad" + + +def test_refresh_clears_real_theta_cache(synthetic_shot): + """Regression: refresh() (fired after every fetch job) cleared every lru_cache + except _real_theta, so a re-pull with more probes kept serving the stale + channel->theta map and poloidal nodes 422'd until server restart.""" + nodes._real_theta(str(synthetic_shot)) + assert nodes._real_theta.cache_info().currsize > 0 + nodes.refresh() + assert nodes._real_theta.cache_info().currsize == 0 + + +def test_channel_usage_tags_plasma_and_qs_channels(synthetic_shot): + """Regression: channel_usage marked kappa/ip/bt (and the QS fit array) as + "unused / droppable" although the theta* correction and the QS helicity + consume them — trimming a pull per that endpoint silently degraded physics.""" + usage = nodes.channel_usage(synthetic_shot) + unused = set(usage["unused"]) + present = {u["name"] for u in usage["used"]} + for nm in ("kappa", "ip", "bt"): + assert nm not in unused, f"{nm} reported droppable but analyses consume it" + assert nm in present + # the QS midplane array must carry a QS role, not sit in unused + roles = {u["name"]: u["roles"] for u in usage["used"]} + assert any("QS fit array" in r for rs in roles.values() for r in rs) diff --git a/tests/test_qs_bridge.py b/tests/test_qs_bridge.py index 3475670..d175934 100644 --- a/tests/test_qs_bridge.py +++ b/tests/test_qs_bridge.py @@ -58,8 +58,11 @@ def test_fit_quality_statuses_use_contract_vocabulary(fit_ds): 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) - assert np.all(np.isfinite(sigma)) + for series in node["series"]: + lower = np.asarray(series["lower"], dtype=float) + upper = np.asarray(series["upper"], dtype=float) + assert np.all(np.isfinite(lower)) and np.all(np.isfinite(upper)) + assert np.all(upper >= np.asarray(series["y"], dtype=float)) def test_svd_energy_node_is_monotonic_and_bounded(fit_ds): @@ -82,9 +85,14 @@ def test_svd_condition_node_is_finite_and_positive(fit_ds): def test_sigma_override_changes_amplitude_uncertainty(synthetic_shot): default_fit = nodes._prep_qs_ds(synthetic_shot, {}).fit overridden_fit = nodes._prep_qs_ds(synthetic_shot, {"sigma": "1.0"}).fit - default_sigma = qs_bridge.fit_to_amplitude_node(default_fit)["meta"]["sigma"] - overridden_sigma = qs_bridge.fit_to_amplitude_node(overridden_fit)["meta"]["sigma"] - assert np.mean(overridden_sigma) > np.mean(default_sigma) + + def band_width(fit): + node = qs_bridge.fit_to_amplitude_node(fit) + return np.mean( + [np.mean(np.asarray(s["upper"]) - np.asarray(s["lower"])) for s in node["series"]] + ) + + assert band_width(overridden_fit) > band_width(default_fit) def test_fit_basis_param_reaches_fit(synthetic_shot): From 885e994882026ff00540aa6249bfc57fa4859f41 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 10 Jul 2026 01:37:13 -0400 Subject: [PATCH 10/16] gui: honest error/staleness handling, retryable fetches, device-snap fixes Frontend correctness fixes from the review, verified with Playwright against the live app: - useNode's stale-while-revalidate is now param-scoped only: a machine switch drops the stale node instead of rendering the previous shot's data under the new shot's labels (and export descriptors). Adds a retryKey so an identical-key fetch can re-run: the QS Plot button and a new Retry button recover from transient failures (previously impossible without jiggling a parameter). - QS tab: only a 422 is diagnosed as "shot lacks the QS array"; other failures show a retryable error instead of a confidently wrong banner. Per-node failures surface in an error strip and in the section placeholders instead of eternal "loading..."; the mode-map panel shows the server error instead of fake "Computing..." progress. Removed the hardcoded shot-specific 3140 ms cursor write; the Plot button's dirty highlight compares params by value. - PullControl snaps per-device defaults whenever the store's device changes - the left rail's picker previously left backend "remote" (and a DIII-D window) active for tree devices while the select displayed "mdsthin". - SensorsTab resets its per-device set overrides on machine change (they are keyed by set name; carrying them over rendered a permanently empty sensor scene). - Both heavy tabs use selective store subscriptions instead of re-rendering every Plotly panel on each credential-field keystroke. Co-Authored-By: Claude Fable 5 --- gui/web/src/components/PullControl.tsx | 25 +++- .../components/tabs/QuasiStationaryTab.tsx | 120 +++++++++++++----- gui/web/src/components/tabs/RotatingTab.tsx | 17 ++- gui/web/src/components/tabs/SensorsTab.tsx | 9 +- gui/web/src/lib/useNode.ts | 32 +++-- 5 files changed, 147 insertions(+), 56 deletions(-) diff --git a/gui/web/src/components/PullControl.tsx b/gui/web/src/components/PullControl.tsx index 9cb7e80..9da30b4 100644 --- a/gui/web/src/components/PullControl.tsx +++ b/gui/web/src/components/PullControl.tsx @@ -61,8 +61,7 @@ export default function PullControl() { // map, so it uses mdsthin over a NARROW window (its raw signals are ~15 MHz and // seconds long — a wide window is gigabytes). An NSTX-style tree device requires a // named sensor set; KSTAR's transport defaults to its declared arrays when none is - // chosen. Called from the device onChange handler — NOT a synchronous effect - // (react-hooks/set-state-in-effect). + // chosen. function snapDeviceDefaults(d: DeviceInfo) { const tree = d.access === "mdsplus_tree"; setBackend(d.remote_capable ? "remote" : "mdsthin"); @@ -72,6 +71,22 @@ export default function PullControl() { if (d.default_shot != null) setShot(String(d.default_shot)); } + // The store owns `device`, and the left rail's Device picker (App.tsx) changes it + // WITHOUT going through this component's select — snapping only in the select's + // onChange left the backend/window/sensor-set at the previous device's values (a + // rail switch to a tree device then POSTed backend "remote" to a device with no + // cluster while the select displayed "mdsthin"). Snap once per device id, from + // wherever the change originated. + const snappedDevice = useRef(null); + useEffect(() => { + if (!device || snappedDevice.current === device) return; + snappedDevice.current = device; + const d = devices.find((x) => x.id === device); + // eslint-disable-next-line react-hooks/set-state-in-effect -- per-device defaults must follow an external (rail) device change + if (d) snapDeviceDefaults(d); + // eslint-disable-next-line react-hooks/exhaustive-deps + }, [device, devices]); + // close any open pull stream when the component unmounts useEffect(() => () => esRef.current?.close(), []); // Device is owned by the store; derive its capabilities. A tree device (NSTX/KSTAR) @@ -150,11 +165,7 @@ export default function PullControl() {

Pull a shot (live)

{devices.length > 0 && (