From c935105bece01dbf6731dd36b9e39e31795e7507 Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 1 Jul 2026 15:15:04 -0400 Subject: [PATCH 1/3] core+service+web: add MAC-based poloidal-m estimator with confidence MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Wire the Olofsson MAC (eigspec eq 9) into a poloidal mode-number estimator that answers even/odd (1/1 vs 2/1) with a probability, an independent cross-check on poloidal_phase_fit's slope fit. - core.mode_shape: mac_m_probability() Monte-Carlos the MAC over the per-probe phase σ to yield P(m); poloidal_m_spectrum() does the poloidal-specific prep — match the *phase pattern* (unit phasors, since a mode number is a phase winding and the per-probe amplitude spans 3-4 decades of probe sensitivity), SNR-gating out dead/saturated probes. - contracts: a `bar` kind (both sides) for the P(m) spectrum. - service: poloidal_m_spectrum node; verdict gated on the absolute MAC so a low match reads as "no clean poloidal mode", not a spurious m/n. - web: MAC+confidence tile (P(m) bars, nominal MAC overlay, verdict banner coloured by quality) in the rotating tab. On real DIII-D 174446 this resolves clean 3/2 and 1/3 modes (MAC 0.6-0.77) and honestly flags mode-free slices. Co-Authored-By: Claude Opus 4.8 --- gui/web/src/components/tabs/RotatingTab.tsx | 43 +++++++ gui/web/src/lib/NodeView.tsx | 36 ++++++ gui/web/src/lib/contract.ts | 14 +++ src/magnetics/core/contracts.py | 20 ++++ src/magnetics/core/mode_shape.py | 126 ++++++++++++++++++++ src/magnetics/service/nodes.py | 73 ++++++++++++ tests/test_mode_shape.py | 88 ++++++++++++++ tests/test_nodes.py | 21 ++++ tests/test_service_api.py | 2 +- 9 files changed, 422 insertions(+), 1 deletion(-) diff --git a/gui/web/src/components/tabs/RotatingTab.tsx b/gui/web/src/components/tabs/RotatingTab.tsx index 22de6fc..3c85241 100644 --- a/gui/web/src/components/tabs/RotatingTab.tsx +++ b/gui/web/src/components/tabs/RotatingTab.tsx @@ -201,6 +201,11 @@ export default function RotatingTab({ machine }: { machine: string }) { // of phase_fit. 422s when the shot lacks the poloidal array; the panel hides then. const { node: poloidalPhaseFitNode } = useNode(machine, "poloidal_phase_fit", { time: cursorMs }); + // MAC-based poloidal m with a Monte-Carlo confidence: an independent, shape-based + // cross-check on the slope fit above, and the even/odd (1/1 vs 2/1) verdict. 422s + // when the shot lacks the poloidal array; the panel hides then. + const { node: mSpectrumNode } = useNode(machine, "poloidal_m_spectrum", { time: cursorMs }); + // One Mirnov probe's raw dB/dt time series in a 4 ms window around the cursor. const { node: rawTraceNode } = useNode(machine, "raw_trace", { time: cursorMs }); @@ -1467,6 +1472,44 @@ export default function RotatingTab({ machine }: { machine: string }) { )} + {mSpectrumNode && mSpectrumNode.kind === "bar" && (() => { + const m = shapeMeta(mSpectrumNode) as Record | undefined; + const pEven = Number(m?.p_even ?? 0), pOdd = Number(m?.p_odd ?? 0); + const parityHi = pEven >= pOdd ? pEven : pOdd; + const parityName = pEven >= pOdd ? "even-m" : "odd-m"; + const quality = String(m?.quality ?? ""); + // colour the verdict by how well the shape matches a *pure* mode: a weak MAC + // means the printed m is untrustworthy, so it must not read as a firm answer. + const verdictColor = quality === "weak" ? "var(--warn, #ffb454)" + : quality === "marginal" ? "var(--text)" : "var(--accent)"; + return ( +
+

+ Poloidal m — MAC + Confidence + + {" · shape-based m (eigspec eq 9), P(m) from a Monte-Carlo over per-probe σ"} + +

+ {/* verdict banner — the plain-language even/odd (1/1 vs 2/1) answer */} +
+ {String(m?.verdict ?? "")} + + P({parityName}) = {parityHi.toFixed(2)} + {" · P(even) "}{pEven.toFixed(2)}{" / P(odd) "}{pOdd.toFixed(2)} + {" · MAC "}{Number(m?.mac_best ?? 0).toFixed(2)} + + + {Number(m?.n_used ?? m?.n_probes ?? 0)}/{Number(m?.n_probes ?? 0)} probes (SNR-gated) + {m?.used_theta_star ? ` · θ* (κ=${m?.kappa})` : " · geometric θ (no κ)"} + {m?.n_draws ? ` · ${Number(m.n_draws)} draws` : ""} + +
+
+ +
+
+ ); + })()} {analysisCard("Mode Persistence", modeTrackNode, "line", 200, shapeMeta(modeTrackNode)?.dominant_n != null ? `shape similarity to the dominant mode vs time (1 = persists) · n≈${shapeMeta(modeTrackNode)!.dominant_n}` diff --git a/gui/web/src/lib/NodeView.tsx b/gui/web/src/lib/NodeView.tsx index ac0b090..e72f9dc 100644 --- a/gui/web/src/lib/NodeView.tsx +++ b/gui/web/src/lib/NodeView.tsx @@ -206,6 +206,42 @@ export default function NodeView({ node, height }: { node: Node; height?: number ); } + case "bar": { + // Categorical spectrum (e.g. P(m) over poloidal mode number). The highlighted + // category (best-fit mode) is drawn in the accent colour, the rest muted; an + // optional secondary series (nominal MAC) rides on top as a marker line. + const xs = node.x.map((v) => String(v)); + const colors = node.x.map((v) => + node.highlight != null && v === node.highlight ? "#ff5cad" : "#4aa3ff", + ); + const traces: Partial[] = [ + { + type: "bar", x: xs, y: node.y, marker: { color: colors }, + name: node.axes.y, hovertemplate: `%{x}: %{y:.3f}`, + } as Partial, + ]; + if (node.secondary) { + traces.push({ + type: "scatter", mode: "lines+markers", x: xs, y: node.secondary.y, + name: node.secondary.name, yaxis: "y2", + marker: { color: "#ffb454", size: 6 }, line: { color: "#ffb454", width: 1, dash: "dot" }, + hovertemplate: `%{x}: %{y:.3f}`, + } as Partial); + } + const layout: Partial = { + ...axisLayout(node.axes), + showlegend: false, + bargap: 0.35, + }; + if (node.secondary) { + layout.yaxis2 = { + overlaying: "y", side: "right", range: [0, 1], showgrid: false, + title: { text: node.secondary.name, font: { size: 9 } }, tickfont: { size: 8 }, + } as Partial; + } + return ; + } + default: { // A kind outside the known union (backend shipped a new node before the // contract was updated, or a malformed payload). Render a placeholder — diff --git a/gui/web/src/lib/contract.ts b/gui/web/src/lib/contract.ts index 927b5cc..da06408 100644 --- a/gui/web/src/lib/contract.ts +++ b/gui/web/src/lib/contract.ts @@ -76,6 +76,19 @@ export interface LineNode { meta?: Record; } +/** A categorical bar chart — e.g. the MAC / probability spectrum over mode number. + * `x` are the categories (mode numbers), `y` the bar heights; `highlight` marks the + * best-fit category and `secondary` is an optional parallel series (nominal MAC). */ +export interface BarNode { + kind: "bar"; + x: (number | string)[]; + y: number[]; + axes: Axes; + highlight?: number | string | null; + secondary?: { name: string; y: number[] } | null; + meta?: Record; +} + /** A plasma equilibrium slice — normalized poloidal flux ψ_N(R,Z) + boundary. * Mirrors an EFIT gEQDSK slice: `psi_n` is row-major [z][r], 0 at the magnetic * axis and 1 at the last closed flux surface. Time-parametrized (one slice). */ @@ -107,6 +120,7 @@ export type Node = | HeatmapNode | Scatter2DNode | LineNode + | BarNode | EquilibriumNode | MetricsNode; diff --git a/src/magnetics/core/contracts.py b/src/magnetics/core/contracts.py index d5c80b0..bd559b6 100644 --- a/src/magnetics/core/contracts.py +++ b/src/magnetics/core/contracts.py @@ -66,6 +66,26 @@ def line(series, axes, *, meta=None) -> dict: return _clean({"kind": "line", "series": series, "axes": axes, "meta": meta}) +def bar(x, y, axes, *, highlight=None, secondary=None, meta=None) -> dict: + """A categorical bar chart — e.g. the MAC / probability spectrum over mode number. + + ``x`` are the categories (mode numbers), ``y`` the bar heights. ``highlight`` marks + one category (the best-fit mode) for emphasis; ``secondary`` is an optional parallel + series shown as an overlaid marker (e.g. the nominal MAC behind the probability). + """ + return _clean( + { + "kind": "bar", + "x": x, + "y": y, + "axes": axes, + "highlight": highlight, + "secondary": secondary, + "meta": meta, + } + ) + + def equilibrium( r, z, psi_n, boundary, axis, *, levels=None, time_ms=0.0, axes=None, meta=None ) -> dict: diff --git a/src/magnetics/core/mode_shape.py b/src/magnetics/core/mode_shape.py index 2a8d70f..cbe330d 100644 --- a/src/magnetics/core/mode_shape.py +++ b/src/magnetics/core/mode_shape.py @@ -310,6 +310,132 @@ def mac_n_spectrum( return ns, macs, int(ns[int(np.argmax(macs))]) +@dataclass(slots=True) +class MacMResult: + kind: str # "mac_m_spectrum" + m_values: NDArray[np.integer] # candidate poloidal mode numbers + mac_nominal: NDArray[np.floating] # MAC of the measured shape vs each e^{−imθ} template + p_by_m: NDArray[np.floating] # P(argmax MAC = this m) over the Monte-Carlo draws + best_m: int # argmax of the nominal (noise-free) MAC spectrum + p_best: float # posterior probability mass on ``best_m`` + p_even: float # Σ P over even |m| (m=0 counted even) → even-parity confidence + p_odd: float # Σ P over odd |m| → odd-parity confidence + n_draws: int # Monte-Carlo draws used (0 → deterministic, no per-probe σ) + + +def _mac_spectrum_matrix(templates_conj, tt_norm, z): + """MAC of one complex shape ``z`` against every template row at once. + + ``templates_conj[m] = conj(e^{−imθ})`` and ``tt_norm[m] = template†template`` are + precomputed so a Monte-Carlo loop is a single matmul per draw, not a Python loop + over modes. Returns the MAC value per template row (each in [0, 1]).""" + num = np.abs(templates_conj @ z) ** 2 # |vᵢ†z|² per row + den = tt_norm * float(np.real(np.vdot(z, z))) # (vᵢ†vᵢ)(z†z) + return np.where(den > 0.0, num / np.where(den > 0.0, den, 1.0), 0.0) + + +def mac_m_probability( + angle_deg: NDArray[np.floating], + complex_shape: NDArray[np.complexfloating], + *, + value_noise: NDArray[np.floating] | None = None, + m_range: tuple[int, int] = (-6, 6), + n_draws: int = 400, + seed: int = 0, +) -> MacMResult: + """Shape-based poloidal mode number m with a Monte-Carlo confidence (eigspec eq 9). + + Matches the measured (φ-detrended) poloidal shape vector against ideal templates + e^{−imθ} — the poloidal twin of :func:`mac_n_spectrum` — then, given the per-probe + 1σ ``value_noise`` (from :func:`shape_noise`), draws ``n_draws`` complex-Gaussian + realizations of the shape and records how often each m wins the MAC. That yields a + posterior P(m) rather than a bare peak, so an even/odd (e.g. 1/1 vs 2/1) call comes + with a probability instead of a single integer. A fixed ``seed`` keeps it + deterministic and cacheable. With no ``value_noise`` it degrades to the nominal peak. + """ + theta = np.deg2rad(np.asarray(angle_deg, dtype=np.float64)) + z = np.asarray(complex_shape, dtype=np.complex128) + ms = np.arange(m_range[0], m_range[1] + 1) + templates_conj = np.exp(1j * np.outer(ms, theta)) # conj(e^{−imθ}) = e^{+imθ} + tt_norm = np.full(ms.size, float(theta.size)) # |e^{±imθ}|² summed = n_probes + + mac_nominal = _mac_spectrum_matrix(templates_conj, tt_norm, z) + best_idx = int(np.argmax(mac_nominal)) + best_m = int(ms[best_idx]) + + if value_noise is None or n_draws <= 0: + p_by_m = (np.arange(ms.size) == best_idx).astype(np.float64) + draws = 0 + else: + sig = np.asarray(value_noise, dtype=np.float64) + rng = np.random.default_rng(seed) + counts = np.zeros(ms.size, dtype=np.float64) + for _ in range(n_draws): + noise = (rng.standard_normal(z.size) + 1j * rng.standard_normal(z.size)) * ( + sig / np.sqrt(2.0) + ) + counts[int(np.argmax(_mac_spectrum_matrix(templates_conj, tt_norm, z + noise)))] += 1.0 + p_by_m = counts / float(n_draws) + draws = int(n_draws) + + even = (np.abs(ms) % 2) == 0 + return MacMResult( + kind="mac_m_spectrum", + m_values=ms, + mac_nominal=mac_nominal, + p_by_m=p_by_m, + best_m=best_m, + p_best=float(p_by_m[best_idx]), + p_even=float(p_by_m[even].sum()), + p_odd=float(p_by_m[~even].sum()), + n_draws=draws, + ) + + +def poloidal_m_spectrum( + angle_deg: NDArray[np.floating], + phase_deg: NDArray[np.floating], + amplitude: NDArray[np.floating], + phase_error_deg: NDArray[np.floating] | None = None, + *, + snr_gate: float = 0.05, + m_range: tuple[int, int] = (-6, 6), + n_draws: int = 400, + seed: int = 0, +) -> tuple[MacMResult, int, int]: + """Poloidal mode number from the array's *phase pattern*, SNR-gated (eigspec eq 9). + + A mode number is a phase winding (δB ∝ e^{i(mθ−nφ)}), so identify m from the + *unit-amplitude* phasor pattern e^{iφ_k}: on a poloidal array the per-probe amplitude + is set by probe sensitivity and flux-surface geometry — spanning 3–4 decades — not by + the mode, so an amplitude-weighted MAC is captured by one or two large probes. Probes + below ``snr_gate``·max|z| are dropped (dead/saturated channels contribute only random + phase); the per-probe phase σ (rad) drives the Monte-Carlo confidence. Falls back to + the full set when too few probes clear the gate. Returns (result, n_used, n_total). + """ + amp = np.abs(np.asarray(amplitude, dtype=np.float64)) + ang = np.asarray(angle_deg, dtype=np.float64) + ph = np.deg2rad(np.asarray(phase_deg, dtype=np.float64)) + n_total = int(amp.size) + + keep = amp >= snr_gate * amp.max() if amp.max() > 0 else np.ones(amp.size, dtype=bool) + if int(keep.sum()) < 4: # too few survive the gate — keep everything rather than fail + keep = np.ones(amp.size, dtype=bool) + + z_unit = np.exp(1j * ph[keep]) # phase pattern only — amplitude discarded + noise = None + if phase_error_deg is not None: + s = np.deg2rad(np.asarray(phase_error_deg, dtype=np.float64))[keep] + ok = np.isfinite(s) & (s > 0) + fill = float(np.median(s[ok])) if np.any(ok) else 0.1 + noise = np.where(ok, s, fill) + + res = mac_m_probability( + ang[keep], z_unit, value_noise=noise, m_range=m_range, n_draws=n_draws, seed=seed + ) + return res, int(keep.sum()), n_total + + # --------------------------------------------------------------------------- # Time-resolved mode tracking (eigspec figure 9) # --------------------------------------------------------------------------- diff --git a/src/magnetics/service/nodes.py b/src/magnetics/service/nodes.py index c5810f3..abbf826 100644 --- a/src/magnetics/service/nodes.py +++ b/src/magnetics/service/nodes.py @@ -916,6 +916,78 @@ def _poloidal_phase_fit(shot, params=None) -> dict: ) +def _poloidal_m_spectrum(shot, params=None) -> dict: + """Shape-based poloidal mode number m via the Modal Assurance Criterion (eigspec + eq 9), with a Monte-Carlo confidence over the per-probe phase σ. + + Independent cross-check on ``poloidal_phase_fit``'s slope fit: match the measured + (φ-detrended, θ*-corrected) poloidal *phase pattern* against ideal e^{−imθ} templates + — after SNR-gating dead probes (see :func:`mode_shape.poloidal_m_spectrum`) — and, + propagating each probe's phase σ, report P(m). That gives an even/odd (e.g. 1/1 vs + 2/1) call with a probability rather than a bare integer. Bars are P(m); the nominal + MAC rides behind as a secondary marker, and ``quality`` gates the verdict on how well + the shape matches a *pure* mode at all.""" + arr, mode, f_khz, t0_ms, kappa = _poloidal_mode(shot, params) + res, n_used, n_total = mode_shape.poloidal_m_spectrum( + mode.toroidal_angle, mode.phase, mode.amplitude, mode.phase_error + ) + + # toroidal n (when a time cursor is set) to name the m/n mode in the verdict + n_tor = None + if t0_ms is not None: + try: + n_tor = _toroidal_n(str(shot), t0_ms * 1e-3, f_khz) + except ValueError, KeyError: + n_tor = None + + best_m = res.best_m + parity = "even" if best_m % 2 == 0 else "odd" + p_parity = res.p_even if best_m % 2 == 0 else res.p_odd + # absolute MAC of the winner: how well the measured shape actually matches a *pure* + # e^{−imθ} mode. A confident P(m) on top of a low MAC just means "least-bad of poor + # fits" (no clean poloidal mode, or |m| beyond what this sparse array resolves), so + # gate the verdict on it rather than printing a spuriously precise m/n. + mac_best = float(res.mac_nominal[int(np.argmax(res.mac_nominal))]) + quality = "clean" if mac_best >= 0.6 else "marginal" if mac_best >= 0.3 else "weak" + mode_label = f"m/n = {abs(best_m)}/{abs(n_tor)}" if n_tor else f"|m| = {abs(best_m)}" + if quality == "weak": + verdict = f"no clean poloidal mode (best {mode_label}, MAC {mac_best:.2f})" + else: + verdict = f"{mode_label} (P={p_parity:.2f} {parity}-m, MAC {mac_best:.2f})" + + ms = [int(m) for m in res.m_values] + return contracts.bar( + ms, + [round(float(p), 4) for p in res.p_by_m], + {"x": "poloidal mode number m", "y": "P(m) — MAC posterior"}, + highlight=best_m, + secondary={"name": "MAC (nominal)", "y": [round(float(v), 4) for v in res.mac_nominal]}, + meta={ + "best_m": best_m, + "n": n_tor, + "verdict": verdict, + "mac_best": round(mac_best, 3), + "quality": quality, + "p_best": round(float(res.p_best), 3), + "p_even": round(float(res.p_even), 3), + "p_odd": round(float(res.p_odd), 3), + "p_by_m": {int(m): round(float(p), 3) for m, p in zip(res.m_values, res.p_by_m)}, + "n_draws": res.n_draws, + "n_probes": n_total, + "n_used": n_used, + "kappa": round(float(kappa), 3) if kappa is not None else None, + "used_theta_star": kappa is not None, + "f_kHz": f_khz, + "t0_ms": t0_ms, + "shot": str(shot), + "note": "MAC of the φ-detrended poloidal phase pattern vs e^{−imθ} templates " + f"({n_used}/{n_total} probes above the SNR gate); P(m) from a Monte-Carlo over " + "the per-probe phase σ. 'quality' reflects the absolute MAC — a low value " + "means no pattern matches a pure mode, so the m is not to be trusted.", + }, + ) + + def _gp_shape(mode): """GP-smoothed complex mode shape with Tier-1-seeded heteroscedastic noise.""" z = mode_shape.shape_vector(mode.phase, mode.amplitude) @@ -1336,6 +1408,7 @@ def _fit_residuals(shot, params=None) -> dict: "toroidal_stripes": _toroidal_stripes, "poloidal_stripes": _poloidal_stripes, "poloidal_phase_fit": _poloidal_phase_fit, + "poloidal_m_spectrum": _poloidal_m_spectrum, "raw_trace": _raw_trace, } diff --git a/tests/test_mode_shape.py b/tests/test_mode_shape.py index 9d73086..e9d799c 100644 --- a/tests/test_mode_shape.py +++ b/tests/test_mode_shape.py @@ -8,8 +8,10 @@ gp_mode_shape, gp_periodic_fit, mac, + mac_m_probability, mac_n_spectrum, mode_pattern_2d, + poloidal_m_spectrum, periodic_kernel, shape_vector, ) @@ -183,6 +185,92 @@ def test_agrees_with_cross_phase_fit(self): assert best == fit.n +class TestMacMProbability: + def _shape(self, theta_deg, m_true): + # φ-detrended poloidal shape: phase = −m·θ (what _poloidal_mode produces) + return shape_vector(-m_true * theta_deg, np.ones_like(theta_deg)) + + def test_peaks_and_parity_on_clean_mode(self): + theta = np.array([0.0, 40.0, 80.0, 130.0, 170.0, 220.0, 270.0, 320.0]) + for m_true, parity in ((1, "odd"), (2, "even"), (3, "odd")): + res = mac_m_probability( + theta, self._shape(theta, m_true), value_noise=np.full(theta.size, 0.05) + ) + assert abs(res.best_m) == m_true + assert res.mac_nominal.max() > 0.99 + hi = res.p_even if parity == "even" else res.p_odd + assert hi > 0.9 # a clean mode is confidently the right parity + + def test_probabilities_form_a_distribution(self): + theta = np.array([0.0, 45.0, 90.0, 150.0, 210.0, 300.0]) + res = mac_m_probability( + theta, self._shape(theta, 2), value_noise=np.full(theta.size, 0.1), n_draws=300 + ) + assert abs(res.p_by_m.sum() - 1.0) < 1e-9 + assert abs(res.p_even + res.p_odd - 1.0) < 1e-9 + assert res.n_draws == 300 + + def test_noise_spreads_the_posterior(self): + # a sparse array + heavy noise must NOT stay pinned at p=1 on one m + theta = np.array([0.0, 30.0, 70.0, 110.0]) + res = mac_m_probability( + theta, self._shape(theta, 2), value_noise=np.full(theta.size, 0.9), seed=1 + ) + assert res.p_best < 1.0 # honest uncertainty, not overconfident + + def test_deterministic_without_noise(self): + theta = np.array([0.0, 45.0, 90.0, 150.0, 210.0, 300.0]) + res = mac_m_probability(theta, self._shape(theta, 2), value_noise=None) + assert res.best_m == 2 and res.n_draws == 0 + assert res.p_by_m.max() == 1.0 # collapses to the nominal peak + + def test_seed_is_reproducible(self): + theta = np.array([0.0, 40.0, 95.0, 160.0, 250.0]) + args = dict(value_noise=np.full(theta.size, 0.3), seed=7, n_draws=200) + a = mac_m_probability(theta, self._shape(theta, 1), **args) + b = mac_m_probability(theta, self._shape(theta, 1), **args) + np.testing.assert_array_equal(a.p_by_m, b.p_by_m) + + +class TestPoloidalMSpectrum: + def test_phase_pattern_ignores_amplitude_spread(self): + # a clean m=2 phase winding with wild per-probe amplitudes (3 decades) — the + # amplitude-weighted MAC would be captured by the big probe; the phase-pattern + # one must still nail m=2. + theta = np.array([0.0, 40.0, 80.0, 130.0, 170.0, 220.0, 270.0, 320.0]) + phase = -2.0 * theta # deg + amp = np.array([1.0, 1000.0, 2.0, 0.5, 800.0, 1.0, 3.0, 0.7]) + res, n_used, n_total = poloidal_m_spectrum(theta, phase, amp) + assert abs(res.best_m) == 2 and res.mac_nominal.max() > 0.99 + assert n_total == 8 and n_used == 8 # nothing below the 5% gate here + + def test_snr_gate_drops_dead_probes(self): + # 5 clean m=1 probes (unevenly spaced, so no m=1↔m=−3 aliasing) + 4 dead ones + # (tiny amplitude, random phase). The gate must drop the dead ones so m=1 wins. + clean = np.array([0.0, 50.0, 110.0, 200.0, 300.0]) + dead = np.array([25.0, 150.0, 240.0, 330.0]) + theta = np.concatenate([clean, dead]) + phase = np.concatenate([-1.0 * clean, np.array([12.0, 200.0, 77.0, 300.0])]) + amp = np.concatenate([np.ones(5), np.full(4, 1e-4)]) + res, n_used, n_total = poloidal_m_spectrum(theta, phase, amp) + assert n_total == 9 and n_used == 5 # only the 5 real probes survive the gate + assert abs(res.best_m) == 1 and res.mac_nominal.max() > 0.99 + + def test_falls_back_when_too_few_survive(self): + # if fewer than 4 clear the gate, keep everything rather than fail + amp = np.array([1.0, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6]) + theta = np.array([0.0, 60.0, 120.0, 180.0, 240.0, 300.0]) + res, n_used, n_total = poloidal_m_spectrum(theta, -theta, amp) + assert n_used == 6 # gate would leave 1 → fall back to all 6 + + def test_carries_monte_carlo_confidence(self): + theta = np.array([0.0, 45.0, 90.0, 150.0, 210.0, 300.0]) + res, _, _ = poloidal_m_spectrum( + theta, -2.0 * theta, np.ones(6), np.full(6, 20.0), n_draws=200 + ) + assert res.n_draws == 200 and abs(res.p_by_m.sum() - 1.0) < 1e-9 + + class TestNoiseCalibration: def test_shape_noise_combines_phase_and_amplitude(self): from magnetics.core.mode_shape import shape_noise diff --git a/tests/test_nodes.py b/tests/test_nodes.py index 7786212..a08bfa5 100644 --- a/tests/test_nodes.py +++ b/tests/test_nodes.py @@ -220,6 +220,27 @@ def test_poloidal_phase_fit_node(): assert "m_fit" in n["meta"] +def test_poloidal_m_spectrum_node(): + shot = _first_shot() + try: + n = nodes.build_node(shot, "poloidal_m_spectrum", {"time": 3000}) + except Exception as e: # noqa: BLE001 — shot may lack the poloidal array + pytest.skip(f"no poloidal array in this shot: {e}") + assert n["kind"] == "bar" + assert len(n["x"]) == len(n["y"]) # one probability per candidate m + m = n["meta"] + assert "best_m" in m and "verdict" in m + assert 0.0 <= m["mac_best"] <= 1.0 + assert m["quality"] in {"clean", "marginal", "weak"} + assert 0 < m["n_used"] <= m["n_probes"] # SNR gate keeps a subset of the array + # P(m) is a proper distribution over the candidates (draws + per-probe σ present) + assert m["n_draws"] > 0 + assert abs(sum(n["y"]) - 1.0) < 1e-6 + assert 0.0 <= m["p_even"] <= 1.0 and 0.0 <= m["p_odd"] <= 1.0 + assert abs(m["p_even"] + m["p_odd"] - 1.0) < 1e-6 + assert n["secondary"]["y"] and len(n["secondary"]["y"]) == len(n["x"]) # nominal MAC + + def test_unknown_node_raises(): shot = _first_shot() with pytest.raises(KeyError): diff --git a/tests/test_service_api.py b/tests/test_service_api.py index e7aa69c..68960f1 100644 --- a/tests/test_service_api.py +++ b/tests/test_service_api.py @@ -16,7 +16,7 @@ client = TestClient(app) _GRID_KINDS = {"heatmap", "contour"} -_VALID_KINDS = {"contour", "heatmap", "scatter2d", "line", "metrics", "equilibrium"} +_VALID_KINDS = {"contour", "heatmap", "scatter2d", "line", "bar", "metrics", "equilibrium"} def test_machines_lists_the_synthetic_shot(synthetic_shot): From 3399fc96aba2ec34e894329a6b84ae4b2a2d8ee5 Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 1 Jul 2026 16:54:00 -0400 Subject: [PATCH 2/3] data: fetch + read the EFIT02 q-profile per shot MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add EFIT 2-D profile support so the q-profile (q vs ψ_N, in time) can be pulled and read back, for the rotating mode-number anchor. - signals: TREE_PROFILES / tree_profiles_for — q_profile from EFIT02 only (MSE-constrained; no efit01 fallback, so a non-MSE q is never silently substituted — the profile is simply omitted if efit02 lacks it). - toksearch: a Profile record + _mdsthin_tree_profiles (2-D fetch, orientation-robust, ψ_N = uniform gEQDSK grid), wired additively into the mdsthin backend and written as its own h5 group; the existing scalar/ pointname path is untouched. - h5source: load_q_profile() reader; channel_names() skips profile groups so they never masquerade as (data, time) channels. - synthetic_shot: a physical synthetic q-profile (q0≈1.05→q_edge≈4.5) so the read/anchor chain is exercised in CI without any tokamak data. The live mdsthin fetch needs a re-pull to validate against real EFIT02; toksearch/remote backends don't yet fetch profiles (follow-up). Co-Authored-By: Claude Opus 4.8 --- src/magnetics/data/fetch/toksearch.py | 111 +++++++++++++++++++++++++- src/magnetics/data/h5source.py | 22 ++++- src/magnetics/data/signals.py | 32 ++++++++ tests/synthetic_shot.py | 13 ++- 4 files changed, 173 insertions(+), 5 deletions(-) diff --git a/src/magnetics/data/fetch/toksearch.py b/src/magnetics/data/fetch/toksearch.py index f4ec417..5149e66 100644 --- a/src/magnetics/data/fetch/toksearch.py +++ b/src/magnetics/data/fetch/toksearch.py @@ -196,6 +196,82 @@ def __init__(self, name, time=None, data=None, ok=False, error=""): self.error = error +# --- EFIT profile record (2-D: value vs ψ_N, evolving in time) ---------------- +class Profile: + """One fetched EFIT profile: time (ms, float64) + psi_n axis + data[ntime, npsi].""" + + __slots__ = ("name", "time", "psi", "data", "ok", "error") + + def __init__(self, name, time=None, psi=None, data=None, ok=False, error=""): + self.name = name + self.time = time + self.psi = psi + self.data = data + self.ok = ok + self.error = error + + +def _mdsthin_tree_profiles(connect, shot, tree_profiles, tmin, tmax, progress): + """Fetch EFIT 2-D profiles (e.g. q(ψ_N, t)) on one connection (openTree + node). + + Mirrors :func:`_mdsthin_tree_channels` but for a 2-D signal: reads the array plus + its time axis (``dim_of``), orients it so axis 0 is time (EFIT gEQDSK QPSI comes as + either [ntime, nψ] or its transpose depending on the tree), and builds the ψ_N axis + as the uniform gEQDSK flux grid linspace(0, 1, nψ). Window applied client-side. Each + candidate (tree, node) is tried in order; the first usable 2-D array wins. Failures + are non-fatal — the profile is just marked missing.""" + out: list[Profile] = [] + if not tree_profiles: + return out + conn = connect() + for name, candidates in tree_profiles.items(): + pf = Profile(name, ok=False, error="no data") + for tree, node in candidates: + try: + conn.openTree(tree, shot) + except Exception as exc: + pf.error = f"open {tree}: {exc}" + continue + try: + raw = np.asarray(conn.get(node).data(), dtype=np.float64) + if raw.ndim != 2 or raw.size == 0: + pf.error = f"{tree}{node}: not 2-D (shape={raw.shape})" + continue + t = np.atleast_1d(conn.get(f"dim_of({node})").data()).astype(np.float64) + # orient so axis 0 = time (QPSI ships either way across EFIT trees) + if raw.shape[0] == t.size: + data = raw + elif raw.shape[1] == t.size: + data = raw.T + else: + pf.error = f"{tree}{node}: time dim {t.size} matches neither axis {raw.shape}" + continue + if not np.all(np.isfinite(t)): + pf.error = f"{tree}{node}: non-finite time" + continue + sel = np.ones(t.size, dtype=bool) + if tmin is not None: + sel &= t >= tmin + if tmax is not None: + sel &= t <= tmax + t, data = t[sel], data[sel] + psi = np.linspace(0.0, 1.0, data.shape[1]) # uniform gEQDSK flux grid + if t.size >= 1 and data.shape[1] >= 2: + pf = Profile( + name, + t.astype(np.float64, copy=False), + psi.astype(np.float64, copy=False), + data.astype(np.float32, copy=False), + ok=True, + ) + break + except Exception as exc: + pf.error = f"{tree}{node}: {exc}" + out.append(pf) + progress(1.0, f"profile:{name} ({'ok' if pf.ok else 'missing'})") + return out + + def _reduce(t: np.ndarray, y: np.ndarray, tmin, tmax, stride: int): """Client-side window + decimate fallback (used after a full fetch).""" if tmin is not None or tmax is not None: @@ -389,6 +465,7 @@ def _fetch_mdsthin( batch_size, progress, tree_signals=None, + tree_profiles=None, ssh_env=None, per_channel=False, ): @@ -551,6 +628,7 @@ def connect(): results = run_all(connect, f"{len(batches)} batches x{batch_size} ({workers} conns)") results += _mdsthin_tree_channels(connect, shot, tree_signals, tmin, tmax, progress) + results += _mdsthin_tree_profiles(connect, shot, tree_profiles, tmin, tmax, progress) return results # Off-network: one SSH tunnel, then a few TCP mdsip conns through it. with _ssh_tunnel(username, gateway, mds_host, mds_port, env=ssh_env) as lport: @@ -560,6 +638,7 @@ def connect(): results = run_all(connect, f"{len(batches)} batches x{batch_size} via tunnel") results += _mdsthin_tree_channels(connect, shot, tree_signals, tmin, tmax, progress) + results += _mdsthin_tree_profiles(connect, shot, tree_profiles, tmin, tmax, progress) return results @@ -876,6 +955,10 @@ def _write_h5( channel_geometry = channel_geometry or {} comp = None if compression in (None, "none") else compression + # EFIT 2-D profiles (q(ψ_N, t), ...) ride in the same results list but are stored + # as their own group, not as (data, time) channels — split them out first. + profiles = [c for c in channels if isinstance(c, Profile)] + channels = [c for c in channels if not isinstance(c, Profile)] missing = [c for c in channels if not c.ok] # A channel can fetch data samples yet come back with a degenerate time axis @@ -965,6 +1048,22 @@ def _write_h5( if v is not None: g.attrs[k] = float(v) + # EFIT profiles: one group per profile with data[ntime, nψ] + its own time (ms) + # and psi (ψ_N) axes. Stored outside the channel/timebase machinery so profile + # arrays never masquerade as (data, time) channels downstream. + for pf in profiles: + if not pf.ok: + continue + if pf.name in h5: # re-fetched: replace + del h5[pf.name] + pg = h5.create_group(pf.name) + pg.attrs["kind"] = "efit_profile" + pg.create_dataset("data", data=pf.data, compression=comp, chunks=True if comp else None) + pg.create_dataset("time", data=pf.time, compression=comp, chunks=True if comp else None) + pg.create_dataset("psi", data=pf.psi, compression=comp, chunks=True if comp else None) + pg.attrs["time_units"] = "ms" + pg.attrs["axis"] = "psi_n" + # Union new channels with whatever the file already recorded; a name that # is now fetched is removed from the missing list. fetched = _attr_names(h5, "channels_fetched") | {c.name for c in got} @@ -1176,6 +1275,10 @@ def fetch_shot( tree_signals = ms.tree_signals_for(analysis) label = analysis + # EFIT02 q-profile (and any other 2-D profiles) for this analysis — fetched only on + # the mdsthin backend; empty for analyses that don't request it. + tree_profiles = ms.tree_profiles_for(analysis) + # Shot-aware pointname resolution: map canonical ids -> the pointnames valid # at THIS shot, dropping channels the shot can't have so we never query them. shot_i = int(shot) @@ -1212,12 +1315,13 @@ def fetch_shot( ) else: merge = True - before = len(pointnames) + len(tree_signals) + before = len(pointnames) + len(tree_signals) + len(tree_profiles) pointnames = [p for p in pointnames if p not in existing] tree_signals = {k: v for k, v in tree_signals.items() if k not in existing} - n_skipped = before - (len(pointnames) + len(tree_signals)) + tree_profiles = {k: v for k, v in tree_profiles.items() if k not in existing} + n_skipped = before - (len(pointnames) + len(tree_signals) + len(tree_profiles)) - if merge and not pointnames and not tree_signals: + if merge and not pointnames and not tree_signals and not tree_profiles: sys.stderr.write(f"All {n_skipped} requested signals already in {out}; nothing to fetch.\n") return out @@ -1294,6 +1398,7 @@ def fetch_shot( batch_size=batch_size, progress=progress, tree_signals=tree_signals, + tree_profiles=tree_profiles, ssh_env=ssh_env, per_channel=per_channel, ) diff --git a/src/magnetics/data/h5source.py b/src/magnetics/data/h5source.py index 1aa903a..8aa9a12 100644 --- a/src/magnetics/data/h5source.py +++ b/src/magnetics/data/h5source.py @@ -128,7 +128,27 @@ def device_id(shot: str | int) -> str: def channel_names(shot: str | int) -> list[str]: with h5py.File(shot_file(shot), "r") as h5: - return [k for k in h5.keys() if k != "_timebases"] + # Exclude the timebase store and any EFIT profile groups (q(ψ_N, t), …): those + # are not (data, time) channels, so array-builders must never see them. + return [ + k + for k in h5.keys() + if k != "_timebases" and _attr_str(h5[k].attrs.get("kind"), "") != "efit_profile" + ] + + +def load_q_profile(shot: str | int): + """The EFIT q-profile as (time_ms float64, psi_n float64, q[ntime, npsi] float32), + or None if the shot was pulled without it (older pulls, or no EFIT02 q available).""" + with h5py.File(shot_file(shot), "r") as h5: + if "q_profile" not in h5: + return None + g = h5["q_profile"] + return ( + np.asarray(g["time"][:], dtype=np.float64), + np.asarray(g["psi"][:], dtype=np.float64), + np.asarray(g["data"][:], dtype=np.float32), + ) def load_channel(shot: str | int, name: str): diff --git a/src/magnetics/data/signals.py b/src/magnetics/data/signals.py index c8ea4f6..c2d79d2 100644 --- a/src/magnetics/data/signals.py +++ b/src/magnetics/data/signals.py @@ -301,6 +301,28 @@ "both": ["kappa"], } +# EFIT PROFILES (2-D: value vs normalized flux ψ_N, evolving in time) -- unlike the +# scalar TREE_SIGNALS above, these are q(ψ_N, t)-style arrays. Fetched from EFIT02 +# specifically: it is the MSE-constrained reconstruction, so its q-profile (hence the +# rational-surface / mode-number anchoring built on it) is trustworthy in the core, +# where a magnetics-only efit01 q can be off by enough to move a rational surface. We +# do NOT fall back to efit01 here: a silently-substituted non-MSE q would defeat the +# point, so if efit02 is absent the profile is simply omitted (analysis degrades to the +# unanchored MAC). The ψ_N axis is the uniform gEQDSK flux grid (linspace 0..1). +TREE_PROFILES: dict[str, list[tuple[str, str]]] = { + "q_profile": [ + ("efit02", r"\top.results.geqdsk:qpsi"), + ("efit02", r"\efit_geqdsk:qpsi"), + ("efit02", r"\qpsi"), + ], +} + +ANALYSIS_TREE_PROFILES: dict[str, list[str]] = { + "quasi-stationary": ["q_profile"], + "rotating": ["q_profile"], + "both": ["q_profile"], +} + def signals_for(analysis: str) -> list[str]: """Return the de-duplicated, ordered pointname list for an analysis type.""" @@ -335,6 +357,16 @@ def tree_signals_for(analysis: str) -> dict[str, list[tuple[str, str]]]: return {name: TREE_SIGNALS[name] for name in names} +def tree_profiles_for(analysis: str) -> dict[str, list[tuple[str, str]]]: + """Return {profile_name: [(tree, node), ...]} EFIT 2-D profiles for an analysis. + + Like :func:`tree_signals_for` but for ψ_N-resolved profiles (q(ψ_N, t)); fetched + from EFIT02 only (MSE-constrained). Absent for an unknown analysis rather than + raising — a missing q-profile is non-fatal (the mode fit degrades to unanchored). + """ + return {name: TREE_PROFILES[name] for name in ANALYSIS_TREE_PROFILES.get(analysis, [])} + + def decimate_allowed(analysis: str) -> bool: """Whether server-side decimation is permitted for this analysis type.""" if analysis not in REDUCTION: diff --git a/tests/synthetic_shot.py b/tests/synthetic_shot.py index ae01e8c..96d50a8 100644 --- a/tests/synthetic_shot.py +++ b/tests/synthetic_shot.py @@ -25,7 +25,7 @@ import numpy as np from magnetics.data import device_geom, devices, diiid, signals -from magnetics.data.fetch.toksearch import Channel, _write_h5, resolve_sensor_set +from magnetics.data.fetch.toksearch import Channel, Profile, _write_h5, resolve_sensor_set class _Mode(NamedTuple): @@ -112,6 +112,17 @@ def _channels(shot: int, *, include_qs: bool) -> list[Channel]: chans.append(Channel("ip", t_ms.copy(), ip, ok=True)) chans.append(Channel("bt", t_ms.copy(), np.full(n, 2.0, np.float32), ok=True)) chans.append(Channel("kappa", t_ms.copy(), np.full(n, 1.8, np.float32), ok=True)) + + if include_qs: + # A synthetic EFIT02-style q-profile: monotonic q0≈1.05 → q_edge≈4.5 on a + # uniform ψ_N grid, ~constant over a handful of coarse EFIT time slices. Gives + # the mode-number anchor real rational surfaces (q=2 at m/n=6/3, etc.) to test + # against — no tokamak data, just a physical shape. + q_t = np.linspace(0.0, _DURATION * 1e3, 5) # 5 EFIT slices (ms) + psi = np.linspace(0.0, 1.0, 65) + q1d = 1.05 + 3.45 * psi**2 # q(0)=1.05, q(1)=4.5 + q2d = np.tile(q1d, (q_t.size, 1)).astype(np.float32) # [ntime, npsi], steady + chans.append(Profile("q_profile", q_t, psi, q2d, ok=True)) return chans From b84ae4d20ea7c7fafc8a3b911932a4acb191081c Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 1 Jul 2026 16:54:17 -0400 Subject: [PATCH 3/3] rotating: anchor poloidal m to the EFIT02 q-profile + flag aliases MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The poloidal MAC winner is easily an alias on the sparse, uneven poloidal array (a physical high-m mode folds down to a spurious small one), so the bare fit could report an unphysical mode like m/n = 1/3 (q ≈ 0.33, inside the q ≈ 1 axis) with false confidence. Fix it three ways: - core.equilibrium (new): resonant_surfaces() + anchor_poloidal_m() — for the reliable toroidal n, keep only m whose q = m/n has a real rational surface (q_min ≤ m/n ≤ q_max) and, if the raw MAC winner is unphysical, swap in the best-MAC physical alias. Returns the q surface + its ψ_N. - core.mode_shape: mac_m_probability now also reports the MAC margin (best − next-best) and the alias set within a margin of the best. - service.nodes: poloidal_m_spectrum reads the q-profile (_q_at) and anchors m to it; a margin-aware quality ladder — weak < ambiguous (alias near-tie, no q to break it) < marginal < clean < q-anchored — so a near-tie reads as "ambiguous: |m| ∈ {…}" instead of a false-confident integer. - web: verdict banner shows the q surface / ψ_N when anchored, the alias set + margin otherwise, coloured by trust. On the synthetic shot the ground-truth 2/1 is recovered from an aliased MAC (q=2.00 @ ψ_N=0.52); real shots without a q-profile now honestly report the ambiguity (e.g. 174446 @ 3000 ms → "|m| ∈ {1,2,6}") instead of "1/3". Co-Authored-By: Claude Opus 4.8 --- gui/web/src/components/tabs/RotatingTab.tsx | 34 +++-- src/magnetics/core/equilibrium.py | 149 ++++++++++++++++++++ src/magnetics/core/mode_shape.py | 18 +++ src/magnetics/service/nodes.py | 112 ++++++++++++--- tests/test_equilibrium.py | 81 +++++++++++ tests/test_nodes.py | 30 +++- 6 files changed, 393 insertions(+), 31 deletions(-) create mode 100644 src/magnetics/core/equilibrium.py create mode 100644 tests/test_equilibrium.py diff --git a/gui/web/src/components/tabs/RotatingTab.tsx b/gui/web/src/components/tabs/RotatingTab.tsx index 2ae2a0c..274aec0 100644 --- a/gui/web/src/components/tabs/RotatingTab.tsx +++ b/gui/web/src/components/tabs/RotatingTab.tsx @@ -1485,10 +1485,14 @@ export default function RotatingTab({ machine }: { machine: string }) { const parityHi = pEven >= pOdd ? pEven : pOdd; const parityName = pEven >= pOdd ? "even-m" : "odd-m"; const quality = String(m?.quality ?? ""); - // colour the verdict by how well the shape matches a *pure* mode: a weak MAC - // means the printed m is untrustworthy, so it must not read as a firm answer. - const verdictColor = quality === "weak" ? "var(--warn, #ffb454)" + // colour the verdict by trust: q-anchored (physical rational surface) is the + // strongest; weak/ambiguous are cautions (the printed m may be a sampling alias) + // so they must not read as a firm answer. + const verdictColor = quality === "q-anchored" ? "var(--good, #54e08a)" + : quality === "weak" || quality === "ambiguous" ? "var(--warn, #ffb454)" : quality === "marginal" ? "var(--text)" : "var(--accent)"; + const qAnchored = m?.q_anchored === true; + const aliasMs = Array.isArray(m?.alias_ms) ? (m!.alias_ms as number[]) : []; return (

@@ -1497,18 +1501,26 @@ export default function RotatingTab({ machine }: { machine: string }) { {" · shape-based m (eigspec eq 9), P(m) from a Monte-Carlo over per-probe σ"}

- {/* verdict banner — the plain-language even/odd (1/1 vs 2/1) answer */} + {/* verdict banner — the plain-language m/n answer */}
{String(m?.verdict ?? "")} + {qAnchored ? ( + + EFIT02 q-anchored: q={Number(m?.q_surface ?? 0).toFixed(2)} + {m?.psi_n != null ? ` @ ψ_N=${Number(m.psi_n).toFixed(2)}` : ""} + {m?.q_corrected ? ` · raw MAC m=${Number(m?.raw_m ?? 0)} was an alias` : ""} + + ) : ( + + P({parityName}) = {parityHi.toFixed(2)} + {aliasMs.length ? ` · aliases m∈{${[Number(m?.raw_m ?? m?.best_m), ...aliasMs].join(", ")}}` : ""} + {` · margin ${Number(m?.margin ?? 0).toFixed(2)}`} + + )} - P({parityName}) = {parityHi.toFixed(2)} - {" · P(even) "}{pEven.toFixed(2)}{" / P(odd) "}{pOdd.toFixed(2)} - {" · MAC "}{Number(m?.mac_best ?? 0).toFixed(2)} - - - {Number(m?.n_used ?? m?.n_probes ?? 0)}/{Number(m?.n_probes ?? 0)} probes (SNR-gated) + MAC {Number(m?.mac_best ?? 0).toFixed(2)} + {" · "}{Number(m?.n_used ?? m?.n_probes ?? 0)}/{Number(m?.n_probes ?? 0)} probes (SNR-gated) {m?.used_theta_star ? ` · θ* (κ=${m?.kappa})` : " · geometric θ (no κ)"} - {m?.n_draws ? ` · ${Number(m.n_draws)} draws` : ""}
diff --git a/src/magnetics/core/equilibrium.py b/src/magnetics/core/equilibrium.py new file mode 100644 index 0000000..a187a21 --- /dev/null +++ b/src/magnetics/core/equilibrium.py @@ -0,0 +1,149 @@ +"""q-profile physics: rational surfaces and mode-number anchoring. + +A magnetic island / tearing mode lives on a *rational* flux surface where the +safety factor q = m/n is a ratio of the poloidal (m) and toroidal (n) mode numbers. +The toroidal n from a well-sampled toroidal array is reliable; the poloidal m from a +sparse, unevenly-sampled poloidal array is easily *aliased* (a physical m folds down +to a spurious small one). The q-profile breaks that degeneracy: for the measured n, +only integer m with a real q = m/n surface in the plasma (q_min ≤ m/n ≤ q_max) are +physically possible, so an aliased m/n below q_min (e.g. m/n = 1/3, q ≈ 0.33, inside +the q ≈ 1 axis) can be rejected and replaced by the physical alias. + +Pure numpy over arrays: no device specifics, no I/O. +""" + +from __future__ import annotations + +from dataclasses import dataclass + +import numpy as np +from numpy.typing import NDArray + + +@dataclass(slots=True) +class RationalSurface: + m: int # poloidal mode number of the resonance + n: int # toroidal mode number (the anchor) + q: float # q = m/n at the surface + psi_n: float # normalized flux ψ_N of the (outermost) crossing, in [0, 1] + + +@dataclass(slots=True) +class MAnchorResult: + chosen_m: int # physically-anchored poloidal mode number (signed, MAC's sign kept) + q_surface: float | None # q = |chosen_m|/|n| if it resonates, else None + psi_n: float | None # ψ_N of the resonant surface, if any + corrected: bool # True if the raw MAC winner was unphysical and got replaced + raw_m: int # the raw MAC-argmax m, before the q-anchor + allowed_ms: list[int] # |m| that have a q = m/n surface (within the MAC alias set) + q_min: float + q_max: float + + +def _crossings(psi_n: NDArray[np.floating], q: NDArray[np.floating], target: float) -> list[float]: + """ψ_N where the (possibly non-monotonic) q-profile crosses ``target`` q, by linear + interpolation across each sign change of q − target. Empty if never reached.""" + d = np.asarray(q, dtype=np.float64) - float(target) + psi = np.asarray(psi_n, dtype=np.float64) + out: list[float] = [] + sign_change = np.flatnonzero(np.diff(np.sign(d)) != 0) + for i in sign_change: + d0, d1 = d[i], d[i + 1] + if d1 != d0: + out.append(float(psi[i] - d0 * (psi[i + 1] - psi[i]) / (d1 - d0))) + return out + + +def q_range(q_psi: NDArray[np.floating]) -> tuple[float, float]: + """(q_min, q_max) over the profile, ignoring NaNs.""" + q = np.asarray(q_psi, dtype=np.float64) + return float(np.nanmin(q)), float(np.nanmax(q)) + + +def resonant_surfaces( + q_psi: NDArray[np.floating], + psi_n: NDArray[np.floating], + n: int, + *, + m_max: int = 15, +) -> list[RationalSurface]: + """All q = m/n rational surfaces present in the profile for toroidal number ``n``. + + Scans |m| = 1…``m_max``; a surface exists when q = m/|n| lies within [q_min, q_max]. + The ψ_N reported is the *outermost* crossing (largest ψ_N) — the mode's resonance is + conventionally the outer one for the higher-order rationals. Returns them sorted by m. + """ + na = abs(int(n)) + if na == 0: + return [] + qmin, qmax = q_range(q_psi) + out: list[RationalSurface] = [] + for m in range(1, m_max + 1): + qt = m / na + if qmin <= qt <= qmax: + xs = _crossings(psi_n, q_psi, qt) + psi_loc = max(xs) if xs else float("nan") + out.append(RationalSurface(m=m, n=na, q=qt, psi_n=psi_loc)) + return out + + +def anchor_poloidal_m( + m_values: NDArray[np.integer], + mac_values: NDArray[np.floating], + n: int, + q_psi: NDArray[np.floating], + psi_n: NDArray[np.floating], + *, + alias_margin: float = 0.1, + m_max: int = 15, +) -> MAnchorResult: + """Pick the physical poloidal m from the MAC spectrum using the q-profile. + + The raw MAC argmax is easily an alias on a sparse poloidal array. Among the MAC + *alias set* — the candidates within ``alias_margin`` MAC of the best — keep only those + whose |m|/|n| has a real q surface, and choose the highest-MAC survivor. If the raw + winner is itself physical, nothing changes; if it is not (e.g. m/n below q_min) it is + replaced and ``corrected`` is set. Sign of m follows the MAC (helicity), the q test + uses |m|. + """ + ms = np.asarray(m_values) + macs = np.asarray(mac_values, dtype=np.float64) + na = abs(int(n)) + qmin, qmax = q_range(q_psi) + + best_idx = int(np.argmax(macs)) + raw_m = int(ms[best_idx]) + + def resonates(m_abs: int) -> bool: + return na > 0 and 1 <= m_abs <= m_max and qmin <= m_abs / na <= qmax + + # alias set: candidates within the MAC margin of the best (excluding q=0/trivial) + alias_mask = macs >= macs[best_idx] - alias_margin + allowed_ms = sorted({abs(int(m)) for m in ms[alias_mask] if resonates(abs(int(m)))}) + + chosen_idx = best_idx + corrected = False + if not resonates(abs(raw_m)): + # raw winner is unphysical → take the best-MAC alias that does resonate + candidates = [i for i in np.argsort(macs)[::-1] if resonates(abs(int(ms[i])))] + if candidates: + chosen_idx = int(candidates[0]) + corrected = True + + chosen_m = int(ms[chosen_idx]) + q_surface = abs(chosen_m) / na if resonates(abs(chosen_m)) else None + psi_loc: float | None = None + if q_surface is not None: + xs = _crossings(psi_n, q_psi, q_surface) + psi_loc = max(xs) if xs else None + + return MAnchorResult( + chosen_m=chosen_m, + q_surface=q_surface, + psi_n=psi_loc, + corrected=corrected, + raw_m=raw_m, + allowed_ms=allowed_ms, + q_min=qmin, + q_max=qmax, + ) diff --git a/src/magnetics/core/mode_shape.py b/src/magnetics/core/mode_shape.py index cbe330d..637abfd 100644 --- a/src/magnetics/core/mode_shape.py +++ b/src/magnetics/core/mode_shape.py @@ -321,6 +321,8 @@ class MacMResult: p_even: float # Σ P over even |m| (m=0 counted even) → even-parity confidence p_odd: float # Σ P over odd |m| → odd-parity confidence n_draws: int # Monte-Carlo draws used (0 → deterministic, no per-probe σ) + margin: float # best MAC − next-best MAC (small ⇒ an alias near-tie, not a clean pick) + alias_ms: list[int] # other m whose MAC is within ``alias_margin`` of the best def _mac_spectrum_matrix(templates_conj, tt_norm, z): @@ -342,6 +344,7 @@ def mac_m_probability( m_range: tuple[int, int] = (-6, 6), n_draws: int = 400, seed: int = 0, + alias_margin: float = 0.1, ) -> MacMResult: """Shape-based poloidal mode number m with a Monte-Carlo confidence (eigspec eq 9). @@ -352,6 +355,10 @@ def mac_m_probability( posterior P(m) rather than a bare peak, so an even/odd (e.g. 1/1 vs 2/1) call comes with a probability instead of a single integer. A fixed ``seed`` keeps it deterministic and cacheable. With no ``value_noise`` it degrades to the nominal peak. + + Also reports the MAC ``margin`` (best − next-best) and the ``alias_ms`` set (other m + within ``alias_margin`` of the best): on a sparse poloidal array several m fit almost + equally, and a small margin flags that the winner is an alias tie, not a clean pick. """ theta = np.deg2rad(np.asarray(angle_deg, dtype=np.float64)) z = np.asarray(complex_shape, dtype=np.complex128) @@ -378,6 +385,15 @@ def mac_m_probability( p_by_m = counts / float(n_draws) draws = int(n_draws) + # margin to the next-best distinct m, and the alias set within `alias_margin` of best + others = np.delete(mac_nominal, best_idx) + margin = float(mac_nominal[best_idx] - others.max()) if others.size else 1.0 + alias_ms = [ + int(m) + for i, m in enumerate(ms) + if i != best_idx and mac_nominal[i] >= mac_nominal[best_idx] - alias_margin + ] + even = (np.abs(ms) % 2) == 0 return MacMResult( kind="mac_m_spectrum", @@ -385,6 +401,8 @@ def mac_m_probability( mac_nominal=mac_nominal, p_by_m=p_by_m, best_m=best_m, + margin=margin, + alias_ms=alias_ms, p_best=float(p_by_m[best_idx]), p_even=float(p_by_m[even].sum()), p_odd=float(p_by_m[~even].sum()), diff --git a/src/magnetics/service/nodes.py b/src/magnetics/service/nodes.py index 39aab00..4201cbb 100644 --- a/src/magnetics/service/nodes.py +++ b/src/magnetics/service/nodes.py @@ -18,7 +18,7 @@ import numpy as np -from ..core import contracts, geometry, mode_shape, qs_bridge, spectral +from ..core import contracts, equilibrium, geometry, mode_shape, qs_bridge, spectral from ..data import device_geom, diiid_geometry, h5source logger = logging.getLogger(__name__) @@ -180,6 +180,33 @@ def _kappa_at(shot, t0_ms=None): return k if 1.0 <= k <= 4.0 else None +def _q_at(shot, t0_ms=None): + """The EFIT02 q-profile q(ψ_N) at the cursor time as (psi_n, q), interpolated across + EFIT slices and clamped to the reconstruction's time range. None when the shot was + pulled without a q-profile (older pulls, or no EFIT02) — the poloidal m fit then runs + unanchored.""" + try: + prof = h5source.load_q_profile(shot) + except KeyError, OSError: + return None + if prof is None: + return None + t, psi, q = prof # t[ntime] ms, psi[npsi], q[ntime, npsi] + q = np.asarray(q, dtype=float) + if q.ndim != 2 or q.shape[0] == 0: + return None + tt = np.asarray(t, dtype=float) + if t0_ms is None or tt.size < 2: + qcol = np.nanmedian(q, axis=0) + else: + t0 = float(np.clip(t0_ms, tt.min(), tt.max())) + qcol = np.array([float(np.interp(t0, tt, q[:, j])) for j in range(q.shape[1])]) + good = np.isfinite(qcol) + if int(good.sum()) < 2: + return None + return np.asarray(psi, dtype=float)[good], qcol[good] + + @lru_cache(maxsize=16) def _stack_cached(shot, names): """Cached channel load for a fixed (shot, names) — HDF5 reads are the dominant @@ -1033,18 +1060,62 @@ def _poloidal_m_spectrum(shot, params=None) -> dict: except ValueError, KeyError: n_tor = None - best_m = res.best_m - parity = "even" if best_m % 2 == 0 else "odd" - p_parity = res.p_even if best_m % 2 == 0 else res.p_odd - # absolute MAC of the winner: how well the measured shape actually matches a *pure* - # e^{−imθ} mode. A confident P(m) on top of a low MAC just means "least-bad of poor - # fits" (no clean poloidal mode, or |m| beyond what this sparse array resolves), so - # gate the verdict on it rather than printing a spuriously precise m/n. + best_m = res.best_m # raw MAC argmax (alias-prone on a sparse array) mac_best = float(res.mac_nominal[int(np.argmax(res.mac_nominal))]) - quality = "clean" if mac_best >= 0.6 else "marginal" if mac_best >= 0.3 else "weak" - mode_label = f"m/n = {abs(best_m)}/{abs(n_tor)}" if n_tor else f"|m| = {abs(best_m)}" + + # q-profile anchor: the poloidal array is sparse and uneven, so the MAC winner is + # often an alias — a physical high-m mode folds down to a spurious small one. With the + # EFIT02 (MSE-constrained) q-profile, keep only m whose q = m/|n| has a real rational + # surface, and if the raw winner is unphysical (e.g. m/n = 1/3, q ≈ 0.33 inside the + # q ≈ 1 axis) swap in the best-MAC physical alias. n from the toroidal array anchors it. + qprof = _q_at(str(shot), t0_ms) if n_tor else None + anchor = None + if qprof is not None and n_tor: + anchor = equilibrium.anchor_poloidal_m( + res.m_values, res.mac_nominal, n_tor, q_psi=qprof[1], psi_n=qprof[0] + ) + # hoist the anchor fields into plain locals (None when unanchored) so the verdict + # logic and meta don't repeatedly re-narrow the optional dataclass. + chosen_m = anchor.chosen_m if anchor is not None else best_m + q_surface = anchor.q_surface if anchor is not None else None + q_psi_n = anchor.psi_n if anchor is not None else None + q_corrected = anchor.corrected if anchor is not None else False + raw_m = anchor.raw_m if anchor is not None else best_m + allowed_ms = list(anchor.allowed_ms) if anchor is not None else None + + parity = "even" if chosen_m % 2 == 0 else "odd" + p_parity = res.p_even if chosen_m % 2 == 0 else res.p_odd + q_ok = q_surface is not None + aliased = res.margin < 0.08 or bool(res.alias_ms) + + # quality ladder: weak (no pure-mode match at all) < ambiguous (alias near-tie, no q + # to break it) < marginal < clean. A q-anchor that finds a physical surface overrides + # the alias ambiguity — that's the whole point of bringing the equilibrium in. + if mac_best < 0.3: + quality = "weak" + elif q_ok: + quality = "q-anchored" + elif aliased: + quality = "ambiguous" + elif mac_best >= 0.6: + quality = "clean" + else: + quality = "marginal" + + mode_label = f"m/n = {abs(chosen_m)}/{abs(n_tor)}" if n_tor else f"|m| = {abs(chosen_m)}" + alias_set = sorted({abs(chosen_m), *(abs(a) for a in res.alias_ms)}) if quality == "weak": verdict = f"no clean poloidal mode (best {mode_label}, MAC {mac_best:.2f})" + elif q_surface is not None: # q-anchored + loc = f" @ ψ_N={q_psi_n:.2f}" if q_psi_n is not None else "" + verdict = f"{mode_label} — resonant q={q_surface:.2f}{loc}" + if q_corrected: + verdict += f" (raw MAC m={raw_m} was a sampling alias)" + elif quality == "ambiguous": + verdict = ( + f"ambiguous: |m| ∈ {alias_set} (MAC margin {res.margin:.2f}) — " + "sampling-limited; no EFIT02 q-profile to disambiguate" + ) else: verdict = f"{mode_label} (P={p_parity:.2f} {parity}-m, MAC {mac_best:.2f})" @@ -1053,14 +1124,22 @@ def _poloidal_m_spectrum(shot, params=None) -> dict: ms, [round(float(p), 4) for p in res.p_by_m], {"x": "poloidal mode number m", "y": "P(m) — MAC posterior"}, - highlight=best_m, + highlight=chosen_m, secondary={"name": "MAC (nominal)", "y": [round(float(v), 4) for v in res.mac_nominal]}, meta={ - "best_m": best_m, + "best_m": chosen_m, + "raw_m": best_m, "n": n_tor, "verdict": verdict, - "mac_best": round(mac_best, 3), "quality": quality, + "mac_best": round(mac_best, 3), + "margin": round(float(res.margin), 3), + "alias_ms": [int(a) for a in res.alias_ms], + "q_anchored": q_ok, + "q_surface": round(float(q_surface), 3) if q_surface is not None else None, + "q_corrected": bool(q_corrected), + "psi_n": round(float(q_psi_n), 3) if q_psi_n is not None else None, + "allowed_ms": allowed_ms, "p_best": round(float(res.p_best), 3), "p_even": round(float(res.p_even), 3), "p_odd": round(float(res.p_odd), 3), @@ -1074,9 +1153,10 @@ def _poloidal_m_spectrum(shot, params=None) -> dict: "t0_ms": t0_ms, "shot": str(shot), "note": "MAC of the φ-detrended poloidal phase pattern vs e^{−imθ} templates " - f"({n_used}/{n_total} probes above the SNR gate); P(m) from a Monte-Carlo over " - "the per-probe phase σ. 'quality' reflects the absolute MAC — a low value " - "means no pattern matches a pure mode, so the m is not to be trusted.", + f"({n_used}/{n_total} probes above SNR gate); P(m) from a Monte-Carlo over the " + "per-probe phase σ. When an EFIT02 q-profile is present the physical m is the " + "MAC alias with a real q=m/n surface; otherwise the bare MAC winner (alias-prone " + "on a sparse array) with a margin/alias-set caveat.", }, ) diff --git a/tests/test_equilibrium.py b/tests/test_equilibrium.py new file mode 100644 index 0000000..978b0c3 --- /dev/null +++ b/tests/test_equilibrium.py @@ -0,0 +1,81 @@ +"""Tests for magnetics.core.equilibrium — q-profile rational surfaces + m anchoring.""" + +import numpy as np + +from magnetics.core.equilibrium import ( + anchor_poloidal_m, + q_range, + resonant_surfaces, +) + +# a monotonic, physical q-profile: q0 ≈ 1 on axis rising to q_edge ≈ 4.5 +PSI = np.linspace(0.0, 1.0, 129) +QPSI = 1.0 + 3.5 * PSI**2 # q(0)=1.0, q(1)=4.5 + + +class TestQRange: + def test_range(self): + lo, hi = q_range(QPSI) + assert abs(lo - 1.0) < 1e-9 and abs(hi - 4.5) < 1e-9 + + def test_ignores_nan(self): + q = QPSI.copy() + q[10] = np.nan + lo, hi = q_range(q) + assert np.isfinite(lo) and np.isfinite(hi) + + +class TestResonantSurfaces: + def test_n1_surfaces_are_the_integer_q_values(self): + # for n=1, q=m/1=m → surfaces at q=2,3,4 (q=1 is the axis endpoint, q=5>qmax) + surfs = resonant_surfaces(QPSI, PSI, n=1) + ms = [s.m for s in surfs] + assert 2 in ms and 3 in ms and 4 in ms + assert 5 not in ms # q=5 exceeds q_edge=4.5 + + def test_n3_includes_q2_at_m6(self): + # n=3 → q=m/3; q=2 is the m=6 surface (a 2/1-helicity resonance) + surfs = resonant_surfaces(QPSI, PSI, n=3) + by_m = {s.m: s for s in surfs} + assert 6 in by_m and abs(by_m[6].q - 2.0) < 1e-9 + assert 0.0 <= by_m[6].psi_n <= 1.0 # located on the grid + # q=1/3 (m=1) is below q_min=1 → no surface + assert 1 not in by_m + + def test_psi_location_matches_profile(self): + # q=2.0 for n=1 (m=2): q=1+3.5ψ²=2 → ψ=sqrt(1/3.5) + surfs = resonant_surfaces(QPSI, PSI, n=1) + s2 = next(s for s in surfs if s.m == 2) + assert abs(s2.psi_n - np.sqrt(1.0 / 3.5)) < 2e-2 + + def test_zero_n_is_empty(self): + assert resonant_surfaces(QPSI, PSI, n=0) == [] + + +class TestAnchorPoloidalM: + def test_rejects_unphysical_alias_and_picks_physical(self): + # MAC winner m=1 (q=1/3, unphysical for n=3), runner-up m=6 (q=2, physical) + ms = np.array([-1, 6, -2, 4]) + macs = np.array([0.61, 0.575, 0.558, 0.49]) + res = anchor_poloidal_m(ms, macs, n=3, q_psi=QPSI, psi_n=PSI) + assert res.raw_m == -1 # what the bare MAC would have said + assert res.corrected is True + assert abs(res.chosen_m) == 6 and abs(res.q_surface - 2.0) < 1e-9 + assert res.psi_n is not None and 0.0 <= res.psi_n <= 1.0 + assert 6 in res.allowed_ms and 1 not in res.allowed_ms + + def test_keeps_physical_winner_untouched(self): + # winner m=6 already resonates → no correction + ms = np.array([6, -1, 3]) + macs = np.array([0.8, 0.5, 0.4]) + res = anchor_poloidal_m(ms, macs, n=3, q_psi=QPSI, psi_n=PSI) + assert res.corrected is False and res.chosen_m == 6 + + def test_no_physical_alias_leaves_raw(self): + # tiny plasma q∈[3,4]; for n=3 that needs m∈[9,12]; candidates {1,2} don't resonate + q = 3.0 + 1.0 * PSI + ms = np.array([1, 2]) + macs = np.array([0.5, 0.4]) + res = anchor_poloidal_m(ms, macs, n=3, q_psi=q, psi_n=PSI) + assert res.corrected is False and res.chosen_m == 1 # nothing physical to swap to + assert res.allowed_ms == [] diff --git a/tests/test_nodes.py b/tests/test_nodes.py index a08bfa5..2bd8987 100644 --- a/tests/test_nodes.py +++ b/tests/test_nodes.py @@ -220,10 +220,10 @@ def test_poloidal_phase_fit_node(): assert "m_fit" in n["meta"] -def test_poloidal_m_spectrum_node(): - shot = _first_shot() +def test_poloidal_m_spectrum_node(synthetic_shot): + shot = synthetic_shot # the full DIII-D synthetic (poloidal array + EFIT02 q-profile) try: - n = nodes.build_node(shot, "poloidal_m_spectrum", {"time": 3000}) + n = nodes.build_node(shot, "poloidal_m_spectrum", {"time": 50}) except Exception as e: # noqa: BLE001 — shot may lack the poloidal array pytest.skip(f"no poloidal array in this shot: {e}") assert n["kind"] == "bar" @@ -231,13 +231,35 @@ def test_poloidal_m_spectrum_node(): m = n["meta"] assert "best_m" in m and "verdict" in m assert 0.0 <= m["mac_best"] <= 1.0 - assert m["quality"] in {"clean", "marginal", "weak"} + assert m["quality"] in {"clean", "marginal", "weak", "ambiguous", "q-anchored"} assert 0 < m["n_used"] <= m["n_probes"] # SNR gate keeps a subset of the array # P(m) is a proper distribution over the candidates (draws + per-probe σ present) assert m["n_draws"] > 0 assert abs(sum(n["y"]) - 1.0) < 1e-6 assert 0.0 <= m["p_even"] <= 1.0 and 0.0 <= m["p_odd"] <= 1.0 assert abs(m["p_even"] + m["p_odd"] - 1.0) < 1e-6 + # the synthetic fixture ships an EFIT02-style q-profile → the fit is q-anchored to a + # real rational surface (the ground-truth 2/1 emerges even though MAC aliases it) + assert m["q_anchored"] is True + assert m["q_surface"] is not None and m["quality"] == "q-anchored" + assert "resonant q" in m["verdict"] + assert m["margin"] is not None and isinstance(m["alias_ms"], list) + + +def test_poloidal_m_spectrum_falls_back_without_q_profile(synthetic_shot, monkeypatch): + """With no q-profile the fit runs unanchored: quality is never 'q-anchored' and a + near-tie is reported as 'ambiguous' with an alias set, not a false-confident m.""" + shot = synthetic_shot + monkeypatch.setattr(nodes, "_q_at", lambda *a, **k: None) + try: + n = nodes.build_node(shot, "poloidal_m_spectrum", {"time": 50}) + except Exception as e: # noqa: BLE001 + pytest.skip(f"no poloidal array in this shot: {e}") + m = n["meta"] + assert m["q_anchored"] is False and m["q_surface"] is None + assert m["quality"] in {"clean", "marginal", "weak", "ambiguous"} + if m["quality"] == "ambiguous": + assert "∈" in m["verdict"] # reports the alias set, not one confident m assert n["secondary"]["y"] and len(n["secondary"]["y"]) == len(n["x"]) # nominal MAC