Bugfix/bandpass fix

docs/specs/quasistationary-mode-analysis.md

@@ -0,0 +1,116 @@
+# Quasi-stationary mode analysis (SLCONTOUR-style) — pipeline details
+
+The **quasi-stationary** path analyses locked and slowly-rotating MHD modes by fitting the *spatial*
+pattern of the perturbed field measured by a sensor **array** (VISION.md §4.1). At each time slice it
+decomposes δB(φ, θ) into a small set of toroidal/poloidal harmonics by SVD-conditioned least squares,
+and reports the design-matrix **condition number K** (the central trust metric — warn K > 10, error
+K > 20), per-coefficient error bars, and reduced χ².
+
+It runs entirely off the project-canonical data — per-shot signals in `data/datafile/shot_<n>.h5`
+and device geometry in `data/device/<device>.json` — with no MDSplus or external framework.
+
+The fit basis (cylindrical geometry) is
+
+```
+δB(φ, θ) = Σ_nm  b_nm · exp(i(nφ + mθ))
+```
+
+> **Scope.** This is the spatial fit of sensor *arrays*. The rotating-mode spectrogram / MODESPEC
+> path (`magnetics.core.spectral`) and 3D coil-current fits are separate and out of scope here.
+
+## The pipeline (`magnetics.core.qs_*`)
+
+All modules live in `src/magnetics/core/` under the unified `qs_` prefix:
+
+| module | role |
+|---|---|
+| `qs_io_data.load_shot` | read `shot_<n>.h5` signals and join the per-channel device geometry → the `raw` / `plasma` Datasets |
+| `qs_prep.prepare` | trim (channels + time), optional integrate (bdot→B), causal band/high/low-pass, detrend, SVD-condition the data matrix |
+| `qs_fit.fit` | build the design matrix from the basis, SVD it (K + error bars), least-squares fit every time slice → the `fit` Dataset |
+| `qs_run.run_steps` | the `load → prep → fit` orchestrator; returns a `MagneticsRun` bundling `raw` / `prepared` / `plasma` / `fit` (+ `condition_number`) |
+| `qs_device` | device-JSON readers — `sensor_geometry` (base + derived `theta`/`*_end1/2`), `resolve_channel_filter`, `list_sensor_subsets`, `load_wall` — **delegating to the `data/` layer** (`data.devices`, `data.diiid_geometry`) |
+| `qs_bridge` | adapt the `fit` Dataset → GUI `kind`-nodes — the **production/service** output path |
+| `qs_plots` | standalone **matplotlib** plots (sensor map, signals, SVD diagnostics, fit quality, mode amp/phase, φ-vs-time contour) — for notebooks/offline use, **not** wired into the GUI |
+
+There are two consumers of the `fit` Dataset. The **service** builds JSON nodes from it via
+`qs_bridge` (served by `service/nodes.py` → the GUI's QS tab). The **notebook/offline** path renders
+it directly with `qs_plots`. They share the same physics and the same reconstruction sign convention
+(below); `qs_bridge` does not import `qs_plots`.
+
+## Data inputs
+
+### Per-shot signals — `data/datafile/shot_<n>.h5`
+
+One HDF5 file per shot (the PTDATA fetch output): a group per channel with `data` + `time`
+(hard-linked into a shared `_timebases` group), **all time in milliseconds**. Channels sample at
+several rates (integrated probes ~50 kHz, bdots ~200 kHz, coils/`ip`/`bt` ~20 kHz). Root attrs
+include `device`, `shot`, `tmin`/`tmax` (ms).
+
+`qs_io_data.load_shot` turns this into the Datasets the pipeline expects:
+
+- `raw` — `signal(channel, time)` (time in **seconds**), `signal_sigma` (a constant 2e-5 T), and the
+  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).
+- `coupling` — `None`; the new files carry no DC vacuum-coupling matrix, so `qs_prep`'s DC
+  compensation (`dc_comp=True`) is unavailable.
+
+### Device geometry — `data/device/<device>.json`
+
+The canonical device description (`diiid.json` for DIII-D): `R0`, the shot-segmented `first_wall`
+outline, the per-sensor base geometry (`r, z, phi, tilt, length, delta_phi, na, pair`), and named
+`sensor_sets`. `qs_device` reads it through the shared `data.devices` resolvers (so the QS pipeline
+and the fetcher never disagree about a shot's geometry) and adds the QS-specific derived `theta` and
+`*_end1/2` sensor-end coordinates in `sensor_geometry`.
+
+## Quick start
+
+The project is a `uv` package rooted at the repo; import `magnetics.core.qs_*` directly (no
+`sys.path` juggling):
+
+```python
+from magnetics.core.qs_run import run_steps
+from magnetics.core import qs_plots as plots
+
+# load → prep → fit the LFS-midplane toroidal Bp array for n = 1,2,3
+r = run_steps(199749, channel_filter="Bp_LFS_midplane", ns=(1, 2, 3), ms=(0,),
+              time_trim=(3.3, 3.5), prep_kwargs=dict(cutoff_hz=(5, 250), energy=0.98))
+
+print("condition number K =", r.condition_number)
+plots.plot_fit_modes(r.fit)                 # amplitude & phase of each mode vs time
+plots.plot_slice(r.fit, fix_coord="theta")  # the SLCONTOUR φ-vs-time contour
+```
+
+`channel_filter` accepts a regex, a list of regexes, or a friendly subset name — with or without
+underscores (e.g. `"Bp_LFS_midplane"` or `"Bp LFS midplane"`, `"Bp_All"`, `"All_3D_Coils"`). List
+them all with `qs_io_data.available_subsets("DIII-D")`. `time_trim` must fall inside the shot file's
+window (shot 199749 spans **3.3–3.5 s**).
+
+The worked, shot-configurable notebook is [`examples/example_magnetics.ipynb`](../examples/example_magnetics.ipynb)
+(needs a fetched `shot_<n>.h5`; the repo ships no tokamak data).
+
+## Reconstruction sign convention
+
+The fit basis is `exp(+i(nφ + mθ))`, but each complex basis column is split into two *real*
+least-squares columns and the complex coefficient is reassembled as `b = x_r + i·x_i`. Because of
+that split, the field that reproduces the fitted signal is
+
+```
+δB(φ, θ) = Re Σ_nm  b_nm · exp(−i(nφ + mθ))
+```
+
+i.e. **reconstruction uses `exp(−i(…))`** even though the fit basis is `+i`. This is honoured by
+`qs_bridge._reconstruct_grid` / `fit_to_phi_t_node` and `qs_plots.plot_slice` (and matches the DIII-D
+SLCONTOUR reference). Do not "align" the reconstruction to `+i` — that mirrors the toroidal phase and
+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.
+- 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,
+  and internal/external source separation (needs Br, which a Bp-only array lacks).

src/magnetics/_slcontour/example_magnetics.ipynb → examples/example_magnetics.ipynb

@@ -4,51 +4,15 @@
    "cell_type": "markdown",
    "id": "0",
    "metadata": {},
-   "source": [
-    "# 3D magnetics spatial fit — DIII-D (configurable shot)\n",
-    "\n",
-    "A worked example of the **SLCONTOUR-style quasi-stationary spatial fit** (VISION.md §4.1),\n",
-    "run **locally** — off OMFIT and off MDSplus — using the ported OMFIT magnetics scripts in this\n",
-    "directory.  It reads the per-shot signals from `data/datafile/shot_<shot>.h5` and the device\n",
-    "geometry from `data/device/diiid.json`.\n",
-    "\n",
-    "**Pipeline:** `load → prep → fit → plot`, the local translation of the OMFIT\n",
-    "`fetch → prep → fit → plot` workflow:\n",
-    "\n",
-    "| local module | OMFIT script |\n",
-    "|---|---|\n",
-    "| `io_data.load_shot` | `SCRIPTS/fetch_magnetics.py` + `init_magnetics.py` (signals + derived geometry) |\n",
-    "| `prep.prepare` | `SCRIPTS/prep_magnetics.py` |\n",
-    "| `fit.fit` | `SCRIPTS/fit_magnetics.py` |\n",
-    "| `run.run_steps` | `SCRIPTS/run_magnetics.py` |\n",
-    "| `plots.*` | `PLOTS/plot_magnetics_*.py` |\n",
-    "\n",
-    "The run is driven entirely by the **Parameters** cell below — change `SHOT`, `CHANNEL_FILTER`,\n",
-    "`TIME_TRIM`, or the mode lists to analyse a different shot or array.  The default analyses the\n",
-    "LFS-midplane toroidal Bp array (`Bp_LFS_midplane`) of shot **199749** for toroidal mode numbers\n",
-    "n = 1, 2, 3."
-   ]
+   "source": "# 3D magnetics spatial fit — DIII-D (configurable shot)\n\nA worked example of the **SLCONTOUR-style quasi-stationary spatial fit** (VISION.md §4.1),\nrun **locally** — off MDSplus — with the `magnetics.core.qs_*` pipeline. It reads the per-shot\nsignals from `data/datafile/shot_<shot>.h5` and the device geometry from `data/device/diiid.json`.\n\n**Pipeline:** `load → prep → fit → plot`:\n\n| module | role |\n|---|---|\n| `qs_io_data.load_shot` | read `shot_<n>.h5` signals + join the device-JSON geometry (derived `*_end1/2` ends) |\n| `qs_prep.prepare` | trim, integrate, causal band-pass, detrend, SVD-condition |\n| `qs_fit.fit` | build the basis matrix + SVD-conditioned least-squares modal fit |\n| `qs_run.run_steps` | the load → prep → fit orchestrator |\n| `qs_plots.*` | sensor map, signals, SVD diagnostics, fit quality, mode amp/phase, φ-vs-time contour |\n\n(The production/GUI path adapts the same `qs_fit` Dataset through `qs_bridge`; `qs_plots` here is\nthe standalone-matplotlib view.)\n\nThe run is driven entirely by the **Parameters** cell below — change `SHOT`, `CHANNEL_FILTER`,\n`TIME_TRIM`, or the mode lists to analyse a different shot or array.  The default analyses the\nLFS-midplane toroidal Bp array (`Bp_LFS_midplane`) of shot **199749** for toroidal mode numbers\nn = 1, 2, 3."
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "1",
    "metadata": {},
    "outputs": [],
-   "source": [
-    "import sys\n",
-    "sys.path.insert(0, '.')  # make the local modules importable\n",
-    "\n",
-    "import numpy as np\n",
-    "import matplotlib.pyplot as plt\n",
-    "\n",
-    "from io_data import load_shot, valid_channels, available_subsets\n",
-    "from omfit_compat import resolve_channel_filter\n",
-    "from run import run_steps\n",
-    "import plots\n",
-    "\n",
-    "%matplotlib inline"
-   ]
+   "source": "import numpy as np\nimport matplotlib.pyplot as plt\n\nfrom magnetics.core.qs_io_data import load_shot, valid_channels, available_subsets\nfrom magnetics.core.qs_device import resolve_channel_filter\nfrom magnetics.core.qs_run import run_steps\nfrom magnetics.core import qs_plots as plots\n\n%matplotlib inline"
   },
   {
    "cell_type": "markdown",
@@ -69,33 +33,26 @@
    "outputs": [],
    "source": [
     "# ── Run parameters — change these to analyse a different shot/array ──\n",
-    "SHOT           = 199749                 # data/datafile/shot_<SHOT>.h5\n",
-    "CHANNEL_FILTER = 'Bp_LFS_midplane'      # friendly name (see available_subsets) or a regex\n",
-    "TIME_TRIM      = (3.3, 3.5)             # seconds; must lie inside the shot file's window\n",
-    "NS, MS         = (1, 2, 3), (0,)        # toroidal / poloidal mode numbers\n",
+    "SHOT = 199749  # data/datafile/shot_<SHOT>.h5\n",
+    "CHANNEL_FILTER = \"Bp_LFS_midplane\"  # friendly name (see available_subsets) or a regex\n",
+    "TIME_TRIM = (3.3, 3.5)  # seconds; must lie inside the shot file's window\n",
+    "NS, MS = (1, 2, 3), (0,)  # toroidal / poloidal mode numbers\n",
     "\n",
     "# detrend: 'none' | 'baseline' | 'linear' | 'endpoints'.  DETREND_BAND is the\n",
     "# sub-window used to estimate the trend — here the first 0.01 s of TIME_TRIM.\n",
-    "DETREND_TYPE   = 'baseline'\n",
-    "DETREND_BAND   = (TIME_TRIM[0], TIME_TRIM[0] + 0.01)\n",
+    "DETREND_TYPE = \"baseline\"\n",
+    "DETREND_BAND = (TIME_TRIM[0], TIME_TRIM[0] + 0.01)\n",
     "\n",
-    "PREP_KWARGS    = dict(cutoff_hz=(5, 250), energy=0.98,\n",
-    "                      detrend_type=DETREND_TYPE, detrend_band=DETREND_BAND)"
+    "PREP_KWARGS = dict(\n",
+    "    cutoff_hz=(5, 250), energy=0.98, detrend_type=DETREND_TYPE, detrend_band=DETREND_BAND\n",
+    ")"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "4",
    "metadata": {},
-   "source": [
-    "## 1. Load the shot (the \"fetch\" + \"init\" step)\n",
-    "\n",
-    "`load_shot` reads the raw sensor signals from `data/datafile/shot_<SHOT>.h5` and joins the\n",
-    "per-channel geometry — including the derived sensor-end coordinates (`*_end1/2`) — from\n",
-    "`data/device/diiid.json`.  All sensor channels are interpolated onto a single common time axis\n",
-    "(in seconds); the global `ip`/`bt` traces become the `plasma` Dataset.  The new files carry no\n",
-    "DC-coupling matrix, so `coupling` is `None`."
-   ]
+   "source": "## 1. Load the shot\n\n`load_shot` reads the raw sensor signals from `data/datafile/shot_<SHOT>.h5` and joins the\nper-channel geometry — including the derived sensor-end coordinates (`*_end1/2`) — from\n`data/device/diiid.json`.  All sensor channels are interpolated onto a single common time axis\n(in seconds); the global `ip`/`bt` traces become the `plasma` Dataset.  The new files carry no\nDC-coupling matrix, so `coupling` is `None`."
   },
   {
    "cell_type": "code",
@@ -105,19 +62,19 @@
    "outputs": [],
    "source": [
     "sd = load_shot(SHOT)\n",
-    "print(f'shot {sd.shot}  device {sd.device}')\n",
-    "print('RAW channels:', sd.raw.sizes['channel'], ' time samples:', sd.raw.sizes['time'])\n",
-    "print(f'file window: {float(sd.raw.time[0]):.4f}–{float(sd.raw.time[-1]):.4f} s')\n",
-    "print('helicity:', sd.plasma.attrs['helicity'])\n",
+    "print(f\"shot {sd.shot}  device {sd.device}\")\n",
+    "print(\"RAW channels:\", sd.raw.sizes[\"channel\"], \" time samples:\", sd.raw.sizes[\"time\"])\n",
+    "print(f\"file window: {float(sd.raw.time[0]):.4f}–{float(sd.raw.time[-1]):.4f} s\")\n",
+    "print(\"helicity:\", sd.plasma.attrs[\"helicity\"])\n",
     "\n",
     "# channels in the selected array that carry good data\n",
     "good = valid_channels(sd.raw, CHANNEL_FILTER, sd.device)\n",
-    "print(f'\\n{CHANNEL_FILTER}: {len(good)} valid sensors  -> resolves |n| <= {(len(good)-1)//2}')\n",
+    "print(f\"\\n{CHANNEL_FILTER}: {len(good)} valid sensors  -> resolves |n| <= {(len(good) - 1) // 2}\")\n",
     "print(good)\n",
     "\n",
     "# the plot helpers match channel names with re.match, so turn the friendly\n",
     "# filter into an explicit regex of its sensor names\n",
-    "ARRAY_REGEX = '|'.join(resolve_channel_filter(CHANNEL_FILTER, sd.device))"
+    "ARRAY_REGEX = \"|\".join(resolve_channel_filter(CHANNEL_FILTER, sd.device))"
    ]
   },
   {
@@ -139,7 +96,7 @@
    "outputs": [],
    "source": [
     "for name, sensors in available_subsets(sd.device).items():\n",
-    "    print(f'{name:24s} {sensors}')"
+    "    print(f\"{name:24s} {sensors}\")"
    ]
   },
   {
@@ -160,10 +117,10 @@
    "outputs": [],
    "source": [
     "fig, ax = plt.subplots(1, 2, figsize=(12, 5))\n",
-    "plots.plot_sensors(sd.raw, ARRAY_REGEX, geometry='rz', ax=ax[0])\n",
-    "ax[0].set_title('Cross-section (R–z)')\n",
-    "plots.plot_sensors(sd.raw, ARRAY_REGEX, geometry='cylindrical', ax=ax[1])\n",
-    "ax[1].set_title('Unrolled (phi–theta)')\n",
+    "plots.plot_sensors(sd.raw, ARRAY_REGEX, geometry=\"rz\", ax=ax[0])\n",
+    "ax[0].set_title(\"Cross-section (R–z)\")\n",
+    "plots.plot_sensors(sd.raw, ARRAY_REGEX, geometry=\"cylindrical\", ax=ax[1])\n",
+    "ax[1].set_title(\"Unrolled (phi–theta)\")\n",
     "plt.tight_layout()"
    ]
   },
@@ -187,20 +144,22 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sets = ['Bp_LFS_midplane', 'Bp_LFS_R+1', 'Bp_LFS_R-1', 'Bp_LFS_R+2', 'Bp_LFS_R-2', 'Bp_HFS_Sensors']\n",
+    "sets = [\"Bp_LFS_midplane\", \"Bp_LFS_R+1\", \"Bp_LFS_R-1\", \"Bp_LFS_R+2\", \"Bp_LFS_R-2\", \"Bp_HFS_Sensors\"]\n",
     "colors = plt.cm.tab10.colors\n",
     "\n",
     "fig, ax = plt.subplots(1, 2, figsize=(13, 5))\n",
     "handles = []\n",
     "for i, name in enumerate(sets):\n",
-    "    regex = '|'.join(resolve_channel_filter(name, sd.device))  # friendly name -> regex of its sensors\n",
-    "    plots.plot_sensors(sd.raw, regex, geometry='cylindrical', ax=ax[0], color=colors[i])\n",
-    "    plots.plot_sensors(sd.raw, regex, geometry='rz',          ax=ax[1], color=colors[i])\n",
+    "    regex = \"|\".join(\n",
+    "        resolve_channel_filter(name, sd.device)\n",
+    "    )  # friendly name -> regex of its sensors\n",
+    "    plots.plot_sensors(sd.raw, regex, geometry=\"cylindrical\", ax=ax[0], color=colors[i])\n",
+    "    plots.plot_sensors(sd.raw, regex, geometry=\"rz\", ax=ax[1], color=colors[i])\n",
     "    handles.append(plt.Line2D([], [], color=colors[i], label=name))\n",
     "\n",
-    "ax[0].set_title('Unrolled (phi–theta)')\n",
-    "ax[1].set_title('Cross-section (R–z)')\n",
-    "ax[0].legend(handles=handles, fontsize=8, loc='upper right')\n",
+    "ax[0].set_title(\"Unrolled (phi–theta)\")\n",
+    "ax[1].set_title(\"Cross-section (R–z)\")\n",
+    "ax[0].legend(handles=handles, fontsize=8, loc=\"upper right\")\n",
     "plt.tight_layout()"
    ]
   },
