Measurement-referenced chi2 coil filter and device registry - #49
Open
d-burg wants to merge 16 commits into
Open
Conversation
WIP checkpoint before end-to-end validation on a real archive.
The docstring claimed |I_base|/|I_meas| is the per-coil turns factor. It is not. Across 26 slices of DIII-D 174823 that ratio is stable only for F6A/F6B/F7A/F7B (54-56, 8-12% spread) and varies 59-295% for 16 of 24 coils: the free-boundary solve fits the BOUNDARY and lands on one of many coil sets consistent with it. The metric is unaffected -- sigma^base is formed per slice from that slice's own baseline and measurement, so z_i is (fractional drift)/(measured fractional precision) and stays self-consistent. Only the stated justification was wrong.
coil_sigma_in_base_units rescales by the INSTANTANEOUS |I_base|/|I_meas|. When a coil's measured current passes near zero that ratio explodes and the coil gets an absurdly loose tolerance, even though the reconstruction's own current is no less determined at that instant. Measured on DIII-D 174823 (worst per-slice sigma / typical sigma): F4B 9.7x ECOILB 8.4x ECOILA 8.2x F5B 4.2x F8A 3.0x F5A 2.7x The concrete case: F5B at t=5.706 s has I_meas = -67.8 A -- above the 50 A floor, so not caught -- against a 2.3e4 A-t baseline, giving sigma = 2426 A-t. coil_sigma_fixed takes the conversion factor as a robust median over the discharge, using only slices where the coil carries enough current for the ratio to mean anything, so sigma is constant in time. Typical values are unchanged (within a few percent on 174823; F6A/F6B/F7A/F7B identical to 3 s.f.) -- only the outliers go away. 4 tests. test_immune_to_a_near_zero_crossing reproduces the real F5B numbers and asserts the 50 A floor does NOT exclude it, which is why the floor alone is not a fix. 404 tests pass.
…eans the same across DAQ epochs pf_active data_error_upper is 10 digitizer LSB (OMAS digi_error) and is ~8x larger on pre-2017 hardware, so a dd-referenced cut accepts everything on e.g. 150000. sigma_ref='d3d' uses the 2017 table chi2_max=4.0 was calibrated against. Default unchanged.
…default; dd_path optional The machine tolerance from EFIT's own calculated-minus-measured F-coil residuals (72 coil-shots, FWTFC=0). The digitizer table stays as the strict 'd3d' option; the dd data_error_upper as None.
…lt; legacy 2% rule selectable via filtering.coil_filter='legacy'
…llback
- bouquet/devices.py: DeviceSpec registry (DIII-D only today): coil signature,
tolerance model floor+fraction with provenance, turns, families, digitizer
table, VSC pair. detect_device() matches the exact mesh coil-set signature.
- resolve_coil_sigma(): explicit filtering.coil_sigma ({floor,fraction} model,
per-coil table, or callable) > device model (named or detected) > raises
CoilSigmaUnavailable.
- Bouquet.filter(): chi2 default; on CoilSigmaUnavailable warns loudly and
uses the legacy rule; records coil_filter_used.
- filter_coil_chi2 stamps coil_filter / coil_sigma_model on the scan group;
dd-referenced sigma_ref modes kept, explicit and requiring dd_path.
- BouquetConfig.device optional.
…35%; 825 A-t before shot 165000); rms model kept as named option Per-shot fits over 497 DIII-D shots (105 CTM + 392 IBS): the offset-removed fraction is ~0.3% in every era, the floor is a DAQ-era effect. The baseline already absorbs each coil's systematic offset, so the draw tolerance is the random part. filtering.coil_sigma='rms_incl_offset' restores the 4-shot rms model. Shot inferred from source/geqdsk name/header for the era lookup.
… |z| <= 5 guard Per-coil random floors with the 0.35% fraction removed in quadrature, median over shots per era, clipped at 100 A-t (modern) / 250 A-t (pre-2014); coils outside the table keep the era floor. F6A/F6B (~800 A-t) and F9A (~580 A-t) are the only modern coils with a real floor. z_max guards against one coil at 7 sigma hiding behind seventeen quiet ones.
…emoved std) and the untested symmetry assumption
…avy tail, drift-dominated std, coils not independent)
…e does not absorb EFIT's offset (25-30 sigma away along the coil null space); std is used because it is the fit-independent part
…ribution (DIII-D q95: chi2/nu <= 6.1, worst |z| <= 6.3); generic 4 / 5 when no device calibration applies Both thresholds default to None -> device acceptance when sigma is the device model, else GENERIC_ACCEPTANCE. z_max=False disables the guard. Provenance carried in the device entry and stamped on the archive.
There was a problem hiding this comment.
🟡 Changes recommended
The new chi² filter has at least one confirmed functional bug (flat/legacy archive layout handling) and a confirmed acceptance-threshold behavior mismatch with the documented intent that can silently change filtering outcomes.
Once you've addressed the issues Copilot identified, you can request another Copilot review.
Pull request overview
This PR replaces the legacy per-draw coil-current band filter with a measurement-referenced per-draw χ²/ν + worst-coil |z| criterion, backed by a device registry (for DIII-D: calibrated coil-current tolerance model + acceptance thresholds) and provenance stamping into the archive.
Changes:
- Add
bouquet.coil_spec(χ² metric, sigma resolution, device/DD-referenced sigma utilities) andbouquet.devices(device registry + tolerance/acceptance calibration). - Implement
filter_coil_chi2(andmeasured_coil_currents) as a new postprocessing coil filter, and make it the default inBouquet.filter(). - Add tests covering sigma transfer, χ² behavior, device detection/acceptance, provenance stamping, and shot inference.
File summaries
| File | Description |
|---|---|
| tests/test_golden_bouquet.py | Adds integration-style tests asserting chi² is the default filter and that provenance + fallback behavior work. |
| tests/test_coil_spec.py | Adds unit tests for χ² scoring, sigma models/resolution, device registry behavior, and acceptance logic. |
| bouquet/run.py | Switches Bouquet.filter() to default to the chi² coil filter with fallback to legacy and adds shot inference helper. |
| bouquet/filtering.py | Adds measured_coil_currents and filter_coil_chi2, including stamping filter provenance into scan groups. |
| bouquet/devices.py | Introduces a device registry (DeviceSpec) and DIII-D calibration data (sigma + acceptance + signatures). |
| bouquet/config.py | Adds config knobs for selecting chi² vs legacy filtering and for providing/overriding sigma/acceptance. |
| bouquet/coil_spec.py | Implements core χ² metric and sigma-resolution utilities (explicit, per-coil, callable, device model, dd-referenced). |
| bouquet/init.py | Exposes the new device/sigma/filter APIs at package top-level. |
Review details
- Files reviewed: 8/8 changed files
- Comments generated: 5
- Review effort level: Lite
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
Comment on lines
+335
to
+342
| with h5py.File(h5path, "r") as hf: | ||
| grp = hf["scan"][str(sv)] | ||
| if "_baseline" not in grp: | ||
| continue | ||
| bn = [x.decode() if isinstance(x, bytes) else str(x) | ||
| for x in grp["_baseline"]["coil_names"][()]] | ||
| baseline = dict(zip(bn, np.asarray( | ||
| grp["_baseline"]["coil_currents"][()], dtype=float).tolist())) |
Comment on lines
+352
to
+356
| acc = dict(GENERIC_ACCEPTANCE); acc_src = "generic" | ||
| if model.get("kind") == "device" and model.get("model") == "random": | ||
| dacc = get_device(model["device"]).acceptance | ||
| if dacc: | ||
| acc.update({k: dacc[k] for k in ("chi2_max", "z_max") if k in dacc}); acc_src = "device q%g" % (100 * dacc.get("quantile", float("nan"))) |
Comment on lines
+1405
to
1409
| self._selection = {"coil": coil_summary, "boundary": bnd_summary, | ||
| "coil_filter_used": coil_filter_used} | ||
| if plot: | ||
| self._selection["figures"] = (coil_fig, bnd_fig) | ||
| self._print_generation_summary(coil_summary, bnd_summary) |
Comment on lines
+503
to
+505
| # "chi2" -> measurement-referenced chi2/nu <= chi2_max with the per-coil | ||
| # sigma from coil_sigma_ref ("efit" = EFIT-residual floor+fraction, | ||
| # no dd needed; "d3d" = digitizer table, needs a dd) [DEFAULT] |
Comment on lines
+66
to
+69
| # coil-current tolerance model sigma_i = hypot(floor, fraction*|I_i|), baseline units. | ||
| # This is the RANDOM part of the reconstruction's coil-current residual (per-shot | ||
| # systematic offsets are already absorbed by the baseline fit, so a draw about the | ||
| # baseline must not re-explore them). The floor may depend on the DAQ era: give |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Replaces the legacy
|dI/I| <= 2 %draw filter with a per-draw χ² test against a per-coil coil-current tolerance, plus a device registry.Tolerance. σ_i² = floor_i² + (0.35 %·|I_i|)² in ampere-turns. For DIII-D the numbers are the offset-removed std of EFIT's own calculated-minus-measured coil current over the flat-top (fit weights zero), from 497 shots: the fraction is the same in every era, the floor is DAQ-era dependent (825 A-t before shot 165000, 325 after) with per-coil floors. The per-shot mean offset is excluded on purpose (it is a property of the reconstruction's null-space choice, not of the machine);
filtering.coil_sigma="rms_incl_offset"selects the looser model. Users may also give{"floor","fraction"}, a per-coil table, or a callable.Acceptance. χ²/ν ≤ 6.1 and worst-coil |z| ≤ 6.3 for DIII-D: the 95th percentile of what real flat-top slices score under the same σ, i.e. a 5 % false-rejection rate on real machine states. Generic 4 / 5 when no device calibration applies; both overridable;
z_max=Falsedisables the guard.Device handling.
bouquet/devices.py: device named in config or detected from the exact mesh coil signature; no signature and no explicit σ →Bouquet.filter()warns loudly and applies the legacy rule.Provenance. Model and thresholds stamped on the archive scan group. Nothing is fetched at run time.
Checked and documented in the module docstring: the residual is symmetric/Gaussian in the core with a heavy far tail; ~63 % of its variance is slow drift (kept); the 18 coil residuals are correlated (n_eff ≈ 4). An unregularised baseline sits ~10 kA-t from the measured currents along the coil null space, which this filter cannot see — addressed by the coil-reg targets PR stacked on this one.
🤖 Generated with Claude Code