Merge branch 'bd-j:main' into main

Author: wangbingjie

doc/prospector-beta_priors.rst

@@ -2,13 +2,19 @@
 Prospector-beta Priors
 ==============
 
-This model is intended for fitting galaxy photometry where the redshift is unknown.
-The priors encode empirical constraints of redshifts, masses, and star formation histories in the galaxy population.
+This model is intended for fitting galaxy photometry where the redshift is
+unknown. The priors encode empirical constraints of redshifts, masses, and star
+formation histories in the galaxy population.
 
-N.B.: Please make sure to update to the Prospector-beta version post commit `09a83f2 <https://github.com/bd-j/prospector/commit/09a83f28cae3bcc0f0397b3a0b8d85aa4f96bf12>`_, merged on May 19, 2023. This is a major update to the SFH(M, z) prior, ensuring the expectation values of SFRs are well-behaved over the full prior space.
+N.B.: Please make sure to update to the Prospector-beta version post commit
+`09a83f2
+<https://github.com/bd-j/prospector/commit/09a83f28cae3bcc0f0397b3a0b8d85aa4f96bf12>`_,
+merged on May 19, 2023. This is a major update to the SFH(M, z) prior, ensuring
+the expectation values of SFRs are well-behaved over the full prior space.
 
-A set of prospector parameters implementing the full set of priors is available as the ``"beta"`` entry
-of :py:class:`prospect.models.templates.TemplateLibrary`.
+A set of prospector parameters implementing the full set of priors is available
+as the ``"beta"`` entry of
+:py:class:`prospect.models.templates.TemplateLibrary`.
 
 Additionally we provide different combinations of the priors for flexibility, which includes the following:
 
@@ -20,43 +26,71 @@ Additionally we provide different combinations of the priors for flexibility, wh
 * ``PhiSFHfixZred``  : same as above, but keeping zred fixed to a user-specified value, 'zred', during fitting
 * ``NzSFH``          : number density + mass function + mass-met + SFH(M, z); this is the full set of Prospector-beta priors.
 
-We describe each of the priors briefly below. Please cite `wang23 <https://ui.adsabs.harvard.edu/abs/2023ApJ...944L..58W/abstract>`_ and the relevant papers if any of the priors are used.
+We describe each of the priors briefly below. Please cite `wang23
+<https://ui.adsabs.harvard.edu/abs/2023ApJ...944L..58W/abstract>`_ and the
+relevant papers if any of the priors are used.
 
 
 Stellar Mass Function
------------
+---------------------
 
 Two options are available, the choice of which depends on the given scientific question.
 The relevant data files are ``prior_data/pdf_of_z_l20.txt`` & ``prior_data/pdf_of_z_l20t18.txt``.
 These mass functions can also be replaced by supplying new data files to ``prior_data/``
 
 1. ``"const_phi = True"``
 
-The mass functions between 0.2 ≤ z ≤ 3.0 are taken from `leja20 <https://ui.adsabs.harvard.edu/abs/2020ApJ...893..111L/abstract>`_. Outside this redshift range, we adopt a nearest-neighbor solution, i.e., the z = 0.2 and z = 3 mass functions.
+The mass functions between 0.2 ≤ z ≤ 3.0 are taken from `leja20
+<https://ui.adsabs.harvard.edu/abs/2020ApJ...893..111L/abstract>`_. Outside this
+redshift range, we adopt a nearest-neighbor solution, i.e., the z = 0.2 and z =
+3 mass functions.
 
 2. ``"const_phi = False"``
 
-The mass functions are switched to those from `tacchella18 <https://ui.adsabs.harvard.edu/abs/2018ApJ...868...92T/abstract>`_ between 4 < z < 12, with the 3 < z < 4 transition from `leja20 <https://ui.adsabs.harvard.edu/abs/2020ApJ...893..111L/abstract>`_ managed with a smoothly-varying average in number density space. We use the z = 12 mass function for z > 12.
+The mass functions are switched to those from `tacchella18
+<https://ui.adsabs.harvard.edu/abs/2018ApJ...868...92T/abstract>`_ between 4 < z
+< 12, with the 3 < z < 4 transition from `leja20
+<https://ui.adsabs.harvard.edu/abs/2020ApJ...893..111L/abstract>`_ managed with
+a smoothly-varying average in number density space. We use the z = 12 mass
+function for z > 12.
 
 
 Galaxy Number Density
 -----------
 
