spergel accuracy

[beckermr]

We seem to have a systematic deviation of our spergel profiles form the galsim ones. This looks like a bug to be honest.

cc @ismael-mendoza @jecampagne

[beckermr]

The k-value differences are due to galsim using fast, approximate exponential, log, and power functions:

https://github.com/GalSim-developers/GalSim/blob/ee177aa50424dbbc5cec4a71cadd6297109c2aa1/src/SBSpergel.cpp#L31-L35

I do not think we want to replicate this in JAX-GalSim, as we should let the XLA backend determine how to compute these special functions.

[beckermr]

@ismael-mendoza Did your code render in real space without a pixel, or in Fourier space with a pixel? If it is in Fourier space, then we should document the k-space behavior, but there is nothing to actually do.

I will look into the real space behavior next.

[ismael-mendoza]

thanks Matt! The behavior I was observing was with FFT drawing with pixel. Here is the test in question in case it helps:

params = {
    "flux_b": 1,
    "hlr_b": 0.4,
    "q_b": 1.0,
    "beta": 0.0,
}
nu = -0.6
slen = 101
psf = galsim.Gaussian(half_light_radius=0.7, flux=1.0)
xpsf = jgs.Gaussian(half_light_radius=0.7, flux=1.0)
fft_size = 256  # checked that increasing this does not change residual

# galsim
bulge = galsim.Spergel(nu=nu, flux=params["flux_b"], half_light_radius=params["hlr_b"])
bulge = bulge.shear(q=params["q_b"], beta=params["beta"] * galsim.degrees)
gal = galsim.Convolve([bulge, psf])
arr_galsim = gal.drawImage(nx=slen, ny=slen, scale=0.2, dtype=np.float64).array
# print(gal.getGoodImageSize(0.2)) ==> result is less than 101

# jax galsim
gsparams = jgs.GSParams(minimum_fft_size=fft_size, maximum_fft_size=fft_size)
bulge = jgs.Spergel(nu=nu, flux=params["flux_b"], half_light_radius=params["hlr_b"])
bulge = bulge.shear(q=params["q_b"], beta=params["beta"] * jgs.degrees)
gal = jgs.Convolve([bulge, xpsf]).withGSParams(gsparams)
arr_jgs = gal.drawImage(nx=slen, ny=slen, scale=0.2, dtype=jnp.float64).array

# comparison
res = arr_galsim - arr_jgs
vmax = np.abs(res).max()
vmin = -vmax

plt.imshow(res, cmap="RdBu_r", vmin=vmin, vmax=vmax)
plt.colorbar()

[beckermr]

Ok. I finally found the issue. I put in an approximation for the half-light radius of the Spergel profile. This is used to compute the scale radius for the profile. Even though that approximation is good to about 1e-8 or so, the residual inaccuracies cause the features above and thus the rendering differences. Without this approximation the profile is too slow when initialized with a half-light radius. So that leaves us in a bit of a bind, since I cannot seem to find anything more accurate that is as fast.