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: Gaussian
using HCIDatasets: BetaPictoris
using Plots

# convenience function for plotting
function imshow(data; kwargs...)
    xlim = extrema(axes(data, 2))
    ylim = extrema(axes(data, 1))
    heatmap(data; xlim=xlim, ylim=ylim, aspect_ratio=1, kwargs...)
end

# get a PSF from HCIDatasets.jl;
# you may be prompted to download the file
psf = BetaPictoris[:psf]

imshow(psf)

using LossFunctions

# generative model
function model(X::AbstractVector{T}) where T
    position = @view X[1:2] # x, y position
    fwhm     = @view X[3:4] # fwhm_x, fwhm_y
    amp      =       X[5]   # amplitude
    return amp * Gaussian{T}(position, fwhm)
end

# objective function
function loss(X::AbstractVector{T}, target) where T
    # cheap way to enforce positivity
    all(>(0), X) || return T(Inf)
    # get generative model
    m = model(X)
    # l2-distance loss (χ² loss) (LossFunctions.jl)
    stamp = @view m[axes(target)...]
    return value(L2DistLoss(), target, stamp, AggMode.Sum())
end

# params are [x, y, fwhm_x, fwhm_y, amp]
test_params = Float32[20, 20, 5, 5, 1]
loss(test_params, psf)

11.579905f0

The objective function can then be used with an optimization library like Optim.jl to find best-fitting parameters

using Optim

# Fit our data using test_params as a starting point
# uses Nelder-Mead optimization
res = optimize(P -> loss(P, psf), test_params)

 * Status: success

 * Candidate solution
    Final objective value:     1.026439e-02

 * Found with
    Algorithm:     Nelder-Mead

 * Convergence measures
    √(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    375
    f(x) calls:    645

# utilize automatic differentiation (AD) to enable
# advanced algorithms, like LBFGS
res_ad = optimize(P -> loss(P, psf), test_params, LBFGS(); autodiff=:forward)

 * Status: success

 * Candidate solution
    Final objective value:     1.026439e-02

 * Found with
    Algorithm:     L-BFGS

 * Convergence measures
    |x - x'|               = 9.54e-07 ≰ 0.0e+00
    |x - x'|/|x'|          = 4.76e-08 ≰ 0.0e+00
    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
    |g(x)|                 = 7.22e-08 ≰ 1.0e-08

 * Work counters
    Seconds run:   1  (vs limit Inf)
    Iterations:    12
    f(x) calls:    52
    ∇f(x) calls:   52

we can see which result has the better loss, and then use the generative model to create a model that we can use elsewhere

best_res = minimum(res) < minimum(res_ad) ? res : res_ad
best_fit_params = Optim.minimizer(best_res)

5-element Array{Float32,1}:
 20.018291
 20.009983
  4.767066
  4.8366284
  0.10024304

synth_psf = model(best_fit_params)

plot(
    imshow(psf, title="Data"),
    plot(synth_psf, axes(psf); title="Model"),
    cbar=false,
    ticks=false,
    layout=2,
    size=(600, 300)
)