@@ -228,14 +187,17 @@
     "r = run_steps(\n",
     "    SHOT,\n",
     "    channel_filter=CHANNEL_FILTER,\n",
-    "    ns=NS, ms=MS,\n",
+    "    ns=NS,\n",
+    "    ms=MS,\n",
     "    time_trim=TIME_TRIM,\n",
     "    prep_kwargs=PREP_KWARGS,\n",
     ")\n",
-    "print(f'\\ncondition number K = {r.condition_number:.2f}')\n",
-    "print(f'mean reduced chi^2 = {float(np.nanmean(r.fit[\"red_chi_sq\"].values)):.1f}')\n",
-    "print(f'data-matrix SVD effective rank (98% energy) = {r.fit.attrs[\"signal_effective_rank\"]} '\n",
-    "      f'of {r.fit.sizes[\"channel\"]} sensors')"
+    "print(f\"\\ncondition number K = {r.condition_number:.2f}\")\n",
+    "print(f\"mean reduced chi^2 = {float(np.nanmean(r.fit['red_chi_sq'].values)):.1f}\")\n",
+    "print(\n",
+    "    f\"data-matrix SVD effective rank (98% energy) = {r.fit.attrs['signal_effective_rank']} \"\n",
+    "    f\"of {r.fit.sizes['channel']} sensors\"\n",
+    ")"
    ]
   },
   {
@@ -347,23 +309,14 @@
    "outputs": [],
    "source": [
     "# fix_value is the poloidal angle of the slice; the LFS-midplane array sits at theta ~ 0\n",
-    "plots.plot_slice(r.fit, fix_coord='theta', fix_value=0.0)"
+    "plots.plot_slice(r.fit, fix_coord=\"theta\", fix_value=0.0)"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "24",
    "metadata": {},
-   "source": [
-    "---\n",
-    "### Recap\n",
-    "\n",
-    "We reproduced the VISION §4.1 SLCONTOUR outputs locally from the ported OMFIT scripts, driven\n",
-    "off the `data/datafile/shot_<SHOT>.h5` signals and `data/device/diiid.json` geometry: sensor map,\n",
-    "conditioned signals, data/design-matrix SVD diagnostics, fit quality, **amplitude & phase of each\n",
-    "(n, m) mode vs time**, and the **phi-vs-time contour**. To analyse a different shot or array, edit\n",
-    "the **Parameters** cell (`SHOT`, `CHANNEL_FILTER`, `TIME_TRIM`, `NS`/`MS`)."
-   ]
+   "source": "---\n### Recap\n\nWe reproduced the VISION §4.1 SLCONTOUR outputs locally with the `magnetics.core.qs_*` pipeline,\ndriven off the `data/datafile/shot_<SHOT>.h5` signals and `data/device/diiid.json` geometry: sensor\nmap, conditioned signals, data/design-matrix SVD diagnostics, fit quality, **amplitude & phase of\neach (n, m) mode vs time**, and the **phi-vs-time contour**. To analyse a different shot or array,\nedit the **Parameters** cell (`SHOT`, `CHANNEL_FILTER`, `TIME_TRIM`, `NS`/`MS`)."
   }
  ],
  "metadata": {
@@ -387,4 +340,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 5
-}
+}