GitHub

Using spectral models

In this page you'll find how to use the spectral model library and how to define your own models. Using the model library is as easy as invoking or composing models from the Model Index. For example:

model = PowerLaw()
┌ PowerLaw
│      K ->  1 ± 0.1  ∈ [ 0, Inf ]   FREE
│      a ->  2 ± 0.2  ∈ [ 0, Inf ]   FREE
└

In the output of the REPL we see the model name, and it's two parameters, with information about those parameters, such as the current value, the associated error (10% by defaul), the minimum and maximum values, and whether the parameter is frozen or not.

Note

See FitParam for full details about fitable parameters.

The parameters can be tweaked by accessing the fields

model.K.value = 2.0
model.K.frozen = true
model
┌ PowerLaw
│      K ->  2                       FROZEN
│      a ->  2 ± 0.2  ∈ [ 0, Inf ]   FREE
└

We can invoke the model on a domain in the following way

domain = collect(range(0.1, 10.0, 100))
invokemodel(domain, model)
99-element Vector{Float64}:
 10.0
  3.333333333333333
  1.666666666666667
  1.0
  0.6666666666666665
  0.4761904761904763
  0.3571428571428572
  0.2777777777777777
  0.22222222222222232
  0.18181818181818188
  ⋮
  0.002388915432393668
  0.002337540906965918
  0.002287805994051678
  0.00223964165733484
  0.00219298245614033
  0.0021477663230240474
  0.002103934357248094
  0.0020614306328591847
  0.0020202020202020055
Note

By default, models are implemented to accept a single input vector with all of the low and high bin edges, and return a flux array with the flux in each energy bin. As such, it is here the case that:

length(flux) == length(energy) - 1

Models need not be defined as such, however. See AbstractLayout for more.

Models can be composed together following the Model algebra. That means to expressive a photoelectric absorption component acting on the power law we can write

model2 = PhotoelectricAbsorption() * model
┌ CompositeModel with 2 model components:
│      m1 * a1
│ Model key and parameters:
│    a1 => PowerLaw
│       K_1 ->  2                       FROZEN
│       a_1 ->  2 ± 0.2  ∈ [ 0, Inf ]   FREE
│    m1 => PhotoelectricAbsorption
│      ηH_1 ->  1 ± 0.1  ∈ [ 0, Inf ]   FREE
└

The parameters of this CompositeModel are are copied from the expression. This means we can modify the K_1 parameter in model2 without having to worry that we are changing model.K:

model2.K_1.frozen = false
model2
┌ CompositeModel with 2 model components:
│      m1 * a1
│ Model key and parameters:
│    a1 => PowerLaw
│       K_1 ->  2 ± 0.1  ∈ [ 0, Inf ]   FREE
│       a_1 ->  2 ± 0.2  ∈ [ 0, Inf ]   FREE
│    m1 => PhotoelectricAbsorption
│      ηH_1 ->  1 ± 0.1  ∈ [ 0, Inf ]   FREE
└ 
model
┌ PowerLaw
│      K ->  2                       FROZEN
│      a ->  2 ± 0.2  ∈ [ 0, Inf ]   FREE
└

Composite models have the same methods as single models. This means we can invoke a model in the same way

invokemodel(domain, model2)
99-element view(::Matrix{Float64}, :, 1) with eltype Float64:
 2.0751585013952092e-263
 4.495947705380242e-37
 6.169182651763269e-13
 4.3452153392915345e-6
 0.00071255337382327
 0.0002950219828216392
 0.0025894772776999717
 0.005476504002915039
 0.008240937581191064
 0.014887524257980244
 ⋮
 0.002352926658712353
 0.0023033682709992184
 0.002255340978829189
 0.002208782849254055
 0.002163635423211736
 0.0021198431935328848
 0.0020773528501903667
 0.0020361145193630906
 0.001996079508045848

Defining new models

To define your own model, you need to tell the package what the model parameters are and how to invoke the model. This is all done by creating a struct which subtypes AbstractSpectralModel.

Let's create a new Additive spectral model:

Base.@kwdef struct MyModel{T} <: AbstractSpectralModel{T,Additive}
    K::T = FitParam(2.0)
    p::T = FitParam(3.0)
end

# implementing a dummy add operation this function can do anything it likes, but
# must write the output into `output` and ideally should be thread safe
function SpectralFitting.invoke!(output, input, model::MyModel)
    SpectralFitting.finite_diff_kernel!(output, input) do E
        E + model.p
    end
end

Here we used the utility method SpectralFitting.finite_diff_kernel! to ensure the additive model is appropriately scaled across the bin width.

Note that Additive models do not need to use the normalization parameter K themselves. This is because when we use invokemodel these sorts of translations are automatically applied, for compatability with external models.

Our model is now ready to use

model = MyModel()
┌ MyModel
│      K ->  2 ± 0.2  ∈ [ 0, Inf ]   FREE
│      p ->  3 ± 0.3  ∈ [ 0, Inf ]   FREE
└ 
domain = collect(range(0.1, 10.0, 100))
invokemodel(domain, model)
99-element Vector{Float64}:
 0.20000000000000018
 0.1999999999999993
 0.20000000000000018
 0.20000000000000018
 0.20000000000000018
 0.20000000000000018
 0.1999999999999993
 0.20000000000000018
 0.20000000000000018
 0.1999999999999993
 ⋮
 0.1999999999999993
 0.20000000000000284
 0.1999999999999993
 0.1999999999999993
 0.1999999999999993
 0.1999999999999993
 0.20000000000000284
 0.1999999999999993
 0.1999999999999993
