lalsimutils: phi12, chi_p_vec, vectorized in-plane spin coordinates (O4d port of #203)

.travis/test-core-units.sh

@@ -160,6 +160,8 @@ FILES=(
   "$C/test/test_tracer_placement_gp.py"
   "$C/test/waveforms/test_uv_symmetry.py"
   "$C/test/test_complex_overlap_interpolate_max.py"
+  # -- coordinates: vectorized in-plane spin / ring coordinates agree with extract_param
+  "$C/test/test_ring_coordinates.py"
 )
 
 # A manifest entry that stops existing is a SILENT no-op: the gate keeps passing while
@@ -313,6 +315,13 @@ done
 #            RIFT_COREUNIT_PYTHON pointed at the IGWN interpreter: per-file 600 over 52 files,
 #            junit 603 collected / 590 passed / 13 skipped / 0 failed.
 #
+#   613/600  test_ring_coordinates.py added (7 tests: vectorized in-plane spin and ring
+#            coordinates, phi12, chi_p_vec, object-array input, Kerr rule), merged with
+#            test_complex_overlap_interpolate_max.py above.  MEASURED on CIT (ldas-grid; numpy
+#            backend) 2026-09-30, IGWN conda python 3.11 / lal 7.7.0, with RIFT_COREUNIT_PYTHON
+#            pointed at the IGWN interpreter: junit 616 collected / 603 passed / 13 skipped /
+#            0 failed, of which 3 are subtests.
+#
 # 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.
 # DO NOT RAISE THESE TO THE RUNNER'S NUMBERS.  The GitHub runner reports 350 collected / 338
@@ -327,8 +336,14 @@ 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.
 # The detector-network sky mapping adds two passing tests and no skips.
-# test_complex_overlap_interpolate_max.py adds 8 passing tests and no skips.
-EXPECTED_TESTS=600
+# test_ring_coordinates.py adds two more, no skips.  RE-MEASURED with it on CIT (ldas-grid;
+# `import cupy` FAILS there, numpy backend) 2026-09-30, IGWN conda python 3.11 / lal 7.7.0:
+# per-file 600, junit 603 collected / 590 passed / 13 skipped / 0 failed, of which 3 are
+# subtests.  So the plugin-free floors are 600/587; the old 592/579 sat 6 below the tree.
+# Review of #377 added five more ring-coordinate tests, no skips (605/592 measured).
+# test_complex_overlap_interpolate_max.py (#375) adds 8 passing tests and no skips.
+# Merged with #375: 613/600 (see the table above).
+EXPECTED_TESTS=613
 # 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
@@ -337,7 +352,7 @@ EXPECTED_TESTS=600
 # (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=587
+EXPECTED_PASSED=600
 MAX_SKIPPED=13
 
 # The floors must be INTEGERS, and this is checked rather than assumed.  `[ 347 -lt FOO ]` does

MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py

@@ -418,7 +418,7 @@ def lsu_StringFromPNOrder(order):
 
 # so far, used for puffball, to prevent insanity (infinite growth) and/or death to downselect
 #   - note we also provide for extrinsic: RA (phi), phiref, psi, just in case we need it in the future
-periodic_params = {'phi1':2*np.pi, 'phi2':2*np.pi, 'phiref':2*np.pi, 'psi':np.pi, 'meanPerAno':2*np.pi, 'phi':2*np.pi, 'phiJL':2*np.pi, 'psiJ':2*np.pi}
+periodic_params = {'phi1':2*np.pi, 'phi2':2*np.pi, 'phi12':2*np.pi, 'phiref':2*np.pi, 'psi':np.pi, 'meanPerAno':2*np.pi, 'phi':2*np.pi, 'phiJL':2*np.pi, 'psiJ':2*np.pi}
 
 tex_dictionary  = {
  "mtot": r'$M$',
@@ -441,6 +441,7 @@ def lsu_StringFromPNOrder(order):
   "DeltaOverM2_perp" : r"$\Delta_\perp$",
   "DeltaOverM2_L" : r"$\Delta_{||}$",
   "SOverM2_perp" : r"$S_\perp$",
+  "chi_p_vec" : r"$\chi_{p,{\rm vec}}$",
   "SOverM2_L" : r"$S_{||}$",
   "eta": r"$\eta$",
   "chi_eff": r"$\chi_{eff}$",
@@ -1390,6 +1391,20 @@ def extract_param(self,p):
             S2p = (m2**2 * chi2)[:2]
             Sp = np.max([np.linalg.norm( A1*S1p), np.linalg.norm(A2*S2p)])
             return Sp/(A1*m1**2)  # divide by term for *larger* BH
+        if p == 'phi12':
+            # azimuth of spin 2's in-plane component relative to spin 1's, in [0, 2 pi), L frame.
+            # Undefined if either in-plane component vanishes; 0 is returned then (same as the vectorized path).
+            if np.hypot(self.s1x, self.s1y) == 0 or np.hypot(self.s2x, self.s2y) == 0:
+                return 0.
+            val = np.mod(np.arctan2(self.s2y, self.s2x) - np.arctan2(self.s1y, self.s1x), 2*np.pi)
+            return 0. if val >= 2*np.pi else val   # np.mod of a tiny negative number rounds to 2 pi
+        if p == 'chi_p_vec':
+            # vector-sum (ring) analogue of chi_p: |A1 S1perp + A2 S2perp| / (A1 m1^2), same A1, A2 as chi_p.
+            # chi_p keeps the larger of the two terms; this keeps their vector sum, so it depends on phi12
+            # (in-plane spins that cancel give a small value).  L frame.
+            q = self.m2/self.m1
+            A1 = (2+ 3.*q/2); A2 = (2+3./(2*q))
+            return np.abs( (self.s1x + 1j*self.s1y) + (A2/A1)*q**2*(self.s2x + 1j*self.s2y) )
         if p == 'chi_pavg':
             if (abs(self.s1x) < 1e-4 and abs(self.s1y) < 1e-4 and abs(self.s2x) < 1e-4 and abs(self.s2y) < 1e-4):
                 chipavg = 0.0
@@ -5918,6 +5933,7 @@ 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) ) )
+    kerr_violation_ring = None   # set by the vectorized in-plane block, which can end the conversion early
     # 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:
