Using a galsim object for a PSF instead of a Gaussian Mixture

[s-Sayan]

Hi Erin,

First of all, thank you for developing such a flexible and generalizable codebase. I’ve implemented and validated ngmix within the SuperBIT pipeline (link), where we currently use the em5 model for PSF fitting. One example of the fit is shown below.

My main concern is related to PSF modeling with a gaussian mixture. Since SuperBIT is a balloon-borne experiment, we don’t have atmospheric seeing contributing to PSF blurring. The PSF is not undersampled (pixel scale ~0.141", PSF FWHM ~0.47", the stability was an issue and the telescope jitter was much higher than we thought). We perform metacalibration using the psf="dilate" option (as you previously suggested) and include PSF shear terms to correct additive biases. However, it still seems that the Gaussian mixture approximation may not capture important PSF features, particularly for exposures with highly elliptical or unstable PSFs. In those cases, we observe PSF-to-galaxy ellipticity leakage and have had to discard some exposures for metacalibration (see our recent paper: Lensing in the Blue III).

Eric Huff mentioned that ngmix may support using a GalSim object directly for the PSF, instead of a Gaussian mixture approximation. Although this would likely be slower than em5, the relatively small number of objects per SuperBIT target makes it potentially feasible for us.

I tried to locate this functionality in the codebase but wasn’t able to find it. If this option exists, could you please point me to how to use it (or whether it requires a custom implementation)?

Thanks again for your help, and for all your work on ngmix!

Best regards,
Sayan

[esheldon]

You can use a galsim fitter with a real PSF image. See ngmix.fitting.galsim_fitters

[s-Sayan]

Hi Erin,

Thanks a lot for the guidance! I am trying to use ngmix.fitting.GalsimFitter in the MetacalBootstrapper, but I'm not sure I'm doing it right — one quick hint is that it runs much faster than the EM fitter I've been using (but my understanding is it should be much slower).

My understanding is that I don't need to set psf_runner anymore. Could you please check the following implementation and let me know if it's correct?

types = ['noshear', '1p', '1m', '2p', '2m', '1p_psf', '1m_psf', '2p_psf', '2m_psf']

guesser = ngmix.guessers.TFluxAndPriorGuesser(
    rng=rng, T=T_guess, flux=flux_guess, prior=prior,
)

fitter = ngmix.fitting.GalsimFitter(model='gauss', prior=prior, fit_pars=lm_pars)

runner = ngmix.runners.Runner(
    fitter=fitter, guesser=guesser, ntry=20,
)

boot = ngmix.metacal.MetacalBootstrapper(
    runner=runner, psf_runner=None,
    rng=rng,
    psf='dilate',
    step=0.01,
    types=types,
)

resdict, obsdict = boot.go(obslist)

The obslist is an ngmix.ObsList with 36 observations, each with its own PSF, created from a MEDS file.

For comparison, here is what I currently do with the EM fitter (heavily inspired by examples in your repo). Please also let me know if this is correct as well:

em_pars = {'tol': 1.0e-6, 'maxiter': 50000}
types = ['noshear', '1p', '1m', '2p', '2m', '1p_psf', '1m_psf', '2p_psf', '2m_psf']

fitter = ngmix.fitting.Fitter(model='gauss', prior=prior, fit_pars=lm_pars)

guesser = ngmix.guessers.TPSFFluxAndPriorGuesser(
    rng=rng, T=T_guess, prior=prior,
)

runner = ngmix.runners.Runner(
    fitter=fitter, guesser=guesser, ntry=20,
)

psf_fitter = ngmix.em.EMFitter(maxiter=em_pars['maxiter'], tol=em_pars['tol'])
psf_guesser = ngmix.guessers.GMixPSFGuesser(rng=rng, ngauss=5)
psf_runner = ngmix.runners.PSFRunner(
    fitter=psf_fitter, guesser=psf_guesser, ntry=20,
)

boot = ngmix.metacal.MetacalBootstrapper(
    runner=runner, psf_runner=psf_runner,
    rng=rng,
    psf='dilate',
    step=0.01,
    types=types,
)

resdict, obsdict = boot.go(obslist)

[esheldon]

I have not used the galsim fitters with metacal bootstrapper.

I think it might make sense to do a very simple simulation test to see if you are getting the result you expect. You could start with round gaussian psf and object (before shear), with very small noise.

As for speed, it could be faster in some cases. Here you are no longer fitting a 5-gaussian model for every PSF, so maybe that was dominating the runtime. I'm not sure.

[s-Sayan]

Yeah, great idea! It seems to overestimate the pre-PSF object size. For example, inputting a Gaussian with σ = 0.4" (T = 2σ² = 0.32 arcsec²) convolved with a PSF of FWHM = 0.5", I get back T ≈ 0.47 arcsec² instead.

The prior is simple (link). Here's a minimal reproducer — let me know if I'm doing something wrong:

flux = 1
sigma = 0.4
scale = 0.141
npix = 53
nse_sd = 1e-4
seed = 42
psf_fwhm = 0.5

