diff --git a/.github/workflows/smoke-test-head.yml b/.github/workflows/smoke-test-head.yml index 6a4db18..a109b90 100644 --- a/.github/workflows/smoke-test-head.yml +++ b/.github/workflows/smoke-test-head.yml @@ -35,6 +35,10 @@ jobs: run: | rail clone-source --package-file rail_packages.yml rail install --from-source --package-file rail_packages.yml + - name: Set up fsps + run: | + git clone "https://github.com/cconroy20/fsps.git" "/opt/hostedtoolcache/Python/fsps" + echo "SPS_HOME=/opt/hostedtoolcache/Python/fsps" >> $GITHUB_ENV - name: Run notebooks run: | rail render-nb --skip examples/creation_examples/plotting_interface_skysim_cosmoDC2_COSMOS2020_demo.ipynb examples/${{ matrix.notebook-dir }}_examples/*.ipynb diff --git a/.github/workflows/smoke-test.yml b/.github/workflows/smoke-test.yml index 9ab6071..a1750cc 100644 --- a/.github/workflows/smoke-test.yml +++ b/.github/workflows/smoke-test.yml @@ -35,6 +35,10 @@ jobs: - name: Run unit tests with pytest run: | python -m pytest tests + - name: Set up fsps + run: | + git clone "https://github.com/cconroy20/fsps.git" "/opt/hostedtoolcache/Python/fsps" + echo "SPS_HOME=/opt/hostedtoolcache/Python/fsps" >> $GITHUB_ENV - name: Run notebooks run: | rail render-nb --skip examples/creation_examples/plotting_interface_skysim_cosmoDC2_COSMOS2020_demo.ipynb examples/${{ matrix.notebook-dir }}_examples/*.ipynb diff --git a/examples/creation_examples/fsps_sed_demo.ipynb b/examples/creation_examples/fsps_sed_demo.ipynb new file mode 100644 index 0000000..c562b84 --- /dev/null +++ b/examples/creation_examples/fsps_sed_demo.ipynb @@ -0,0 +1,540 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Using rail_fsps to generate galaxy rest-frame spectral energy distributions and compute apparent magnitudes" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "author: Luca Tortorelli, Josue de Santiago, Eric Charles\n", + " \n", + "last run successfully: Aug 2nd, 2023\n", + " \n", + "This notebook demonstrates how to use rail_fsps to generate galaxy rest-frame spectral energy distributions (SEDs) with FSPS and how to compute apparent magnitudes from them.\n", + "\n", + "In order to run this notebook you need to have FSPS and Python-FSPS installed. The easiest way to do this is the following (first line applies only in case you already installed it via pip):\n", + "\n", + " pip uninstall fsps\n", + " git clone --recursive https://github.com/dfm/python-fsps.git\n", + " cd python-fsps\n", + " python -m pip install .\n", + " export SPS_HOME=$(pwd)/src/fsps/libfsps" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "from rail.core.data import TableHandle\n", + "from rail.creation.engines.fsps_sed_modeler import *\n", + "import rail.fsps\n", + "import h5py\n", + "import numpy as np" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We'll start by setting up the Rail data store. RAIL uses ceci, which is designed for pipelines rather than interactive notebooks, the data store will work around that and enable us to use data interactively. See the rail/examples/goldenspike/goldenspike.ipynb example notebook for more details on the Data Store." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from rail.core.stage import RailStage\n", + "DS = RailStage.data_store\n", + "DS.__class__.allow_overwrite = True\n", + "RAIL_FSPS_DIR = os.path.abspath(os.path.join(os.path.dirname(rail.fsps.__file__), '..', '..'))\n", + "default_rail_fsps_files_folder = os.path.join(RAIL_FSPS_DIR, 'rail', 'examples_data', 'creation_data', 'data',\n", + " 'fsps_default_data')\n", + "input_file = os.path.join(default_rail_fsps_files_folder, 'input_galaxy_properties_fsps.hdf5')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We generate some mock input data for the sed modeler class that we store in an hdf5 file." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "n_galaxies = 10\n", + "\n", + "redshifts = np.linspace(0.1,1,num=n_galaxies)\n", + "zmet = np.full(n_galaxies, 1, dtype=int)\n", + "stellar_metallicity = np.full(n_galaxies, 0.0) # log10(Z/Zsun)\n", + "pmetals = np.full(n_galaxies, 2.0)\n", + "stellar_velocity_dispersion = np.full(n_galaxies, 100.)\n", + "gas_ionization = np.full(n_galaxies, -1)\n", + "gas_metallicity = np.full(n_galaxies, 0.0)\n", + "tau = np.full(n_galaxies, 1.0)\n", + "const = np.full(n_galaxies, 0.0)\n", + "sf_start = np.full(n_galaxies, 0.0)\n", + "sf_trunc = np.full(n_galaxies, 0.0)\n", + "stellar_age = np.full(n_galaxies, 2.)\n", + "fburst = np.full(n_galaxies, 0.0)\n", + "tburst = np.full(n_galaxies, 11.0)\n", + "sf_slope = np.full(n_galaxies, 0.0)\n", + "dust1_birth_cloud = np.full(n_galaxies, 0.1)\n", + "dust2_diffuse = np.full(n_galaxies, 0.1)\n", + "dust_index = np.full(n_galaxies, -0.7)\n", + "dust_calzetti_modifier = np.full(n_galaxies, -1.)\n", + "mwr = np.full(n_galaxies, 3.1)\n", + "uvb = np.full(n_galaxies, 1.0)\n", + "wgp1 = np.full(n_galaxies, 1)\n", + "dust_gamma = np.full(n_galaxies, 0.01)\n", + "dust_umin = np.full(n_galaxies, 1.0)\n", + "dust_qpah = np.full(n_galaxies, 3.5)\n", + "f_agn = np.full(n_galaxies, 0.01)\n", + "tau_agn = np.full(n_galaxies, 10.0)\n", + "\n", + "gal_t_table = np.linspace(0.05, 13.8, 100) # age of the universe in Gyr\n", + "gal_sfr_table = np.random.uniform(0, 10, gal_t_table.size) # SFR in Msun/yr\n", + "tabulated_sfh = np.full((n_galaxies, 2, len(gal_sfr_table)), [gal_t_table,gal_sfr_table])\n", + "\n", + "wave_lsf = np.linspace(3000, 10000, 2000)\n", + "sigma_lsf = np.full_like(wave_lsf, 0.5)\n", + "tabulated_lsf = np.full((n_galaxies, 2, len(wave_lsf)), [wave_lsf, sigma_lsf])\n", + "\n", + "with h5py.File(input_file, 'w') as h5table:\n", + " data = h5table.create_group('model')\n", + " data.create_dataset(name='redshifts', data=redshifts)\n", + " data.create_dataset(name='zmet', data=zmet)\n", + " data.create_dataset(name='stellar_metallicity', data=stellar_metallicity)\n", + " data.create_dataset(name='pmetals', data=pmetals)\n", + " data.create_dataset(name='stellar_velocity_dispersion', data=stellar_velocity_dispersion)\n", + " data.create_dataset(name='gas_ionization', data=gas_ionization)\n", + " data.create_dataset(name='gas_metallicity', data=gas_metallicity)\n", + " data.create_dataset(name='tau', data=tau)\n", + " data.create_dataset(name='const', data=const)\n", + " data.create_dataset(name='sf_start', data=sf_start)\n", + " data.create_dataset(name='sf_trunc', data=sf_trunc)\n", + " data.create_dataset(name='stellar_age', data=stellar_age)\n", + " data.create_dataset(name='fburst', data=fburst)\n", + " data.create_dataset(name='tburst', data=tburst)\n", + " data.create_dataset(name='sf_slope', data=sf_slope)\n", + " data.create_dataset(name='dust1_birth_cloud', data=dust1_birth_cloud)\n", + " data.create_dataset(name='dust2_diffuse', data=dust2_diffuse)\n", + " data.create_dataset(name='dust_index', data=dust_index)\n", + " data.create_dataset(name='dust_calzetti_modifier', data=dust_calzetti_modifier)\n", + " data.create_dataset(name='mwr', data=mwr)\n", + " data.create_dataset(name='uvb', data=uvb)\n", + " data.create_dataset(name='wgp1', data=wgp1)\n", + " data.create_dataset(name='dust_gamma', data=dust_gamma)\n", + " data.create_dataset(name='dust_umin', data=dust_umin)\n", + " data.create_dataset(name='dust_qpah', data=dust_qpah)\n", + " data.create_dataset(name='f_agn', data=f_agn)\n", + " data.create_dataset(name='tau_agn', data=tau_agn)\n", + " data.create_dataset(name='tabulated_lsf', data=tabulated_lsf)\n", + " data.create_dataset(name='tabulated_sfh', data=tabulated_sfh)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here we load the data to the data store to use it later, and we also set the redshifts in which the output will be generated." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "trainFile = os.path.join(input_file)\n", + "training_data = DS.read_file(\"training_data\", TableHandle, trainFile)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's create an FSPSSedModeler class object. The latter has a number of configuration parameters to set. All the parameters have default values. Therefore, if the user is not sure about a particular value for a certain parameter, the latter can be left at default. A short description of each parameter can be found in src/rail/creation/engines/fsps_sed_modeler.py. \n", + "\n", + "We run the FSPSSedModeler in sequential mode. Note that each galaxy spectrum takes a few seconds to generate, so it is advisable to proceed in this way only for a limited sample of objects or for testing the code." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fspssedmodeler = FSPSSedModeler.make_stage(chunk_size=10, hdf5_groupname='model', name='FSPSSedModeler',\n", + " compute_vega_mags=False,vactoair_flag=False,\n", + " zcontinuous=1, add_agb_dust_model=True,\n", + " add_dust_emission=True, add_igm_absorption=True,\n", + " add_neb_emission=True, add_neb_continuum=True,\n", + " add_stellar_remnants=True, compute_light_ages=False,\n", + " nebemlineinspec=True, smooth_velocity=True,\n", + " smooth_lsf=False, cloudy_dust=False,\n", + " agb_dust=1.0, tpagb_norm_type=2, dell=0.0,\n", + " delt=0.0, redgb=1.0, agb=1.0, fcstar=1.0, sbss=0.0,\n", + " fbhb=0.0, pagb=1.0, redshifts_key='redshifts',\n", + " zmet_key='zmet', stellar_metallicities_key='stellar_metallicity',\n", + " pmetals_key='pmetals', imf_type=1, imf_upper_limit=120.,\n", + " imf_lower_limit=0.08, imf1=1.3, imf2=2.3,imf3=2.3,vdmc=0.08,\n", + " mdave=0.5,evtype=-1,use_wr_spectra=1,logt_wmb_hot=0.0,masscut=150.0,\n", + " velocity_dispersions_key='stellar_velocity_dispersion',min_wavelength=3000,\n", + " max_wavelength=10000,gas_ionizations_key='gas_ionization',\n", + " gas_metallicities_key='gas_metallicity',igm_factor=1.0,sfh_type=3,\n", + " tau_key='tau',const_key='const',sf_start_key='sf_start',\n", + " sf_trunc_key='sf_trunc',stellar_ages_key='stellar_age',\n", + " fburst_key='fburst',tburst_key='tburst',sf_slope_key='sf_slope',\n", + " dust_type=2,dust_tesc=7.0,dust_birth_cloud_key='dust1_birth_cloud',\n", + " dust_diffuse_key='dust2_diffuse',dust_clumps=-99,frac_nodust= 0.0,\n", + " frac_obrun=0.0,dust_index_key='dust_index',\n", + " dust_powerlaw_modifier_key='dust_calzetti_modifier',mwr_key='mwr',\n", + " uvb_key='uvb',wgp1_key='wgp1',wgp2=1,wgp3=1,\n", + " dust_emission_gamma_key='dust_gamma',dust_emission_umin_key='dust_umin',\n", + " dust_emission_qpah_key='dust_qpah',fraction_agn_bol_lum_key='f_agn',\n", + " agn_torus_opt_depth_key='tau_agn',tabulated_sfh_key='tabulated_sfh',\n", + " tabulated_lsf_key='tabulated_lsf',physical_units=False,\n", + " restframe_wave_key='restframe_wavelengths',\n", + " restframe_sed_key='restframe_seds')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here we show the example where we provide tabulated star-formation histories to generate the final SED. In this case, FSPS outputs the emission per total stellar mass. We call the fit_model() function to generate the model SEDs." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fspssedmodel = fspssedmodeler.fit_model(training_data)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fspssedmodel.data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We plot the first spectrum to check that the SED generation worked correctly." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, + "outputs": [], + "source": [ + "with h5py.File('model_FSPSSedModeler.hdf5','r') as h5table:\n", + " for key in h5table.keys():\n", + " print(key)\n", + " redshifts = h5table['redshifts'][()]\n", + " restframe_seds = h5table['restframe_seds'][()]\n", + " restframe_wavelengths = h5table['restframe_wavelengths'][()]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "plt.clf()\n", + "plt.plot(restframe_wavelengths[0], restframe_seds[0], lw=2, color='black')\n", + "plt.xlim(3000, 10000)\n", + "plt.xlabel(r'wavelength [$\\AA$]')\n", + "plt.ylabel(r'luminosity density [$\\mathrm{Lsun \\ Hz^{-1}}$]')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Running it with ceci" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "import ceci\n", + "pipe = ceci.Pipeline.interactive()\n", + "stages = [fspssedmodeler]\n", + "for stage in stages:\n", + " pipe.add_stage(stage)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "help(pipe.initialize)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pipe.initialize(dict(training_data=trainFile), dict(output_dir='./temp_output_rail_fsps', log_dir='./logs_rail_fsps',\n", + " resume=False, nprocess=2), None)\n", + "pipe.save('./temp_output_rail_fsps/pipe_saved.yml')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import ceci\n", + "pr = ceci.Pipeline.read('./temp_output_rail_fsps/pipe_saved.yml')\n", + "pr.run()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import tables_io\n", + "rest_frame_sed_models = tables_io.read('temp_output_rail_fsps/model_FSPSSedModeler.hdf5')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "plt.clf()\n", + "plt.plot(rest_frame_sed_models['restframe_wavelengths'][0], rest_frame_sed_models['restframe_seds'][1],\n", + " lw=2, color='black')\n", + "plt.xlim(3000, 10000)\n", + "plt.xlabel(r'wavelength [$\\AA$]')\n", + "plt.ylabel(r'luminosity density [$\\mathrm{Lsun \\ Hz^{-1}}$]')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we use these rest-frame SEDs to generate LSST photometry at user-provided redshifts.\n", + "For this we need the FSPSPhotometryCreator class. The class parameters are:\n", + "- redshifts_key: redshift keyword for the Hdf5 input table.\n", + "- restframe_sed_key: rest-frame sed keyword for the Hdf5 input table.\n", + "- restframe_wave_key: wavelength keyword for the Hdf5 input table.\n", + "- apparent_mags_key: apparent magnitudes keyword for the Hdf5 output table.\n", + "- filter_folder: path to the folder where filter bands are stored.\n", + "- instrument_name: name of the instrument for which we want to compute the magnitudes. The syntax of filenames should be of type instrument_band_transmission.h5. The wavelengths should be in units of Angstrom.\n", + "- wavebands: comma-separated list of filter bands.\n", + "- filter_wave_key: wavelength keyword in the hdf5 table of filter bands.\n", + "- filter_transm_key: transmission keyword in the hdf5 table of filter bands.\n", + "- Om0, Ode0, w0, wa, h: cosmological parameters for a w0waCDM cosmology\n", + "- use_planck_cosmology: True to overwrite the cosmological parameters with Planck15 cosmology model from Astropy\n", + "- physical_units: same meaning as above." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, + "outputs": [], + "source": [ + "from rail.core.stage import RailStage\n", + "import os\n", + "DS = RailStage.data_store\n", + "DS.__class__.allow_overwrite = True\n", + "from rail.core.data import TableHandle\n", + "import rail.fsps" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, + "outputs": [], + "source": [ + "trainFile = os.path.join('model_FSPSSedModeler.hdf5')\n", + "fsps_model = DS.read_file(\"fsps_model\", TableHandle, trainFile)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fspssedmodeler.get_handle('model').path" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from rail.creation.engines.fsps_photometry_creator import *\n", + "\n", + "RAIL_FSPS_DIR = os.path.abspath(os.path.join(os.path.dirname(rail.fsps.__file__), '..', '..'))\n", + "default_rail_fsps_filter_folder = os.path.join(RAIL_FSPS_DIR, 'rail', 'examples_data', 'creation_data', 'data',\n", + " 'fsps_default_data', 'filters')\n", + "\n", + "fspsphotometrycreator = FSPSPhotometryCreator.make_stage(name=\"FSPSPhotometryCreator\", \n", + " model=trainFile,\n", + " redshifts_key='redshifts',\n", + " redshift_key='redshifts',\n", + " restframe_sed_key='restframe_seds',\n", + " restframe_wave_key='restframe_wavelengths',\n", + " apparent_mags_key='apparent_mags',\n", + " filter_folder=default_rail_fsps_filter_folder,\n", + " instrument_name='lsst', wavebands='u,g,r,i,z,y',\n", + " filter_wave_key='wave', filter_transm_key='transmission',\n", + " Om0=0.3, Ode0=0.7, w0=-1, wa=0.0, h=0.7,\n", + " use_planck_cosmology=True, physical_units=False)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We call the sample method to generate LSST photometry once the class has been initialised. The output is a table in Hdf5 format, with three columns: sequential ID, redshifts, apparent AB magnitudes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fspsphotometry = fspsphotometrycreator.sample(input_data=fsps_model)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fspsphotometry.data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The RAIL stages can be chained together conveniently using Ceci. The following is an example of how the FSPSSedGenerator and FSPSPhotometryCreator can be run as part of the pipeline Ceci stages." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pipe = ceci.Pipeline.interactive()\n", + "stages = [fspssedmodeler, fspsphotometrycreator]\n", + "for stage in stages:\n", + " pipe.add_stage(stage)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pipe.initialize(dict(training_data=input_file, fsps_model=trainFile), dict(output_dir='./temp_output_rail_fsps', log_dir='./logs', resume=False, nprocess=2),'./temp_output_rail_fsps/pipe_saved_config.yml')\n", + "pipe.save('./temp_output_rail_fsps/pipe_2_saved.yml')\n", + "pr = ceci.Pipeline.read('./temp_output_rail_fsps/pipe_2_saved.yml')\n", + "pr.run()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The creation of galaxy SEDs is a computationally intensive process. To speed things up, it is convenient to parallelize the process using MPI. To do that, one needs to set the parameters into the configuration file pipe_saved_config.yml. An example command to be run in the command line with n_cores is:\n", + "\n", + "mpiexec -n n_cores --mpi python3 -m rail SedGenerator --input=input_galaxy_properties_fsps.h5 --name=sed_generator_test --config=pipe_saved_config.yml --output=output_sed_generator_test.h5\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "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": 4 +}