diff --git a/.travis/test-core-units.sh b/.travis/test-core-units.sh index 1114c955d..6f4733768 100755 --- a/.travis/test-core-units.sh +++ b/.travis/test-core-units.sh @@ -81,6 +81,8 @@ FILES=( "$C/test/test_srate_resample_time_marginalization.py" "$C/test/test_vectorized_lal_tools_split.py" "$C/test/test_noloop_accumulator_shapes.py" + # -- detector-network sky coordinates: direction and inverse round trip + "$C/test/test_sky_rotations.py" # The Bilby-convention noise evidence: its normalization is checked against the # 4/T sum |d|^2/S formula through the real ComplexIP, so a factor of two regresses # loudly instead of shifting every reported log evidence by a plausible amount. @@ -317,8 +319,8 @@ done # direction that matters: 350 >= 347 passes today, and if pytest-subtests ever leaves the # 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=590 +# The detector-network sky mapping adds two passing tests and no skips. +EXPECTED_TESTS=592 # 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 13 skips are # environment legs -- cupy in test_seeding_reproducibility, device legs in @@ -327,7 +329,7 @@ EXPECTED_TESTS=590 # (mcsamplerNFlow is an optional dependency and is absent from the IGWN environment), and # the xfail in test_uv_symmetry. test_eos_portfolio_sampler.py adds 12 tests and # test_cip_portfolio_members.py 4, none of them skips. -EXPECTED_PASSED=577 +EXPECTED_PASSED=579 MAX_SKIPPED=13 # The floors must be INTEGERS, and this is checked rather than assumed. `[ 347 -lt FOO ]` does diff --git a/MonteCarloMarginalizeCode/Code/RIFT/misc/sky_rotations.py b/MonteCarloMarginalizeCode/Code/RIFT/misc/sky_rotations.py index 4459c4480..897f683f9 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/misc/sky_rotations.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/misc/sky_rotations.py @@ -37,6 +37,20 @@ def assign_sky_frame(det0,det1,theEpochFiducial): frmInverse= np.asarray(np.matrix(frm).I) # Create an orthonormal frame to undo the transform above +def physical_to_network(theta, phi, frame=None, xpy=np): + """Map physical equatorial angles into the assigned network frame.""" + if frame is None: + frame = frm + return lalsimutils.polar_angles_in_frame_alt(frame, theta, phi, xpy=xpy) + + +def network_to_physical(theta, phi, inverse_frame=None, xpy=np): + """Map sampled network-frame angles back to physical equatorial angles.""" + if inverse_frame is None: + inverse_frame = frmInverse + return lalsimutils.polar_angles_in_frame_alt(inverse_frame, theta, phi, xpy=xpy) + + # USE INSTEAD # functools(lalsimutils.polar_angles_in_frame_alt,frm) def rotate_sky_forwards_scalar(theta,phi,frm=frm): # When theta=0 we are describing the coordinats of the zhat direction in the vecZ frame diff --git a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode index 6b9667527..cf175e6a7 100755 --- a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode +++ b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode @@ -2245,9 +2245,13 @@ if opts.internal_sky_network_coordinates: opts.internal_sky_network_coordinates = False # stop using this code else: sky_rotations.assign_sky_frame(ifo_list[0], ifo_list[1], fiducial_epoch) - frm = identity_convert_togpu(sky_rotations.frm) - my_rotation = functools.partial(lalsimutils.polar_angles_in_frame_alt,frm,xpy=xpy_default) - my_rotation_cpu = functools.partial(lalsimutils.polar_angles_in_frame_alt,sky_rotations.frm,xpy=np) + # The integration variables are angles in the network-aligned frame. The + # likelihood and output XML require physical equatorial angles, so apply + # the inverse frame here. Using ``frm`` would instead apply the + # physical->network transform a second time. + frm_inverse = identity_convert_togpu(sky_rotations.frmInverse) + my_rotation = functools.partial(sky_rotations.network_to_physical,inverse_frame=frm_inverse,xpy=xpy_default) + my_rotation_cpu = functools.partial(sky_rotations.network_to_physical,inverse_frame=sky_rotations.frmInverse,xpy=np) if opts.internal_sky_network_coordinates and (opts.limit_right_ascension or opts.limit_declination): # The sampled RA/dec live in the network-aligned frame, so a truth-centered sky box given in # equatorial coordinates would silently select the wrong patch of sky. Fail loudly. diff --git a/MonteCarloMarginalizeCode/Code/test/test_sky_rotations.py b/MonteCarloMarginalizeCode/Code/test/test_sky_rotations.py new file mode 100644 index 000000000..875c3abff --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_sky_rotations.py @@ -0,0 +1,45 @@ +"""Regression tests for the detector-network sky coordinate frame.""" + +import numpy as np +import lal + +from RIFT import lalsimutils +from RIFT.misc import sky_rotations + + +def _angular_separation(theta_a, phi_a, theta_b, phi_b): + """Return the angle between two directions without RA wrap ambiguity.""" + a = lalsimutils.nhat(theta_a, phi_a) + b = lalsimutils.nhat(theta_b, phi_b) + return np.arccos(np.clip(np.dot(a, b), -1.0, 1.0)) + + +def test_physical_to_network_polar_angle_is_baseline_projection(): + """Network-frame cos(theta) is n dot the detector baseline at the epoch.""" + epoch = lal.LIGOTimeGPS(1389783918.0) + sky_rotations.assign_sky_frame("H1", "L1", epoch) + + theta = np.deg2rad(90.0 - 42.0195) + phi = np.deg2rad(167.5254) + theta_network, _ = sky_rotations.physical_to_network(theta, phi, xpy=np) + + expected = np.dot(lalsimutils.nhat(theta, phi), sky_rotations.vecZnew) + assert np.isclose(np.cos(theta_network), expected, rtol=0.0, atol=1e-12) + + +def test_network_to_physical_round_trip_uses_inverse_frame(): + """The inverse frame maps sampled network angles back to equatorial sky.""" + epoch = lal.LIGOTimeGPS(1389783918.0) + sky_rotations.assign_sky_frame("H1", "L1", epoch) + + for dec_deg, ra_deg in ((42.0195, 167.5254), (14.7999, 135.8866)): + theta = np.deg2rad(90.0 - dec_deg) + phi = np.deg2rad(ra_deg) + theta_network, phi_network = sky_rotations.physical_to_network( + theta, phi, xpy=np + ) + theta_out, phi_out = sky_rotations.network_to_physical( + theta_network, phi_network, xpy=np + ) + + assert _angular_separation(theta, phi, theta_out, phi_out) < 1e-7