Pairwise-difference basis for differential magnetic sensors

CLAUDE.md

@@ -39,7 +39,7 @@ DIII-D shots:
 - **QS live (Day-3 night):** `qs_fit` / `phi_t` / `fit_quality` / `chi_sq_t` / sensor-map / signal
   nodes serve the **real** SLCONTOUR fit. Shots pulled rotating-only (no Bp LFS midplane array)
   return a clean 422 and the QS tab shows a "no quasi-stationary array" banner. Remaining fidelity
-  gap: the data layer's per-sensor σ / helicity (fit currently uses a constant σ).
+  gap: the data layer's per-sensor σ (the fit uses a constant σ; helicity is computed from Ip·Bt).
 - **Geometry shot-indexed (Day-3 night):** `data/device/diiid.json` sensor availability + positions
   are now segmented back to shot 124400 (legacy dense set) / 151593 (3D-upgrade). The Sensors tab
   renders wall + vacuum vessel + perturbation coils + saddle loops (2D honoring each loop's tilt).
@@ -71,13 +71,13 @@ rather than fabricating data in the GUI. Keep `contracts.py` and `contract.ts` i
   label → live via `usingLiveBackend()`; build the Sensors geometry view).
 - **Slow Rollers + Meg (quasi-stationary):** the real `qs_fit` node (K / χ² / modes) + the QS GUI
   tab are **live** (Day-3 night). Remaining: finish the pure `core/quasistationary.py` port and
