From 1b679032294001b2cd1d5e2420fb30ba68d51c43 Mon Sep 17 00:00:00 2001 From: Yun Date: Tue, 8 Sep 2026 20:05:12 +0900 Subject: [PATCH 1/2] Harmonize toroidal mode-number sign convention across magnetics routines (fix #638) - Update toroidal_mode_analysis in vaft/process/magnetics.py to compute n_raw = -phase / float(phase_geometry), matching toroidal_phase_fit_at_time and SXR rank_toroidal_mode_numbers. - Explicitly document right-handed cylindrical coordinates (R, phi, Z) and positive n propagation in module and function docstrings. - Update test_mirnov_pipeline.py synthetic phase lag for expected_n = 2. - Add comprehensive test suite in test/test_process_magnetics.py asserting exact agreement between pair analysis and array fitting across positive, negative, zero, and concurrent modes. --- test/test_mirnov_pipeline.py | 2 +- test/test_process_magnetics.py | 181 +++++++++++++++++++++++++++++++++ vaft/process/magnetics.py | 37 ++++--- 3 files changed, 207 insertions(+), 13 deletions(-) create mode 100644 test/test_process_magnetics.py diff --git a/test/test_mirnov_pipeline.py b/test/test_mirnov_pipeline.py index 20ba33c63..52f95ee83 100644 --- a/test/test_mirnov_pipeline.py +++ b/test/test_mirnov_pipeline.py @@ -141,7 +141,7 @@ def test_toroidal_mode_analysis_recovers_phase_mode(): phase_geometry = np.pi / 6 expected_n = 2 signal_a = np.sin(2.0 * np.pi * 1_000.0 * time) - signal_b = np.sin(2.0 * np.pi * 1_000.0 * time + expected_n * phase_geometry) + signal_b = np.sin(2.0 * np.pi * 1_000.0 * time - expected_n * phase_geometry) result = toroidal_mode_analysis( signal_a, diff --git a/test/test_process_magnetics.py b/test/test_process_magnetics.py new file mode 100644 index 000000000..5b806e3ed --- /dev/null +++ b/test/test_process_magnetics.py @@ -0,0 +1,181 @@ +"""Unit tests for toroidal mode analysis and phase fitting sign conventions (Issue #638). + +Verifies that: +1. `toroidal_mode_analysis` and `toroidal_phase_fit_at_time` yield identical mode + numbers for identical synthetic signals across positive, negative, and zero modes. +2. Positive n corresponds to co-current / +phi propagation (downstream phase lag). +3. Negative n corresponds to counter-current / -phi propagation (downstream phase lead). +4. Geometry sign inversions (Delta_phi < 0) are handled consistently. +5. Edge cases including phase wrapping and multiple simultaneous modes are recovered. +""" + +from __future__ import annotations + +import numpy as np +import pytest + +from vaft.process.magnetics import ( + toroidal_mode_analysis, + toroidal_phase_fit_at_time, +) + + +@pytest.mark.parametrize("n_expected", [1, 2, -1, -2, 0]) +def test_toroidal_mode_analysis_and_fit_agree(n_expected: int): + """Both entry points must yield the exact same mode number for identical signals.""" + sample_rate = 50_000.0 + n_samples = 4096 + time = np.arange(n_samples, dtype=float) / sample_rate + frequency = 5_000.0 + initial_phase = 0.35 + + # 4 sensors spaced around the torus + angles = np.deg2rad([0.0, 30.0, 60.0, 90.0]) + phase_geometry = float(angles[1] - angles[0]) + + # Under canonical convention: phase(phi) = initial_phase - n * phi + phases = initial_phase - float(n_expected) * angles + signals = np.vstack( + [np.sin(2.0 * np.pi * frequency * time + phi) for phi in phases] + ) + + # 1. Pair analysis on probe 0 and probe 1 + res_pair = toroidal_mode_analysis( + signals[0], + signals[1], + sample_rate=sample_rate, + phase_geometry=phase_geometry, + peak_threshold=0.05, + nperseg=1024, + ) + assert len(res_pair.n) > 0, "Expected at least one coherent peak in pair analysis" + # Find peak closest to target frequency + freq_idx = int(np.argmin(np.abs(res_pair.frequency - frequency))) + n_pair = int(res_pair.n[freq_idx]) + + # 2. Wrapped phase array fit across all 4 probes + res_fit = toroidal_phase_fit_at_time( + time, + signals, + angles, + center_time=0.04, + sample_rate=sample_rate, + window_size=512, + frequencies=[frequency], + candidate_n=tuple(range(-4, 5)), + ) + assert len(res_fit.modes) == 1, "Expected one fitted mode" + n_fit = res_fit.modes[0].n + + # Assert both routines match the expected physical mode number and each other + assert n_pair == n_expected, f"Pair analysis got n={n_pair}, expected {n_expected}" + assert n_fit == n_expected, f"Array fit got n={n_fit}, expected {n_expected}" + assert n_pair == n_fit, f"Pair analysis (n={n_pair}) != Array fit (n={n_fit})" + + +def test_co_current_and_counter_current_propagation(): + """Verify co-current (+phi) is n > 0 and counter-current (-phi) is n < 0.""" + sample_rate = 100_000.0 + time = np.arange(4096, dtype=float) / sample_rate + freq = 8_000.0 + dphi = np.pi / 4 # 45 degrees + + # Co-current wave propagating in +phi: downstream probe lags upstream probe + # signal_upstream = sin(omega * t) + # signal_downstream = sin(omega * t - delta) with delta > 0 + sig_up = np.sin(2.0 * np.pi * freq * time) + sig_down_lag = np.sin(2.0 * np.pi * freq * time - 2 * dphi) # n = +2 + res_co = toroidal_mode_analysis( + sig_up, + sig_down_lag, + sample_rate=sample_rate, + phase_geometry=dphi, + peak_threshold=0.05, + nperseg=1024, + ) + assert 2 in set(res_co.n.astype(int)) + + # Counter-current wave propagating in -phi: downstream probe leads upstream probe + sig_down_lead = np.sin(2.0 * np.pi * freq * time + 2 * dphi) # n = -2 + res_counter = toroidal_mode_analysis( + sig_up, + sig_down_lead, + sample_rate=sample_rate, + phase_geometry=dphi, + peak_threshold=0.05, + nperseg=1024, + ) + assert -2 in set(res_counter.n.astype(int)) + + +def test_negative_phase_geometry(): + """When coil B is at smaller toroidal angle than coil A (dphi < 0), n is preserved.""" + sample_rate = 50_000.0 + time = np.arange(4096, dtype=float) / sample_rate + freq = 6_000.0 + phi_a = np.pi / 3 + phi_b = 0.0 + dphi = phi_b - phi_a # -pi/3 < 0 + expected_n = 2 + + sig_a = np.sin(2.0 * np.pi * freq * time - expected_n * phi_a) + sig_b = np.sin(2.0 * np.pi * freq * time - expected_n * phi_b) + + res = toroidal_mode_analysis( + sig_a, + sig_b, + sample_rate=sample_rate, + phase_geometry=dphi, + peak_threshold=0.05, + nperseg=1024, + ) + assert expected_n in set(res.n.astype(int)) + + +def test_multiple_simultaneous_modes(): + """Verify both routines on a composite signal containing two distinct mode frequencies.""" + sample_rate = 100_000.0 + time = np.arange(8192, dtype=float) / sample_rate + f1, n1 = 4_000.0, 1 + f2, n2 = 12_000.0, 3 + + angles = np.deg2rad([0.0, 30.0, 60.0, 90.0]) + dphi = float(angles[1] - angles[0]) + + signals = np.vstack( + [ + np.sin(2.0 * np.pi * f1 * time - n1 * phi) + + 0.8 * np.sin(2.0 * np.pi * f2 * time - n2 * phi) + for phi in angles + ] + ) + + # Pair analysis + res_pair = toroidal_mode_analysis( + signals[0], + signals[1], + sample_rate=sample_rate, + phase_geometry=dphi, + peak_threshold=0.05, + nperseg=2048, + ) + # Check that both n1 and n2 are present at their respective frequencies + idx1 = int(np.argmin(np.abs(res_pair.frequency - f1))) + idx2 = int(np.argmin(np.abs(res_pair.frequency - f2))) + assert int(res_pair.n[idx1]) == n1 + assert int(res_pair.n[idx2]) == n2 + + # Array fit + res_fit = toroidal_phase_fit_at_time( + time, + signals, + angles, + center_time=0.04, + sample_rate=sample_rate, + window_size=1024, + frequencies=[f1, f2], + candidate_n=tuple(range(-4, 5)), + ) + fit_modes = {round(m.frequency, -2): m.n for m in res_fit.modes} + assert fit_modes[round(f1, -2)] == n1 + assert fit_modes[round(f2, -2)] == n2 diff --git a/vaft/process/magnetics.py b/vaft/process/magnetics.py index e08ddc476..6e0584d06 100644 --- a/vaft/process/magnetics.py +++ b/vaft/process/magnetics.py @@ -45,9 +45,17 @@ which is the one place in the module where a convention is decided by the signal rather than declared. -**The two mode-number entry points disagree on the sign of n**, one carrying a -minus that the other does not; see #638. Compare results from them only after -reading that issue. +**Toroidal mode numbers adhere to standard right-handed cylindrical coordinates.** +Coordinates (R, phi, Z) are oriented such that phi increases counter-clockwise +when viewed from above. A positive toroidal mode number (n > 0) denotes a +perturbation propagating in the positive phi direction (co-current / toroidal +direction). Across toroidally separated sensors, the signal phase varies as +theta(phi) = theta_0 - n * phi (negative phase slope d(theta)/d(phi) = -n). For +two coils separated by toroidal angle Delta_phi = phase_geometry > 0, the +cross-spectral phase of the downstream coil relative to the upstream coil is +Delta_theta = -n * Delta_phi, so the mode number is n = -Delta_theta / Delta_phi. +Both :func:`toroidal_mode_analysis` and :func:`toroidal_phase_fit_at_time` +share this convention. Provenance ---------- @@ -1225,14 +1233,18 @@ def toroidal_mode_analysis( 1. Estimate the cross-spectral density and the coherence of the pair. 2. Keep frequencies whose coherence exceeds the significance level. 3. Find the peaks among them above the threshold. - 4. Divide each peak's phase by the toroidal separation to get its mode number. + 4. Divide the negative of each peak's phase by the toroidal separation to get its mode number. Convention ---------- - **The mode number is the phase over the separation, with no minus sign.** - :func:`toroidal_phase_fit_at_time` fits a model that carries one, so **the two - entry points report opposite signs for the same physical mode**; tracked in - #638. Compare a result from one against the other only after reading it. + **The mode number carries a minus sign**: ``n_raw = -phase / phase_geometry``. + Under standard right-handed cylindrical coordinates (R, phi, Z) with phi + counter-clockwise from above, a mode propagating in the co-current (+phi) + direction has phase decreasing with increasing phi (theta(phi) = theta_0 - + n * phi, so d(theta)/d(phi) = -n). A sensor at Delta_phi > 0 lags the + reference sensor by Delta_theta = -n * Delta_phi, so n = -Delta_theta / + Delta_phi. This agrees with the model fitted by + :func:`toroidal_phase_fit_at_time`. The coherence threshold is the 95 percent significance level for a magnitude-squared coherence averaged over the given number of sensors, so it @@ -1273,7 +1285,7 @@ def toroidal_mode_analysis( frequencies, cross_power = csd(a, b, fs=sample_rate, nperseg=segment) _, coherence_values = coherence(a, b, fs=sample_rate, nperseg=segment) phase = np.angle(cross_power) - n_raw = phase / float(phase_geometry) + n_raw = -phase / float(phase_geometry) n_rounded = np.round(n_raw) power_abs = np.abs(cross_power) @@ -1381,9 +1393,10 @@ def toroidal_phase_fit_at_time( Convention ---------- **The model carries a minus sign**: the fitted phase decreases with increasing - toroidal angle for a positive mode number. :func:`toroidal_mode_analysis` uses - the opposite sense, so **the two report opposite signs for the same physical - mode**; tracked in #638. + toroidal angle for a positive mode number (``fitted = intercept - n * toroidal_angle``). + Under standard right-handed cylindrical coordinates (R, phi, Z), this + corresponds to a perturbation propagating in the positive phi (co-current) + direction, in agreement with :func:`toroidal_mode_analysis`. Phases are wrapped to a single turn and the intercept is a circular mean, so the fit is insensitive to where the branch cut falls, which a plain average From 1d2c428cc728efcc024f8cee389559b6b82057bd Mon Sep 17 00:00:00 2001 From: Yun Date: Tue, 8 Sep 2026 20:11:20 +0900 Subject: [PATCH 2/2] Address adversarial review feedback for toroidal mode sign harmonization - Update user guide formula in docs/_guide/Processing.md to include leading minus. - Update candidate_n default in vaft/plot/mirnov.py to range(-6, 7). - Update comment in test/test_process_docstrings.py reflecting harmonized status. - Clarify phase_geometry docstring (phi_b - phi_a) and use np.isclose in magnetics.py. - Tighten test assertions with frequency indexing and add tests for descending angles, noise, and near-Nyquist cases in test_process_magnetics.py. --- docs/_guide/Plotting.md | 2 +- docs/_guide/Processing.md | 2 +- test/test_process_docstrings.py | 4 +- test/test_process_magnetics.py | 96 +++++++++++++++++++++++++++++++-- vaft/plot/mirnov.py | 2 +- vaft/process/magnetics.py | 4 +- 6 files changed, 100 insertions(+), 10 deletions(-) diff --git a/docs/_guide/Plotting.md b/docs/_guide/Plotting.md index b4edcf6a9..9fea1e2bb 100644 --- a/docs/_guide/Plotting.md +++ b/docs/_guide/Plotting.md @@ -417,7 +417,7 @@ your own subplot grids. Everything after the first argument is keyword-only. | `mirnov_signal` | `(ods, channels=None, *, probe_group='b_field_pol_probe', time_range=None, preprocess=False, gains=None, ax=None, show=True)` | | `mirnov_spectrogram` | `(ods, channel=0, *, probe_group='b_field_pol_probe', time_range=None, preprocess=True, gain=None, sample_rate=None, window_size=500, time_resolution=1, max_frequency=None, cmap='hot_r', ax=None, show=True, return_result=False)` | | `toroidal_mode_spectrum` | `(ods, channel_pair=(65, 67), *, probe_group='b_field_pol_probe', time_range=None, preprocess=True, gains=None, phase_geometry=np.pi/6, peak_threshold=0.1, sample_rate=None, axes=None, show=True, return_result=False)` | -| `toroidal_phase_mode_fit` | `(ods, center_time, *, channels=(64, 65, 66, 67), probe_group='b_field_pol_probe', time_range=None, frequencies=None, num_modes=2, candidate_n=tuple(range(0, 7)), window_size=500, preprocess=True, gains=None, sample_rate=None, peak_threshold=0.1, ax=None, show=True, save_path=None, return_result=False)` | +| `toroidal_phase_mode_fit` | `(ods, center_time, *, channels=(64, 65, 66, 67), probe_group='b_field_pol_probe', time_range=None, frequencies=None, num_modes=2, candidate_n=tuple(range(-6, 7)), window_size=500, preprocess=True, gains=None, sample_rate=None, peak_threshold=0.1, ax=None, show=True, save_path=None, return_result=False)` | ```python import matplotlib.pyplot as plt diff --git a/docs/_guide/Processing.md b/docs/_guide/Processing.md index e330f68af..1e673ab9f 100644 --- a/docs/_guide/Processing.md +++ b/docs/_guide/Processing.md @@ -622,7 +622,7 @@ best = fit.modes[0] # sorted by amp print(best.frequency, best.n, best.rms_error) ``` -$n$ is recovered as $\arg \mathrm{CSD}(a,b) / \Delta\phi$, peak-picked on $\lvert \mathrm{CSD} \rvert$ +$n$ is recovered as $-\arg \mathrm{CSD}(a,b) / \Delta\phi$, peak-picked on $\lvert \mathrm{CSD} \rvert$ and filtered by a coherence threshold. The plot module wraps all of this against an ODS — this is the path the fluctuation notebook takes on shot 44740: diff --git a/test/test_process_docstrings.py b/test/test_process_docstrings.py index ed22932a8..f0e781bdf 100644 --- a/test/test_process_docstrings.py +++ b/test/test_process_docstrings.py @@ -187,8 +187,8 @@ # and is never converted "pedestal_top", # magnetics / electromagnetics / fluctuation (#418): integration sign, - # shot-era baselines, per-unit-current responses, and the two mode-number - # entry points that disagree on the sign of n (#638) + # shot-era baselines, per-unit-current responses, and the toroidal mode-number + # entry points harmonized under standard right-handed coordinates (#638) "analyze_fluctuation_spectrum", "b_field_pol_probe_field", "calc_grid", diff --git a/test/test_process_magnetics.py b/test/test_process_magnetics.py index 5b806e3ed..c9f1dfdad 100644 --- a/test/test_process_magnetics.py +++ b/test/test_process_magnetics.py @@ -93,7 +93,8 @@ def test_co_current_and_counter_current_propagation(): peak_threshold=0.05, nperseg=1024, ) - assert 2 in set(res_co.n.astype(int)) + idx_co = int(np.argmin(np.abs(res_co.frequency - freq))) + assert int(res_co.n[idx_co]) == 2 # Counter-current wave propagating in -phi: downstream probe leads upstream probe sig_down_lead = np.sin(2.0 * np.pi * freq * time + 2 * dphi) # n = -2 @@ -105,7 +106,8 @@ def test_co_current_and_counter_current_propagation(): peak_threshold=0.05, nperseg=1024, ) - assert -2 in set(res_counter.n.astype(int)) + idx_counter = int(np.argmin(np.abs(res_counter.frequency - freq))) + assert int(res_counter.n[idx_counter]) == -2 def test_negative_phase_geometry(): @@ -129,7 +131,95 @@ def test_negative_phase_geometry(): peak_threshold=0.05, nperseg=1024, ) - assert expected_n in set(res.n.astype(int)) + idx = int(np.argmin(np.abs(res.frequency - freq))) + assert int(res.n[idx]) == expected_n + + +def test_decreasing_angle_array_fit(): + """Array fit with angles listed in descending order recovers identical n.""" + sample_rate = 50_000.0 + time = np.arange(4096, dtype=float) / sample_rate + freq = 7_000.0 + expected_n = 3 + angles = np.deg2rad([90.0, 60.0, 30.0, 0.0]) + phases = 0.2 - expected_n * angles + signals = np.vstack([np.sin(2.0 * np.pi * freq * time + p) for p in phases]) + + res = toroidal_phase_fit_at_time( + time, + signals, + angles, + center_time=0.04, + sample_rate=sample_rate, + window_size=512, + frequencies=[freq], + candidate_n=tuple(range(-5, 6)), + ) + assert len(res.modes) == 1 + assert res.modes[0].n == expected_n + + +def test_near_nyquist_boundary(): + """Mode number near Nyquist (n=5 with dphi=pi/6, phase separation 5pi/6) is recovered.""" + sample_rate = 100_000.0 + time = np.arange(4096, dtype=float) / sample_rate + freq = 10_000.0 + expected_n = 5 + dphi = np.pi / 6 # 30 deg; Nyquist is n=6 + + sig_a = np.sin(2.0 * np.pi * freq * time) + sig_b = np.sin(2.0 * np.pi * freq * time - expected_n * dphi) + + res = toroidal_mode_analysis( + sig_a, + sig_b, + sample_rate=sample_rate, + phase_geometry=dphi, + peak_threshold=0.05, + nperseg=1024, + ) + idx = int(np.argmin(np.abs(res.frequency - freq))) + assert int(res.n[idx]) == expected_n + + +def test_noisy_mode_signal_recovery(): + """Additive Gaussian noise does not perturb rounded mode recovery.""" + rng = np.random.default_rng(42) + sample_rate = 100_000.0 + time = np.arange(8192, dtype=float) / sample_rate + freq = 8_000.0 + expected_n = 2 + angles = np.deg2rad([0.0, 30.0, 60.0, 90.0]) + dphi = float(angles[1] - angles[0]) + + pure_signals = np.vstack( + [np.sin(2.0 * np.pi * freq * time - expected_n * phi) for phi in angles] + ) + noise = 0.15 * rng.standard_normal(pure_signals.shape) + signals = pure_signals + noise + + res_pair = toroidal_mode_analysis( + signals[0], + signals[1], + sample_rate=sample_rate, + phase_geometry=dphi, + peak_threshold=0.05, + nperseg=2048, + ) + idx_pair = int(np.argmin(np.abs(res_pair.frequency - freq))) + assert int(res_pair.n[idx_pair]) == expected_n + + res_fit = toroidal_phase_fit_at_time( + time, + signals, + angles, + center_time=0.04, + sample_rate=sample_rate, + window_size=1024, + frequencies=[freq], + candidate_n=tuple(range(-4, 5)), + ) + assert res_fit.modes[0].n == expected_n def test_multiple_simultaneous_modes(): diff --git a/vaft/plot/mirnov.py b/vaft/plot/mirnov.py index a03421e4b..5200df320 100644 --- a/vaft/plot/mirnov.py +++ b/vaft/plot/mirnov.py @@ -418,7 +418,7 @@ def toroidal_phase_mode_fit( time_range: tuple[float, float] | None = None, frequencies: Sequence[float] | None = None, num_modes: int = 2, - candidate_n: Sequence[int] = tuple(range(0, 7)), + candidate_n: Sequence[int] = tuple(range(-6, 7)), window_size: int = 500, preprocess: bool = True, gains: Any = None, diff --git a/vaft/process/magnetics.py b/vaft/process/magnetics.py index 6e0584d06..a78df969d 100644 --- a/vaft/process/magnetics.py +++ b/vaft/process/magnetics.py @@ -1215,7 +1215,7 @@ def toroidal_mode_analysis( sample_rate : float, optional Sample rate [Hz]. phase_geometry : float, optional - Toroidal separation of the two coils [rad]. + Toroidal angle of the second coil minus the first coil, phi_b - phi_a [rad]. peak_threshold : float, optional Fraction of the maximum a peak must reach to be reported [-]. sensor_count : int, optional @@ -1278,7 +1278,7 @@ def toroidal_mode_analysis( if a.size < 2: empty = np.array([]) return ToroidalModeResult(empty, empty, empty, empty, empty, empty, empty.astype(int), empty, empty, empty) - if phase_geometry == 0: + if np.isclose(float(phase_geometry), 0.0): raise ValueError("phase_geometry must be non-zero.") segment = min(a.size, int(nperseg) if nperseg is not None else 256)