gwastro
Extreme discontinuities in optimized_match that are not present in match
When calculating the mismatch between waves, discontinuities in the match value can occur when using pycbc.filter.matchedfilter.optimized_match that are not present when using match. Specifically, the research I'm doing is investigating lensed vs unlensed waves and certain combinations of chirp mass and time delay cause a spike in the mismatch (1-match). An example is shown in the following graph:
This is not an isolated issue, it occurs in multiple regions of phase-space, as seen by the sharp red-orange regions in the following contour plot:
Or zoomed-in:
Here's a link to a pkl file containing the waveform strain for the range of waves from the first plot for testing purposes. Usage is as follows:
with open('data.pkl', 'rb') as file:
td_list, lensed_waveform_list, unlensed_waveform_list = pickle.load(file)
Each item in lensed_waveform_list and unlensed_waveform_list is a FrequencySeries object of the strain of the wave. The PSD is calculated using the following function:
def Sn(f_arr, f_min=20, delta_f=0.25, frequencySeries=True):
"""
Calculates the power spectral density of the aLIGO noise curve based on arXiv:0903.0338.
Parameters
----------
f_arr : np.ndarray
The frequency array.
f_min : float, optional
The minimum frequency. Defaults to 20 Hz.
delta_f : float, optional
The frequency step size. Defaults to 0.25 Hz.
frequencySeries : bool, optional
If True, returns a FrequencySeries object. Defaults to True.
Returns
-------
np.ndarray or FrequencySeries
The power spectral density of the aLIGO noise curve.
"""
Sn_val = np.zeros_like(f_arr)
for i in range(len(f_arr)):
if f_arr[i] < f_min:
Sn_val[i] = np.inf
else:
S0 = 1e-49
f0 = 215
Sn_temp = (
np.power(f_arr[i] / f0, -4.14)
- 5 * np.power(f_arr[i] / f0, -2)
+ 111
* (
(1 - np.power(f_arr[i] / f0, 2) + 0.5 * np.power(f_arr[i] / f0, 4))
/ (1 + 0.5 * np.power(f_arr[i] / f0, 2))
)
)
Sn_val[i] = Sn_temp * S0
if frequencySeries:
return FrequencySeries(Sn_val, delta_f=delta_f)
return Sn_val
@bmck039 Thank you for reporting this. It definitely seems odd. That function uses an optimizer to try to fit to subsample precision. Perhaps that process is failing?
https://github.com/gwastro/pycbc/blob/master/pycbc/filter/matchedfilter.py#L2075
It doesn't look like we are checking the result of the minimizer to see if it has properly converged. I don't know if @jacopok would be able to look at this, but if you are comfortable modifying the code, perhaps to check the convergence result. One may need to pass in the ability to set the minimizer options to at least give some way to avoid failure modes (assuming this is the problem to begin with).
We are happy to accept contributions here if you have ideas to address this.
Yeah, it's the optimizer messing up - but there's an easy solution.
This is what happened in your test case - one might wonder, how is it possible that the optimizer found a local minimum outside its bounds? It turns out that the brent optimization method does not enforce its bounds, even if they are explicitly given.
So, the solution is simply to switch to the bounded method, don't know why I didn't do it from the beginning. Then your plot will look like this:
Would you be willing to make the PR @bmck039 ?
Essentially you'd need to change the optimizer call to
res = minimize_scalar(
to_minimize,
method="bounded",
bounds=(-delta_t, delta_t)
)
in the code @ahnitz linked to.
It would be nice to have an additional test in the automated test suite checking that optimized_match always returns a higher value than match (which will be true as long as the bounds are respected) - maybe hardcoding the waveforms at index 545 in your example, the ones for which the bug happens?
