From 8c7b8d34d8bdd1835464960838571cd18050f9d3 Mon Sep 17 00:00:00 2001
From: Aditya Vijaykumar <vijaykumar.aditya@gmail.com>
Date: Mon, 16 May 2022 12:46:17 +0530
Subject: [PATCH] changes to source.py
---
AUTHORS.md | 1 +
bilby/gw/source.py | 189 +++++++++++++++++++++++++++++++++++++++
test/gw/source_test.py | 194 +++++++++++++++++++++++++++++++++++++++++
3 files changed, 384 insertions(+)
diff --git a/AUTHORS.md b/AUTHORS.md
index 98fae5bce..cdb687ce1 100644
--- a/AUTHORS.md
+++ b/AUTHORS.md
@@ -40,6 +40,7 @@ Karl Wette
Katerina Chatziioannou
Kaylee de Soto
Khun Sang Phukon
+Kruthi Krishna
Kshipraa Athar
Leslie Wade
Liting Xiao
diff --git a/bilby/gw/source.py b/bilby/gw/source.py
index 4d7793b92..f32af4bf2 100644
--- a/bilby/gw/source.py
+++ b/bilby/gw/source.py
@@ -828,6 +828,195 @@ def _base_waveform_frequency_sequence(
return dict(plus=h_plus.data.data, cross=h_cross.data.data)
+def lal_binary_black_hole_relativebinning(
+ frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
+ phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, fiducial, **kwargs):
+
+ waveform_kwargs = dict(
+ waveform_approximant='IMRPhenomPv2', reference_frequency=50.0,
+ minimum_frequency=20.0, maximum_frequency=frequency_array[-1],
+ catch_waveform_errors=False, pn_spin_order=-1, pn_tidal_order=-1,
+ pn_phase_order=-1, pn_amplitude_order=0)
+ waveform_kwargs.update(kwargs)
+
+ if fiducial == 1:
+ return _base_lal_cbc_fd_waveform(
+ frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
+ luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
+ a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, phi_jl=phi_jl,
+ phi_12=phi_12, lambda_1=0.0, lambda_2=0.0, **waveform_kwargs)
+
+ else:
+ return _base_relativebinning_waveform(
+ frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
+ luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
+ a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, phi_jl=phi_jl,
+ phi_12=phi_12, lambda_1=0.0, lambda_2=0.0, **waveform_kwargs)
+
+
+def lal_binary_neutron_star_relativebinning(
+ frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
+ phi_12, a_2, tilt_2, phi_jl, lambda_1, lambda_2, theta_jn, phase,
+ fiducial, **kwargs):
+
+ waveform_kwargs = dict(
+ waveform_approximant='IMRPhenomPv2_NRTidal', reference_frequency=50.0,
+ minimum_frequency=20.0, maximum_frequency=frequency_array[-1],
+ catch_waveform_errors=False, pn_spin_order=-1, pn_tidal_order=-1,
+ pn_phase_order=-1, pn_amplitude_order=0)
+ waveform_kwargs.update(kwargs)
+
+ if fiducial == 1:
+
+ return _base_lal_cbc_fd_waveform(
+ frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
+ luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
+ a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, phi_12=phi_12,
+ phi_jl=phi_jl, lambda_1=lambda_1, lambda_2=lambda_2, **waveform_kwargs)
+
+ else:
+ return _base_relativebinning_waveform(
+ frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
+ luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
+ a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, phi_jl=phi_jl,
+ phi_12=phi_12, lambda_1=lambda_1, lambda_2=lambda_2, **waveform_kwargs)
+
+
+def _base_relativebinning_waveform(
+ frequency_array, mass_1, mass_2, luminosity_distance, theta_jn, phase,
+ a_1=0.0, a_2=0.0, tilt_1=0.0, tilt_2=0.0, phi_12=0.0, phi_jl=0.0,
+ lambda_1=0.0, lambda_2=0.0, **waveform_kwargs):
+ """ Generate a cbc waveform model using lalsimulation
+
+ Parameters
+ ----------
+ frequency_array: array_like
+ The frequencies at which we want to calculate the strain. Ignored
+ for the relative binning model.
+ mass_1: float
+ The mass of the heavier object in solar masses
+ mass_2: float
+ The mass of the lighter object in solar masses
+ luminosity_distance: float
+ The luminosity distance in megaparsec
+ a_1: float
+ Dimensionless primary spin magnitude
+ tilt_1: float
+ Primary tilt angle
+ phi_12: float
+ Azimuthal angle between the component spins
+ a_2: float
+ Dimensionless secondary spin magnitude
+ tilt_2: float
+ Secondary tilt angle
+ phi_jl: float
+ Azimuthal angle between the total and orbital angular momenta
+ theta_jn: float
+ Orbital inclination
+ phase: float
+ The phase at coalescence
+ lambda_1: float
+ Tidal deformability of the more massive object
+ lambda_2: float
+ Tidal deformability of the less massive object
+ kwargs: dict
+ Optional keyword arguments
+
+ Waveform arguments
+ ------------------
+ Non-sampled extra data used in the source model calculation
+ frequency_bin_edges: np.array
+
+ Returns
+ -------
+ dict: A dictionary with the plus and cross polarisation strain modes
+ """
+ import lal
+ import lalsimulation as lalsim
+
+ waveform_approximant = waveform_kwargs['waveform_approximant']
+ reference_frequency = waveform_kwargs['reference_frequency']
+ catch_waveform_errors = waveform_kwargs['catch_waveform_errors']
+ pn_spin_order = waveform_kwargs['pn_spin_order']
+ pn_tidal_order = waveform_kwargs['pn_tidal_order']
+ pn_phase_order = waveform_kwargs['pn_phase_order']
+ pn_amplitude_order = waveform_kwargs['pn_amplitude_order']
+ waveform_dictionary = waveform_kwargs.get(
+ 'lal_waveform_dictionary', lal.CreateDict()
+ )
+
+ frequency_bin_edges = waveform_kwargs['frequency_bin_edges']
+
+ approximant = lalsim_GetApproximantFromString(waveform_approximant)
+
+ luminosity_distance = luminosity_distance * 1e6 * utils.parsec
+ mass_1 = mass_1 * utils.solar_mass
+ mass_2 = mass_2 * utils.solar_mass
+
+ iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = bilby_to_lalsimulation_spins(
+ theta_jn=theta_jn, phi_jl=phi_jl, tilt_1=tilt_1, tilt_2=tilt_2,
+ phi_12=phi_12, a_1=a_1, a_2=a_2, mass_1=mass_1, mass_2=mass_2,
+ reference_frequency=reference_frequency, phase=phase)
+
+ lalsim.SimInspiralWaveformParamsInsertPNSpinOrder(
+ waveform_dictionary, int(pn_spin_order))
+ lalsim.SimInspiralWaveformParamsInsertPNTidalOrder(
+ waveform_dictionary, int(pn_tidal_order))
+ lalsim.SimInspiralWaveformParamsInsertPNPhaseOrder(
+ waveform_dictionary, int(pn_phase_order))
+ lalsim.SimInspiralWaveformParamsInsertPNAmplitudeOrder(
+ waveform_dictionary, int(pn_amplitude_order))
+ lalsim_SimInspiralWaveformParamsInsertTidalLambda1(
+ waveform_dictionary, lambda_1)
+ lalsim_SimInspiralWaveformParamsInsertTidalLambda2(
+ waveform_dictionary, lambda_2)
+
+ for key, value in waveform_kwargs.items():
+ func = getattr(lalsim, "SimInspiralWaveformParamsInsert" + key, None)
+ if func is not None:
+ func(waveform_dictionary, value)
+
+ if waveform_kwargs.get('numerical_relativity_file', None) is not None:
+ lalsim.SimInspiralWaveformParamsInsertNumRelData(
+ waveform_dictionary, waveform_kwargs['numerical_relativity_file'])
+
+ if ('mode_array' in waveform_kwargs) and waveform_kwargs['mode_array'] is not None:
+ mode_array = waveform_kwargs['mode_array']
+ mode_array_lal = lalsim.SimInspiralCreateModeArray()
+ for mode in mode_array:
+ lalsim.SimInspiralModeArrayActivateMode(mode_array_lal, mode[0], mode[1])
+ lalsim.SimInspiralWaveformParamsInsertModeArray(waveform_dictionary, mode_array_lal)
+
+ try:
+ h_binned_plus, h_binned_cross = lalsim_SimInspiralChooseFDWaveformSequence(
+ phase, mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y,
+ spin_2z, reference_frequency, luminosity_distance, iota,
+ waveform_dictionary, approximant, frequency_bin_edges)
+ except Exception as e:
+ if not catch_waveform_errors:
+ raise
+ else:
+ EDOM = (e.args[0] == 'Internal function call failed: Input domain error')
+ if EDOM:
+ failed_parameters = dict(mass_1=mass_1, mass_2=mass_2,
+ spin_1=(spin_1x, spin_2y, spin_1z),
+ spin_2=(spin_2x, spin_2y, spin_2z),
+ luminosity_distance=luminosity_distance,
+ iota=iota, phase=phase,
+ )
+ logger.warning("Evaluating the waveform failed with error: {}\n".format(e) +
+ "The parameters were {}\n".format(failed_parameters) +
+ "Likelihood will be set to -inf.")
+ return None
+ else:
+ raise
+
+ waveform_polarizations = dict(
+ plus=h_binned_plus.data.data, cross=h_binned_cross.data.data)
+
+ return waveform_polarizations
+
+
def sinegaussian(frequency_array, hrss, Q, frequency, **kwargs):
r"""
A frequency-domain sine-Gaussian burst source model.
diff --git a/test/gw/source_test.py b/test/gw/source_test.py
index 627278a33..f3fcc1201 100644
--- a/test/gw/source_test.py
+++ b/test/gw/source_test.py
@@ -423,5 +423,199 @@ class TestBNSfreqseq(unittest.TestCase):
self.assertLess(diff / norm, 1e-5)
+class TestRelbinBBH(unittest.TestCase):
+ def setUp(self):
+ self.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.waveform_kwargs_fiducial = dict(
+ waveform_approximant="IMRPhenomPv2",
+ reference_frequency=50.0,
+ minimum_frequency=20.0,
+ catch_waveform_errors=True,
+ fiducial=True,
+ )
+ self.waveform_kwargs_binned = dict(
+ waveform_approximant="IMRPhenomPv2",
+ reference_frequency=50.0,
+ minimum_frequency=20.0,
+ catch_waveform_errors=True,
+ fiducial=False,
+ frequency_bin_edges=np.arange(20, 1500, 50)
+ )
+ self.frequency_array = bilby.core.utils.create_frequency_series(2048, 4)
+ self.bad_parameters = copy(self.parameters)
+ self.bad_parameters["mass_1"] = -30.0
+
+ def tearDown(self):
+ del self.parameters
+ del self.waveform_kwargs_fiducial
+ del self.waveform_kwargs_binned
+ del self.frequency_array
+ del self.bad_parameters
+
+ def test_relbin_fiducial_bbh_works_runs_valid_parameters(self):
+ self.parameters.update(self.waveform_kwargs_fiducial)
+ self.assertIsInstance(
+ bilby.gw.source.lal_binary_black_hole_relativebinning(
+ self.frequency_array, **self.parameters
+ ),
+ dict,
+ )
+
+ def test_relbin_binned_bbh_works_runs_valid_parameters(self):
+ self.parameters.update(self.waveform_kwargs_binned)
+ self.assertIsInstance(
+ bilby.gw.source.lal_binary_black_hole_relativebinning(
+ self.frequency_array, **self.parameters
+ ),
+ dict,
+ )
+
+ def test_waveform_error_catching_fiducial(self):
+ self.bad_parameters.update(self.waveform_kwargs_fiducial)
+ self.assertIsNone(
+ bilby.gw.source.lal_binary_black_hole_relativebinning(
+ self.frequency_array, **self.bad_parameters
+ )
+ )
+
+ def test_waveform_error_catching_binned(self):
+ self.bad_parameters.update(self.waveform_kwargs_binned)
+ self.assertIsNone(
+ bilby.gw.source.lal_binary_black_hole_relativebinning(
+ self.frequency_array, **self.bad_parameters
+ )
+ )
+
+ def test_waveform_error_raising_fiducial(self):
+ raise_error_parameters = copy(self.bad_parameters)
+ raise_error_parameters.update(self.waveform_kwargs_fiducial)
+ raise_error_parameters["catch_waveform_errors"] = False
+ with self.assertRaises(Exception):
+ bilby.gw.source.lal_binary_black_hole_relativebinning(
+ self.frequency_array, **raise_error_parameters
+ )
+
+ def test_waveform_error_raising_binned(self):
+ raise_error_parameters = copy(self.bad_parameters)
+ raise_error_parameters.update(self.waveform_kwargs_binned)
+ raise_error_parameters["catch_waveform_errors"] = False
+ with self.assertRaises(Exception):
+ bilby.gw.source.lal_binary_black_hole_relativebinning(
+ self.frequency_array, **raise_error_parameters
+ )
+
+ def test_relbin_bbh_fails_without_fiducial_option(self):
+ with self.assertRaises(KeyError):
+ bilby.gw.source.lal_binary_black_hole_relativebinning(
+ self.frequency_array, **self.parameters
+ )
+
+ def test_relbin_bbh_xpprecession_version(self):
+ self.parameters.update(self.waveform_kwargs_fiducial)
+ self.parameters["waveform_approximant"] = "IMRPhenomXP"
+
+ # Test that we can modify the XP precession version
+ out_v223 = bilby.gw.source.lal_binary_black_hole(
+ self.frequency_array, PhenomXPrecVersion=223, **self.parameters
+ )
+ out_v102 = bilby.gw.source.lal_binary_black_hole(
+ self.frequency_array, PhenomXPrecVersion=102, **self.parameters
+ )
+ self.assertFalse(np.all(out_v223["plus"] == out_v102["plus"]))
+
+
+class TestRelbinBNS(unittest.TestCase):
+ def setUp(self):
+ self.parameters = dict(
+ mass_1=1.4,
+ mass_2=1.4,
+ luminosity_distance=400.0,
+ a_1=0.4,
+ a_2=0.3,
+ tilt_1=0.2,
+ tilt_2=1.7,
+ phi_jl=0.2,
+ phi_12=0.9,
+ theta_jn=1.7,
+ phase=0.0,
+ lambda_1=100.0,
+ lambda_2=100.0,
+ )
+ self.waveform_kwargs_fiducial = dict(
+ waveform_approximant="IMRPhenomPv2_NRTidal",
+ reference_frequency=50.0,
+ minimum_frequency=20.0,
+ fiducial=True,
+ )
+ self.waveform_kwargs_binned = dict(
+ waveform_approximant="IMRPhenomPv2_NRTidal",
+ reference_frequency=50.0,
+ minimum_frequency=20.0,
+ fiducial=False,
+ frequency_bin_edges=np.arange(20, 1500, 50),
+ )
+ self.frequency_array = bilby.core.utils.create_frequency_series(2048, 4)
+
+ def tearDown(self):
+ del self.parameters
+ del self.waveform_kwargs_fiducial
+ del self.waveform_kwargs_binned
+ del self.frequency_array
+
+ def test_relbin_fiducial_bns_runs_with_valid_parameters(self):
+ self.parameters.update(self.waveform_kwargs_fiducial)
+ self.assertIsInstance(
+ bilby.gw.source.lal_binary_neutron_star_relativebinning(
+ self.frequency_array, **self.parameters
+ ),
+ dict,
+ )
+
+ def test_relbin_binned_bns_runs_with_valid_parameters(self):
+ self.parameters.update(self.waveform_kwargs_binned)
+ self.assertIsInstance(
+ bilby.gw.source.lal_binary_neutron_star_relativebinning(
+ self.frequency_array, **self.parameters
+ ),
+ dict,
+ )
+
+ def test_relbin_bns_fails_without_fiducial_option(self):
+ with self.assertRaises(KeyError):
+ bilby.gw.source.lal_binary_neutron_star_relativebinning(
+ self.frequency_array, **self.parameters
+ )
+
+ def test_fiducial_fails_without_tidal_parameters(self):
+ self.parameters.pop("lambda_1")
+ self.parameters.pop("lambda_2")
+ self.parameters.update(self.waveform_kwargs_fiducial)
+ with self.assertRaises(TypeError):
+ bilby.gw.source.lal_binary_neutron_star_relativebinning(
+ self.frequency_array, **self.parameters
+ )
+
+ def test_binned_fails_without_tidal_parameters(self):
+ self.parameters.pop("lambda_1")
+ self.parameters.pop("lambda_2")
+ self.parameters.update(self.waveform_kwargs_binned)
+ with self.assertRaises(TypeError):
+ bilby.gw.source.lal_binary_neutron_star_relativebinning(
+ self.frequency_array, **self.parameters
+ )
+
+
if __name__ == "__main__":
unittest.main()
--
GitLab