1.0 release

.github/workflows/CI.yml

@@ -63,7 +63,7 @@ jobs:
       - uses: julia-actions/cache@v2	  
       - uses: julia-actions/julia-buildpkg@v1
       - run: julia --color=yes -e 'using Pkg; Pkg.add("Aqua")'
-      - run: julia --project=@. --color=yes -e 'using Aqua, StarFormationHistories; Aqua.test_all(StarFormationHistories; ambiguities=false, piracies=(treat_as_own=[mean],))'
+      - run: julia --project=@. --color=yes -e 'using Aqua, StarFormationHistories; Aqua.test_all(StarFormationHistories; ambiguities=true, piracies=(treat_as_own=[mean],))'
 
   docs:
     name: Documentation

Project.toml

@@ -1,7 +1,7 @@
 name = "StarFormationHistories"
 uuid = "774d3ad4-2d55-48fa-ab42-467705ebff24"
 authors = ["cgarling <[email protected]>"]
-version = "0.2.0"
+version = "1.0.0"
 
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
@@ -26,6 +26,7 @@ Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+TaskLocalValues = "ed4db957-447d-4319-bfb6-7fa9ae7ecf34"
 Trapz = "592b5752-818d-11e9-1e9a-2b8ca4a44cd1"
 VectorizationBase = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
 
@@ -66,6 +67,7 @@ SpecialFunctions = "1.2, 2"
 StableRNGs = "1"
 StaticArrays = "1"
 StatsBase = "0.32, 0.33, 0.34"
+TaskLocalValues = "0.1"
 Test = "<0.0.1, 1"
 Trapz = "2"
 TypedTables = "1"

README.md

