API reference

Must support the same API, but may also have some query methods for specific internals.

abstract type AbstractStatistic

Abstract supertype for all statistic measures.

ParameterTriple(model_index, component, parameter)
ParameterTriple(info_tuple)

Used to unambiguously denote a given parameter in a fitting problem.

abstract type AbstractTableModel{T,K} <: AbstractSpectralModel{T,K} end

Abstract type representing table models, i.e. those models that interpolate or load data from a table.

First field in the struct must be table. See PhotoelectricAbsorption for an example implementation.

Constant

• value: Constant value.

ModelPatchView{P<:ParameterPatchView}

Similar to ParameterPatchView but for accessing the parameter components by symbol. Completes the ability to modify parameters as

pv = ParameterPatchView(...)
pv.a1.K = pv.a2.K * 3 + 2

TableModelInterpolation
TableModelInterpolation(tmd::TableModelData)

Wraps a TableModelData and augments it with a multi-linear interpolation cache.

This can then be used with interpolate_table! to interpolate the table data.

abstract type AbstractSpectralModelImplementation end

Details about the implementation are represented by this abstract type, used in the trait pattern. Concrete types are

• JuliaImplementation(): means the model is implemented in the Julia programming language.
• XSPECImplementation(): means the model is vendored from the XSPEC model library.

ParameterPatchView{V<:AbstractVector}

Part of ParameterPatch.

Used to create a view on a set of parameters that can be manipulated in an intuitive way. Defines property setting / getting methods so that the parameters associated with a particular model can be reached

p = ParameterPatchView(...)
p.a1 # returns a ModelPatchView

The possible properties that can be access are defined in a symbol vector.

Every "model" accessed this way will return a view on the parameters via ModelPatchView.

TableGridData

Used to contain table grid data in the interpolation scheme.

TableModelData

A container for data read in from an XSPEC-style Xtable table model. This structure currently assumes a number of things:

• The energy bins are strictly contiguous (E_lo[i+1] = E_hi[i])
• The parameter grid is rectilinear

• energy_bins: Energy bins used for all of the grid entries.

• params: The parameter axes, with the range of tabulated values.

• grids: All grids, laid out in the same way as the parameter axes.

FittableMultiDataset(d1, d2...)

A thin wrapper representing multiple datasets.

struct BinnedData{Tag,V,D} <: AbstractDataset

function BinnedData(
    domain,
    codomain;
    tag = nothing,
    domain_variance = nothing,
    codomain_variance = nothing,
    user_data = nothing,
    name = "[no-name]",
)
    @assert length(domain) == length(codomain) + 1 "Data does not look like it is binned!"
    BinnedData{tag}(domain, codomain, domain_variance, codomain_variance, name, user_data)
end

A datastructure for contiguously binned data (see ContiguouslyBinned). This is the binned analogue to InjectiveData.

The parameters have the following meaning:

• Tag: an optional user-defined tag that can be used to overwrite method dispatches (default: nothing).
• V: the type of the domain, co-domain, and variance vectors (default: inferred from arguments).
• D: the optional additional user-defined data type (default: nothing).

• domain: The domain of the data as pairs of low and high bins, such that E_i and E_i+1 are the low and high bins of element i respectively.

• codomain: The co-domain (objective) of the dataset.

• domain_variance: The variance for each bin edge. Todo: maybe remove?

• codomain_variance: The variance for the value in each bin.

• name: Dataset name.

• user_data: Optional user specific data that can be kept along with the dataset.

FittingProblem

A struct representing a combination of models and dataset to be fit.

FittableMultiModel(m1, m2...)

A thin wrapper representing multiple models.

Base.copy(m::AbstractTableModel)

Create a copy of an AbstractTableModel. This will copy all fields except the table field, which is assumed to be a constant set of values that can be shared by multiple copies.

When this is not the case, the user should redefine Base.copy for their particular table model to copy the table as needed.

Base.copy(m::AsConvolution)

Creates a copy of an AsConvolution wrapped model. Will make a deepcopy of the cache to eliminate possible thread contention, but does not copy the domain.

_accumulated_indices(items)

items is a tuple or vector of lengths n1, n2, ...

Returns a tuple or array with same length as items, which gives the index boundaries of an array with size n1 + n2 + ....

Bind the symbols of last(pair) in all models indexes by first(pair).

Assumes A is binned in X1 and kernel is binned in X2. Output will also be binned on X1

