Final Diffstarpop model

README.md

@@ -26,9 +26,17 @@ $ conda create -c conda-forge -n diffit python=3.11 numpy numba flake8 pytest ja
 Data for this project can be found [at this URL](https://portal.nersc.gov/project/hacc/aphearin/diffstar_data/).
 
 ## Scripts and demo notebooks
+
+The `demo_diffstar_sfh.ipynb` notebook in the `docs` folder illustrates how to use the Diffstar model for individual SFH, and how to generate the SFH of a population of galaxies using DiffstarPop.
+
+The `demo_diffmahpop_diffstarpop_sfh.ipynb` notebook illustrates how to generate a subhalo catalog, and how to generate SFHs for each halo using parameters that reproduce UniverseMachine, IllustrisTNG or Galacticus.
+
+See `diffstar_fitting_script_umachine_mgash.py` for an example of how to fit the SFHs of a large number of simulated galaxies in parallel with mpi4py.
+
 The `diffstar_fitter_demo.ipynb` notebook demonstrates how to fit the SFH of a simulated galaxy with a diffstar approximation.
 
-See `history_fitting_script.py` for an example of how to fit the SFHs of a large number of simulated galaxies in parallel with mpi4py.
+
+See `fit_mstar_ssfr_pdfs_mgash.py` for an example of how to use DiffstarPop to fit a set of Mstar and sSFR PDFs, and `measure_smhm_smdpl_script_mpi_mgash.py` for an example of how to generate the target data from a set of Diffstar fits.
 
 ## Citing diffstar
 [The Diffstar paper](https://arxiv.org/abs/2205.04273) has been published in [Monthly Notices of the Royal Astronomical Society](https://academic.oup.com/mnras/article-abstract/518/1/562/6795944?redirectedFrom=fulltext). Citation information for the paper can be found at [this ADS link](https://ui.adsabs.harvard.edu/abs/2023MNRAS.518..562A/abstract), copied below for convenience:

diffstar/diffstarpop/kernels/__init__.py

@@ -1,4 +1,3 @@
-"""
-"""
+""" """
 
 # flake8: noqa

diffstar/diffstarpop/kernels/params/__init__.py

@@ -7,3 +7,8 @@
     DiffstarPop_Params_Diffstarfits_mgash,
     DiffstarPop_UParams_Diffstarfits_mgash,
 )
+
+from .params_diffstarpopfits_mgash import (
+    DiffstarPop_Params_Diffstarpopfits_mgash,
+    DiffstarPop_UParams_Diffstarpopfits_mgash,
+)

diffstar/diffstarpop/kernels/params/params_diffstarfits_mgash_galacticus_in_plus_ex_situ.py

@@ -1,104 +1,112 @@
-""""""
-
 from collections import OrderedDict, namedtuple
 
-from ..defaults_mgash import DiffstarPopParams, get_unbounded_diffstarpop_params
+from ..defaults_mgash import (
+    DiffstarPopParams,
+    get_unbounded_diffstarpop_params,
+)
+
 from ..satquenchpop_model import DEFAULT_SATQUENCHPOP_PARAMS
 
 SFH_PDF_QUENCH_MU_PDICT = OrderedDict(
     [
-        ("mean_ulgm_mseq_xtp", 12.780),
-        ("mean_ulgm_mseq_ytp", 11.668),
-        ("mean_ulgm_mseq_lo", 0.748),
+        ("mean_ulgm_mseq_xtp", 12.029),
+        ("mean_ulgm_mseq_ytp", 11.300),
+        ("mean_ulgm_mseq_lo", 2.812),
         ("mean_ulgm_mseq_hi", 0.400),
-        ("mean_ulgy_mseq_int", -9.299),
-        ("mean_ulgy_mseq_slp", 0.409),
+        ("mean_ulgy_mseq_xtp", 12.344),
+        ("mean_ulgy_mseq_ytp", -9.296),
+        ("mean_ulgy_mseq_lo", 0.902),
+        ("mean_ulgy_mseq_hi", 0.303),
         ("mean_ul_mseq_int", -0.639),
         ("mean_ul_mseq_slp", 0.988),
-        ("mean_uh_mseq_int", -0.607),
-        ("mean_uh_mseq_slp", -0.312),
-        ("mean_ulgm_qseq_xtp", 13.252),
-        ("mean_ulgm_qseq_ytp", 11.971),
-        ("mean_ulgm_qseq_lo", 0.639),
-        ("mean_ulgm_qseq_hi", 0.208),
-        ("mean_ulgy_qseq_int", -9.243),
-        ("mean_ulgy_qseq_slp", 0.660),
-        ("mean_ul_qseq_int", -0.758),
-        ("mean_ul_qseq_slp", 0.421),
-        ("mean_uh_qseq_int", -0.661),
-        ("mean_uh_qseq_slp", -0.334),
-        ("mean_uqt_int", 1.065),
-        ("mean_uqt_slp", -0.003),
-        ("mean_uqs_int", -0.352),
-        ("mean_uqs_slp", -0.273),
-        ("mean_udrop_int", -1.938),
-        ("mean_udrop_slp", 0.551),
-        ("mean_urej_int", -0.620),
-        ("mean_urej_slp", -0.206),
+        ("mean_uh_mseq_int", -0.391),
+        ("mean_uh_mseq_slp", 0.161),
+        ("mean_ulgm_qseq_xtp", 13.078),
+        ("mean_ulgm_qseq_ytp", 11.903),
+        ("mean_ulgm_qseq_lo", 0.794),
+        ("mean_ulgm_qseq_hi", 0.230),
+        ("mean_ulgy_qseq_xtp", 12.112),
+        ("mean_ulgy_qseq_ytp", -9.421),
+        ("mean_ulgy_qseq_lo", 0.725),
+        ("mean_ulgy_qseq_hi", 0.356),
+        ("mean_ul_qseq_int", -0.714),
+        ("mean_ul_qseq_slp", 0.616),
+        ("mean_uh_qseq_int", -0.582),
+        ("mean_uh_qseq_slp", -0.115),
+        ("mean_uqt_xtp", 12.922),
+        ("mean_uqt_ytp", 1.047),
+        ("mean_uqt_lo", -0.100),
+        ("mean_uqt_hi", -0.100),
+        ("mean_uqs_int", -0.298),
+        ("mean_uqs_slp", -0.411),
+        ("mean_udrop_int", -1.953),
+        ("mean_udrop_slp", 0.586),
+        ("mean_urej_int", -0.853),
+        ("mean_urej_slp", 0.329),
     ]
 )
 
 SFH_PDF_QUENCH_COV_MS_BLOCK_PDICT = OrderedDict(
     [
-        ("std_ulgm_mseq_int", 0.355),
-        ("std_ulgm_mseq_slp", -0.015),
-        ("std_ulgy_mseq_int", 0.277),
-        ("std_ulgy_mseq_slp", 0.103),
-        ("std_ul_mseq_int", 2.221),
-        ("std_ul_mseq_slp", 0.900),
-        ("std_uh_mseq_int", 0.651),
-        ("std_uh_mseq_slp", -0.155),
-        ("std_ulgm_qseq_int", 0.313),
-        ("std_ulgm_qseq_slp", -0.022),
-        ("std_ulgy_qseq_int", 0.257),
-        ("std_ulgy_qseq_slp", 0.075),
-        ("std_ul_qseq_int", 2.241),
-        ("std_ul_qseq_slp", -0.130),
-        ("std_uh_qseq_int", 0.482),
-        ("std_uh_qseq_slp", 0.007),
+        ("std_ulgm_mseq_int", 0.552),
+        ("std_ulgm_mseq_slp", -0.433),
+        ("std_ulgy_mseq_int", 0.298),
+        ("std_ulgy_mseq_slp", 0.046),
+        ("std_ul_mseq_int", 2.469),
+        ("std_ul_mseq_slp", 0.677),
+        ("std_uh_mseq_int", 0.705),
+        ("std_uh_mseq_slp", -0.287),
+        ("std_ulgm_qseq_int", 0.492),
+        ("std_ulgm_qseq_slp", -0.293),
+        ("std_ulgy_qseq_int", 0.264),
+        ("std_ulgy_qseq_slp", 0.063),
+        ("std_ul_qseq_int", 2.263),
+        ("std_ul_qseq_slp", -0.165),
+        ("std_uh_qseq_int", 0.528),
+        ("std_uh_qseq_slp", -0.067),
     ]
 )
 
 SFH_PDF_QUENCH_COV_Q_BLOCK_PDICT = OrderedDict(
     [
-        ("std_uqt_int", 0.067),
-        ("std_uqt_slp", 0.013),
+        ("std_uqt_int", 0.075),
+        ("std_uqt_slp", -0.008),
         ("std_uqs_int", 0.900),
         ("std_uqs_slp", 0.900),
-        ("std_udrop_int", 0.589),
-        ("std_udrop_slp", 0.049),
-        ("std_urej_int", 0.867),
-        ("std_urej_slp", 0.049),
+        ("std_udrop_int", 0.600),
+        ("std_udrop_slp", 0.020),
+        ("std_urej_int", 0.979),
+        ("std_urej_slp", -0.205),
     ]
 )
 
 SFH_PDF_FRAC_QUENCH_PDICT = OrderedDict(
     [
         ("frac_quench_cen_x0_tpeak", 7.000),
         ("frac_quench_cen_k_tpeak", 2.000),
-        ("frac_quench_cen_x0_ylotpeak", 11.100),
-        ("frac_quench_cen_x0_yhitpeak", 13.008),
+        ("frac_quench_cen_x0_ylotpeak", 11.750),
+        ("frac_quench_cen_x0_yhitpeak", 12.965),
         ("frac_quench_cen_ylo_ylotpeak", 0.990),
-        ("frac_quench_cen_ylo_yhitpeak", 0.446),
+        ("frac_quench_cen_ylo_yhitpeak", 0.625),
         ("frac_quench_cen_k", 3.848),
         ("frac_quench_cen_yhi", 0.971),
         ("frac_quench_sat_x0_tpeak", 7.000),
         ("frac_quench_sat_k_tpeak", 2.000),
-        ("frac_quench_sat_x0_ylotpeak", 11.100),
-        ("frac_quench_sat_x0_yhitpeak", 13.008),
+        ("frac_quench_sat_x0_ylotpeak", 11.750),
+        ("frac_quench_sat_x0_yhitpeak", 12.965),
         ("frac_quench_sat_ylo_ylotpeak", 0.990),
-        ("frac_quench_sat_ylo_yhitpeak", 0.446),
+        ("frac_quench_sat_ylo_yhitpeak", 0.625),
         ("frac_quench_sat_k", 3.848),
         ("frac_quench_sat_yhi", 0.971),
     ]
 )
 DELTA_UQT_PDICT = OrderedDict(
     [
         ("delta_uqt_x0", 1.001),
-        ("delta_uqt_k", 0.836),
-        ("delta_uqt_ylo", -0.583),
-        ("delta_uqt_yhi", -0.006),
-        ("delta_uqt_slope", -0.017),
+        ("delta_uqt_k", 0.714),
+        ("delta_uqt_ylo", -0.611),
+        ("delta_uqt_yhi", 0.002),
+        ("delta_uqt_slope", -0.024),
     ]
 )
 SFH_PDF_QUENCH_PDICT = SFH_PDF_FRAC_QUENCH_PDICT.copy()
@@ -112,7 +120,6 @@
 _UPNAMES = ["u_" + key for key in QseqParams._fields]
 QseqUParams = namedtuple("QseqUParams", _UPNAMES)
 
-
 DIFFSTARFITS_GALACTICUS_INPLUSEX_DIFFSTARPOP_PARAMS = DiffstarPopParams(
     *SFH_PDF_QUENCH_PARAMS, *DEFAULT_SATQUENCHPOP_PARAMS
 )

diffstar/diffstarpop/kernels/params/params_diffstarfits_mgash_galacticus_in_situ.py

@@ -1,104 +1,112 @@
-""""""
-
 from collections import OrderedDict, namedtuple
 
-from ..defaults_mgash import DiffstarPopParams, get_unbounded_diffstarpop_params
+from ..defaults_mgash import (
+    DiffstarPopParams,
+    get_unbounded_diffstarpop_params,
+)
+
 from ..satquenchpop_model import DEFAULT_SATQUENCHPOP_PARAMS
 
 SFH_PDF_QUENCH_MU_PDICT = OrderedDict(
     [
-        ("mean_ulgm_mseq_xtp", 13.561),
-        ("mean_ulgm_mseq_ytp", 12.181),
-        ("mean_ulgm_mseq_lo", 0.628),
-        ("mean_ulgm_mseq_hi", -0.819),
-        ("mean_ulgy_mseq_int", -9.663),
-        ("mean_ulgy_mseq_slp", 0.105),
+        ("mean_ulgm_mseq_xtp", 13.189),
+        ("mean_ulgm_mseq_ytp", 12.079),
+        ("mean_ulgm_mseq_lo", 0.783),
+        ("mean_ulgm_mseq_hi", -0.001),
+        ("mean_ulgy_mseq_xtp", 11.995),
+        ("mean_ulgy_mseq_ytp", -9.781),
+        ("mean_ulgy_mseq_lo", 0.927),
+        ("mean_ulgy_mseq_hi", 0.091),
         ("mean_ul_mseq_int", -0.715),
         ("mean_ul_mseq_slp", 1.717),
-        ("mean_uh_mseq_int", -0.392),
-        ("mean_uh_mseq_slp", 0.043),
-        ("mean_ulgm_qseq_xtp", 12.728),
-        ("mean_ulgm_qseq_ytp", 11.931),
-        ("mean_ulgm_qseq_lo", 0.769),
-        ("mean_ulgm_qseq_hi", 0.123),
-        ("mean_ulgy_qseq_int", -9.647),
-        ("mean_ulgy_qseq_slp", 0.382),
-        ("mean_ul_qseq_int", -0.378),
-        ("mean_ul_qseq_slp", 1.047),
-        ("mean_uh_qseq_int", -0.196),
-        ("mean_uh_qseq_slp", 0.157),
-        ("mean_uqt_int", 1.034),
-        ("mean_uqt_slp", -0.096),
-        ("mean_uqs_int", 0.040),
-        ("mean_uqs_slp", -0.173),
-        ("mean_udrop_int", -1.963),
-        ("mean_udrop_slp", 0.746),
-        ("mean_urej_int", -0.524),
-        ("mean_urej_slp", 0.164),
+        ("mean_uh_mseq_int", -0.253),
+        ("mean_uh_mseq_slp", 0.364),
+        ("mean_ulgm_qseq_xtp", 13.473),
+        ("mean_ulgm_qseq_ytp", 12.089),
+        ("mean_ulgm_qseq_lo", 0.404),
+        ("mean_ulgm_qseq_hi", -0.228),
+        ("mean_ulgy_qseq_xtp", 12.303),
+        ("mean_ulgy_qseq_ytp", -9.608),
+        ("mean_ulgy_qseq_lo", 0.446),
+        ("mean_ulgy_qseq_hi", 0.017),
+        ("mean_ul_qseq_int", -0.388),
+        ("mean_ul_qseq_slp", 0.954),
+        ("mean_uh_qseq_int", -0.153),
+        ("mean_uh_qseq_slp", 0.284),
+        ("mean_uqt_xtp", 12.746),
+        ("mean_uqt_ytp", 1.035),
+        ("mean_uqt_lo", -0.100),
+        ("mean_uqt_hi", -0.150),
+        ("mean_uqs_int", 0.057),
+        ("mean_uqs_slp", -0.196),
+        ("mean_udrop_int", -1.917),
+        ("mean_udrop_slp", 0.684),
+        ("mean_urej_int", -0.749),
+        ("mean_urej_slp", 0.472),
     ]
 )
 
 SFH_PDF_QUENCH_COV_MS_BLOCK_PDICT = OrderedDict(
     [
-        ("std_ulgm_mseq_int", 0.405),
-        ("std_ulgm_mseq_slp", -0.094),
-        ("std_ulgy_mseq_int", 0.229),
-        ("std_ulgy_mseq_slp", -0.006),
-        ("std_ul_mseq_int", 2.073),
-        ("std_ul_mseq_slp", 0.374),
-        ("std_uh_mseq_int", 1.176),
-        ("std_uh_mseq_slp", -0.200),
-        ("std_ulgm_qseq_int", 0.420),
-        ("std_ulgm_qseq_slp", -0.019),
-        ("std_ulgy_qseq_int", 0.227),
-        ("std_ulgy_qseq_slp", -0.014),
-        ("std_ul_qseq_int", 2.174),
-        ("std_ul_qseq_slp", -0.137),
-        ("std_uh_qseq_int", 0.761),
-        ("std_uh_qseq_slp", -0.141),
+        ("std_ulgm_mseq_int", 0.527),
+        ("std_ulgm_mseq_slp", -0.248),
+        ("std_ulgy_mseq_int", 0.247),
+        ("std_ulgy_mseq_slp", -0.031),
+        ("std_ul_mseq_int", 2.156),
+        ("std_ul_mseq_slp", 0.260),
+        ("std_uh_mseq_int", 1.112),
+        ("std_uh_mseq_slp", -0.113),
+        ("std_ulgm_qseq_int", 0.434),
+        ("std_ulgm_qseq_slp", -0.097),
+        ("std_ulgy_qseq_int", 0.229),
+        ("std_ulgy_qseq_slp", -0.018),
+        ("std_ul_qseq_int", 2.179),
+        ("std_ul_qseq_slp", -0.124),
+        ("std_uh_qseq_int", 0.770),
+        ("std_uh_qseq_slp", -0.177),
     ]
 )
 
 SFH_PDF_QUENCH_COV_Q_BLOCK_PDICT = OrderedDict(
     [
-        ("std_uqt_int", 0.083),
-        ("std_uqt_slp", 0.070),
+        ("std_uqt_int", 0.110),
+        ("std_uqt_slp", 0.032),
         ("std_uqs_int", 0.900),
-        ("std_uqs_slp", 0.002),
-        ("std_udrop_int", 0.607),
-        ("std_udrop_slp", 0.097),
-        ("std_urej_int", 0.893),
-        ("std_urej_slp", 0.260),
+        ("std_uqs_slp", 0.120),
+        ("std_udrop_int", 0.639),
+        ("std_udrop_slp", 0.053),
+        ("std_urej_int", 1.118),
+        ("std_urej_slp", -0.044),
     ]
 )
 
 SFH_PDF_FRAC_QUENCH_PDICT = OrderedDict(
     [
         ("frac_quench_cen_x0_tpeak", 7.000),
         ("frac_quench_cen_k_tpeak", 2.000),
-        ("frac_quench_cen_x0_ylotpeak", 11.708),
-        ("frac_quench_cen_x0_yhitpeak", 13.900),
+        ("frac_quench_cen_x0_ylotpeak", 11.692),
+        ("frac_quench_cen_x0_yhitpeak", 12.963),
         ("frac_quench_cen_ylo_ylotpeak", 0.990),
-        ("frac_quench_cen_ylo_yhitpeak", 0.413),
+        ("frac_quench_cen_ylo_yhitpeak", 0.525),
         ("frac_quench_cen_k", 3.848),
         ("frac_quench_cen_yhi", 0.971),
         ("frac_quench_sat_x0_tpeak", 7.000),
         ("frac_quench_sat_k_tpeak", 2.000),
-        ("frac_quench_sat_x0_ylotpeak", 11.708),
-        ("frac_quench_sat_x0_yhitpeak", 13.900),
+        ("frac_quench_sat_x0_ylotpeak", 11.692),
+        ("frac_quench_sat_x0_yhitpeak", 12.963),
         ("frac_quench_sat_ylo_ylotpeak", 0.990),
-        ("frac_quench_sat_ylo_yhitpeak", 0.413),
+        ("frac_quench_sat_ylo_yhitpeak", 0.525),
         ("frac_quench_sat_k", 3.848),
         ("frac_quench_sat_yhi", 0.971),
     ]
 )
 DELTA_UQT_PDICT = OrderedDict(
     [
-        ("delta_uqt_x0", 1.977),
-        ("delta_uqt_k", 1.334),
-        ("delta_uqt_ylo", -0.367),
-        ("delta_uqt_yhi", -0.008),
-        ("delta_uqt_slope", -0.013),
+        ("delta_uqt_x0", 1.659),
+        ("delta_uqt_k", 0.939),
+        ("delta_uqt_ylo", -0.466),
+        ("delta_uqt_yhi", 0.000),
+        ("delta_uqt_slope", -0.022),
     ]
 )
 SFH_PDF_QUENCH_PDICT = SFH_PDF_FRAC_QUENCH_PDICT.copy()
@@ -112,7 +120,6 @@
 _UPNAMES = ["u_" + key for key in QseqParams._fields]
 QseqUParams = namedtuple("QseqUParams", _UPNAMES)
 
-
 DIFFSTARFITS_GALACTICUS_IN_DIFFSTARPOP_PARAMS = DiffstarPopParams(
     *SFH_PDF_QUENCH_PARAMS, *DEFAULT_SATQUENCHPOP_PARAMS
 )