@@ -7,10 +7,10 @@
 [![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
 
 
-This package implements methods for modelling observed Hess diagrams (which are just binned color-magnitude diagrams, also called CMDs) and using them to fit astrophysical star formation histories (SFHs). A paper describing the implemented methodologies is currently in preparation. This package also provides utilities for simulating CMDs given input SFHs and photometric error and completeness functions, which can be useful for planning observations and writing proposals. Please see our documentation (linked in the badge above) for more information. A rendered Jupyter notebook with example usage of this package is available [here](https://nbviewer.org/github/cgarling/StarFormationHistories.jl/blob/main/examples/fitting1.ipynb). This package does not currently ship with stellar models or isochrones, these are expected to be provided by the user and can be sourced from online resources like [the CMD webform](http://stev.oapd.inaf.it/cgi-bin/cmd) for PARSEC models.
+This package implements methods for modelling observed Hess diagrams (which are just binned color-magnitude diagrams, also called CMDs) and using them to fit astrophysical star formation histories (SFHs). The paper describing the implemented methodologies has been accepted to ApJS (see below). This package also provides utilities for simulating CMDs given input SFHs and photometric error and completeness functions, which can be useful for planning observations and writing proposals. Please see our documentation (linked in the badge above) for more information. A rendered Jupyter notebook with example usage of this package is available [here](https://nbviewer.org/github/cgarling/StarFormationHistories.jl/blob/main/examples/fitting1.ipynb). This package does not currently ship with stellar models or isochrones, these are expected to be provided by the user and can be sourced from online resources like [the CMD webform](http://stev.oapd.inaf.it/cgi-bin/cmd) for PARSEC models.
 
 ## Paper
-The paper describing the methodologies used in this package is available on ArXiv [here](http://arxiv.org/abs/2407.19534). We ask that published work using this package cite this paper.
+The initial paper describing the methodologies used in this package is available on ArXiv [here](http://arxiv.org/abs/2407.19534) with permanent DOI identifier [here](https://doi.org/10.3847/1538-4365/adbb64). We ask that published work using this package cite this paper.
 
 ## Installation
 

examples/fitting1.ipynb

@@ -3968,7 +3968,7 @@
     "# Walker initial positions sampled from a multivariate normal distribution with means equal to the BFGS MLE\n",
     "# and identity covariance matrix. The `max` call makes sure none of the cofficients have initial conditions less than 0. \n",
     "mcmc_x0 = [max.(0.0, rand(MvNormal(bfgs_result.mle.μ, I))) for i in 1:nwalkers]\n",
-    "mcmc_result = mcmc_sample(free_templates, data, mcmc_x0, nwalkers, nsteps; nthin=5); # nthin=5 to only save every 5 steps"
+    "mcmc_result = mcmc_sample(free_templates, data, mcmc_x0, nsteps; nthin=5); # nthin=5 to only save every 5 steps"
    ]
   },
   {

examples/log_amr/log_amr_example.jl

@@ -1,4 +1,6 @@
-import StarFormationHistories: Z_from_MH, MH_from_Z, calculate_coeffs_logamr, calculate_αβ_logamr
+# import StarFormationHistories: Z_from_MH, MH_from_Z, calculate_coeffs_logamr, calculate_αβ_logamr
+using StarFormationHistories: Z_from_MH, MH_from_Z, LogarithmicAMR, GaussianDispersion,
+    calculate_coeffs, fittable_params
 using Plots
 gr()
 
@@ -10,16 +12,18 @@ MH = repeat(unique_MH; outer=length(unique_logAge))
 
 # Set coefficients for the logarithmic age-metallicity relation
 T_max::Float64 = 13.7 # exp10(maximum(unique_logAge)-9)
-α::Float64, β::Float64 = calculate_αβ_logamr( (-0.8, 0.0), (-2.5, 13.7), T_max, Z_from_MH )
+MH_model = LogarithmicAMR((-0.8, 0.0), (-2.5, 13.7), T_max)
+α, β = fittable_params(MH_model)
 σ::Float64 = 0.2 # Metallicity spread at fixed logAge; Units of dex
-# <Z>(unique_logAge) for plotting
-plot_Z = α .* (T_max .- exp10.(unique_logAge)./1e9) .+ β
-
+disp_model = GaussianDispersion(σ)
 # Calculate the relative weights such that the sum of all coefficients
-# across a single logAge entry logAge[i] is 1. 
-relative_weights = calculate_coeffs_logamr(ones(length(unique_logAge)), logAge, MH, T_max, α, β, σ)
+# across a single logAge entry logAge[i] is 1.
+relative_weights = calculate_coeffs(MH_model, disp_model, ones(length(unique_logAge)), logAge, MH)
 # Reshape vector into matrix for display
 relweights_matrix = reshape(relative_weights,(length(unique_MH), length(unique_logAge)))
+# <Z>(unique_logAge) for plotting
+plot_Z = α .* (T_max .- exp10.(unique_logAge)./1e9) .+ β
+
 # xbins = vcat(unique_MH .- step(unique_MH)/2,last(unique_MH) + step(unique_MH)/2)
 # ybins = vcat(unique_logAge .- step(unique_logAge)/2,last(unique_logAge) + step(unique_logAge)/2)
 l = @layout [ a{0.33h} ; b{0.33h}; c{0.33h} ]

examples/templates/smooth_template.jl

@@ -21,7 +21,7 @@ savefig = ("DOCSBUILD" in keys(ENV)) && (ENV["DOCSBUILD"] == "true")
 
 # Load example isochrone
 # Path is relative to location of script, so use @__DIR__
-isochrone, mag_names = readdlm(joinpath(@__DIR__, "../../data/isochrone.txt"), ' ',
+isochrone, mag_names = readdlm(joinpath(@__DIR__, "..", "..", "data/isochrone.txt"), ' ',
                                Float64, '\n'; header=true)
 # Unpack
 m_ini = isochrone[:,1]

src/StarFormationHistories.jl

@@ -21,6 +21,7 @@ import ProgressMeter
 using QuadGK: quadgk # For general mean(imf::UnivariateDistribution{Continuous}; kws...)
 using Random: AbstractRNG, default_rng
 using Roots: find_zero # For mass_limits in simulate.jl
+using TaskLocalValues: TaskLocalValue # For threadsafe buffers
 import Trapz: trapz
 using SpecialFunctions: erf
 # LoopVectorization has a specialfunctions extension that provides erf(x::AbstractSIMD)

src/fitting/hierarchical/construct_x0_mdf.jl

@@ -0,0 +1,118 @@
+"""
+    x0::Vector = construct_x0_mdf(logAge::AbstractVector{T},
+                                  [ cum_sfh, ]
+                                  T_max::Number;
+                                  normalize_value::Number = one(T)) where T <: Number
+
+Generates a vector of initial stellar mass normalizations for input to [`fit_sfh`](@ref) and similar methods with a total stellar mass of `normalize_value`. The `logAge` vector must contain the `log10(Age [yr])` of each isochrone that you are going to input as models. If `cum_sfh` is not provided, a constant star formation rate is assumed. For the purposes of computing the constant star formation rate, the provided `logAge` are treated as left-bin edges, with the final right-bin edge being `T_max`, which has units of Gyr. For example, you might have `logAge=[6.6, 6.7, 6.8]` in which case a final logAge of 6.9 would give equal bin widths (in log-space). In this case you would set `T_max = exp10(6.9) / 1e9 ≈ 0.0079` so that the width of the final bin for the star formation rate calculation has the same `log10(Age [yr])` step as the other bins.
+
+A desired cumulative SFH vector `cum_sfh::AbstractVector{<:Number}` can be provided as the second argument, which should correspond to a lookback time vector `unique(logAge)`. You can also provide `cum_sfh` as a length-2 indexable (e.g., a length-2 `Vector{Vector{<:Number}}`) with the first element containing a list of `log10(Age [yr])` values and the second element containing the cumulative SFH values at those values. This cumulative SFH is then interpolated onto the `logAge` provided in the first argument. This method should be used when you want to define the cumulative SFH on a different age grid from the `logAge` you provide in the first argument. The examples below demonstrate these use cases.
+
+The difference between this function and [`StarFormationHistories.construct_x0`](@ref) is that this function generates an `x0` vector that is of length `length(unique(logage))` (that is, a single normalization factor for each unique entry in `logAge`) while [`StarFormationHistories.construct_x0`](@ref) returns an `x0` vector that is of length `length(logAge)`; that is, a normalization factor for every entry in `logAge`. The order of the coefficients is such that the coefficient `x[i]` corresponds to the entry `unique(logAge)[i]`. 
+
+# Notes
+
+# Examples -- Constant SFR
+```jldoctest; setup = :(import StarFormationHistories: construct_x0_mdf, construct_x0)
+julia> isapprox( construct_x0_mdf([9.0, 8.0, 7.0], 10.0; normalize_value=5.0),
+                 [4.504504504504504, 0.4504504504504504, 0.04504504504504504] )
+true
+
+julia> isapprox( construct_x0_mdf(repeat([9.0, 8.0, 7.0, 8.0]; inner=3), 10.0; normalize_value=5.0),
+                 [4.504504504504504, 0.4504504504504504, 0.04504504504504504] )
+true
+
+julia> isapprox( construct_x0_mdf(repeat([9.0, 8.0, 7.0, 8.0]; outer=3), 10.0; normalize_value=5.0),
+                 construct_x0([9.0, 8.0, 7.0], 10.0; normalize_value=5.0) )
+true
+```
+
+# Examples -- Input Cumulative SFH defined on same `logAge` grid
+```jldoctest; setup = :(import StarFormationHistories: construct_x0_mdf)
+julia> isapprox( construct_x0_mdf([9.0, 8.0, 7.0], [0.9009, 0.99099, 1.0], 10.0; normalize_value=5.0),
+                 [4.5045, 0.4504, 0.0450]; atol=1e-3 )
+true
+
+julia> isapprox( construct_x0_mdf([9.0, 8.0, 7.0], [0.1, 0.5, 1.0], 10.0; normalize_value=5.0),
+                 [0.5, 2.0, 2.5] )
+true
+
+julia> isapprox( construct_x0_mdf([7.0, 8.0, 9.0], [1.0, 0.5, 0.1], 10.0; normalize_value=5.0),
+                 [2.5, 2.0, 0.5] )
+true
+```
+
+# Examples -- Input Cumulative SFH with separate `logAge` grid
+```jldoctest; setup = :(import StarFormationHistories: construct_x0_mdf)
+julia> isapprox( construct_x0_mdf([9.0, 8.0, 7.0],
+                                  [[9.0, 8.0, 7.0], [0.9009, 0.99099, 1.0]], 10.0; normalize_value=5.0),
+                 construct_x0_mdf([9.0, 8.0, 7.0], [0.9009, 0.99099, 1.0], 10.0; normalize_value=5.0) )
+true
+
+julia> isapprox( construct_x0_mdf([9.0, 8.0, 7.0],
+                                  [[9.0, 8.5, 8.25, 7.0], [0.9009, 0.945945, 0.9887375, 1.0]], 10.0; normalize_value=5.0),
+                 construct_x0_mdf([9.0, 8.0, 7.0], [0.9009, 0.99099, 1.0], 10.0; normalize_value=5.0) )
+true
+```
+"""
+function construct_x0_mdf(logAge::AbstractVector{T}, T_max::Number; normalize_value::Number=one(T)) where T <: Number
+    minlog, maxlog = extrema(logAge)
+    max_logAge = log10(T_max) + 9 # T_max in units of Gyr
+    @assert max_logAge > maxlog   # max_logAge has to be greater than the maximum of logAge vector
+    sfr = normalize_value / (exp10(max_logAge) - exp10(minlog)) # Average SFR / yr
+    unique_logAge = unique(logAge)
+    idxs = sortperm(unique_logAge)
+    sorted_ul = vcat(unique_logAge[idxs], max_logAge)
+    dt = diff(exp10.(sorted_ul))
+    return [ begin
+                idx = findfirst(Base.Fix1(==, sorted_ul[i]), unique_logAge) # findfirst(==(sorted_ul[i]), unique_logAge)
+                sfr * dt[idx]
+             end for i in eachindex(unique_logAge) ]
+end
+
+function construct_x0_mdf(logAge::AbstractVector{<:Number}, cum_sfh::AbstractVector{T},
+                          T_max::Number; normalize_value::Number=one(T)) where T <: Number
+    cmin, cmax = extrema(cum_sfh)
+    @assert cmin ≥ zero(T)
+    !isapprox(cmax, one(T)) && @warn "Maximum of `cum_sfh` argument is $cmax which is not approximately equal to 1."
+    maximum(cum_sfh) 
+    minlog, maxlog = extrema(logAge)
+    max_logAge = log10(T_max) + 9 # T_max in units of Gyr
+    @assert max_logAge > maxlog   # max_logAge has to be greater than the maximum of logAge vector
+    unique_logAge = unique(logAge)
+    @assert length(unique_logAge) == length(cum_sfh)
+
+    idxs = sortperm(unique_logAge)
+    sorted_cum_sfh = vcat(cum_sfh[idxs], zero(T))
+    # Test that cum_sfh is properly monotonic
+    @assert(all(sorted_cum_sfh[i] ≤ sorted_cum_sfh[i-1] for i in eachindex(sorted_cum_sfh)[2:end]),
+            "Provided `cum_sfh` must be monotonically increasing as `logAge` decreases.")
+    sorted_ul = vcat(unique_logAge[idxs], max_logAge)
+    return [ begin
+                idx = findfirst(Base.Fix1(==, sorted_ul[i]), unique_logAge) # findfirst(==(sorted_ul[i]), unique_logAge)
+                (normalize_value * (sorted_cum_sfh[idx] - sorted_cum_sfh[idx+1]))
+             end for i in eachindex(unique_logAge) ]
+end
+
+# For providing cum_sfh as a length-2 indexable with cum_sfh[1] = log(age) values
+# and cum_sfh[2] = cumulative SFH values at those log(age) values.
+# This interpolates the provided pair onto the provided logAge in first argument.
+function construct_x0_mdf(logAge::AbstractVector{T}, cum_sfh_vec,
+                          T_max::Number; normalize_value::Number=one(T)) where T <: Number
+    @assert(length(cum_sfh_vec) == 2,
+            "`cum_sfh` must either be a vector of numbers or vector containing two vectors that \
+             define a cumulative SFH with a different log(age) discretization than the provided \
+             `logAge` argument.")
+    # Extract cum_sfh info and concatenate with T_max where cum_sfh=0 by definition
+    @assert(maximum(last(cum_sfh_vec)) ≤ one(T),
+            "Maximum of cumulative SFH must be less than or equal to one when passing a custom \
+             cumulative SFH to `construct_x0_mdf`.")
+    idxs = sortperm(first(cum_sfh_vec))
+    cum_sfh_la = vcat(first(cum_sfh_vec)[idxs], log10(T_max) + 9)
+    cum_sfh_in = vcat(last(cum_sfh_vec)[idxs], zero(T))
+    # Construct interpolant and evaluate at unique(logAge)
+    itp = extrapolate(interpolate((cum_sfh_la,), cum_sfh_in, Gridded(Linear())), Flat())
+    cum_sfh = itp.(unique(logAge))
+    # Feed into above method
+    return construct_x0_mdf(logAge, cum_sfh, T_max; normalize_value=normalize_value)
+end

src/fitting/hierarchical/hierarchical_models.jl

@@ -16,12 +16,11 @@ julia> nparams(LinearAMR(1.0, 1.0), GaussianDispersion(0.2))
 """
 nparams(models...) = mapreduce(nparams, +, models)
 
+include("construct_x0_mdf.jl")  # utility function to set initial guess x0
 include("transformations.jl")   # Variable transformations
 include("dispersion_models.jl") # AbstractDispersionModel and subtypes
 include("bfgs_result.jl")       # BFGSResult and CompositeBFGSResult types
 include("fixed_amr.jl")         # Fit under a fixed AMR -- deprecate?
-include("linear_amr/linear_amr.jl")
-include("log_amr/log_amr.jl")
 include("amr.jl")               # Age-metallicity relations
 include("mzr.jl")               # Mass-metallicity relations
 include("generic_fitting.jl")   # Fitting and sampling functions for both AMR and MZR