From d91d413f3f719acd6ecd1d2621ad2ace90691ae0 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Wed, 12 Mar 2025 13:07:06 +0000 Subject: [PATCH 1/4] Add new functionality for finding and reading data --- bilby/gw/detector/strain_data.py | 81 ++++++++++++++++++++++++++++++++ 1 file changed, 81 insertions(+) diff --git a/bilby/gw/detector/strain_data.py b/bilby/gw/detector/strain_data.py index bca7acced..d4aed9697 100644 --- a/bilby/gw/detector/strain_data.py +++ b/bilby/gw/detector/strain_data.py @@ -832,3 +832,84 @@ def check_frequency(self, freq): if notch.check_frequency(freq): return True return False + + +def resample_with_gwpy(data, sampling_frequency): + data = data.resample(sampling_frequency) + + +def resample_with_lal(data, sampling_frequency): + import lal + from gwpy.timeseries import TimeSeries + lal_timeseries = data.to_lal() + lal.ResampleREAL8TimeSeries( + lal_timeseries, float(1 / sampling_frequency) + ) + return TimeSeries( + lal_timeseries.data.data, + epoch=lal_timeseries.epoch, + dt=lal_timeseries.deltaT, + ) + + +RESAMPLING_FUNCTIONS = dict( + gwpy=resample_with_gwpy, + lal=resample_with_lal +) + + +def resample_timeseries(data, sampling_frequency, resampling_method="lal"): + + if data.sample_rate.value == sampling_frequency: + logger.info("Sample rate matches data no resampling") + elif resampling_method in RESAMPLING_FUNCTIONS: + logger.info(f"Resampling data to sampling_frequency {sampling_frequency} using {resampling_method}") + return RESAMPLING_FUNCTIONS[resampling_method](data, sampling_frequency) + else: + raise ValueError("Resampling method {resampling_method} not implemented") + + +def find_and_read_data(start, end, ifo, frametype, channel, find_url_kwargs=None, read_kwargs=None, + sampling_frequency=None, resampling_method="lal", dtype="float64"): + """ + Find data using gw_data_find and then read it with gwpy. This assumes the data + exists on the provided host. Usually this is locally. + + Parameters + ---------- + + start, end: float + The GPS start and end time + ifo: str [H1, L1, V1] + The detector to use + frametype: str + The frametype to search for: not including the prepended detector name, e.g. "HOFT_C00_AR" + channel: str + The channel within the frame to read: not include the detector name, e.g. "GDS-CALIB_STRAIN_AR" + find_url_kwargs: dict + A dictionary of kwargs to pass to `gwdatafind.find_urls()` + read_kwargs: dict + A dictionary of kwargs to pass to `gwpy.timeseries.TimeSeries.read()` + + Returns + ------- + data: gwpy.timeseries.TimeSeries + The gwpy timeseries of the data + + """ + from gwpy.timeseries import TimeSeries + from gwdatafind import find_urls + + frametype_with_ifo = f"{ifo}_{frametype}" + channel_with_ifo = f"{ifo}:{channel}" + single_letter_ifo = ifo[0] + + urls = find_urls(single_letter_ifo, frametype_with_ifo, start, end, **find_url_kwargs) + + type_kwargs = dict(dtype=dtype, subok=True, copy=False) + if len(urls) == 0: + raise ValueError(f"No data found for {ifo} {frametype} {start} {end}") + + data = TimeSeries.read(urls, channel_with_ifo, start=start, end=end, **read_kwargs) .astype(**type_kwargs) + data = resample_timeseries(data, sampling_frequency, resampling_method) + return data From 3f00ea7deeb49ec0faec201b3bb858834a2efe0e Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Wed, 12 Mar 2025 13:12:34 +0000 Subject: [PATCH 2/4] Update the docstrings --- bilby/gw/detector/strain_data.py | 57 +++++++++++++++++++++++++++++++- 1 file changed, 56 insertions(+), 1 deletion(-) diff --git a/bilby/gw/detector/strain_data.py b/bilby/gw/detector/strain_data.py index d4aed9697..5db38703f 100644 --- a/bilby/gw/detector/strain_data.py +++ b/bilby/gw/detector/strain_data.py @@ -835,10 +835,42 @@ def check_frequency(self, freq): def resample_with_gwpy(data, sampling_frequency): + """ + Resample a GWPy TimeSeries to a new sampling frequency. + + Parameters: + ---------- + data : gwpy.timeseries.TimeSeries + The input time series data to be resampled. + sampling_frequency : float + The target sampling frequency (Hz) for resampling. + + Returns: + ------- + gwpy.timeseries.TimeSeries + A new TimeSeries object resampled to the desired frequency. + + """ data = data.resample(sampling_frequency) def resample_with_lal(data, sampling_frequency): + """ + Resample a GWPy TimeSeries using LAL's ResampleREAL8TimeSeries function. + + Parameters: + ---------- + data : gwpy.timeseries.TimeSeries + The input time series data to be resampled. + sampling_frequency : float + The target sampling frequency (Hz) for resampling. + + Returns: + ------- + gwpy.timeseries.TimeSeries + A new TimeSeries object resampled to the desired frequency. + + """ import lal from gwpy.timeseries import TimeSeries lal_timeseries = data.to_lal() @@ -859,7 +891,31 @@ def resample_with_lal(data, sampling_frequency): def resample_timeseries(data, sampling_frequency, resampling_method="lal"): + """ + Resample a time series to a specified sampling frequency using a chosen method. + + Parameters: + ---------- + data : gwpy.timeseries.TimeSeries + The input time series data to be resampled. + sampling_frequency : float + The target sampling frequency (Hz) for resampling. + resampling_method : str, optional + The resampling method to use. Defaults to "lal". + Must be one of the methods defined in + `bilby.gw.detector.strain_data.RESAMPLING_FUNCTIONS`. + + Returns: + ------- + gwpy.timeseries.TimeSeries + A new TimeSeries object resampled to the desired frequency. + + Raises: + ------ + ValueError + If the specified resampling method is not implemented. + """ if data.sample_rate.value == sampling_frequency: logger.info("Sample rate matches data no resampling") elif resampling_method in RESAMPLING_FUNCTIONS: @@ -877,7 +933,6 @@ def find_and_read_data(start, end, ifo, frametype, channel, find_url_kwargs=None Parameters ---------- - start, end: float The GPS start and end time ifo: str [H1, L1, V1] From 03d6f4d7470acc0bbf48e54cdd4172ea5ecd5894 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Wed, 26 Aug 2026 17:06:53 +0100 Subject: [PATCH 3/4] Address review comments AI summary: 1. resample_with_gwpy now accepts **kwargs passed through to .resample(), and fixed the bug where it discarded the result instead of returning it (it was returning None). 2. resample_with_lal accepts **kwargs for interface consistency, raising a clear ValueError if any are passed (LAL's resampler doesn't support extra options). 3. resample_timeseries passes **kwargs through to whichever resampling function is selected, and now returns data unchanged when no resampling is needed (previously returned None in that branch too). 4. find_and_read_data gained a resample_kwargs parameter so resampling kwargs can flow all the way through; fixed the docstring to clarify ifo isn't limited to H1/L1/V1 (K1, A1 work fine via ifo[0]); fixed the "not include" grammar; and fixed a latent bug where find_url_kwargs/read_kwargs defaulting to None would crash on **None. 5. Added a TestResampling and TestFindAndReadData test class covering the resampling functions and find_and_read_data (mocking gwdatafind.find_urls and TimeSeries.read), addressing his request for test coverage. --- bilby/gw/detector/strain_data.py | 43 +++++++--- test/gw/detector/strain_data_test.py | 122 +++++++++++++++++++++++++++ 2 files changed, 152 insertions(+), 13 deletions(-) diff --git a/bilby/gw/detector/strain_data.py b/bilby/gw/detector/strain_data.py index 140df848e..cd4fdcd34 100644 --- a/bilby/gw/detector/strain_data.py +++ b/bilby/gw/detector/strain_data.py @@ -871,7 +871,7 @@ def check_frequency(self, freq): return False -def resample_with_gwpy(data, sampling_frequency): +def resample_with_gwpy(data, sampling_frequency, **kwargs): """ Resample a GWPy TimeSeries to a new sampling frequency. @@ -881,6 +881,8 @@ def resample_with_gwpy(data, sampling_frequency): The input time series data to be resampled. sampling_frequency : float The target sampling frequency (Hz) for resampling. + **kwargs: + Additional keyword arguments passed to `gwpy.timeseries.TimeSeries.resample`. Returns: ------- @@ -888,10 +890,10 @@ def resample_with_gwpy(data, sampling_frequency): A new TimeSeries object resampled to the desired frequency. """ - data = data.resample(sampling_frequency) + return data.resample(sampling_frequency, **kwargs) -def resample_with_lal(data, sampling_frequency): +def resample_with_lal(data, sampling_frequency, **kwargs): """ Resample a GWPy TimeSeries using LAL's ResampleREAL8TimeSeries function. @@ -901,6 +903,8 @@ def resample_with_lal(data, sampling_frequency): The input time series data to be resampled. sampling_frequency : float The target sampling frequency (Hz) for resampling. + **kwargs: + Not used, present for interface compatibility with `resample_with_gwpy`. Returns: ------- @@ -908,6 +912,8 @@ def resample_with_lal(data, sampling_frequency): A new TimeSeries object resampled to the desired frequency. """ + if kwargs: + raise ValueError(f"resample_with_lal does not support additional kwargs: {kwargs}") import lal from gwpy.timeseries import TimeSeries lal_timeseries = data.to_lal() @@ -927,7 +933,7 @@ def resample_with_lal(data, sampling_frequency): ) -def resample_timeseries(data, sampling_frequency, resampling_method="lal"): +def resample_timeseries(data, sampling_frequency, resampling_method="lal", **kwargs): """ Resample a time series to a specified sampling frequency using a chosen method. @@ -941,6 +947,8 @@ def resample_timeseries(data, sampling_frequency, resampling_method="lal"): The resampling method to use. Defaults to "lal". Must be one of the methods defined in `bilby.gw.detector.strain_data.RESAMPLING_FUNCTIONS`. + **kwargs: + Additional keyword arguments passed to the chosen resampling function. Returns: ------- @@ -955,15 +963,17 @@ def resample_timeseries(data, sampling_frequency, resampling_method="lal"): """ if data.sample_rate.value == sampling_frequency: logger.info("Sample rate matches data no resampling") + return data.copy() elif resampling_method in RESAMPLING_FUNCTIONS: logger.info(f"Resampling data to sampling_frequency {sampling_frequency} using {resampling_method}") - return RESAMPLING_FUNCTIONS[resampling_method](data, sampling_frequency) + return RESAMPLING_FUNCTIONS[resampling_method](data, sampling_frequency, **kwargs) else: - raise ValueError("Resampling method {resampling_method} not implemented") + raise ValueError(f"Resampling method {resampling_method} not implemented") def find_and_read_data(start, end, ifo, frametype, channel, find_url_kwargs=None, read_kwargs=None, - sampling_frequency=None, resampling_method="lal", dtype="float64"): + sampling_frequency=None, resampling_method="lal", resample_kwargs=None, + dtype="float64"): """ Find data using gw_data_find and then read it with gwpy. This assumes the data exists on the provided host. Usually this is locally. @@ -972,16 +982,19 @@ def find_and_read_data(start, end, ifo, frametype, channel, find_url_kwargs=None ---------- start, end: float The GPS start and end time - ifo: str [H1, L1, V1] - The detector to use + ifo: str + The detector to use, e.g. "H1", "L1", "V1", "K1", "A1" frametype: str The frametype to search for: not including the prepended detector name, e.g. "HOFT_C00_AR" channel: str - The channel within the frame to read: not include the detector name, e.g. "GDS-CALIB_STRAIN_AR" + The channel within the frame to read: not including the detector name, e.g. "GDS-CALIB_STRAIN_AR" find_url_kwargs: dict A dictionary of kwargs to pass to `gwdatafind.find_urls()` read_kwargs: dict A dictionary of kwargs to pass to `gwpy.timeseries.TimeSeries.read()` + resample_kwargs: dict + A dictionary of kwargs to pass to the resampling function, see + `bilby.gw.detector.strain_data.resample_timeseries`. Returns ------- @@ -992,16 +1005,20 @@ def find_and_read_data(start, end, ifo, frametype, channel, find_url_kwargs=None from gwpy.timeseries import TimeSeries from gwdatafind import find_urls + find_url_kwargs = find_url_kwargs or dict() + read_kwargs = read_kwargs or dict() + resample_kwargs = resample_kwargs or dict() + frametype_with_ifo = f"{ifo}_{frametype}" channel_with_ifo = f"{ifo}:{channel}" single_letter_ifo = ifo[0] urls = find_urls(single_letter_ifo, frametype_with_ifo, start, end, **find_url_kwargs) - type_kwargs = dict(dtype=dtype, subok=True, copy=False) if len(urls) == 0: raise ValueError(f"No data found for {ifo} {frametype} {start} {end}") - data = TimeSeries.read(urls, channel_with_ifo, start=start, end=end, **read_kwargs) .astype(**type_kwargs) - data = resample_timeseries(data, sampling_frequency, resampling_method) + type_kwargs = dict(dtype=dtype, subok=True, copy=False) + data = TimeSeries.read(urls, channel_with_ifo, start=start, end=end, **read_kwargs).astype(**type_kwargs) + data = resample_timeseries(data, sampling_frequency, resampling_method, **resample_kwargs) return data diff --git a/test/gw/detector/strain_data_test.py b/test/gw/detector/strain_data_test.py index 0f82a40a2..f60e05608 100644 --- a/test/gw/detector/strain_data_test.py +++ b/test/gw/detector/strain_data_test.py @@ -3,8 +3,15 @@ import numpy as np import scipy.signal +from gwpy.timeseries import TimeSeries import bilby +from bilby.gw.detector.strain_data import ( + find_and_read_data, + resample_timeseries, + resample_with_gwpy, + resample_with_lal, +) class TestInterferometerStrainData(unittest.TestCase): @@ -420,5 +427,120 @@ def test_init_fail(self): bilby.gw.detector.strain_data.NotchList([(30, 20, 20)]) +class TestResampling(unittest.TestCase): + def setUp(self): + self.sampling_frequency = 512 + self.data = TimeSeries( + np.random.normal(0, 1, self.sampling_frequency * 4), + sample_rate=self.sampling_frequency, + epoch=0, + ) + + def test_resample_with_gwpy(self): + new_data = resample_with_gwpy(self.data, self.sampling_frequency / 2) + self.assertIsInstance(new_data, TimeSeries) + self.assertEqual(new_data.sample_rate.value, self.sampling_frequency / 2) + + def test_resample_with_lal(self): + new_data = resample_with_lal(self.data, self.sampling_frequency / 2) + self.assertIsInstance(new_data, TimeSeries) + self.assertEqual(new_data.sample_rate.value, self.sampling_frequency / 2) + + def test_resample_with_lal_unsupported_kwargs(self): + with self.assertRaises(ValueError): + resample_with_lal(self.data, self.sampling_frequency / 2, window="hann") + + def test_resample_timeseries_same_rate(self): + new_data = resample_timeseries(self.data, self.sampling_frequency) + self.assertIsNot(new_data, self.data) + self.assertTrue(np.array_equal(new_data.value, self.data.value)) + + def test_resample_timeseries_lal(self): + new_data = resample_timeseries( + self.data, self.sampling_frequency / 2, resampling_method="lal" + ) + self.assertEqual(new_data.sample_rate.value, self.sampling_frequency / 2) + + def test_resample_timeseries_gwpy(self): + new_data = resample_timeseries( + self.data, self.sampling_frequency / 2, resampling_method="gwpy" + ) + self.assertEqual(new_data.sample_rate.value, self.sampling_frequency / 2) + + def test_resample_timeseries_kwargs_passed_through(self): + new_data = resample_timeseries( + self.data, + self.sampling_frequency / 2, + resampling_method="gwpy", + window="hann", + ) + self.assertEqual(new_data.sample_rate.value, self.sampling_frequency / 2) + + def test_resample_timeseries_invalid_method(self): + with self.assertRaises(ValueError): + resample_timeseries( + self.data, self.sampling_frequency / 2, resampling_method="not-a-method" + ) + + +class TestFindAndReadData(unittest.TestCase): + def setUp(self): + self.start = 0 + self.end = 4 + self.ifo = "H1" + self.frametype = "HOFT_C00_AR" + self.channel = "GDS-CALIB_STRAIN_AR" + self.sampling_frequency = 512 + self.data = TimeSeries( + np.random.normal(0, 1, self.sampling_frequency * 4), + sample_rate=self.sampling_frequency, + epoch=self.start, + ) + + @mock.patch("gwdatafind.find_urls") + @mock.patch("gwpy.timeseries.TimeSeries.read") + def test_find_and_read_data(self, mock_read, mock_find_urls): + mock_find_urls.return_value = ["file://fake/H-H1_HOFT_C00_AR-0-4.gwf"] + mock_read.return_value = self.data + + data = find_and_read_data( + self.start, self.end, self.ifo, self.frametype, self.channel, + sampling_frequency=self.sampling_frequency, + ) + + mock_find_urls.assert_called_once_with( + "H", f"{self.ifo}_{self.frametype}", self.start, self.end + ) + mock_read.assert_called_once_with( + mock_find_urls.return_value, + f"{self.ifo}:{self.channel}", + start=self.start, + end=self.end, + ) + self.assertEqual(data.sample_rate.value, self.sampling_frequency) + + @mock.patch("gwdatafind.find_urls") + def test_find_and_read_data_no_urls_raises(self, mock_find_urls): + mock_find_urls.return_value = [] + with self.assertRaises(ValueError): + find_and_read_data( + self.start, self.end, self.ifo, self.frametype, self.channel, + sampling_frequency=self.sampling_frequency, + ) + + @mock.patch("gwdatafind.find_urls") + @mock.patch("gwpy.timeseries.TimeSeries.read") + def test_find_and_read_data_resamples(self, mock_read, mock_find_urls): + mock_find_urls.return_value = ["file://fake/H-H1_HOFT_C00_AR-0-4.gwf"] + mock_read.return_value = self.data + + data = find_and_read_data( + self.start, self.end, self.ifo, self.frametype, self.channel, + sampling_frequency=self.sampling_frequency / 2, + ) + + self.assertEqual(data.sample_rate.value, self.sampling_frequency / 2) + + if __name__ == "__main__": unittest.main() From 0584cd9268455ccfd1c112dd69ca741b9370cd4d Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Wed, 26 Aug 2026 21:44:38 +0100 Subject: [PATCH 4/4] Address review comments from copilot --- bilby/gw/detector/strain_data.py | 20 +++++++++++++++----- test/gw/detector/strain_data_test.py | 12 ++++++++++++ 2 files changed, 27 insertions(+), 5 deletions(-) diff --git a/bilby/gw/detector/strain_data.py b/bilby/gw/detector/strain_data.py index cd4fdcd34..d4b60f7f8 100644 --- a/bilby/gw/detector/strain_data.py +++ b/bilby/gw/detector/strain_data.py @@ -914,8 +914,11 @@ def resample_with_lal(data, sampling_frequency, **kwargs): """ if kwargs: raise ValueError(f"resample_with_lal does not support additional kwargs: {kwargs}") - import lal - from gwpy.timeseries import TimeSeries + try: + import lal + from gwpy.timeseries import TimeSeries + except ModuleNotFoundError: + raise ModuleNotFoundError("Cannot resample with lal: lal and gwpy are required") lal_timeseries = data.to_lal() lal.ResampleREAL8TimeSeries( lal_timeseries, float(1 / sampling_frequency) @@ -992,6 +995,9 @@ def find_and_read_data(start, end, ifo, frametype, channel, find_url_kwargs=None A dictionary of kwargs to pass to `gwdatafind.find_urls()` read_kwargs: dict A dictionary of kwargs to pass to `gwpy.timeseries.TimeSeries.read()` + sampling_frequency: float, optional + The target sampling frequency (Hz) to resample the data to. If not + given (default), no resampling is performed. resample_kwargs: dict A dictionary of kwargs to pass to the resampling function, see `bilby.gw.detector.strain_data.resample_timeseries`. @@ -1002,8 +1008,11 @@ def find_and_read_data(start, end, ifo, frametype, channel, find_url_kwargs=None The gwpy timeseries of the data """ - from gwpy.timeseries import TimeSeries - from gwdatafind import find_urls + try: + from gwpy.timeseries import TimeSeries + from gwdatafind import find_urls + except ModuleNotFoundError: + raise ModuleNotFoundError("Cannot find and read data: gwpy and gwdatafind are required") find_url_kwargs = find_url_kwargs or dict() read_kwargs = read_kwargs or dict() @@ -1020,5 +1029,6 @@ def find_and_read_data(start, end, ifo, frametype, channel, find_url_kwargs=None type_kwargs = dict(dtype=dtype, subok=True, copy=False) data = TimeSeries.read(urls, channel_with_ifo, start=start, end=end, **read_kwargs).astype(**type_kwargs) - data = resample_timeseries(data, sampling_frequency, resampling_method, **resample_kwargs) + if sampling_frequency is not None: + data = resample_timeseries(data, sampling_frequency, resampling_method, **resample_kwargs) return data diff --git a/test/gw/detector/strain_data_test.py b/test/gw/detector/strain_data_test.py index f60e05608..7ae6eaac8 100644 --- a/test/gw/detector/strain_data_test.py +++ b/test/gw/detector/strain_data_test.py @@ -541,6 +541,18 @@ def test_find_and_read_data_resamples(self, mock_read, mock_find_urls): self.assertEqual(data.sample_rate.value, self.sampling_frequency / 2) + @mock.patch("gwdatafind.find_urls") + @mock.patch("gwpy.timeseries.TimeSeries.read") + def test_find_and_read_data_no_sampling_frequency_skips_resampling(self, mock_read, mock_find_urls): + mock_find_urls.return_value = ["file://fake/H-H1_HOFT_C00_AR-0-4.gwf"] + mock_read.return_value = self.data + + data = find_and_read_data( + self.start, self.end, self.ifo, self.frametype, self.channel, + ) + + self.assertEqual(data.sample_rate.value, self.sampling_frequency) + if __name__ == "__main__": unittest.main()