Parameters
One of the core interface that SpectralFitting.jl provides are the model abstractions, and a method for controlling their parameters. Models are defined with a generic field type, so that they can be converted to primative types (e.g. Float64 or Float32) during evaluation. But whilst defining a model or a problem, the field type of all parameters is FitParam
.
SpectralFitting.FitParam
— TypeFitParam(
value;
# default error is 10%
error = 0.1 * value,
frozen = false,
lower_limit = 0,
upper_limit = Inf
)
A parameter of an AbstractSpectralModel
. The fittable parameter can be used to control how a model will be invoked during fit, e.g. what rangers a parameter can take, whether they should be fronze, and so forth.
Be aware, not all fitting algorithms support "boxed" parameters, where the lower and upper limits are taken into account.
To link parameters together, see bind!
or ParameterPatch
.
Parameters can be modified in a model in the conventional way:
model = PowerLaw()
# set the value
model.K = 2.0
# set the error
model.K.error = 0.3
# set the lower bound
model.K.lower_limit = 1.0
Accessing the values in a generic function should use get_value
.
To iterate over all parameters of a model, use parameter_vector
.
Fields
value
: The value (i.e. mean) of the parameter.error
: The absolute (±, or standard deviation) error on the parameter.lower_limit
: The lower limit the parameter can take.upper_limit
: The upper limit the parameter can take.frozen
: Is this a frozen parameter during fits?patched
: Is this parameter patched in some way?
A FitParam
is a fittable (or frozen) parameter of an AbstractSpectralModel
. Every model will instantiate with FitParam
as their parameter types so that they may be modified, bounded, frozen, or otherwise, as desired before a fit is attempted.
The following methods should be preferred over direct field access for a FitParam
:
Parameter binding
When performing a fit, it is desireable to bind certain parameters together. This ensures that they will have the same value; for example, if you were fitting two simultaneous datasets with two PowerLaw
models, you may want to have different normalisations of the model components, but enforce the power law index to be the same. To achieve this, SpectralFitting has the bind!
function that applies to your FittingProblem
.
SpectralFitting.bind!
— Functionbind!(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 paramter 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 paramters 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'sK
: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)
SpectralFitting.bindall!
— Functionbindall!(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))
Bindings are treated not as specific to the model but specific to the FittingProblem
. This is because you may want to use the same model for multiple different datasets, and have slightly different binding requirements for each one (e.g. depending on the instruments you are using). If you do need the same binding applied to two different problems, you can do that with
append!(prob1.bindings, prob2.bindings)
Caution however, this will only make sense if you are using precisely the same model in both problems.
Let's try it out. We'll generate some arbitrary powerlaw spectra with different normalisations and fit them simultaneously.
using SpectralFitting, Plots
energy = collect(range(0.1, 10.0, 100))
# two different models with different normalisations
model1 = PowerLaw(K = FitParam(100.0), a = FitParam(1.2))
model2 = PowerLaw(K = FitParam(300.0), a = FitParam(2.0))
data1 = simulate(energy, model1, var = 1e-3, seed = 42)
data2 = simulate(energy, model2, var = 1e-3, seed = 42)
plot(data1, xscale = :log10, yscale = :log10)
plot!(data2, xscale = :log10, yscale = :log10)
Now we want to fit a single powerlaw model to both of these spectra simultaneously, but with the powerlaw index fixed to be the same in both models.
model = PowerLaw()
prob = FittingProblem(model => data1, model => data2)
bindall!(prob, :a)
prob
┌ FittingProblem:
│ . Models : 2
│ . Datasets : 2
│ Parameter Summary:
│ . Total : 4
│ . Frozen : 0
│ . Bound : 1
│ . Free : 3
└
We can get a better look at our model configuration by using the details
method:
details(prob)
┌ Models:
│ Model 1: PowerLaw
│ K -> 1 ± 0.1 ∈ [ 0, Inf ] FREE
│ a -> 2 ± 0.2 ∈ [ 0, Inf ] FREE
│
│ Model 2: PowerLaw
│ K -> 1 ± 0.1 ∈ [ 0, Inf ] FREE
│ a -> Bound to Model 1 -> a
└
In this printout we see that the a
parameter of Model 2
is bound to the a
parameter of Model 1
.
result = fit(prob, LevenbergMarquadt())
plot(data1, xscale = :log10, yscale = :log10)
plot!(data2, xscale = :log10, yscale = :log10)
plot!(result[1])
plot!(result[2])
Note that this fit is bad, because the underlying data have different power law indices, but our fit is required to enforce the models to have the same value. If we release this requirement, the fit will be much better, but the models will be entirely independent.
prob = FittingProblem(model => data1, model => data2)
result = SpectralFitting.fit(prob, LevenbergMarquadt())
plot(data1, xscale = :log10, yscale = :log10)
plot!(data2, xscale = :log10, yscale = :log10)
plot!(result[1])
plot!(result[2])
Parameter patching
The parameter patching interface is still experimental and will likely change in subsesquent versions of SpectralFitting.jl
Sometimes simply linking values is not sufficient, and you need to express a complex relationship between parameters. This is where ParameterPatch
, a type of AbstractModelWrapper
is useful.
Consider the following
using SpectralFitting, Plots
real_model = GaussianLine(K = FitParam(2.0)) + GaussianLine(μ = FitParam(2.0), K = FitParam(4.0))
energy = collect(range(0.1, 10.0, 100))
data = simulate(energy, real_model; var = 4e-4, seed = 42)
plot(data)
We can fit this simple dataset quite easily, but suppose we wanted to constrain the solution to have the normalisation of one model component be exactly twice that of the other. To achieve this, we can use a patch:
function my_patch!(p)
p.a2.K = p.a1.K * 2
end
# reset the parameters so the fit starts "fresh"
model = GaussianLine() + GaussianLine(μ = FitParam(3.0))
# Wrap the model with a parameter patch
patched_model = ParameterPatch(model; patch = my_patch!)
# be sure to freeze any parameter you are planning to overwrite in a patch
patched_model.a2.K.frozen = true
patched_model
┌ ParameterPatch with 2 model components:
│ a1 + a2
│ Model key and parameters:
│ a1 => GaussianLine
│ K -> 1 ± 0.1 ∈ [ 0, Inf ] FREE
│ μ -> 6.4 ± 0.6 ∈ [ 0, Inf ] FREE
│ σ -> 1 ± 0.1 ∈ [ 0, Inf ] FREE
│ a2 => GaussianLine
│ K -> 1 FROZEN/PATCHED
│ μ -> 3 ± 0.3 ∈ [ 0, Inf ] FREE
│ σ -> 1 ± 0.1 ∈ [ 0, Inf ] FREE
└
This can then be fit as usual:
prob = FittingProblem(patched_model => data)
result = fit(prob, LevenbergMarquadt())
plot!(result)
To apply the result with a parameter patch back on a model, use update_model!
or apply_patch!
update_model!(patched_model, result)
patched_model
┌ ParameterPatch with 2 model components:
│ a1 + a2
│ Model key and parameters:
│ a1 => GaussianLine
│ K -> 1.97 ± 0.1 ∈ [ 0, Inf ] FREE
│ μ -> 6.34 ± 0.6 ∈ [ 0, Inf ] FREE
│ σ -> 1.05 ± 0.1 ∈ [ 0, Inf ] FREE
│ a2 => GaussianLine
│ K -> 3.93 FROZEN/PATCHED
│ μ -> 2.01 ± 0.3 ∈ [ 0, Inf ] FREE
│ σ -> 0.961 ± 0.1 ∈ [ 0, Inf ] FREE
└