From 06764daf12a55077562c15cb28c2228654ae6d16 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Fri, 7 Aug 2026 06:49:06 -0700 Subject: [PATCH] =?UTF-8?q?lalsimutils:=20fix=20spurious=20(l,=C2=B1m)=20a?= =?UTF-8?q?symmetry=20in=20hlmoft=20FD-mode=20conditioning?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit SimInspiralChooseFDModes returns modes on an ascending two-sided grid [-fNyq, +fNyq] with DC at the center index (odd length TDlen+1), but the FD high-pass window in hlmoft assigned frequencies via evaluate_fvals, which assumes RIFT's reversed packing and, at odd length, is offset by deltaF/2 (f_assumed = -f_true + deltaF/2). Because (l,m) and (l,-m) modes occupy opposite signs of f, the 'symmetric' window was shifted by one bin between the members of each pair -- exactly in the taper band containing the |h(f)| maximum -- violating the conjugate-pair identity h_{l,-m}(f) = (-1)^l conj(h_{lm}(-f)) for every ChooseFDModes-consumed model (IMRPhenomXHM/XPHM/XPNR/XO4a/HM), e.g. 1.02% (2,±2) amplitude asymmetry for an exactly-symmetric aligned-spin XHM input. Build the window on the true frequency grid (exactly even in f), and zero the surviving -fNyq bin when the resize truncates its +fNyq partner, so all conditioning commutes with complex conjugation. Verification (nonprec control, M=80, fmin=fref=20, deltaF=1/16): - mode-level TD conjugate-pair residual: 1.02e-2 -> 2.5e-16 (XHM); XPNR and HM machine-exact; XPHM reduced to its intrinsic model asymmetry; precessing XPHM amplitude identity at 4.2e-16. - U/V diagnostics (test/waveforms/uv_parity_diagnostics.py): D2 1.04e-9 -> 1.28e-17, D3 1.7e-9 -> 1.4e-16. (The corrupted band lay mostly below fmin, so U/V leakage was ~1e-9, not the mode-level 1e-2; the fix removes the interface floor entirely.) - Conditioned waveform change: 1.1e-2 full-band, 99% of it below fmin in the discarded turn-on band; 1.6e-3 in-band. Co-Authored-By: Claude Fable 5 --- .../Code/RIFT/lalsimutils.py | 33 +++++++++++++++++-- 1 file changed, 30 insertions(+), 3 deletions(-) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py index fd1859d4a..db5f791fe 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py @@ -3275,7 +3275,11 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil # but that is very difficult to do because modes of different 'm' and thus typical frequency generally mix # Note also that unless the segment length is large, this is often surprisingly few frequency bins for tapering if not(no_condition): - our_fvals = evaluate_fvals(hlmsdict[(2,2)]) + # lal_convention=True is REQUIRED here: ChooseFDModes returns the ascending + # center-DC grid, and the default (RIFT reversed) convention would shift the + # high-pass window by one bin between (l,m) and (l,-m), breaking the + # conjugate-pair identity at the percent level. See evaluate_fvals docs. + our_fvals = evaluate_fvals(hlmsdict[(2,2)], lal_convention=True) vectaper_symmetric = np.ones(len(our_fvals)) indx_below = np.logical_and(np.abs(our_fvals)=P.fmin*fd_standoff_factor) vectaper_symmetric[indx_below] = 0.5 + 0.5*np.cos(np.pi* (np.abs(our_fvals[indx_below])/P.fmin - 1)/(1-fd_standoff_factor)) @@ -3300,7 +3304,15 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil indx_crit = TDlen - ntaper for mode in hlmsdict: + npts_fd_pre_resize = hlmsdict[mode].data.length hlmsdict[mode] = lal.ResizeCOMPLEX16FrequencySeries(hlmsdict[mode],0, TDlen) + if npts_fd_pre_resize > TDlen: + # The resize above truncated the +fNyq bin from the two-sided grid but kept + # its -fNyq partner (index 0). Zero it so the truncation commutes with the + # conjugate-pair (f -> -f) reflection for models with support at Nyquist. + # Tied to the truncation itself, NOT to no_condition: an asymmetric + # truncation would reintroduce (l,±m) asymmetry even on the raw path. + hlmsdict[mode].data.data[0] = 0 hlmsT[mode] = DataInverseFourier(hlmsdict[mode]) # Phase factors: see crazy conventions in https://git.ligo.org/lscsoft/lalsuite/-/blob/master/lalsimulation/lib/LALSimInspiral.c if True: #P.approx == lalIMRPhenomXHM or P.approx == lalIMRPhenomHM: @@ -4901,13 +4913,26 @@ def psd_windowing_factor(window_shape, TDlen): def evaluate_tvals(lal_tseries): return float(lal_tseries.epoch) +lal_tseries.deltaT*np.arange(lal_tseries.data.length) -def evaluate_fvals(lal_2sided_fseries): +def evaluate_fvals(lal_2sided_fseries, lal_convention=False): r""" - evaluate_fvals(lal_2sided_fseries) + evaluate_fvals(lal_2sided_fseries, lal_convention=False) Associates frequencies with a 2sided lal complex array. Compare with 'self.longweights' code Done by HAND in PrecessingOrbitModesOfFrequency Manually *reverses* convention re sign of \omega used in lal! + lal_convention=False (default): RIFT's own two-sided packing (arrays produced by + DataFourier/complex_hoff): even length, REVERSED, f[k] = deltaF*(npts/2 - k). + WARNING: for odd-length arrays this assigns f on a grid offset by deltaF/2 -- + it is only correct for the even-length RIFT packing. + lal_convention=True: the packing of two-sided series returned directly by LAL + generators (e.g. SimInspiralChooseFDModes): ASCENDING [-fNyq, ..., 0, ..., +fNyq] + with DC at index npts//2 (odd length TDlen+1; even length after truncation keeps + DC at npts//2). f[k] = deltaF*(k - npts//2), exact for both parities of npts. + Use this - not the default - on ChooseFDModes output: the default's reversal + + half-bin offset shifts any |f|-symmetric window by one bin between (l,m) and + (l,-m) modes (which live at opposite signs of f), breaking the conjugate-pair + identity h_{l,-m}(f) = (-1)^l conj(h_{lm}(-f)) at the percent level. + Notes: a) XXXFrequencySeries have an f0 and a deltaF, and *logically* they should run from f0....f0+N*df b) I will always use COMPLEX16FrequencySeries, which should run from -fNyq...fNyq-df @@ -4939,6 +4964,8 @@ def evaluate_fvals(lal_2sided_fseries): """ npts = lal_2sided_fseries.data.length df = lal_2sided_fseries.deltaF + if lal_convention: + return df*(np.arange(npts) - npts//2) fvals = np.zeros(npts) # https://www.lsc-group.phys.uwm.edu/daswg/projects/lal/nightly/docs/html/group___time_freq_f_f_t__h.html # https://www.lsc-group.phys.uwm.edu/daswg/projects/lal/nightly/docs/html/_time_freq_f_f_t_8h.html