Generating a weighted Monte Carlo lightcone of Diffsky galaxies

This notebook demonstrates how to generate a lightcone of diffsky galaxies with SEDs, star formation histories, dust, and other properties.

import numpy as np
from matplotlib import pyplot as plt

Download stellar population synthesis data

The very first thing we will do in this notebook is download some LSST-like transmission curves, and also some SEDs of simple stellar populations (SSPs) used in the SED modeling. You can skip this step of you have already cached previously-stored SSP data as described in the DSPS Quickstart Guide.

Note that for the demonstration purposes of these docs, we use a “sparse” set of SSP SEDs, but high-res SSPs should be used for more accurate SEDs.

import os
os.makedirs('./drn_dsps_temp', exist_ok=True)

filter_names = ('u', 'g', 'r', 'i', 'z', 'y')
base_url = "https://portal.nersc.gov/project/hacc/aphearin/DSPS_data/filters/lsst_{}_transmission.h5"

for filt in filter_names:
    ! wget -q -P ./drn_dsps_temp {base_url.format(filt)}

! wget -q -P ./drn_dsps_temp https://portal.nersc.gov/project/hacc/aphearin/DSPS_data/fsps_v0.4.7_mist_c3k_a_kroupa_wNE_logGasU-2.0_logGasZ0.0.sparse.h5

The load_ssp_templates function in DSPS loads an hdf5 file and packs it into a namedtuple with the expected field names.

from dsps.data_loaders import load_ssp_templates
ssp_data = load_ssp_templates(
    fn="drn_dsps_temp/fsps_v0.4.7_mist_c3k_a_kroupa_wNE_logGasU-2.0_logGasZ0.0.sparse.h5")

print(ssp_data._fields)

('ssp_lgmet', 'ssp_lg_age_gyr', 'ssp_wave', 'ssp_flux', 'ssp_emline_wave', 'ssp_emline_luminosity')

The load_transmission_curve function in DSPS loads an hdf5 file and packs it into a namedtuple with the expected field names.

from dsps.data_loaders import load_transmission_curve
fn_pat = "drn_dsps_temp/lsst_{}_transmission.h5"
filter_names = ('u', 'g', 'r', 'i', 'z', 'y')
fn_list = [fn_pat.format(x) for x in filter_names]
tcurves = [load_transmission_curve(fn) for fn in fn_list]

print(tcurves[0]._fields)

('wave', 'transmission')

Generate the halo lightcone and additional data needed to model/compute SEDs

First specify halo lightcone specs

from diffsky.experimental.lightcone_generators import weighted_lc_photdata

num_halos = 500
z_min, z_max = 0.1, 2.0
lgmp_min, lgmp_max = 10.5, 15.0
sky_area_degsq = 10.0

Get a random number seed from JAX

from jax import random as jran
ran_key = jran.key(0)

Define a redshift table used for photometry interpolation

n_z_phot_table = 15
z_phot_table = np.linspace(z_min, z_max, n_z_phot_table)

Generate the lightcone of dark matter halos and additional SPS data

halo_lc_data = (num_halos, z_min, z_max, lgmp_min, lgmp_max, sky_area_degsq)
phot_data = (ssp_data, tcurves, z_phot_table)

ran_key, lc_halo_key = jran.split(ran_key, 2)
args = (lc_halo_key, *halo_lc_data, *phot_data)

lc_data = weighted_lc_photdata(*args)

Populate the lightcone with diffsky galaxies

Now we will pass the lc_data to the diffsky SED kernels to populate the halo lightcone with galaxies.

from diffsky.experimental.mc_phot import mc_lc_phot

mc_merge = 0
ran_key, sed_key = jran.split(ran_key, 2)
phot_info, __, __ = mc_lc_phot(sed_key, lc_data, mc_merge)
print(phot_info._fields)

('obs_mags', 't_table', 'lgmcrit', 'lgy_at_mcrit', 'indx_lo', 'indx_hi', 'lg_qt', 'qlglgdt', 'lg_drop', 'lg_rejuv', 'sfh_table', 'logsm_obs', 'logssfr_obs', 'mc_sfh_type', 'ssp_weights', 'lgmet_weights', 'lgfburst', 'lgyr_peak', 'lgyr_max', 'av', 'delta', 'funo', 'dust_frac_trans', 'ssp_photflux_table', 'frac_ssp_errors', 'wave_eff_galpop', 'obs_mags_ms', 'obs_mags_q', 'obs_mags_bursty', 'frac_q', 'obs_mags_weighted', 'diffstar_info_ms', 'diffstar_info_q', 'burstiness_info_ms', 'burstiness_info_q', 'logsm_obs_in_situ', 'obs_mags_in_situ', 'p_merge', 'uran_pmerge')

Choosing an alternative set of diffsky model parameters

There are several choices for diffsky parameters that are based on ongoing work calibrating Diffsky to the COSMOS and FENIKS-UDS datasets:

COSMOS calibrated diffsky model parameters

Below we show how to check which COSMOS parameters are available, and then select the cosmos_260210 calibration to populate the lightcone.

