Multi-component PSF models?

[entaylor]

I'm trying to get a pipeline running for PSF characterisation and then PSF photometry to calibrate a giant set of VISTA VIRCam images. Of course once you start looking, the PSF is more complicated than previously realised. A 2D Moffat profile does a pretty good job, but there is a concentrated central core that a single Moffat profile misses.

Is it feasible/easy to construct a 2-component PSF model; e.g. a double 2D Moffat or a 2D Moffat + 1D or 2D Gaussian?

I have tried a bunch of stuff to refine and validate how i am doing the PSF modelling to make sure, e.g., that the fits aren't being dragged off to a wide-winged solution by bad background subtraction; i think the PSF really just does have significant departures from an ordinary Moffat shape. I could try a 'non-parameteric' pixelated psf model, but the data are shallow and i typically only have a few dozen high S/N stars per chip, which i think is not sufficient to avoid excess noise in the PSF model. I could try something like modelling every star as a Gaussian that is constrained to have the same shape parameters for all stars, but then i would be convolving the Gaussian (source component) with the Moffat (PSF component), which i don't think is what i want to do.

I've attached a few diagnostic plots to show the PSFs that i'm considering. The first and second panel show the image and post-fit residuals on the same arcsinh gray-scale. In each, the red contours illustrate the fit (2D Moffat; thanks again for that, btw!) PSF in log-10 spaced contours. The right panel shows the pixel-by-pixel profile for each star: the black shows the image, the red show the PSF fit, the blue shows the difference between the two. In each panel, you can see that the best fit Moffat profile is low by about 30-50 % for the central few pixels, and then it overshoots slightly at intermediate radii to compensate for the max L/min Chi2. These plots make me think that i can get a really excellent PSF description with a 2 component model.

[ConnorStoneAstro]

Hi @entaylor very cool analysis! Yes it is possible to have multi component PSF models. Looking through the documentation I realize its a major oversight that I don't mention this! But a psf group model should behave basically the same as a group model, but it is intended to hold psf model type objects. Currently the only thing in the docs is the source code: https://astrophot.readthedocs.io/en/latest/_modules/astrophot/models/group_psf_model.html#PSF_Group_Model

I will certainly need to add that to the PSF tutorials. I think I can add it at the bottom of the PSF Model Basics tutorial since it is likely something many others will want to do!

Let me know if you run into any hurdles. I'll make sure to at least add a basic example within the next couple weeks so lets keep this issue open until I get that done!

[ConnorStoneAstro]

Oh, one thing I can pre-empt that you will need to do. By default the I0 parameter is fixed for moffat PSF models, you will need to set it free for one of the components. It is unnecessary for single PSF models since they are always normalized to a total flux of 1. But with two components, you need some way to give the two components a relative brightness, so one of the I0 parameters should be un-fixed. You don't need both of them free since only the relative value is what matters. In principle you can also free the position of the models to handle cases where the two components are not perfectly aligned, but I doubt that will be necessary.

[ConnorStoneAstro]

Ok, this is why I need to make a tutorial haha. I also remembered that for each psf model object, you will need to set normalize_psf = False otherwise they will be individually normalized to a total flux of 1, rather than adding together first. Currently the psf group model doesn't normalize the total flux to 1, which is something I can work on while building the demo (though it isn't critical).

[entaylor]

Connor, you are incredible. Thank you for this lightning fast response. I have something running, but i think i haven't quite got it 100% right yet. With luck i'll have some time to work more on it tomorrow, and i hope i will be able to share with you some nice results. Thank you again!

[entaylor]

Actually ... i couldn't keep myself from tinkering a little bit further still ... i think i'm really close, but i haven't quite got the fixed vs free handling of the I0s right yet. More on this to come in the next day or three, though! : )

[ConnorStoneAstro]

Nice! And that's my bad for using the wrong word, you need to set locked = False so in the parameters dictionary when creating the object you would include arguments like this:

parameters = {
    "I0": {"value": 0.0, "locked": False},
},
normalize_psf = False,

I have been meaning to fix the confusing fixed for models and locked for parameters for a while, but that would be a breaking change. I am planning a big breaking update for when I start my next postdoc where I will have more dedicated time to work on AstroPhot. So if you think of anything that could use some polishing, let me know and I'll see if I can fix it in the new version!

