diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json
index 3278ac7..164a0c1 100644
--- a/dev/.documenter-siteinfo.json
+++ b/dev/.documenter-siteinfo.json
@@ -1 +1 @@
-{"documenter":{"julia_version":"1.10.3","generation_timestamp":"2024-06-04T17:16:40","documenter_version":"1.4.1"}}
\ No newline at end of file
+{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-06-05T13:56:49","documenter_version":"1.4.1"}}
\ No newline at end of file
diff --git a/dev/binaries/index.html b/dev/binaries/index.html
index ef05b5c..7757d2c 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 Tuesday 4 June 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 Wednesday 5 June 2024. Using Julia version 1.10.4.
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.
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 Tuesday 4 June 2024. Using Julia version 1.10.3.