FEAT: add option to crop whitened time series#1030
Conversation
f39871f to
5fbc454
Compare
| @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 | ||
| ) | ||
|
|
There was a problem hiding this comment.
We should probably have a utility function that does this.
0b74e2b to
0972d62
Compare
GregoryAshton
left a comment
There was a problem hiding this comment.
Overall looks good and I'm supportive, just a few things I wanted to check.
b8f4b97 to
c7e762a
Compare
| else: | ||
| return self.whiten_and_crop(frequency_series=frequency_series) | ||
|
|
||
| def whiten_and_crop(self, frequency_series : np.array) -> np.array: |
There was a problem hiding this comment.
It seems like this procedure is more fundamental than suggested by its being a method of the interferometer class. The nice thing about using the gwutils SNR functions was that you could calculate SNRs using them without necessarily initializing an interferometer object. Is it possible to pull this out to the utils and just call it here?
| return np.nan_to_num(whitened) * frequency_mask * (4 / duration)**0.5 | ||
|
|
||
|
|
||
| def whiten_and_crop(frequency_series, amplitude_spectral_density, frequency_mask, time_mask, duration): |
There was a problem hiding this comment.
Should the name state what input it accepts? E.g. like 'frequency_domain_whiten' to make it clear it is acting on the frequency-domain input?
There was a problem hiding this comment.
I can make this change, but it's pretty clear from the signature and docstring.
GregoryAshton
left a comment
There was a problem hiding this comment.
A few more comments after another look.
| * 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 |
There was a problem hiding this comment.
Why does this use both np.sqrt and **0.5. It would probably be easier to read if it was consistent.
There was a problem hiding this comment.
I think I was trying to avoid explicitly using numpy, but that's unavoidable since we need an fft and sqrt is sometimes meaningfully faster.
| * self.time_mask.mean()**0.5 | |
| * np.sqrt(self.time_mask.mean()) |
| 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: |
There was a problem hiding this comment.
What happens (or can) crop_duration be "None" or None etc? This would silently fail and then trigger the else clause, likely also erroring, but could be confusing?
There was a problem hiding this comment.
This would fail in the crop_duration.setter in the InterferometerStrainData with a reasonably explicit error message.
| ) | ||
| _mask = interferometer.frequency_mask | ||
|
|
||
| if 'recalib_index' in parameters: |
There was a problem hiding this comment.
Why does the signal no longer require masking here?
There was a problem hiding this comment.
Unclear, is this tested?
Co-authored-by: Colm Talbot <talbotcolm@gmail.com>
962bc37 to
85294b4
Compare
Refactor likelihood tests to store results in variables before assertions.
asb5468
left a comment
There was a problem hiding this comment.
A few more general questions/comments:
- In https://arxiv.org/pdf/2508.11091, there are two options presented for how to treat the template, "Signal models can then either be generated in the frequency domain with frequency spacing 1/T and whitened before subtracting from the whitened and cropped frequency-domain data or generated with the full duration T + 2∆ and processed using the same procedure as the data." Do I understand correctly that this MR has implemented the first option?
- The code itself seems fine, but it's hard to tell definitively whether there are any unintended consequences, probably because I'm not familiar with all the review tests, and right now they are failing. Have you done any more extensive test runs that you could share, like a GW150914 analysis with the new likelihood with no cropping to make sure the results are unchanged?
| 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() |
There was a problem hiding this comment.
Why does the whitened signal need to be transposed? In the function below, template_template_inner_product, it is not.
There was a problem hiding this comment.
I suspect this is a typo, it likely wouldn't do anything for one dimensional arrays which I think we have an implicit expectation this is, but I'll remove it.
| * 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 |
There was a problem hiding this comment.
Where does this extra factor of self.time_mask.mean()**0.5 come from? Can you point to a specific equation in https://arxiv.org/pdf/2508.11091?
There was a problem hiding this comment.
Note for myself later, it is (10), this is
| 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 |
There was a problem hiding this comment.
The docstring here should add the units for this argument. I'm assuming it's seconds?
| @property | ||
| def crop_duration(self): | ||
| """ | ||
| The duration of data to crop at the beginning/end of the segment |
There was a problem hiding this comment.
Add the units of this property to the docstring. Also, crop_duration and cropped_duration are really similar and I'm leery of them getting mixed up. Perhaps crop_duration could be renamed to crop_time?
| cropped_frequencies | ||
| ) | ||
|
|
||
| self._frequency_mask_updated = True |
There was a problem hiding this comment.
It seems like self._frequency_mask_updated now controls whether both frequency_mask and cropped_frequency_mask are updated. Is there ever a situation where a user would want to update one but not the other?
There was a problem hiding this comment.
I don't think so, I feel like we should require that they are kept synchronized.
| signal[_mask] *= self.calibration_draws[interferometer.name][int(parameters['recalib_index'])] | ||
| signal *= self.calibration_draws[interferometer.name][int(parameters['recalib_index'])] |
There was a problem hiding this comment.
Why is this no longer being applied only to the masked frequencies?
| 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[:-1], axis=0).T |
There was a problem hiding this comment.
Why is the 4/T normalization no longer applied here?
There was a problem hiding this comment.
Ah, is it because it's already applied within interferometer.whitened_frequency_domain_strain and whitened_signal?
| with self.assertRaises(ValueError): | ||
| self.ifo.inject_signal(injection_polarizations=None, parameters=None) | ||
|
|
||
| def test_optimal_snr_squared(self): |
This PR implements the proposed modification to the likelihood from https://iopscience.iop.org/article/10.1088/1361-6382/ae1ac7/meta.
The main new feature is adding a
crop_durationto theInterferometerStrainDataand some accompanying utility code.I modified how the
GravitationalWaveTransientandBasicGravitationalWaveTransientcall the data. I did not have to change the reference likelihood tests, so there is good evidence that this doesn't impact the numerical values when not using the cropping.I suspect there is a performance penalty due to repeatedly whitening the data. There should be a way to cache this calculation, it just requires some cleverness about when the time/frequency masks might have changed.
I also haven't validated how this interacts with time/calibration marginalization.