From eb0f34e70ae75eecee7ed30def3debc1f90faa2a Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Tue, 30 Jan 2024 10:30:23 +0100 Subject: [PATCH] Updated tutorial on fitting a grid of model spectra --- docs/tutorials/fitting_model_spectra.ipynb | 29 +++++++++++----------- species/fit/fit_model.py | 3 ++- species/util/fit_util.py | 2 +- 3 files changed, 18 insertions(+), 16 deletions(-) diff --git a/docs/tutorials/fitting_model_spectra.ipynb b/docs/tutorials/fitting_model_spectra.ipynb index 8f1d9f53..7be9e498 100644 --- a/docs/tutorials/fitting_model_spectra.ipynb +++ b/docs/tutorials/fitting_model_spectra.ipynb @@ -11,7 +11,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "In this tutorial, we will fit spectra and photometric fluxes of [beta Pic b](http://exoplanet.eu/catalog/beta_pic_b/) with the synthetic spectra from a radiative-convective equilibrium atmosphere model. Here we will use [DRIFT-PHOENIX](https://leap2010blog.wordpress.com/2014/07/03/drift-phoenix-atmosphere-models-creating-new-worlds/) but there are also [several other models](https://species.readthedocs.io/en/latest/overview.html#supported-data) supported by `species`." + "In this tutorial, we will fit spectra and photometric fluxes of [beta Pic b](https://exoplanet.eu/catalog/beta_pic_b--521/) with the synthetic spectra from a radiative-convective equilibrium atmosphere model. Here we will use [DRIFT-PHOENIX](https://leap2010blog.wordpress.com/2014/07/03/drift-phoenix-atmosphere-models-creating-new-worlds/) but there are also [several other models](https://species.readthedocs.io/en/latest/overview.html#supported-data) supported by `species`." ] }, { @@ -58,6 +58,7 @@ "from species.fit.fit_model import FitModel\n", "from species.read.read_model import ReadModel\n", "from species.plot.plot_mcmc import plot_posterior\n", + "from species.plot.plot_spectrum import plot_spectrum\n", "from species.util.box_util import update_objectbox\n", "from species.util.fit_util import get_residuals, multi_photometry" ] @@ -66,7 +67,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We initiate *species* by running the [SpeciesInit](https://species.readthedocs.io/en/latest/species.core.html#species.core.init.SpeciesInit) class. By doing so, both the configuration file and the HDF5 database are created in the working folder." + "We initiate *species* by running the [SpeciesInit](https://species.readthedocs.io/en/latest/species.core.html#species.core.species_init.SpeciesInit) class. By doing so, both the configuration file and the HDF5 database are created in the working folder." ] }, { @@ -239,7 +240,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We will also add [photometric data](https://github.com/tomasstolker/species/blob/master/species/data/companion_data.json) of beta Pic b that have been adopted in `species` from various references. To get an overview of all photometric data that were stored in the database, we can use the [list_companions](https://species.readthedocs.io/en/latest/species.data.html#species.data.database.Database.list_companions) method of the [Database](https://species.readthedocs.io/en/latest/species.data.html#species.data.database.Database) object (i.e. `database.list_companions()` in this case)." + "We will also add [photometric data](https://github.com/tomasstolker/species/blob/master/species/data/companion_data/companion_data.json) of beta Pic b that have been adopted in `species` from various references. To get an overview of all photometric data that were stored in the database, we can use the [list_companions](https://species.readthedocs.io/en/latest/species.data.html#species.data.database.Database.list_companions) method of the [Database](https://species.readthedocs.io/en/latest/species.data.html#species.data.database.Database) object (i.e. `database.list_companions()` in this case)." ] }, { @@ -494,11 +495,11 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "To fit the data with the grid of model spectra (i.e. *grid retrieval*), we start by creating an instance of [FitModel](https://species.readthedocs.io/en/latest/species.analysis.html#species.analysis.fit_model.FitModel), which requires the database tag of the object (`beta Pic b`) and the atmosphere model (`drift-phoenix`).\n", + "To fit the data with the grid of model spectra (i.e. *grid retrieval*), we start by creating an instance of [FitModel](https://species.readthedocs.io/en/latest/species.fit.html#species.fit.fit_model.FitModel), which requires the database tag of the object (`beta Pic b`) and the atmosphere model (`drift-phoenix`).\n", "\n", "The dictionary of the `bounds` parameter contains the prior boundaries for the various parameters. By default, the parameters of the atmosphere model are set to the full range of the available spectra (e.g. $\\log(g)$ and $\\mathrm{[Fe/H]}$ in the example below).\n", "\n", - "In addition to the parameters of the atmosphere model, several optional parameters can be fitted, for example to account for calibration systematics, extinction, and excess emission from a disk. More details on the various parameters can be found in the [API documentation](https://species.readthedocs.io/en/latest/species.analysis.html#species.analysis.fit_model.FitModel) of `FitModel`. As example, we fit a scaling parameter for the GPI $H$ band spectrum.\n", + "In addition to the parameters of the atmosphere model, several optional parameters can be fitted, for example to account for calibration systematics, extinction, and excess emission from a disk. More details on the various parameters can be found in the [API documentation](https://species.readthedocs.io/en/latest/species.fit.html#species.fit.fit_model.FitModel) of `FitModel`. As example, we fit a scaling parameter for the GPI $H$ band spectrum.\n", "\n", "The arguments of `inc_phot` and `inc_spec` are either a boolean or a list of photometry and spectra that were previously stored in the database for the object beta Pic b. We will use all the spectra but only the photometric data in the $L$ and $M$ bands.\n", "\n", @@ -525,7 +526,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We now create the instance of [FitModel](https://species.readthedocs.io/en/latest/species.analysis.html#species.analysis.fit_model.FitModel). By doing so, the grid of model spectra will be interpolated to the wavelengths of the data." + "We now create the instance of [FitModel](https://species.readthedocs.io/en/latest/species.fit.html#species.fit.fit_model.FitModel). By doing so, the grid of model spectra will be interpolated to the wavelengths of the data." ] }, { @@ -587,21 +588,21 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We are now ready to sample the posterior distribution by either using [MultiNest](https://johannesbuchner.github.io/PyMultiNest/index.html) or [UltraNest](https://johannesbuchner.github.io/UltraNest/index.html) with the [run_multinest](https://species.readthedocs.io/en/latest/species.analysis.html#species.analysis.fit_model.FitModel.run_multinest)) and [run_ultranest](https://species.readthedocs.io/en/latest/species.analysis.html#species.analysis.fit_model.FitModel.run_ultranest)) methods of `FitModel`. Both are nested sampling algorithms which are powerful in sampling multi-modal distributions and will estimate the marginalized likelihood (i.e. *model evidence*), which enables pair-wise model comparison through the Bayes factor." + "We are now ready to sample the posterior distribution by either using [MultiNest](https://johannesbuchner.github.io/PyMultiNest/index.html) or [UltraNest](https://johannesbuchner.github.io/UltraNest/index.html) with the [run_multinest](https://species.readthedocs.io/en/latest/species.fit.html#species.fit.fit_model.FitModel.run_multinest) and [run_ultranest](https://species.readthedocs.io/en/latest/species.fit.html#species.fit.fit_model.FitModel.run_ultranest) methods of `FitModel`. Both are nested sampling algorithms which are powerful in sampling multi-modal distributions and will estimate the marginalized likelihood (i.e. *model evidence*), which enables pair-wise model comparison through the Bayes factor." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "In this example, we will use the [run_multinest](https://species.readthedocs.io/en/latest/species.analysis.html#species.analysis.fit_model.FitModel.run_multinest) method since `UltraNest` is somewhat more straightforward to install than `MultiNest`. We specify the database tag were the posterior samples will be stored, the number of live points, the output folder where `MultiNest` will write its data, and we use the mass of beta Pic b from [Gravity Collaboration et al. 2020](https://ui.adsabs.harvard.edu/abs/2020A%26A...633A.110G/abstract) as Gaussian prior on $\\log(g)$." + "In this example, we will use the [run_multinest](https://species.readthedocs.io/en/latest/species.analysis.html#species.analysis.fit_model.FitModel.run_multinest) method. We specify the database tag were the posterior samples will be stored, the number of live points, the output folder where `MultiNest` will write its data, and we use the mass of beta Pic b from [Gravity Collaboration et al. 2020](https://ui.adsabs.harvard.edu/abs/2020A%26A...633A.110G/abstract) as Gaussian prior on $\\log(g)$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "To speed up the computation, it is possible to run the nested sampling in parallel (e.g. with `mpirun`) to benefit from the multiprocessing support by `UltraNest` and `MultiNest`. In that case it is important that any functions of `species` that will write to the [Database](https://species.readthedocs.io/en/latest/species.data.html#species.data.database.Database) will be commented out since simultaneous writing to the HDF5 database by different processes is not possible. It is therefore recommended to first add all the required data to the database and then only run [SpeciesInit](https://species.readthedocs.io/en/latest/species.core.html#species.core.init.SpeciesInit), [FitModel](https://species.readthedocs.io/en/latest/species.analysis.html#species.analysis.fit_model.FitModel), and the sampler ([run_multinest](https://species.readthedocs.io/en/latest/species.analysis.html#species.analysis.fit_model.FitModel.run_multinest) or [run_ultranest](https://species.readthedocs.io/en/latest/species.analysis.html#species.analysis.fit_model.FitModel.run_ultranest)) in parallel with MPI." + "To speed up the computation, it is possible to run the nested sampling in parallel (e.g. with `mpirun`) to benefit from the multiprocessing support by `UltraNest` and `MultiNest`. In that case it is important that any functions of `species` that will write to the [Database](https://species.readthedocs.io/en/latest/species.data.html#species.data.database.Database) will be commented out since simultaneous writing to the HDF5 database by different processes is not possible. It is therefore recommended to first add all the required data to the database and then only run [SpeciesInit](https://species.readthedocs.io/en/latest/species.core.html#species.core.species_init.SpeciesInit), [FitModel](https://species.readthedocs.io/en/latest/species.fit.html#species.fit.fit_model.FitModel), and the sampler ([run_multinest](https://species.readthedocs.io/en/latest/species.fit.html#species.fit.fit_model.FitModel.run_multinest) or [run_ultranest](https://species.readthedocs.io/en/latest/species.fit.html#species.fit.fit_model.FitModel.run_ultranest)) in parallel with MPI." ] }, { @@ -664,7 +665,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The samples from the parameter estimation have been stored in the database. We can now run the [plot_posterior](https://species.readthedocs.io/en/latest/species.plot.html#species.plot.plot_mcmc.plot_posterior) function to plot the 1D and 2D projections of the posterior distributions by making use of [corner.py](https://corner.readthedocs.io). Here we specify the database tag with the results from [run_ultranest](https://species.readthedocs.io/en/latest/species.analysis.html#species.analysis.fit_model.FitModel.run_ultranest) and we also include the derived posterior for the bolometric luminosity and mass." + "The samples from the parameter estimation have been stored in the database. We can now run the [plot_posterior](https://species.readthedocs.io/en/latest/species.plot.html#species.plot.plot_mcmc.plot_posterior) function to plot the 1D and 2D projections of the posterior distributions by making use of [corner.py](https://corner.readthedocs.io). Here we specify the database tag with the results from [run_multinest](https://species.readthedocs.io/en/latest/species.fit.html#species.fit.fit_model.FitModel.run_multinest) and we also include the derived posterior for the bolometric luminosity and mass." ] }, { @@ -923,7 +924,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We then run the [ObjectBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.ObjectBox) through the [update_objectbox](https://species.readthedocs.io/en/latest/species.util.html#species.util.read_util.update_objectbox) function with the best-fit parameters. This function applies the flux scaling and could also inflate the errors in case the `model_param` dictionary contains error inflation parameters." + "We then run the [ObjectBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.ObjectBox) through the [update_objectbox](https://species.readthedocs.io/en/latest/species.util.html#species.util.box_util.update_objectbox) function with the best-fit parameters. This function applies the flux scaling and could also inflate the errors in case the `model_param` dictionary contains error inflation parameters." ] }, { @@ -947,7 +948,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Next, we calculate the residuals of the best-fit model with [get_residuals](https://species.readthedocs.io/en/latest/species.util.html#species.util.phot_util.get_residuals). This function will interpolate the grid at the specified parameters, calculate synthetic photometry for the filters that are provided as argument of `inc_phot`, and smooth and resample the spectra to the resolution and wavelengths of the spectra that are provided as argument of `inc_spec` (all spectra of the [ObjectBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.ObjectBox) in this case). The returned residuals will be stored in a [ResidualsBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.ResidualsBox)." + "Next, we calculate the residuals of the best-fit model with [get_residuals](https://species.readthedocs.io/en/latest/species.util.html#species.util.fit_util.get_residuals). This function will interpolate the grid at the specified parameters, calculate synthetic photometry for the filters that are provided as argument of `inc_phot`, and smooth and resample the spectra to the resolution and wavelengths of the spectra that are provided as argument of `inc_spec` (all spectra of the [ObjectBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.ObjectBox) in this case). The returned residuals will be stored in a [ResidualsBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.ResidualsBox)." ] }, { @@ -988,7 +989,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Finally, we compute synthetic photometry from the best-fit model spectrum with the [multi_photometry](https://species.readthedocs.io/en/latest/species.util.html#species.util.phot_util.multi_photometry) function. Since we will plot all available photometric data of the [ObjectBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.ObjectBox), we will calculate synthetic photometry for all the names of the `filters` attribute of the [ObjectBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.ObjectBox)." + "Finally, we compute synthetic photometry from the best-fit model spectrum with the [multi_photometry](https://species.readthedocs.io/en/latest/species.util.html#species.util.fit_util.multi_photometry) function. Since we will plot all available photometric data of the [ObjectBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.ObjectBox), we will calculate synthetic photometry for all the names of the `filters` attribute of the [ObjectBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.ObjectBox)." ] }, { @@ -1051,7 +1052,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Now that we have gather all the `Box` objects with the results, we can pass them as list to the `boxes` parameter of the [plot_spectrum](https://species.readthedocs.io/en/latest/species.plot.html#species.plot.plot_spectrum.plot_spectrum) function. The [ResidualsBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.ResidualsBox) is separately provided as argument of `residuals`. We also include a list of filter names for which the transmission profiles will be plotted. The arguments of `residuals` and `filters` can also be set to `None` to not include these data in the plot.\n", + "Now that we have gathered all the `Box` objects with the results, we can pass them as list to the `boxes` parameter of the [plot_spectrum](https://species.readthedocs.io/en/latest/species.plot.html#species.plot.plot_spectrum.plot_spectrum) function. The [ResidualsBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.ResidualsBox) is separately provided as argument of `residuals`. We also include a list of filter names for which the transmission profiles will be plotted. The arguments of `residuals` and `filters` can also be set to `None` to not include these data in the plot.\n", "\n", "The somewhat complex part of [plot_spectrum](https://species.readthedocs.io/en/latest/species.plot.html#species.plot.plot_spectrum.plot_spectrum) is the optional `plot_kwargs` parameter. The argument is a list with the same length as the list of `boxes`. Each item contains a dictionary with keyword arguments that can be used to fine-tune the styling of the plot. The item of the [SynphotBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.SynphotBox) can be set to `None` because it is automatically adjusted to the styles that are used for the [ObjectBox](https://species.readthedocs.io/en/latest/species.core.html#species.core.box.ObjectBox)." ] diff --git a/species/fit/fit_model.py b/species/fit/fit_model.py index b75601a4..7285ae40 100644 --- a/species/fit/fit_model.py +++ b/species/fit/fit_model.py @@ -173,7 +173,8 @@ def __init__( linearly sampled (``flux_scaling``) or logarithmically sampled (``log_flux_scaling``). Additionally, it is also possible to fit a flux offset (``flux_offset``), - which adds a constant flux to the model spectrum. + which adds a constant flux (in W m-2 um-1) to the + model spectrum. Blackbody parameters (with ``model='planck'``): diff --git a/species/util/fit_util.py b/species/util/fit_util.py index 1fcf3446..4660a147 100644 --- a/species/util/fit_util.py +++ b/species/util/fit_util.py @@ -211,7 +211,7 @@ def get_residuals( provided. radtrans : ReadRadtrans, None Instance of :class:`~species.read.read_radtrans.ReadRadtrans`. - Only required with ``spectrum='petitradtrans'`. Make sure that + Only required with ``spectrum='petitradtrans'``. Make sure that the ``wavel_range`` of the ``ReadRadtrans`` instance is sufficiently broad to cover all the photometric and spectroscopic data of ``inc_phot`` and ``inc_spec``. Not used