-This prior informs the model about the survey volume being probed. It is sensitive to the mass-completeness limit of the data. We provide a default setting derived from a mock JWST catalog, which is contained in ``prior_data/mc_from_mocks.txt``.
-In practice one would likely need to obtain the mass-completeness limits from using SED-modeling heuristics based on the flux-completeness limits in a given catalog.
+This prior informs the model about the survey volume being probed. It is
+sensitive to the mass-completeness limit of the data. We provide a default
+setting derived from a mock JWST catalog, which is contained in
+``prior_data/mc_from_mocks.txt``.
+In practice one would likely need to obtain
+the mass-completeness limits from using SED-modeling heuristics based on the
+flux-completeness limits in a given catalog.
 
 
 Dynamic Star-formation History
 -----------
 
-The SFH is described non-parametrically as in Prospector-alpha; the number of age bins is set by ``"nbins_sfh"``.
+The SFH is described non-parametrically as in Prospector-alpha; the number of
+age bins is set by ``"nbins_sfh"``.
 
-In contrast to the null expectation assumed in Prospector-alpha, the expectation value in each age bin is matched to the cosmic star formation rate densities in `behroozi19 <https://ui.adsabs.harvard.edu/abs/2019MNRAS.488.3143B/abstract>`_, while the distribution about the mean remains to be the Student’s-t distribution. The sigma of the Student’s-t distribution is set by ``"logsfr_ratio_tscale"``, and the range is clipped to be within ``"logsfr_ratio_mini"`` and ``"logsfr_ratio_maxi"``.
+In contrast to the null expectation assumed in Prospector-alpha, the expectation
+value in each age bin is matched to the cosmic star formation rate densities in
+`behroozi19 <https://ui.adsabs.harvard.edu/abs/2019MNRAS.488.3143B/abstract>`_,
+while the distribution about the mean remains to be the Student's-t
+distribution. The sigma of the Student's-t distribution is set by
+``"logsfr_ratio_tscale"``, and the range is clipped to be within
+``"logsfr_ratio_mini"`` and ``"logsfr_ratio_maxi"``.
 
-A simple mass dependence on SFH is further introduced by shifting the start of the age bins as a function of mass. This SFH prior effectively encodes an expectation that high-mass galaxies form earlier, and low-mass galaxies form later, than naive expectations from the cosmic SFRD.
+A simple mass dependence on SFH is further introduced by shifting the start of
+the age bins as a function of mass. This SFH prior effectively encodes an
+expectation that high-mass galaxies form earlier, and low-mass galaxies form
+later, than naive expectations from the cosmic SFRD.
 
 
-Stellar Mass–Stellar Metallicity
+Stellar Mass-Stellar Metallicity
 -----------
 
-This is the stellar mass–stellar metallicity relationship measured from the SDSS in `gallazzi05 <https://ui.adsabs.harvard.edu/abs/2005MNRAS.362...41G/abstract>`_, introduced in `leja19 <https://ui.adsabs.harvard.edu/abs/2019ApJ...876....3L/abstract>`_.
+This is the stellar mass-stellar metallicity relationship measured from the SDSS
+in `gallazzi05
+<https://ui.adsabs.harvard.edu/abs/2005MNRAS.362...41G/abstract>`_, introduced
+in `leja19 <https://ui.adsabs.harvard.edu/abs/2019ApJ...876....3L/abstract>`_.

doc/sfhs.rst

@@ -77,16 +77,17 @@ Continuity SFH
 ^^^^^^^^^^^^^^
 See `leja19 <https://ui.adsabs.harvard.edu/abs/2019ApJ...876....3L/abstract>`_,
 `johnson21 <https://ui.adsabs.harvard.edu/abs/2021ApJS..254...22J/abstract>`_
-for more details. A set of propsector parameters implementing this treatment
-with 3 bins is available as the ``"continuity_sfh"`` entry of
+for more details. A basic set of prospector parameters implementing this
+treatment with 3 bins is available as the ``"continuity_sfh"`` entry of
 :py:class:`prospect.models.templates.TemplateLibrary`.
 
 In this parameterization, the SFR of each bin is derived from sampling a vector
 of parameters describing the *ratio* of SFRs in adjacent temporal bins.  By
 default, a Student-t prior distribution (like a Gaussian but with heavier tails)
 is placed on the log of these ratios.  This results in a prior SFH that tends
 toward constant SFR, and down-weights drmamtic changes in the SFR between