Does the background subtraction and returns units of counts. That means we have multiplied through by a factor $t_D$ relative to the reference equation (2.3) in the XSPEC manual.

allocate_model_output(model::AbstractSpectralModel, domain::AbstractVector)

Allocate the output space for the AbstractSpectralModel for a given domain. The output type will be the promoted element type of the model and domain.

Uses construct_objective_cache to construct the appropriate layout.

apply_patch!(model::ParameterPatch)

Apply a patch to the model (i.e. use the patch! function to update the model parameters).

bind!(prob::FittingProblem, pairs[, pairs...])

Bind parameters together within a FittingProblem. Parameters bound together will be mandated to have same value during the fit.

The binding must be specified in double or triple selector, which follow the format:

pair := (model_index, :component_name, :parameter_symbol)
     := (model_index, :parameter_symbol)

The model_index is the index of the model in a multi-fit problem, i.e. 1, 2, and so on.

The component name is a CompositeModel model name, e.g. :a1, or :c3. This can be omitted if the model is not a CompositeModel.

The parameter symbol is a symbol representing the field of the parameter in the model. That is, :K or :log10Flux.

Bindings are specified using a chain of pairs (root) => (target) [=> (target)]. The root parameter is kept as is, and all subsequent parameters are bound to the root. Multiple chains of pairs may be specified in a single call to bind!, or, alternatively, multiple bindings may be specified with successive calls to bind!.

Bindings can be inspected with details.

See also bindall!.

Examples

• Bind model 1's K parameter to model 2's second additive model's K:

  bind!(prob, (1, :K) => (2, :a2, :K))

• Bind model 3's :a2.K parameter to model4's :m3.L and model 6's :a1.a:

  bind!(prob, (3, :a2, :K) => (4, :m3, :K) => (6, :a1, :a))

Consider the following two models

model1 = PhotoelectricAbsorption() * (BlackBody() + PowerLaw())
model2 = PhotoelectricAbsorption() * (PowerLaw() + PowerLaw())

prob = FittingProblem(model1 => data1, model2 => data2)

# Bind the power law indices in the two models
bindall!(prob, :a)

# Bind the normalisation of powerlaws in the 2nd model:
bind!(prob, (2, :a1, :K) => (2, :a2, :K))

# To inspect the overall bindings.
details(prob)

bindall!(prob::FittingProblem, item[, item...])

Bind a common parameter across all models. The item is used to select the parameter to bind, and may either be a single symbol, or a model-symbol double.

Examples

# bind parameter `a` in all models
bindall!(prob, :a)

# bind parameter `K` in component `a3` in all models
bindall!(prob, (:a3, :K))

# multiple simultaneously
bindall!(prob, :E, (:a2, :K))

common_support(x, y)

Find the common AbstractDataLayout of x and y, following the ordering of preferred_support.

count_error(k, σ)

Gives the error on k mean photon counts to a significance of σ standard deviations about the mean.

Derived from likelihood of binomial distributions being the beta function.

details(prob::FittingProblem)

Show details about the fitting problem, including the specific model parameters that are bound together.

SpectralFitting.download_all_model_data()

Downloads all model data for the models currently registered with SpectralFitting.register_model_data. Calls SpectralFitting.download_model_data to perform the download.

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.

error_statistic(::AbstractDataset)

Should return an ErrorStatistics describing which error statistic this data uses.

If undefined for a derived type, returns ErrorStatistics.Unknown.

finite_diff_kernel!(f::Function, flux, energy)

Calculates the finite difference of the function f over the energy bin between the high and low bin edges, via

\[c_i = f(E_{i,\text{high}}) - f(E_{i,\text{low}}),\]

similar to evaluating the limits of the integral between $E_{i,\text{high}}$ and $E_{i,\text{low}}$.

This utility function is primarily used for Additive models to ensure the flux per bin is normalised for the energy over the bin.

folded_energy(response::ResponseMatrix)

Get the contiguously binned energy corresponding to the output (folded) domain of the response matrix. That is, the channel energies as used by the spectrum.

get_invoke_wrapper(config::FittingConfig)

Creates a wrapper around the config that can be used to invoke the model by only passing the free parameters as arguments:

f = get_invoke_wrapper(config)
f(arg1, arg2, arg3, ...)

See also: get_invoke_wrapper_single.

get_invoke_wrapper_single(config::FittingConfig)

Similar to get_invoke_wrapper but calls only on the output of the model to unpack the result tuple.

get_model_data_paths(M::Type)

