diff --git a/bilby/gw/source.py b/bilby/gw/source.py index 97ad51b51..f09e0b16c 100644 --- a/bilby/gw/source.py +++ b/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), diff --git a/pixi.toml b/pixi.toml index faad40e3b..4234b3d85 100644 --- a/pixi.toml +++ b/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 = "*" diff --git a/pyproject.toml b/pyproject.toml index 167dfedd8..8ae2e4139 100644 --- a/pyproject.toml +++ b/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/*"] diff --git a/test/gw/likelihood/roq_test.py b/test/gw/likelihood/roq_test.py index 34fa5504f..412181b9e 100644 --- a/test/gw/likelihood/roq_test.py +++ b/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) diff --git a/test/gw/plot_test.py b/test/gw/plot_test.py index 32313c9d4..bfcb1b9f9 100644 --- a/test/gw/plot_test.py +++ b/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), diff --git a/test/gw/result_test.py b/test/gw/result_test.py index 56e190449..01d46b34b 100644 --- a/test/gw/result_test.py +++ b/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), diff --git a/test/gw/source_test.py b/test/gw/source_test.py index 8295e675a..427935b64 100644 --- a/test/gw/source_test.py +++ b/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() diff --git a/test/gw/utils_test.py b/test/gw/utils_test.py index 9b6815264..c4601eec8 100644 --- a/test/gw/utils_test.py +++ b/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) diff --git a/test/integration/make_standard_data.py b/test/integration/make_standard_data.py index 4f5d4073d..c24f3e985 100644 --- a/test/integration/make_standard_data.py +++ b/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, diff --git a/test/integration/sample_from_the_prior_test.py b/test/integration/sample_from_the_prior_test.py index 7fce8b162..0882e352f 100644 --- a/test/integration/sample_from_the_prior_test.py +++ b/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, ) diff --git a/test/integration/test_waveforms.py b/test/integration/test_waveforms.py index c564c1052..9b026ff65 100644 --- a/test/integration/test_waveforms.py +++ b/test/integration/test_waveforms.py @@ -46,20 +46,20 @@ def setUp(self): lambda_2=1500, ) - def test_IMRPhenomPv2(self): - waveform_approximant = "IMRPhenomPv2" + def test_IMRPhenomXPHM(self): + waveform_approximant = "IMRPhenomXPHM" self.run_for_approximant(waveform_approximant, source="bbh") - def test_IMRPhenomD(self): - waveform_approximant = "IMRPhenomD" + def test_IMRPhenomXAS(self): + waveform_approximant = "IMRPhenomXAS" self.run_for_approximant(waveform_approximant, source="bbh") - def test_IMRPhenomPv2_NRTidal(self): - waveform_approximant = "IMRPhenomPv2_NRTidal" + def test_IMRPhenomPv2_NRTidalv2(self): + waveform_approximant = "IMRPhenomPv2_NRTidalv2" self.run_for_approximant(waveform_approximant, source="bns") - def test_IMRPhenomD_NRTidal(self): - waveform_approximant = "IMRPhenomD_NRTidal" + def test_IMRPhenomD_NRTidalv2(self): + waveform_approximant = "IMRPhenomD_NRTidalv2" self.run_for_approximant(waveform_approximant, source="bns") def test_TaylorF2(self):