-adjacent bins.  The overall normalization is provided by the ``logmass``
+adjacent bins.  The width of the distribution can be ajusted to produce smoother
+or burstier SFHs.The overall normalization is provided by the ``logmass``
 parameter.
 
 In detail, the SFR in each timetime is computed as
@@ -140,6 +141,18 @@ required to model post-starburst SFHs. A set of prospector parameters
 implementing this treatment is available as the ``"continuity_psb_sfh"`` entry
 of :py:class:`prospect.models.templates.TemplateLibrary`
 
+
+'Stochastic' SFH
+^^^^^^^^^^^^^^^^
+This SFH (hyper-)prior uses the power spectrum of SFH fluctuations -- the
+parameters of which can be sampled -- to determine the covariance matrix between
+(adjacent and non-adjacent) temporal bins of SFR. See `Wan et al. 24 <https://ui.adsabs.harvard.edu/abs/2024arXiv240414494W/abstract>`_ for
+details.  This prior is adapted from the Extended Regulator model developed in
+`Caplar & Tacchella (2019) <https://ui.adsabs.harvard.edu/abs/2019MNRAS.487.3845C/abstract>`_ and `Tacchella, Forbes & Caplar (2020) <https://ui.adsabs.harvard.edu/abs/2020MNRAS.497..698T/abstract>`_ , in
+conjunction with the GP implementation of `Iyer & Speagle et al. (2024) <https://ui.adsabs.harvard.edu/abs/2024ApJ...961...53I/abstract>`_ taken from `this module <https://github.com/kartheikiyer/GP-SFH>`_ .
+
+
+
 Dirichlet SFH
 ^^^^^^^^^^^^^
 See `leja17 <https://ui.adsabs.harvard.edu/abs/2017ApJ...837..170L/abstract>`_,

prospect/models/hyperparam_transforms.py