Load the model data using the Base.load function and return the corresponding structures.

See also: get_model_data_paths

get_model_data_paths(M::Type)

Get the absolute file paths to the local data files.

See also: get_model_data

implementation(model::AbstractSpectralModel)
implementation(::Type{<:AbstractSpectralModel})

Get the AbstractSpectralModelImplementation for a given AbstractSpectralModel or model type.

This is used primarily to learn what optimizations we can do with a model, for example propagating auto-diff gradients through a model or arbitrary precision numbers.

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.

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

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)

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.

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)

make_domain_variance(layout::AbstractDataLayout, dataset::AbstractDataset)

Make the variance vector associated with the domain.

make_label(d::AbstractDataset)

Return a string that gives a descriptive label for this dataset.

make_model_domain(layout::AbstractDataLayout, dataset::AbstractDataset)

Returns the array used as the domain for the modelling. This is paired with make_domain_variance

make_objective(layout::AbstractDataLayout, dataset::AbstractDataset)

Returns the array used as the target for model fitting. The array must correspond to the data AbstractDataLayout specified by the layout parameter.

In as far as it can be guaranteed, the memory in the returned array will not be mutated by any fitting procedures.

Domain for this objective should be returned by make_model_domain.

make_objective_variance(layout::AbstractDataLayout, dataset::AbstractDataset)

Make the variance vector associated with each objective point.

make_output_domain(layout::AbstractDataLayout, dataset::AbstractDataset)

Returns the array used as the output domain. That is, in cases where the model input and output map to different domains, the input domain is said to be the model domain, the input domain is said to be the model domain.

The distinction is mainly used for the purposes of simulating data and for visualising data.

make_surrogate_harness(
    model::M,
    lowerbounds::T,
    upperbounds::T;
    optimization_samples = 200,
    seed_samples = 50,
    S::Type = RadialBasis,
    sample_type = SobolSample(),
    verbose = false,
)

Creates and optimizes a surrogate model of type S for model, using wrap_model_as_objective and optimize_accuracy! for optimization_samples iterations. Model is initially seeded with seed_samples points prior to optimization.

modelkind(M::Type{<:AbstractSpectralModel})
modelkind(::AbstractSpectralModel)

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

objective_transformer(layout::AbstractDataLayout, dataset::AbstractDataset)

Used to transform the output of the model onto the output domain. For spectral fitting, this is the method used to do response folding and bin masking.

If none provided, uses the _DEFAULT_TRANSFORMER

function _DEFAULT_TRANSFORMER()
    function _transformer!!(domain, objective)
        objective
    end
    function _transformer!!(output, domain, objective)
        @. output = objective
    end
    _transformer!!
end

operation_symbol(::AbstractCompositeOperator)
operation_symbol(::Type{<:AbstractCompositeOperator})

Obtain the model symbol from a given AbstractCompositeOperator.

optimize_accuracy!(
    surr::AbstractSurrogate,
    obj::Function,
    lb,
    ub;
    sample_type::SamplingAlgorithm = SobolSample(),
    maxiters = 200,
    N_truth = 5000,
    verbose = false,
)

Improve accuracy (faithfullness) of the surrogate model in recreating the objective function.

Samples a new space of N_truth points between lb and ub, and calculates the objective function obj at each. Finds the point with largest MSE between surrogate and objective, and adds the point to the surrogate pool. Repeats maxiters times, adding maxiters points to surrogate model.

Optionally print to stdout the MSE and iteration count with verbose = true.

Note that upper- and lower-bounds should be in the form (E, params...), where E is a single energy and params are the model parameters.

parameter_count(m::AbstractSpectralModel)

