Advanced PSF modeling

Ideally we always have plenty of well separated bright, but not oversaturated, stars to use to construct a PSF model. These models are incredibly important for certain science objectives that rely on precise shape measurements and not just total light measures. Here we demonstrate some of the special capabilities AstroPhot has to handle challenging scenarios where a good PSF model is needed but there are only very faint stars, poorly placed stars, or even no stars to work with!

import astrophot as ap
import numpy as np
import torch
from astropy.io import fits
import matplotlib.pyplot as plt
from time import time
%matplotlib inline

Making a PSF model

Before we can optimize a PSF model, we need to make the model and get some starting parameters. If you already have a good guess at some starting parameters then you can just enter them yourself, however if you don’t then AstroPhot provides another option; if you have an empirical PSF estimate (a stack of a few stars from the field), then you can have a PSF model initialize itself on the empirical PSF just like how other AstroPhot models can initialize themselves on target images. Let’s see how that works!

# First make a mock empirical PSF image
np.random.seed(124)
psf = ap.utils.initialize.gaussian_psf(2., 101, 0.5)
psf += 1e-7
psf += np.random.normal(scale = psf / 4)
psf[psf < 0] = ap.utils.initialize.gaussian_psf(2., 101, 0.5)[psf < 0]

psf_target = ap.image.PSF_Image(
    data=psf,
    pixelscale=0.5,
)

# To ensure the PSF has a normalized flux of 1, we call
psf_target.normalize()

fig, ax = plt.subplots()
ap.plots.psf_image(fig, ax, psf_target)
ax.set_title("mock empirical PSF")
plt.show()

# Now we initialize on the image
psf_model = ap.models.AstroPhot_Model(
    name = "init psf",
    model_type = "moffat psf model",
    target = psf_target,
)

psf_model.initialize()

fig, ax = plt.subplots()
ap.plots.psf_image(fig, ax, psf_model)
ax.set_title("PSF model fit to mock empirical PSF")
plt.show()

That’s pretty good! it doesn’t need to be perfect, so this is already in the right ballpark, just based on the size of the main light concentration. For the examples below, we will just start with some simple given initial parameters, but for real analysis this is quite handy.

PSF modeling without stars

Can it be done? Let’s see!

# Lets make some data that we need to fit

true_psf = ap.utils.initialize.moffat_psf(
    2., # n                                !!!!! Take note, we want to get n = 2. !!!!!!
    3., # Rd                               !!!!! Take note, we want to get Rd = 3.!!!!!!
    51, # pixels
    1.  # pixelscale
)

target = ap.image.Target_Image(
    data = torch.zeros(100,100),
    pixelscale = 1.,
    psf = true_psf,
)

true_model = ap.models.AstroPhot_Model(
    name = "true model",
    model_type = "sersic galaxy model",
    target = target,
    parameters = {
        "center": [50.,50.],
        "q": 0.4,
        "PA": np.pi/3,
        "n": 2,
        "Re": 25,
        "Ie": 1,
    },
    psf_mode = "full",
)

# use the true model to make some data
sample = true_model()
torch.manual_seed(61803398)
target.data = sample.data + torch.normal(torch.zeros_like(sample.data), 0.1)
target.variance = 0.01*torch.ones_like(sample.data)

fig, ax = plt.subplots(1,2, figsize = (16,7))
ap.plots.model_image(fig, ax[0], true_model)
ap.plots.target_image(fig, ax[1], target)
ax[0].set_title("true sersic+psf model")
ax[1].set_title("mock observed data")
plt.show()

# Now we will try and fit the data using just a plain sersic

# Here we set up a sersic model for the galaxy
plain_galaxy_model = ap.models.AstroPhot_Model(
    name = "galaxy model",
    model_type = "sersic galaxy model",
    target = target,
)

# Let AstroPhot determine its own intial parameters, so it has to start with whatever it decides automatically,
# just like a real fit.
plain_galaxy_model.initialize()

result = ap.fit.LM(plain_galaxy_model, verbose = 1).fit()
print(result.message)

