Low-Level Functions

 composite!(composite::AbstractMatrix{<:Number},
            coeffs::AbstractVector{<:Number},
            models::AbstractVector{T}) where T <: AbstractMatrix{<:Number}

Updates the composite matrix in place with the linear combination of sum( coeffs .* models ); this is equation 1 in Dolphin 2002, $m_i = \sum_j \, r_j \, c_{i,j}$.

Examples

julia> C = zeros(5,5);
julia> models = [rand(size(C)...) for i in 1:5];
julia> coeffs = rand(length(models));
julia> composite!(C, coeffs, models);
julia> C ≈ sum( coeffs .* models)
true

composite!(composite::AbstractVector{<:Number},
           coeffs::AbstractVector{<:Number},
           models::AbstractMatrix{<:Number})

Updates the composite vector with the matrix-vector product of models * coeffs. This is equation 1 in Dolphin 2002, $m_i = \sum_j \, r_j \, c_{i,j}$.

Examples

julia> hist_size = (5,10);
julia> models = reduce(hcat,rand(prod(hist_size)) for i in 1:20);
julia> coeffs = rand(length(axes(models,2)));
julia> C = zeros(length(axes(models,1)));
julia> composite!(C, coeffs, models);
julia> C ≈ models * coeffs
true

Notes

While the other call signature for this function more closely mirrors the natural data structure for Hess diagrams (2D matrices for composite and each entry in models), this method operates on the same data but flattened. Thus composite becomes a vector rather than a matrix and models becomes a single matrix rather than a vector of matrices. The method StarFormationHistories.stack_models is provided to stack the models into this format. This data layout enables us to use the highly optimized LinearAlgebra.mul! function to perform the matrix-vector product which typically achieves >30% speedup relative to the more natural formulation. Additionally, as mul! will typically call to a BLAS matrix-vector product function like gemv! for our use-case, we can switch out Julia's default OpenBLAS at runtime for other BLAS libraries with Julia bindings like MKL and Apple Accelerate, enabling even greater performance improvements.

loglikelihood(composite::AbstractArray{<:Number}, data::AbstractArray{<:Number})

Returns the logarithm of the Poisson likelihood ratio given by equation 10 in Dolphin 2002,

\[\text{ln} \, \mathscr{L} = \sum_i -m_i + n_i \times \left( 1 - \text{ln} \, \left( \frac{n_i}{m_i} \right) \right)\]

with composite being the complex Hess model diagram $m_i$ (see StarFormationHistories.composite!) and data being the observed Hess diagram $n_i$.

Performance Notes

• ~18.57 μs for composite=Matrix{Float64}(undef,99,99) and data=similar(composite).
• ~20 μs for composite=Matrix{Float64}(undef,99,99) and data=Matrix{Int64}(undef,99,99).
• ~9.3 μs for composite=Matrix{Float32}(undef,99,99) and data=similar(composite).
• ~9.6 μs for composite=Matrix{Float32}(undef,99,99) and data=Matrix{Int64}(undef,99,99).

loglikelihood(coeffs::AbstractVector{<:Number},
              models::AbstractVector{T},
              data::AbstractMatrix{<:Number}) where T <: AbstractMatrix{<:Number}
loglikelihood(coeffs::AbstractVector{<:Number},
              models::AbstractMatrix{<:Number},
              data::AbstractVector{<:Number})

Returns the logarithm of the Poisson likelihood ratio, but constructs the complex Hess diagram model as sum(coeffs .* models) rather than taking composite directly as an argument. Second call signature supports the flattened formats for models and data. See the notes for the flattened call signature of StarFormationHistories.composite! for more details.

∇loglikelihood(model::AbstractArray{<:Number},
               composite::AbstractArray{<:Number},
               data::AbstractArray{<:Number})

Returns the partial derivative of the logarithm of the Poisson likelihood ratio (StarFormationHistories.loglikelihood) with respect to the coefficient $r_j$ on the provided model. If the complex Hess diagram model is $m_i = \sum_j \, r_j \, c_{i,j}$, then model is $c_{i,j}$, and this function computes the partial derivative of $\text{log} \, \mathscr{L}$ with respect to the coefficient $r_j$. This is given by equation 21 in Dolphin 2002,

\[\frac{\partial \, \text{log} \, \mathscr{L}}{\partial \, r_j} = \sum_i c_{i,j} \left( \frac{n_i}{m_i} - 1 \right)\]

where $n_i$ is bin $i$ of the observed Hess diagram data.

Performance Notes

• ~4.1 μs for model, composite, data all being Matrix{Float64}(undef,99,99).
• ~1.3 μs for model, composite, data all being Matrix{Float32}(undef,99,99).

