diff --git a/AUTHORS.md b/AUTHORS.md
index afa1592dcb8ad3a95769e962cbf0c0ab63ac9c47..7bf380bd69fde382046bfcc339fef7294cb75e12 100644
--- a/AUTHORS.md
+++ b/AUTHORS.md
@@ -9,6 +9,7 @@ Aditya Vijaykumar
Andrew Kim
Andrew Miller
Antoni Ramos-Buades
+Apratim Ganguly
Avi Vajpeyi
Bruce Edelman
Carl-Johan Haster
@@ -40,6 +41,7 @@ Karl Wette
Katerina Chatziioannou
Kaylee de Soto
Khun Sang Phukon
+Kruthi Krishna
Kshipraa Athar
Leslie Wade
Liting Xiao
diff --git a/bilby/gw/conversion.py b/bilby/gw/conversion.py
index c49a7217bdca5a80cf97325012d9cb291b5545f4..3d207e4833ee7facb3df9ee25791bbc0ebb5e1de 100644
--- a/bilby/gw/conversion.py
+++ b/bilby/gw/conversion.py
@@ -1407,11 +1407,9 @@ def compute_snrs(sample, likelihood, npool=1):
if likelihood is not None:
if isinstance(sample, dict):
likelihood.parameters.update(sample)
- signal_polarizations =\
- likelihood.waveform_generator.frequency_domain_strain(sample)
+ signal_polarizations = likelihood.waveform_generator.frequency_domain_strain(sample)
for ifo in likelihood.interferometers:
- per_detector_snr = likelihood.calculate_snrs(
- signal_polarizations, ifo)
+ per_detector_snr = likelihood.calculate_snrs(signal_polarizations, ifo)
sample['{}_matched_filter_snr'.format(ifo.name)] =\
per_detector_snr.complex_matched_filter_snr
sample['{}_optimal_snr'.format(ifo.name)] = \
diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py
index 3eea4bcf0dc0a00094f27a506f21d588991c383d..9f8c29147aec9fab09f7a7dc63882588c44fd2e0 100644
--- a/bilby/gw/detector/interferometer.py
+++ b/bilby/gw/detector/interferometer.py
@@ -285,7 +285,7 @@ class Interferometer(object):
else:
return 0
- def get_detector_response(self, waveform_polarizations, parameters):
+ def get_detector_response(self, waveform_polarizations, parameters, frequencies=None):
""" Get the detector response for a particular waveform
Parameters
@@ -294,11 +294,21 @@ class Interferometer(object):
polarizations of the waveform
parameters: dict
parameters describing position and time of arrival of the signal
-
+ frequencies: array-like, optional
+ The frequency values to evaluate the response at. If
+ not provided, the response is computed using
+ :code:`self.frequency_array`. If the frequencies are
+ specified, no frequency masking is performed.
Returns
=======
array_like: A 3x3 array representation of the detector response (signal observed in the interferometer)
"""
+ if frequencies is None:
+ frequencies = self.frequency_array[self.frequency_mask]
+ mask = self.frequency_mask
+ else:
+ mask = np.ones(len(frequencies), dtype=bool)
+
signal = {}
for mode in waveform_polarizations.keys():
det_response = self.antenna_response(
@@ -308,9 +318,7 @@ class Interferometer(object):
parameters['psi'], mode)
signal[mode] = waveform_polarizations[mode] * det_response
- signal_ifo = sum(signal.values())
-
- signal_ifo *= self.strain_data.frequency_mask
+ signal_ifo = sum(signal.values()) * mask
time_shift = self.time_delay_from_geocenter(
parameters['ra'], parameters['dec'], parameters['geocent_time'])
@@ -320,12 +328,11 @@ class Interferometer(object):
dt_geocent = parameters['geocent_time'] - self.strain_data.start_time
dt = dt_geocent + time_shift
- signal_ifo[self.strain_data.frequency_mask] = signal_ifo[self.strain_data.frequency_mask] * np.exp(
- -1j * 2 * np.pi * dt * self.strain_data.frequency_array[self.strain_data.frequency_mask])
+ signal_ifo[mask] = signal_ifo[mask] * np.exp(-1j * 2 * np.pi * dt * frequencies)
- signal_ifo[self.strain_data.frequency_mask] *= self.calibration_model.get_calibration_factor(
- self.strain_data.frequency_array[self.strain_data.frequency_mask],
- prefix='recalib_{}_'.format(self.name), **parameters)
+ signal_ifo[mask] *= self.calibration_model.get_calibration_factor(
+ frequencies, prefix='recalib_{}_'.format(self.name), **parameters
+ )
return signal_ifo
diff --git a/bilby/gw/likelihood/__init__.py b/bilby/gw/likelihood/__init__.py
index 88fb527ad579513e39c4cda69e0d20a9c1d38e84..d752f77a915a54000d6a65e6e66b3b027828665f 100644
--- a/bilby/gw/likelihood/__init__.py
+++ b/bilby/gw/likelihood/__init__.py
@@ -2,6 +2,7 @@ from .base import GravitationalWaveTransient
from .basic import BasicGravitationalWaveTransient
from .roq import BilbyROQParamsRangeError, ROQGravitationalWaveTransient
from .multiband import MBGravitationalWaveTransient
+from .relative import RelativeBinningGravitationalWaveTransient
from ..source import lal_binary_black_hole
from ..waveform_generator import WaveformGenerator
diff --git a/bilby/gw/likelihood/base.py b/bilby/gw/likelihood/base.py
index 985fd94b6e925c322b00e8e33e0f4365739555d5..53621643b6d68b1fef504f0ff89b56b92ac74a4c 100644
--- a/bilby/gw/likelihood/base.py
+++ b/bilby/gw/likelihood/base.py
@@ -265,8 +265,10 @@ class GravitationalWaveTransient(Likelihood):
The bilby interferometer object
"""
- signal = interferometer.get_detector_response(
- waveform_polarizations, self.parameters)
+ signal = self._compute_full_waveform(
+ signal_polarizations=waveform_polarizations,
+ interferometer=interferometer,
+ )
_mask = interferometer.frequency_mask
if 'recalib_index' in self.parameters:
@@ -614,8 +616,10 @@ class GravitationalWaveTransient(Likelihood):
for ifo in self.interferometers:
ifo_length = len(ifo.frequency_domain_strain)
mask = ifo.frequency_mask
- signal = ifo.get_detector_response(
- signal_polarizations, self.parameters)
+ signal = self._compute_full_waveform(
+ signal_polarizations=signal_polarizations,
+ interferometer=ifo,
+ )
signal_long[:ifo_length] = signal
data[:ifo_length] = np.conj(ifo.frequency_domain_strain)
psd[:ifo_length][mask] = ifo.power_spectral_density_array[mask]
@@ -703,6 +707,23 @@ class GravitationalWaveTransient(Likelihood):
h_inner_h += per_detector_snr.optimal_snr_squared
return d_inner_h, h_inner_h
+ def _compute_full_waveform(self, signal_polarizations, interferometer):
+ """
+ Project the waveform polarizations against the interferometer
+ response. This is useful for likelihood classes that don't
+ use the full frequency array, e.g., the relative binning
+ likelihood.
+
+ Parameters
+ ==========
+ signal_polarizations: dict
+ Dictionary containing the waveform evaluated at
+ :code:`interferometer.frequency_array`.
+ interferometer: bilby.gw.detector.Interferometer
+ Interferometer to compute the response with respect to.
+ """
+ return interferometer.get_detector_response(signal_polarizations, self.parameters)
+
def generate_phase_sample_from_marginalized_likelihood(
self, signal_polarizations=None):
r"""
diff --git a/bilby/gw/likelihood/relative.py b/bilby/gw/likelihood/relative.py
new file mode 100644
index 0000000000000000000000000000000000000000..ddb4d3ee872b3f2c43efb181f4b3fdaba4b3db08
--- /dev/null
+++ b/bilby/gw/likelihood/relative.py
@@ -0,0 +1,395 @@
+import numpy as np
+from scipy.optimize import differential_evolution
+
+from .base import GravitationalWaveTransient
+from ...core.utils import logger
+from ...core.prior.base import Constraint
+from ...core.prior import DeltaFunction
+from ..utils import noise_weighted_inner_product
+
+
+class RelativeBinningGravitationalWaveTransient(GravitationalWaveTransient):
+ """A gravitational-wave transient likelihood object which uses the relative
+ binning procedure to calculate a fast likelihood. See Zackay et al.
+ arXiv1806.08792
+
+ Parameters
+ ----------
+ interferometers: list, bilby.gw.detector.InterferometerList
+ A list of `bilby.detector.Interferometer` instances - contains the
+ detector data and power spectral densities
+ waveform_generator: `bilby.waveform_generator.WaveformGenerator`
+ An object which computes the frequency-domain strain of the signal,
+ given some set of parameters
+ fiducial_parameters: dict, optional
+ A starting guess for initial parameters of the event for finding the
+ maximum likelihood (fiducial) waveform.
+ parameter_bounds: dict, optional
+ Dictionary of bounds (lists) for the initial parameters when finding
+ the initial maximum likelihood (fiducial) waveform.
+ distance_marginalization: bool, optional
+ If true, marginalize over distance in the likelihood.
+ This uses a look up table calculated at run time.
+ The distance prior is set to be a delta function at the minimum
+ distance allowed in the prior being marginalised over.
+ time_marginalization: bool, optional
+ If true, marginalize over time in the likelihood.
+ This uses a FFT to calculate the likelihood over a regularly spaced
+ grid.
+ In order to cover the whole space the prior is set to be uniform over
+ the spacing of the array of times.
+ If using time marginalisation and jitter_time is True a "jitter"
+ parameter is added to the prior which modifies the position of the
+ grid of times.
+ phase_marginalization: bool, optional
+ If true, marginalize over phase in the likelihood.
+ This is done analytically using a Bessel function.
+ The phase prior is set to be a delta function at phase=0.
+ priors: dict, optional
+ If given, used in the distance and phase marginalization.
+ distance_marginalization_lookup_table: (dict, str), optional
+ If a dict, dictionary containing the lookup_table, distance_array,
+ (distance) prior_array, and reference_distance used to construct
+ the table.
+ If a string the name of a file containing these quantities.
+ The lookup table is stored after construction in either the
+ provided string or a default location:
+ '.distance_marginalization_lookup_dmin{}_dmax{}_n{}.npz'
+ jitter_time: bool, optional
+ Whether to introduce a `time_jitter` parameter. This avoids either
+ missing the likelihood peak, or introducing biases in the
+ reconstructed time posterior due to an insufficient sampling frequency.
+ Default is False, however using this parameter is strongly encouraged.
+ reference_frame: (str, bilby.gw.detector.InterferometerList, list), optional
+ Definition of the reference frame for the sky location.
+ - "sky": sample in RA/dec, this is the default
+ - e.g., "H1L1", ["H1", "L1"], InterferometerList(["H1", "L1"]):
+ sample in azimuth and zenith, `azimuth` and `zenith` defined in the
+ frame where the z-axis is aligned the the vector connecting H1
+ and L1.
+ time_reference: str, optional
+ Name of the reference for the sampled time parameter.
+ - "geocent"/"geocenter": sample in the time at the Earth's center,
+ this is the default
+ - e.g., "H1": sample in the time of arrival at H1
+ chi: float, optional
+ Tunable parameter which limits the perturbation of alpha when setting
+ up the bin range. See https://arxiv.org/abs/1806.08792.
+ epsilon: float, optional
+ Tunable parameter which limits the differential phase change in each
+ bin when setting up the bin range. See https://arxiv.org/abs/1806.08792.
+
+ Returns
+ -------
+ Likelihood: `bilby.core.likelihood.Likelihood`
+ A likelihood object, able to compute the likelihood of the data given
+ some model parameters.
+
+ Notes
+ -----
+ The relative binning likelihood does not currently support calibration marginalization.
+ """
+
+ def __init__(self, interferometers,
+ waveform_generator,
+ fiducial_parameters=None,
+ parameter_bounds=None,
+ maximization_kwargs=None,
+ update_fiducial_parameters=False,
+ distance_marginalization=False,
+ time_marginalization=False,
+ phase_marginalization=False,
+ priors=None,
+ distance_marginalization_lookup_table=None,
+ jitter_time=True,
+ reference_frame="sky",
+ time_reference="geocenter",
+ chi=1,
+ epsilon=0.5):
+
+ super(RelativeBinningGravitationalWaveTransient, self).__init__(
+ interferometers=interferometers,
+ waveform_generator=waveform_generator,
+ distance_marginalization=distance_marginalization,
+ phase_marginalization=phase_marginalization,
+ time_marginalization=time_marginalization,
+ priors=priors,
+ distance_marginalization_lookup_table=distance_marginalization_lookup_table,
+ jitter_time=jitter_time,
+ reference_frame=reference_frame,
+ time_reference=time_reference)
+
+ if fiducial_parameters is None:
+ logger.info("Drawing fiducial parameters from prior.")
+ fiducial_parameters = priors.sample()
+ fiducial_parameters["fiducial"] = 0
+ if self.time_marginalization:
+ fiducial_parameters["geocent_time"] = interferometers.start_time
+ if self.distance_marginalization:
+ fiducial_parameters["luminosity_distance"] = self._ref_dist
+ if self.phase_marginalization:
+ fiducial_parameters["phase"] = 0.0
+ self.fiducial_parameters = fiducial_parameters
+ self.chi = chi
+ self.epsilon = epsilon
+ self.gamma = np.array([-5 / 3, -2 / 3, 1, 5 / 3, 7 / 3])
+ self.maximum_frequency = waveform_generator.frequency_array[-1]
+ self.fiducial_waveform_obtained = False
+ self.check_if_bins_are_setup = False
+ self.fiducial_polarizations = None
+ self.per_detector_fiducial_waveforms = dict()
+ self.per_detector_fiducial_waveform_points = dict()
+ self.bin_freqs = dict()
+ self.bin_inds = dict()
+ self.bin_widths = dict()
+ self.bin_centers = dict()
+ self.set_fiducial_waveforms(self.fiducial_parameters)
+ logger.info("Initial fiducial waveforms set up")
+ self.setup_bins()
+ self.compute_summary_data()
+ logger.info("Summary Data Obtained")
+
+ if update_fiducial_parameters:
+ # write a check to make sure prior is not None
+ logger.info("Using scipy optimization to find maximum likelihood parameters.")
+ self.parameters_to_be_updated = [key for key in priors if not isinstance(
+ priors[key], (DeltaFunction, Constraint, float, int))]
+ logger.info(f"Parameters over which likelihood is maximized: {self.parameters_to_be_updated}")
+ if parameter_bounds is None:
+ logger.info("No parameter bounds were given. Using priors instead.")
+ self.parameter_bounds = self.get_bounds_from_priors(priors)
+ else:
+ self.parameter_bounds = self.get_parameter_list_from_dictionary(parameter_bounds)
+ self.fiducial_parameters = self.find_maximum_likelihood_parameters(
+ self.parameter_bounds, maximization_kwargs=maximization_kwargs)
+ self.parameters.update(self.fiducial_parameters)
+ logger.info(f"Fiducial likelihood: {self.log_likelihood_ratio():.2f}")
+
+ def __repr__(self):
+ return self.__class__.__name__ + '(interferometers={},\n\twaveform_generator={},\n\fiducial_parameters={},' \
+ .format(self.interferometers, self.waveform_generator, self.fiducial_parameters)
+
+ def setup_bins(self):
+ frequency_array = self.waveform_generator.frequency_array
+ gamma = self.gamma[:, np.newaxis]
+ maximum_frequency = frequency_array[0]
+ minimum_frequency = frequency_array[-1]
+ for interferometer in self.interferometers:
+ maximum_frequency = max(maximum_frequency, interferometer.maximum_frequency)
+ minimum_frequency = min(minimum_frequency, interferometer.minimum_frequency)
+ maximum_frequency = min(maximum_frequency, self.maximum_frequency)
+ frequency_array_useful = frequency_array[
+ (frequency_array >= minimum_frequency)
+ & (frequency_array <= maximum_frequency)
+ ]
+
+ d_alpha = self.chi * 2 * np.pi / np.abs(
+ (minimum_frequency ** gamma) * np.heaviside(-gamma, 1)
+ - (maximum_frequency ** gamma) * np.heaviside(gamma, 1)
+ )
+ d_phi = np.sum(
+ np.sign(gamma) * d_alpha * frequency_array_useful ** gamma,
+ axis=0
+ )
+ d_phi_from_start = d_phi - d_phi[0]
+ self.number_of_bins = int(d_phi_from_start[-1] // self.epsilon)
+ self.bin_freqs = np.zeros(self.number_of_bins + 1)
+ self.bin_inds = np.zeros(self.number_of_bins + 1, dtype=int)
+
+ for i in range(self.number_of_bins + 1):
+ bin_index = np.where(d_phi_from_start >= ((i / self.number_of_bins) * d_phi_from_start[-1]))[0][0]
+ bin_freq = frequency_array_useful[bin_index]
+ self.bin_freqs[i] = bin_freq
+ self.bin_inds[i] = np.where(frequency_array >= bin_freq)[0][0]
+
+ logger.debug(
+ f"Set up {self.number_of_bins} bins "
+ f"between {minimum_frequency} Hz and {maximum_frequency} Hz"
+ )
+ self.waveform_generator.waveform_arguments["frequency_bin_edges"] = self.bin_freqs
+ self.bin_widths = self.bin_freqs[1:] - self.bin_freqs[:-1]
+ self.bin_centers = (self.bin_freqs[1:] + self.bin_freqs[:-1]) / 2
+ for interferometer in self.interferometers:
+ name = interferometer.name
+ self.per_detector_fiducial_waveform_points[name] = (
+ self.per_detector_fiducial_waveforms[name][self.bin_inds]
+ )
+
+ def set_fiducial_waveforms(self, parameters):
+ _fiducial = parameters["fiducial"]
+ parameters["fiducial"] = 1
+ self.fiducial_polarizations = self.waveform_generator.frequency_domain_strain(
+ parameters)
+
+ maximum_nonzero_index = np.where(self.fiducial_polarizations["plus"] != 0j)[0][-1]
+ logger.debug(f"Maximum Nonzero Index is {maximum_nonzero_index}")
+ maximum_nonzero_frequency = self.waveform_generator.frequency_array[maximum_nonzero_index]
+ logger.debug(f"Maximum Nonzero Frequency is {maximum_nonzero_frequency}")
+ self.maximum_frequency = maximum_nonzero_frequency
+
+ if self.fiducial_polarizations is None:
+ raise ValueError(f"Cannot compute fiducial waveforms for {parameters}")
+
+ for interferometer in self.interferometers:
+ logger.debug(f"Maximum Frequency is {interferometer.maximum_frequency}")
+ wf = interferometer.get_detector_response(self.fiducial_polarizations, parameters)
+ wf[interferometer.frequency_array > self.maximum_frequency] = 0
+ self.per_detector_fiducial_waveforms[interferometer.name] = wf
+
+ parameters["fiducial"] = _fiducial
+
+ def find_maximum_likelihood_parameters(self, parameter_bounds,
+ iterations=5, maximization_kwargs=None):
+ if maximization_kwargs is None:
+ maximization_kwargs = dict()
+ self.parameters.update(self.fiducial_parameters)
+ self.parameters["fiducial"] = 0
+ updated_parameters_list = self.get_parameter_list_from_dictionary(self.fiducial_parameters)
+ old_fiducial_ln_likelihood = self.log_likelihood_ratio()
+ logger.info(f"Fiducial ln likelihood ratio: {old_fiducial_ln_likelihood:.2f}")
+ for it in range(iterations):
+ logger.info(f"Optimizing fiducial parameters. Iteration : {it + 1}")
+ output = differential_evolution(
+ self.lnlike_scipy_maximize,
+ bounds=parameter_bounds,
+ x0=updated_parameters_list,
+ **maximization_kwargs,
+ )
+ updated_parameters_list = output['x']
+ updated_parameters = self.get_parameter_dictionary_from_list(updated_parameters_list)
+ self.parameters.update(updated_parameters)
+ self.set_fiducial_waveforms(updated_parameters)
+ self.setup_bins()
+ self.compute_summary_data()
+ new_fiducial_ln_likelihood = self.log_likelihood_ratio()
+ logger.info(f"Fiducial ln likelihood ratio: {new_fiducial_ln_likelihood:.2f}")
+ if new_fiducial_ln_likelihood - old_fiducial_ln_likelihood < 0.1:
+ break
+ old_fiducial_ln_likelihood = new_fiducial_ln_likelihood
+
+ logger.info("Fiducial waveforms updated")
+ logger.info("Summary Data updated")
+ return updated_parameters
+
+ def lnlike_scipy_maximize(self, parameter_list):
+ self.parameters.update(self.get_parameter_dictionary_from_list(parameter_list))
+ return -self.log_likelihood_ratio()
+
+ def get_parameter_dictionary_from_list(self, parameter_list):
+ parameter_dictionary = dict(zip(self.parameters_to_be_updated, parameter_list))
+ excluded_parameter_keys = set(self.fiducial_parameters) - set(self.parameters_to_be_updated)
+ for key in excluded_parameter_keys:
+ parameter_dictionary[key] = self.fiducial_parameters[key]
+ return parameter_dictionary
+
+ def get_parameter_list_from_dictionary(self, parameter_dict):
+ return [parameter_dict[k] for k in self.parameters_to_be_updated]
+
+ def get_bounds_from_priors(self, priors):
+ bounds = []
+ for key in self.parameters_to_be_updated:
+ bounds.append([priors[key].minimum, priors[key].maximum])
+ return bounds
+
+ def compute_summary_data(self):
+ summary_data = dict()
+
+ for interferometer in self.interferometers:
+ mask = interferometer.frequency_mask
+ masked_frequency_array = interferometer.frequency_array[mask]
+ masked_bin_inds = []
+ for edge in self.bin_freqs:
+ index = np.where(masked_frequency_array == edge)[0][0]
+ masked_bin_inds.append(index)
+ masked_strain = interferometer.frequency_domain_strain[mask]
+ masked_h0 = self.per_detector_fiducial_waveforms[interferometer.name][mask]
+ masked_psd = interferometer.power_spectral_density_array[mask]
+ duration = interferometer.duration
+ a0, b0, a1, b1 = np.zeros((4, self.number_of_bins), dtype=complex)
+
+ for i in range(self.number_of_bins):
+ start_idx = masked_bin_inds[i]
+ end_idx = masked_bin_inds[i + 1]
+ start = masked_frequency_array[start_idx]
+ stop = masked_frequency_array[end_idx]
+ idxs = slice(start_idx, end_idx)
+
+ strain = masked_strain[idxs]
+ h0 = masked_h0[idxs]
+ psd = masked_psd[idxs]
+
+ frequencies = masked_frequency_array[idxs]
+ central_frequency = (start + stop) / 2
+ delta_frequency = frequencies - central_frequency
+
+ a0[i] = noise_weighted_inner_product(h0, strain, psd, duration)
+ b0[i] = noise_weighted_inner_product(h0, h0, psd, duration)
+ a1[i] = noise_weighted_inner_product(h0, strain * delta_frequency, psd, duration)
+ b1[i] = noise_weighted_inner_product(h0, h0 * delta_frequency, psd, duration)
+
+ summary_data[interferometer.name] = (a0, a1, b0, b1)
+
+ self.summary_data = summary_data
+
+ def compute_waveform_ratio_per_interferometer(self, waveform_polarizations, interferometer):
+ name = interferometer.name
+ strain = interferometer.get_detector_response(
+ waveform_polarizations=waveform_polarizations,
+ parameters=self.parameters,
+ frequencies=self.bin_freqs,
+ )
+ reference_strain = self.per_detector_fiducial_waveform_points[name]
+ waveform_ratio = strain / reference_strain
+
+ r0 = (waveform_ratio[1:] + waveform_ratio[:-1]) / 2
+ r1 = (waveform_ratio[1:] - waveform_ratio[:-1]) / self.bin_widths
+
+ return [r0, r1]
+
+ def _compute_full_waveform(self, signal_polarizations, interferometer):
+ fiducial_waveform = self.per_detector_fiducial_waveforms[interferometer.name]
+ r0, r1 = self.compute_waveform_ratio_per_interferometer(
+ waveform_polarizations=signal_polarizations,
+ interferometer=interferometer,
+ )
+ f = interferometer.frequency_array
+ duplicated_r0, duplicated_r1, duplicated_fm = np.zeros((3, f.shape[0]), dtype=complex)
+
+ for i in range(self.number_of_bins):
+ idxs = slice(self.bin_inds[i], self.bin_inds[i + 1])
+ duplicated_fm[idxs] = self.bin_centers[i]
+ duplicated_r0[idxs] = r0[i]
+ duplicated_r1[idxs] = r1[i]
+
+ full_waveform_ratio = duplicated_r0 + duplicated_r1 * (f - duplicated_fm)
+ return fiducial_waveform * full_waveform_ratio
+
+ def calculate_snrs(self, waveform_polarizations, interferometer):
+ r0, r1 = self.compute_waveform_ratio_per_interferometer(
+ waveform_polarizations=waveform_polarizations,
+ interferometer=interferometer,
+ )
+ a0, a1, b0, b1 = self.summary_data[interferometer.name]
+ d_inner_h = np.sum(a0 * np.conjugate(r0) + a1 * np.conjugate(r1))
+ h_inner_h = np.sum(b0 * np.abs(r0) ** 2 + 2 * b1 * np.real(r0 * np.conjugate(r1)))
+ optimal_snr_squared = h_inner_h
+ complex_matched_filter_snr = d_inner_h / (optimal_snr_squared ** 0.5)
+
+ if self.time_marginalization:
+ full_waveform = self._compute_full_waveform(
+ signal_polarizations=waveform_polarizations,
+ interferometer=interferometer,
+ )
+ d_inner_h_array = 4 / self.waveform_generator.duration * np.fft.fft(
+ full_waveform[0:-1]
+ * interferometer.frequency_domain_strain.conjugate()[0:-1]
+ / interferometer.power_spectral_density_array[0:-1])
+
+ else:
+ d_inner_h_array = None
+
+ return self._CalculatedSNRs(
+ d_inner_h=d_inner_h, optimal_snr_squared=optimal_snr_squared,
+ complex_matched_filter_snr=complex_matched_filter_snr,
+ d_inner_h_array=d_inner_h_array, optimal_snr_squared_array=None,
+ d_inner_h_squared_tc_array=None)
diff --git a/bilby/gw/source.py b/bilby/gw/source.py
index 4d7793b92fe8124e6e0c3f4397df1bfacc3ca91d..e13abb9e21935f7637f26bef261d96731bec982f 100644
--- a/bilby/gw/source.py
+++ b/bilby/gw/source.py
@@ -444,6 +444,78 @@ def binary_neutron_star_roq(
phi_12=phi_12, lambda_1=lambda_1, lambda_2=lambda_2, **waveform_kwargs)
+def lal_binary_black_hole_relative_binning(
+ 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):
+ """ Source model to go with RelativeBinningGravitationalWaveTransient likelihood.
+
+ All parameters are the same as in the usual source models, except `fiducial`
+
+ fiducial: float
+ If fiducial=1, waveform evaluated on the full frequency grid is returned.
+ If fiducial=0, waveform evaluated at waveform_kwargs["frequency_bin_edges"]
+ is returned.
+ """
+
+ 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:
+ waveform_kwargs["frequencies"] = waveform_kwargs.pop("frequency_bin_edges")
+ return _base_waveform_frequency_sequence(
+ 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_relative_binning(
+ 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):
+ """ Source model to go with RelativeBinningGravitationalWaveTransient likelihood.
+
+ All parameters are the same as in the usual source models, except `fiducial`
+
+ fiducial: float
+ If fiducial=1, waveform evaluated on the full frequency grid is returned.
+ If fiducial=0, waveform evaluated at waveform_kwargs["frequency_bin_edges"]
+ is returned.
+ """
+
+ 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:
+ waveform_kwargs["frequencies"] = waveform_kwargs.pop("frequency_bin_edges")
+ return _base_waveform_frequency_sequence(
+ 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_roq_waveform(
frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
phi_12, a_2, tilt_2, lambda_1, lambda_2, phi_jl, theta_jn, phase,
diff --git a/examples/gw_examples/data_examples/GW190425.py b/examples/gw_examples/data_examples/GW190425.py
index b575a4ed88842cae9f0b4683878e310e8197e654..8d9fb02af9e3de6e32b56f5c50fd0eeaa73d0809 100644
--- a/examples/gw_examples/data_examples/GW190425.py
+++ b/examples/gw_examples/data_examples/GW190425.py
@@ -1,11 +1,11 @@
#!/usr/bin/env python
"""
-Tutorial to demonstrate running parameter estimation on GW190425
+Tutorial to demonstrate running parameter estimation on GW190425.
This example estimates all 17 parameters of the binary neutron star system using
-commonly used prior distributions. This will take several hours to run. The
-data is obtained using gwpy, see [1] for information on how to access data on
-the LIGO Data Grid instead.
+commonly used prior distributions. We shall use the relative binning likelihood.
+This will take around an hour to run. The data is obtained using gwpy,
+see [1] for information on how to access data on the LIGO Data Grid instead.
[1] https://gwpy.github.io/docs/stable/timeseries/remote-access.html
"""
@@ -16,7 +16,7 @@ logger = bilby.core.utils.logger
outdir = "outdir"
label = "GW190425"
-# Note you can get trigger times using the gwosc package, e.g.:
+# Note you can get trigger times using the gwosc package, e.g.
# > from gwosc import datasets
# > datasets.event_gps("GW190425")
trigger_time = 1240215503.0
@@ -29,6 +29,29 @@ post_trigger_duration = 2 # Time between trigger time and end of segment
end_time = trigger_time + post_trigger_duration
start_time = end_time - duration
+# The fiducial parameters are taken to me the max likelihood sample from the
+# posterior sample release of LIGO-Virgo
+# https://www.gw-openscience.org/eventapi/html/O3_Discovery_Papers/GW190425/
+
+fiducial_parameters = {
+ "a_1": 0.018,
+ "a_2": 0.016,
+ "chirp_mass": 1.48658,
+ "dec": 0.438,
+ "geocent_time": 1240215503.039,
+ "lambda_1": 446.941,
+ "lambda_2": 43.386,
+ "luminosity_distance": 206.751,
+ "mass_ratio": 0.8955,
+ "phase": 3.0136566567608765,
+ "phi_12": 4.319,
+ "phi_jl": 5.07,
+ "psi": 0.281,
+ "ra": 4.2,
+ "theta_jn": 0.185,
+ "tilt_1": 0.879,
+ "tilt_2": 0.514,
+}
psd_duration = 1024
psd_start_time = start_time - psd_duration
psd_end_time = start_time
@@ -64,6 +87,7 @@ ifo_list.plot_data(outdir=outdir, label=label)
# You can overwrite this using the syntax below in the file,
# or choose a fixed value by just providing a float value as the prior.
priors = bilby.gw.prior.BBHPriorDict(filename="GW190425.prior")
+priors["fiducial"] = 0
# Add the geocent time prior
priors["geocent_time"] = bilby.core.prior.Uniform(
@@ -76,7 +100,7 @@ priors["geocent_time"] = bilby.core.prior.Uniform(
# the waveform approximant and reference frequency and a parameter conversion
# which allows us to sample in chirp mass and ratio rather than component mass
waveform_generator = bilby.gw.WaveformGenerator(
- frequency_domain_source_model=bilby.gw.source.lal_binary_neutron_star,
+ frequency_domain_source_model=bilby.gw.source.lal_binary_neutron_star_relative_binning,
parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_neutron_star_parameters,
waveform_arguments={
"waveform_approximant": "IMRPhenomPv2_NRTidalv2",
@@ -87,13 +111,14 @@ waveform_generator = bilby.gw.WaveformGenerator(
# In this step, we define the likelihood. Here we use the standard likelihood
# function, passing it the data and the waveform generator.
# Note, phase_marginalization is formally invalid with a precessing waveform such as IMRPhenomPv2
-likelihood = bilby.gw.likelihood.GravitationalWaveTransient(
+likelihood = bilby.gw.likelihood.RelativeBinningGravitationalWaveTransient(
ifo_list,
waveform_generator,
priors=priors,
- time_marginalization=True,
- phase_marginalization=False,
+ time_marginalization=False,
+ phase_marginalization=True,
distance_marginalization=True,
+ fiducial_parameters=fiducial_parameters,
)
# Finally, we run the sampler. This function takes the likelihood and prior
@@ -108,6 +133,6 @@ result = bilby.run_sampler(
check_point_delta_t=600,
check_point_plot=True,
npool=1,
- conversion_function=bilby.gw.conversion.generate_all_bbh_parameters,
+ conversion_function=bilby.gw.conversion.generate_all_bns_parameters,
)
result.plot_corner()
diff --git a/examples/gw_examples/injection_examples/relative_binning.py b/examples/gw_examples/injection_examples/relative_binning.py
new file mode 100644
index 0000000000000000000000000000000000000000..5a848f003b0a88262d0380272c85b624ef306111
--- /dev/null
+++ b/examples/gw_examples/injection_examples/relative_binning.py
@@ -0,0 +1,162 @@
+#!/usr/bin/env python
+"""
+Tutorial to demonstrate running parameter estimation on a reduced parameter
+space for an injected signal, using the relative binning likelihood.
+
+This example estimates the masses using a uniform prior in both component masses
+and distance using a uniform in comoving volume prior on luminosity distance
+between luminosity distances of 100Mpc and 5Gpc, the cosmology is Planck15.
+"""
+
+import bilby
+import numpy as np
+from tqdm.auto import trange
+
+# Set the duration and sampling frequency of the data segment that we're
+# going to inject the signal into
+duration = 4.0
+sampling_frequency = 2048.0
+minimum_frequency = 20
+
+# Specify the output directory and the name of the simulation.
+outdir = "outdir"
+label = "relative"
+bilby.core.utils.setup_logger(outdir=outdir, label=label)
+
+# Set up a random seed for result reproducibility. This is optional!
+np.random.seed(88170235)
+
+# We are going to inject a binary black hole waveform. We first establish a
+# dictionary of parameters that includes all of the different waveform
+# parameters, including masses of the two black holes (mass_1, mass_2),
+# spins of both black holes (a, tilt, phi), etc.
+injection_parameters = dict(
+ mass_1=36.0,
+ mass_2=29.0,
+ a_1=0.4,
+ a_2=0.3,
+ tilt_1=0.5,
+ tilt_2=1.0,
+ phi_12=1.7,
+ phi_jl=0.3,
+ luminosity_distance=2000.0,
+ theta_jn=0.4,
+ psi=2.659,
+ phase=1.3,
+ geocent_time=1126259642.413,
+ ra=1.375,
+ dec=-1.2108,
+ fiducial=1,
+)
+
+# Fixed arguments passed into the source model
+waveform_arguments = dict(
+ waveform_approximant="IMRPhenomXP",
+ reference_frequency=50.0,
+ minimum_frequency=minimum_frequency,
+)
+
+# Create the waveform_generator using a LAL BinaryBlackHole source function
+waveform_generator = bilby.gw.WaveformGenerator(
+ duration=duration,
+ sampling_frequency=sampling_frequency,
+ # frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
+ frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole_relative_binning,
+ parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,
+ waveform_arguments=waveform_arguments,
+)
+
+# Set up interferometers. In this case we'll use two interferometers
+# (LIGO-Hanford (H1), LIGO-Livingston (L1). These default to their design
+# sensitivity
+ifos = bilby.gw.detector.InterferometerList(["H1", "L1"])
+ifos.set_strain_data_from_power_spectral_densities(
+ sampling_frequency=sampling_frequency,
+ duration=duration,
+ start_time=injection_parameters["geocent_time"] - 2,
+)
+ifos.inject_signal(
+ waveform_generator=waveform_generator, parameters=injection_parameters
+)
+
+# Set up a PriorDict, which inherits from dict.
+# By default we will sample all terms in the signal models. However, this will
+# take a long time for the calculation, so for this example we will set almost
+# all of the priors to be equall to their injected values. This implies the
+# prior is a delta function at the true, injected value. In reality, the
+# sampler implementation is smart enough to not sample any parameter that has
+# a delta-function prior.
+# The above list does *not* include mass_1, mass_2, theta_jn and luminosity
+# distance, which means those are the parameters that will be included in the
+# sampler. If we do nothing, then the default priors get used.
+priors = bilby.gw.prior.BBHPriorDict()
+for key in [
+ "a_1",
+ "a_2",
+ "tilt_1",
+ "tilt_2",
+ "phi_12",
+ "phi_jl",
+ "psi",
+ "ra",
+ "dec",
+ "geocent_time",
+ "phase",
+]:
+ priors[key] = injection_parameters[key]
+priors["fiducial"] = 0
+
+# Perform a check that the prior does not extend to a parameter space longer than the data
+priors.validate_prior(duration, minimum_frequency)
+
+# Initialise the likelihood by passing in the interferometer data (ifos) and
+# the waveform generator
+likelihood = bilby.gw.likelihood.RelativeBinningGravitationalWaveTransient(
+ interferometers=ifos,
+ waveform_generator=waveform_generator,
+ priors=priors,
+ distance_marginalization=True,
+ fiducial_parameters=injection_parameters,
+)
+
+# Run sampler. In this case we're going to use the `dynesty` sampler
+result = bilby.run_sampler(
+ likelihood=likelihood,
+ priors=priors,
+ sampler="nestle",
+ npoints=100,
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label=label,
+ conversion_function=bilby.gw.conversion.generate_all_bbh_parameters,
+)
+
+alt_waveform_generator = bilby.gw.WaveformGenerator(
+ duration=duration,
+ sampling_frequency=sampling_frequency,
+ frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
+ # frequency_domain_source_model=lal_binary_black_hole_relative_binning,
+ parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters,
+ waveform_arguments=waveform_arguments,
+)
+alt_likelihood = bilby.gw.likelihood.GravitationalWaveTransient(
+ interferometers=ifos,
+ waveform_generator=alt_waveform_generator,
+)
+likelihood.distance_marginalization = False
+weights = list()
+for ii in trange(len(result.posterior)):
+ parameters = dict(result.posterior.iloc[ii])
+ likelihood.parameters.update(parameters)
+ alt_likelihood.parameters.update(parameters)
+ weights.append(
+ alt_likelihood.log_likelihood_ratio() - likelihood.log_likelihood_ratio()
+ )
+weights = np.exp(weights)
+print(
+ f"Reweighting efficiency is {np.mean(weights)**2 / np.mean(weights**2) * 100:.2f}%"
+)
+print(f"Binned vs unbinned log Bayes factor {np.log(np.mean(weights)):.2f}")
+
+# Make a corner plot.
+# result.plot_corner()
diff --git a/test/gw/likelihood/relative_binning_test.py b/test/gw/likelihood/relative_binning_test.py
new file mode 100644
index 0000000000000000000000000000000000000000..2f1352e4174bad7cafe3854644968b91807c40b9
--- /dev/null
+++ b/test/gw/likelihood/relative_binning_test.py
@@ -0,0 +1,166 @@
+import unittest
+from copy import deepcopy
+
+import bilby
+import numpy as np
+from parameterized import parameterized
+
+
+class TestRelativeBinningLikelihood(unittest.TestCase):
+ def setUp(self):
+ duration = 16
+ fmin = 20
+ sampling_frequency = 8192
+ self.test_parameters = dict(
+ chirp_mass=13,
+ mass_ratio=0.5,
+ a_1=0.3,
+ a_2=0.4,
+ tilt_1=1.0,
+ tilt_2=0.2,
+ phi_12=1.0,
+ phi_jl=2.0,
+ luminosity_distance=2000.0,
+ theta_jn=0.4,
+ psi=0.659,
+ phase=1.3,
+ geocent_time=1187008882,
+ ra=1.3,
+ dec=-1.2,
+ )
+
+ ifos = bilby.gw.detector.InterferometerList(["H1", "L1", "V1"])
+ # np.random.seed(170817)
+ ifos.set_strain_data_from_power_spectral_densities(
+ sampling_frequency=sampling_frequency, duration=duration,
+ start_time=self.test_parameters['geocent_time'] - duration + 2.
+ )
+ for ifo in ifos:
+ ifo.minimum_frequency = fmin
+
+ spline_calibration_nodes = 10
+ self.calibration_parameters = {}
+ for ifo in ifos:
+ ifo.calibration_model = bilby.gw.calibration.CubicSpline(
+ prefix=f"recalib_{ifo.name}_",
+ minimum_frequency=ifo.minimum_frequency,
+ maximum_frequency=ifo.maximum_frequency,
+ n_points=spline_calibration_nodes
+ )
+ for i in range(spline_calibration_nodes):
+ self.test_parameters[f"recalib_{ifo.name}_amplitude_{i}"] = 0
+ self.test_parameters[f"recalib_{ifo.name}_phase_{i}"] = 0
+ # Calibration errors of 5% in amplitude and 5 degrees in phase
+ self.calibration_parameters[f"recalib_{ifo.name}_amplitude_{i}"] = \
+ np.random.normal(loc=0, scale=0.05)
+ self.calibration_parameters[f"recalib_{ifo.name}_phase_{i}"] = \
+ np.random.normal(loc=0, scale=5 * np.pi / 180)
+
+ priors = bilby.gw.prior.BBHPriorDict()
+ priors.pop("mass_1")
+ priors.pop("mass_2")
+ # priors["chirp_mass"] = bilby.core.prior.Uniform(1.29, 1.31)
+ priors["chirp_mass"] = bilby.core.prior.Uniform(12, 14)
+ priors["mass_ratio"] = bilby.core.prior.Uniform(0.125, 1)
+ priors["geocent_time"] = bilby.core.prior.Uniform(
+ self.test_parameters['geocent_time'] - 0.1,
+ self.test_parameters['geocent_time'] + 0.1)
+ priors["fiducial"] = 0
+ self.priors = priors
+
+ approximant = "IMRPhenomXP"
+ non_bin_wfg = bilby.gw.WaveformGenerator(
+ duration=duration, sampling_frequency=sampling_frequency,
+ frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
+ waveform_arguments=dict(
+ reference_frequency=fmin, minimum_frequency=fmin, approximant=approximant)
+ )
+ bin_wfg = bilby.gw.waveform_generator.WaveformGenerator(
+ duration=duration, sampling_frequency=sampling_frequency,
+ frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole_relative_binning,
+ waveform_arguments=dict(
+ reference_frequency=fmin, approximant=approximant, minimum_frequency=fmin)
+ )
+ ifos.inject_signal(
+ parameters=self.test_parameters,
+ waveform_generator=non_bin_wfg,
+ raise_error=False,
+ )
+ self.ifos = ifos
+
+ self.non_bin = bilby.gw.likelihood.GravitationalWaveTransient(
+ interferometers=ifos, waveform_generator=deepcopy(non_bin_wfg),
+ priors=priors.copy()
+ )
+ self.binned = bilby.gw.likelihood.RelativeBinningGravitationalWaveTransient(
+ interferometers=ifos, waveform_generator=deepcopy(bin_wfg),
+ fiducial_parameters=self.test_parameters,
+ priors=priors.copy(),
+ epsilon=0.05,
+ # chi=0.2,
+ )
+ self.non_bin.parameters.update(self.test_parameters)
+ self.reference_ln_l = self.non_bin.log_likelihood_ratio()
+ self.bin_wfg = bin_wfg
+ self.priors = priors
+
+ def tearDown(self):
+ del (
+ self.non_bin,
+ self.binned,
+ )
+
+ def test_matches_non_binned_many(self):
+ for _ in range(100):
+ self.non_bin.parameters.update(self.priors.sample())
+ self.binned.parameters.update(self.priors.sample())
+ regular_ln_l = self.non_bin.log_likelihood_ratio()
+ binned_ln_l = self.binned.log_likelihood_ratio()
+ if regular_ln_l - self.reference_ln_l < -20:
+ continue
+ self.assertLess(
+ abs(regular_ln_l - binned_ln_l)
+ / self.reference_ln_l,
+ 1.5e-2
+ )
+
+ @parameterized.expand([(False, ), (True, )])
+ def test_matches_non_binned(self, add_cal_errors):
+ self.non_bin.parameters.update(self.test_parameters)
+ self.binned.parameters.update(self.test_parameters)
+ if add_cal_errors:
+ self.non_bin.parameters.update(self.calibration_parameters)
+ self.binned.parameters.update(self.calibration_parameters)
+ self.assertLess(
+ abs(self.non_bin.log_likelihood_ratio() - self.binned.log_likelihood_ratio())
+ / self.reference_ln_l,
+ 1.5e-2
+ )
+
+ def test_optimization_gives_good_match(self):
+ fiducial_parameters = self.test_parameters.copy()
+ fiducial_parameters["chirp_mass"] *= 0.99
+ priors = self.priors.copy()
+ for key in [
+ "ra", "dec", "geocent_time", "phase", "psi", "theta_jn", "luminosity_distance",
+ "a_1", "a_2", "tilt_1", "tilt_2", "phi_12", "phi_jl",
+ ]:
+ priors[key] = self.test_parameters[key]
+ binned = bilby.gw.likelihood.RelativeBinningGravitationalWaveTransient(
+ interferometers=self.ifos, waveform_generator=deepcopy(self.bin_wfg),
+ priors=priors,
+ fiducial_parameters=fiducial_parameters,
+ epsilon=0.05,
+ update_fiducial_parameters=True,
+ )
+ self.non_bin.parameters.update(self.test_parameters)
+ binned.parameters.update(self.test_parameters)
+ self.assertLess(
+ abs(self.non_bin.log_likelihood_ratio() - binned.log_likelihood_ratio())
+ / self.reference_ln_l,
+ 1.5e-2
+ )
+
+
+if __name__ == "__main__":
+ unittest.main()
diff --git a/test/gw/likelihood_test.py b/test/gw/likelihood_test.py
index 29962de2b96b5e2cf1eb0a5e97e906d7a7f871a2..e7e707c287aeb473cf3bb0528d9c78e876a39dd4 100644
--- a/test/gw/likelihood_test.py
+++ b/test/gw/likelihood_test.py
@@ -2,6 +2,7 @@ import itertools
import os
import pytest
import unittest
+from copy import deepcopy
from itertools import product
from parameterized import parameterized
@@ -433,7 +434,7 @@ class TestMarginalizations(unittest.TestCase):
The `time_jitter` parameter makes this a weaker dependence during sampling.
"""
_parameters = product(
- ["regular", "roq"],
+ ["regular", "roq", "relbin"],
["luminosity_distance", "geocent_time", "phase"],
[True, False],
[True, False],
@@ -465,6 +466,7 @@ class TestMarginalizations(unittest.TestCase):
ra=1.375,
dec=-1.2108,
time_jitter=0,
+ fiducial=0,
)
self.interferometers = bilby.gw.detector.InterferometerList(["H1"])
@@ -524,6 +526,18 @@ class TestMarginalizations(unittest.TestCase):
self.roq_linear_matrix_file = f"{roq_dir}/B_linear.npy"
self.roq_quadratic_matrix_file = f"{roq_dir}/B_quadratic.npy"
+ self.relbin_waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
+ duration=self.duration,
+ sampling_frequency=self.sampling_frequency,
+ frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole_relative_binning,
+ start_time=1126259640,
+ waveform_arguments=dict(
+ reference_frequency=20.0,
+ minimum_frequency=20.0,
+ approximant="IMRPhenomPv2",
+ )
+ )
+
def tearDown(self):
del self.duration
del self.sampling_frequency
@@ -542,7 +556,7 @@ class TestMarginalizations(unittest.TestCase):
def likelihood_kwargs(self, kind, time_marginalization, phase_marginalization, distance_marginalization, priors):
if priors is None:
- priors = self.priors.copy()
+ priors = deepcopy(self.priors)
if distance_marginalization and phase_marginalization:
lookup = TestMarginalizations.lookup_phase
elif distance_marginalization:
@@ -566,6 +580,9 @@ class TestMarginalizations(unittest.TestCase):
))
if os.path.exists(self.__class__.path_to_roq_weights):
kwargs["weights"] = self.__class__.path_to_roq_weights
+ elif kind == "relbin":
+ kwargs["fiducial_parameters"] = deepcopy(self.parameters)
+ kwargs["waveform_generator"] = self.relbin_waveform_generator
return kwargs
def get_likelihood(
@@ -583,6 +600,10 @@ class TestMarginalizations(unittest.TestCase):
cls_ = bilby.gw.likelihood.GravitationalWaveTransient
elif kind == "roq":
cls_ = bilby.gw.likelihood.ROQGravitationalWaveTransient
+ elif kind == "relbin":
+ cls_ = bilby.gw.likelihood.RelativeBinningGravitationalWaveTransient
+ kwargs["epsilon"] = 0.3
+ self.parameters["fiducial"] = 0
else:
raise ValueError(f"kind {kind} not understood")
like = cls_(**kwargs)
@@ -591,6 +612,10 @@ class TestMarginalizations(unittest.TestCase):
like.parameters = self.parameters.copy()
if time_marginalization:
like.parameters["geocent_time"] = self.interferometers.start_time
+ if distance_marginalization:
+ like.parameters["luminosity_distance"] = like._ref_dist
+ if phase_marginalization:
+ like.parameters["phase"] = 0.0
return like
def _template(self, marginalized, non_marginalized, key, prior=None, values=None):
@@ -630,7 +655,8 @@ class TestMarginalizations(unittest.TestCase):
key=key,
)
- def test_time_marginalisation_full_segment(self):
+ @parameterized.expand(["regular", "relbin"])
+ def test_time_marginalisation_full_segment(self, kind):
"""
Test time marginalised likelihood matches brute force version over
just part of a segment.
@@ -642,15 +668,15 @@ class TestMarginalizations(unittest.TestCase):
)
priors["geocent_time"] = prior
self._template(
- self.get_likelihood("regular", time_marginalization=True, priors=priors.copy()),
- self.get_likelihood("regular", priors=priors.copy()),
+ self.get_likelihood(kind, time_marginalization=True, priors=priors.copy()),
+ self.get_likelihood(kind, priors=priors.copy()),
key="geocent_time",
values=self.waveform_generator.time_array,
prior=prior,
)
@parameterized.expand(
- itertools.product(["regular", "roq"], *itertools.repeat([True, False], 3)),
+ itertools.product(["regular", "roq", "relbin"], *itertools.repeat([True, False], 3)),
name_func=lambda func, num, param: (
f"{func.__name__}_{num}__{param.args[0]}_" + "_".join([
["D", "P", "T"][ii] for ii, val
@@ -832,19 +858,6 @@ class TestROQLikelihood(unittest.TestCase):
self.roq.parameters["geocent_time"] = -5
self.assertEqual(self.roq.log_likelihood_ratio(), np.nan_to_num(-np.inf))
- def test_phase_marginalisation_roq(self):
- """Test phase marginalised likelihood matches brute force version"""
- self.non_roq_phase.parameters = self.test_parameters.copy()
- self.roq_phase.parameters = self.test_parameters.copy()
- self.assertLess(
- abs(
- self.non_roq_phase.log_likelihood_ratio()
- - self.roq_phase.log_likelihood_ratio()
- )
- / self.non_roq_phase.log_likelihood_ratio(),
- 1e-3,
- )
-
def test_create_roq_weights_with_params(self):
roq = bilby.gw.likelihood.ROQGravitationalWaveTransient(
interferometers=self.ifos,
diff --git a/test/gw/source_test.py b/test/gw/source_test.py
index 627278a337df155d09a82a342e3e80bc4bf6608f..8a4b7625b4ed2c92990f36fab855de77f13b71a5 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_relative_binning(
+ 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_relative_binning(
+ 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_relative_binning(
+ 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_relative_binning(
+ 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_relative_binning(
+ 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_relative_binning(
+ self.frequency_array, **raise_error_parameters
+ )
+
+ def test_relbin_bbh_fails_without_fiducial_option(self):
+ with self.assertRaises(TypeError):
+ bilby.gw.source.lal_binary_black_hole_relative_binning(
+ 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_relative_binning(
+ self.frequency_array, PhenomXPrecVersion=223, **self.parameters
+ )
+ out_v102 = bilby.gw.source.lal_binary_black_hole_relative_binning(
+ 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_relative_binning(
+ 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_relative_binning(
+ self.frequency_array, **self.parameters
+ ),
+ dict,
+ )
+
+ def test_relbin_bns_fails_without_fiducial_option(self):
+ with self.assertRaises(TypeError):
+ bilby.gw.source.lal_binary_neutron_star_relative_binning(
+ 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_relative_binning(
+ 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_relative_binning(
+ self.frequency_array, **self.parameters
+ )
+
+
if __name__ == "__main__":
unittest.main()