Generating a 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.

Before generating these docs, these data were already been downloaded, and so we will just load them directly. See the DSPS docs for download instructions.

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

Note: This tutorial uses unrealistically sparse SSP data to minimize memory demands of generating the documentation. You should use full-resolution SSP SEDs in real applications.

from diffsky.data_loaders import load_ssp_templates
ssp_data = load_ssp_templates(
    drn="dsps_drn",
    bn="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.

You can use any transmission curves to compute Diffsky photometry. You just need to create a namedtuple storing ‘wave’, ‘transmission’ of each transmission curve.

from dsps.data_loaders import load_transmission_curve
bn_base_pat = "lsst_{}*"
filter_names = ('u', 'g', 'r', 'i', 'z', 'y')
bn_pat_list = [bn_base_pat.format(x) for x in filter_names]
tcurves = [load_transmission_curve(bn_pat=pat, drn="dsps_drn/filters") for pat in bn_pat_list]

print(tcurves[0]._fields)

('wave', 'transmission')

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

First specify halo lightcone specs:

• ran_key: Random number seed.

• z_min, z_max: the redshift range

• lgmp_min should be small enough to resolve the faintest galaxies in your population (values of 10-11 are typical for cosmological surveys).

• lgmsub_min is the lowest mass of the subhalos. In most cases, this should typically be set equal to lgmp_min.

• sky_area_degsq is the sky area in square degrees.

from diffsky.experimental.mc_lightcone_generators import mc_lc_photdata

z_min, z_max = 0.1, 0.15
lgmp_min, lgmp_sub_min = 11.0, 11.0
sky_area_degsq = 1.0

Note: This tutorial demonstrates how to use the mc_lc_photdata and mc_lc_phot functions, which is currently imported from within diffsky.experimental. Code located within diffsky.experimental may not have a stable API, and will be relocated elsewhere in the library after development of this feature settles down. This tutorial is built dynamically with each PR, and so the syntax demonstrated here is always up-to-date with the main branch of the code.

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

ran_key, lc_halo_key = jran.split(ran_key, 2)
args = (lc_halo_key,
    z_min,
    z_max,
    lgmp_min,
    lgmp_sub_min,
    sky_area_degsq,
    ssp_data,
    tcurves,
    z_phot_table,)

lc_data = mc_lc_photdata(*args)
print(lc_data._fields)
print(f"NUmber of halos = {lc_data.logmp_obs.size}")

('cen_weight', 'z_obs', 't_obs', 'logmp_obs', 'mah_params', 'logmp0', 't_table', 'ssp_data', 'precomputed_ssp_mag_table', 'z_phot_table', 'wave_eff_table', 'sat_weight', 't_infall', 'logmp_infall', 'logmhost_infall', 'is_central', 'halo_indx', 'halo_weight', 'precomputed_ssp_linelum_cgs_table', 'line_wave_table')
NUmber of halos = 499

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.

For Monte Carlo realizations, we set mc_merge = 1: satellite galaxies that are either merged with their associated central (p_merge=1, or they continue to exist as they orbit within the host halo p_merge=0. The photometry and stellar mass of merged satellites and their centrals has been self-consistently calculated by the mc_lc_phot function, and so you can treat merged satellites as if they do not exist by simply masking them out of your analysis.

from diffsky.experimental.mc_phot import mc_lc_phot

mc_merge = 1  # select descrete realization of merging PDF for Monte Carlo lightcones
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')

phot_info.ssp_weights.shape

(499, 6, 8)

print(np.unique(phot_info.p_merge))
print(np.mean(phot_info.p_merge>0))

[0.       0.999999]
0.3466934

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. Below we show how to check which parameters are available.

from diffsky.param_utils import DIFFSKY_FIT_PARAMS
print(DIFFSKY_FIT_PARAMS.keys())

dict_keys(['c260710', 'feniks_260617', 'sdss_feniks_260701', 'sdss_feniks_hizels_260710'])

cosmos_param_collection = DIFFSKY_FIT_PARAMS['c260710']
feniks_param_collection = DIFFSKY_FIT_PARAMS['sdss_feniks_hizels_260710']

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)

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 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.p_merge==0], phot_info.logssfr_obs[phot_info.p_merge==0], 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[phot_info.p_merge==0], ri[phot_info.p_merge==0], 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)
ylim = ax.set_ylim(-0.8, 0.8)
iz = phot_info.obs_mags[:, 3]-phot_info.obs_mags[:, 4]

msk_cens = lc_data.is_central==1
msk_sats = (lc_data.is_central==0) & (phot_info.p_merge==0)
__=ax.scatter(lc_data.z_obs[msk_cens], iz[msk_cens], s=1, label=r'${\rm centrals}$')
__=ax.scatter(lc_data.z_obs[msk_sats], iz[msk_sats], s=1, label=r'${\rm satellites}$')
xlabel = ax.set_xlabel(r'${\rm redshift}$')
ylabel = ax.set_ylabel(r'${\rm i-z}$')
leg = ax.legend(markerscale=4)