from diffsky.param_utils.cosmos_calibrations import COSMOS_PARAM_FITS
print(COSMOS_PARAM_FITS.keys())

dict_keys(['cosmos_260105', 'cosmos_260120_UM', 'cosmos_260210', 'cosmos_260215'])

from diffsky.param_utils.get_mock_params import get_param_collection_for_mock
cosmos_param_collection = get_param_collection_for_mock(cosmos_fit='cosmos_260210')

Using cosmos_fit = cosmos_260210

FENIKS-UDS calibrated diffsky model parameters

Below we show how to load a FENIKS-UDS based calibration feniks_260617:

from diffsky.param_utils.get_calib_params import get_calib_params
feniks_param_collection = get_calib_params(
    calibration_dir="feniks_calibrations", calibration_name="feniks_260617"
)

We can generate the photometry using these alternative set of diffsky model parameters as below

phot_info_feniks, _, _ = mc_lc_phot(sed_key, lc_data, mc_merge, param_collection = feniks_param_collection)

For the weighted lightcone, each halo has multiplicity according to its abundance in the volume

This cen_weight and sat_weight columns need to be taken into account when predicting summary statistics from the weighted lightcone. The product of these two weights supplies the total weight of each galaxy.

gal_weight = lc_data.cen_weight * lc_data.sat_weight

fig, ax = plt.subplots(1, 1)
__=ax.loglog()
__=ax.scatter(10**lc_data.logmp_obs, gal_weight, s=1)
xlabel = ax.set_xlabel(r'$M_{\rm halo}\ [M_{\odot}]$')
ylabel = ax.set_ylabel(r'$N_{\rm halos}$')

Calculate the halo mass function, accounting for halo weights

The unweighted version uniformly spans \(\log_{10}M_{\rm halo}\). The weighted version has the expected Schechter-type shape of the HMF.

fig, ax = plt.subplots(1, 1)
yscale = ax.set_yscale('log')
__=ax.hist(lc_data.logmp_obs, bins=100, alpha=0.7, label=r'${\rm unweighted}$')
__=ax.hist(lc_data.logmp_obs, bins=100, weights=gal_weight, alpha=0.7, label=r'${\rm weighted}$')
xlabel = ax.set_xlabel(r'$\log_{10}M_{\rm halo}/M_{\odot}$')
ylabel = ax.set_ylabel(r'$N_{\rm halos}$')
leg = ax.legend()

Visually inspect the diversity of SFHs

fig, ax = plt.subplots(1, 1)
yscale = ax.set_yscale('log')
ylim = ax.set_ylim(8e-4, 5e2)
xscale = ax.set_xscale('log')
xlim = ax.set_xlim(1, 15)

n_plot = 5
for i in range(n_plot):
    __=ax.plot(lc_data.t_table, phot_info.sfh_table[i, :])

xlabel = ax.set_xlabel(r'${\rm cosmic\ time\ [Gyr]}$')
ylabel = ax.set_ylabel(r'${\rm SFR\ [M_{\odot}/yr]}$')

Visually inspect the sSFR PDF

Note that the plot below shows the PDF for all galaxies/halos in the lightcone, without accounting for the weights

fig, ax = plt.subplots(1, 1)
xlim = ax.set_xlim(-15, -7)
yscale = ax.set_yscale('log')
__=ax.hist(phot_info.logssfr_obs, bins=70, alpha=0.7, weights=gal_weight, density=True)
xlabel = ax.set_xlabel(r'${\rm log_{10}(sSFR)}$')

Visually inspect star-forming sequence

fig, ax = plt.subplots(1, 1)
ylim = ax.set_ylim(-13.5, -7.5)
xlim = ax.set_xlim(7.5, 13)
__=ax.scatter(phot_info.logsm_obs, phot_info.logssfr_obs, s=1)
xlabel = ax.set_xlabel(r'${\rm log_{10}(M_{\star})}$')
ylabel = ax.set_ylabel(r'${\rm log_{10}(sSFR)}$')

Plot color-color diagram

Note that the plot below shows the PDF for all galaxies/halos in the lightcone, without accounting for the weights, and without selecting a particular galaxy sample of interest.

fig, ax = plt.subplots(1, 1)

ri = phot_info.obs_mags[:, 2]-phot_info.obs_mags[:, 3]
iz = phot_info.obs_mags[:, 3]-phot_info.obs_mags[:, 4]
__=ax.scatter(iz, ri, s=1)
xlabel = ax.set_xlabel(r'${\rm i-z}$')
ylabel = ax.set_ylabel(r'${\rm r-i}$')

Plot color–redshift relation

Note that the plot below shows the PDF for all galaxies/halos in the lightcone, without accounting for the weights, and without selecting a particular galaxy sample of interest.

fig, ax = plt.subplots(1, 1)

iz = phot_info.obs_mags[:, 3]-phot_info.obs_mags[:, 4]
__=ax.scatter(lc_data.z_obs, iz, s=1)
xlabel = ax.set_xlabel(r'${\rm redshift}$')
ylabel = ax.set_ylabel(r'${\rm i-z}$')

Docs about accounting for the weights, computing high-res SEDs, and more coming soon…