TST: change all likelihood and waveform tests to use more recent waveforms (#791)

bilby/gw/source.py

@@ -17,6 +17,35 @@
 Unused waveform_kwargs: {waveform_kwargs}
 """
 
+# LAL error-code messages that indicate the requested parameters are outside
+# the model's domain of validity (e.g. a sampler proposing an unphysical
+# point) rather than a genuine bug. Different approximants surface an
+# out-of-domain point through different XLAL error codes:
+#   XLAL_EDOM   -> "Input domain error"            (e.g. IMRPhenomPv2/XP)
+#   XLAL_EFUNC  -> "Internal function call failed"  (e.g. IMRPhenomXPHM via
+#                                                    the frequency sequence)
+# When ``catch_waveform_errors`` is enabled these are converted to a -inf
+# likelihood instead of being raised.
+#
+# XLAL_EINVAL ("Invalid argument") is deliberately not included: LAL raises
+# it when, e.g., precessing spins are passed to an aligned-spin approximant,
+# which is a configuration error that should not be silently ignored.
+WAVEFORM_DOMAIN_ERROR_MESSAGES = (
+    "Input domain error",
+    "Internal function call failed",
+)
+
+
+def _is_waveform_domain_error(exception):
+    """Return True if ``exception`` looks like an out-of-domain waveform error.
+
+    See :data:`WAVEFORM_DOMAIN_ERROR_MESSAGES` for the recognised messages.
+    """
+    if not exception.args:
+        return False
+    message = str(exception.args[0])
+    return any(msg in message for msg in WAVEFORM_DOMAIN_ERROR_MESSAGES)
+
 
 def gwsignal_binary_black_hole(frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
                                phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, **kwargs):
@@ -217,12 +246,7 @@ def gwsignal_binary_black_hole(frequency_array, mass_1, mass_2, luminosity_dista
         if not catch_waveform_errors:
             raise
         else:
-            EDOM = (
-                "Internal function call failed: Input domain error" in e.args[0]
-            ) or "Input domain error" in e.args[
-                0
-            ]
-            if EDOM:
+            if _is_waveform_domain_error(e):
                 failed_parameters = dict(mass_1=mass_1, mass_2=mass_2,
                                          spin_1=(spin_1x, spin_1y, spin_1z),
                                          spin_2=(spin_2x, spin_2y, spin_2z),
@@ -645,8 +669,7 @@ def _base_lal_cbc_fd_waveform(
         if not catch_waveform_errors:
             raise
         else:
-            EDOM = (e.args[0] == 'Internal function call failed: Input domain error')
-            if EDOM:
+            if _is_waveform_domain_error(e):
                 failed_parameters = dict(mass_1=mass_1, mass_2=mass_2,
                                          spin_1=(spin_1x, spin_1y, spin_1z),
                                          spin_2=(spin_2x, spin_2y, spin_2z),
@@ -1131,8 +1154,7 @@ def _base_waveform_frequency_sequence(
         if not catch_waveform_errors:
             raise
         else:
-            EDOM = (e.args[0] == 'Internal function call failed: Input domain error')
-            if EDOM:
+            if _is_waveform_domain_error(e):
                 failed_parameters = dict(mass_1=mass_1, mass_2=mass_2,
                                          spin_1=(spin_1x, spin_1y, spin_1z),
                                          spin_2=(spin_2x, spin_2y, spin_2z),

pixi.toml

@@ -47,7 +47,7 @@ nestle-bilby = "*"
 zeus-bilby = "*"
 
 [feature.dynesty.pypi-dependencies]
-dynesty = {git = "https://github.com/joshspeagle/dynesty.git", branch = "master"}
+dynesty = { git = "https://github.com/joshspeagle/dynesty.git", branch = "master" }
 
 [feature.gw.dependencies]
 astropy = ">=5.0.0"
@@ -98,7 +98,7 @@ schwimmbad = "*"
 scikit-learn = "*"
 
 [feature.test.pypi-dependencies]
-bilby = {path = ".", editable = true}
+bilby = { path = ".", editable = true }
 
 [feature.torch.dependencies]
 pytorch = "*"

pyproject.toml

@@ -8,7 +8,7 @@ requires = [
 
 [project]
 authors = [
-  {name = "Bilby Developers"}
+  { name = "Bilby Developers" }
 ]
 classifiers = [
   "Development Status :: 5 - Production/Stable",
@@ -22,13 +22,13 @@ classifiers = [
 ]
 description = "A user-friendly Bayesian inference library"
 dynamic = ["dependencies", "optional-dependencies", "version"]
-license = {file = "LICENSE.md"}
+license = { file = "LICENSE.md" }
 maintainers = [
-  {name = "Greg Ashton"},
-  {name = "Sylvia Biscoveanu"},
-  {name = "Colm Talbot"},
-  {name = "Aditya Vijaykumar"},
-  {name = "Michael J. Williams"}
+  { name = "Greg Ashton" },
+  { name = "Sylvia Biscoveanu" },
+  { name = "Colm Talbot" },
+  { name = "Aditya Vijaykumar" },
+  { name = "Michael J. Williams" }
 ]
 name = "bilby"
 readme = "README.rst"
@@ -122,18 +122,18 @@ packages = [
 ]
 
 [tool.setuptools.dynamic]
-dependencies = {file = ["requirements.txt"]}
+dependencies = { file = ["requirements.txt"] }
 
 [tool.setuptools.dynamic.optional-dependencies]
-all = {file = [
+all = { file = [
   "gw_requirements.txt",
   "jax_requirements.txt",
   "mcmc_requirements.txt",
   "optional_requirements.txt"
-]}
-gw = {file = ["gw_requirements.txt"]}
-jax = {file = ["jax_requirements.txt"]}
-mcmc = {file = ["mcmc_requirements.txt"]}
+] }
+gw = { file = ["gw_requirements.txt"] }
+jax = { file = ["jax_requirements.txt"] }
+mcmc = { file = ["mcmc_requirements.txt"] }
 
 [tool.setuptools.package-data]
 "bilby.gw" = ["prior_files/*"]

test/gw/likelihood/roq_test.py

@@ -188,9 +188,9 @@ def test_create_roq_weights_with_params(self):
             priors=self.priors,
         )
         roq_llr = roq.log_likelihood_ratio(self.test_parameters)
-        self.assertEqual(
-            roq_llr,
-            self.roq.log_likelihood_ratio(self.test_parameters)
+        self.assertAlmostEqual(
+            float(roq_llr),
+            float(self.roq.log_likelihood_ratio(self.test_parameters)),
         )
         self.assertEqual(aac.get_namespace(roq_llr), self.xp)
 

test/gw/plot_test.py

@@ -21,7 +21,7 @@ def setUp(self):
                 frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
                 time_domain_source_model=None,
                 waveform_arguments=dict(
-                    reference_frequency=20.0, waveform_approximant="IMRPhenomPv2"
+                    reference_frequency=20.0, waveform_approximant="IMRPhenomXPHM"
                 ),
                 interferometers=dict(
                     H1=dict(optimal_SNR=1, parameters=injection_parameters),

test/gw/result_test.py

@@ -27,7 +27,7 @@ def setUp(self):
                 time_marginalization=True,
                 frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
                 waveform_arguments=dict(
-                    reference_frequency=20.0, waveform_approximant="IMRPhenomPv2"
+                    reference_frequency=20.0, waveform_approximant="IMRPhenomXPHM"
                 ),
                 interferometers=dict(
                     H1=dict(optimal_SNR=1, parameters=injection_parameters),

test/gw/source_test.py

@@ -8,6 +8,7 @@
 
 import numpy as np
 from copy import copy
+from parameterized import parameterized
 
 
 class TestLalBBH(unittest.TestCase):
@@ -26,7 +27,7 @@ def setUp(self):
             phase=0.0,
         )
         self.waveform_kwargs = dict(
-            waveform_approximant="IMRPhenomPv2",
+            waveform_approximant="IMRPhenomXPHM",
             reference_frequency=50.0,
             minimum_frequency=20.0,
             catch_waveform_errors=True,
@@ -208,7 +209,7 @@ def setUp(self):
             lambda_2=100.0,
         )
         self.waveform_kwargs = dict(
-            waveform_approximant="IMRPhenomPv2_NRTidal",
+            waveform_approximant="IMRPhenomPv2_NRTidalv2",
             reference_frequency=50.0,
             minimum_frequency=20.0,
         )
@@ -487,7 +488,7 @@ def setUp(self):
         self.full_frequencies_to_sequence = self.frequency_array >= self.minimum_frequency
         self.frequencies = self.frequency_array[self.full_frequencies_to_sequence]
         self.waveform_kwargs = dict(
-            waveform_approximant="IMRPhenomPv2_NRTidal",
+            waveform_approximant="IMRPhenomPv2_NRTidalv2",
             reference_frequency=50.0,
         )
 
@@ -546,14 +547,14 @@ def setUp(self):
             phase=0.0,
         )
         self.waveform_kwargs_fiducial = dict(
-            waveform_approximant="IMRPhenomPv2",
+            waveform_approximant="IMRPhenomXPHM",
             reference_frequency=50.0,
             minimum_frequency=20.0,
             catch_waveform_errors=True,
             fiducial=True,
         )
         self.waveform_kwargs_binned = dict(
-            waveform_approximant="IMRPhenomPv2",
+            waveform_approximant="IMRPhenomXPHM",
             reference_frequency=50.0,
             minimum_frequency=20.0,
             catch_waveform_errors=True,
@@ -682,13 +683,13 @@ def setUp(self):
             lambda_2=100.0,
         )
         self.waveform_kwargs_fiducial = dict(
-            waveform_approximant="IMRPhenomPv2_NRTidal",
+            waveform_approximant="IMRPhenomPv2_NRTidalv2",
             reference_frequency=50.0,
             minimum_frequency=20.0,
             fiducial=True,
         )
         self.waveform_kwargs_binned = dict(
-            waveform_approximant="IMRPhenomPv2_NRTidal",
+            waveform_approximant="IMRPhenomPv2_NRTidalv2",
             reference_frequency=50.0,
             minimum_frequency=20.0,
             fiducial=False,
@@ -769,5 +770,114 @@ def test_binned_fails_without_tidal_parameters(self):
             )
 
 
+PRECESSING_APPROXIMANTS = ["IMRPhenomPv2", "IMRPhenomXP", "IMRPhenomXPHM"]
+ALIGNED_APPROXIMANTS = ["IMRPhenomXAS", "IMRPhenomD", "IMRPhenomHM"]
+SOURCE_MODELS = ["lal", "frequency_sequence", "relbin_fiducial", "relbin_binned"]
+
+
+class TestWaveformErrorCatchingByApproximant(unittest.TestCase):
+    """
+    Check which LAL errors ``catch_waveform_errors`` converts to ``None``.
+
+    Out-of-domain parameters (here a negative mass) should be caught for all
+    approximants and source models. Precessing spins passed to an aligned-spin
+    approximant are a configuration error and should always be raised, even
+    with ``catch_waveform_errors=True``.
+    """
+
+    def setUp(self):
+        self.precessing_parameters = dict(
+            mass_1=30.0,
+            mass_2=30.0,
+            luminosity_distance=400.0,
+            a_1=0.4,
+            tilt_1=0.2,
+            phi_12=1.0,
+            a_2=0.8,
+            tilt_2=2.7,
+            phi_jl=2.9,
+            theta_jn=0.3,
+            phase=0.0,
+        )
+        self.aligned_parameters = copy(self.precessing_parameters)
+        self.aligned_parameters.update(tilt_1=0.0, tilt_2=0.0, phi_12=0.0, phi_jl=0.0)
+        self.frequency_array = bilby.core.utils.create_frequency_series(2048, 4)
+
+    def call_source(self, source_model, approximant, parameters, catch_waveform_errors=True):
+        kwargs = dict(
+            waveform_approximant=approximant,
+            reference_frequency=50.0,
+            catch_waveform_errors=catch_waveform_errors,
+            **parameters,
+        )
+        if source_model == "lal":
+            return bilby.gw.source.lal_binary_black_hole(
+                self.frequency_array, minimum_frequency=20.0, **kwargs
+            )
+        elif source_model == "frequency_sequence":
+            frequencies = self.frequency_array[self.frequency_array >= 20.0]
+            return bilby.gw.source.binary_black_hole_frequency_sequence(
+                self.frequency_array, frequencies=frequencies, **kwargs
+            )
+        elif source_model == "relbin_fiducial":
+            return bilby.gw.source.lal_binary_black_hole_relative_binning(
+                self.frequency_array, minimum_frequency=20.0, fiducial=True, **kwargs
+            )
+        elif source_model == "relbin_binned":
+            return bilby.gw.source.lal_binary_black_hole_relative_binning(
+                self.frequency_array,
+                minimum_frequency=20.0,
+                fiducial=False,
+                frequency_bin_edges=np.arange(20, 1500, 50),
+                **kwargs,
+            )
+        raise ValueError(f"Unknown source model {source_model}")
+
+    @parameterized.expand(
+        [(model, approx) for model in SOURCE_MODELS for approx in PRECESSING_APPROXIMANTS]
+    )
+    def test_precessing_negative_mass_is_caught(self, source_model, approximant):
+        parameters = copy(self.precessing_parameters)
+        parameters["mass_1"] = -30.0
+        self.assertIsNone(self.call_source(source_model, approximant, parameters))
+
+    @parameterized.expand(
+        [(model, approx) for model in SOURCE_MODELS for approx in ALIGNED_APPROXIMANTS]
+    )
+    def test_aligned_negative_mass_is_caught(self, source_model, approximant):
+        parameters = copy(self.aligned_parameters)
+        parameters["mass_1"] = -30.0
+        self.assertIsNone(self.call_source(source_model, approximant, parameters))
+
+    @parameterized.expand(
+        [(model, approx) for model in SOURCE_MODELS
+         for approx in PRECESSING_APPROXIMANTS + ALIGNED_APPROXIMANTS]
+    )
+    def test_negative_mass_raises_when_not_catching(self, source_model, approximant):
+        parameters = copy(self.aligned_parameters)
+        parameters["mass_1"] = -30.0
+        with self.assertRaises(RuntimeError):
+            self.call_source(
+                source_model, approximant, parameters, catch_waveform_errors=False
+            )
+
+    @parameterized.expand(
+        [(model, approx) for model in SOURCE_MODELS for approx in ALIGNED_APPROXIMANTS]
+    )
+    def test_aligned_approximant_with_precessing_spins_raises(self, source_model, approximant):
+        with self.assertRaises(RuntimeError) as context:
+            self.call_source(source_model, approximant, self.precessing_parameters)
+        self.assertIn("Invalid argument", str(context.exception))
+
+    @parameterized.expand(
+        [(model, approx) for model in SOURCE_MODELS for approx in PRECESSING_APPROXIMANTS]
+    )
+    def test_precessing_approximant_with_precessing_spins_works(self, source_model, approximant):
+        self.assertIsInstance(
+            self.call_source(source_model, approximant, self.precessing_parameters),
+            dict,
+        )
+
+
 if __name__ == "__main__":
     unittest.main()

test/gw/utils_test.py

@@ -203,9 +203,12 @@ def test_get_approximant(self):
             gwutils.lalsim_GetApproximantFromString(10)
 
     def test_lalsim_SimInspiralChooseFDWaveform(self):
+        # Masses are passed to LALSimulation in SI units, so use solar masses
+        # here. IMRPhenomXPHM (unlike IMRPhenomPv2) cannot handle unphysical
+        # sub-solar-mass values and aborts the process for them.
         a = gwutils.lalsim_SimInspiralChooseFDWaveform(
-            35.2,
-            20.4,
+            35.2 * lal.MSUN_SI,
+            20.4 * lal.MSUN_SI,
             0.1,
             0.2,
             0.2,
@@ -223,7 +226,7 @@ def test_lalsim_SimInspiralChooseFDWaveform(self):
             1000,
             20,
             None,
-            lalsim.IMRPhenomPv2,
+            lalsim.IMRPhenomXPHM,
         )
         self.assertEqual(len(a), 2)
         self.assertEqual(type(a[0]), lal.COMPLEX16FrequencySeries)

test/integration/make_standard_data.py

@@ -22,7 +22,7 @@
     luminosity_distance=100.0,
     theta_jn=0.4,
     phase=1.3,
-    waveform_approximant="IMRPhenomPv2",
+    waveform_approximant="IMRPhenomXPHM",
     reference_frequency=50.0,
     ra=1.375,
     dec=-1.2108,

test/integration/sample_from_the_prior_test.py

@@ -42,7 +42,7 @@ def test_fifteen_dimensional_cbc(self):
         bilby.core.utils.random.seed(8817021)
 
         waveform_arguments = dict(
-            waveform_approximant="IMRPhenomPv2",
+            waveform_approximant="IMRPhenomXPHM",
             reference_frequency=50.0,
             minimum_frequency=20.0,
         )