From f0a293a4547149d6ebe3326863c85e8df2099ce9 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Thu, 13 Aug 2026 10:30:15 -0700 Subject: [PATCH] GWSignal: fix TEOBResumSDALI global sign (psi was displaced by pi/2) gwsignal's TEOBResumSDALI GenerateTDModes output relates to the same plugin's polarizations by hp - i hc = - nu M G/(c^2 D) * sum_lm h_lm {-2}Y_lm(iota, 0) Measured one (l,|m|) at a time via ModeArray, so that a single complex scalar is exact and an azimuth convention cannot masquerade as a sign: arg(c) = 180 deg with residual 0.0000 for m = 1, 2, 3 and 4 alike, and |c| matching nu M G/(c^2 D) to six digits. A constant arg across m is a global sign, not a phase convention. hlmoft applied the magnitude but not the minus sign. This negates the rescaling and changes nothing else. Other approximants are untouched: RIFT's SEOBNRv5EHM mode sum already reproduces the gwsignal polarizations exactly, every per-m coefficient 1.000 at 0.0 deg. A global sign on h is exactly psi -> psi + pi/2, so this was not visible as a bad fit. It displaced psi by a quarter turn in every RIFT TEOBResumSDALI posterior in the eccentric-PE task force comparison while leaving sky location, distance, masses, inclination and the peak likelihood in agreement with bilby and dingo. Confirmed end to end: before, the template RIFT labels psi_R matches the gwsignal/bilby template at psi_R + 1.03 x pi/2; after, |dpsi| <= 0.10 x pi/2. Same answer in H1, L1 and V1, at phi_ref = 0 and 0.3. Deliberately NOT addressed here: TEOB's modes also do not want the exp(i m phi_shift) that hlmoft applies, which leaves RIFT's reported `phase` for this approximant offset by pi/2 from the polarization path. That is a phase-convention difference rather than a sign error, it is unchanged by this commit, and it is left for separate review since the shift is load-bearing for every other approximant. The regression test asserts the sign, not the convention. Fitting a free complex coefficient c_m per azimuthal index, the templates are correct iff c_m = s exp(i m delta) with s = +1; a common delta is only a relabelling of phi_ref, but s = -1 is the bug. One m cannot separate them -- -exp(i m delta) is reabsorbed into delta, which is why the quadrupole alone never saw this -- so the test uses arg(c_4) - 2 arg(c_2) == 0 if s = +1, == pi if s = -1 which is invariant under any choice of delta. It reads +178.7 deg on the unpatched tree and passes after the fix, and it will keep passing if the phase convention is revisited later. The test skips where the TEOBResumSDALI plugin is absent. This branch has no CI workflow, so run it manually where the plugin exists. Co-Authored-By: Claude Opus 5 --- .../Code/RIFT/physics/GWSignal.py | 14 +- .../Code/test/test_gwsignal_teob_mode_sign.py | 173 ++++++++++++++++++ 2 files changed, 186 insertions(+), 1 deletion(-) create mode 100644 MonteCarloMarginalizeCode/Code/test/test_gwsignal_teob_mode_sign.py diff --git a/MonteCarloMarginalizeCode/Code/RIFT/physics/GWSignal.py b/MonteCarloMarginalizeCode/Code/RIFT/physics/GWSignal.py index d6998bfa6..ec1666556 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/physics/GWSignal.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/physics/GWSignal.py @@ -162,7 +162,19 @@ def hlmoft(P, Lmax=2,approx_string=None,no_trust_align_method=None,internal_phas .value ) for mode in hlmT: - hlmT[mode].data.data = distance_rescaling*hlmT[mode].data.data + # NOTE THE SIGN. gwsignal's TEOBResumSDALI modes are MINUS the + # polarization convention, not plus. Getting this wrong is not + # visible as a bad fit: a global sign on h is exactly + # psi -> psi + pi/2, so it silently displaces the polarization + # angle by a quarter turn and leaves every other parameter, and + # the peak likelihood, looking fine. + # + # Separately, and NOT corrected here: TEOB's modes also do not + # want the exp(i m phi_shift) applied above, which leaves RIFT's + # reported `phase` for this approximant offset by pi/2 from the + # polarization path. That is a phase-convention difference, not + # a sign error, and it is unchanged by this fix. + hlmT[mode].data.data = -distance_rescaling*hlmT[mode].data.data return hlmT diff --git a/MonteCarloMarginalizeCode/Code/test/test_gwsignal_teob_mode_sign.py b/MonteCarloMarginalizeCode/Code/test/test_gwsignal_teob_mode_sign.py new file mode 100644 index 000000000..f668642a5 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_gwsignal_teob_mode_sign.py @@ -0,0 +1,173 @@ +"""Regression guard: RIFT's gwsignal TEOBResumSDALI templates must not carry +a global sign flip relative to the same plugin's polarizations. + +Background +---------- +RIFT reads TEOBResumSDALI through ``GenerateTDModes``; bilby and dingo read it +through ``GenerateFDWaveform`` (polarizations). For this one plugin the modes +are *minus* the polarization convention, and dimensionless rather than in +strain units. ``GWSignal.hlmoft`` rescales by ``nu M G/(c^2 D)``; that factor +must carry the minus sign. + +Getting it wrong is invisible in a fit. The antenna patterns obey +``F+(psi+pi/2) = -F+(psi)`` and ``Fx(psi+pi/2) = -Fx(psi)`` in every detector +for every mode, so a global sign on h is *exactly* ``psi -> psi + pi/2`` and +nothing else: sky location, distance, masses, inclination and the peak +likelihood all still agree. It displaced psi by a quarter turn in every RIFT +TEOBResumSDALI posterior in the eccentric-PE task force comparison. + +What this test asserts, and what it deliberately does not +-------------------------------------------------------- +It tests the *physics* (is there a sign error?), not the phase convention. + +Fitting a free complex coefficient ``c_m`` per azimuthal index m, RIFT's +templates are correct iff ``c_m = s * exp(i m delta)`` with ``s = +1``: a +common ``delta`` is just a relabelling of ``phi_ref``, but ``s = -1`` is the +psi bug. For a single m the two are indistinguishable -- ``-exp(i m delta)`` +can be reabsorbed into ``delta`` -- which is exactly why the quadrupole alone +cannot see this. Two m values sharing one ``delta`` separate them: + + arg(c_4) - 2 arg(c_2) == 0 if s = +1, == pi if s = -1 + +That combination is invariant under any choice of ``delta``, so this test +passes for either TEOB phase convention and fails only on a true sign error. +It is therefore agnostic about whether the ``exp(i m phi_shift)`` in +``hlmoft`` is right for TEOB -- a separate, open question about making +``phase`` line up across approximants. + +A whole-waveform scalar fit would NOT work here: the buggy and correct forms +differ per mode by ``-exp(i m pi/2)``, which is ``+1`` at ``m = +-2``, so the +dominant quadrupole agrees either way and carries nearly all the power. + +Not run in CI: the runners have no TEOBResumSDALI plugin, so this skips there. +Run it where the plugin is installed: + + python -m pytest -q MonteCarloMarginalizeCode/Code/test/test_gwsignal_teob_mode_sign.py +""" +import numpy as np +import pytest + +lal = pytest.importorskip("lal") +lalsim = pytest.importorskip("lalsimulation") +u = pytest.importorskip("astropy.units") + +gws = pytest.importorskip("lalsimulation.gwsignal") +wfm = pytest.importorskip("lalsimulation.gwsignal.core.waveform") + +import RIFT.lalsimutils as lsu # noqa: E402 +import RIFT.physics.GWSignal as rgws # noqa: E402 + +M1, M2 = 80.0, 40.0 # q = 2, so odd-m modes are alive +DIST, IOTA = 1500.0, np.pi / 4 +DELTAT, DELTAF = 1.0 / 4096, 1.0 / 16 +F22 = FREF = 20.0 + + +def _pdict(): + return { + "mass1": M1 * u.solMass, "mass2": M2 * u.solMass, + "spin1x": 0.0 * u.dimensionless_unscaled, + "spin1y": 0.0 * u.dimensionless_unscaled, + "spin1z": 0.0 * u.dimensionless_unscaled, + "spin2x": 0.0 * u.dimensionless_unscaled, + "spin2y": 0.0 * u.dimensionless_unscaled, + "spin2z": 0.0 * u.dimensionless_unscaled, + "deltaT": DELTAT * u.s, "f22_start": F22 * u.Hz, + "f22_ref": FREF * u.Hz, "phi_ref": 0.0 * u.rad, + "distance": DIST * u.Mpc, "inclination": IOTA * u.rad, + "eccentricity": 0.0 * u.dimensionless_unscaled, + "longAscNodes": 0.0 * u.rad, "meanPerAno": 0.0 * u.rad, + "condition": 0, + } + + +def _P(): + P = lsu.ChooseWaveformParams() + P.m1, P.m2 = M1 * lal.MSUN_SI, M2 * lal.MSUN_SI + P.s1x = P.s1y = P.s1z = P.s2x = P.s2y = P.s2z = 0.0 + P.dist = DIST * 1e6 * lal.PC_SI + P.incl = IOTA + P.phiref = P.psi = 0.0 + P.fmin, P.fref = F22, FREF + P.deltaT, P.deltaF = DELTAT, DELTAF + P.eccentricity, P.meanPerAno = 0.0, 0.0 + P.taper = lalsim.SIM_INSPIRAL_TAPER_NONE + P.approx = lalsim.IMRPhenomXPHM # unused; approx_string drives gwsignal + return P + + +def _arr(x): + return np.asarray(x.value if hasattr(x, "value") else x) + + +def _times(x): + t = x.times + return np.asarray(t.value if hasattr(t, "value") else t) + + +def _regrid(ts, a, grid): + return a[np.clip(np.searchsorted(ts, grid), 0, len(a) - 1)] + + +def _coeffs_per_m(approx): + """Free complex coefficient per azimuthal index, fitting the gwsignal + polarizations with RIFT's own modes.""" + gen = gws.models.gwsignal_get_waveform_generator(approx) + hp, hc = wfm.GenerateTDWaveform(_pdict(), gen) + t_pol = _times(hp) + + hlm = rgws.hlmoft(_P(), Lmax=4, approx_string=approx) + key0 = sorted(hlm)[0] + t_rift = (float(hlm[key0].epoch) + + np.arange(hlm[key0].data.length) * hlm[key0].deltaT) + + grid = np.arange(max(t_rift[0], t_pol[0]), min(t_rift[-1], t_pol[-1]), + DELTAT) + target = (_regrid(t_pol, _arr(hp), grid) + - 1j * _regrid(t_pol, _arr(hc), grid)) + + ms = sorted({k[1] for k in hlm}) + cols = [] + for m in ms: + b = np.zeros(len(t_rift), dtype=complex) + for k in hlm: + if k[1] == m: + b += hlm[k].data.data * lal.SpinWeightedSphericalHarmonic( + IOTA, 0.0, -2, k[0], k[1]) + cols.append(_regrid(t_rift, b, grid)) + + A = np.array(cols).T + keep = np.abs(target) > 0.01 * np.abs(target).max() + c, *_ = np.linalg.lstsq(A[keep], target[keep], rcond=None) + return dict(zip(ms, c)) + + +def _sign_invariant(c): + """arg(c_4) - 2 arg(c_2), in radians, folded to [-pi, pi]. + + 0 => no global sign error (any phi_ref convention). +-pi => psi is + displaced by pi/2. + """ + z = c[4] * np.conj(c[2]) ** 2 + return np.angle(z / np.abs(z)) + + +@pytest.mark.parametrize("approx", ["TEOBResumSDALI", "SEOBNRv5EHM"]) +def test_no_global_sign_error_vs_polarizations(approx): + try: + c = _coeffs_per_m(approx) + except Exception as exc: # plugin absent / model failed + pytest.skip("%s unavailable: %r" % (approx, exc)) + + for m in (2, 4): + assert m in c, "need m=%d to separate a sign from a phase" % m + + d = _sign_invariant(c) + assert abs(d) < np.radians(45), ( + "%s: global sign error -- RIFT's templates are -1 x the gwsignal " + "polarizations, i.e. psi is displaced by pi/2.\n" + " arg(c_4) - 2 arg(c_2) = %+.1f deg (expected ~0, pi means the bug)\n" + " per-m coefficients: %s" + % (approx, np.degrees(d), + " ".join("m=%+d: %.3f/%+.1fd" % (m, abs(v), np.degrees(np.angle(v))) + for m, v in sorted(c.items()))))