Why & how

SpectralFitting.jl is a package for fitting models to spectral data, similar to XSPEC, Sherpa or ISIS.

The rationale for this package is to provide a unanimous interface for different model libraries, and to leverage advancements in the computional methods that are available in Julia, including the rich statistics ecosystem, with automatic-differentiation and speed.

Longer term ambitions include

• Multi-wavelength fits
• Radiative transfer embedded into the package
• Spectral and timing fits

SpectralFitting aims to provide highly optimised and flexible fitting algorithms, along with a library of spectral models, for use in any field of Astronomy that concerns itself with spectral data.

Rewriting model calls during invocation

SpectralFitting.jl tries to optimise model invocation through source-rewriting. For compatibility with the XSPEC model library, this is achieved through aggressive pre-allocation and shuffling of output vectors. All XSPEC models output the result of their calculation through side effects into a flux array passed as an argument, and therefore each model invocation requires its own output before addition or multiplication of fluxes may occur.

Principally, combining several models together would look like this:

energy = collect(range(0.1, 20.0, 100))

flux1 = invokemodel(energy, PowerLaw())
flux2 = invokemodel(energy, PowerLaw(a=FitParam(3.0)))
flux3 = invokemodel(energy, BlackBody())
flux4 = invokemodel(energy, PhotoelectricAbsorption())

total_flux = @. flux4 * (flux1 + flux2 + flux3)
sum(total_flux)

6.392479254813619

But these operations could also be performed in a different order:

flux1 = invokemodel(energy, PowerLaw())
flux2 = invokemodel(energy, PowerLaw(a=FitParam(3.0)))
total_flux = @. flux1 + flux2

flux3 = invokemodel(energy, BlackBody())
@. total_flux = total_flux + flux3

flux4 = invokemodel(energy, PhotoelectricAbsorption())
@. total_flux = total_flux * flux4
sum(total_flux)

6.392479254813619

Doing so would allow us to only pre-allocate 2 flux arrays, instead of 4 when using the in-place variants:

fluxes = zeros(Float64, (length(energy) - 1, 2))
flux1, flux2 = eachcol(fluxes)

invokemodel!(flux1, energy, PowerLaw())
invokemodel!(flux2, energy, PowerLaw(a=FitParam(3.0)))
@. flux1 = flux1 + flux2

invokemodel!(flux2, energy, BlackBody())
@. flux1 = flux1 + flux2

invokemodel!(flux2, energy, PhotoelectricAbsorption())
@. flux1 = flux1 * flux2
sum(flux1)

6.392479254813619

It is precisely this re-writing that SpectralFitting performs via @generated functions. We can inspect the code used to generate the invocation body after defining a CompositeModel:

model = PhotoelectricAbsorption() * (
    PowerLaw() + PowerLaw(a=FitParam(3.0)) + BlackBody()
)

params = get_value.(SpectralFitting.parameter_tuple(model))
SpectralFitting.Reflection.assemble_composite_model_call(typeof(model), typeof(params))

quote
    #= /home/runner/work/SpectralFitting.jl/SpectralFitting.jl/src/reflection.jl:345 =#
    #= /home/runner/work/SpectralFitting.jl/SpectralFitting.jl/src/reflection.jl:345 =# @fastmath begin
            #= /home/runner/work/SpectralFitting.jl/SpectralFitting.jl/src/reflection.jl:346 =#
            #= /home/runner/work/SpectralFitting.jl/SpectralFitting.jl/src/reflection.jl:346 =# @inbounds let objective1 = view(objectives, :, 1), objective2 = view(objectives, :, 2), objective3 = view(objectives, :, 3)
                    #= /home/runner/work/SpectralFitting.jl/SpectralFitting.jl/src/reflection.jl:347 =#
                    var"##table#554" = getfield(getfield(model, :left), :table)
                    #= /home/runner/work/SpectralFitting.jl/SpectralFitting.jl/src/reflection.jl:348 =#
                    var"##K#547" = parameters[1]
                    var"##kT#548" = parameters[2]
                    var"##K#549" = parameters[3]
                    var"##a#550" = parameters[4]
                    var"##K#551" = parameters[5]
                    var"##a#552" = parameters[6]
                    var"##ηH#553" = parameters[7]
                    #= /home/runner/work/SpectralFitting.jl/SpectralFitting.jl/src/reflection.jl:349 =#
                    invokemodel!(objective1, domain, (BlackBody){Float64}(var"##K#547", var"##kT#548"))
                    invokemodel!(objective2, domain, (PowerLaw){Float64}(var"##K#549", var"##a#550"))
                    invokemodel!(objective3, domain, (PowerLaw){Float64}(var"##K#551", var"##a#552"))
                    #= /home/runner/work/SpectralFitting.jl/SpectralFitting.jl/src/reflection.jl:286 =# @__dot__ objective2 = objective2 + objective3
                    #= /home/runner/work/SpectralFitting.jl/SpectralFitting.jl/src/reflection.jl:286 =# @__dot__ objective1 = objective1 + objective2
                    invokemodel!(objective2, domain, (PhotoelectricAbsorption){Interpolations.Extrapolation{Float64, 1, Interpolations.GriddedInterpolation{Float64, 1, Vector{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}}}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Interpolations.Line{Nothing}}, Float64}(var"##table#554", var"##ηH#553"))
                    #= /home/runner/work/SpectralFitting.jl/SpectralFitting.jl/src/reflection.jl:286 =# @__dot__ objective1 = objective1 * objective2
                    #= /home/runner/work/SpectralFitting.jl/SpectralFitting.jl/src/reflection.jl:350 =#
                    return objective1
                end
        end
end

This generated function also takes care of some other things for us, such as unpacking parameters (optionally unpacking frozen parameters separately), and ensuring any closure are passed to invokemodel if a model needs them (e.g., SurrogateSpectralModel).

This is achieved by moving as much information as possible about the model and its construction to its type, such that all of the invocation and parameter unpacking may be inferred at compile time.