diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json
index b55a97b..bcf72b6 100644
--- a/dev/.documenter-siteinfo.json
+++ b/dev/.documenter-siteinfo.json
@@ -1 +1 @@
-{"documenter":{"julia_version":"1.10.3","generation_timestamp":"2024-05-20T10:36:38","documenter_version":"1.4.1"}}
\ No newline at end of file
+{"documenter":{"julia_version":"1.10.3","generation_timestamp":"2024-05-20T11:24:33","documenter_version":"1.4.1"}}
\ No newline at end of file
diff --git a/dev/binaries/index.html b/dev/binaries/index.html
index dc4f65a..7d972ed 100644
--- a/dev/binaries/index.html
+++ b/dev/binaries/index.html
@@ -1,2 +1,2 @@
-Multi-Star Systems · StarFormationHistories.jl
StarFormationHistories.AbstractBinaryModel is the abstract supertype for all types that are used to model multi-star systems in the package. All concrete subtypes must implement the StarFormationHistories.sample_system method and the Base.length method, which should return an integer indicating the number of stars per system that can be sampled by the model; this is equivalent to the length of the mass vector returned by sample_system.
Simulates the effects of non-interacting, unresolved stellar companions on stellar photometry. Implementation depends on the choice of binarymodel.
Arguments
imf: an object implementing rand(imf) to draw a random mass for a single star or a stellar system (depends on choice of binarymodel)
rng::AbstractRNG: the random number generator to use when sampling stars
binarymodel::StarFormationHistories.AbstractBinaryModel: an instance of a binary model that determines which implementation will be used; currently provided options are NoBinaries, RandomBinaryPairs, and BinaryMassRatio.
Returns
masses::SVector{N,eltype(imf)}: the masses of the individual stars sampled in the system in descending order where N is the maximum number of stars that can be sampled by the provided binarymodel as given by Base.length(binarymodel).
The RandomBinaryPairs type takes one argument 0 <= fraction::Real <= 1 that denotes the number fraction of binaries (e.g., 0.3 for 30% binary fraction) and will sample binaries as random pairs of two stars drawn from the same single-star IMF. This model will ONLY generate up to one additional star – it will not generate any 3+ star systems. This model typically incurs a 10–20% speed penalty relative to NoBinaries.
The BinaryMassRatio type takes two arguments; the binary fraction 0 <= fraction::Real <= 1 and a continuous univariate distribution qdist from which to sample binary mass ratios, defined as the ratio of the secondary mass to the primary mass: $q = \text{M}_S / \text{M}_P$. The provided qdist must have the proper support of (minimum(qdist) >= 0) & (maximum(qdist) <= 1); users may find the Distributions.truncated method useful for enforcing this support on more general distributions. The default qdist is a uniform distribution from 0.1 to 1, which appears to give reasonably good agreement to observations (see, e.g., Goodwin 2013).
StarFormationHistories.AbstractBinaryModel is the abstract supertype for all types that are used to model multi-star systems in the package. All concrete subtypes must implement the StarFormationHistories.sample_system method and the Base.length method, which should return an integer indicating the number of stars per system that can be sampled by the model; this is equivalent to the length of the mass vector returned by sample_system.
Simulates the effects of non-interacting, unresolved stellar companions on stellar photometry. Implementation depends on the choice of binarymodel.
Arguments
imf: an object implementing rand(imf) to draw a random mass for a single star or a stellar system (depends on choice of binarymodel)
rng::AbstractRNG: the random number generator to use when sampling stars
binarymodel::StarFormationHistories.AbstractBinaryModel: an instance of a binary model that determines which implementation will be used; currently provided options are NoBinaries, RandomBinaryPairs, and BinaryMassRatio.
Returns
masses::SVector{N,eltype(imf)}: the masses of the individual stars sampled in the system in descending order where N is the maximum number of stars that can be sampled by the provided binarymodel as given by Base.length(binarymodel).
The RandomBinaryPairs type takes one argument 0 <= fraction::Real <= 1 that denotes the number fraction of binaries (e.g., 0.3 for 30% binary fraction) and will sample binaries as random pairs of two stars drawn from the same single-star IMF. This model will ONLY generate up to one additional star – it will not generate any 3+ star systems. This model typically incurs a 10–20% speed penalty relative to NoBinaries.
The BinaryMassRatio type takes two arguments; the binary fraction 0 <= fraction::Real <= 1 and a continuous univariate distribution qdist from which to sample binary mass ratios, defined as the ratio of the secondary mass to the primary mass: $q = \text{M}_S / \text{M}_P$. The provided qdist must have the proper support of (minimum(qdist) >= 0) & (maximum(qdist) <= 1); users may find the Distributions.truncated method useful for enforcing this support on more general distributions. The default qdist is a uniform distribution from 0.1 to 1, which appears to give reasonably good agreement to observations (see, e.g., Goodwin 2013).
We have constructed an example Jupyter notebook that highlights some common use cases supported by the package. It is recommended that you read some of the background in our main documentation before or concurrently with the example. The notebook is available in the repository as examples/fitting1.ipynb and can be viewed in rendered form at this link. It relies on a custom isochrone file which can be made available upon request to the package author.
Settings
This document was generated with Documenter.jl version 1.4.1 on Monday 20 May 2024. Using Julia version 1.10.3.
We have constructed an example Jupyter notebook that highlights some common use cases supported by the package. It is recommended that you read some of the background in our main documentation before or concurrently with the example. The notebook is available in the repository as examples/fitting1.ipynb and can be viewed in rendered form at this link. It relies on a custom isochrone file which can be made available upon request to the package author.
Settings
This document was generated with Documenter.jl version 1.4.1 on Monday 20 May 2024. Using Julia version 1.10.3.
diff --git a/dev/fitting/fitting_intro/index.html b/dev/fitting/fitting_intro/index.html
index 1cb3e41..7815c97 100644
--- a/dev/fitting/fitting_intro/index.html
+++ b/dev/fitting/fitting_intro/index.html
@@ -1,7 +1,7 @@
Background and Template Construction · StarFormationHistories.jl
In the classic formulation of star formation history fitting from resolved-star photometry (Dolphin 2002), an observed color-magnitude diagram (CMD) is binned into a 2-D histogram known as a Hess diagram. Such a CMD and Hess diagram pair is shown below.
The representation of the observations as a Hess diagram allows one to apply Poisson statistics, specifically the Poisson likelihood ratio (Equations 7–10 in Dolphin 2002), to model the observations. As the CMD of a complex stellar population is simply the sum of the CMDs of its sub-populations, one need only prepare a number of templates for each simple stellar population (SSP) which may make up the complex population in question and model the observed Hess diagram as a linear combination of these templates. Keeping the same notation as Dolphin 2002 (Equation 1), the complex model Hess diagram is simply
\[m_i = \sum_j \, r_j \, c_{i,j}\]
where $m_i$ is the value of the complex model in bin $i$, $c_{i,j}$ is the value of simple template $j$ in bin $i$, and $r_j$ is the multiplicative coefficient determining how significant template $j$ is to the complex population. In Dolphin 2002, he normalizes the templates to identical star formation rates (SFRs) and so the $r_j$ are SFRs as well. In this package, we prefer to normalize our templates to identical population stellar masses, so our $r_j$ are stellar masses, but the principal is the same.
Construction of the templates is, however, not a trivial exercise. Ideally, a template constructed to represent a particular SSP would accurately reflect the expectation of how such a population would be observed. Thus, these templates must be adjusted for photometric error, incompleteness, and other effects such as those caused by unresolved binary- and multi-star systems. Observational effects such as photometric error and incompleteness are best measured from artificial star tests (ASTs). It is worth mentioning that ASTs can often return "best case" results, as they typically neglect systematics like uncertainty in the point-spread model used for the photometry; as such it is sometimes necessary to add a systematic error floor to photometric error results from ASTs.
Such templates can be constructed by sampling many mock stars from an initial mass function (IMF), interpolating their absolute magnitudes from an isochrone of the relevant SSP, and "mock observing" them by applying photometric error and completeness functions (for example, by looking up the $1\sigma$ photometric error and completeness value from a catalog of artificial stars). Such Monte Carlo templates can be slow to construct and exhibit Poisson shot-noise, requiring a statistical data–data comparison rather than a model–data comparison. Thus this method is non-optimal from both a practical and statistical perspective.
It is better to form what Dolphin 2002 calls a "blurred isochrone;" in this form of template, the SSP isochrone is first interpolated in initial stellar mass to improve the point density along the isochrone. The number of interpolated points is generally a function of the size of the bins in the Hess diagram and the observational error; more points are required as the bin size or photometric errors become smaller. These points are then weighted according to the IMF and the photometric completeness, and this weight is distributed into the Hess diagram following the photometric error distribution determined by similar artificial stars. Dolphin 2002 also mentions interpolating across stellar age/metallicity when constructing such templates; for example, for an SSP with an age of 1 Gyr and a metallicity of [M/H]=-1.0, you could interpolate the isochrones to introduce a Gaussian metallicity spread of 0.05 dex or an age spread of 100 Myr. The general effects of this form of interpolation is to broaden the model templates, particularly features that are very sharp in true SSP models. We neglect this form of interpolation in our implementation as it adds significant complexity and requires users to provide more information about the isochrones that are providing. Such widening of the individual templates is most impactful when photometric errors in the observational data are low (perhaps <0.10 mag).
While the above description summarizes the necessary components for constructing such a blurred isochrone, it can be a bit difficult to figure out how best to actually construct them. Specifically there are many ways that one could implement the observational effects of photometric error and incompleteness. We provide a method partial_cmd_smooth to construct such templates under the assumption of Gaussian photometric error distributions, which is often a good approximation in the high-completeness regime. This method makes use of user-defined functions for the mean photometric error and completeness as a function of magnitude and filter, such that these can be defined in a number of ways; for example, as direct lookups from a large table of ASTs or as simple function evaluations of analytic approximations or fits to the ASTs.
This method begins by interpolating the provided SSP isochrone to increase point density. For every such point with $i$ band apparent magnitude $m_i$, it calls a user-defined function to estimate the $1\sigma$ photometric error as $\sigma_i = f_i(m_i)$. The photometric error on the x-axis color for the Hess diagram is estimated from the individual-band $\sigma_i$. These errors are used to define an asymmetric Gaussian kernel for each point in the interpolated isochrone. This kernel describes the shape of the probability distribution for where in the Hess diagram the isochrone point would be observed. However, it also must be normalized (weighted) according to the IMF and observational completeness functions.
Assume that the vector of initial stellar masses for the points in the interpolated isochrone are $m_i$ and that they are sorted such that $m_i < m_{i+1}$. The IMF weight on point $m_i$ can be approximated as the number fraction of stars born between $m_i$ and $m_{i+1}$ divided by the mean mass per star born $\langle m \rangle$, such that the weight effectively represents the number of stars expected to be born with masses between $m_i$ and $m_{i+1}$ per solar mass of star formation:
The numerator can either be calculated as the difference in the cumulative distribution function across the bin or approximated efficiently via the trapezoidal rule. The denominator is a function only of the IMF and need only be calculated once. Multiplying this weight by the probability of detection in the relevant bands gives the final weight.
Below we show a comparison of a smooth Hess diagram template constructed with partial_cmd_smooth with a Monte Carlo realization created with generate_stars_mass and mock-observed with model_cmd. These use an SSP isochrone of age 10 Gyr and metallicity [M/H] of -2 from PARSEC with identical observational error and completeness models. For the provided stellar mass of $10^7 \, \text{M}_\odot$, the Monte Carlo model is fairly well-sampled but still noticably noisy in regions of the Hess diagram that are less well-populated.
result::StatsBase.Histogram = partial_cmd_smooth( m_ini::AbstractVector{<:Number}, mags::AbstractVector{<:AbstractVector{<:Number}}, mag_err_funcs, y_index, color_indices, imf, completeness_funcs=[one for i in mags]; dmod::Number=0, normalize_value::Number=1, mean_mass=mean(imf), edges=nothing, xlim=nothing, ylim=nothing, nbins=nothing, xwidth=nothing, ywidth=nothing )
Main function for generating template Hess diagrams from a simple stellar population of stars from an isochrone, including photometric error and completeness.
Arguments
m_ini::AbstractVector{<:Number} is a vector containing the initial stellar masses of the stars from the isochrone.
mags::AbstractVector{<:AbstractVector{<:Number}} is a vector of vectors. Each constituent vector with index i should have length(mags[i]) == length(m_ini), representing the magnitudes of the isochrone stars in each of the magnitudes considered. In most cases, mags should contain 2 (if y-axis mag is also involved in the x-axis color) or 3 vectors.
mag_err_funcs must be an indexable object (e.g., a vector or tuple) that contains callables (e.g., a Function) to compute the 1σ photometric errors in the filters provided in mags. Each callable must take a single argument and return a Number. The length mag_err_funcs must be equal to the length of mags.
y_index gives a valid index (e.g., an Int or CartesianIndex) into mags for the filter you want to have on the y-axis of the Hess diagram. For example, if the mags argument contains the B and V band magnitudes as mags=[B, V] and you want V on the y-axis, you would set y_index as 2.
color_indices is a length-2 indexable object giving the indices into mags that are to be used to compute the x-axis color. For example, if the mags argument contains the B and V band magnitudes as mags=[B, V], and you want B-V to be the x-axis color, then color_indices should be [1,2] or (1,2) or similar.
imf is a callable that takes an initial stellar mass as its sole argument and returns the (properly normalized) probability density of your initial mass function model. All the models from InitialMassFunctions.jl are valid for imf.
completeness_functions must be an indexable object (e.g., a vector or tuple) that contains callables (e.g., a Function) to compute the single-filter completeness fractions as a function of magnitude. Each callable in this argument must correspond to the matching filter provided in mags.
Keyword Arguments
dmod::Number=0 distance modulus in magnitudes to apply to the input mags.
normalize_value::Number=1 gives the total stellar mass of the population you wish to model.
mean_mass::Number gives the expectation value for a random star drawn from your provided imf. This will be computed for you if your provided imf is a valid continuous, univariate Distributions.Distribution object.
edges is a tuple of vector-like objects 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 overrides the following keyword arguments that offer other ways to specify the extent of the Hess diagram.
xlim; a length-2 indexable object (e.g., a vector or tuple) 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 is not provided.
ylim; as xlim but for the y-axis corresponding to the provided mags array. Example [25.0, 20.0]. This is only used if edges is not provided.
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 is not provided.
xwidth; the bin width along the x-axis for the colors array. This is only used if edges and nbins are not provided. Example: 0.1.
ywidth; as xwidth but for the y-axis corresponding to the provided mags array. Example: 0.1.
Returns
This method returns the Hess diagram as a StatsBase.Histogram; you should refer to the StatsBase documentation for more information. In short, if the output of this method is result, then the Hess diagram represented as a Matrix is available as result.weights (this is what you would want for fit_templates and similar functions) and the edges of the histogram are available as result.edges.
We note that in many cases it can also be helpful to add in a foreground/background template that models contamination of the Hess diagram from stars not in your population of interest – this is often done using observations of parallel fields though there are several other possible methods.
Photometric catalogs can be processed into Hess diagrams meeting our formatting requirements with the method bin_cmd.
The numerator can either be calculated as the difference in the cumulative distribution function across the bin or approximated efficiently via the trapezoidal rule. The denominator is a function only of the IMF and need only be calculated once. Multiplying this weight by the probability of detection in the relevant bands gives the final weight.
Below we show a comparison of a smooth Hess diagram template constructed with partial_cmd_smooth with a Monte Carlo realization created with generate_stars_mass and mock-observed with model_cmd. These use an SSP isochrone of age 10 Gyr and metallicity [M/H] of -2 from PARSEC with identical observational error and completeness models. For the provided stellar mass of $10^7 \, \text{M}_\odot$, the Monte Carlo model is fairly well-sampled but still noticably noisy in regions of the Hess diagram that are less well-populated.
result::StatsBase.Histogram = partial_cmd_smooth( m_ini::AbstractVector{<:Number}, mags::AbstractVector{<:AbstractVector{<:Number}}, mag_err_funcs, y_index, color_indices, imf, completeness_funcs=[one for i in mags]; dmod::Number=0, normalize_value::Number=1, mean_mass=mean(imf), edges=nothing, xlim=nothing, ylim=nothing, nbins=nothing, xwidth=nothing, ywidth=nothing )
Main function for generating template Hess diagrams from a simple stellar population of stars from an isochrone, including photometric error and completeness.
Arguments
m_ini::AbstractVector{<:Number} is a vector containing the initial stellar masses of the stars from the isochrone.
mags::AbstractVector{<:AbstractVector{<:Number}} is a vector of vectors. Each constituent vector with index i should have length(mags[i]) == length(m_ini), representing the magnitudes of the isochrone stars in each of the magnitudes considered. In most cases, mags should contain 2 (if y-axis mag is also involved in the x-axis color) or 3 vectors.
mag_err_funcs must be an indexable object (e.g., a vector or tuple) that contains callables (e.g., a Function) to compute the 1σ photometric errors in the filters provided in mags. Each callable must take a single argument and return a Number. The length mag_err_funcs must be equal to the length of mags.
y_index gives a valid index (e.g., an Int or CartesianIndex) into mags for the filter you want to have on the y-axis of the Hess diagram. For example, if the mags argument contains the B and V band magnitudes as mags=[B, V] and you want V on the y-axis, you would set y_index as 2.
color_indices is a length-2 indexable object giving the indices into mags that are to be used to compute the x-axis color. For example, if the mags argument contains the B and V band magnitudes as mags=[B, V], and you want B-V to be the x-axis color, then color_indices should be [1,2] or (1,2) or similar.
imf is a callable that takes an initial stellar mass as its sole argument and returns the (properly normalized) probability density of your initial mass function model. All the models from InitialMassFunctions.jl are valid for imf.
completeness_functions must be an indexable object (e.g., a vector or tuple) that contains callables (e.g., a Function) to compute the single-filter completeness fractions as a function of magnitude. Each callable in this argument must correspond to the matching filter provided in mags.
Keyword Arguments
dmod::Number=0 distance modulus in magnitudes to apply to the input mags.
normalize_value::Number=1 gives the total stellar mass of the population you wish to model.
mean_mass::Number gives the expectation value for a random star drawn from your provided imf. This will be computed for you if your provided imf is a valid continuous, univariate Distributions.Distribution object.
edges is a tuple of vector-like objects 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 overrides the following keyword arguments that offer other ways to specify the extent of the Hess diagram.
xlim; a length-2 indexable object (e.g., a vector or tuple) 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 is not provided.
ylim; as xlim but for the y-axis corresponding to the provided mags array. Example [25.0, 20.0]. This is only used if edges is not provided.
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 is not provided.
xwidth; the bin width along the x-axis for the colors array. This is only used if edges and nbins are not provided. Example: 0.1.
ywidth; as xwidth but for the y-axis corresponding to the provided mags array. Example: 0.1.
Returns
This method returns the Hess diagram as a StatsBase.Histogram; you should refer to the StatsBase documentation for more information. In short, if the output of this method is result, then the Hess diagram represented as a Matrix is available as result.weights (this is what you would want for fit_templates and similar functions) and the edges of the histogram are available as result.edges.
We note that in many cases it can also be helpful to add in a foreground/background template that models contamination of the Hess diagram from stars not in your population of interest – this is often done using observations of parallel fields though there are several other possible methods.
Photometric catalogs can be processed into Hess diagrams meeting our formatting requirements with the method bin_cmd.
Returns a StatsBase.Histogram type containing the Hess diagram from the provided x-axis photometric colors and y-axis photometric magnitudes mags. These must all be vectors equal in length. You can either specify the bin edges directly via the edges keyword (e.g., edges = (range(-0.5, 1.6, length=100), range(17.0, 26.0, length=100))), or you can set the x- and y-limits via xlim and ylim and the number of bins as nbins, or you can omit nbins and instead pass the bin width in the x and y directions, xwidth and ywidth. See below for more info on the keyword arguments. To plot this with PyPlot you should do plt.imshow(result.weights', origin="lower", ...).
Keyword Arguments
weights::AbstractVector{<:Number} is a array of length equal to colors and mags that contains the probabilistic weights associated with each point. This is passed to StatsBase.fit as StatsBase.Weights(weights).
edges is a tuple of vector-like objects 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, weights is the only other keyword that will be read; edges supercedes the other construction methods.
xlim; a length-2 indexable object (e.g., a vector or tuple) 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 is not provided.
ylim; as xlim but for the y-axis corresponding to the provided mags array. Example [25.0, 20.0]. This is only used if edges is not provided.
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 is not provided.
xwidth; the bin width along the x-axis for the colors array. This is only used if edges and nbins are not provided. Example: 0.1.
ywidth; as xwidth but for the y-axis corresponding to the provided mags array. Example: 0.1.
It is expected that the user will typically have model templates stored as two-dimensional matrices as these are the obvious choice for representing a binned two-dimensional histogram. We fully support supplying the list of model templates as a list of matrices (e.g., a Vector{Matrix{<:Number}}) to the fitting functions discussed below. The important computational kernels composite! and ∇loglikelihood! have custom loops for these input types.
However, additional optimizations are possible by flattening the data. By flattening each matrix in the list of model templates into a column vector and concatenating them such that the list of model templates becomes a single matrix, we can compute the complex model Hess diagram as a single matrix-vector product rather than using a custom loop. The same optimization can be made when computing the gradient of the loglikelihood (discussed more below). The majority of all computation for the fitting methods below is spent in these two functions, so optimizing their performance translates directly to improved fitting runtimes. With this flattened memory arrangement we can use the highly optimized LinearAlgbera.mul! method to do the in-place matrix-vector product. This will typically be translated into a call to a BLAS function like gemv!. As such, we can benefit from Julia's ability to switch BLAS implementations at runtime to use Intel's Math Kernel Library, Apple's Accelerate, and others.
Most of the fitting methods below support both the natural and flattened data layouts. We provide the stack_models method to produce the optimized layout for the list of model templates.
Transforms a vector of matrices into a single matrix, with each matrix from models being transcribed into a single column in the output matrix. This data layout enables more efficient calculations in some of our internal functions like composite! and ∇loglikelihood!.
Examples
julia> stack_models([rand(5,5) for i in 1:10])
+ ywidth = nothing)
Returns a StatsBase.Histogram type containing the Hess diagram from the provided x-axis photometric colors and y-axis photometric magnitudes mags. These must all be vectors equal in length. You can either specify the bin edges directly via the edges keyword (e.g., edges = (range(-0.5, 1.6, length=100), range(17.0, 26.0, length=100))), or you can set the x- and y-limits via xlim and ylim and the number of bins as nbins, or you can omit nbins and instead pass the bin width in the x and y directions, xwidth and ywidth. See below for more info on the keyword arguments. To plot this with PyPlot you should do plt.imshow(result.weights', origin="lower", ...).
Keyword Arguments
weights::AbstractVector{<:Number} is a array of length equal to colors and mags that contains the probabilistic weights associated with each point. This is passed to StatsBase.fit as StatsBase.Weights(weights).
edges is a tuple of vector-like objects 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, weights is the only other keyword that will be read; edges supercedes the other construction methods.
xlim; a length-2 indexable object (e.g., a vector or tuple) 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 is not provided.
ylim; as xlim but for the y-axis corresponding to the provided mags array. Example [25.0, 20.0]. This is only used if edges is not provided.
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 is not provided.
xwidth; the bin width along the x-axis for the colors array. This is only used if edges and nbins are not provided. Example: 0.1.
ywidth; as xwidth but for the y-axis corresponding to the provided mags array. Example: 0.1.
It is expected that the user will typically have model templates stored as two-dimensional matrices as these are the obvious choice for representing a binned two-dimensional histogram. We fully support supplying the list of model templates as a list of matrices (e.g., a Vector{Matrix{<:Number}}) to the fitting functions discussed below. The important computational kernels composite! and ∇loglikelihood! have custom loops for these input types.
However, additional optimizations are possible by flattening the data. By flattening each matrix in the list of model templates into a column vector and concatenating them such that the list of model templates becomes a single matrix, we can compute the complex model Hess diagram as a single matrix-vector product rather than using a custom loop. The same optimization can be made when computing the gradient of the loglikelihood (discussed more below). The majority of all computation for the fitting methods below is spent in these two functions, so optimizing their performance translates directly to improved fitting runtimes. With this flattened memory arrangement we can use the highly optimized LinearAlgbera.mul! method to do the in-place matrix-vector product. This will typically be translated into a call to a BLAS function like gemv!. As such, we can benefit from Julia's ability to switch BLAS implementations at runtime to use Intel's Math Kernel Library, Apple's Accelerate, and others.
Most of the fitting methods below support both the natural and flattened data layouts. We provide the stack_models method to produce the optimized layout for the list of model templates.
Transforms a vector of matrices into a single matrix, with each matrix from models being transcribed into a single column in the output matrix. This data layout enables more efficient calculations in some of our internal functions like composite! and ∇loglikelihood!.
Examples
julia> stack_models([rand(5,5) for i in 1:10])
25×10 Matrix{Float64}:
-...
While generally the methods using BLAS routines offer significant performance improvements, there is a caveat when multithreading from within Julia. By default Julia will allow BLAS to use multiple threads even if Julia itself is started with a single thread (i.e., by running julia -t 1). BLAS threads do not compose with Julia threads. That is, if you start Julia with N>1 threads (julia -t N) and write a threaded workload where each Julia thread is doing BLAS operations concurrently, you can easily oversubscribe the CPU. Specific recommendations vary depending on BLAS vendor (see this page and the linked discourse threads), but generally this package is in the regime of doing many small calculations that do not individually benefit much from BLAS threading (e.g., performance for OpenBLAS with 8 threads is only ~2x the 1 thread performance). As such it is often sufficient to set BLAS to use a single thread (via LinearAlgbera.BLAS.set_num_threads(1) or environment variables; see above link).
Settings
This document was generated with Documenter.jl version 1.4.1 on Monday 20 May 2024. Using Julia version 1.10.3.
While generally the methods using BLAS routines offer significant performance improvements, there is a caveat when multithreading from within Julia. By default Julia will allow BLAS to use multiple threads even if Julia itself is started with a single thread (i.e., by running julia -t 1). BLAS threads do not compose with Julia threads. That is, if you start Julia with N>1 threads (julia -t N) and write a threaded workload where each Julia thread is doing BLAS operations concurrently, you can easily oversubscribe the CPU. Specific recommendations vary depending on BLAS vendor (see this page and the linked discourse threads), but generally this package is in the regime of doing many small calculations that do not individually benefit much from BLAS threading (e.g., performance for OpenBLAS with 8 threads is only ~2x the 1 thread performance). As such it is often sufficient to set BLAS to use a single thread (via LinearAlgbera.BLAS.set_num_threads(1) or environment variables; see above link).
Settings
This document was generated with Documenter.jl version 1.4.1 on Monday 20 May 2024. Using Julia version 1.10.3.
