The bug on the backend of `ifft' function from pycbc/fft/fftw.py

[BenjaminDbb]

Recently, we use the same pycbc code to perform ifft with the GW waveform in the different computers, then we obtain different results.

As we check the details in the inner code of pycbc, we find that, when pycbc calls the ifft from the npfft.py, the result is correct, while calling the ifft from fftw.py, the result is wrong. (In default, the macOS system calls the ifft in fftw.py, some Linux systems call the ifft in the npfft.py)

ifft from fftw.py:

def ifft(invec, outvec, prec, itype, otype):
    theplan, destroy = plan(len(outvec), invec.dtype, outvec.dtype, FFTW_BACKWARD,
                            get_measure_level(),(check_aligned(invec.data) and check_aligned(outvec.data)),
                   _scheme.mgr.state.num_threads, (invec.ptr == outvec.ptr))
    execute(theplan, invec, outvec)
    destroy(theplan)

The bug is, the ifft function from fftw.py changes the values of invec. It makes no sense.

The minimal code is

import pylab
from pycbc import types, fft, waveform

delta_t = 1 / 4096
delta_f = 1.0 / 4
f_lower = 40

hp_tilde, _ = waveform.get_fd_waveform(approximant='IMRPhenomA', mass1=6, mass2=6, delta_f=delta_f, f_lower=f_lower)

tlen = int(1.0 / delta_t / delta_f)
hp_tilde.resize(tlen / 2 + 1)

# perform the first ifft
hp = types.TimeSeries(types.zeros(tlen), delta_t=delta_t)
fft.ifft(hp_tilde, hp)

# perform the second ifft (same operation as the first one)
hp_new = types.TimeSeries(types.zeros(tlen), delta_t=delta_t)
fft.ifft(hp_tilde, hp_new)

pylab.plot(hp_new.sample_times, hp_new, label='Second ifft', color='red')  
pylab.plot(hp.sample_times, hp, label='First ifft', color='blue')  

pylab.ylabel('Strain')
pylab.xlabel('t')
pylab.legend()
pylab.show()

The results:

The correct result (using the backend ofifft from npfft.py; the value of hp_tilde is not changed after 'First ifft'):

The wrong result (using the backend ofifft in the fftw.py; the value of hp_tilde is changed after 'First ifft')

[ahnitz]

@josh-willis Can you look at this? It seem quite worrying.

[ahnitz]

@josh-willis I can reproduce this with the following as well (essentially, the same as posted, just minor simplifications and switched to PhenomD as PhenomA has backward Fourier conventions).

import pylab
from pycbc import types, fft, waveform

delta_t = 1 / 4096
delta_f = 1.0 / 4
f_lower = 40

hp_tilde, _ = waveform.get_fd_waveform(approximant='IMRPhenomD', mass1=6,
                                       mass2=6, delta_f=delta_f, f_lower=f_lower)

tlen = (len(hp_tilde) -1) * 2

hp=types.TimeSeries(types.zeros(tlen), delta_t=delta_t)
hp_new=hp.copy()

fft.ifft(hp_tilde, hp)
fft.ifft(hp_tilde, hp_new)

pylab.plot(hp_new.sample_times, hp_new, label='Second ifft', color='red', alpha=0.5)  
pylab.plot(hp.sample_times, hp, label='First ifft', color='blue', alpha=0.5)  

pylab.ylabel('Strain')
pylab.xlabel('t')
pylab.legend()
pylab.show()

[ahnitz]

@josh-willis However, the following does give the correct result on the same environment.

import pylab
from pycbc import types, fft, waveform

delta_t = 1 / 4096
delta_f = 1.0 / 4
f_lower = 40

hp_tilde, _ = waveform.get_fd_waveform(approximant='IMRPhenomD', mass1=6,
                                       mass2=6, delta_f=delta_f, f_lower=f_lower)

hp = hp_tilde.to_timeseries()
hp_new = hp_tilde.to_timeseries()

pylab.plot(hp_new.sample_times, hp_new, label='Second ifft', color='red', alpha=0.5)  
pylab.plot(hp.sample_times, hp, label='First ifft', color='blue', alpha=0.5)  

pylab.ylabel('Strain')
pylab.xlabel('t')
pylab.legend()
pylab.show()

[josh-willis]

The behavior observed is in fact the expected behavior for the FFTW backend, when performing a complex-to-real transformation. Per the FFTW documentation:

• FFTW_PRESERVE_INPUT: specifies that an out-of-place transform must not change its input array. This is ordinarily the default, except for c2r and hc2r (i.e. complex-to-real) transforms for which FFTW_DESTROY_INPUT is the default. In the latter cases, passing FFTW_PRESERVE_INPUT will attempt to use algorithms that do not destroy the input, at the expense of worse performance; for multi-dimensional c2r transforms, however, no input-preserving algorithms are implemented and the planner will return NULL if one is requested.

I may consider changing the behavior for the function interface (but not the class) in future development, but in the meantime, you should either make a copy of input arrays or series, or choose the numpy FFT backend if you do not care about performance, and want more intuitive behavior. This can be done via a command line interface (--fft-backends 'numpy') or through the module interface, which is (deliberately) not exposed at top level. You would need to call fft.backend_support.set_backend(['numpy']) (note you are providing a single-element list, whose only entry is the string 'numpy'; in general, you provide a preference-ordered list from among the available backends).

Though this can be confusing for novice users, we do place a high priority on having the default behavior be the one that gives good performance, which in particular means not silently disabling the performance features of the underlying numeric libraries that we call. As I said, I will consider changing the behavior for the function interface, though it is not likely to be a high priority given time constraints. So I will leave this issue open for now, but I am removing the 'bug' and 'urgent' labels.

[BenjaminDbb]

@josh-willis Thanks for your patient explanation. It really helps me!