[ConnorStoneAstro]

Hi @entaylor I put together a quick initial version of the psf group model tutorial to help give you an idea of what things should look like. Along the way I am doing some minor fixes (like making the psf group model psfs normalized). The tutorial is far from finished, but I'm sending you the preview of what it will look like so hopefully it can help you!

https://astrophot--253.org.readthedocs.build/en/253/tutorials/PSFModelling.html#multi-component-psf-models-psf-group-model

[entaylor]

Hi again ... i've had a little time to try to work this through carefully, and i still have something going awry that i can't quite diagnose. For completeness, i'm using astrophot v0.16.9 for what follows.

Below is the code that i have, closely following your (very helpful!) new documentation, for a small test psf_data.npz.

The difficulty i have is that i think i have not got the component vs group normalisation issues working correctly. The output results from the model looks like:

psf star:
center: [0.0, 0.0] +- [0.1, 0.1] [arcsec], locked
n: 1.7357481961978787 +- 0.003602003871550701 [none], limits: (0.1, 10.0)
Rd: 1.1905451097711186 +- 0.0038431652596877548 [arcsec], limits: (0.0, None)
**I0: -0.4144167529393543 +- 5938.661913045064 [log10(flux/arcsec^2)]**
q: 0.8910732629316971 +- 0.0007274614521312713 [b/a], limits: (0.0, 1.0)
PA: 0.09075197937667052 +- 0.003860378652493767 [radians], limits: (0.0, 3.141592653589793), cyclic
center: [0.0, 0.0] +- [0.1, 0.1] [arcsec], locked
n: 9.256166020731236 +- 1.183105575872009e+18 [none], limits: (0.1, 10.0)
Rd: 0.09969753145832971 +- 2614508676407232.5 [arcsec], limits: (0.0, None)
**I0: -9.001603209348824 +- 3.886882872358984e+16 [log10(flux/arcsec^2)]**
q: 0.002826848212200517 +- 112186900232551.8 [b/a], limits: (0.0, 1.0)
PA: 0.8928834226567623 +- 459494881142817.1 [radians], limits: (0.0, 3.141592653589793), cyclic
center: [20.552189247462, 20.49030332918398] +- [0.0005033964604892446, 0.00048695407016306877] [arcsec]
**flux: -0.3904478587001068 +- 5938.661913045068 [log10(flux)]**

Note the fluxes appear to be unconstrained?

The code i have used for this follows.

import astrophot
import numpy as np
from matplotlib import pyplot as plt
import torch

# load some PSF data
data = np.load( 'psf_data.npz' )
psf_img = data[ 'psf_img' ]
psf_var = data[ 'psf_var']

# create a PSF_Image target
psf_image_target = astrophot.image.PSF_Image( 
    data=psf_img, 
    pixelscale=1. )

# create two Moffat components
moffat1 = astrophot.models.AstroPhot_Model(
    name='psf moffat1', 
    model_type='moffat2d psf model',
    target=psf_image_target,
    parameters={'I0' : {'value':-0.1, 'locked':False} },
    normalize_psf=False, 
)
moffat1.initialize()

moffat2 = astrophot.models.AstroPhot_Model(
    name='psf moffat2', 
    model_type='moffat2d psf model',
    target=psf_image_target,
    parameters={'I0' : {'value':-2, 'locked':False} },
    normalize_psf=False, 
)
moffat2.initialize()

# combine the two into a single group model
psf_group_model = astrophot.models.AstroPhot_Model(
    name="double moffat",
    model_type="psf group model",
    target=psf_image_target,
    models=[moffat1, moffat2],
    normalize_psf=True,
)
psf_group_model.initialize()

fig, ax = plt.subplots(1, 2, figsize=(13, 5))
# astrophot.plots.psf_image(fig, ax[0], psf_group_model)
ax[0].set_title("PSF model initialization")
astrophot.plots.residual_image(fig, ax[1], psf_group_model, 
                               normalize_residuals=torch.tensor(psf_var))
ax[1].set_title("initial residuals")
plt.show()

