Merge pull request #130 from fjebaker/fergus/doc-binding

Document parameter binding

docs/make.jl

@@ -25,12 +25,14 @@ makedocs(
             "Composite models" => "models/composite-models.md",
             "Surrogate models" => "models/surrogate-models.md",
         ],
-        # "Parameters" => "parameters.md",
+        "Fitting" => [
+            "parameters.md",
+            # "fitting.md",
+        ],
         "Datasets" => [
             "Using datasets" => "datasets/datasets.md",
             "Mission support" => "datasets/mission-support.md",
         ],
-        # "Fitting" => "fitting.md",
         "Why & How" => "why-and-how.md",
         "Reference" => "reference.md",
     ],

docs/src/parameters.md

@@ -1 +1,72 @@
-# Parameters
+# Parameters
+
+## 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`](@ref) 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!`](@ref) function that applies to your [`FittingProblem`](@ref).
+
+```@docs
+bind!
+```
+
+!!! note
+    Bindings are treated not as specific to the model but specific to the [`FittingProblem`](@ref). 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
+    ```julia
+    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.
+```julia
+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(1.22))
+
+data1 = simulate(energy, model1, var = 1e-3)
+data2 = simulate(energy, model2, var = 1e-3)
+
+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.
+```julia
+model = PowerLaw()
+prob = FittingProblem(model => data1, model => data2)
+
+bind!(prob, :a)
+prob
+```
+
+We can get a better look at our model configuration by using the [`details`](@ref) method:
+```julia
+details(prob)
+```
+
+In this printout we see that the `a` parameter of `Model 2` is bound to the `a` parameter of `Model 1`.
+
+```julia
+result = SpectralFitting.fit(prob, LevenbergMarquadt())
+
+plot(data1, xscale = :log10, yscale = :log10)
+plot!(data2, xscale = :log10, yscale = :log10)
+plot!(result[1])
+plot!(result[2])
+```
+
+Note that these fits are not perfect, because the underlying data have subtly 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 better, but the models will be entirely independent.
+
+```julia
+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])
+```

src/SpectralFitting.jl

@@ -71,6 +71,7 @@ include("datasets/datasets.jl")
 include("datasets/binning.jl")
 include("datasets/grouping.jl")
 include("datasets/injectivedata.jl")
+include("datasets/binneddata.jl")
 
 include("model-data-io.jl")
 

src/abstract-models.jl

@@ -328,7 +328,7 @@ end
 
 # printing
 
-function _printinfo(io::IO, m::M) where {M<:AbstractSpectralModel}
+function _printinfo(io::IO, m::M; bindings = nothing) where {M<:AbstractSpectralModel}
     param_tuple = parameter_named_tuple(m)
     params = [String(s) => p for (s, p) in zip(keys(param_tuple), param_tuple)]
     basename = Base.typename(M).name
@@ -343,9 +343,14 @@ function _printinfo(io::IO, m::M) where {M<:AbstractSpectralModel}
         length("$s")
     end
 
-    for (s, param) in params
+    for (i, (s, param)) in enumerate(params)
         free = param isa FitParam ? !isfrozen(param) : true
-        _print_param(io, free, s, param, param_offset, q1, q2, q3, q4)
+        val, binding = if !isnothing(bindings) && !isempty(bindings)
+            get(bindings, i, param => nothing)
+        else
+            param, nothing
+        end
+        _print_param(io, free, s, val, param_offset, q1, q2, q3, q4; binding)
     end
 end
 

src/datasets/binneddata.jl

@@ -0,0 +1,78 @@
+struct BinnedData{V} <: AbstractDataset
+    domain::V
+    codomain::V
+    domain_variance::Union{Nothing,V}
+    codomain_variance::Union{Nothing,V}
+    name::String
+end
+
+function BinnedData(
+    domain,
+    codomain;
+    domain_variance = nothing,
+    codomain_variance = nothing,
+    name = "[no-name]",
+)
+    @assert length(domain) == length(codomain) + 1 "Data does not look like it is binned!"
+    BinnedData(domain, codomain, domain_variance, codomain_variance, name)
+end
+
+supports(::Type{<:BinnedData}) = (ContiguouslyBinned(), OneToOne())
+
+bin_widths(dataset::BinnedData) = diff(dataset.domain)
+objective_units(::BinnedData) = u"counts / (s * keV)"
+spectrum_energy(dataset::BinnedData) =
+    @views (dataset.domain[1:end-1] .+ dataset.domain[2:end]) ./ 2
+
+make_model_domain(::ContiguouslyBinned, dataset::BinnedData) = dataset.domain
+make_objective(::ContiguouslyBinned, dataset::BinnedData) = dataset.codomain
+
+make_model_domain(::OneToOne, dataset::BinnedData) = dataset.domain[1:end-1]
+make_objective(::OneToOne, dataset::BinnedData) = dataset.codomain
+
+make_output_domain(layout::AbstractLayout, dataset::BinnedData) =
+    make_model_domain(layout, dataset)
+
+function make_objective_variance(::AbstractLayout, dataset::BinnedData{V})::V where {V}
+    if !isnothing(dataset.codomain_variance)
+        dataset.codomain_variance
+    else
+        # TODO: i dunno just something
+        ones(eltype(V), length(dataset.codomain))
+    end
+end
+
+function objective_transformer(::AbstractLayout, dataset::BinnedData)
+    function _transformer!!(domain, objective)
+        @views objective
+    end
+    function _transformer!!(output, domain, objective)
+        @. output = objective
+    end
+    _transformer!!
+end
+
+make_label(dataset::BinnedData) = dataset.name
+
+function _printinfo(io::IO, data::BinnedData)
+    println(
+        io,
+        Crayons.Crayon(foreground = :cyan),
+        "BinnedData",
+        Crayons.Crayon(reset = true),
+        " with ",
+        Crayons.Crayon(foreground = :cyan),
+        length(data.domain),
+        Crayons.Crayon(reset = true),
+        " data points:",
+    )
+    dmin, dmax = prettyfloat.(extrema(data.domain))
+    comin, comax = prettyfloat.(extrema(data.codomain))
+    descr = """  Name                  : $(data.name)
+      . Domain (min/max)    : ($dmin, $dmax)
+      . Codomain (min/max)  : ($comin, $comax)
+    """
+    print(io, descr)
+end
+
+export BinnedData

