gwastro
Numpy 2 leads to numerical underflows in PyCBC Live
Numpy 2 introduced changes to data type promotion, in particular this one:
np.float32(3) + 3. now returns a float32 when it previously returned a float64.
This introduces a subtle bug in PyCBC Live. PSD estimates are internally represented as single-precision float arrays, which is possible because the GW data is scaled by the "dynamic range factor" (see FindChirp paper for details) before doing math on it. In a couple places, however, we convert the PSD estimate from the internal representation back to "physical" values by doing this:
actual_psd = single_precision_psd / DYN_RANGE_FAC ** 2
where single_precision_psd is typically around 1e-4, and DYN_RANGE_FAC is ~6e20. This line used to work before Numpy 2 because single_precision_psd would be promoted to double precision before being divided. In Numpy 2, the promotion no longer happens, so the division underflows and we get an array of (single-precision) zeroes.
In practice this causes the checks in the example script to fail (well done, checks!) and it causes the PSD estimates in the LIGOLW XML files to be all zeroes (which would break whatever reads them later, like BAYESTAR).
I think this went unnoticed because Numpy 2 was never pulled in any of the CI tests due to other dependencies. Thanks to Thomas for spotting it while trying the PyCBC Live example on his laptop :)
@titodalcanton As a first step should we force numpy 2.0 to be used for tests? Although as you suggest, this probably is incompatible with some other libraries.
Apart from fixing this right away (Thomas is working on a PR) I am inclined to just let the tests use whatever they have to, and rather add even more strict checks to the example. As soon as the dependencies change and we start getting Numpy 2, we will notice immediately.