-  wire it in place of the `_slcontour` reference pipeline (#40); consume real per-sensor σ +
-  helicity once the data layer provides them (the fit currently uses a constant σ).
+  wire it in place of the `_slcontour` reference pipeline (#40); consume real per-sensor σ once the
+  data layer provides it (the fit currently uses a constant σ; helicity is computed from Ip·Bt).
 - **Data Streamers:** the DIII-D geometry table is now **shot-indexed** (`diiid.json` segmented to
   124400 / 151593; the cosmetic θ in `magnetics/data/diiid.py` — now a thin shim over the
   device-agnostic `data/device_geom.py` — is superseded by the real device table). NSTX/NSTX-U
   fetch + node rendering landed on `feature/nstxu-data-fetch` (see the NSTX status bullet above).
-  Remaining: give `h5source` per-sensor σ + helicity (last QS-fidelity gap); a `DataSource`
+  Remaining: give `h5source` per-sensor σ (last QS-fidelity gap; helicity now computed from Ip·Bt); a `DataSource`
   abstraction with an array cache; populate the shot-segmented legacy-NSTX tree/wall from the
   `config_hf/hn.mm` files; import other-device geometry the same way.
 - **Structural cleanup (LANDED — PR #41, Day-2 night):** the project was hoisted to the repo root

docs/specs/quasistationary-mode-analysis.md

@@ -52,7 +52,8 @@ include `device`, `shot`, `tmin`/`tmax` (ms).
   per-channel geometry joined from the device JSON, including the derived `*_end1/2` coordinates. All
   sensor channels are interpolated onto one common time axis (the densest native grid, clipped to the
   window they share).
-- `plasma` — `Ip`, `Bt` from the `ip`/`bt` channels (time in **ms**), `helicity` attr (default −1).
+- `plasma` — `Ip`, `Bt` from the `ip`/`bt` channels (time in **ms**), and a `helicity` attr computed
+  as sign(Bt·Ip) from those traces (−1 fallback when they're absent).
 - `coupling` — `None`; the new files carry no DC vacuum-coupling matrix, so `qs_prep`'s DC
   compensation (`dc_comp=True`) is unavailable.
 
@@ -108,8 +109,8 @@ flips helicity. A regression test in `tests/test_qs_bridge.py` pins this.
 ## Notes & limitations
 
 - Reduced χ² typically runs above 1 because the constant 2e-5 T sensor σ is optimistic relative to
-  higher-n structure not in the basis; residuals stay small versus the signals. Per-sensor σ /
-  helicity from the data layer is an open item.
+  higher-n structure not in the basis; residuals stay small versus the signals. Per-sensor σ from the
+  data layer is an open item (helicity is now computed from the shot's Ip·Bt in `qs_io_data`).
 - The pure-numpy relatives of `qs_fit` (a device-agnostic `core` port) and the `data/`-layer
   set-flattening dedup are open follow-ups.
 - Out of scope for this path: the rotating-mode spectrogram/MODESPEC analysis, 3D coil-current fits,

src/magnetics/core/qs_device.py

@@ -32,22 +32,16 @@
 from ..data import devices as _devices
 from ..data import diiid_geometry as _diiid_geo
 
-# ``data.devices.load_device`` resolves ``<name>.json`` by lowercasing the name; a few
-# display names don't slugify by lowercasing alone (``"DIII-D"`` → ``diiid.json``, not
-# ``diii-d.json``). Preserve the alias the old shim carried.
-_DEVICE_ALIAS = {"DIII-D": "diiid"}
-
-
-def _device_key(device: str) -> str:
-    """Map a device display name to its ``data/device/<key>.json`` stem."""
-    if device in _DEVICE_ALIAS:
-        return _DEVICE_ALIAS[device]
-    return re.sub(r"[^a-z0-9]", "", str(device).lower())
-
 
 def load_device(device: str = "DIII-D") -> dict:
-    """The parsed device JSON dict, via the data layer (with QS name aliasing)."""
-    return _devices.load_device(_device_key(device))
+    """The parsed device JSON dict, via the data layer's name/id resolver.
+
+    ``devices.resolve_device_id`` maps a display name (e.g. ``"DIII-D"``) or a config
+    id to the JSON file stem using each device file's ``name`` field — so no device
+    alias is hardcoded here. Falls back to a slugified name if it can't be resolved.
+    """
+    dev_id = _devices.resolve_device_id(device) or re.sub(r"[^a-z0-9]", "", str(device).lower())
+    return _devices.load_device(dev_id)
 
 
 # ── device-name comparison ─────────────────────────────────────────────────────
@@ -132,7 +126,9 @@ def sensor_geometry(device="DIII-D", shot=None):
     Base fields (``r, z, phi, tilt, length, delta_phi, na``, plus ``pair``) come from
     the data layer's shot-aware resolver (:func:`devices.geometry_nearest`); this
     function adds the derived poloidal angle ``theta`` and the sensor-end coordinates
-    ``{r,z,phi,theta}_end1/2`` that the SLCONTOUR fit needs. A sensor is modelled as a
+    ``{r,z,phi,theta}_end1/2`` that the SLCONTOUR fit needs, plus the paired sensor's
+    ``pair_{phi,theta,z}_end1/2`` (for the pairwise-difference basis; NaN when unpaired).
+    A sensor is modelled as a
     straight segment of physical ``length`` centred at ``(r, z)`` tilted by ``tilt``
     (degrees) in the poloidal plane and spanning ``delta_phi`` (degrees) toroidally;
     angles are measured about the magnetic axis at ``(R0, 0)``.
@@ -169,6 +165,20 @@ def sensor_geometry(device="DIII-D", shot=None):
     theta_end1 = np.degrees(np.arctan2(z_end1, r_end1 - R0))
     theta_end2 = np.degrees(np.arctan2(z_end2, r_end2 - R0))
 
+    # ── paired-sensor ends (for the SLCONTOUR pairwise-difference basis) ──────────
+    # A differential sensor (pair != "None") reports field(X) - field(pair); the fit
+    # differences the basis at X and at its pair. Resolve each channel's pair to its
+    # row and gather the pair's *_end1/2 here (NaN when unpaired/unresolved) so the
+    # pair geometry rides the normal geometry-attach path into the fit.
+    name_to_row = {c: i for i, c in enumerate(channels)}
+    pair_row = np.array([name_to_row.get(str(p), -1) if str(p) != "None" else -1 for p in pairs])
+
+    def _by_pair(arr):
+        out = np.full(len(channels), np.nan)
+        valid = pair_row >= 0
+        out[valid] = arr[pair_row[valid]]
+        return out
+
     ds = xr.Dataset(
         {
             **{k: ("channel", cols[k]) for k in base},
@@ -181,6 +191,12 @@ def sensor_geometry(device="DIII-D", shot=None):
             "phi_end2": ("channel", phi_end2),
             "theta_end1": ("channel", theta_end1),
             "theta_end2": ("channel", theta_end2),
+            "pair_phi_end1": ("channel", _by_pair(phi_end1)),
+            "pair_phi_end2": ("channel", _by_pair(phi_end2)),
+            "pair_theta_end1": ("channel", _by_pair(theta_end1)),
+            "pair_theta_end2": ("channel", _by_pair(theta_end2)),
+            "pair_z_end1": ("channel", _by_pair(z_end1)),
+            "pair_z_end2": ("channel", _by_pair(z_end2)),
         },
         coords={"channel": channels},
     )

src/magnetics/core/qs_fit.py

@@ -33,8 +33,6 @@
 import numpy as np
 import xarray as xr
 
-from .qs_device import is_device
-
 logger = logging.getLogger(__name__)
 
 
@@ -251,16 +249,32 @@ def _printv(*a):
                 fill = np.nanmean(sigma) if np.isfinite(np.nanmean(sigma)) else 1.0
             sigma[bad] = fill
 
+    # ── paired-sensor differencing (SLCONTOUR pairwise-difference basis) ────────
+    # A differential sensor (pair != "None") reports field(X) - field(pair), so the
+    # design-matrix column is differenced at X and at its pair. The pair ends ride in
+    # on the geometry (qs_device.sensor_geometry); NaN marks an unpaired channel.
+    if f"pair_{xkey}_end1" in ds:
+        px1 = ds[f"pair_{xkey}_end1"].values
+        px2 = ds[f"pair_{xkey}_end2"].values
+        py1 = ds[f"pair_{ykey}_end1"].values
+        py2 = ds[f"pair_{ykey}_end2"].values
+    else:
+        px1 = px2 = py1 = py2 = np.full_like(x1, np.nan)
+    has_pair = np.isfinite(px1)
+    # OMFIT sources the pair's signal_sigma (fit_magnetics.py:136-141: all-NaN -> own
+    # sigma, then NaN -> nanmean). Our sigma is the constant per-shot value for every
+    # channel, so pair sigma == own sigma; revisit when real per-sensor sigma lands.
+    psigma = sigma
+
     # ── basis-specific mode / centre setup ────────────────────────────────────
     is_sinusoidal = fit_basis.startswith("sinusoidal")
     is_gaussian = fit_basis.startswith("gaussian")
 
     if is_sinusoidal:
         ms = np.atleast_1d(ms)
         ns = np.atleast_1d(ns)
-        dp = None  # paired sensors not present in this dataset
-        if dp is None and 0 not in ns:
-            _printv("WARNING: Sensors are not paired! Consider including n=0.")
+        if (~has_pair).any() and 0 not in ns:
+            _printv("WARNING: Some sensors are not paired! Consider including n=0.")
         if not np.all(ms == 0) and not np.any(np.sign(ms) == helicity):
             ms = ms * -1
             _printv(f"WARNING: Flipping sign of m to conform to helicity {helicity:+}")
@@ -271,11 +285,6 @@ def _printv(*a):
         _ncycle = _mcycle = 0
         _nepsilon = _mepsilon = np.inf
 
-        if is_device(ds.attrs.get("device", ""), "DIII-D"):
-            if any(re.match(k, c) for k in ("C.*", "IL.*", "IU.*") for c in channels):
-                if fit_basis != "sinusoidal-point":
-                    logger.warning("sinusoidal-point basis is used by DIII-D 3D coil operators")
-
     elif is_gaussian:
         # Derive RBF centres from sensor extent (mirrors OMFIT gaussian setup block)
         x_all = np.hstack((x1, x2))
@@ -325,6 +334,9 @@ def _printv(*a):
     if is_sinusoidal:
         for n, m in nms:
             fmn = form_basis_function(n, m, x1, x2, y1, y2, fit_basis) / sigma
+            if has_pair.any():
+                fmn_pair = form_basis_function(n, m, px1, px2, py1, py2, fit_basis) / psigma
+                fmn = np.where(has_pair, (fmn - fmn_pair) / 2.0, fmn)
             if np.allclose(fmn.imag, 0):
                 A_cols.append(fmn.real)
                 ncomp.append(1)
@@ -345,6 +357,14 @@ def _printv(*a):
                         form_basis_function(n, m, x1 + x0, x2 + x0, y1 + y0, y2 + y0, fit_basis)
                         / sigma
                     )
+                    if has_pair.any():
+                        fmn_pair = (
+                            form_basis_function(
+                                n, m, px1 + x0, px2 + x0, py1 + y0, py2 + y0, fit_basis
+                            )
+                            / psigma
+                        )
+                        fmn = np.where(has_pair, (fmn - fmn_pair) / 2.0, fmn)
                     A_cols = A_cols[:-2] + [fmn.real]
                     ncomp[-1] = 1
     else:  # gaussian — basis functions are always real
@@ -365,6 +385,24 @@ def _printv(*a):
                 )
                 / sigma
             )
+            if has_pair.any():
+                fmn_pair = (
+                    form_basis_function(
+                        n,
+                        m,
+                        px1,
+                        px2,
+                        py1,
+                        py2,
+                        fit_basis,
+                        ncycle=_ncycle,
+                        mcycle=_mcycle,
+                        nepsilon=_nepsilon,
+                        mepsilon=_mepsilon,
+                    )
+                    / psigma
+                )
+                fmn = np.where(has_pair, (fmn - fmn_pair) / 2.0, fmn)
             A_cols.append(fmn)
             ncomp.append(1)
 
@@ -448,6 +486,7 @@ def _printv(*a):
         fit_mcycle=_mcycle,
         fit_neps=_nepsilon,
         fit_meps=_mepsilon,
+        n_paired=int(has_pair.sum()),
     )
     _condition_warning(raw_cn)
     return ds

src/magnetics/core/qs_io_data.py

@@ -82,8 +82,10 @@ def load_shot(shot, data_root=DATAFILE_ROOT, helicity=-1):
 
     :param shot: shot number (int/str) or a path to a ``shot_<n>.h5`` file.
     :param data_root: directory holding the ``shot_<n>.h5`` files.
-    :param helicity: field/current helicity convention used to orient mode signs
-        (default ``-1``; the new data files do not store it).
+    :param helicity: **fallback** helicity only. The real value (sign of Bt·Ip, used
+        to orient the poloidal-mode sign) is computed from the shot's own ``ip``/``bt``
+        traces; the default ``-1`` is used only when those are missing (the data files
+        do not store helicity directly).
     :return: :class:`ShotData` (``coupling`` is ``None`` — the new files carry no
         DC vacuum-coupling matrix).
     """
@@ -160,11 +162,27 @@ def _build_raw(names, sigs, times, geo, shot_no, device):
     return raw
 
 
+def _compute_helicity(sigs, fallback=-1):
+    """Field/current helicity = sign(median Bt) / sign(median Ip) (±1), from the shot's
+    own ``bt``/``ip`` traces (matches OMFIT ``init_magnetics``) — it orients the fitted
+    poloidal-mode sign. Returns ``fallback`` when either trace is absent or its median
+    sign is 0 (indeterminate)."""
+    if "bt" in sigs and "ip" in sigs:
+        sbt = np.sign(np.nanmedian(np.asarray(sigs["bt"], dtype=float)))
+        sip = np.sign(np.nanmedian(np.asarray(sigs["ip"], dtype=float)))
+        if sbt != 0 and sip != 0:
+            return int(sbt / sip)
+    return int(fallback)
+
+
 def _build_plasma(sigs, times, helicity):
-    """Assemble the ``plasma`` Dataset (Ip/Bt on a shared ms time base)."""
+    """Assemble the ``plasma`` Dataset (Ip/Bt on a shared ms time base).
+
+    The ``helicity`` attr is computed from the ``ip``/``bt`` medians
+    (:func:`_compute_helicity`); the ``helicity`` argument is only the fallback."""
     if not sigs:
         plasma = xr.Dataset(coords={"time": np.array([], dtype=float)})
-        plasma.attrs["helicity"] = int(helicity)
+        plasma.attrs["helicity"] = _compute_helicity(sigs, helicity)
         return plasma
 
     base = "ip" if "ip" in times else next(iter(times))
@@ -181,7 +199,7 @@ def _build_plasma(sigs, times, helicity):
         data_vars[var] = ("time", y)
 
     plasma = xr.Dataset(data_vars, coords={"time": t_ms})
-    plasma.attrs["helicity"] = int(helicity)
+    plasma.attrs["helicity"] = _compute_helicity(sigs, helicity)
     return plasma
 
 

tests/synthetic_shot.py

@@ -99,13 +99,27 @@ def _channels(shot: int, *, include_qs: bool) -> list[Channel]:
         names += _MPID_POLOIDAL
     names = list(dict.fromkeys(names))  # dedup, keep order
 
+    from magnetics.data import devices
+
+    dev = devices.load_device("diiid")
+
     chans: list[Channel] = []
     for nm in names:
         phi = diiid.phi_of(nm, shot)
         if phi is None:
             continue  # not modeled at this shot → leave out (fetcher would skip it)
         theta = diiid.real_theta_of(nm, shot) or 0.0
-        chans.append(Channel(nm, t_ms.copy(), _sensor_signal(t_s, phi, theta), ok=True))
+        sig = _sensor_signal(t_s, phi, theta)
+        # Differential probes (MPID/ISLD/ESLD) physically report field(X) − field(pair);
+        # subtract the pair's field so the signal matches the pairwise-difference basis
+        # the QS fit builds for paired sensors (see qs_fit / qs_device).
+        pair = (devices.geometry_nearest(dev, nm, shot) or {}).get("pair")
+        if pair and str(pair) != "None":
+            pphi = diiid.phi_of(pair, shot)
+            if pphi is not None:
+                ptheta = diiid.real_theta_of(pair, shot) or 0.0
+                sig = sig - _sensor_signal(t_s, pphi, ptheta)
+        chans.append(Channel(nm, t_ms.copy(), sig, ok=True))
 
     # Plasma traces (routed to the plasma Dataset / κ path; not sensor channels).
     ip = (1.0e6 * t_s / t_s[-1]).astype(np.float32)  # 0 → 1 MA ramp