Chi^2/DoF: 66.40813954210755, L: 1.0
Chi^2/DoF: 31.807673267257808, L: 0.3333333333333333
Chi^2/DoF: 25.371868148829822, L: 0.1111111111111111
Chi^2/DoF: 23.536846683022027, L: 0.012345679012345678
Chi^2/DoF: 23.522809952008853, L: 1.693508780843029e-05
Chi^2/DoF: 23.522745912573384, L: 5.645029269476763e-06
Final Chi^2/DoF: 23.522745279428626, L: 5.645029269476763e-06. Converged: success
success

# The shape of the residuals here shows that there is still missing information; this is of course
# from the missing PSF convolution to blur the model. In fact, the shape of those residuals is very
# commonly seen in real observed data (ground based) when it is fit without accounting for PSF blurring.
fig, ax = plt.subplots(1,2, figsize = (16,7))
ap.plots.model_image(fig, ax[0], plain_galaxy_model)
ap.plots.residual_image(fig, ax[1], plain_galaxy_model)
ax[0].set_title("fitted sersic only model")
ax[1].set_title("residuals")
plt.show()

# Now we will try and fit the data with a sersic model and a "live" psf

# Here we create a target psf model which will determine the specs of our live psf model
psf_target = ap.image.PSF_Image(
    data = np.zeros((51,51)),
    pixelscale = target.pixelscale,
)

# Here we create a moffat model for the PSF. Note that this is just a regular AstroPhot model that we have chosen
# to be a moffat, really any model can be used. To make it suitable as a PSF we will need to apply some very
# specific settings.
live_psf_model = ap.models.AstroPhot_Model(
    name = "psf",
    model_type = "moffat psf model",
    target = psf_target,
    parameters = {
        "n": 1., # True value is 2.
        "Rd": 2., # True value is 3.
    },
)

# Here we set up a sersic model for the galaxy
live_galaxy_model = ap.models.AstroPhot_Model(
    name = "galaxy model",
    model_type = "sersic galaxy model",
    target = target,
    psf_mode = "full",
    psf = live_psf_model, # Here we bind the PSF model to the galaxy model, this will add the psf_model parameters to the galaxy_model
)

live_psf_model.initialize()
live_galaxy_model.initialize()

result = ap.fit.LM(live_galaxy_model, verbose = 1).fit()
result.update_uncertainty()

Chi^2/DoF: 1115.302684858344, L: 1.0
Chi^2/DoF: 770.6161754844939, L: 0.3333333333333333
Chi^2/DoF: 18.38798157962339, L: 0.1111111111111111
Chi^2/DoF: 5.574257574641829, L: 0.037037037037037035
Chi^2/DoF: 1.9165797920121488, L: 0.012345679012345678
Chi^2/DoF: 1.2494113233360604, L: 0.004115226337448559
Chi^2/DoF: 0.9991516375298524, L: 0.0013717421124828531
Chi^2/DoF: 0.9929959054956735, L: 1.881676423158921e-06
Chi^2/DoF: 0.9929955231494145, L: 6.27225474386307e-07
Final Chi^2/DoF: 0.9929955231490599, L: 2.323057312541878e-08. Converged: success

print("fitted n for moffat PSF: ", live_galaxy_model["psf:n"].value.item(), "we were hoping to get 2!")
print("fitted Rd for moffat PSF: ", live_galaxy_model["psf:Rd"].value.item(), "we were hoping to get 3!")
print(live_galaxy_model.parameters)

fitted n for moffat PSF:  1.995246419826767 we were hoping to get 2!
fitted Rd for moffat PSF:  2.972340397385543 we were hoping to get 3!
galaxy model:
center: [0.0, 0.0] +- [0.1, 0.1] [arcsec], locked
n: 1.995246419826767 +- 0.012799404244268937 [none], limits: (0.1, 10.0)
Rd: 2.972340397385543 +- 0.012574434548764288 [arcsec], limits: (0.0, None)
I0: 0.0 [log10(flux/arcsec^2)], locked
center: [50.002139693462674, 49.998723271753754] +- [0.0015637301450824137, 0.001167111468237485] [arcsec]
q: 0.4006732240751394 +- 0.0007932093092872861 [b/a], limits: (0.0, 1.0)
PA: 1.0468726647855944 +- 0.0001814039797940729 [radians], limits: (0.0, 3.141592653589793), cyclic
n: 1.9926241231069444 +- 0.0058373483591286045 [none], limits: (0.36, 8.0)
Re: 24.952783031471295 +- 0.020947532222143964 [arcsec], limits: (0.0, None)
Ie: 1.0011860195123792 +- 0.0007932478943630854 [log10(flux/arcsec^2)]