@@ -0,0 +1,183 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+"""hyperparam_transforms.py -- This module contains parameter transformations that are
+used in the stochastic SFH prior.
+
+These are taken from the implementation of
+https://ui.adsabs.harvard.edu/abs/2024ApJ...961...53I/abstract given in
+https://github.com/kartheikiyer/GP-SFH
+
+They can be used as ``"depends_on"`` entries in parameter specifications.
+"""
+
+import numpy as np
+from astropy.cosmology import FlatLambdaCDM
+
+__all__ = ["get_sfr_covar", "sfr_covar_to_sfr_ratio_covar"]
+
+
+# --------------------------------------
+# --- Functions/transforms for stochastic SFH prior ---
+# --------------------------------------
+
+
+# Creating a base class that simplifies a lot of things.
+# The way this is set up, you can pass a kernel as an argument
+# to compute the covariance matrix and draw samples from it.
+class simple_GP_sfh():
+
+    """
+    A class that creates and holds information about a specific
+    kernel, and can generate samples from it.
+
+    From https://github.com/kartheikiyer/GP-SFH
+
+    Attributes
+    ----------
+    tarr: fiducial time array used to draw samples
+    kernel: accepts an input function as an argument,
+        of the format:
+
+            def kernel_function(delta_t, **kwargs):
+                ... function interior ...
+                return kernel_val[array of len(delta_t)]
+
+    Methods
+    -------
+    get_covariance_matrix
+        [although this has double for loops for maximum flexibility
+        with generic kernel functions, it only has to be computed once,
+        which makes drawing random samples super fast once it's computed.]
+    sample_kernel
+    plot_samples
+    plot_kernel
+    [to-do] condition on data
+
+    """
+
+    def __init__(self, sp = 'none', cosmo = FlatLambdaCDM(H0=70, Om0=0.3), zval = 0.1):
+
+
+        self.kernel = []
+        self.covariance_matrix = []
+        self.zval = zval
+        self.sp = sp
+        self.cosmo = cosmo
+        self.get_t_univ()
+        self.get_tarr()
+
+
+    def get_t_univ(self):
+
+        self.t_univ = self.cosmo.age(self.zval).value
+        return
+
+    def get_tarr(self, n_tarr = 1000):
+
+        self.get_t_univ()
+        if n_tarr > 1:
+            self.tarr = np.linspace(0,self.t_univ, n_tarr)
+        elif n_tarr < 1:
+            self.tarr = np.arange(0,self.t_univ, n_tarr)
+        else:
+            raise('Undefined n_tarr: expected int or float.')
+        return
+
+
+    def get_covariance_matrix(self, show_prog = True, **kwargs):
+        """
+        Evaluate covariance matrix with a particular kernel
+        """
+
+        cov_matrix = np.zeros((len(self.tarr),len(self.tarr)))
+
+        if show_prog == True:
+            iterrange = range(len(cov_matrix))
+        else:
+            iterrange = range(len(cov_matrix))
+        for i in iterrange:
+            for j in range(len(cov_matrix)):
+                    cov_matrix[i,j] = self.kernel(self.tarr[i] - self.tarr[j], **kwargs)
+
+        return cov_matrix
+
+
+
+def extended_regulator_model_kernel_paramlist(delta_t, kernel_params, base_e_to_10 = False):
+    """
+    A basic implementation of the regulator model kernel, with five parameters:
+    kernel_params = [sigma, tau_eq, tau_in, sigma_gmc, tau_gmc]
+    sigma: \sigma, the amount of overall variance
+    tau_eq: equilibrium timescale
+    tau_x: inflow correlation timescale (includes 2pi factor)
+    sigma_gmc: gmc variability
+    tau_l: cloud lifetime
+
+    from https://github.com/kartheikiyer/GP-SFH
+    """
+
+    sigma, tau_eq, tau_in, sigma_gmc, tau_gmc = kernel_params
+
+    if base_e_to_10 == True:
+        # in TCF20, this is defined in base e, so convert to base 10
+        sigma = sigma*np.log10(np.e)
+        sigma_gmc = sigma_gmc*np.log10(np.e)
+
+    tau = np.abs(delta_t)
+
+    if tau_in == tau_eq:
+        c_reg = sigma**2 * (1 + tau/tau_eq) * (np.exp(-tau/tau_eq))
+    else:
+        c_reg = sigma**2 / (tau_in - tau_eq) * (tau_in*np.exp(-tau/tau_in) - tau_eq*np.exp(-tau/tau_eq))
+
+    c_gmc = sigma_gmc**2 * np.exp(-tau/tau_gmc)
+
+    kernel_val = (c_reg + c_gmc)
+    return kernel_val
+
+
+
+def get_sfr_covar(psd_params, agebins=[], **extras):
+    """Caluclates SFR covariance matrix for a given set of PSD parameters and agebins
+    PSD parameters must be in the order: [sigma_reg, tau_eq, tau_in, sigma_dyn, tau_dyn]
+
+    from https://github.com/kartheikiyer/GP-SFH
+
+    Returns
+    -------
+    covar_matrix: (Nbins, Nbins)-dim array of covariance values for SFR
+    """
+
+    bincenters = np.array([np.mean(agebins[i]) for i in range(len(agebins))])
+    bincenters = (10**bincenters)/1e9
+    case1 = simple_GP_sfh()
+    case1.tarr = bincenters
+    case1.kernel = extended_regulator_model_kernel_paramlist
+    covar_matrix = case1.get_covariance_matrix(kernel_params = psd_params, show_prog=False)
+
+    return covar_matrix
+
+
+def sfr_covar_to_sfr_ratio_covar(covar_matrix):
+    """Caluclates log SFR ratio covariance matrix from SFR covariance matrix
+
+    from https://github.com/kartheikiyer/GP-SFH
+
+    Returns
+    -------
+    sfr_ratio_covar: (Nbins-1, Nbins-1)-dim array of covariance values for log SFR
+    """
+
+    dim = covar_matrix.shape[0]
+
+    sfr_ratio_covar = []
+
+    for i in range(dim-1):
+        row = []
+        for j in range(dim-1):
+            cov = covar_matrix[i][j] - covar_matrix[i+1][j] - covar_matrix[i][j+1] + covar_matrix[i+1][j+1]
+            row.append(cov)
+        sfr_ratio_covar.append(row)
+
+    return np.array(sfr_ratio_covar)