Skip to content

Commit

Permalink
Updated tutorial on fitting a grid of model spectra
Browse files Browse the repository at this point in the history
  • Loading branch information
@tomasstolker
tomasstolker committed Jan 30, 2024
1 parent 4804241 commit eb0f34e
Show file tree
Hide file tree
Showing 3 changed files with 18 additions and 16 deletions.
29 changes: 15 additions & 14 deletions docs/tutorials/fitting_model_spectra.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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`."
]
},
{
Expand Down Expand Up @@ -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"
]
Expand All @@ -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."
]
},
{
Expand Down Expand Up @@ -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)."
]
},
{
Expand Down Expand Up @@ -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",
Expand All @@ -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."
]
},
{
Expand Down Expand Up @@ -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."
]
},
{
Expand Down Expand Up @@ -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."
]
},
{
Expand Down Expand Up @@ -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."
]
},
{
Expand All @@ -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)."
]
},
{
Expand Down Expand Up @@ -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)."
]
},
{
Expand Down Expand Up @@ -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)."
]
Expand Down
3 changes: 2 additions & 1 deletion species/fit/fit_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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'``):
Expand Down
2 changes: 1 addition & 1 deletion species/util/fit_util.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit eb0f34e

Please sign in to comment.