From 380c9548a2943ef385a325cdfa467fd94b9b55dc Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Sun, 4 Oct 2026 16:43:12 -0700 Subject: [PATCH] lalsimutils: convert_waveform_coordinates applies source_redshift to input masses With source_redshift=z, the vectorized branches computed mu1/mu2 from the source-frame mc. On rift_O4c the final rescale also scaled a ROW instead of a column (fixed on rift_O4d by da75550b8). Now the low-level mass inputs (mc, mc_ecc, m1, m2, mtot) are multiplied by 1+z once, on a copy, at entry. The per-row fallback rescales P.m1, P.m2 only when no input was rescaled. Test: vectorized output vs the per-row extract_param path at z=0 and 0.3. Co-Authored-By: Claude Opus 5.5 --- .travis/test-core-units.sh | 8 +- .../Code/RIFT/lalsimutils.py | 23 +++-- ...est_convert_coordinates_source_redshift.py | 84 +++++++++++++++++++ 3 files changed, 104 insertions(+), 11 deletions(-) create mode 100644 MonteCarloMarginalizeCode/Code/test/test_convert_coordinates_source_redshift.py diff --git a/.travis/test-core-units.sh b/.travis/test-core-units.sh index 2355f30be..c25d3b35b 100755 --- a/.travis/test-core-units.sh +++ b/.travis/test-core-units.sh @@ -111,6 +111,8 @@ FILES=( "$C/test/test_lisa_ini_contract.py" "$C/test/test_tracer_placement_gp.py" "$C/test/waveforms/test_uv_symmetry.py" + # -- coordinate conversion (source-frame CIP input -> detector-frame fit coordinates) + "$C/test/test_convert_coordinates_source_redshift.py" ) # A manifest entry that stops existing is a SILENT no-op: the gate keeps passing while @@ -163,6 +165,8 @@ done # fed odeint's float probe to len(x) pdfs; the t_ref wiring in all three ILE drivers) # 378/366 + test_response_order.py (8 tests: SNR tightening, Halton independence, # compound-axis semantics, reference resolution/tail charging, and bank preflight) +# 393/381 + test_convert_coordinates_source_redshift.py (13 tests: vectorized +# convert_waveform_coordinates vs per-row extract_param at source_redshift 0 and 0.3) # # RAISE these when files are added: a floor left at the old value passes while covering less, # which is the failure this gate exists to catch. @@ -178,12 +182,12 @@ done # runner's closure the count falls back to 347 and still passes. Pinning 350 would turn an # unrelated dependency change into a red gate. # Two XML/grid template-finalization regressions, with no added skips. -EXPECTED_TESTS=380 +EXPECTED_TESTS=393 # Outcomes, not just exit status: a collection floor cannot see a test that collects, runs and # asserts nothing, and a pytest.skip can quietly absorb a lost gate. The 12 skips are # environment legs -- cupy in test_seeding_reproducibility, device legs in # test_dslice_device_native, and the xfail in test_uv_symmetry. -EXPECTED_PASSED=368 +EXPECTED_PASSED=381 MAX_SKIPPED=12 # The floors must be INTEGERS, and this is checked rather than assumed. `[ 347 -lt FOO ]` does diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py index 0c2f48111..b9035b890 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py @@ -5919,6 +5919,17 @@ def convert_waveform_coordinates(x_in,coord_names=['mc', 'eta'],low_level_coord_ - source_redshift: if nonzero, convert m1 -> m1 (1+z)=m_z, as fit is done in the detector frame. We are **assuming source-frame sampling** """ x_out = np.zeros( (len(x_in), len(coord_names) ) ) + # Source-frame input: move the mass coordinates to the detector frame ONCE, here, so every vectorized branch below + # and the per-row fallback see detector-frame masses. Other mass-dimensional inputs (e.g. mu1, mu2) do not scale + # by (1+z); for those, the per-row fallback applies the redshift to P.m1, P.m2 instead. + redshift_applied_to_input = False + if source_redshift: + mass_names_in = [p for p in ['mc', 'mc_ecc', 'm1', 'm2', 'mtot'] if p in low_level_coord_names] + if mass_names_in: + x_in = np.array(x_in, dtype=float) # copy: never rescale the caller's array + for p in mass_names_in: + x_in[:, low_level_coord_names.index(p)] *= (1+source_redshift) + redshift_applied_to_input = True # Check for trivial identity transformations and do those by direct copy, then remove those from the list of output coord names coord_names_reduced = coord_names.copy() for p in low_level_coord_names: @@ -6341,13 +6352,6 @@ def convert_waveform_coordinates(x_in,coord_names=['mc', 'eta'],low_level_coord_ coord_names_reduced.remove('DeltaLambdaTilde') - # perform any mass conversions needed, so output in detector frame given input in source frame - if source_redshift: - for name in ['mc', 'm1', 'm2']: - if name in coord_names: - indx_name = coord_names.index(name) - x_out[:,indx_name] *= (1+source_redshift) - # return if we don't need to do any more conversions (e.g., if we only have --parameter specification) if len(coord_names_reduced)<1: return x_out @@ -6361,8 +6365,9 @@ def convert_waveform_coordinates(x_in,coord_names=['mc', 'eta'],low_level_coord_ if low_level_coord_names[indx] != 'chi_pavg': P.assign_param( low_level_coord_names[indx], x_in[indx_out,indx]) # Apply redshift: assume input is source-frame mass, convert m1 -> m1(1+z) = m1_z, as fit used detector frame - P.m1 = P.m1*(1+source_redshift) - P.m2 = P.m2*(1+source_redshift) + if source_redshift and not redshift_applied_to_input: + P.m1 = P.m1*(1+source_redshift) + P.m2 = P.m2*(1+source_redshift) for indx in np.arange(len(coord_names_reduced)): p = coord_names_reduced[indx] indx_p_out= coord_names.index(p) diff --git a/MonteCarloMarginalizeCode/Code/test/test_convert_coordinates_source_redshift.py b/MonteCarloMarginalizeCode/Code/test/test_convert_coordinates_source_redshift.py new file mode 100644 index 000000000..76551c149 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_convert_coordinates_source_redshift.py @@ -0,0 +1,84 @@ +"""convert_waveform_coordinates(..., source_redshift=z) must return detector-frame values. + +Reference: the per-row ChooseWaveformParams path the function itself falls back to +(assign source-frame low-level coordinates, scale P.m1, P.m2 by 1+z, extract_param). +Masses are in Msun, as CIP passes them. +""" +import numpy as np +import pytest + +from RIFT import lalsimutils as lsu + +N = 50 + + +def _draw(low_level_coord_names, rng): + ranges = { + 'mc': (5., 40.), 'delta_mc': (0.05, 0.85), 'eta': (0.1, 0.249), + 's1z': (-0.9, 0.9), 's2z': (-0.9, 0.9), + 'chi1': (0., 0.99), 'chi2': (0., 0.99), + 'cos_theta1': (-1., 1.), 'cos_theta2': (-1., 1.), + 'phi1': (0., 2*np.pi), 'phi2': (0., 2*np.pi), + 's1z_bar': (-0.9, 0.9), 's2z_bar': (-0.9, 0.9), + 'chi1_perp_bar': (0., 0.9), 'chi2_perp_bar': (0., 0.9), + 'lambda1': (0., 1000.), 'lambda2': (0., 1000.), + } + return np.column_stack([rng.uniform(*ranges[p], N) for p in low_level_coord_names]) + + +def _per_row(x_in, coord_names, low_level_coord_names, z): + out = np.zeros((len(x_in), len(coord_names))) + for i, row in enumerate(x_in): + P = lsu.ChooseWaveformParams() + for p, val in zip(low_level_coord_names, row): + P.assign_param(p, val) + P.m1 *= 1 + z + P.m2 *= 1 + z + out[i] = [P.extract_param(p) for p in coord_names] + return out + + +CASES = [ + # spherical spins: vectorized mu1, mu2 (the case PR #204's RF features use) + (['mc', 'delta_mc', 'mu1', 'mu2', 'xi', 'chiMinus', 's1x', 's1y', 's2x', 's2y'], + ['mc', 'delta_mc', 'chi1', 'cos_theta1', 'phi1', 'chi2', 'cos_theta2', 'phi2']), + # aligned cartesian: vectorized mu1, mu2 + (['mu1', 'mu2', 'delta_mc', 'chiMinus', 'mc'], ['mc', 'delta_mc', 's1z', 's2z']), + # vectorized m1, m2 from mc, eta + (['mc', 'm1', 'm2', 'eta', 'xi'], ['mc', 'eta', 's1z', 's2z']), + # pseudo-cylindrical spins: vectorized mu1, mu2. s*z_bar precedes chi*_perp_bar because + # assign_param('chi1_perp_bar') reads the current s1z; the reverse order builds a different spin. + (['mu1', 'mu2', 'delta_mc', 'xi', 'chiMinus', 'chi_p'], + ['mc', 'delta_mc', 's1z_bar', 'chi1_perp_bar', 'phi1', 's2z_bar', 'chi2_perp_bar', 'phi2']), + # tidal: vectorized mu1, mu2 next to LambdaTilde + (['mu1', 'mu2', 'delta_mc', 'LambdaTilde', 'DeltaLambdaTilde'], + ['mc', 'delta_mc', 's1z', 's2z', 'lambda1', 'lambda2']), + # mtot and q go through the per-row fallback: the redshift must be applied once, not twice + (['mc', 'mtot', 'q', 'm1'], ['mc', 'delta_mc', 's1z', 's2z']), +] + + +@pytest.mark.parametrize("z", [0.0, 0.3]) +@pytest.mark.parametrize("coord_names,low_level_coord_names", CASES) +def test_vectorized_matches_per_row(coord_names, low_level_coord_names, z): + rng = np.random.default_rng(20261004) + x_in = _draw(low_level_coord_names, rng) + x_in_before = x_in.copy() + got = lsu.convert_waveform_coordinates(x_in, coord_names=coord_names, + low_level_coord_names=low_level_coord_names, + source_redshift=z) + ref = _per_row(x_in, coord_names, low_level_coord_names, z) + np.testing.assert_array_equal(x_in, x_in_before) # caller's array is not rescaled in place + np.testing.assert_allclose(got, ref, rtol=1e-9, atol=1e-12, + err_msg=f"z={z} coords={coord_names}") + + +def test_mass_columns_scale_on_every_row(): + rng = np.random.default_rng(1) + low = ['mc', 'eta', 's1z', 's2z'] + x_in = _draw(low, rng) + names = ['mc', 'm1', 'm2'] + z0 = lsu.convert_waveform_coordinates(x_in, coord_names=names, low_level_coord_names=low) + z1 = lsu.convert_waveform_coordinates(x_in, coord_names=names, low_level_coord_names=low, + source_redshift=0.3) + np.testing.assert_allclose(z1, 1.3*z0, rtol=1e-12)