Loss of precision

[fredRos]

Loss of precision

[fredRos]

Oh and the disagreement is stronger with increasing image size. Either because prob. of mismatch is constant and we just have more trials, or it is really due to summing more terms in a different way

[mtazzari]

I am looking into this right now. I did not realize that sample() and chi2() tests were not passing already at the commit of 19 days ago.
I am puzzled from the fact that all the steps pass if tested individually, but sample() and chi2() - which do all the steps - fail.

[fredRos]

I think the difference is really due to shift FFT shift operation for single precision. Using atol just masks this problem because it allows tests to pass when numbers happen to be close to zero and if subsequent operations move values away from zero, our stricter rtol causes the test to fail.

conclusion

numpy and FFTW vary quite a bit in single precision and there is nothing we can do about it. We should recommend users to use double for scientific work. Single precision at your own risk!

Investigation of discrepancies

def test_shift_fft_shift(size, complex_type, rtol, atol, acc_lib):

    reference_image = create_reference_image(size=size, kernel='gaussian', dtype=complex_type)
    cpu_shift_fft_shift = np.fft.fftshift(np.fft.fft2(np.fft.fftshift(reference_image)))

    # create a copy of reference_image because galario makes in-place FFT
    ref_complex = reference_image.copy()
    acc_lib.fftshift_fft2d_fftshift(ref_complex)

    print()
    print('%.15e' % cpu_shift_fft_shift.real[0,0])
    print('%.15e' % ref_complex.real[0,0])
    print('---')

Let's test for different image sizes

size = 1024

# single precision
-1.599041643203236e-03
-1.560211181640625e-03
---
.
# double
-1.599940267100308e-03
-1.599940267100308e-03
---

So in double precision, all digits agree but in single precision, only 2 digits match. This relative mismatch of ~2.5 % gets propagated. Numpy's single fft is good to 3 digits vs double. I know that shift just copies values so the question is why FFTW (and cuFFT) are so far off in single precision.

size = 1000

1.779161626470635e-03
1.781748607754707e-03
---
.
1.780879036946570e-03
1.780879036946813e-03

Now numpy a bit more off in single. Double results not identical.

size = 1500

3.309356397994316e-04
3.327933081891388e-04
---
F
3.327185188833359e-04
3.327185189048832e-04

And now numpy is far off in single.