Fix AV network sky-frame inverse mapping

.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

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

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.

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