From 5014f4d29aa31a3906758b0fd1d7418e8bc1975b Mon Sep 17 00:00:00 2001 From: Colm Talbot Date: Tue, 23 Dec 2025 11:59:13 -0500 Subject: [PATCH 01/16] FEAT: add option to crop whitened time series --- bilby/core/series.py | 4 ++ bilby/gw/detector/interferometer.py | 92 +++++++++++++++++++------ bilby/gw/detector/strain_data.py | 65 ++++++++++++++++- bilby/gw/likelihood/base.py | 54 ++++++--------- bilby/gw/likelihood/basic.py | 14 ++-- test/core/series_test.py | 4 ++ test/gw/detector/interferometer_test.py | 21 ++++-- test/gw/detector/strain_data_test.py | 19 +++++ test/gw/likelihood_test.py | 33 +++++++-- 9 files changed, 235 insertions(+), 71 deletions(-) diff --git a/bilby/core/series.py b/bilby/core/series.py index ba1d0ffcb..fd7d36f93 100644 --- a/bilby/core/series.py +++ b/bilby/core/series.py @@ -128,3 +128,7 @@ def start_time(self): def start_time(self, start_time): self._start_time = start_time self._time_array_updated = False + + @property + def end_time(self) + return self.start_time + self.duration diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py index 94267d30c..8d7a96d73 100644 --- a/bilby/gw/detector/interferometer.py +++ b/bilby/gw/detector/interferometer.py @@ -9,6 +9,7 @@ from ...core import utils from ...core.utils import docstring, logger, PropertyAccessor, safe_file_dump +from ...core.utils.series import infft, nfft from ...core.utils.env import string_to_boolean from .. import utils as gwutils from .calibration import Recalibrate @@ -43,6 +44,9 @@ class Interferometer(object): minimum_frequency = PropertyAccessor('strain_data', 'minimum_frequency') maximum_frequency = PropertyAccessor('strain_data', 'maximum_frequency') frequency_mask = PropertyAccessor('strain_data', 'frequency_mask') + time_mask = PropertyAccessor('strain_data', 'time_mask') + crop_duration = PropertyAccessor('strain_data', 'crop_duration') + cropped_duration = PropertyAccessor('strain_data', 'cropped_duration') frequency_domain_strain = PropertyAccessor('strain_data', 'frequency_domain_strain') time_domain_strain = PropertyAccessor('strain_data', 'time_domain_strain') @@ -599,10 +603,7 @@ def optimal_snr_squared(self, signal): ======= float: The optimal signal to noise ratio possible squared """ - return gwutils.optimal_snr_squared( - signal=signal[self.strain_data.frequency_mask], - power_spectral_density=self.power_spectral_density_array[self.strain_data.frequency_mask], - duration=self.strain_data.duration) + return (abs(self.whiten_frequency_series(signal))**2).sum() def inner_product(self, signal): """ @@ -616,11 +617,9 @@ def inner_product(self, signal): ======= float: The optimal signal to noise ratio possible squared """ - return gwutils.noise_weighted_inner_product( - aa=signal[self.strain_data.frequency_mask], - bb=self.strain_data.frequency_domain_strain[self.strain_data.frequency_mask], - power_spectral_density=self.power_spectral_density_array[self.strain_data.frequency_mask], - duration=self.strain_data.duration) + whitened_signal = self.whiten_frequency_series(signal) + whitened_data = self.whitened_frequency_domain_strain + return (whitened_signal.T * whitened_data.conj()).sum() def template_template_inner_product(self, signal_1, signal_2): """A noise weighted inner product between two templates, using this ifo's PSD. @@ -636,11 +635,9 @@ def template_template_inner_product(self, signal_1, signal_2): ======= float: The noise weighted inner product of the two templates """ - return gwutils.noise_weighted_inner_product( - aa=signal_1[self.strain_data.frequency_mask], - bb=signal_2[self.strain_data.frequency_mask], - power_spectral_density=self.power_spectral_density_array[self.strain_data.frequency_mask], - duration=self.strain_data.duration) + whitened_1 = self.whiten_frequency_series(signal_1) + whitened_2 = self.whiten_frequency_series(signal_2) + return (whitened_1 * whitened_2.conj()).sum() def matched_filter_snr(self, signal): """ @@ -655,11 +652,7 @@ def matched_filter_snr(self, signal): complex: The matched filter signal to noise ratio """ - return gwutils.matched_filter_snr( - signal=signal[self.strain_data.frequency_mask], - frequency_domain_strain=self.strain_data.frequency_domain_strain[self.strain_data.frequency_mask], - power_spectral_density=self.power_spectral_density_array[self.strain_data.frequency_mask], - duration=self.strain_data.duration) + return self.inner_product(signal) / self.optimal_snr_squared(signal)**0.5 def whiten_frequency_series(self, frequency_series : np.array) -> np.array: """Whitens a frequency series with the noise properties of the detector @@ -680,7 +673,62 @@ def whiten_frequency_series(self, frequency_series : np.array) -> np.array: frequency_series : np.array The frequency series, whitened by the ASD """ - return frequency_series / (self.amplitude_spectral_density_array * np.sqrt(self.duration / 4)) + if self.crop_duration == 0: + return np.nan_to_num(frequency_series / ( + self.amplitude_spectral_density_array + * np.sqrt(self.duration / 4) + * self.frequency_mask + )) + else: + return self._whiten_and_crop(frequency_series=frequency_series) + + def _whiten_and_crop(self, frequency_series : np.array) -> np.array: + """ + Whitens a frequency series with the noise properties and applies + the time mask to the whitened time-domain strain. + + First, we naively whiten the data in the frequency domain and apply + our frequency mask :math:`\\tilde{m}(f)` + + .. math:: + \\tilde{a}_w(f) = \\tilde{a}(f) \\tilde{m}(f) / \\sqrt{S_n(f)}. + + We then inverse discrete Fourier transform to get the whitened + time-domain strain :math:`w(t)` and apply the time-domain mask + :math:`m(t)`. We then perform a final discrete Fourier transform. + + .. math:: + \\bar{a}_{w}(f) = \\mathcal{F}(\\mathcal{iF}(\\tilde{a}_w(f))[m(t)]) + + Finally, we normalize the whitened strain by a factor :math:`N` + + .. math:: + N = \\sqrt{\\frac{4}{T_{c}}} + + where :math:`T_{c}` is the cropped duration such that + + .. math:: + Var(n) = \\frac{1}{N} \\sum_{k=0}^N n_W(f_k)n_W^*(f_k) = 2. + + Where the factor of two is due to the independent real and imaginary + components. + + + Parameters + ========== + frequency_series : np.array + The frequency series, whitened by the ASD + """ + initial_white = np.nan_to_num( + frequency_series + * self.frequency_mask + / self.amplitude_spectral_density_array + ) + time_series = infft(initial_white, self.sampling_frequency) + cropped_time_series = time_series[self.time_mask] + cropped_white, _ = nfft(cropped_time_series, self.sampling_frequency) + + return cropped_white * (4 / (self.duration - 2 * self.cropped_duration))**0.5 def get_whitened_time_series_from_whitened_frequency_series( self, @@ -718,7 +766,9 @@ def get_whitened_time_series_from_whitened_frequency_series( whitened_time_series = ( np.fft.irfft(whitened_frequency_series) - * np.sqrt(np.sum(self.frequency_mask)) / frequency_window_factor + * np.sqrt(np.sum(self.frequency_mask)) + / frequency_window_factor + * self.time_mask.mean()**0.5 ) return whitened_time_series diff --git a/bilby/gw/detector/strain_data.py b/bilby/gw/detector/strain_data.py index bca7acced..27ca5e03f 100644 --- a/bilby/gw/detector/strain_data.py +++ b/bilby/gw/detector/strain_data.py @@ -12,11 +12,12 @@ class InterferometerStrainData(object): duration = PropertyAccessor('_times_and_frequencies', 'duration') sampling_frequency = PropertyAccessor('_times_and_frequencies', 'sampling_frequency') start_time = PropertyAccessor('_times_and_frequencies', 'start_time') + end_time = PropertyAccessor('_times_and_frequencies', 'end_time') frequency_array = PropertyAccessor('_times_and_frequencies', 'frequency_array') time_array = PropertyAccessor('_times_and_frequencies', 'time_array') def __init__(self, minimum_frequency=0, maximum_frequency=np.inf, - roll_off=0.2, notch_list=None): + roll_off=0.2, notch_list=None, crop_duration=0): """ Initiate an InterferometerStrainData object The initialised object contains no data, this should be added using one @@ -33,6 +34,11 @@ def __init__(self, minimum_frequency=0, maximum_frequency=np.inf, This corresponds to alpha * duration / 2 for scipy tukey window. notch_list: bilby.gw.detector.strain_data.NotchList A list of notches + crop_duration: float | tuple + The duration of data to crop at the beginning/end of the segment + to avoid whitening artifacts. If a float, that duration is excluded + at each end, if a tuple, this specifies the truncation duration + at the beginning and end. """ @@ -41,11 +47,14 @@ def __init__(self, minimum_frequency=0, maximum_frequency=np.inf, self.notch_list = notch_list self.roll_off = roll_off self.window_factor = 1 + self._crop_duration = crop_duration self._times_and_frequencies = CoupledTimeAndFrequencySeries() self._frequency_mask_updated = False self._frequency_mask = None + self._time_mask_updated = False + self._time_mask = None self._frequency_domain_strain = None self._time_domain_strain = None self._channel = None @@ -135,6 +144,33 @@ def notch_list(self, notch_list): raise ValueError("notch_list {} not understood".format(notch_list)) self._frequency_mask_updated = False + @property + def crop_duration(self): + """ + The duration of data to crop at the beginning/end of the segment + to avoid conditioning artifacts. If a float, that duration is + excluded at each end, if a tuple, this specifies the truncation + duration at the beginning and end. + """ + return self._crop_duration + + @crop_duration.setter + def crop_duration(self, crop_duration): + if not isinstance(self.crop_duration, (float, int, list, tuple)): + raise TypeError(f"Invalid crop specification {self.crop_duration}") + self._crop_duration = crop_duration + self._time_mask_updated = False + + @property + def cropped_duration(self): + """ + The duration after applying the time-domain mask. + """ + if isinstance(self.crop_duration, (float, int)): + return self.duration - 2 * self.crop_duration + else: + return self.duration - sum(self.crop_duration[:2]) + @property def frequency_mask(self): """ Masking array for limiting the frequency band. @@ -159,6 +195,33 @@ def frequency_mask(self, mask): self._frequency_mask = mask self._frequency_mask_updated = True + @property + def time_mask(self): + """ Masking array for cropping corrupted data at the edges. + + Returns + ======= + mask: np.ndarray + An array of boolean values + """ + if not self._time_mask_updated: + if isinstance(self.crop_duration, (tuple, list)): + crop_start, crop_end = self.crop_duration + elif isinstance(self.crop_duration, (float, int)): + crop_start = crop_end = self.crop_duration + + time_array = self._times_and_frequencies.time_array + mask = ((time_array > self.start_time + crop_start) & + (time_array <= self.end_time - crop_end)) + self._time_mask = mask + self._time_mask_updated = True + return self._time_mask + + @time_mask.setter + def time_mask(self, mask): + self._time_mask = mask + self._time_mask_updated = True + @property def alpha(self): return 2 * self.roll_off / self.duration diff --git a/bilby/gw/likelihood/base.py b/bilby/gw/likelihood/base.py index b88b788ac..53c014a87 100644 --- a/bilby/gw/likelihood/base.py +++ b/bilby/gw/likelihood/base.py @@ -281,60 +281,51 @@ def calculate_snrs(self, waveform_polarizations, interferometer, return_array=Tr _mask = interferometer.frequency_mask if 'recalib_index' in parameters: - signal[_mask] *= self.calibration_draws[interferometer.name][int(parameters['recalib_index'])] + signal *= self.calibration_draws[interferometer.name][int(parameters['recalib_index'])] - d_inner_h = interferometer.inner_product(signal=signal) - optimal_snr_squared = interferometer.optimal_snr_squared(signal=signal) + whitened_signal = interferometer.whiten_frequency_series(signal) + + d_inner_h = (interferometer.whitened_frequency_domain_strain.conjugate() * whitened_signal).sum() + optimal_snr_squared = (abs(whitened_signal)**2).sum() complex_matched_filter_snr = d_inner_h / (optimal_snr_squared**0.5) d_inner_h_array = None optimal_snr_squared_array = None - normalization = 4 / self.waveform_generator.duration - if return_array is False: d_inner_h_array = None optimal_snr_squared_array = None elif self.time_marginalization and self.calibration_marginalization: d_inner_h_integrand = np.tile( - interferometer.frequency_domain_strain.conjugate() * signal / - interferometer.power_spectral_density_array, (self.number_of_response_curves, 1)).T + interferometer.whitened_frequency_domain_strain.conjugate() * whitened_signal, + (self.number_of_response_curves, 1) + ).T - d_inner_h_integrand[_mask] *= self.calibration_draws[interferometer.name].T + d_inner_h_integrand *= self.calibration_draws[interferometer.name].T - d_inner_h_array = 4 / self.waveform_generator.duration * np.fft.fft( - d_inner_h_integrand[0:-1], axis=0 - ).T + d_inner_h_array = np.fft.fft(d_inner_h_integrand[0:-1], axis=0).T - optimal_snr_squared_integrand = ( - normalization * np.abs(signal)**2 / interferometer.power_spectral_density_array - ) + optimal_snr_squared_integrand = np.abs(whitened_signal)**2 optimal_snr_squared_array = np.dot( - optimal_snr_squared_integrand[_mask], + optimal_snr_squared_integrand, self.calibration_abs_draws[interferometer.name].T ) elif self.time_marginalization and not self.calibration_marginalization: - d_inner_h_array = normalization * np.fft.fft( - signal[0:-1] - * interferometer.frequency_domain_strain.conjugate()[0:-1] - / interferometer.power_spectral_density_array[0:-1] + d_inner_h_array = np.fft.fft( + whitened_signal + * interferometer.whitened_frequency_domain_strain.conjugate() ) elif self.calibration_marginalization and ('recalib_index' not in parameters): d_inner_h_integrand = ( - normalization * - interferometer.frequency_domain_strain.conjugate() * signal - / interferometer.power_spectral_density_array + interferometer.whitened_frequency_domain_strain.conjugate() * whitened_signal ) - d_inner_h_array = np.dot(d_inner_h_integrand[_mask], self.calibration_draws[interferometer.name].T) + d_inner_h_array = np.dot(d_inner_h_integrand, self.calibration_draws[interferometer.name].T) - optimal_snr_squared_integrand = ( - normalization * np.abs(signal)**2 / interferometer.power_spectral_density_array - ) + optimal_snr_squared_integrand = np.abs(whitened_signal)**2 optimal_snr_squared_array = np.dot( - optimal_snr_squared_integrand[_mask], self.calibration_abs_draws[interferometer.name].T ) @@ -391,12 +382,9 @@ def priors(self, priors): def _calculate_noise_log_likelihood(self): log_l = 0 for interferometer in self.interferometers: - mask = interferometer.frequency_mask - log_l -= noise_weighted_inner_product( - interferometer.frequency_domain_strain[mask], - interferometer.frequency_domain_strain[mask], - interferometer.power_spectral_density_array[mask], - self.waveform_generator.duration) / 2 + log_l -= ( + abs(interferometer.whitened_frequency_domain_strain)**2 + ).sum() / 2 return float(np.real(log_l)) def noise_log_likelihood(self): diff --git a/bilby/gw/likelihood/basic.py b/bilby/gw/likelihood/basic.py index b2f04eb69..64c46bb82 100644 --- a/bilby/gw/likelihood/basic.py +++ b/bilby/gw/likelihood/basic.py @@ -43,9 +43,9 @@ def noise_log_likelihood(self): """ log_l = 0 for interferometer in self.interferometers: - log_l -= 2. / self.waveform_generator.duration * np.sum( - abs(interferometer.frequency_domain_strain) ** 2 / - interferometer.power_spectral_density_array) + log_l -= ( + abs(interferometer.whitened_frequency_domain_strain)**2 + ).sum() / 2 return log_l.real def log_likelihood(self, parameters=None): @@ -87,8 +87,8 @@ def log_likelihood_interferometer(self, waveform_polarizations, signal_ifo = interferometer.get_detector_response( waveform_polarizations, parameters) - log_l = - 2. / self.waveform_generator.duration * np.vdot( - interferometer.frequency_domain_strain - signal_ifo, - (interferometer.frequency_domain_strain - signal_ifo) / - interferometer.power_spectral_density_array) + residual = interferometer.frequency_domain_strain - signal_ifo + white_residual = interferometer.whiten_frequency_series(residual) + + log_l = - (abs(white_residual)**2).sum() / 2 return log_l.real diff --git a/test/core/series_test.py b/test/core/series_test.py index bf1b19c43..abea7e646 100644 --- a/test/core/series_test.py +++ b/test/core/series_test.py @@ -42,6 +42,10 @@ def test_sampling_from_init(self): def test_start_time_from_init(self): self.assertEqual(self.start_time, self.series.start_time) + def test_end_time(self): + expected = self.start_time + self.duration + self.assertEqual(expected, self.series.end_time) + def test_frequency_array_type(self): self.assertIsInstance(self.series.frequency_array, np.ndarray) diff --git a/test/gw/detector/interferometer_test.py b/test/gw/detector/interferometer_test.py index cb1666320..ce69bb878 100644 --- a/test/gw/detector/interferometer_test.py +++ b/test/gw/detector/interferometer_test.py @@ -6,6 +6,7 @@ import lalsimulation import pytest from shutil import rmtree +from parameterized import parameterized_class import numpy as np @@ -596,6 +597,7 @@ def test_time_delay_vs_lal(self): @pytest.mark.flaky(reruns=3, only_rerun=["AssertionError"]) +@parameterized_class(("crop_duration",), [(0,), (4,), (16,)]) class TestInterferometerWhitenedStrain(unittest.TestCase): def setUp(self): self.duration = 64 @@ -603,6 +605,7 @@ def setUp(self): self.ifo = bilby.gw.detector.get_empty_interferometer('H1') self.ifo.set_strain_data_from_power_spectral_density( sampling_frequency=self.sampling_frequency, duration=self.duration) + self.ifo.crop_duration = self.crop_duration self.waveform_generator = bilby.gw.waveform_generator.WaveformGenerator( duration=self.duration, sampling_frequency=self.sampling_frequency, @@ -646,9 +649,17 @@ def _check_time_series_whiteness(self, time_series): std = np.std(time_series) self.assertAlmostEqual(std, 1, places=2) + @property + def frequency_mask(self): + frequencies = bilby.core.utils.series.create_frequency_series( + duration=self.ifo.cropped_duration, sampling_frequency=self.sampling_frequency + ) + return (self.ifo.minimum_frequency <= frequencies) & ( + frequencies <= self.ifo.maximum_frequency + ) + def test_frequency_domain_whitened_strain(self): - mask = self.ifo.frequency_mask - white = self.ifo.whitened_frequency_domain_strain[mask] + white = self.ifo.whitened_frequency_domain_strain[self.frequency_mask] self._check_frequency_series_whiteness(white) def test_time_domain_whitened_strain(self): @@ -666,8 +677,10 @@ def test_frequency_domain_noise_and_signal_whitening(self): ) # Whiten the template whitened_signal_ifo = self.ifo.whiten_frequency_series(signal_ifo) - mask = self.ifo.frequency_mask - white = self.ifo.whitened_frequency_domain_strain[mask] - whitened_signal_ifo[mask] + white = ( + self.ifo.whitened_frequency_domain_strain[self.frequency_mask] + - whitened_signal_ifo[self.frequency_mask] + ) self._check_frequency_series_whiteness(white) def test_time_domain_noise_and_signal_whitening(self): diff --git a/test/gw/detector/strain_data_test.py b/test/gw/detector/strain_data_test.py index 0f82a40a2..e04cfd8e0 100644 --- a/test/gw/detector/strain_data_test.py +++ b/test/gw/detector/strain_data_test.py @@ -77,6 +77,25 @@ def test_notches_frequency_mask(self): idxs = (freqs > 100) * (freqs < 101) self.assertTrue(len(freqs[idxs]) == 0) + def test_time_mask(self): + strain_data = bilby.gw.detector.InterferometerStrainData( + minimum_frequency=20, maximum_frequency=512, crop_duration=0.1) + strain_data.set_from_time_domain_strain( + time_domain_strain=np.random.normal(0, 1, 4096), + time_array=np.arange(0, 4, 4 / 4096), + ) + + # Test from init + times = strain_data.time_array[strain_data.time_mask] + self.assertTrue(all(times > 0.1)) + self.assertTrue(all(times <= 3.9)) + + # Test from update + strain_data.crop_duration = (0.5, 1) + times = strain_data.time_array[strain_data.time_mask] + self.assertTrue(all(times > 0.5)) + self.assertTrue(all(times <= 3)) + def test_set_data_fails(self): with mock.patch("bilby.core.utils.create_frequency_series") as m: m.return_value = [1, 2, 3] diff --git a/test/gw/likelihood_test.py b/test/gw/likelihood_test.py index 25a6a7b11..5f0fc6820 100644 --- a/test/gw/likelihood_test.py +++ b/test/gw/likelihood_test.py @@ -61,11 +61,27 @@ def test_noise_log_likelihood(self): -4014.1787704539474, self.likelihood.noise_log_likelihood(), 3 ) + def test_noise_log_likelihood_with_cropping(self): + """Test noise log likelihood matches precomputed value""" + self.interferometers[0].crop_duration = 1.0 + self.likelihood.noise_log_likelihood() + self.assertAlmostEqual( + -1991.4986550018828, self.likelihood.noise_log_likelihood(), 3 + ) + def test_log_likelihood(self): """Test log likelihood matches precomputed value""" self.likelihood.log_likelihood() self.assertAlmostEqual(self.likelihood.log_likelihood(), -4032.4397343470005, 3) + def test_log_likelihood_with_cropping(self): + """Test noise log likelihood matches precomputed value""" + self.interferometers[0].crop_duration = 1.0 + self.likelihood.log_likelihood() + self.assertAlmostEqual( + -2009.7686313437, self.likelihood.log_likelihood(), 3 + ) + def test_log_likelihood_ratio(self): """Test log likelihood ratio returns the correct value""" self.assertAlmostEqual( @@ -152,14 +168,21 @@ def test_log_likelihood(self): self.assertAlmostEqual(self.likelihood.log_likelihood(), -4032.4397343470005, 3) - def test_log_likelihood_ratio(self): - """Test log likelihood ratio returns the correct value""" + def test_noise_log_likelihood_with_cropping(self): + """Test noise log likelihood matches precomputed value""" + self.interferometers[0].crop_duration = 1.0 + self.likelihood.noise_log_likelihood() self.assertAlmostEqual( - self.likelihood.log_likelihood() - self.likelihood.noise_log_likelihood(), - self.likelihood.log_likelihood_ratio(), - 3, + -1991.4986550018828, self.likelihood.noise_log_likelihood(), 3 ) + def test_log_likelihood_with_cropping(self): + """Test log likelihood matches precomputed value""" + self.interferometers[0].crop_duration = 1.0 + self.likelihood.log_likelihood() + self.assertAlmostEqual(self.likelihood.log_likelihood(), + -2009.7686313436998, 3) + def test_likelihood_zero_when_waveform_is_none(self): """Test log likelihood returns np.nan_to_num(-np.inf) when the waveform is None""" From c8453d0161202608c2c01e2cb98bd06a8cbda102 Mon Sep 17 00:00:00 2001 From: Colm Talbot Date: Tue, 23 Dec 2025 13:11:26 -0500 Subject: [PATCH 02/16] FMT: fix precommits --- bilby/core/series.py | 2 +- bilby/gw/detector/strain_data.py | 2 +- bilby/gw/likelihood/base.py | 3 +-- 3 files changed, 3 insertions(+), 4 deletions(-) diff --git a/bilby/core/series.py b/bilby/core/series.py index fd7d36f93..e607b5f57 100644 --- a/bilby/core/series.py +++ b/bilby/core/series.py @@ -130,5 +130,5 @@ def start_time(self, start_time): self._time_array_updated = False @property - def end_time(self) + def end_time(self): return self.start_time + self.duration diff --git a/bilby/gw/detector/strain_data.py b/bilby/gw/detector/strain_data.py index 27ca5e03f..5d3545828 100644 --- a/bilby/gw/detector/strain_data.py +++ b/bilby/gw/detector/strain_data.py @@ -221,7 +221,7 @@ def time_mask(self): def time_mask(self, mask): self._time_mask = mask self._time_mask_updated = True - + @property def alpha(self): return 2 * self.roll_off / self.duration diff --git a/bilby/gw/likelihood/base.py b/bilby/gw/likelihood/base.py index 53c014a87..b0fdb8b4d 100644 --- a/bilby/gw/likelihood/base.py +++ b/bilby/gw/likelihood/base.py @@ -11,7 +11,7 @@ from ...core.prior import Interped, Prior, Uniform, DeltaFunction from ..detector import InterferometerList, get_empty_interferometer, calibration from ..prior import BBHPriorDict, Cosmological -from ..utils import noise_weighted_inner_product, zenith_azimuth_to_ra_dec, ln_i0 +from ..utils import zenith_azimuth_to_ra_dec, ln_i0 class GravitationalWaveTransient(Likelihood): @@ -278,7 +278,6 @@ def calculate_snrs(self, waveform_polarizations, interferometer, return_array=Tr interferometer=interferometer, parameters=parameters, ) - _mask = interferometer.frequency_mask if 'recalib_index' in parameters: signal *= self.calibration_draws[interferometer.name][int(parameters['recalib_index'])] From d102673ec5a49e3771c31663820006100d5c2d6d Mon Sep 17 00:00:00 2001 From: Colm Talbot Date: Tue, 23 Dec 2025 13:17:09 -0500 Subject: [PATCH 03/16] Expose whiten and crop method --- bilby/gw/detector/interferometer.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py index 8d7a96d73..d9b5b67eb 100644 --- a/bilby/gw/detector/interferometer.py +++ b/bilby/gw/detector/interferometer.py @@ -680,9 +680,9 @@ def whiten_frequency_series(self, frequency_series : np.array) -> np.array: * self.frequency_mask )) else: - return self._whiten_and_crop(frequency_series=frequency_series) + return self.whiten_and_crop(frequency_series=frequency_series) - def _whiten_and_crop(self, frequency_series : np.array) -> np.array: + def whiten_and_crop(self, frequency_series : np.array) -> np.array: """ Whitens a frequency series with the noise properties and applies the time mask to the whitened time-domain strain. From e6c6e3c9b85d67408f6089e9408255ad76eedaa0 Mon Sep 17 00:00:00 2001 From: Colm Talbot Date: Tue, 20 Jan 2026 16:21:32 +0000 Subject: [PATCH 04/16] BUG: squash some bugs --- bilby/gw/detector/interferometer.py | 5 ++--- bilby/gw/likelihood/base.py | 12 +++++++----- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py index d9b5b67eb..8678d1ce1 100644 --- a/bilby/gw/detector/interferometer.py +++ b/bilby/gw/detector/interferometer.py @@ -677,8 +677,7 @@ def whiten_frequency_series(self, frequency_series : np.array) -> np.array: return np.nan_to_num(frequency_series / ( self.amplitude_spectral_density_array * np.sqrt(self.duration / 4) - * self.frequency_mask - )) + )) * self.frequency_mask else: return self.whiten_and_crop(frequency_series=frequency_series) @@ -728,7 +727,7 @@ def whiten_and_crop(self, frequency_series : np.array) -> np.array: cropped_time_series = time_series[self.time_mask] cropped_white, _ = nfft(cropped_time_series, self.sampling_frequency) - return cropped_white * (4 / (self.duration - 2 * self.cropped_duration))**0.5 + return cropped_white * (4 / self.cropped_duration)**0.5 def get_whitened_time_series_from_whitened_frequency_series( self, diff --git a/bilby/gw/likelihood/base.py b/bilby/gw/likelihood/base.py index b0fdb8b4d..8d01fe690 100644 --- a/bilby/gw/likelihood/base.py +++ b/bilby/gw/likelihood/base.py @@ -301,30 +301,32 @@ def calculate_snrs(self, waveform_polarizations, interferometer, return_array=Tr (self.number_of_response_curves, 1) ).T - d_inner_h_integrand *= self.calibration_draws[interferometer.name].T + d_inner_h_integrand[interferometer.frequency_mask] *= self.calibration_draws[interferometer.name].T - d_inner_h_array = np.fft.fft(d_inner_h_integrand[0:-1], axis=0).T + d_inner_h_array = np.fft.fft(d_inner_h_integrand[:-1], axis=0).T optimal_snr_squared_integrand = np.abs(whitened_signal)**2 optimal_snr_squared_array = np.dot( - optimal_snr_squared_integrand, + optimal_snr_squared_integrand[interferometer.frequency_mask], self.calibration_abs_draws[interferometer.name].T ) elif self.time_marginalization and not self.calibration_marginalization: - d_inner_h_array = np.fft.fft( + d_inner_h_integrand = ( whitened_signal * interferometer.whitened_frequency_domain_strain.conjugate() ) + d_inner_h_array = np.fft.fft(d_inner_h_integrand[:-1]) elif self.calibration_marginalization and ('recalib_index' not in parameters): d_inner_h_integrand = ( interferometer.whitened_frequency_domain_strain.conjugate() * whitened_signal - ) + )[interferometer.frequency_mask] d_inner_h_array = np.dot(d_inner_h_integrand, self.calibration_draws[interferometer.name].T) optimal_snr_squared_integrand = np.abs(whitened_signal)**2 optimal_snr_squared_array = np.dot( + optimal_snr_squared_integrand[interferometer.frequency_mask], self.calibration_abs_draws[interferometer.name].T ) From 45fc89374229c3a2464231ec55035f72e6dbfd38 Mon Sep 17 00:00:00 2001 From: Colm Talbot Date: Tue, 20 Jan 2026 16:46:07 +0000 Subject: [PATCH 05/16] TST: remove now obsolete test --- test/gw/detector/interferometer_test.py | 21 --------------------- 1 file changed, 21 deletions(-) diff --git a/test/gw/detector/interferometer_test.py b/test/gw/detector/interferometer_test.py index ce69bb878..a2f784078 100644 --- a/test/gw/detector/interferometer_test.py +++ b/test/gw/detector/interferometer_test.py @@ -303,27 +303,6 @@ def test_inject_signal_raises_value_error(self): with self.assertRaises(ValueError): self.ifo.inject_signal(injection_polarizations=None, parameters=None) - def test_optimal_snr_squared(self): - """ - Merely checks parameters are given in the right order and the frequency - mask is applied. - """ - with mock.patch("bilby.gw.utils.noise_weighted_inner_product") as m: - m.side_effect = lambda a, b, c, d: [a, b, c, d] - signal = np.ones_like(self.ifo.power_spectral_density_array) - mask = self.ifo.frequency_mask - expected = [ - signal[mask], - signal[mask], - self.ifo.power_spectral_density_array[mask], - self.ifo.strain_data.duration, - ] - actual = self.ifo.optimal_snr_squared(signal=signal) - self.assertTrue(np.array_equal(expected[0], actual[0])) - self.assertTrue(np.array_equal(expected[1], actual[1])) - self.assertTrue(np.array_equal(expected[2], actual[2])) - self.assertEqual(expected[3], actual[3]) - def test_template_template_inner_product(self): signal_1 = np.ones_like(self.ifo.power_spectral_density_array) signal_2 = np.ones_like(self.ifo.power_spectral_density_array) * 2 From 0e185179bcb46dc9a3ee1c39261865bcaf8f63cc Mon Sep 17 00:00:00 2001 From: Colm Talbot Date: Fri, 23 Jan 2026 11:36:52 -0500 Subject: [PATCH 06/16] FEAT: add whitened frequeny mask --- bilby/gw/detector/interferometer.py | 1 + bilby/gw/detector/strain_data.py | 30 +++++++++++++++++++++---- test/gw/detector/interferometer_test.py | 15 +++---------- 3 files changed, 30 insertions(+), 16 deletions(-) diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py index 8678d1ce1..f14e58566 100644 --- a/bilby/gw/detector/interferometer.py +++ b/bilby/gw/detector/interferometer.py @@ -47,6 +47,7 @@ class Interferometer(object): time_mask = PropertyAccessor('strain_data', 'time_mask') crop_duration = PropertyAccessor('strain_data', 'crop_duration') cropped_duration = PropertyAccessor('strain_data', 'cropped_duration') + cropped_frequency_mask = PropertyAccessor('strain_data', 'cropped_frequency_mask') frequency_domain_strain = PropertyAccessor('strain_data', 'frequency_domain_strain') time_domain_strain = PropertyAccessor('strain_data', 'time_domain_strain') diff --git a/bilby/gw/detector/strain_data.py b/bilby/gw/detector/strain_data.py index 5d3545828..b26d409af 100644 --- a/bilby/gw/detector/strain_data.py +++ b/bilby/gw/detector/strain_data.py @@ -181,20 +181,42 @@ def frequency_mask(self): An array of boolean values """ if not self._frequency_mask_updated: - frequency_array = self._times_and_frequencies.frequency_array + self._update_frequency_mask() + return self._frequency_mask + + def _update_frequency_mask(self): + def calculate_frequency_mask(frequency_array): mask = ((frequency_array >= self.minimum_frequency) & (frequency_array <= self.maximum_frequency)) for notch in self.notch_list: mask[notch.get_idxs(frequency_array)] = False - self._frequency_mask = mask - self._frequency_mask_updated = True - return self._frequency_mask + return mask + + self._frequency_mask = calculate_frequency_mask( + self._times_and_frequencies.frequency_array + ) + + cropped_frequencies = bilby.core.utils.series.create_frequency_series( + duration=self.cropped_duration, + sampling_frequency=self.sampling_frequency + ) + self._cropped_frequency_mask = calculate_frequency_mask( + cropped_frequencies + ) + + self._frequency_mask_updated = True @frequency_mask.setter def frequency_mask(self, mask): self._frequency_mask = mask self._frequency_mask_updated = True + @property + def cropped_frequency_mask(self): + if not self._frequency_mask_updated: + self._update_frequency_mask + return self._cropped_frequency_mask + @property def time_mask(self): """ Masking array for cropping corrupted data at the edges. diff --git a/test/gw/detector/interferometer_test.py b/test/gw/detector/interferometer_test.py index a2f784078..6de3bca42 100644 --- a/test/gw/detector/interferometer_test.py +++ b/test/gw/detector/interferometer_test.py @@ -628,17 +628,8 @@ def _check_time_series_whiteness(self, time_series): std = np.std(time_series) self.assertAlmostEqual(std, 1, places=2) - @property - def frequency_mask(self): - frequencies = bilby.core.utils.series.create_frequency_series( - duration=self.ifo.cropped_duration, sampling_frequency=self.sampling_frequency - ) - return (self.ifo.minimum_frequency <= frequencies) & ( - frequencies <= self.ifo.maximum_frequency - ) - def test_frequency_domain_whitened_strain(self): - white = self.ifo.whitened_frequency_domain_strain[self.frequency_mask] + white = self.ifo.whitened_frequency_domain_strain[self.ifo.cropped_frequency_mask] self._check_frequency_series_whiteness(white) def test_time_domain_whitened_strain(self): @@ -657,8 +648,8 @@ def test_frequency_domain_noise_and_signal_whitening(self): # Whiten the template whitened_signal_ifo = self.ifo.whiten_frequency_series(signal_ifo) white = ( - self.ifo.whitened_frequency_domain_strain[self.frequency_mask] - - whitened_signal_ifo[self.frequency_mask] + self.ifo.whitened_frequency_domain_strain[self.ifo.cropped_frequency_mask] + - whitened_signal_ifo[self.ifo.cropped_frequency_mask] ) self._check_frequency_series_whiteness(white) From bd9e29e6060498b4552fd1fe7ac7706f31681f3b Mon Sep 17 00:00:00 2001 From: Colm Talbot Date: Fri, 23 Jan 2026 11:53:32 -0500 Subject: [PATCH 07/16] TYPO: fix imports --- bilby/gw/detector/strain_data.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/bilby/gw/detector/strain_data.py b/bilby/gw/detector/strain_data.py index b26d409af..aeb6beb7a 100644 --- a/bilby/gw/detector/strain_data.py +++ b/bilby/gw/detector/strain_data.py @@ -1,7 +1,7 @@ import numpy as np from ...core import utils -from ...core.series import CoupledTimeAndFrequencySeries +from ...core.series import create_frequency_series, CoupledTimeAndFrequencySeries from ...core.utils import logger, PropertyAccessor from .. import utils as gwutils @@ -196,7 +196,7 @@ def calculate_frequency_mask(frequency_array): self._times_and_frequencies.frequency_array ) - cropped_frequencies = bilby.core.utils.series.create_frequency_series( + cropped_frequencies = create_frequency_series( duration=self.cropped_duration, sampling_frequency=self.sampling_frequency ) From 7d6a3992a35dfd32b20d31a80f8cfeb8935efc18 Mon Sep 17 00:00:00 2001 From: Colm Talbot Date: Fri, 23 Jan 2026 12:01:49 -0500 Subject: [PATCH 08/16] BUG: fix imports --- bilby/gw/detector/strain_data.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/bilby/gw/detector/strain_data.py b/bilby/gw/detector/strain_data.py index aeb6beb7a..9ce05df0b 100644 --- a/bilby/gw/detector/strain_data.py +++ b/bilby/gw/detector/strain_data.py @@ -1,8 +1,9 @@ import numpy as np from ...core import utils -from ...core.series import create_frequency_series, CoupledTimeAndFrequencySeries +from ...core.series import CoupledTimeAndFrequencySeries from ...core.utils import logger, PropertyAccessor +from ...core.utils.series import create_frequency_series from .. import utils as gwutils From 06adb1ed7cb31142fe58eef0e3200bbfe9cc8393 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Thu, 29 Jan 2026 14:21:25 +0000 Subject: [PATCH 09/16] Update bilby/gw/detector/interferometer.py Co-authored-by: Colm Talbot --- bilby/gw/detector/interferometer.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py index f14e58566..4de9ff0af 100644 --- a/bilby/gw/detector/interferometer.py +++ b/bilby/gw/detector/interferometer.py @@ -616,7 +616,9 @@ def inner_product(self, signal): Returns ======= - float: The optimal signal to noise ratio possible squared + float: + The noise-weighted inner product between the passed signal + and the data stored in the :code:`Interferometer`. """ whitened_signal = self.whiten_frequency_series(signal) whitened_data = self.whitened_frequency_domain_strain From 1e7086ec359959e5f8548d9bc6a65516fed018c3 Mon Sep 17 00:00:00 2001 From: Colm Talbot Date: Thu, 5 Feb 2026 14:37:58 +0000 Subject: [PATCH 10/16] DOC: fix docstring --- bilby/gw/detector/interferometer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py index 4de9ff0af..6b4b304ce 100644 --- a/bilby/gw/detector/interferometer.py +++ b/bilby/gw/detector/interferometer.py @@ -602,7 +602,7 @@ def optimal_snr_squared(self, signal): Returns ======= - float: The optimal signal to noise ratio possible squared + float: The optimal signal-to-noise ratio squared of the signal """ return (abs(self.whiten_frequency_series(signal))**2).sum() From 730a12ed03a6ed789b7e21b22d15cb8879e1264c Mon Sep 17 00:00:00 2001 From: Colm Talbot Date: Fri, 6 Feb 2026 10:18:29 -0500 Subject: [PATCH 11/16] REFACTOR: move whitening functions into gwutils --- bilby/gw/detector/interferometer.py | 68 +++++--------------- bilby/gw/utils.py | 97 +++++++++++++++++++++++++++++ 2 files changed, 111 insertions(+), 54 deletions(-) diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py index 6b4b304ce..915ca486f 100644 --- a/bilby/gw/detector/interferometer.py +++ b/bilby/gw/detector/interferometer.py @@ -657,7 +657,7 @@ def matched_filter_snr(self, signal): """ return self.inner_product(signal) / self.optimal_snr_squared(signal)**0.5 - def whiten_frequency_series(self, frequency_series : np.array) -> np.array: + def whiten_frequency_series(self, frequency_series: np.array) -> np.array: """Whitens a frequency series with the noise properties of the detector .. math:: @@ -677,60 +677,20 @@ def whiten_frequency_series(self, frequency_series : np.array) -> np.array: The frequency series, whitened by the ASD """ if self.crop_duration == 0: - return np.nan_to_num(frequency_series / ( - self.amplitude_spectral_density_array - * np.sqrt(self.duration / 4) - )) * self.frequency_mask + return gwutils.frequency_domain_whiten( + frequency_series=frequency_series, + asd_array=self.amplitude_spectral_density_array, + frequency_mask=self.frequency_mask, + duration=duration, + ) else: - return self.whiten_and_crop(frequency_series=frequency_series) - - def whiten_and_crop(self, frequency_series : np.array) -> np.array: - """ - Whitens a frequency series with the noise properties and applies - the time mask to the whitened time-domain strain. - - First, we naively whiten the data in the frequency domain and apply - our frequency mask :math:`\\tilde{m}(f)` - - .. math:: - \\tilde{a}_w(f) = \\tilde{a}(f) \\tilde{m}(f) / \\sqrt{S_n(f)}. - - We then inverse discrete Fourier transform to get the whitened - time-domain strain :math:`w(t)` and apply the time-domain mask - :math:`m(t)`. We then perform a final discrete Fourier transform. - - .. math:: - \\bar{a}_{w}(f) = \\mathcal{F}(\\mathcal{iF}(\\tilde{a}_w(f))[m(t)]) - - Finally, we normalize the whitened strain by a factor :math:`N` - - .. math:: - N = \\sqrt{\\frac{4}{T_{c}}} - - where :math:`T_{c}` is the cropped duration such that - - .. math:: - Var(n) = \\frac{1}{N} \\sum_{k=0}^N n_W(f_k)n_W^*(f_k) = 2. - - Where the factor of two is due to the independent real and imaginary - components. - - - Parameters - ========== - frequency_series : np.array - The frequency series, whitened by the ASD - """ - initial_white = np.nan_to_num( - frequency_series - * self.frequency_mask - / self.amplitude_spectral_density_array - ) - time_series = infft(initial_white, self.sampling_frequency) - cropped_time_series = time_series[self.time_mask] - cropped_white, _ = nfft(cropped_time_series, self.sampling_frequency) - - return cropped_white * (4 / self.cropped_duration)**0.5 + return gwutils.whiten_and_crop( + frequency_series=frequency_series, + asd_array=asd_array, + frequency_mask=frequency_mask, + time_mask=time_mask, + duration=duration, + ) def get_whitened_time_series_from_whitened_frequency_series( self, diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py index f1f4c0291..9597458be 100644 --- a/bilby/gw/utils.py +++ b/bilby/gw/utils.py @@ -228,6 +228,103 @@ def overlap(signal_a, signal_b, power_spectral_density=None, delta_frequency=Non return sum(integral).real +def frequency_domain_whiten(frequency_series, amplitude_spectral_density, frequency_mask, duration): + """ + Whitens a frequency series with the provided asd in the frequency domain + + .. math:: + \\tilde{a}_w(f) = \\tilde{a}(f) \\sqrt{\\frac{4}{T S_n(f)}} + + Such that + + .. math:: + Var(n) = \\frac{1}{N} \\sum_{k=0}^N n_W(f_k)n_W^*(f_k) = 2 + + Where the factor of two is due to the independent real and imaginary + components. + + Parameters + ========== + frequency_series : np.ndarray + The frequency series to whiten + amplitude_spectral_density: np.ndarray + The amplitude spectral density array to use for whitening + frequency_mask: np.ndarray + A boolean mask to apply to the whitened strain + duration: float + The effective duration of the passed frequency series + + Returns + ======= + np.ndarray + The whitened frequency series + """ + whitened = frequency_series / amplitude_spectral_density + return np.nan_to_num(whitened) * frequency_mask * (self.duration / 4)**0.5 + + +def whiten_and_crop(self, frequency_series, amplitude_spectral_density, frequency_mask, time_mask, duration): + """ + Whitens a frequency series with the noise properties and applies + the time mask to the whitened time-domain strain [1]. + + First, we naively whiten the data in the frequency domain and apply + our frequency mask :math:`\\tilde{m}(f)` + + .. math:: + \\tilde{a}_w(f) = \\tilde{a}(f) \\tilde{m}(f) / \\sqrt{S_n(f)}. + + We then inverse discrete Fourier transform to get the whitened + time-domain strain :math:`w(t)` and apply the time-domain mask + :math:`m(t)`. We then perform a final discrete Fourier transform. + + .. math:: + \\bar{a}_{w}(f) = \\mathcal{F}(\\mathcal{iF}(\\tilde{a}_w(f))[m(t)]) + + Finally, we normalize the whitened strain by a factor :math:`N` + + .. math:: + N = \\sqrt{\\frac{4}{T_{c}}} + + where :math:`T_{c}` is the cropped duration such that + + .. math:: + Var(n) = \\frac{1}{N} \\sum_{k=0}^N n_W(f_k)n_W^*(f_k) = 2. + + Where the factor of two is due to the independent real and imaginary + components. + + + Parameters + ========== + frequency_series : np.ndarray + The frequency series to whiten + amplitude_spectral_density: np.ndarray + The amplitude spectral density array to use for whitening + frequency_mask: np.ndarray + A boolean mask to apply to the whitened strain + time_mask: np.ndarray + A boolean mask used to crop the whitened time series + duration: float + The effective duration of the passed frequency series + + Returns + ======= + np.ndarray + The whitened frequency series + + .. [1] `C. Talbot *et al.* 2025 *Class. Quantum Grav.* **42** 235023 `_ + """ + whitened = frequency_series / amplitude_spectral_density + initial_white = np.nan_to_num(whitened) * frequency_mask + time_series = np.fft.irfft(initial_white) + cropped_time_series = time_series[time_mask] + cropped_white, _ = np.fft.rfft(cropped_time_series) + cropped_duration = duration * len(cropped_time_series) / len(time_series) + + return cropped_white * (4 / cropped_duration)**0.5 + + def zenith_azimuth_to_theta_phi(zenith, azimuth, ifos): """ Convert from the 'detector frame' to the Earth frame. From 032763b60d5bb11e6a6bbae3bfe5988483b9bf5e Mon Sep 17 00:00:00 2001 From: Colm Talbot Date: Fri, 6 Feb 2026 10:22:07 -0500 Subject: [PATCH 12/16] BUG: pre-commit fixes --- bilby/gw/detector/interferometer.py | 13 ++++++------- bilby/gw/utils.py | 2 +- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py index 915ca486f..db24c8c58 100644 --- a/bilby/gw/detector/interferometer.py +++ b/bilby/gw/detector/interferometer.py @@ -9,7 +9,6 @@ from ...core import utils from ...core.utils import docstring, logger, PropertyAccessor, safe_file_dump -from ...core.utils.series import infft, nfft from ...core.utils.env import string_to_boolean from .. import utils as gwutils from .calibration import Recalibrate @@ -679,17 +678,17 @@ def whiten_frequency_series(self, frequency_series: np.array) -> np.array: if self.crop_duration == 0: return gwutils.frequency_domain_whiten( frequency_series=frequency_series, - asd_array=self.amplitude_spectral_density_array, + amplitude_spectral_density=self.amplitude_spectral_density_array, frequency_mask=self.frequency_mask, - duration=duration, + duration=self.duration, ) else: return gwutils.whiten_and_crop( frequency_series=frequency_series, - asd_array=asd_array, - frequency_mask=frequency_mask, - time_mask=time_mask, - duration=duration, + amplitude_spectral_density=self.amplitude_spectral_density_array, + frequency_mask=self.frequency_mask, + time_mask=self.time_mask, + duration=self.duration, ) def get_whitened_time_series_from_whitened_frequency_series( diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py index 9597458be..f927e272d 100644 --- a/bilby/gw/utils.py +++ b/bilby/gw/utils.py @@ -260,7 +260,7 @@ def frequency_domain_whiten(frequency_series, amplitude_spectral_density, freque The whitened frequency series """ whitened = frequency_series / amplitude_spectral_density - return np.nan_to_num(whitened) * frequency_mask * (self.duration / 4)**0.5 + return np.nan_to_num(whitened) * frequency_mask * (duration / 4)**0.5 def whiten_and_crop(self, frequency_series, amplitude_spectral_density, frequency_mask, time_mask, duration): From f3bec335d26a3d986c2e22517675f5f32ca54940 Mon Sep 17 00:00:00 2001 From: Colm Talbot Date: Fri, 6 Feb 2026 11:13:39 -0500 Subject: [PATCH 13/16] BUG: typo fix --- bilby/gw/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py index f927e272d..9b46dcbc7 100644 --- a/bilby/gw/utils.py +++ b/bilby/gw/utils.py @@ -260,7 +260,7 @@ def frequency_domain_whiten(frequency_series, amplitude_spectral_density, freque The whitened frequency series """ whitened = frequency_series / amplitude_spectral_density - return np.nan_to_num(whitened) * frequency_mask * (duration / 4)**0.5 + return np.nan_to_num(whitened) * frequency_mask * (4 / duration)**0.5 def whiten_and_crop(self, frequency_series, amplitude_spectral_density, frequency_mask, time_mask, duration): From a52830d52da0c309c7a97f8e6035bd7ebb978b62 Mon Sep 17 00:00:00 2001 From: Colm Talbot Date: Fri, 6 Feb 2026 13:54:39 -0500 Subject: [PATCH 14/16] TYPO: remove old argument --- bilby/gw/utils.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py index 9b46dcbc7..5b31e8236 100644 --- a/bilby/gw/utils.py +++ b/bilby/gw/utils.py @@ -263,10 +263,10 @@ def frequency_domain_whiten(frequency_series, amplitude_spectral_density, freque return np.nan_to_num(whitened) * frequency_mask * (4 / duration)**0.5 -def whiten_and_crop(self, frequency_series, amplitude_spectral_density, frequency_mask, time_mask, duration): +def whiten_and_crop(frequency_series, amplitude_spectral_density, frequency_mask, time_mask, duration): """ Whitens a frequency series with the noise properties and applies - the time mask to the whitened time-domain strain [1]. + the time mask to the whitened time-domain strain[1]_. First, we naively whiten the data in the frequency domain and apply our frequency mask :math:`\\tilde{m}(f)` @@ -313,7 +313,8 @@ def whiten_and_crop(self, frequency_series, amplitude_spectral_density, frequenc np.ndarray The whitened frequency series - .. [1] `C. Talbot *et al.* 2025 *Class. Quantum Grav.* **42** 235023 `_ + .. [1] C. Talbot et al. 2025 + `Class. Quantum Grav. 42 235023 `_ """ whitened = frequency_series / amplitude_spectral_density initial_white = np.nan_to_num(whitened) * frequency_mask From 85294b4682c3afb23b29af0c1bdfe519aaec3aea Mon Sep 17 00:00:00 2001 From: Colm Talbot Date: Sat, 7 Feb 2026 11:14:20 -0500 Subject: [PATCH 15/16] TYPO: typo fixes in whitening function --- bilby/gw/utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py index 5b31e8236..77906527b 100644 --- a/bilby/gw/utils.py +++ b/bilby/gw/utils.py @@ -313,14 +313,14 @@ def whiten_and_crop(frequency_series, amplitude_spectral_density, frequency_mask np.ndarray The whitened frequency series - .. [1] C. Talbot et al. 2025 + .. [1] C. Talbot et al. 2025 `Class. Quantum Grav. 42 235023 `_ """ whitened = frequency_series / amplitude_spectral_density initial_white = np.nan_to_num(whitened) * frequency_mask time_series = np.fft.irfft(initial_white) cropped_time_series = time_series[time_mask] - cropped_white, _ = np.fft.rfft(cropped_time_series) + cropped_white = np.fft.rfft(cropped_time_series) cropped_duration = duration * len(cropped_time_series) / len(time_series) return cropped_white * (4 / cropped_duration)**0.5 From f885be2ca3b03bcf570ae1ed57332b3d1b0ba305 Mon Sep 17 00:00:00 2001 From: Colm Talbot Date: Fri, 15 May 2026 10:48:02 -0400 Subject: [PATCH 16/16] TEST: add back in removed test --- test/gw/likelihood_test.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/gw/likelihood_test.py b/test/gw/likelihood_test.py index 8d4ba2006..ee16ac4dd 100644 --- a/test/gw/likelihood_test.py +++ b/test/gw/likelihood_test.py @@ -166,6 +166,14 @@ def test_log_likelihood(self): self.assertAlmostEqual(self.likelihood.log_likelihood(self.parameters), -4032.4397343470005, 3) + def test_log_likelihood_ratio(self): + """Test log likelihood ratio returns the correct value""" + self.assertAlmostEqual( + self.likelihood.log_likelihood(self.parameters) - self.likelihood.noise_log_likelihood(), + self.likelihood.log_likelihood_ratio(self.parameters), + 3, + ) + def test_noise_log_likelihood_with_cropping(self): """Test noise log likelihood matches precomputed value""" self.interferometers[0].crop_duration = 1.0