@@ -6131,6 +6147,35 @@ def convert_waveform_coordinates(x_in,coord_names=['mc', 'eta'],low_level_coord_
             x_out[:,indx_pout_s2y] = x_in[:,indx_chi2]*sintheta2*sinphi2
             coord_names_reduced.remove('s2x')
             coord_names_reduced.remove('s2y')
+        # in-plane magnitudes, relative azimuth, and ring coordinates, vectorized.  L frame only:
+        # for any other spin_convention these fall through to extract_param, as before.
+        ring_names = ['chi1_perp', 'chi2_perp', 'phi12', 'SOverM2_perp', 'DeltaOverM2_perp', 'chi_p_vec']
+        if spin_convention == "L" and any(p in coord_names_reduced for p in ring_names):
+            # CIP's default sampler passes an object array of python floats; ufuncs need a float array
+            xf = np.asarray(x_in, dtype=float)
+            m1f = np.asarray(m1_vals, dtype=float); m2f = np.asarray(m2_vals, dtype=float)
+            indx_phi1 = low_level_coord_names.index('phi1')
+            indx_phi2 = low_level_coord_names.index('phi2')
+            chi1_perp = xf[:,indx_chi1]*np.sqrt(1-xf[:,indx_ct1]**2)
+            chi2_perp = xf[:,indx_chi2]*np.sqrt(1-xf[:,indx_ct2]**2)
+            v1 = chi1_perp*np.exp(1j*xf[:,indx_phi1])
+            v2 = chi2_perp*np.exp(1j*xf[:,indx_phi2])
+            mtot_vals = m1f + m2f
+            q_vals = m2f/m1f
+            A1 = 2 + 1.5*q_vals; A2 = 2 + 1.5/q_vals
+            phi12 = np.where((chi1_perp > 0) & (chi2_perp > 0), np.mod(xf[:,indx_phi2] - xf[:,indx_phi1], 2*np.pi), 0.)
+            ring_vals = {'chi1_perp': chi1_perp, 'chi2_perp': chi2_perp,
+                         'phi12': np.where(phi12 >= 2*np.pi, 0., phi12),
+                         'SOverM2_perp': np.abs(v1*m1f**2 + v2*m2f**2)/mtot_vals**2,
+                         'DeltaOverM2_perp': np.abs(v1*m1f - v2*m2f)/mtot_vals,
+                         'chi_p_vec': np.abs(v1 + (A2/A1)*q_vals**2*v2)}
+            for p in ring_names:
+                if p in coord_names_reduced:
+                    x_out[:,coord_names.index(p)] = ring_vals[p]
+                    coord_names_reduced.remove(p)
+            if enforce_kerr:
+                # same rule as the per-row fallthrough below, which this block can bypass
+                kerr_violation_ring = (xf[:,indx_chi1] > 1) | (xf[:,indx_chi2] > 1)
 
     # Spin pseudo-cylindrical coordinate names, standard framing
     if  ('s1z_bar' in low_level_coord_names) and ('phi1' in low_level_coord_names)  and ('s2z_bar' in low_level_coord_names) and ('phi2' in low_level_coord_names) and ('mc' in low_level_coord_names) and ('eta' in low_level_coord_names or 'delta_mc' in low_level_coord_names):
@@ -6348,6 +6393,8 @@ def convert_waveform_coordinates(x_in,coord_names=['mc', 'eta'],low_level_coord_
                 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 kerr_violation_ring is not None:
+        x_out[kerr_violation_ring] = -np.inf
     if len(coord_names_reduced)<1:
         return x_out
 

MonteCarloMarginalizeCode/Code/test/test_ring_coordinates.py

@@ -0,0 +1,101 @@
+"""phi12 and the ring coordinates (chi1_perp, chi2_perp, SOverM2_perp, DeltaOverM2_perp, chi_p_vec).
+
+The vectorized spherical-spin path of convert_waveform_coordinates must agree with
+ChooseWaveformParams.extract_param, and must not fall through to the per-row loop.
+"""
+
+import numpy as np
+import lal
+
+from RIFT import lalsimutils
+
+RING = ['chi1_perp', 'chi2_perp', 'phi12', 'SOverM2_perp', 'DeltaOverM2_perp', 'chi_p_vec']
+LOW = ['mc', 'delta_mc', 'chi1', 'cos_theta1', 'phi1', 'chi2', 'cos_theta2', 'phi2']
+
+
+def _draws(n=300, seed=4):
+    rng = np.random.default_rng(seed)
+    return np.column_stack([rng.uniform(5, 30, n), rng.uniform(0.01, 0.8, n), rng.uniform(0, 0.99, n),
+                            rng.uniform(-1, 1, n), rng.uniform(0, 2 * np.pi, n), rng.uniform(0, 0.99, n),
+                            rng.uniform(-1, 1, n), rng.uniform(0, 2 * np.pi, n)])
+
+
+def _params(row):
+    mc, dmc, c1, ct1, p1, c2, ct2, p2 = row
+    m1, m2 = lalsimutils.m1m2(mc, 0.25 * (1 - dmc ** 2))
+    P = lalsimutils.ChooseWaveformParams()
+    P.m1, P.m2 = m1 * lal.MSUN_SI, m2 * lal.MSUN_SI
+    s1 = c1 * np.sqrt(1 - ct1 ** 2)
+    s2 = c2 * np.sqrt(1 - ct2 ** 2)
+    P.s1x, P.s1y, P.s1z = s1 * np.cos(p1), s1 * np.sin(p1), c1 * ct1
+    P.s2x, P.s2y, P.s2z = s2 * np.cos(p2), s2 * np.sin(p2), c2 * ct2
+    return P
+
+
+def test_vectorized_matches_extract_param(capsys):
+    x = _draws()
+    y = lalsimutils.convert_waveform_coordinates(x, coord_names=RING, low_level_coord_names=LOW)
+    assert "Fallthrough" not in capsys.readouterr().out
+    for i, row in enumerate(x):
+        P = _params(row)
+        for j, name in enumerate(RING):
+            ref = P.extract_param(name)
+            d = abs(y[i, j] - ref)
+            if name == 'phi12':
+                d = min(d, 2 * np.pi - d)
+            assert d < 1e-9, (name, y[i, j], ref)
+
+
+def test_chi_p_vec_limits():
+    # one in-plane spin: chi_p_vec equals chi1_perp.  Antiparallel in-plane spins of the weighted size cancel.
+    P = lalsimutils.ChooseWaveformParams()
+    P.m1, P.m2 = 10 * lal.MSUN_SI, 5 * lal.MSUN_SI
+    P.s1x, P.s1y, P.s2x, P.s2y = 0.4, 0.0, 0.0, 0.0
+    assert abs(P.extract_param('chi_p_vec') - 0.4) < 1e-12
+    q = 0.5
+    A1 = 2 + 1.5 * q
+    A2 = 2 + 1.5 / q
+    P.s2x = -0.4 / ((A2 / A1) * q ** 2)
+    assert P.extract_param('chi_p_vec') < 1e-12
+    assert abs(P.extract_param('phi12') - np.pi) < 1e-12
+
+
+def _one(*row):
+    return np.array([row], dtype=float)
+
+
+def test_phi12_sign_both_paths():
+    # spin 2 sixty degrees ahead of spin 1: phi12 = pi/3, not 5 pi/3 (a sign flip in both paths fails here)
+    row = (10., 0.3, 0.5, 0., 0.2, 0.5, 0., 0.2 + np.pi / 3)
+    y = lalsimutils.convert_waveform_coordinates(_one(*row), coord_names=['phi12'], low_level_coord_names=LOW)
+    assert abs(y[0, 0] - np.pi / 3) < 1e-12
+    assert abs(_params(row).extract_param('phi12') - np.pi / 3) < 1e-12
+
+
+def test_object_array_input():
+    # CIP's default sampler (adaptive_cartesian) passes an object array of python floats
+    x = _draws(20)
+    y = lalsimutils.convert_waveform_coordinates(x.astype(object), coord_names=RING, low_level_coord_names=LOW)
+    y_ref = lalsimutils.convert_waveform_coordinates(x, coord_names=RING, low_level_coord_names=LOW)
+    assert np.allclose(np.asarray(y, dtype=float), y_ref, atol=1e-12)
+
+
+def test_phi12_zero_inplane_spin():
+    # phi12 is undefined without an in-plane component; both paths return 0
+    for row in [(10., 0.3, 0.0, 0.5, 1.5, 0.5, 0.2, 2.5), (10., 0.3, 0.5, 0.2, 1.5, 0.5, -1.0, 2.5)]:
+        y = lalsimutils.convert_waveform_coordinates(_one(*row), coord_names=['phi12'], low_level_coord_names=LOW)
+        assert y[0, 0] == 0.
+        assert _params(row).extract_param('phi12') == 0.
+
+
+def test_enforce_kerr_in_ring_block():
+    # the vectorized block can end the conversion early; it applies the fallthrough's Kerr rule itself
+    x = np.vstack([_one(10., 0.3, 1.2, 0.2, 1.5, 0.5, 0.2, 2.5), _one(10., 0.3, 0.5, 0.2, 1.5, 0.5, 0.2, 2.5)])
+    y = lalsimutils.convert_waveform_coordinates(x, coord_names=['mc', 'delta_mc', 'chi1_perp', 'chi_p_vec'],
+                                                 low_level_coord_names=LOW, enforce_kerr=True)
+    assert np.all(y[0] == -np.inf) and np.all(np.isfinite(y[1]))
+
+
+def test_chi_p_vec_not_assignable():
+    # grid readers call assign_param on every column named in valid_params; chi_p_vec is derived only
+    assert 'chi_p_vec' not in lalsimutils.valid_params