gal = galsim.Gaussian(sigma=sigma, flux=flux)
psf = galsim.Gaussian(fwhm=psf_fwhm)
obj = galsim.Convolve(gal, psf)

gal_image = gal.drawImage(nx=npix, ny=npix, scale=scale).array
psf_image = psf.drawImage(nx=npix, ny=npix, scale=scale).array
obj_im = obj.drawImage(nx=npix, ny=npix, scale=scale).array

rng = np.random.RandomState(seed=seed)
nse_im = rng.normal(size=obj_im.shape, scale=nse_sd)
image = obj_im + nse_im

cen = npix // 2
jac = ngmix.jacobian.DiagonalJacobian(scale=scale, row=cen, col=cen)

psf_obs = ngmix.Observation(image=psf_image, jacobian=jac)
obj_obs = ngmix.Observation(
    image=image,
    noise=nse_im,
    weight=np.ones_like(nse_im) / nse_sd**2,
    jacobian=jac,
    bmask=np.zeros_like(nse_im, dtype=np.int32),
    ormask=np.zeros_like(nse_im, dtype=np.int32),
    psf=psf_obs,
)

# Initial guess via GaussMom
gm = GaussMom(1.2).go(obj_obs)

guesser = ngmix.guessers.TFluxAndPriorGuesser(
    rng=rng, T=gm['T'], flux=gm['flux'], prior=prior,
)
lm_pars = {'maxfev': 2000, 'xtol': 5.0e-5, 'ftol': 5.0e-5}
fitter = ngmix.fitting.GalsimFitter(model='gauss', prior=prior, fit_pars=lm_pars)
runner = ngmix.runners.Runner(fitter=fitter, guesser=guesser, ntry=20)

resdict = runner.go(obj_obs)
print(resdict['pars'][4])

[s-Sayan]

And just to compare, using the usual Fitter with a PSFRunner via Bootstrapper, I do get the correct pre-PSF size back (~0.32):

fitter = ngmix.fitting.Fitter(model="gauss", prior=prior, fit_pars=lm_pars)
guesser = ngmix.guessers.TPSFFluxAndPriorGuesser(rng=rng, T=T_guess, prior=prior)
runner = ngmix.runners.Runner(fitter=fitter, guesser=guesser, ntry=20)

psf_lm_pars = {'maxfev': 4000, 'xtol': 5.0e-5, 'ftol': 5.0e-5}
psf_fitter = ngmix.fitting.Fitter(model='gauss', fit_pars=psf_lm_pars)
psf_guesser = ngmix.guessers.SimplePSFGuesser(rng=rng)
psf_runner = ngmix.runners.PSFRunner(fitter=psf_fitter, guesser=psf_guesser, ntry=20)

boot = ngmix.bootstrap.Bootstrapper(runner=runner, psf_runner=psf_runner)
resdict = boot.go(obj_obs)
print(resdict['pars'][4])
# >> np.float64(0.32032405584731727)

[esheldon]

the galsim fitter works in half light radius rather than T.

You can check with the input by doing fwhm_meas = resdict['pars'][4] * 2

[s-Sayan]

Okay that matches perfectly, thanks! I also verified this works with MetacalBootstrapper on the same single round object — res['pars'][4] is consistent with the metacal quantities, and R11/R22 are both ~1 as expected (done with 100 objects with different noise):

R11: mean=0.9993, std=0.0059
R22: mean=0.9999, std=0.0059

Using psf_runner=None with psf='dilate':

types = ['noshear', '1p', '1m', '2p', '2m', '1p_psf', '1m_psf', '2p_psf', '2m_psf']
guesser = ngmix.guessers.TFluxAndPriorGuesser(rng=rng, T=T_guess, flux=flux_guess, prior=prior)
fitter = ngmix.fitting.GalsimFitter(model='gauss', prior=prior, fit_pars=lm_pars)
runner = ngmix.runners.Runner(fitter=fitter, guesser=guesser, ntry=20)
boot = ngmix.metacal.MetacalBootstrapper(
    runner=runner, psf_runner=None,
    rng=rng, psf='dilate', step=0.01, types=types,
)
resdict, obsdict = boot.go(obj_obs)

Would you recommend any more robust simple test for this setup, before I implement this in our pipeline check it with image simulations? Also, would GalsimFitter + MetacalBootstrapper work with an ObsList where observations have different PSFs and Jacobians (SuperBIT has a rotational dither between exposures)?
Thanks a lot again -- this has already been super helpful!

[esheldon]

Yes, it should work with an ObsList

You should be able to do a low noise, round objects test with metacal that runs fast and should come out unbiased. I'd recommend that before diving into a big simulation.

[s-Sayan]

Thanks @esheldon! I ran the test you suggested, comparing the GalsimFitter against the same EMFitter (5-gauss PSF) runner as before, with MetacalBootstrapper. Looks consistent with EM.