This is truly remarkable! With no stars available we were still able to extract an accurate PSF from the image! To be fair, this example is essentially perfect for this kind of fitting and we knew the true model types (sersic and moffat) from the start. Still, this is a powerful capability in certain scenarios. For many applications (e.g. weak lensing) it is essential to get the absolute best PSF model possible. Here we have shown that not only stars, but galaxies in the field can be useful tools for measuring the PSF!

fig, ax = plt.subplots(1,2, figsize = (16,7))
ap.plots.model_image(fig, ax[0], live_galaxy_model)
ap.plots.residual_image(fig, ax[1], live_galaxy_model)
ax[0].set_title("fitted sersic + psf model")
ax[1].set_title("residuals")
plt.show()

There are regions of parameter space that are degenerate and so even in this idealized scenario the PSF model can get stuck. If you rerun the notebook with different random number seeds for pytorch you may find some where the optimizer “fails by immobility” this is when it get’s stuck in the parameter space and can’t find any way to improve the likelihood. In fact most of these “fail” fits do return really good values for the PSF model, so keep in mind that the “fail” flag only means the possibility of a truly failed fit. Unfortunatly, detecting convergence is hard.

PSF fitting with a faint star

Fitting a PSF to a galaxy is perhaps not the most stable way to get a good model. However, there is a very common situation where this kind of fitting is quite helpful. Consider the scenario that there is a star, but it is not very bright and it is right next to a galaxy. Now we need to model the galaxy and the star simultaneously, but the galaxy should be convolved with the PSF for the fit to be stable (otherwise you’ll have to do several iterations to converge). If there were many stars you could perhaps just stack a bunch of them and hope the average is close enough, but in this case we don’t have many to work with so we need to squeeze out as much statistical power as possible.

# Lets make some data that we need to fit

true_psf2 = ap.utils.initialize.moffat_psf(
    2., # n                                !!!!! Take note, we want to get n = 2. !!!!!!
    3., # Rd                               !!!!! Take note, we want to get Rd = 3.!!!!!!
    51, # pixels
    1.  # pixelscale
)

target2 = ap.image.Target_Image(
    data = torch.zeros(100,100),
    pixelscale = 1.,
    psf = true_psf,
)

true_galaxy2 = ap.models.AstroPhot_Model(
    name = "true galaxy",
    model_type = "sersic galaxy model",
    target = target2,
    parameters = {
        "center": [50.,50.],
        "q": 0.4,
        "PA": np.pi/3,
        "n": 2,
        "Re": 25,
        "Ie": 1,
    },
    psf_mode = "full",
)
true_star2 = ap.models.AstroPhot_Model(
    name = "true star",
    model_type = "point model",
    target = target2,
    parameters = {
        "center": [70,70],
        "flux": 2.,
    },
)
true_model2 = ap.models.AstroPhot_Model(
    name = "true model",
    model_type = "group model",
    target = target2,
    models = [true_galaxy2, true_star2],
)

# use the true model to make some data
sample2 = true_model2()
torch.manual_seed(1618033988)
target2.data = sample2.data + torch.normal(torch.zeros_like(sample2.data), 0.1)
target2.variance = 0.01 * torch.ones_like(sample2.data)

fig, ax = plt.subplots(1,2, figsize = (16,7))
ap.plots.model_image(fig, ax[0], true_model2)
ap.plots.target_image(fig, ax[1], target2)
ax[0].set_title("true model")
ax[1].set_title("mock observed data")
plt.show()

# Now we will try and fit the data

psf_model2 = ap.models.AstroPhot_Model(
    name = "psf",
    model_type = "moffat psf model",
    target = psf_target,
    parameters = {
        "n": 1., # True value is 2.
        "Rd": 2., # True value is 3.
    },
)

# Here we set up a sersic model for the galaxy
galaxy_model2 = ap.models.AstroPhot_Model(
    name = "galaxy model",
    model_type = "sersic galaxy model",
    target = target,
    psf_mode = "full",
    psf = psf_model2,
)

