Examples
Fitting a PSF
Here is a brief example which shows how to construct a loss function for fitting a PSFModel
to some data.
using PSFModels
using PSFModels: fit
using HCIDatasets: BetaPictoris
using Plots
using Statistics
# convenience function for plotting
function imshow(data; kwargs...)
xlim = extrema(axes(data, 1))
ylim = extrema(axes(data, 2))
heatmap(transpose(data); xlim=xlim, ylim=ylim,
aspect_ratio=1, clims=(1e-5, Inf), kwargs...)
end
# get a PSF from HCIDatasets.jl;
# you may be prompted to download the file
psf = BetaPictoris[:psf]
inds = CartesianIndices(psf)
imshow(psf)
We can fit this data with a variety of models, here showcasing the flexible PSFModels.fit
function.
Gaussian
Using gaussian
params = (x=20, y=20, fwhm=5, amp=0.1)
P_gauss, mod_gauss = fit(gaussian, params, psf)
pairs(P_gauss)
pairs(::NamedTuple) with 4 entries: :x => 20.0191 :y => 20.01 :fwhm => 4.80293 :amp => 0.100216
plot(
imshow(psf, title="Data"),
imshow(mod_gauss.(inds), title="Model"),
cbar=false,
ticks=false,
xlabel="",
ylabel="",
layout=2,
size=(600, 300)
)
and now using a rotated, elliptical Gaussian
params = (x=20, y=20, fwhm=(5, 5), amp=0.1, theta=0)
P_ellip, mod_ellip = fit(gaussian, params, psf)
pairs(P_ellip)
pairs(::NamedTuple) with 5 entries: :x => 20.0183 :y => 20.01 :fwhm => (4.76786, 4.83557) :amp => 0.100245 :theta => 0.0
plot(
imshow(psf, title="Data"),
imshow(mod_ellip.(inds), title="Model"),
cbar=false,
ticks=false,
xlabel="",
ylabel="",
layout=2,
size=(600, 300)
)
Airy disk
Now with airydisk
params = (x=20, y=20, fwhm=5, amp=0.1, ratio=0.3)
P_airy, mod_airy = fit(airydisk, params, psf)
pairs(P_airy)
pairs(::NamedTuple) with 5 entries: :x => 19.9891 :y => 20.0061 :fwhm => 5.05739 :amp => 0.102935 :ratio => 0.375187
plot(
imshow(psf, title="Data"),
imshow(mod_airy.(inds), title="Model"),
cbar=false,
ticks=false,
xlabel="",
ylabel="",
layout=2,
size=(600, 300)
)
Moffat
And finally, with moffat
params = (x=20, y=20, fwhm=(5, 5), amp=0.1, theta=0, alpha=2)
P_moff, mod_moff = fit(moffat, params, psf)
pairs(P_moff)
pairs(::NamedTuple) with 6 entries: :x => 20.0359 :y => 20.0151 :fwhm => (4.96808, 4.9635) :amp => 0.116352 :theta => 0.0 :alpha => 1.46154
plot(
imshow(psf, title="Data"),
imshow(mod_moff.(inds), title="Model"),
cbar=false,
ticks=false,
xlabel="",
ylabel="",
layout=2,
size=(600, 300)
)
Changing optimization parameters
Any keyword arguments get passed on to Optim.optimize
, and you can change the algorithm used with the alg
keyword
# load Optim.jl to use the Newton method
using Optim
params = (x=20, y=20, fwhm=(5, 5), amp=0.1, theta=0, alpha=2)
P_moff, mod_moff = fit(moffat, params, psf; alg=Newton())
pairs(P_moff)
pairs(::NamedTuple) with 6 entries: :x => 20.0359 :y => 20.0151 :fwhm => (4.96808, 4.96349) :amp => 0.116352 :theta => 0.0 :alpha => 1.46154
We can also "freeze" parameters by creating a named tuple and passing it to func_kwargs
params = (;x=10, y=20, fwhm=(5, 5), amp=0.1)
func_kwargs = (;alpha=2)
P_moff2, mod_moff2 = fit(moffat, params, psf; func_kwargs)
pairs(P_moff2)
pairs(::NamedTuple) with 4 entries: :x => 20.0349 :y => 20.0144 :fwhm => (6.51242, 6.51775) :amp => 0.110515
plot(
imshow(psf, title="Data"),
imshow(mod_moff2.(inds), title="Model"),
cbar=false,
ticks=false,
xlabel="",
ylabel="",
layout=2,
size=(600, 300)
)