IntegratedGaussianPRF does not update sigma in residual images, and nor will any model which uses parameters other than x, y, flux

[Onoddil]

Any model fit in PSF photometry, passed through image = subtract_psf(image, self.psf_model, param_table, subshape=self.fitshape) will only have its x and y coordinates and flux updated. subtract_psf -- used to create the fit-subtracted residual images -- will only find the parameter names of x, y, and flux via xname, yname, fluxname = _extract_psf_fitting_names(psf), updating the PSF to subtract accordingly with

getattr(psf, xname).value = row['x_fit']
getattr(psf, yname).value = row['y_fit']
getattr(psf, fluxname).value = row['flux_fit']

thus any extra parameter -- such as sigma in IntegratedGaussianPRF -- will be left at its default values, creating incorrect residual images, affecting any iterative PSF photometry routines.

Minimum working test example:

from photutils import IntegratedGaussianPRF
from photutils.psf import BasicPSFPhotometry, DAOGroup
import matplotlib.pyplot as plt
from astropy.table import Table
from astropy.visualization import simple_norm
import numpy as np

psf = IntegratedGaussianPRF(flux=1, sigma=0.6)
psf.sigma.fixed = True
psf.flux.fixed = True

star_xy_orig = np.random.uniform(-0.5, 0.5, size=2)
xs_orig, ys_orig = np.arange(-12, 12.1, 1), np.arange(-12, 12.1, 1)
xs_orig_mg, ys_orig_mg = np.meshgrid(xs_orig, ys_orig, indexing='xy')
psf_orig = psf.evaluate(xs_orig_mg, ys_orig_mg, psf.flux, star_xy_orig[0], star_xy_orig[1],
                        psf.sigma)

psf_g = IntegratedGaussianPRF(flux=1, sigma=1.2)
psf_g.sigma.fixed = False
psf_g.flux.fixed = False

# drop one from each side of the grid for fitting due to the subtract_psf add_array bug
fitshape_orig = tuple(qq for qq in psf_orig.shape)
daogroup = DAOGroup(crit_separation=8)
basic_phot_orig = BasicPSFPhotometry(group_maker=daogroup, psf_model=psf_g,
                                     fitshape=fitshape_orig, bkg_estimator=None)
y0_orig, x0_orig = tuple((i - 1) / 2 + pp for i, pp in zip(xs_orig_mg.shape,
                                                           star_xy_orig.flatten()))
results = basic_phot_orig(psf_orig, init_guesses=Table(names=['x_0', 'y_0',
                          'flux_0'], data=[[x0_orig], [y0_orig], [1]]))

manual_resid_bad = psf_orig - psf_g.evaluate(xs_orig_mg, ys_orig_mg, results['flux_fit'][0],
                                             star_xy_orig[0], star_xy_orig[1], psf_g.sigma)
manual_resid_good = psf_orig - psf_g.evaluate(xs_orig_mg, ys_orig_mg, results['flux_fit'][0],
                                              star_xy_orig[0], star_xy_orig[1],
                                              results['sigma_fit'][0])
fig, axs = plt.subplots(2, 2, figsize=(9, 9))

for ax, p, t in zip(axs.flatten(), [psf_orig, basic_phot_orig.get_residual_image(),
                                    manual_resid_bad, manual_resid_good],
                    ['PSF Response', r'get\_residual\_image() residual',
                     'Original sigma fit residual', 'Correct sigma fit residual']):
    norm = simple_norm(p, 'log', min_cut=0)
    img = ax.imshow(p, origin='lower', cmap='viridis', norm=norm, aspect='auto')
    cb = plt.colorbar(img, ax=ax, use_gridspec=True)
    cb.set_label(t)
    ax.set_xlabel('x / pixel')
    ax.set_ylabel('y / pixel')

plt.tight_layout()
plt.show()

print(results)