# Let AstroPhot determine its own intial parameters, so it has to start with whatever it decides automatically,
# just like a real fit.
galaxy_model2.initialize()

star_model2 = ap.models.AstroPhot_Model(
    name = "star model",
    model_type = "point model",
    target = target2,
    psf = psf_model2,
    parameters = {
        "center": [70,70], # start the star in roughly the right location
        "flux": 2.5,
    },
)

star_model2.initialize()

full_model2 = ap.models.AstroPhot_Model(
    name = "full model",
    model_type = "group model",
    models = [galaxy_model2, star_model2],
    target = target2,
)

result = ap.fit.LM(full_model2, verbose = 1).fit()

Chi^2/DoF: 1116.2498824099835, L: 1.0
Chi^2/DoF: 793.718437875768, L: 0.3333333333333333
Chi^2/DoF: 14.515879653006385, L: 0.1111111111111111
Chi^2/DoF: 6.287894183775265, L: 0.037037037037037035
Chi^2/DoF: 2.013731912346296, L: 0.012345679012345678
Chi^2/DoF: 1.2291545048429289, L: 0.004115226337448559
Chi^2/DoF: 1.005507474787216, L: 0.0013717421124828531
Chi^2/DoF: 1.0006501100242953, L: 1.881676423158921e-06
Chi^2/DoF: 1.0006498195221003, L: 2.0907515812876902e-07
Final Chi^2/DoF: 1.00064981946267, L: 2.323057312541878e-08. Converged: success

fig, ax = plt.subplots(1,2, figsize = (16,7))
ap.plots.model_image(fig, ax[0], full_model2)
ap.plots.residual_image(fig, ax[1], full_model2)
ax[0].set_title("fitted sersic+star model")
ax[1].set_title("residuals")
plt.show()

print("fitted n for moffat PSF: ", galaxy_model2["psf:n"].value.item(), "we were hoping to get 2!")
print("fitted Rd for moffat PSF: ", galaxy_model2["psf:Rd"].value.item(), "we were hoping to get 3!")

print("---Note that we can just as well look at the star model parameters since they are the same---")
print("fitted n for moffat PSF: ", psf_model2["n"].value.item(), "we were hoping to get 2!")
print("fitted Rd for moffat PSF: ", psf_model2["Rd"].value.item(), "we were hoping to get 3!")

fitted n for moffat PSF:  1.9624007733141735 we were hoping to get 2!
fitted Rd for moffat PSF:  2.9787398126512596 we were hoping to get 3!
---Note that we can just as well look at the star model parameters since they are the same---
fitted n for moffat PSF:  1.9624007733141735 we were hoping to get 2!
fitted Rd for moffat PSF:  2.9787398126512596 we were hoping to get 3!

Note that the fitted moffat parameters aren’t much better than they were earlier when we just fit the galaxy alone. This shows us that extended objects have plenty of constraining power when it comes to PSF fitting, all this information has previously been left on the table! It makes sense that the galaxy dominates the PSF fit here, while the star is very simple to fit, it has much less light than the galaxy in this scenario so the S/N for the galaxy dominates. The reason this works really well is of course that the true data is in fact a sersic model, so we are working in a very idealized scenario. Real world galaxies are not necessarily well described by a sersic, so it is worthwhile to be cautious when doing this kind of fitting. Always make sure the results make sense before storming ahead with galaxy based PSF models, that said the payoff can be well worth it.

PSF fitting for faint stars

Sometimes there are stars available, but they are faint and it is hard to see how a reliable fit could be obtained. We have already seen how faint stars next to galaxies are still viable for PSF fitting. Now we will consider the case of isolated but faint stars. The trick here is that we have a second high resolution image, perhaps in a different band. To perform this fitting we will link up the two bands using joint modelling to constrain the star centers, this will constrain some of the parameters making it easier to fit a PSF model.

# Coming soon

PSF fitting for saturated stars

A saturated star is a bright star, and it’s just begging to be used for modelling a PSF. There’s just one catch, the highest signal to noise region is completely messed up and can’t be used! Traditionally these stars are either ignored, or a two stage fit is performed to get an “inner psf” and an “outer psf” which are then merged. Why not fit the inner and outer PSFs all at once! This can be done with AstroPhot using parameter constraints and masking.

# Coming soon