∇loglikelihood(models::AbstractVector{T},
               composite::AbstractMatrix{<:Number},
               data::AbstractMatrix{<:Number}) where T <: AbstractMatrix{<:Number}
∇loglikelihood(models::AbstractMatrix{<:Number},
               composite::AbstractVector{<:Number},
               data::AbstractVector{<:Number})

Computes the gradient of the logarithm of the Poisson likelihood ratio with respect to the coefficients by calling the single-model ∇loglikelihood for every model in models. Second call signature supports the flattened formats for models and data. See the notes for the flattened call signature of StarFormationHistories.composite! for more details.

∇loglikelihood(coeffs::AbstractVector{<:Number},
               models::AbstractVector{T},
               data::AbstractMatrix{<:Number}) where T <: AbstractMatrix{<:Number}
∇loglikelihood(coeffs::AbstractVector{<:Number},
               models::AbstractMatrix{<:Number},
               data::AbstractVector{<:Number})

Forms the composite matrix from coefficients coeffs and model templates models and returns the gradient of the loglikelihood with respect to the coefficients. Second call signature supports the flattened formats for models and data. See the notes for the flattened call signature of StarFormationHistories.composite! for more details.

 ∇loglikelihood!(G::AbstractVector,
                 composite::AbstractMatrix{<:Number},
                 models::AbstractVector{S},
                 data::AbstractMatrix{<:Number}) where S <: AbstractMatrix{<:Number}

Efficiently computes the gradient of StarFormationHistories.loglikelihood with respect to all coefficients by updating G with the gradient. This will overwrite composite with the result of 1 .- (data ./ composite) so it shouldn't be reused after being passed to this function.

Arguments

• G::AbstractVector is the vector that will be mutated in-place with the computed gradient values.
• models::AbstractVector{<:AbstractMatrix{<:Number}} is the vector of matrices giving the model Hess diagrams.
• composite::AbstractMatrix{<:Number} is a matrix that contains the composite model sum(coeffs .* models).
• data::AbstractMatrix{<:Number} contains the observed Hess diagram that is being fit.

G = ∇loglikelihood!(G::AbstractVector,
                    composite::AbstractVector{<:Number},
                    models::AbstractMatrix{<:Number},
                    data::AbstractVector{<:Number})

Updates and returns G with the gradient of the loglikelihood with respect to all coefficients. This call signature supports the flattened formats for models and data. See the notes for the flattened call signature of StarFormationHistories.composite! for more details.

fg!(F, G, MH_model0::AbstractMetallicityModel,
    disp_model0::AbstractDispersionModel,
    variables::AbstractVector{<:Number},
    models::Union{AbstractMatrix{<:Number},
                  AbstractVector{<:AbstractMatrix{<:Number}}},
    data::Union{AbstractVector{<:Number}, AbstractMatrix{<:Number}},
    composite::Union{AbstractVector{<:Number}, AbstractMatrix{<:Number}},
    logAge::AbstractVector{<:Number},
    metallicities::AbstractVector{<:Number})

Main function that differs between AMR and MZR models that accounts for the difference in the gradient formulations between models. F and G mirror the Optim.jl interface for computing the objective and gradient in a single function to make use of common intermediate computations.

Arguments

• F controls whether the objective is being requested. If !isnothing(F), the negative log likelihood will be returned from fg!.
• G controls whether the gradient of the objective with respect to the variables is being requested. If !isnothing(G), G will be updated in-place with the gradient of the negative log likelihood with respect to the fitting parameters variables.
• MH_model0 is an instance of a concrete subtype of AbstractMetallicityModel (e.g., PowerLawMZR) that specifies the metallicity evolution model to use. The parameters contained in MH_model0 are used as initial parameters to begin the optimization in fit_sfh, but are not used internally in fg! – new instances are constructed from variables instead.
• disp_model0 is an instance of a concrete subtype of AbstractDispersionModel (e.g., GaussianDispersion) that specifies the PDF of the metallicities of stars forming at fixed time. The parameters contained in disp_model0 are used as initial parameters to begin the optimization in fit_sfh, but are not used internally in fg! – new instances are constructed from variables instead.
• variables are the fitting parameters, including free and fixed parameters. This vector is split into stellar mass coefficients R_j, metallicity model parameters, and dispersion model parameters, and so must contain all relevant fittable parameters, even those that are to be fixed during the solve.
• models are the template Hess diagrams for the SSPs that compose the observed Hess diagram.
• data is the Hess diagram for the observed data.
• composite is the pre-allocated array in which to store the complex Hess diagram model. Must have same shape as data.
• logAge::AbstractVector{<:Number} is the vector containing the effective ages of the stellar populations used to create the templates in models, in units of log10(age [yr]). For example, if a population has an age of 1 Myr, its entry in logAge should be log10(10^6) = 6.0.
• metallicities::AbstractVector{<:Number} is the vector containing the effective metallicities of the stellar populations used to create the templates in models. These should be logarithmic abundances like [M/H] or [Fe/H]. There are some notes on the Wikipedia that might be useful.