Return the total number of parameters of the model (i.e. how many fields are FitParam.

preferred_support(x)

Get the preferred AbstractDataLayout of x. If multiple supports are available, the DEFAULT_SUPPORT_ORDERING is followed:

DEFAULT_SUPPORT_ORDERING = (ContiguouslyBinned{Nothing}(nothing), OneToOne{Nothing}(nothing))

preferred_units(::Type{<:AbstractDataset}, s::AbstractStatistic)
preferred_units(x, s::AbstractStatistic)

Get the preferred units that a given dataset would use to fit the AbstractStatistic in. For example, for ChiSquared, the units of the model may be a rate, however for Cash the preferred units might be counts.

Returning nothing from this function implies there is no unit preference.

If undefined for a derived type, returns nothing.

See also support_units.

SpectralFitting.register_model_data(M::Type{<:AbstractSpectralModel}, model_data::ModelDataInfo...)
SpectralFitting.register_model_data(M::Type{<:AbstractSpectralModel}, remote_and_local::Tuple{String,String}...)
SpectralFitting.register_model_data(M::Type{<:AbstractSpectralModel}, filenames::String...)
SpectralFitting.register_model_data(s::Symbol, filenames::String...)

Register filenames as model data associated with the model given by type M or symbol s. This function does not download any files, but rather adds the relevant filenames to a lookup which SpectralFitting.download_model_data consults when invoked, and consequently model data is only downloaded when needed.

Example

# by type
register_model_data(XS_Laor, "ari.mod")
# by symbol
register_model_data(:XS_KyrLine, "KBHline01.fits")

remake_with_parameters(model::AbstractSpectralModel, parameters)

Rebuild the model with all FitParam type parameters replaced with parameters, in the same order as they appear in the fields of the struct.

response_energy(response::ResponseMatrix)

Get the contiguously binned energy corresponding to the input domain of the response matrix. This is equivalent to the model domain.

simulate(model::AbstractSpectralModel; kwargs...)

Returns an InjectiveDataset with a realisation of the model. Can also add noise to the objective, but not to the domain.

The kwargs are:

• seed: if not nothing used to set the PRNG seed.
• var: variance of the data (assuming mean of 0)

support_units(x)

Return the units of a particular layout. If this method returns nothing, assume the layout does not care about the units and handle that information appropriately (throw an error or set defaults).

supports(x::Type)

Used to define whether a given type has support for a specific AbstractDataLayout. Should return a tuple of the supported layouts. This method should be implemented to express new support, not the query method.

To query, there is

supports(layout::AbstractDataLayout, x)::Bool

Example

supports(::Type{typeof(x)}) = (OneToOne(),)
@assert supports(ContiguouslyBinned(), x) == false

unfold(data::SpectralData)

Unfold the dataset using the response and ancillary, as in e.g. ISIS.

Returns a new Spectrum that has had the domain masking applied.

with_units(::AbstractDataLayout, units)

Remake the AbstractDataLayout with the desired units. This may be a no-op if the layout does not care about units, see support_units.

wrap_model_as_objective(model::AbstractSpectralModel; ΔE = 1e-1)
wrap_model_as_objective(M::Type{<:AbstractSpectralModel}; ΔE = 1e-1)

Wrap a spectral model into an objective function for building/optimizing a surrogate model. Returns an anonymous function taking the tuple (E, params...) as the argument, and returning a single flux value.

• https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSmodelDiskbb.html

XS_NeutralHydrogenAbsorption(nH)

• nH: Neutral hydrogen column density (units of 10²² atoms per cm⁻²)

Example

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

            XS_NeutralHydrogenAbsorption
     ┌────────────────────────────────────────┐
   1 │       ...''''''''''''''''''''''''''''''│
     │      .'                                │
     │     :                                  │
     │    :'                                  │
     │    :                                   │
     │   :                                    │
     │   :                                    │
     │   :                                    │
     │  :'                                    │
     │  :                                     │
     │  :                                     │
     │  :                                     │
     │  :                                     │
     │ :                                      │
   0 │.:                                      │
     └────────────────────────────────────────┘
      0                                     20
                       E (keV)

XS_Optxagnf()

• K: Normalisation must be frozen.

• mass: Black hole mass in solar masses.

• Dco: Comoving (proper) distance in Mpc.

• logLoLEdd: The Eddington ratio. log(L/L_Edd).

• astar: The dimensionless black hole spin.

• rcor: The coronal radius in Rg = GM/c^2.

• logrout: The log of the outer radius of the disk in units of Rg.

• kT_e: The electron temperature for the soft Comptonisation component (soft excess), in keV.

• τ: The optical depth of the soft Comptonisation component.

• Gamma: The spectral index of the hard Comptonisation component (‘power law’) which has temperature fixed to 100 keV.

• fpl: the fraction of the power below rcor which is emitted in the hard comptonisation component.

• z: Redshift.

Example

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

                                        XS_Optxagnf
                        ┌────────────────────────────────────────┐
                      0 │                                        │
                        │:                                       │
                        │:                                       │
                        │:                                       │
                        │'.                                      │
                        │ :                                      │
                        │ '.                                     │
   Flux (log scale)     │  :.                                    │
                        │   :                                    │
                        │    '.                                  │
                        │     ''':.....                          │
                        │              '''''''':.............    │
                        │                                    ''''│
                        │                                        │
                    -30 │                                        │
                        └────────────────────────────────────────┘
                         0                                     20
                                          E (keV)

XS_Jet(K, mass, Dco, log_mdot, thetaobs, BulkG, phi, zdiss, B, logPrel, gmin_inj, gbreak, gmax, s1, s2, z)

• K: MUST BE FIXED AT UNITY as the jet spectrum normalisation is set by the relativisitic particle power.

• mass: Black hole mass in solar masses.

• Dco: Comoving (proper) distance in Mpc.

• log_mdot: log(L/L_Edd.

• thetaobs: Inclination angle (deg).

• BulkG: Bulk lorentz factor of the jet.

• phi: Angular size scale (radians) of the jet acceleration region as seen from the black hole.

• zdiss: Vertical distance from the black hole of the jet dissipation region (r_g).

• B: Magnetic field in the jet (Gauss).

• logPrel: Log of the power injected in relativisitic particles (ergs/s).

• gmin_inj: Minimum lorentz factor of the injected electrons.

• gbreak: Lorentz factor of the break in injected electron distribution.

• gmax: Maximum lorentz factor.

• s1: Injected index of the electron distribution below the break.

• s2: Injected index of the electron distribution above the break.

• z: Cosmological redshift corresponding to the comoving distance used above.

Example

using UnicodePlots
energy = 10 .^collect(range(-8.0, 8.0, 100))
m = invokemodel(energy, XS_Jet())
lineplot(energy[1:end-1],m,xscale=:log10,yscale=:log10,xlim=(1e-8,1e8),ylim=(1e-8,1e8),xlabel="Energy (keV)",ylabel="Flux",title="XS_Jet",canvas=DotCanvas)

                        XS_Jet
        ┌────────────────────────────────────────┐
10⁸     │                                        │
        │                                        │
        │                                        │
        │                                        │
        │  :':..                                 │
        │ .'   ''.                               │
        │:        ':.                            │
Flux    │'          ':                           │
        │             '.                         │
        │              :                         │
        │               ''.....                  │
        │                     '''''..            │
        │                            ':..        │
        │                               ':.      │
10⁻⁸    │                                 ''.    │
        └────────────────────────────────────────┘
        10⁻⁸                                 10⁸
                        Energy (keV)

XS_CutOffPowerLaw(K, Γ, Ecut, z)

• K: Normalisation.

• Γ: Photon index.

• Ecut: Cut-off energy (keV).

• z: Redshift.

Example

using SpectralFitting
using UnicodePlots
energy = 10 .^collect(range(-1.0, 2.0, 100))
m = invokemodel(energy, XS_CutOffPowerLaw())
lineplot(energy[1:end-1],m,xscale=:log10,yscale=:log10,xlim=(1e-1,1e2),ylim=(1e-6,1e0),xlabel="Energy (keV)",ylabel="Flux",title="XS_CutOffPowerLaw",canvas=DotCanvas)

                    XS_CutOffPowerLaw
        ┌────────────────────────────────────────┐
10⁰     │:..                                     │
        │   '':...                               │
        │        '':..                           │
        │             '''..                      │
        │                  '':..                 │
        │                      '':.              │
        │                          ':.           │
Flux    │                             '..        │
        │                               ':.      │
        │                                 '.     │
        │                                   :    │
        │                                    :.  │
        │                                     :  │
        │                                      : │
10⁻⁶    │                                       :│
        └────────────────────────────────────────┘
         10⁻¹                                 10²
                        Energy (keV)

XS_Kerrconv()

• Index1: The emissivity index for the inner disk.

• Index2: The emissivity index for the outer disk.

• r_br_g: The break radius separating the inner and outer portions of the disk, in gravitational radii.

• a: The dimensionless spin parameter of the black hole.

• Incl: The disk inclination angle, in degrees. A face-on disk has Incl=0.

• Rin_ms: The inner radius of the disk, in units of the radius of marginal stability.

• Rout_ms: The outer radius of the disk, in units of the radius of marginal stability.

function _safe_ffi_invoke!(output, input, params, ModelType::Type{<:AbstractSpectralModel})

Wrapper to do a foreign function call to an XSPEC model.

function _unsafe_ffi_invoke!(
    output,
    error_vec,
    input,
    params,
    ModelType::Type{<:AbstractSpectralModel},
)

Wrapper to do a foreign function call to an XSPEC model.

See also _safe_ffi_invoke!.