spergel accuracy #257

Activity

1. beckermr commented on Jun 17, 2026

   

   beckermr

   on Jun 17, 2026

   CollaboratorAuthor

   More actions

   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.

   👍React with 👍1Reacted by Ismael Mendoza
2. beckermr commented on Jun 17, 2026

   

   beckermr

   on Jun 17, 2026

   CollaboratorAuthor

   More actions

   @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.

   👍React with 👍1Reacted by Ismael Mendoza
3. ismael-mendoza commented on Jun 17, 2026

   

   ismael-mendoza

   on Jun 17, 2026

   Collaborator

   More actions

   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()

   ❤️React with ❤️1Reacted by Matthew R. Becker
4. beckermr commented on Jun 18, 2026

   

   beckermr

   on Jun 18, 2026

   CollaboratorAuthor

   More actions

   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.
5. beckermr commented on Jul 10, 2026

   

   beckermr

   on Jul 10, 2026

   CollaboratorAuthor

   More actions

   These lines might be related:

   https://github.com/GalSim-developers/GalSim/blob/ee177aa50424dbbc5cec4a71cadd6297109c2aa1/src/SBSpergel.cpp#L31-L35
6. beckermr commented on Jul 15, 2026

   

   beckermr

   on Jul 15, 2026

   CollaboratorAuthor

   More actions

   It is the FFT size used for rendering. See the test I added to #258
7. beckermr commented on Jul 15, 2026

   

   beckermr

   on Jul 15, 2026

   CollaboratorAuthor

   More actions

   It is also the scale radius comp as well to be clear. That is fixed in #258
8. beckermr

   closed this as completedon Jul 19, 2026