Returns

• Negative log likelihood if !isnothing(F).

keep_idx::Vector{Int} =
    truncate_relweights(relweightsmin::Number,
                        relweights::AbstractVector{<:Number},
                        logAge::AbstractVector{<:Number})

Method to truncate an isochrone grid with log10(age [yr]) values logAge and relative weights relweights derived from a hierarchical metallicity evolution model to only include template models with relweights greater than relweightsmin times the maximum relative weight for each unique entry in logAge. The input vectors are the same as those for StarFormationHistories.fixed_amr, which includes more information. Returns a vector of the indices into relweights and logAge of the template models whose relative weights are significant given the provided relweightsmin.

Examples

When using a fixed input age-metallicity relation as enabled by, for example, StarFormationHistories.fixed_amr, only the star formation rate (or total stellar mass) coefficients need to be fit, as the metallicity distribution is no longer a free parameter in the model. As such, the relative weights of each model with identical logAge but different metallicities only need to be computed once at the start of the optimization. As the metallicity distribution is not a free parameter, it is also possible to truncate the list of models to only those that contribute significantly to the final composite model to improve runtime performance. That is what this method does.

A simple isochrone grid will be two-dimensional, encompassing age and metallicity. Consider a subset of the model grid with the same age such that unique(logAge) = [10.0] but a series of different metallicities, metallicities = -2.5:0.25:0. If we model the metallicity distribution function for this age as having a mean [M/H] of -2.0 and a Gaussian spread of 0.2 dex, then the relative weights of these models can be approximated as

import Distributions: Normal, pdf
metallicities = -2.5:0.25:0
relweights = pdf.(Normal(-2.0, 0.2), metallicities)
relweights ./= sum(relweights) # Normalize the relative weights to unity sum

11-element Vector{Float64}:
 0.021919934465195145
 0.2284109622221623
 0.4988954088848224
 0.2284109622221623
 0.021919934465195145
 0.0004409368867815243
 1.8592101580561089e-6
 1.6432188478108614e-9
 3.0442281937632026e-13
 1.1821534989089337e-17
 9.622444440364979e-23

Several of these models with very low relative weights are unlikely to contribute significantly to the final composite model. We can select out only the significant ones with, say, relative weights greater than 10% of the maximum as StarFormationHistories.truncate_relweights(0.1, relweights, fill(10.0,length(metallicities))) which will return indices into relweights whose values are greater than 0.1 * maximum(relweights) = 0.04988954088848224,

3-element Vector{Int64}:
 2
 3
 4

which correspond to relweights[2,3,4] = [ 0.2284109622221623, 0.4988954088848224, 0.2284109622221623 ]. If we use only these 3 templates in the fit, instead of the original 11, we will achieve a speedup of almost 4x with a minor loss in precision which, in most cases, will be less than the numerical uncertainties on the individual star formation rate parameters. However, as fits of these sort are naturally quite fast, we recommend use of this type of truncation only in applications where many fits are required (e.g., Monte Carlo experiments). For most applications, this level of optimization is not necessary.

calculate_edges(edges, xlim, ylim, nbins, xwidth, ywidth)

Function to calculate the bin edges for 2D histograms. Returns (xbins, ybins) with both entries being ranges.

Keyword Arguments

• edges is a tuple of ranges defining the left-side edges of the bins along the x-axis (edges[1]) and the y-axis (edges[2]). Example: (-1.0:0.1:1.5, 22:0.1:27.2). If edges is provided, it will simply be returned.
• xlim is a length-2 indexable object (e.g., a Vector{Float64} or NTuple{2,Float64)) giving the lower and upper bounds on the x-axis corresponding to the provided colors array. Example: [-1.0, 1.5]. This is only used if edges==nothing.
• ylim is like xlim but for the y-axis corresponding to the provided mags array. Example [25, 20]. This is only used if edges==nothing.
• nbins::NTuple{2,<:Integer} is a 2-tuple of integers providing the number of bins to use along the x- and y-axes. This is only used if edges==nothing.
• xwidth is the bin width along the x-axis for the colors array. This is only used if edges==nothing and nbins==nothing. Example: 0.1.
• ywidth is like xwidth but for the y-axis corresponding to the provided mags array. Example: 0.1.