Setup: 10,000 round Gaussian objects, HLR drawn uniformly in 0.1–1 arcsec, flux=1, no shear. Each object's ObsList has 36 observations, each convolved with a different SuperBIT PSF drawn from its FOV position (matching real multi-exposure data; median PSF FWHM ≈ 0.46″) and a different noise realization (noise_sd=1e-4).

Timing: GalSimFitter was ~6.6× faster than EMFitter (1.65 s/obj vs 10.92 s/obj over 10,000 objects).

Checking: HLR recovery, g1/g2 (should be consistent with zero), R11/R22, and a PSF→galaxy ellipticity leakage plot. Results below.

[esheldon]

This looks fine.

But by unbiased I meant that shear is unbiased. Even if there is some bias in the properties you measure, the shear should be well calibrated

[s-Sayan]

That makes more sense, thanks! Following your metacal example, I used the galsim fitter for a Moffat PSF with a slight shear (code):

S/N: 84551.9
R11: 0.923006
m: -0.000105781 +/- 0.000413673 (99.7% conf)
c: 5.28524e-07 +/- 4.41977e-06 (99.7% conf)

Running the identical code with the SuperBIT PSF (code) gives:

S/N: 68002
R11: 0.929858
m: -0.00528693 +/- 0.000396201 (99.7% conf)
c: -6.91925e-05 +/- 3.99757e-06 (99.7% conf)

I noticed your metacal bootstrapper uses psf=gauss. Switching to psf='dilate' blows up m for the superbit psf case:

S/N: 43381.2
R11: 0.944036
m: -0.143113 +/- 0.000325965 (99.7% conf)
c: -0.000921407 +/- 3.30204e-06 (99.7% conf)

But with the EMFitter (code), psf='dilate' does noticeably better:

S/N: 61077.6
R11: 0.957135
m: -0.0650563 +/- 0.00337187 (99.7% conf)
c: -0.00372223 +/- 3.93602e-05 (99.7% conf)

Any thoughts on what's driving the difference?

[esheldon]

The m bias you see might actually be partly or fully from additive bias in the g1 direction (c1). As I recall the superbit PSF is complex and not round.

Note that the c here is only coming from looking at g2 and m is from looking at g1. If there is additive in g1 then it would look like an m.

So dilate leaves in PSF shape that psf='fitgauss' or psf='gauss' would not. Of course we think fitgauss/gauss have different problems due to the very non gaussian nature of the superbit PSF, but they will suffer less additive.

A quick test would be to run without any shear and see if you are still seeing a mean in shear1, shear2.

Even better, you can isolate this by running with a range of shears, say [-0.02, -0.01, 0.0, 0.01, 0.02]

If there is additive it will look like a line with an offset at zero (c), and the slope will be (1 + m)

If additive is confirmed for dilate, which would not be surprising, the next step would be to implement the PSF response and subtract off the additive.

[s-Sayan]

Yes, that makes perfect sense! I just did this, you're absolutely right that the large $m$ was due to mismodelling, i.e. assuming there is no additive bias. Running multipoint fitting for both $g_1$ and $g_2$:

Now, to implement this PSF related correction — if I understand correctly, I compute the mean shear following Huff & Mandelbaum 2017:

$$ \left<\boldsymbol{\gamma}\right> \approx \left<\boldsymbol{R}^{\gamma}\right>^{-1} \Big(\left<\boldsymbol{e}\right> - \left<\boldsymbol{c}^{\text{total}}\right>\Big) $$

where the per-object total additive bias is:

$$ \boldsymbol{c}^{\text{total}} = \boldsymbol{c}^{\gamma} + \boldsymbol{c}^{\text{PSF}} $$

with

$$ \boldsymbol{c}^{\gamma} = \frac{\boldsymbol{e}^+ + \boldsymbol{e}^-}{2} - \boldsymbol{e}_{\text{noshear}}, \qquad \boldsymbol{c}^{\text{PSF}} = \frac{\boldsymbol{e}^{\text{PSF},+} + \boldsymbol{e}^{\text{PSF},-}}{2} - \boldsymbol{e}_{\text{noshear}} $$

So I add 1p_psf, 1m_psf, 2p_psf, 2m_psf to the types in the MetacalBootstrapper, compute the mean $\boldsymbol{c}^{\text{total}}$ per object, subtract it from each shape, then divide by $\boldsymbol{R}^{\gamma}$ — is that the right procedure?

Note: the above result uses the GalSimFitter in the MetacalBootstrapper (with psf='dilate') run on single observations, where every object's PSF is same, rendered from one single position using a SuperBIT .psf file.

[esheldon]

For the second term you need to subtract Rpsf * epsf where Rpsf is the response of the measurement to psf shearing.

Rpsf = (e^psf+ - e^psf-) / (2 delta_gamma_psf)

where here gamma_psf is how much you sheared the PSF (say 0.01)

The idea is this: you know some of the psf epsf leaked into your measurement. But you don't just subtract epsf you need to know how much your measurement responded to that, so it gets weighted Rpsf * epsf