changes to source.py

AUTHORS.md

@@ -40,6 +40,7 @@ Karl Wette
 Katerina Chatziioannou
 Kaylee de Soto
 Khun Sang Phukon
+Kruthi Krishna
 Kshipraa Athar
 Leslie Wade
 Liting Xiao

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.

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()