src/datasets/injectivedata.jl

@@ -1,4 +1,3 @@
-
 struct InjectiveData{V} <: AbstractDataset
     domain::V
     codomain::V

src/fitting/binding.jl

@@ -66,6 +66,56 @@ function _bind_pairs!(prob::FittingProblem, pairs::Vararg{<:Pair{Int,Symbol}})
     push!(prob.bindings, binding)
 end
 
+"""
+    bind!(prob::FittingProblem, pairs::Pair{Int,Symbol}...)
+    bind!(prob::FittingProblem, symbols::Symbol...)
+
+Bind parameters together within a [`FittingProblem`](@ref). Parameters bound
+together will be mandated to have exact same value during the fit.
+
+The binding may either be a single symbol that is present in all models in the
+fitting problem, or a series of pairs `Int => Symbol` which index the specific
+model and parameters to bind together. All bindings specified in a single call
+to `bind!` will be bound together. Multiple bindings are possible with repeated
+call to `bind!`.
+
+- Bind model 1's `K_1` parameter to model 2's `K_3`:
+
+```julia
+bind!(prob, 1 => :K_1, 2 => :K_3)
+```
+
+- Bind model 3's `K_2` parameter to model4's `:L_1` and model 6's `a_3`:
+
+```julia
+bind!(prob, 3 => :K_2, 4 => :L_1, 6 => :a_3)
+```
+
+- Bind the `K_1` parameter across all the models:
+
+```julia
+bind!(prob, :K_1)
+```
+
+## Examples
+
+Consider the following two models
+```julia
+model1 = PhotoelectricAbsorption() * (BlackBody() + PowerLaw())
+model2 = PhotoelectricAbsorption() * (PowerLaw() + PowerLaw())
+
+prob = FittingProblem(model1 => data1, model2 => data2)
+
+# bind the power law indices in the two models
+bind!(prob, :a_1)
+
+# bind the normalisation of powerlaws in the 2nd model:
+bind!(prob, 2 => :K_1, 2 => :K_2)
+```
+
+!!! note
+    Only free parameters can be bound together.
+"""
 bind!(prob::FittingProblem, pairs::Vararg{<:Pair}) = _bind_pairs!(prob, pairs...)
 
 function bind!(prob::FittingProblem, symbs::Vararg{Symbol})

src/plots-recipes.jl

@@ -32,6 +32,7 @@ end
     yerr --> _yerr
     _xerr = SpectralFitting.bin_widths(dataset) ./ 2
     xerr --> _xerr
+    yscale --> :identity
     markerstrokecolor --> :auto
     xlabel --> "Energy (keV)"
     ylabel --> SpectralFitting.objective_units(dataset)

src/simulate.jl

@@ -143,5 +143,35 @@ function simulate(model::AbstractSpectralModel, dataset::AbstractDataset; kwargs
     simulate!(conf; kw...)
 end
 
+"""
+    simulate(model::AbstractSpectralModel; kwargs...)
+
+Returns an [`InjectiveDataset`](@ref) 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)
+"""
+function simulate(
+    domain::AbstractVector,
+    model::AbstractSpectralModel;
+    seed::Union{Nothing,Int} = nothing,
+    simulate_distribution = Distributions.Normal,
+    var = 0.1,
+)
+    _seed::Int = isnothing(seed) ? time_ns() : seed
+    rng = Random.default_rng()
+    Random.seed!(rng, _seed)
+
+    flux = invokemodel(domain, model)
+
+    realisation = map(flux) do f
+        rand(rng, simulate_distribution(f, sqrt(var)))
+    end
+    variances = fill(var, size(realisation))
+    BinnedData(domain, realisation; codomain_variance = variances)
+end
+
 
 export simulate