target_image = astrophot.image.Target_Image(
    data=psf_img, variance=psf_var, pixelscale=1.
)

model = astrophot.models.AstroPhot_Model(
    name="psf star",
    model_type="point model",
    target=target_image,
    psf=psf_group_model,
)
model.initialize()

result = astrophot.fit.LM(model, verbose=0).fit()
result.update_uncertainty()
print( model )

fig, ax = plt.subplots(1, 2, figsize=(13, 5))
# astrophot.plots.psf_image(fig, ax[0], psf_group_model)
ax[0].set_title("PSF model fit to mock empirical PSF")
astrophot.plots.residual_image(fig, ax[1], psf_group_model, 
                               normalize_residuals=torch.tensor(psf_var))
ax[1].set_title("residuals")
plt.show()

One bit of funniness is that if i try the astrophot.plots.psf_image with the group model, i get the error:
AttributeError: 'str' object has no attribute 'pixel_shape'
with the traceback pointing to:

File ~/work/photometry/.venv/lib/python3.9/site-packages/astrophot/image/window_object.py:467, in Window.__eq__(self, other)
    464 @torch.no_grad()
    465 def __eq__(self, other):
    466     return (
--> 467         torch.all(self.pixel_shape == other.pixel_shape)
    468         and torch.all(self.pixelscale == other.pixelscale)
    469         and (self.projection == other.projection)
    470         and (
    471             torch.all(
    472                 self.pixel_to_plane(torch.zeros_like(self.reference_imageij))
    473                 == other.pixel_to_plane(torch.zeros_like(other.reference_imageij))
    474             )
    475         )
    476     )

I will try to keep tinkering on this over the next few days, but posting now in case you can see something obvious that i am missing?

[ConnorStoneAstro]

Hi @entaylor ah yes, looks like you are doing things right, but the fact that the psf group model doesn't normalize it's psf image is causing the flux to be unconstrained. I will treat this as a bug and have it fixed by the end of today! As for the plots.psf_image plotting bug, I know what is causing that and I'll fix it in the same update today.

[ConnorStoneAstro]

Hi @entaylor the new version 0.16.10 should have the appropriate fixes in it for your flux values to make more sense! Now the point model flux will be what you are after, and the two I0 values will only really be constrained in ratio (so you can actually leave one fixed if you like). Let me know if you run into any other problems!

[entaylor]

Hi again Conner — Good news and bad news.

The good news is that i do have successful multiple moffat2d PSF fits! They are glorious!! I'm very excited!!

The bad news is that i think i have a new bug, which may or may not be related to these changes. I'm continuing this conversation here, with apologies if this should be redirected to a separate issue.

I'm uploading a example_data.npz and an ipython notebook to show the problem.

But in a nutshell: I have some data, and some good-ish stars. In a loop, i can fit each star one-by-one to get the best fit moffat2d profile for each. Previously (unless i have inadvertently changed something, but i'm 98% sure i haven't) i could just take the list of all of these model objects, and then fit a single profile to the ensemble of all stars as a group model. I hope that makes sense!

The critical line is the one that constructs the group model:

    model = astrophot.models.AstroPhot_Model(
        name='psf star models',
        model_type='group model',
        psf_mode='full',
        models=models_to_fit,
        target=target
    )

Where previously i think this has worked, now i get the traceback and error:

    283 @psf_mode.setter
    284 def psf_mode(self, value):
    285     self._psf_mode = value
--> 286     for model in self.models.values():
    287         model.psf_mode = value

AttributeError: 'list' object has no attribute 'values'

Am i doing something dumb, or has something changed here?

Thank you again for all your help!

Example data and notebook are attached.

psffitting.zip

[ConnorStoneAstro]

Hi @entaylor that's great to hear about the psf group models!

Indeed this bug was introduced in the previous fix. It's an easy one for me to solve, I just need to initialize the models list before updating parameters like psf_mode which try to update all the models to keep things in sync within the group model. I have already done the fix, but it will take a little while for a new version to be released as I have to run a bunch of checks.

[ConnorStoneAstro]

New version 0.16.11 is live! It has a fix for the bug and so you should be all set to go! Let me know if anything else pops up. I really appreciate everyone who sends me bug reports like this so it can be fixed for everyone!