diff --git a/examples/creation_examples/plotting_interface_skysim_cosmoDC2_COSMOS2020_demo.ipynb b/examples/creation_examples/plotting_interface_skysim_cosmoDC2_COSMOS2020_demo.ipynb new file mode 100644 index 0000000..2983f6c --- /dev/null +++ b/examples/creation_examples/plotting_interface_skysim_cosmoDC2_COSMOS2020_demo.ipynb @@ -0,0 +1,1573 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "c087f952-7aad-4b5e-96f2-8fe4bb72243b", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import GCRCatalogs\n", + "import matplotlib.pyplot as plt\n", + "import pandas as pd\n", + "import corner\n", + "import numpy as np\n", + "import h5py\n", + "from astropy.table import Table\n", + "import os" + ] + }, + { + "cell_type": "markdown", + "id": "6b84d8ef-ea56-4bb8-b15f-5a5b10691d66", + "metadata": {}, + "source": [ + "## Read skysim v3.1.0 and cosmoDC2 v1.1.4 catalogs with the GCRCatalogs class" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3f367960-308e-49a9-8ba1-095c409b7915", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "skysim3 = GCRCatalogs.load_catalog('skysim_v3.1.0')\n", + "cosmodc2 = GCRCatalogs.load_catalog('cosmoDC2_v1.1.4_small')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a0675c78-2fff-4692-bfbb-de25a5a60e62", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "skysim3_allcols = skysim3.list_all_native_quantities()\n", + "cosmodc2_allcols = cosmodc2.list_all_native_quantities()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c12c1711-b8a7-4d60-9d94-4455b37a2490", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "skysim3_relevantcols = ['LSST_obs_g','LSST_obs_i','LSST_obs_r','LSST_obs_u','LSST_obs_y','LSST_obs_z','LSST_rest_G','LSST_rest_I','LSST_rest_R',\n", + " 'LSST_rest_U' ,'LSST_rest_Y','LSST_rest_Z','lg_met_mean','lg_met_scatter','log_sm','log_ssfr','obs_sfr','obs_sm',\n", + " 'redshift','sfr','HSC_obs_g','HSC_obs_i','HSC_obs_r','HSC_obs_y','HSC_obs_z','HSC_rest_G','HSC_rest_I','HSC_rest_R',\n", + " 'HSC_rest_Y','HSC_rest_Z', 'diffmah_early_index', 'lg_met_scatter', \n", + " 'diffstar_u_indx_lo', 'diffstar_u_tau_dep', 'diffmah_mah_logtc', 'diffmah_logmp_fit',\n", + " 'lg_met_mean', 'diffstar_u_qs', 'diffstar_u_q_rejuv', 'diffstar_u_lgmcrit',\n", + " 'diffstar_u_qt', 'diffmah_late_index', 'diffstar_u_q_drop', 'diffstar_u_lgy_at_mcrit',\n", + " 'diffstar_u_indx_hi']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f7f7ed74-1146-4094-b95a-3524e208bb49", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "cosmodc2_relevantcols = ['LSST_filters/magnitude:LSST_u:rest','LSST_filters/magnitude:LSST_g:rest','LSST_filters/magnitude:LSST_r:rest',\n", + " 'LSST_filters/magnitude:LSST_i:rest','LSST_filters/magnitude:LSST_z:rest','LSST_filters/magnitude:LSST_y:rest',\n", + " 'LSST_filters/magnitude:LSST_u:observed','LSST_filters/magnitude:LSST_g:observed','LSST_filters/magnitude:LSST_r:observed',\n", + " 'LSST_filters/magnitude:LSST_i:observed','LSST_filters/magnitude:LSST_z:observed','LSST_filters/magnitude:LSST_y:observed',\n", + " 'baseDC2/sfr','baseDC2/redshift','baseDC2/obs_sfr','baseDC2/obs_sm','baseDC2/sm']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f8769596-f16f-46f5-84ed-e88aa56f0ef0", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "skysim3_quantities = skysim3.get_quantities(skysim3_relevantcols, native_filters=['healpix_pixel == 9816'])\n", + "cosmodc2_quantities = cosmodc2.get_quantities(cosmodc2_relevantcols, native_filters=['healpix_pixel == 9816'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20db1f1b-d877-4134-9991-7b448d09a584", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "skysim3_df = pd.DataFrame(skysim3_quantities)\n", + "cosmodc2_df = pd.DataFrame(cosmodc2_quantities)" + ] + }, + { + "cell_type": "markdown", + "id": "cfdc79d5-b9d7-4b3e-9ae0-b7ebcac8838f", + "metadata": {}, + "source": [ + "## Select a random subset of 100k objects per catalog to facilitate plotting" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "934e777e-232e-4bc8-a59e-0331541aa074", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "skysim3_df_100k = skysim3_df.sample(n=100000)\n", + "cosmodc2_df_100k = cosmodc2_df.sample(n=100000)" + ] + }, + { + "cell_type": "markdown", + "id": "b7090af7-a65b-4bbe-ab17-37ea5219c112", + "metadata": { + "tags": [] + }, + "source": [ + "## Compare physical and photometric properties distribution between skysim v3.1.0 and cosmoDC2 v1.1.4" + ] + }, + { + "cell_type": "markdown", + "id": "13cf0157-1d43-4ab2-8735-0952eb5c5014", + "metadata": {}, + "source": [ + "#### Redshift distribution comparison" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1a4885dc-447f-43e1-8b6a-1f2fc3a68155", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "plt.clf()\n", + "plt.rcParams['figure.figsize']=(10,8)\n", + "plt.rcParams['xtick.labelsize']=25\n", + "plt.rcParams['ytick.labelsize']=25\n", + "plt.hist(cosmodc2_df_100k['baseDC2/redshift'], bins=50, lw=2, ls='solid', color='black', histtype='step',density=True, label='cosmoDC2 v1.1.4')\n", + "plt.hist(skysim3_df_100k['redshift'], bins=50, lw=2, ls='dashed', color='red', histtype='step',density=True, label='skysim v3.1.0')\n", + "plt.legend(loc='upper right', fontsize=20)\n", + "plt.xlim(0,3.5)\n", + "plt.xlabel(r'redshift', fontsize=20)\n", + "plt.ylabel(r'', fontsize=20)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "e2835ec7-5518-41a6-8a9e-0b9afa7164bc", + "metadata": {}, + "source": [ + "#### Stellar mass - Star formation rate plane comparison using values quoted in the catalogs. For skysim, these values refer to the stellar mass and star-formation rate evaluated at the last age bin of the galaxy SFH, i.e. 13.7 Gyr." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7c376cda-c68d-49a7-b2b1-575155ab1695", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "w = np.where((cosmodc2_df_100k['baseDC2/sm'].values>0)&(cosmodc2_df_100k['baseDC2/sfr'].values>0))\n", + "prop_cosmodc2 = np.empty((len(cosmodc2_df_100k['baseDC2/sm'].values[w]), 2))\n", + "prop_cosmodc2[:, 0] = np.log10(cosmodc2_df_100k['baseDC2/sm'].values[w])\n", + "prop_cosmodc2[:, 1] = np.log10(cosmodc2_df_100k['baseDC2/sfr'].values[w])\n", + "\n", + "prop_skysim = np.empty((len(skysim3_df_100k['log_sm']),2))\n", + "prop_skysim[:, 0] = skysim3_df_100k['log_sm'].values\n", + "prop_skysim[:, 1] = np.log10(skysim3_df_100k['sfr'].values)\n", + "\n", + "plt.clf()\n", + "plt.rcParams['figure.figsize']=(10,8)\n", + "plt.rcParams['xtick.labelsize']=15\n", + "plt.rcParams['ytick.labelsize']=15\n", + "\n", + "fig = corner.corner(prop_cosmodc2, hist_kwargs=dict(density=True, linewidth=1.5), #range=extents,\n", + " color='blue', labels=[r'$\\log(M/M_{\\odot})$',r'$\\log(SFR[M_{\\odot}/yr])$'],\n", + " label_kwargs={\"fontsize\":15}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "fig = corner.corner(prop_skysim, fig=fig, hist_kwargs=dict(density=True, linewidth=1.5), #range=extents,\n", + " color='red', labels=[r'$\\log(M/M_{\\odot})$',r'$\\log(SFR[M_{\\odot}/yr])$'],\n", + " label_kwargs={\"fontsize\":15}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "axes = np.array(fig.axes)\n", + "proxy = [plt.Rectangle((0,0),1,1,fc = 'blue'), \n", + " plt.Rectangle((0,0),1,1,fc = 'red')]\n", + "axes[1].legend(proxy, ['cosmoDC2 v1.1.4', 'skysim v3.1.0'],loc='upper right',fontsize=10)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "22d0d3f9-782f-4715-a524-9e450953472a", + "metadata": {}, + "source": [ + "#### Stellar mass - Star formation rate plane comparison using *_obs values, which refers to the UniverseMachine values used in galsampling. We refer to them as _um in the plot." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "523307f6-052a-4f05-a907-7e6cac8bdf18", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "w = np.where((cosmodc2_df_100k['baseDC2/obs_sm'].values>0)&(cosmodc2_df_100k['baseDC2/obs_sfr'].values>0))\n", + "prop_cosmodc2 = np.empty((len(cosmodc2_df_100k['baseDC2/obs_sm'].values[w]), 2))\n", + "prop_cosmodc2[:, 0] = np.log10(cosmodc2_df_100k['baseDC2/obs_sm'].values[w])\n", + "prop_cosmodc2[:, 1] = np.log10(cosmodc2_df_100k['baseDC2/obs_sfr'].values[w])\n", + "\n", + "prop_skysim = np.empty((len(skysim3_df_100k['obs_sm']),2))\n", + "prop_skysim[:, 0] = np.log10(skysim3_df_100k['obs_sm'].values)\n", + "prop_skysim[:, 1] = np.log10(skysim3_df_100k['obs_sfr'].values)\n", + "\n", + "plt.clf()\n", + "plt.rcParams['figure.figsize']=(10,8)\n", + "plt.rcParams['xtick.labelsize']=15\n", + "plt.rcParams['ytick.labelsize']=15\n", + "\n", + "fig = corner.corner(prop_cosmodc2, hist_kwargs=dict(density=True, linewidth=1.5), #range=extents,\n", + " color='blue', labels=[r'$\\log(M/M_{\\odot})_{um}$',r'$\\log(SFR[M_{\\odot}/yr])_{um}$'],\n", + " label_kwargs={\"fontsize\":15}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "fig = corner.corner(prop_skysim, fig=fig, hist_kwargs=dict(density=True, linewidth=1.5), #range=extents,\n", + " color='red', labels=[r'$\\log(M/M_{\\odot})_{um}$',r'$\\log(SFR[M_{\\odot}/yr])_{um}$'],\n", + " label_kwargs={\"fontsize\":15}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "axes = np.array(fig.axes)\n", + "proxy = [plt.Rectangle((0,0),1,1,fc = 'blue'), \n", + " plt.Rectangle((0,0),1,1,fc = 'red')]\n", + "axes[1].legend(proxy, ['cosmoDC2 v1.1.4', 'skysim v3.1.0'],loc='upper right',fontsize=10)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "b42fa42d-1d84-4d60-86e3-709ff0469b05", + "metadata": {}, + "source": [ + "#### Stellar mass - Star formation rate comparison where the quantities are computed for skysim by building the SFH of each galaxy and evaluating them at the time of observation corresponding to the galaxy redshift. This plot uses the SFHs computed with diffmah and diffstar, following the example script kindly provided by Andrew Hearin." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3a577536-2449-483c-b12c-6736a645dfae", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "from diffstar.defaults import DEFAULT_N_STEPS, INDX_K, LGT0, FB, T_BIRTH_MIN\n", + "from dsps.constants import T_TABLE_MIN\n", + "from diffstar.sfh import get_sfh_from_mah_kern\n", + "from dsps.cosmology import age_at_z, DEFAULT_COSMOLOGY\n", + "from dsps.utils import cumulative_mstar_formed\n", + "from jax import vmap\n", + "from jax import jit as jjit\n", + "from jax import numpy as jnp" + ] + }, + { + "cell_type": "markdown", + "id": "532fd8e9-a1db-4ecc-809e-359bccad94a4", + "metadata": {}, + "source": [ + "First we get the relevant parameters to build the SFH from the skysim catalog" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3f5b4b52-07e9-456c-877c-69107ccb1b41", + "metadata": {}, + "outputs": [], + "source": [ + "def get_fit_params(data, mah_keys=None, ms_keys=None, q_keys=None):\n", + " \"\"\"Read the mock galaxy data table and return the diffmah and diffstar fit params\n", + " Parameters\n", + " ----------\n", + " data : astropy Table of length (n_h, )\n", + " Synthetic galaxy catalog\n", + " mah_keys:\n", + " ms_keys:\n", + " q_keys:\n", + " Returns\n", + " -------\n", + " mah_params : ndarray of shape (n_h, 4)\n", + " ms_u_params : ndarray of shape (n_h, 5)\n", + " q_u_params : ndarray of shape (n_h, 4)\n", + " \"\"\"\n", + " if q_keys is None:\n", + " q_keys = DIFFSTAR_U_Q_KEYS\n", + " if ms_keys is None:\n", + " ms_keys = DIFFSTAR_U_MS_KEYS\n", + " if mah_keys is None:\n", + " mah_keys = DIFFMAH_KEYS\n", + " mah_params = np.array([data[key] for key in mah_keys]).T\n", + " u_ms_params = np.array([data[key] for key in ms_keys]).T\n", + " u_q_params = np.array([data[key] for key in q_keys]).T\n", + " return mah_params, u_ms_params, u_q_params\n", + "\n", + "DIFFMAH_KEYS = [\"diffmah_logmp_fit\", \"diffmah_mah_logtc\", \"diffmah_early_index\", \"diffmah_late_index\"]\n", + "_ms_keylist = (\"lgmcrit\", \"lgy_at_mcrit\", \"indx_lo\", \"indx_hi\", \"tau_dep\")\n", + "DIFFSTAR_U_MS_KEYS = [\"diffstar_u_\" + key for key in _ms_keylist]\n", + "DIFFSTAR_U_Q_KEYS = [\"diffstar_u_\" + key for key in (\"qt\", \"qs\", \"q_drop\", \"q_rejuv\")]\n", + "\n", + "mah_params, ms_params, q_params = get_fit_params(skysim3_df_100k,mah_keys=DIFFMAH_KEYS,\n", + " ms_keys=DIFFSTAR_U_MS_KEYS,\n", + " q_keys=DIFFSTAR_U_Q_KEYS)" + ] + }, + { + "cell_type": "markdown", + "id": "7ddcab45-6ea2-45bb-8684-29d817f2dde2", + "metadata": {}, + "source": [ + "Building the SFH using the kernel in diffstar" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d97eb2f1-a0fb-4ba5-9bc5-51f63ed5f35b", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "t_table = np.linspace(T_TABLE_MIN, 10**LGT0, DEFAULT_N_STEPS)\n", + "sfh_from_mah_kern = get_sfh_from_mah_kern(n_steps=DEFAULT_N_STEPS,\n", + " tacc_integration_min=T_BIRTH_MIN,\n", + " tobs_loop='vmap', galpop_loop='vmap')\n", + "sfh = sfh_from_mah_kern(t_table, mah_params, ms_params, q_params, LGT0, FB)" + ] + }, + { + "cell_type": "markdown", + "id": "b49fdf95-0410-42e4-bafe-81dad89b9d27", + "metadata": {}, + "source": [ + "Compute the values of the SFR at the age that the galaxy has when observed at redshift z" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "148eb3a2-877d-4520-bb29-9cfa2a93da24", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "t_obs = age_at_z(skysim3_df_100k['redshift'].values, *DEFAULT_COSMOLOGY)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "31e6a01b-753f-4be2-8e72-3b68a02102f3", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "def multiInterp2(x, xp, fp):\n", + " i = np.arange(x.size)\n", + " j = np.searchsorted(xp, x) - 1\n", + " d = (x - xp[j]) / (xp[j + 1] - xp[j])\n", + " return (1 - d) * fp[i, j] + fp[i, j + 1] * d\n", + "\n", + "sfr = multiInterp2(t_obs, t_table, sfh)" + ] + }, + { + "cell_type": "markdown", + "id": "121fbe0d-530c-49c2-8820-712cff4a1a19", + "metadata": { + "tags": [] + }, + "source": [ + "Compute the cumulative formed stellar mass up to t_obs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8fc4916c-060a-413c-a211-0d28d08fb8de", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "_c = [None, 0]\n", + "_calc_smh_table_vmap = jjit(vmap(cumulative_mstar_formed, in_axes=_c))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5ccc9d4b-7d0b-4392-954c-03c5769c8a43", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "args_smh_table = (t_table,sfh)\n", + "smh_table = _calc_smh_table_vmap(*args_smh_table)\n", + "logsmh_table = jnp.log10(smh_table)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "00a084e9-4223-4f23-a71f-48669d57ba8b", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "log_sm = multiInterp2(t_obs, t_table, logsmh_table)" + ] + }, + { + "cell_type": "markdown", + "id": "0f16ad1a-f067-4def-a85b-c5cf14477bf9", + "metadata": {}, + "source": [ + "Plot the Stellar Mass - Star formation rate plane" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0afb8b80-494f-4167-9602-5f337d19d5f1", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "w = np.where((cosmodc2_df_100k['baseDC2/sm'].values>0)&(cosmodc2_df_100k['baseDC2/sfr'].values>0))\n", + "prop_cosmodc2 = np.empty((len(cosmodc2_df_100k['baseDC2/sm'].values[w]), 2))\n", + "prop_cosmodc2[:, 0] = np.log10(cosmodc2_df_100k['baseDC2/sm'].values[w])\n", + "prop_cosmodc2[:, 1] = np.log10(cosmodc2_df_100k['baseDC2/sfr'].values[w])\n", + "\n", + "prop_skysim = np.empty((len(log_sm),2))\n", + "prop_skysim[:, 0] = log_sm\n", + "prop_skysim[:, 1] = np.log10(sfr)\n", + "\n", + "plt.clf()\n", + "plt.rcParams['figure.figsize']=(10,8)\n", + "plt.rcParams['xtick.labelsize']=15\n", + "plt.rcParams['ytick.labelsize']=15\n", + "\n", + "fig = corner.corner(prop_cosmodc2, hist_kwargs=dict(density=True, linewidth=1.5), #range=extents,\n", + " color='blue', labels=[r'$\\log(M/M_{\\odot})$',r'$\\log(SFR[M_{\\odot}/yr])$'],\n", + " label_kwargs={\"fontsize\":15}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "fig = corner.corner(prop_skysim, fig=fig, hist_kwargs=dict(density=True, linewidth=1.5), #range=extents,\n", + " color='red', labels=[r'$\\log(M/M_{\\odot})$',r'$\\log(SFR[M_{\\odot}/yr])$'],\n", + " label_kwargs={\"fontsize\":15}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "axes = np.array(fig.axes)\n", + "proxy = [plt.Rectangle((0,0),1,1,fc = 'blue'), \n", + " plt.Rectangle((0,0),1,1,fc = 'red')]\n", + "axes[1].legend(proxy, ['cosmoDC2 v1.1.4', 'skysim v3.1.0'],loc='upper right',fontsize=10)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "7d969daa-7cc3-48bc-be32-3e29cbb96374", + "metadata": {}, + "source": [ + "#### LSST rest-frame magnitudes and colours comparison" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7476e411-e0cc-4477-8f60-9f5f2cefefc0", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "wavebands = ['u', 'g', 'r', 'i', 'z', 'y']\n", + "lsst_restframe_mags_cosmodc2 = np.array([cosmodc2_df_100k['LSST_filters/magnitude:LSST_{}:rest'.format(waveband)].values for waveband in wavebands])\n", + "lsst_restframe_mags_skysim = np.array([skysim3_df_100k['LSST_rest_{}'.format(waveband.capitalize())].values for waveband in wavebands])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f4632b1b-cd41-40bd-895b-d8c8a06a4cf2", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "plt.clf()\n", + "plt.rcParams['figure.figsize']=(10,8)\n", + "plt.rcParams['xtick.labelsize']=15\n", + "plt.rcParams['ytick.labelsize']=15\n", + "\n", + "fig = corner.corner(np.swapaxes(lsst_restframe_mags_cosmodc2,0,1), hist_kwargs=dict(density=True, linewidth=1.5), #range=extents,\n", + " color='blue', labels=[r'$u_{rest}$',r'$g_{rest}$',r'$r_{rest}$',r'$i_{rest}$',r'$z_{rest}$',r'$y_{rest}$'],\n", + " label_kwargs={\"fontsize\":30}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "fig = corner.corner(np.swapaxes(lsst_restframe_mags_skysim,0,1), fig=fig, hist_kwargs=dict(density=True, linewidth=1.5), #range=extents,\n", + " color='red', labels=[r'$u_{rest}$',r'$g_{rest}$',r'$r_{rest}$',r'$i_{rest}$',r'$z_{rest}$',r'$y_{rest}$'],\n", + " label_kwargs={\"fontsize\":30}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "axes = np.array(fig.axes)\n", + "proxy = [plt.Rectangle((0,0),1,1,fc = 'blue'), \n", + " plt.Rectangle((0,0),1,1,fc = 'red')]\n", + "axes[1].legend(proxy, ['cosmoDC2 v1.1.4', 'skysim v3.1.0'],loc='upper right',fontsize=10)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e54f64d0-4ca6-4ee1-be99-7a195448f051", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "lsst_restframe_colours_cosmodc2 = np.array([cosmodc2_df_100k['LSST_filters/magnitude:LSST_{}:rest'.format(wavebands[i])].values - cosmodc2_df_100k['LSST_filters/magnitude:LSST_{}:rest'.format(wavebands[i+1])].values for i in range(len(wavebands)-1)])\n", + "lsst_restframe_colours_skysim = np.array([skysim3_df_100k['LSST_rest_{}'.format(wavebands[i].capitalize())].values-skysim3_df_100k['LSST_rest_{}'.format(wavebands[i+1].capitalize())].values for i in range(len(wavebands)-1)])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ad3bae02-f15b-4a06-bf32-3aead7ffa312", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "plt.clf()\n", + "plt.rcParams['figure.figsize']=(10,8)\n", + "plt.rcParams['xtick.labelsize']=15\n", + "plt.rcParams['ytick.labelsize']=15\n", + "\n", + "fig = corner.corner(np.swapaxes(lsst_restframe_colours_cosmodc2,0,1), hist_kwargs=dict(density=True, linewidth=1.5), #range=extents,\n", + " color='blue', labels=[r'$(u-g)_{rest}$',r'$(g-r)_{rest}$',r'$(r-i)_{rest}$',r'$(i-z)_{rest}$',r'$(z-y)_{rest}$'],\n", + " label_kwargs={\"fontsize\":30}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "fig = corner.corner(np.swapaxes(lsst_restframe_colours_skysim,0,1), fig=fig, hist_kwargs=dict(density=True, linewidth=1.5), #range=extents,\n", + " color='red', labels=[r'$(u-g)_{rest}$',r'$(g-r)_{rest}$',r'$(r-i)_{rest}$',r'$(i-z)_{rest}$',r'$(z-y)_{rest}$'],\n", + " label_kwargs={\"fontsize\":30}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "axes = np.array(fig.axes)\n", + "proxy = [plt.Rectangle((0,0),1,1,fc = 'blue'), \n", + " plt.Rectangle((0,0),1,1,fc = 'red')]\n", + "axes[1].legend(proxy, ['cosmoDC2 v1.1.4', 'skysim v3.1.0'],loc='upper right',fontsize=10)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "8ab0af22-0a4f-408a-b42f-34e5901bddf6", + "metadata": {}, + "source": [ + "#### LSST observed-frame magnitudes and colours comparison" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7e8e524c-9013-41a9-ab25-cc6efd7dc24c", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "wavebands = ['u', 'g', 'r', 'i', 'z', 'y']\n", + "lsst_obsframe_mags_cosmodc2 = np.array([cosmodc2_df_100k['LSST_filters/magnitude:LSST_{}:observed'.format(waveband)].values for waveband in wavebands])\n", + "lsst_obsframe_mags_skysim = np.array([skysim3_df_100k['LSST_obs_{}'.format(waveband)].values for waveband in wavebands])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9b8a9b12-86e4-43c6-bf3a-f3b018a891dd", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "plt.clf()\n", + "plt.rcParams['figure.figsize']=(10,8)\n", + "plt.rcParams['xtick.labelsize']=15\n", + "plt.rcParams['ytick.labelsize']=15\n", + "\n", + "fig = corner.corner(np.swapaxes(lsst_obsframe_mags_cosmodc2,0,1), hist_kwargs=dict(density=True, linewidth=1.5), #range=extents,\n", + " color='blue', labels=[r'$u_{obs}$',r'$g_{obs}$',r'$r_{obs}$',r'$i_{obs}$',r'$z_{obs}$',r'$y_{obs}$'],\n", + " label_kwargs={\"fontsize\":30}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "fig = corner.corner(np.swapaxes(lsst_obsframe_mags_skysim,0,1), fig=fig, hist_kwargs=dict(density=True, linewidth=1.5), #range=extents,\n", + " color='red', labels=[r'$u_{obs}$',r'$g_{obs}$',r'$r_{obs}$',r'$i_{obs}$',r'$z_{obs}$',r'$y_{obs}$'],\n", + " label_kwargs={\"fontsize\":30}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "axes = np.array(fig.axes)\n", + "proxy = [plt.Rectangle((0,0),1,1,fc = 'blue'), \n", + " plt.Rectangle((0,0),1,1,fc = 'red')]\n", + "axes[1].legend(proxy, ['cosmoDC2 v1.1.4', 'skysim v3.1.0'],loc='upper right',fontsize=10)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8271ee48-6c22-4eac-b94d-2441a2f9a5c7", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "lsst_obsframe_colours_cosmodc2 = np.array([cosmodc2_df_100k['LSST_filters/magnitude:LSST_{}:observed'.format(wavebands[i])].values - cosmodc2_df_100k['LSST_filters/magnitude:LSST_{}:observed'.format(wavebands[i+1])].values for i in range(len(wavebands)-1)])\n", + "lsst_obsframe_colours_skysim = np.array([skysim3_df_100k['LSST_obs_{}'.format(wavebands[i])].values-skysim3_df_100k['LSST_obs_{}'.format(wavebands[i+1])].values for i in range(len(wavebands)-1)])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b95a43aa-0f6f-4e5b-85e3-133fa7689596", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "plt.clf()\n", + "plt.rcParams['figure.figsize']=(10,8)\n", + "plt.rcParams['xtick.labelsize']=15\n", + "plt.rcParams['ytick.labelsize']=15\n", + "\n", + "fig = corner.corner(np.swapaxes(lsst_obsframe_colours_cosmodc2,0,1), hist_kwargs=dict(density=True, linewidth=1.5), #range=extents,\n", + " color='blue', labels=[r'$(u-g)_{obs}$',r'$(g-r)_{obs}$',r'$(r-i)_{obs}$',r'$(i-z)_{obs}$',r'$(z-y)_{obs}$'],\n", + " label_kwargs={\"fontsize\":30}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "fig = corner.corner(np.swapaxes(lsst_obsframe_colours_skysim,0,1), fig=fig, hist_kwargs=dict(density=True, linewidth=1.5), #range=extents,\n", + " color='red', labels=[r'$(u-g)_{obs}$',r'$(g-r)_{obs}$',r'$(r-i)_{obs}$',r'$(i-z)_{obs}$',r'$(z-y)_{obs}$'],\n", + " label_kwargs={\"fontsize\":30}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "axes = np.array(fig.axes)\n", + "proxy = [plt.Rectangle((0,0),1,1,fc = 'blue'), \n", + " plt.Rectangle((0,0),1,1,fc = 'red')]\n", + "axes[1].legend(proxy, ['cosmoDC2 v1.1.4', 'skysim v3.1.0'],loc='upper right',fontsize=10)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "8cb15051-18ad-4f16-bffc-0e92d66e9069", + "metadata": { + "tags": [] + }, + "source": [ + "#### Observed color-redshift relation comparison" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "46198d27-0dc1-43e0-a454-6533c08f1202", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "colours = [r'$(u-g)_{obs}$',r'$(g-r)_{obs}$',r'$(r-i)_{obs}$',r'$(i-z)_{obs}$',r'$(z-y)_{obs}$']\n", + "for i in range(len(colours)):\n", + " plt.clf()\n", + " plt.rcParams['figure.figsize']=(10,8)\n", + " plt.rcParams['xtick.labelsize']=25\n", + " plt.rcParams['ytick.labelsize']=25\n", + " plt.scatter(cosmodc2_df_100k['baseDC2/redshift'], lsst_obsframe_colours_cosmodc2[i,:], s=1, color='blue', label='cosmoDC2 v1.1.4')\n", + " plt.scatter(skysim3_df_100k['redshift'], lsst_obsframe_colours_skysim[i,:], s=1, color='red', label='skysim v3.1.0')\n", + " plt.legend(loc='upper right', fontsize=20)\n", + " plt.xlim(0,3.5)\n", + " plt.ylim(-1.5,4.5)\n", + " plt.ylabel(colours[i], fontsize=30)\n", + " plt.xlabel(r'redshift', fontsize=20)\n", + " plt.tight_layout()\n", + " plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "6622b0cb-8eac-4b87-bef8-065640d39cdf", + "metadata": {}, + "source": [ + "## Load COSMOS2020 catalog" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e05937ce-1a04-4901-816a-a467712cd35e", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "output_directory = os.getcwd()\n", + "\n", + "if not os.path.isfile(os.path.join(output_directory,'COSMOS2020_CLASSIC_R1_v2.0.fits')):\n", + " os.system('curl -O https://portal.nersc.gov/cfs/lsst/PZ/COSMOS2020_CLASSIC_R1_v2.0.fits '\n", + " '--output-dir {}'.format(output_directory))\n", + "\n", + "cosmos2020 = Table.read(os.path.join(output_directory,'COSMOS2020_CLASSIC_R1_v2.0.fits'), format='fits')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "30624b69-b012-4f3c-994e-03f2fb0e1352", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "mask_bad_values = np.where((cosmos2020['CFHT_u_MAG_AUTO'] > 0) & \n", + " (cosmos2020['HSC_i_MAG_AUTO'] > 0)&\n", + " (cosmos2020['HSC_g_MAG_AUTO'] > 0)&\n", + " (cosmos2020['HSC_r_MAG_AUTO'] > 0) & \n", + " (cosmos2020['HSC_z_MAG_AUTO'] > 0)&\n", + " (cosmos2020['HSC_y_MAG_AUTO'] > 0)&\n", + " (cosmos2020['UVISTA_J_MAG_AUTO'] > 0) &\n", + " (cosmos2020['UVISTA_H_MAG_AUTO'] > 0) & \n", + " (cosmos2020['UVISTA_Ks_MAG_AUTO'] > 0))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e4dd14a6-585c-45bf-9522-6454ef335f6b", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "cosmos2020 = cosmos2020[mask_bad_values]" + ] + }, + { + "cell_type": "markdown", + "id": "aea7d6bc-40b0-40c8-b604-0e81e2155a0b", + "metadata": {}, + "source": [ + "## Since skysim and cosmoDC2 catalogs contain model galaxies with noiseless properties, we compare them against COSMOS2020 observed galaxies by applying the LSST Y1 gold sample cut of 24.1 in the i-band (see DESC SRM v1.0.2, page 48). We use the HSC i-band in COSMOS since the LSST one is not available." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a34df6e0-3ecc-413b-a53f-1a8496b038fc", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "skysim3_df_lsst_y1_all = skysim3_df.loc[skysim3_df['LSST_obs_i']<=24.1]\n", + "cosmodc2_df_lsst_y1_all = cosmodc2_df.loc[cosmodc2_df['LSST_filters/magnitude:LSST_i:observed']<=24.1]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "42cc8668-d57d-48cb-9e65-eef4a663c744", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "mask_cosmos_lsst_y1 = np.where(cosmos2020['HSC_i_MAG_AUTO']<=24.1)\n", + "cosmos2020_lsst_y1_all = cosmos2020[mask_cosmos_lsst_y1]\n", + "rand_idx = np.random.randint(low=0, high=len(cosmos2020_lsst_y1_all),size=100000)\n", + "cosmos2020_lsst_y1 = cosmos2020_lsst_y1_all[rand_idx]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "200713a3-f6fd-4c5f-a46a-6f1d3ed92bc6", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "skysim3_df_lsst_y1 = skysim3_df_lsst_y1_all.sample(n=100000)\n", + "cosmodc2_df_lsst_y1 = cosmodc2_df_lsst_y1_all.sample(n=100000)" + ] + }, + { + "cell_type": "markdown", + "id": "a2e2d8bf-cd2c-45b2-97f0-80cde932eb91", + "metadata": {}, + "source": [ + "## Compare physical properties distribution among skysim v3.1.0, cosmoDC2 v1.1.4 and COSMOS2020, and photometric properties distribution between Skysim and COSMOS2020" + ] + }, + { + "cell_type": "markdown", + "id": "09d0a9b8-2521-4081-8cbd-76edd409b0d7", + "metadata": {}, + "source": [ + "#### Redshift distribution comparison between cosmoDC2 v1.1.4 and skysim v3.1.0 input model redshifts and COSMOS2020 observed redshifts obtained via Eazy photo-z (maximum a-posteriori)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "77f497c4-50b7-4a4f-8d13-c575b4ffc53d", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "plt.clf()\n", + "plt.rcParams['figure.figsize']=(10,8)\n", + "plt.rcParams['xtick.labelsize']=25\n", + "plt.rcParams['ytick.labelsize']=25\n", + "w = np.where((cosmos2020_lsst_y1['ez_z_phot']>0)&(cosmos2020_lsst_y1['ez_z_phot']<20))\n", + "plt.hist(cosmos2020_lsst_y1['ez_z_phot'][w], bins=50, lw=3, ls='solid', color='black', histtype='step',density=True, label='COSMOS2020 (Eazy)')\n", + "plt.hist(cosmodc2_df_lsst_y1['baseDC2/redshift'], bins=50, lw=2, ls='solid', color='blue', histtype='step',density=True, label='cosmoDC2 v1.1.4')\n", + "plt.hist(skysim3_df_lsst_y1['redshift'], bins=50, lw=2, ls='solid', color='red', histtype='step',density=True, label='skysim v3.1.0')\n", + "plt.vlines(np.nanmedian(cosmos2020_lsst_y1['ez_z_phot'][w]),0,1.2, lw=2, ls='dashed', color='black', label=r'$Me(\\bar{z}_{COSMOS2020(Eazy)})=$'+'{:.3f}'.format(np.nanmedian(cosmos2020_lsst_y1['ez_z_phot'][w])))\n", + "plt.vlines(np.nanmedian(cosmodc2_df_lsst_y1['baseDC2/redshift']),0,1.2, lw=2, ls='dashed', color='blue', label=r'$Me(\\bar{z}_{cosmoDC2})=$'+'{:.3f}'.format(np.nanmedian(cosmodc2_df_lsst_y1['baseDC2/redshift'])))\n", + "plt.vlines(np.nanmedian(skysim3_df_lsst_y1['redshift']),0,1.2, lw=2, ls='dashed', color='red', label=r'$Me(\\bar{z}_{skysim})=$'+'{:.3f}'.format(np.nanmedian(skysim3_df_lsst_y1['redshift'])))\n", + "plt.legend(loc='upper right', fontsize=10)\n", + "plt.xlim(0,3.5)\n", + "plt.xlabel(r'redshift', fontsize=20)\n", + "plt.ylabel(r'', fontsize=20)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "d7316a1e-def7-4554-ace7-cfa331d06e0c", + "metadata": {}, + "source": [ + "#### Redshift distribution comparison between cosmoDC2 v1.1.4 and skysim v3.1.0 input model redshifts and COSMOS2020 observed redshifts obtained via LePhare photo-z (median of the likelihood distribution). " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b8d4ffc4-cd78-4e2a-9fac-5e21dfaf3dd4", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "plt.clf()\n", + "plt.rcParams['figure.figsize']=(10,8)\n", + "plt.rcParams['xtick.labelsize']=25\n", + "plt.rcParams['ytick.labelsize']=25\n", + "w = np.where((cosmos2020_lsst_y1['lp_zPDF']>0)&(cosmos2020_lsst_y1['lp_zPDF']<20))\n", + "plt.hist(cosmos2020_lsst_y1['lp_zPDF'][w], bins=50, lw=3, ls='solid', color='black', histtype='step',density=True, label='COSMOS2020 (LePhare)')\n", + "plt.hist(cosmodc2_df_lsst_y1['baseDC2/redshift'], bins=50, lw=2, ls='solid', color='blue', histtype='step',density=True, label='cosmoDC2 v1.1.4')\n", + "plt.hist(skysim3_df_lsst_y1['redshift'], bins=50, lw=2, ls='solid', color='red', histtype='step',density=True, label='skysim v3.1.0')\n", + "plt.vlines(np.nanmedian(cosmos2020_lsst_y1['lp_zPDF'][w]),0,1.2, lw=2, ls='dashed', color='black', label=r'$Me(\\bar{z}_{COSMOS2020(LePhare)})=$'+'{:.3f}'.format(np.nanmedian(cosmos2020_lsst_y1['lp_zPDF'][w])))\n", + "plt.vlines(np.nanmedian(cosmodc2_df_lsst_y1['baseDC2/redshift']),0,1.2, lw=2, ls='dashed', color='blue', label=r'$Me(\\bar{z}_{cosmoDC2})=$'+'{:.3f}'.format(np.nanmedian(cosmodc2_df_lsst_y1['baseDC2/redshift'])))\n", + "plt.vlines(np.nanmedian(skysim3_df_lsst_y1['redshift']),0,1.2, lw=2, ls='dashed', color='red', label=r'$Me(\\bar{z}_{skysim})=$'+'{:.3f}'.format(np.nanmedian(skysim3_df_lsst_y1['redshift'])))\n", + "plt.legend(loc='upper right', fontsize=10)\n", + "plt.xlim(0,3.5)\n", + "plt.xlabel(r'redshift', fontsize=20)\n", + "plt.ylabel(r'', fontsize=20)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "2c0fdabe-8e8d-425a-8fc6-4722a3f3b664", + "metadata": {}, + "source": [ + "#### Stellar mass - Star formation rate plane comparison with Eazy derived physical properties" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7326c019-4fb8-48aa-a92b-e2468513783e", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "mah_params, ms_params, q_params = get_fit_params(skysim3_df_lsst_y1,mah_keys=DIFFMAH_KEYS,\n", + " ms_keys=DIFFSTAR_U_MS_KEYS,\n", + " q_keys=DIFFSTAR_U_Q_KEYS)\n", + "t_table = np.linspace(T_TABLE_MIN, 10**LGT0, DEFAULT_N_STEPS)\n", + "sfh_from_mah_kern = get_sfh_from_mah_kern(n_steps=DEFAULT_N_STEPS,\n", + " tacc_integration_min=T_BIRTH_MIN,\n", + " tobs_loop='vmap', galpop_loop='vmap')\n", + "sfh = sfh_from_mah_kern(t_table, mah_params, ms_params, q_params, LGT0, FB)\n", + "t_obs = age_at_z(skysim3_df_lsst_y1['redshift'].values, *DEFAULT_COSMOLOGY)\n", + "sfr_lsst_y1 = multiInterp2(t_obs, t_table, sfh)\n", + "_c = [None, 0]\n", + "_calc_smh_table_vmap = jjit(vmap(cumulative_mstar_formed, in_axes=_c))\n", + "args_smh_table = (t_table,sfh)\n", + "smh_table = _calc_smh_table_vmap(*args_smh_table)\n", + "logsmh_table = jnp.log10(smh_table)\n", + "log_sm_lsst_y1 = multiInterp2(t_obs, t_table, logsmh_table)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c0e03ed0-7eee-4517-a350-9eb8f31f5f4f", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "w = np.where((cosmos2020_lsst_y1['ez_mass']>0)&(10**(cosmos2020_lsst_y1['ez_sfr'])>0))\n", + "#idxs_cosmos = np.random.randint(0,len(cosmos2020[w]),100000) # 100k random samples\n", + "prop_cosmos = np.empty((len(cosmos2020_lsst_y1['ez_mass'][w]), 2))\n", + "prop_cosmos[:, 0] = cosmos2020_lsst_y1['ez_mass'][w]\n", + "prop_cosmos[:, 1] = cosmos2020_lsst_y1['ez_sfr'][w]\n", + "\n", + "w = np.where((cosmodc2_df_lsst_y1['baseDC2/sm'].values>0)&(cosmodc2_df_lsst_y1['baseDC2/sfr'].values>0))\n", + "prop_cosmodc2 = np.empty((len(cosmodc2_df_lsst_y1['baseDC2/sm'].values[w]), 2))\n", + "prop_cosmodc2[:, 0] = np.log10(cosmodc2_df_lsst_y1['baseDC2/sm'].values[w])\n", + "prop_cosmodc2[:, 1] = np.log10(cosmodc2_df_lsst_y1['baseDC2/sfr'].values[w])\n", + "\n", + "w = np.where((log_sm_lsst_y1>0)&(sfr_lsst_y1>0))\n", + "prop_skysim = np.empty((len(log_sm_lsst_y1[w]),2))\n", + "prop_skysim[:, 0] = log_sm_lsst_y1[w]\n", + "prop_skysim[:, 1] = np.log10(sfr_lsst_y1[w])\n", + "\n", + "plt.clf()\n", + "plt.rcParams['figure.figsize']=(10,8)\n", + "plt.rcParams['xtick.labelsize']=15\n", + "plt.rcParams['ytick.labelsize']=15\n", + "\n", + "fig = corner.corner(prop_cosmos, hist_kwargs=dict(density=True, linewidth=1.5,), #range=extents,\n", + " color='black', labels=[r'$\\log(M/M_{\\odot})_{obs}$',r'$\\log(SFR[M_{\\odot}/yr])_{obs}$'],\n", + " label_kwargs={\"fontsize\":15}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "fig = corner.corner(prop_cosmodc2, fig=fig, hist_kwargs=dict(density=True, linewidth=1.5), #range=extents,\n", + " color='blue', labels=[r'$\\log(M/M_{\\odot})_{obs}$',r'$\\log(SFR[M_{\\odot}/yr])_{obs}$'],\n", + " label_kwargs={\"fontsize\":15}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "fig = corner.corner(prop_skysim, fig=fig, hist_kwargs=dict(density=True, linewidth=1.5), #range=extents,\n", + " color='red', labels=[r'$\\log(M/M_{\\odot})_{obs}$',r'$\\log(SFR[M_{\\odot}/yr])_{obs}$'],\n", + " label_kwargs={\"fontsize\":15}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "axes = np.array(fig.axes)\n", + "proxy = [plt.Rectangle((0,0),1,1,fc = 'black'), \n", + " plt.Rectangle((0,0),1,1,fc = 'blue'),\n", + " plt.Rectangle((0,0),1,1,fc = 'red')]\n", + "axes[1].legend(proxy, ['COSMOS2020 (Eazy)','cosmoDC2 v1.1.4', 'skysim v3.1.0'],loc='upper right',fontsize=10)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "100ba1aa-66fa-4b08-a939-acedc5fe4ed6", + "metadata": {}, + "source": [ + "#### Stellar mass - Star formation rate plane comparison with LePhare derived physical properties" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8bfd84d6-816b-481f-b2ea-5c30ea670785", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "w = np.where((cosmos2020_lsst_y1['lp_mass_med']>0)&(cosmos2020_lsst_y1['lp_SFR_med']>-20))\n", + "#idxs_cosmos = np.random.randint(0,len(cosmos2020[w]),100000) # 100k random samples\n", + "prop_cosmos_lephare = np.empty((len(cosmos2020_lsst_y1['lp_mass_med'][w]), 2))\n", + "prop_cosmos_lephare[:, 0] = cosmos2020_lsst_y1['lp_mass_med'][w]\n", + "prop_cosmos_lephare[:, 1] = cosmos2020_lsst_y1['lp_SFR_med'][w]\n", + "\n", + "plt.clf()\n", + "plt.rcParams['figure.figsize']=(10,8)\n", + "plt.rcParams['xtick.labelsize']=15\n", + "plt.rcParams['ytick.labelsize']=15\n", + "\n", + "fig = corner.corner(prop_cosmos_lephare, hist_kwargs=dict(density=True, linewidth=1.5,), #range=extents,\n", + " color='black', labels=[r'$\\log(M/M_{\\odot})_{obs}$',r'$\\log(SFR[M_{\\odot}/yr])_{obs}$'],\n", + " label_kwargs={\"fontsize\":15}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "fig = corner.corner(prop_cosmodc2, fig=fig, hist_kwargs=dict(density=True, linewidth=1.5), #range=extents,\n", + " color='blue', labels=[r'$\\log(M/M_{\\odot})_{obs}$',r'$\\log(SFR[M_{\\odot}/yr])_{obs}$'],\n", + " label_kwargs={\"fontsize\":15}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "fig = corner.corner(prop_skysim, fig=fig, hist_kwargs=dict(density=True, linewidth=1.5), #range=extents,\n", + " color='red', labels=[r'$\\log(M/M_{\\odot})_{obs}$',r'$\\log(SFR[M_{\\odot}/yr])_{obs}$'],\n", + " label_kwargs={\"fontsize\":15}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "axes = np.array(fig.axes)\n", + "proxy = [plt.Rectangle((0,0),1,1,fc = 'black'), \n", + " plt.Rectangle((0,0),1,1,fc = 'blue'),\n", + " plt.Rectangle((0,0),1,1,fc = 'red')]\n", + "axes[1].legend(proxy, ['COSMOS2020 (LePhare)','cosmoDC2 v1.1.4', 'skysim v3.1.0'],loc='upper right',fontsize=10)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "e841fdab-11e4-4722-8cf6-238c21978ba7", + "metadata": {}, + "source": [ + "#### HSC observed-frame magnitudes and colours comparison between COSMOS2020 and skysim v3.1.0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f646931e-0716-4b12-a21e-2a4fce38cc9b", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "wavebands = ['g', 'r', 'i', 'z', 'y']\n", + "hsc_obsframe_mags_skysim = np.array([skysim3_df_lsst_y1['HSC_obs_{}'.format(waveband)].values for waveband in wavebands])\n", + "hsc_obsframe_colours_skysim = np.array([skysim3_df_lsst_y1['HSC_obs_{}'.format(wavebands[i])].values-skysim3_df_lsst_y1['HSC_obs_{}'.format(wavebands[i+1])].values for i in range(len(wavebands)-1)])\n", + "hsc_obsframe_mags_cosmos = np.array([cosmos2020_lsst_y1['HSC_{}_MAG_AUTO'.format(waveband)] for waveband in wavebands])\n", + "hsc_obsframe_colours_cosmos = np.array([cosmos2020_lsst_y1['HSC_{}_MAG_AUTO'.format(wavebands[i])]-cosmos2020_lsst_y1['HSC_{}_MAG_AUTO'.format(wavebands[i+1])] for i in range(len(wavebands)-1)])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fbd16a39-d120-4d72-a711-4960b59f33ad", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "plt.clf()\n", + "plt.rcParams['figure.figsize']=(10,8)\n", + "plt.rcParams['xtick.labelsize']=15\n", + "plt.rcParams['ytick.labelsize']=15\n", + "\n", + "fig = corner.corner(np.swapaxes(hsc_obsframe_mags_cosmos,0,1), hist_kwargs=dict(density=True, linewidth=1.5), #range=extents,\n", + " color='blue', labels=[r'$g_{obs}$',r'$r_{obs}$',r'$i_{obs}$',r'$z_{obs}$',r'$y_{obs}$'],\n", + " label_kwargs={\"fontsize\":30}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "fig = corner.corner(np.swapaxes(hsc_obsframe_mags_skysim,0,1), fig=fig, hist_kwargs=dict(density=True, linewidth=1.5), #range=extents,\n", + " color='red', labels=[r'$g_{obs}$',r'$r_{obs}$',r'$i_{obs}$',r'$z_{obs}$',r'$y_{obs}$'],\n", + " label_kwargs={\"fontsize\":30}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "axes = np.array(fig.axes)\n", + "proxy = [plt.Rectangle((0,0),1,1,fc = 'blue'), \n", + " plt.Rectangle((0,0),1,1,fc = 'red')]\n", + "axes[1].legend(proxy, ['COSMOS2020', 'skysim v3.1.0'],loc='upper right',fontsize=10)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3b246bc2-1761-4339-bc99-a978e41936a2", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "plt.clf()\n", + "plt.rcParams['figure.figsize']=(10,8)\n", + "plt.rcParams['xtick.labelsize']=15\n", + "plt.rcParams['ytick.labelsize']=15\n", + "\n", + "extents=[(-1.5,2.5),(-1.5,2.5),(-1.5,2.5),(-1.5,2.5)]\n", + "\n", + "fig = corner.corner(np.swapaxes(hsc_obsframe_colours_cosmos,0,1), hist_kwargs=dict(density=True, linewidth=1.5), range=extents,\n", + " color='blue', labels=[r'$(g-r)_{obs}$',r'$(r-i)_{obs}$',r'$(i-z)_{obs}$',r'$(z-y)_{obs}$'],\n", + " label_kwargs={\"fontsize\":30}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "fig = corner.corner(np.swapaxes(hsc_obsframe_colours_skysim,0,1), fig=fig, hist_kwargs=dict(density=True, linewidth=1.5), range=extents,\n", + " color='red', labels=[r'$(g-r)_{obs}$',r'$(r-i)_{obs}$',r'$(i-z)_{obs}$',r'$(z-y)_{obs}$'],\n", + " label_kwargs={\"fontsize\":30}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "axes = np.array(fig.axes)\n", + "proxy = [plt.Rectangle((0,0),1,1,fc = 'blue'), \n", + " plt.Rectangle((0,0),1,1,fc = 'red')]\n", + "axes[1].legend(proxy, ['COSMOS2020', 'skysim v3.1.0'],loc='upper right',fontsize=10)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0187a89d-bd5e-43f0-a115-05320865eaa5", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "colours = [r'$(g-r)_{obs}$',r'$(r-i)_{obs}$',r'$(i-z)_{obs}$',r'$(z-y)_{obs}$']\n", + "for i in range(len(colours)):\n", + " plt.clf()\n", + " plt.rcParams['figure.figsize']=(10,8)\n", + " plt.rcParams['xtick.labelsize']=25\n", + " plt.rcParams['ytick.labelsize']=25\n", + " plt.hist(hsc_obsframe_colours_cosmos[i,:], bins='auto', lw=2, ls='solid', color='blue', histtype='step',density=True, label='COSMOS2020')\n", + " plt.hist(hsc_obsframe_colours_skysim[i,:], bins='auto', lw=2, ls='dashed', color='red', histtype='step',density=True, label='skysim v3.1.0')\n", + " plt.legend(loc='upper right', fontsize=20)\n", + " plt.xlim(-1.5,4.5)\n", + " plt.xlabel(colours[i], fontsize=30)\n", + " plt.ylabel(r'', fontsize=20)\n", + " plt.tight_layout()\n", + " plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "41c9c04f-c3de-4828-9a8c-9013ac25bdc7", + "metadata": {}, + "source": [ + "#### Observed color-redshift relation comparison between COSMOS2020 and skysim v3.1.0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a35ee912-151c-4f66-8a64-394d12d31606", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "colours = [r'$(g-r)_{obs}$',r'$(r-i)_{obs}$',r'$(i-z)_{obs}$',r'$(z-y)_{obs}$']\n", + "for i in range(len(colours)):\n", + " plt.clf()\n", + " plt.rcParams['figure.figsize']=(10,8)\n", + " plt.rcParams['xtick.labelsize']=25\n", + " plt.rcParams['ytick.labelsize']=25\n", + " plt.scatter(cosmos2020_lsst_y1['ez_z_phot'], hsc_obsframe_colours_cosmos[i,:], s=1, color='blue', label='COSMOS2020 (Eazy)')\n", + " plt.scatter(skysim3_df_lsst_y1['redshift'], hsc_obsframe_colours_skysim[i,:], s=1, color='red', label='skysim v3.1.0',alpha=0.5)\n", + " plt.legend(loc='upper right', fontsize=20)\n", + " plt.xlim(0,3.5)\n", + " plt.ylim(-1.5,4.5)\n", + " plt.ylabel(colours[i], fontsize=30)\n", + " plt.xlabel(r'redshift', fontsize=20)\n", + " plt.tight_layout()\n", + " plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "45ceb235-cf56-497b-9d4b-4f722e6e116b", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "colours = [r'$(g-r)_{obs}$',r'$(r-i)_{obs}$',r'$(i-z)_{obs}$',r'$(z-y)_{obs}$']\n", + "for i in range(len(colours)):\n", + " plt.clf()\n", + " plt.rcParams['figure.figsize']=(10,8)\n", + " plt.rcParams['xtick.labelsize']=25\n", + " plt.rcParams['ytick.labelsize']=25\n", + " plt.scatter(cosmos2020_lsst_y1['lp_zPDF'], hsc_obsframe_colours_cosmos[i,:], s=1, color='blue', label='COSMOS2020 (LePhare)')\n", + " plt.scatter(skysim3_df_lsst_y1['redshift'], hsc_obsframe_colours_skysim[i,:], s=1, color='red', label='skysim v3.1.0')\n", + " plt.legend(loc='upper right', fontsize=20)\n", + " plt.xlim(0,3.5)\n", + " plt.ylim(-1.5,4.5)\n", + " plt.ylabel(colours[i], fontsize=30)\n", + " plt.xlabel(r'redshift', fontsize=20)\n", + " plt.tight_layout()\n", + " plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "b857b21e-ea36-432c-9428-98846985f3b7", + "metadata": {}, + "source": [ + "## Load an hdf5 table of rest-frame spectra generated with rail_fsps and rail_dsps using skysim v3.1.0 input properties" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e6161baa-848d-497b-836e-02ccf17418dd", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "output_directory = os.getcwd()\n", + "\n", + "if not os.path.isfile(os.path.join(output_directory,'skysim_v3.1.0_100k.csv')):\n", + " os.system('curl -O https://portal.nersc.gov/cfs/lsst/PZ/skysim_v3.1.0_100k.csv '\n", + " '--output-dir {}'.format(output_directory))\n", + "\n", + "if not os.path.isfile(os.path.join(output_directory,'skysim_v3.1.0_sed_library.h5')):\n", + " os.system('curl -O https://portal.nersc.gov/cfs/lsst/PZ/skysim_v3.1.0_sed_library.h5 '\n", + " '--output-dir {}'.format(output_directory))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "74313e12-110b-4460-a5fc-f38b118875d0", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "skysim_for_seds_df = pd.read_csv(os.path.join(output_directory, 'skysim_v3.1.0_100k.csv'))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d452814c-10b8-4039-bce2-6a8b2b220501", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "with h5py.File(os.path.join(output_directory, 'skysim_v3.1.0_sed_library.h5'), 'r') as h5table:\n", + " idx = h5table['index'][()]\n", + " fsps_ab_magnitudes = h5table['fsps_ab_magnitudes'][()]\n", + " fsps_rest_wavelength_A = h5table['fsps_rest_wavelength_A'][()]\n", + " fsps_rest_Lsolar_Hz = h5table['fsps_rest_Lsolar_Hz'][()]\n", + " dsps_ab_magnitudes = h5table['dsps_ab_magnitudes'][()]\n", + " dsps_rest_wavelength_A = h5table['dsps_rest_wavelength_A'][()]\n", + " dsps_rest_Lsolar_Hz = h5table['dsps_rest_Lsolar_Hz'][()]\n", + " filters_list = [name.decode('utf8') for name in h5table['filters_list'][()]]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3a81b022-a71c-41f5-89a7-a02106bee42c", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "print(filters_list)" + ] + }, + { + "cell_type": "markdown", + "id": "fa534109-e617-46eb-bfd9-9f14be38bfa7", + "metadata": {}, + "source": [ + "#### Comparison of rest-frame spectra generated with FSPS and DSPS on the same input skysim v3.1.0 galaxy properties" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d4c23fe6-4259-4624-85cd-23dc29dfd61f", + "metadata": {}, + "outputs": [], + "source": [ + "plt.clf()\n", + "plt.rcParams['figure.figsize']=(10,8)\n", + "plt.rcParams['xtick.labelsize']=25\n", + "plt.rcParams['ytick.labelsize']=25\n", + "\n", + "mask_fsps = np.where((fsps_rest_wavelength_A>=3500) & (fsps_rest_wavelength_A<=10000))\n", + "mask_dsps = np.where((dsps_rest_wavelength_A>=3500) & (dsps_rest_wavelength_A<=10000))\n", + "\n", + "plt.plot(fsps_rest_wavelength_A[mask_fsps], np.nanmedian(fsps_rest_Lsolar_Hz,axis=0)[mask_fsps], \n", + " lw=3, color='red', label='FSPS')\n", + "plt.plot(dsps_rest_wavelength_A[mask_dsps], np.nanmedian(dsps_rest_Lsolar_Hz,axis=0)[mask_dsps], \n", + " lw=1, color='green', label='DSPS')\n", + "\n", + "plt.legend(loc='upper right', fontsize=20)\n", + "plt.xlabel(r'Wavelength [$\\AA$]', fontsize=20)\n", + "plt.ylabel(r'luminosity density $L_{\\nu} [\\mathrm{L_{\\odot} Hz^{-1}}]$', fontsize=20)\n", + "plt.tight_layout()\n", + "plt.xlim(3000,10000)\n", + "#plt.xscale('log')\n", + "plt.yscale('log')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "68a9cc55-e9f6-4cb4-ace5-953f1a1b47d1", + "metadata": {}, + "outputs": [], + "source": [ + "frac_diff_dspsvsfsps_wrt_fsps = abs((dsps_rest_Lsolar_Hz-fsps_rest_Lsolar_Hz))/fsps_rest_Lsolar_Hz\n", + "\n", + "fifty_percent = np.percentile(frac_diff_dspsvsfsps_wrt_fsps, 50, axis=0,method=\"closest_observation\") \n", + "fifty_percent_idx = np.where(frac_diff_dspsvsfsps_wrt_fsps == fifty_percent)\n", + "\n", + "sixtyeight_percent = np.percentile(frac_diff_dspsvsfsps_wrt_fsps,68,axis=0,method=\"closest_observation\") \n", + "sixtyeight_percent_idx = np.where(frac_diff_dspsvsfsps_wrt_fsps == sixtyeight_percent)\n", + "\n", + "ninetyfive_percent = np.percentile(frac_diff_dspsvsfsps_wrt_fsps,95,axis=0,method=\"closest_observation\") \n", + "ninetyfive_percent_idx = np.where(frac_diff_dspsvsfsps_wrt_fsps == ninetyfive_percent)\n", + "\n", + "ninetynine_percent = np.percentile(frac_diff_dspsvsfsps_wrt_fsps,99,axis=0,method=\"closest_observation\") \n", + "ninetynine_percent_idx = np.where(frac_diff_dspsvsfsps_wrt_fsps == ninetynine_percent)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ef53e377-53a1-421b-a36e-0c0083b29090", + "metadata": {}, + "outputs": [], + "source": [ + "plt.clf()\n", + "plt.rcParams['figure.figsize']=(10,8)\n", + "plt.rcParams['xtick.labelsize']=25\n", + "plt.rcParams['ytick.labelsize']=25\n", + "\n", + "mask_dsps = np.where((dsps_rest_wavelength_A>=3500) & (dsps_rest_wavelength_A<=10000))\n", + "\n", + "plt.fill_between(dsps_rest_wavelength_A[mask_dsps], np.full(len(dsps_rest_wavelength_A[mask_dsps]), 0), \n", + " fifty_percent[mask_dsps], \n", + " facecolor=\"darkred\",label='50%')\n", + "plt.fill_between(dsps_rest_wavelength_A[mask_dsps], np.full(len(dsps_rest_wavelength_A[mask_dsps]), 0), \n", + " sixtyeight_percent[mask_dsps], \n", + " facecolor=\"red\", alpha=0.7,label='68%')\n", + "plt.fill_between(dsps_rest_wavelength_A[mask_dsps], np.full(len(dsps_rest_wavelength_A[mask_dsps]), 0), \n", + " ninetyfive_percent[mask_dsps], \n", + " facecolor=\"salmon\", alpha=0.5,label='95%')\n", + "plt.fill_between(dsps_rest_wavelength_A[mask_dsps], np.full(len(dsps_rest_wavelength_A[mask_dsps]), 0), \n", + " ninetynine_percent[mask_dsps], \n", + " facecolor=\"mistyrose\", alpha=0.5,label='99%')\n", + "\n", + "emission_lines_names = [r'$\\mathrm{OII_{doublet}}$',r'$\\mathrm{H_{\\beta}}$', r'$\\mathrm{OII_1}$+$\\mathrm{OIII_2}$', \n", + " r'$\\mathrm{NII_1}$+$\\mathrm{H_{\\alpha}}$+$\\mathrm{NII_2}$',\n", + " r'$\\mathrm{SII_1}$+$\\mathrm{SII_2}$',r'$\\mathrm{SIII_{12}}$',r'$\\mathrm{SIII_{22}}$']\n", + "emission_lines_wavelengths = [3726,4861, 4983, 6562, 6725, 9068, 9531]\n", + "colours_emlines = ['purple','blue', 'green', 'orange', 'red', 'magenta', 'pink']\n", + "lower_lim_emlines_wave = [3680, 4800, 4950, 6500, 6680, 9018, 9480]\n", + "upper_lim_emlines_wave = [3780, 4900, 5050, 6600, 6780, 9108, 9580]\n", + "colours_emlines = ['purple','blue', 'green', 'orange', 'red', 'magenta','pink']\n", + "y_coord_emlines = [0.07,0.09,0.07,0.09,0.07,0.09,0.07]\n", + "for k in range(len(emission_lines_wavelengths)):\n", + " plt.vlines(lower_lim_emlines_wave[k], ymin=0, ymax=0.3, color='black', lw=2)\n", + " plt.vlines(upper_lim_emlines_wave[k], ymin=0, ymax=0.3, color='black', lw=2)\n", + " plt.text(lower_lim_emlines_wave[k]+50, y_coord_emlines[k], emission_lines_names[k], fontsize=20)\n", + " plt.fill_betweenx(np.arange(0,0.31,0.01), lower_lim_emlines_wave[k], upper_lim_emlines_wave[k],\n", + " alpha=0.3, color=colours_emlines[k])\n", + "\n", + " \n", + "plt.legend(loc='center right', fontsize=20)\n", + "plt.xlabel(r'Wavelength [$\\AA$]', fontsize=20)\n", + "plt.ylabel(r'|DSPS-FSPS|/FSPS', fontsize=20)\n", + "plt.tight_layout()\n", + "plt.xlim(3000,10000)\n", + "plt.yscale('log')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "e4d7bfad-3ca8-430a-8baa-76527463cc58", + "metadata": {}, + "source": [ + "#### Comparison of optical and infrared colours between COSMOS2020 and those obtained via DSPS. We use the same CFHTLS, HSC and VISTA magnitudes as those contained in the COSMOS2020 catalog. We also cut at i=24.1 in the i-band." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f97d57a2-f148-4f7f-b69e-1e8019f835fb", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "skysim_mags_names = ['cfhtls_u','hsc_g', 'hsc_r', 'hsc_i', 'hsc_z', 'hsc_y','vista_vircam_J', 'vista_vircam_H',\n", + " 'vista_vircam_Ks']\n", + "cosmos_mags_names = ['CFHT_u_MAG_AUTO','HSC_g_MAG_AUTO','HSC_r_MAG_AUTO','HSC_i_MAG_AUTO','HSC_z_MAG_AUTO',\n", + " 'HSC_y_MAG_AUTO','UVISTA_J_MAG_AUTO','UVISTA_H_MAG_AUTO','UVISTA_Ks_MAG_AUTO']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5e883f61-73fb-466c-96b3-14518432fc88", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "idx_hsc_i = filters_list.index('hsc_i')\n", + "\n", + "dsps_y1_mask = np.where(dsps_ab_magnitudes[:, idx_hsc_i] <= 24.1)\n", + "\n", + "dsps_y1_ab_magnitudes = dsps_ab_magnitudes[dsps_y1_mask]\n", + "\n", + "idx_rand_cosmos2020_y1 = np.random.randint(low=0, high=len(cosmos2020_lsst_y1_all), size=10000)\n", + "cosmos2020_lsst_y1_red = cosmos2020_lsst_y1_all[idx_rand_cosmos2020_y1]\n", + "\n", + "dsps_colors = np.empty((len(dsps_y1_ab_magnitudes), 8))\n", + "cosmos2020_colors = np.empty((len(cosmos2020_lsst_y1_red), 8))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4a1bfa9a-0937-42e4-983d-912943d6dacd", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "for i in range(len(skysim_mags_names)-1):\n", + " idx_mag1 = filters_list.index(skysim_mags_names[i])\n", + " idx_mag2 = filters_list.index(skysim_mags_names[i+1])\n", + " dsps_colors[:, i] = dsps_y1_ab_magnitudes[:, idx_mag1] - dsps_y1_ab_magnitudes[:, idx_mag2]\n", + " cosmos2020_colors[:,i] = cosmos2020_lsst_y1_red[cosmos_mags_names[i]] - cosmos2020_lsst_y1_red[cosmos_mags_names[i+1]]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f7d648d0-e9e5-41ff-8487-4b391b075ffa", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "plt.clf()\n", + "\n", + "plt.rcParams['xtick.labelsize']=15\n", + "plt.rcParams['ytick.labelsize']=15\n", + "\n", + "extents = [(-1.5,4.), (-1.5,4.), (-1.5,4.), (-1.5,4.),(-1.5,4.)]\n", + "\n", + "fig = corner.corner(cosmos2020_colors[:,:5], hist_kwargs=dict(density=True, linewidth=1.5), range=extents,\n", + " color='blue', labels=colours[:5], label_kwargs={\"fontsize\":30}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "fig = corner.corner(dsps_colors[:,:5], fig=fig, hist_kwargs=dict(density=True, linewidth=1.5), range=extents,\n", + " color='red', labels=colours[:5], label_kwargs={\"fontsize\":30}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + " \n", + "axes = np.array(fig.axes)\n", + "proxy = [plt.Rectangle((0,0),1,1,fc = 'blue'), \n", + " plt.Rectangle((0,0),1,1,fc = 'red')]\n", + "\n", + "axes[3].legend(proxy, ['COSMOS2020','DSPS from skysim v3.1.0'],loc='upper right',fontsize=20)\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "38fad410-59fc-4f83-af41-a86b0871c249", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "plt.clf()\n", + "\n", + "plt.rcParams['xtick.labelsize']=15\n", + "plt.rcParams['ytick.labelsize']=15\n", + "\n", + "extents = [(-1.5,4.), (-1.5,4.), (-1.5,4.)]\n", + "\n", + "fig = corner.corner(cosmos2020_colors[:,5:], hist_kwargs=dict(density=True, linewidth=1.5), range=extents,\n", + " color='blue', labels=colours[5:], label_kwargs={\"fontsize\":30}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + "fig = corner.corner(dsps_colors[:,5:], fig=fig, hist_kwargs=dict(density=True, linewidth=1.5), range=extents,\n", + " color='red', labels=colours[5:], label_kwargs={\"fontsize\":30}, use_math_text=True,\n", + " plot_datapoints=False, levels=(1-np.exp(-0.5),1-np.exp(-2),1-np.exp(-4.5)), show_titles=False,\n", + " title_fmt='.4f')\n", + " \n", + "axes = np.array(fig.axes)\n", + "proxy = [plt.Rectangle((0,0),1,1,fc = 'blue'), \n", + " plt.Rectangle((0,0),1,1,fc = 'red')]\n", + "\n", + "axes[2].legend(proxy, ['COSMOS2020','DSPS from skysim v3.1.0'],loc='upper right',fontsize=20)\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "new_rail_env", + "language": "python", + "name": "new_rail_env" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pyproject.toml b/pyproject.toml index f937fda..0234ed5 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -35,6 +35,10 @@ dev = [ "pz-rail-pipelines[full]", "jupyter", "seaborn", + "corner", + "pandas", + "h5py", + "matplotlib", "coverage", "pytest", "pytest-cov", # Used to report total code coverage