Calculating delta_frequency in source.py can cause inconsistent arrays

This was recently discovered. For certain sets of duration and sampling_frequency there are rounding errors that cause inconsistent shapes between frequency_array and the output modes.

Try running duration=32, sampling_frequency=4400 with an injected lalsim binary neutron star to reproduce this issue.

One solution might be replacing delta_frequency = frequency_array[1] - frequency_array[0] with delta_frequency = np.mean(np.diff(frequency_array))

Hi @Moritz Hübner, thank you for your answers!

Indeed I agree with your comment on the issue: #153 (closed)

In my opinion the first culprit is in the core.utils.create_frequency_series() function where sampling_frequency / 2. is appended at the end of the array. The result is that the last two frequency values of frequency_array are not spaced by delta_freq.

Then in the source.lal_binary_neutron_star() model, we set maximum_frequency = frequency_array[-1], ie maximum_frequency = sampling_frequency / 2.

But from what I understand lalsim computes its number of frequencies (at least for the TF2 waveform) using this formula:

n = (size_t) (f_max / deltaF + 1)

(cf this page line 512: https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/_l_a_l_sim_inspiral_taylor_f2_8c_source.html).

Replacing f_max in lalsim by bilby’s input maximum_frequency it yields:

n = (size_t) (sampling_frequency / 2 / deltaF + 1)

Even if that was to be corrected I think there is another issue: the frequency_array created with np.linspace() is inherently not evenly spaced, to some numerical precision at least, and that ends up in different array length calculation depending on what value we use for delta_freq. This little following code shows it:

import numpy as np
duration = 32
sampling_frequency = 4397
f_max = sampling_frequency / 2
delta_freq = 1 / duration
number_of_frequencies = int(f_max / delta_freq + 1) # = 70353

frequencies = delta_freq * np.linspace(0, number_of_frequencies, number_of_frequencies)

print(delta_freq == frequencies[1]-frequencies[0]) # returns False
print(delta_freq == frequencies[-1]-frequencies[-2]) # returns False
print(frequencies[1]-frequencies[0] == frequencies[-1]-frequencies[-2]) # returns False

# Now if we want to recover the number of frequencies as lalsim does we could get different results:
n1 = int(frequencies[-1]/delta_freq + 1)
n2 = int(frequencies[-1]/(frequencies[1]-frequencies[0]) + 1)
n3 = int(frequencies[-1]/(frequencies[-1]-frequencies[-2]) + 1)
print(n1, n2, n3) # returns 70354 70353 70352

However if we use np.arange() instead of np.linspace(), those discrepancies are corrected:

import numpy as np
duration = 32
sampling_frequency = 4397
f_max = sampling_frequency / 2
delta_freq = 1 / duration
number_of_frequencies = int(f_max / delta_freq + 1) # = 70353

frequencies = np.arange(0, delta_freq*number_of_frequencies, delta_freq)

print(delta_freq == frequencies[1]-frequencies[0]) # returns True
print(delta_freq == frequencies[-1]-frequencies[-2]) # returns True
print(frequencies[1]-frequencies[0] == frequencies[-1]-frequencies[-2]) # returns True

# Now if we want to recover the number of frequencies as lalsim does we could get different results:
n1 = int(frequencies[-1]/delta_freq + 1)
n2 = int(frequencies[-1]/(frequencies[1]-frequencies[0]) + 1)
n3 = int(frequencies[-1]/(frequencies[-1]-frequencies[-2]) + 1)
print(n1, n2, n3) # returns 70353 70353 70353

So I would suggest to correct np.linspace() in core.utils.create_frequency_series() with np.arange().

To be exhaustive (even though I feel I already got everyone bored here…) I should mention that if in my little example you change the values using: duration = 31.91560167202568, sampling_frequency = 4397.2097653746623, some discrepancies reappear even with np.arange() but it keeps delta_freq == frequencies[1]-frequencies[0] so then bilby and lalsim agree on array lengths.

I hope that helps a bit and that I’m not fooling myself !

Nice digging @marc.arene, at an absolute minimum we want to agree with lalsim on the lengths of the array's.

@marc.arene We will probably come up with a lasting solution early next year. I have been working on a MR but I will need to sit down with some people here before it is ready to go. One thing that we realized is that duration and sampling_frequency will have to multiply to an integer number for everything to be consistent. Otherwise it will not be possible create a time array that is both equidistantly spaced and will return the same duration and sampling_frequency when we put it into bilby.core.utils.get_sampling_frequency_and_duration_from_time_array(time_array)

This also relates to to FFTs/IFFTs, so getting it all correct is a bit tricky.

!298 (merged)

closed