Image simple galaxy sheared: Real_space vs FFT

[jecampagne]

Hello, Here is a small snippet done with @aguinot (here on Colab)

import galsim
import matplotlib.pyplot as plt
print("Galsim version:", galsim.__version__)  # 2.4.7

gal_flux=1e5
g1=0.2
g2=0.1
gal_sigma = 1

pixel_scale = 0.2

gal_high_flux = galsim.Gaussian(sigma=gal_sigma, flux=gal_flux).shear(g1=g1, g2=g2)
gal_low_flux = galsim.Gaussian(sigma=gal_sigma, flux=1.).shear(g1=g1, g2=g2)

img_real_hf = gal_high_flux.drawImage(method='real_space', scale=pixel_scale)
img_fft_hf = gal_high_flux.drawImage(method='fft', scale=pixel_scale)

img_real_lf = gal_low_flux.drawImage(method='real_space', scale=pixel_scale)
img_fft_lf = gal_low_flux.drawImage(method='fft', scale=pixel_scale)

plt.figure()
plt.title(f"flux={gal_flux}\ngalsim_fft - galsim_real")
plt.imshow(img_fft_hf.array-img_real_hf.array)
plt.colorbar()
plt.show()
plt.figure()
plt.title(f"flux=1\ngalsim_fft - galsim_real")
plt.imshow(img_fft_lf.array-img_real_lf.array)
plt.colorbar()
plt.show()

We were wandering if these wiggling patterns are expected and what about their magnitude too. Are they due to some code optimisation which makes some hidden truncations/approximations?

[rmjarvis]

Yes, there are various numerical approximations made throughout GalSim, which are usually controllable via GSParams objects. Doing all the math to machine precision is overkill for almost all use cases. Our default accuracy/speed trade-off tries to target plausible accuracy requirements from typical weak lensing surveys (e.g. LSST). But we allow users to tune things differently, either for more accuracy or faster rendering. cf GSParams docs

It's not always obvious a priori which parameter is the most relevant for any particular imaging artifact. Before I tested, I expected maxk_threshold to be the one, since that controls the size of the FFT image that gets made. But in fact, in this case kvalue_accuracy is the relevant one. If you lower that from the default of 1.e-5 to 1.e-8, here is what I get back:

[jecampagne]

Thanks !

[rmjarvis]

BTW, the easiest way to adjust the GSParams of an object that is already made is to use the withGSParams method. E.g.

gal_high_flux = gal_high_flux.withGSParams(kvalue_accuracy=1.e-8)

You can include multiple parameters here. Or do it multiple times with different parameters if you prefer. This can make it easier to comment in or out particular choices to see which ones matter.

If you know what you want in advance, you can also create a gsparams instance and pass it as a parameter to the GSObject (Gaussian in this case) when you initially make it.