Note

To add new XSPEC or foreign function models, see Wrapping new XSPEC models.

Model abstraction

All spectral models are a sub-type of AbstractSpectralModel.

SpectralFitting.AbstractSpectralModel — Type
abstract type AbstractSpectralModel{T,K<:AbstractSpectralModelKind} end

Supertype of all spectral models, tracking the number type T and AbstractSpectralModelKind denoted K.

Implementation

Sub-types must implement the following interface (see the function's documentation for examples):

Usage

The available API for a spectral model is detailed below:

The following query functions exist:

Model reflection is supported by the following functions. These are intended for internal use and are not exported.

The parametric type parameter T is the number type of the model and K defines the AbstractSpectralModelKind.

source
SpectralFitting.invoke! — Function
SpectralFitting.invoke!(output, domain, M::Type{<:AbstractSpectralModel}, params...)

Used to define the behaviour of models. Should calculate the output of the model and write in-place into output. The model parameters are passed in the model structure.

Warning

This function should not be called directly. Use invokemodel instead. invoke! is only to define the model, not to use it. Users should always call models using invokemodel or invokemodel! to ensure normalisations and closures are accounted for.

Example

Base.@kwdef struct MyModel{T} <: AbstractSpectralModel{T,Multiplicative}
    p1::T = FitParam(1.0)
    p2::T = FitParam(2.0)
    p3::T = FitParam(3.0)
end

would have the arguments passed to invoke! as

function SpectralFitting.invoke!(output, domain, model::MyModel)
    # ...
end
source
SpectralFitting.modelkind — Function
modelkind(M::Type{<:AbstractSpectralModel})
modelkind(::AbstractSpectralModel)

Return the kind of model given by M: either Additive, Multiplicative, or Convolutional.

source
SpectralFitting.numbertype — Function
numbertype(::AbstractSpectralModel)

Get the numerical type of the model. This goes through FitParam, so that the number type returned is as close to a primative as possible.

See also paramtype.

Example

numbertype(PowerLaw()) == Float64
source

Model methods

SpectralFitting.invokemodel — Function
invokemodel(domain, model::AbstractSpectralModel)

Invoke the AbstractSpectralModel given by model over the domain domain.

This function will perform any normalisation or post-processing tasks that a specific model kind may require, e.g. multiplying by a normalisation constant for Additive models.

Note

invokemodel allocates the needed output arrays based on the element type of free_params to allow automatic differentation libraries to calculate parameter gradients.

In-place non-allocating variants are the invokemodel! functions.

Example

model = PowerLaw()
domain = collect(range(0.1, 20.0, 100))

invokemodel(domain, model)
source
SpectralFitting.invokemodel! — Function
invokemodel!(output, domain, model)
invokemodel!(output, domain, model, params::AbstractVector)
invokemodel!(output, domain, model, params::ParameterCache)

In-place variant of invokemodel, calculating the output of an AbstractSpectralModel given by model, optionally overriding the parameters using a ParameterCache or an AbstractVector.

The output may not necessarily be a single vector, and one should use allocate_model_output to allocate the output structure.

Example

model = PowerLaw()
domain = collect(range(0.1, 20.0, 100))
output = allocate_model_output(model, domain)
invokemodel!(output, domain, model)
source

Model algebra

Models exist as three different kinds, defined by an AbstractSpectralModelKind trait.

SpectralFitting.Additive — Type
Additive <: AbstractSpectralModelKind
Additive()

Additive models are effectively the sources of photons, and are the principle building blocks of composite models. Every additive model has a normalisation parameter which re-scales the output by a constant factor K.

Note

Defining custom additive models requires special care. See Defining new models.

source
SpectralFitting.Multiplicative — Type
Multiplicative <: AbstractSpectralModelKind
Multiplicative()

Multiplicative models act on Additive models, by element-wise multiplying the output in each domain bin of the additive model by a different factor.

source

Model data availability

Many of the XSPEC implemented models use tabular data, such as FITS, and return results interpolated from these pre-calculated tables. In some cases, these table models have data files that are multiple gigabytes in size, and would be very unwieldy to ship indiscriminantly. SpectralFitting attempts to circumnavigate this bloat by downloading the model data on an ut opus basis.

SpectralFitting.download_model_data — Function
SpectralFitting.download_model_data(model::AbstractSpectralModel; kwargs...)
SpectralFitting.download_model_data(M::Type{<:AbstractSpectralModel}; kwargs...)
SpectralFitting.download_model_data(s::Symbol; kwargs...)

Downloads the model data for a model specified either by model, type M, or symbol s. Datafiles associated with a specific model may be registered using SpectralFitting.register_model_data. The download is currently unconfigurable, but permits slight control via a number of keyword arguments:

  • progress::Bool = true

Display a progress bar for the download.

  • model_source_url::String = "http://www.star.bris.ac.uk/fbaker/XSPEC-model-data"

The source URL used to download the model data.

All standard XSPEC spectral model data is currently being hosted on the University of Bristol astrophysics servers, and should be persistently available to anyone.

source

Special care must be taken if new XSPEC models are wrapped to ensure the data is available. For more on this, see Wrapping new XSPEC models.

Model data may also alternatively be copied in by-hand from a HEASoft XSPEC source directory. In this case, the location to copy the data to may be determined via joinpath(SpectralFitting.LibXSPEC_jll.artifact_dir, "spectral", "modelData").