Small Discrepancy between Real and Fourier Space Shift

[EiffL]

While working on porting new features to jax-galsim, we noticed that the GalSim phase shifts are not exactly equivalent to the real domain equivalent.

Here is an example of what I mean:

import galsim

gal = galsim.Gaussian(sigma=1.).shift([5,5]).withFlux(200)

real_im = gal.drawImage(nx=128, ny=128, scale=0.2, method='real_space')

fft_im = gal.drawImage(nx=128, ny=128, scale=0.2, method='fft')

notebook

So, it's definitely small, but not a nice pattern...

My guess... is that it might have something to do with this: https://github.com/GalSim-developers/GalSim/blob/3a9ab46b3c11333cc5def004cbce3c74f6f186cf/src/SBTransform.cpp#L749 because when we use the exact expression in jax-galsim we recover perfectly the same answer as drawing the image in real space.

Just wanted to report this to see if this is known/not a big deal. Feel free to close this issue if this is the case :-)

[rmjarvis]

I couldn't reproduce the exact result you had. On my laptop, this results in the following:

A little bit smaller residuals, and a different pattern. Apparently the rounding behavior on my laptop is somewhat different than on collab. If you use double precision though, the differences are ~1.e-15.

So I think this isn't a problem with the method per se. Just single-precision floating point errors accumulating. But they are a fair bit higher than the single-precision epsilon value, so I switched the code to use double precision for the phases, which basically fixed it. cf. #1234

[EiffL]

Hum interesting, I wouldn't have guessed platform-specific discrepancies at that level. But ok, if you recover 1e-15 using double precision then you are right it's not related to an approximation in this formula used to compute the phase.

[rmjarvis]

The first equation is nominally exact, but can have rounding errors. The line you pointed to was intended to get back to |z|=1, which it does to the accuracy of the data type. So the errors are just that numerical rounding errors can add up, maybe semi-coherently, rather than purely randomly. And apparently differently on different systems. That surprised me too, FWIW.

[rmjarvis]

Done. #1234