Retrievals: Emission Spectra Retrieval
This retrieval is based on Mollière+20 for HR8799e, and shows a more realistic example of how to set up a retrieval. Emission spectra retrievals, particularly when multiple datasets are included can take a lot of computational resources to run. Therefore, this notebook outlines the setup for such a retrieval, but we advise only running this retrieval on a server or cluster.
For this example, the module poor_mans_noneq_chemistry is used to solve disequilibrium chemistry, and the cloud_cond module is used for condensation. You can import other models from models.py, or write your own function, ensuring that it takes in the same set of inputs and outputs, which are documented in models, see the API documentation.
This model uses a simple adaptive mesh refinement (AMR) algorithm to improve the pressure resolution around the location of the cloud bases.
Getting started
You should already have an installation of pRT on your machine, please see the installation manual if you have not installed it yet. pyMultiNest is also required. Ultranest is required if you want to use that as your sampler rather than pyMultiNest. See the Ultranest documentation for why you might want to choose this method. Using nested sampling rather than MCMC is faster, handles multimodal cases better, and directly provides the Bayesian evidence, which allows for easier model comparison.
In this tutorial, we will outline the process of setting up a RetrievalConfig object, which is the class used to set up a pRT retrieval. The basic process will always be to set up the configuration, and then pass it to the Retrieval class to run the retrieval using pyMultiNest. Several standard plotting outputs will also be produced by the retrieval class. Most of the classes and functions used in this tutorial have more advanced features than what will be explained here, so it’s highly recommended to take a look at the code and API documentation. There should be enough flexibility built in to cover most typical retrieval studies, but if you have feature requests please get in touch, or open an issue on gitlab.
[3]:
# Let's start by importing everything we need.
import os
# To not have numpy start parallelizing on its own
os.environ["OMP_NUM_THREADS"] = "1"
#os.environ["pRT_input_data_path"] = "/path/to/petitRADTRANS/petitRADTRANS/input_data/"
import numpy as np
import matplotlib.pyplot as plt
from petitRADTRANS import Radtrans
from petitRADTRANS import nat_cst as nc
# Import the class used to set up the retrieval.
from petitRADTRANS.retrieval import Retrieval,RetrievalConfig
# Import Prior functions, if necessary.
from petitRADTRANS.retrieval.util import gaussian_prior
# Import atmospheric model function
from petitRADTRANS.retrieval.models import emission_model_diseq
[4]:
##########################
# Define the pRT run setup
##########################
# retrieve.py expects this object to be called RunDefinition.
RunDefinition = RetrievalConfig(retrieval_name = "HR8799e_example", # Give a useful name for your retrieval
run_mode = "retrieval", # 'retrieval' to run, or 'evaluate' to make plots
AMR = True, # Adaptive mesh refinement, slower if True
scattering = True) # Add scattering for emission spectra clouds
For this example we include the GRAVITY data as published in Mollière+20. To reproduce the published results, please also include the archival SPHERE and GPI data from Zurlo et al. (2015) and Greenbaum et al. (2016).
[5]:
##################
# Read in Data
##################
import petitRADTRANS # need to get the name for the example data
path_to_data = petitRADTRANS.__file__.split('__init__.py')[0] # Default location for the example
RunDefinition.add_data('GRAVITY',
f"{path_to_data}retrieval/examples/emission/observations/HR8799e_Spectra.fits",
data_resolution = 500,
model_resolution = 1000,
model_generating_function = emission_model_diseq)
Photometric data
If we want to add photometry, we can do that as well! The photometry file should have the format:
# Name, lower wavelength edge [um], upper wavelength edge [um], flux density [W/m2/micron], flux error [W/m2/micron]
You are required to provide a model function for calculating the spectrum, as with spectral data, but also a photometric transformation function, which is used to convert the model spectrum into synthetic photometry. This would typically make use of the transmission function for a particular filter. We recommend the use of the species package (https://species.readthedocs.io/), in particular the SyntheticPhotometry module to provide these functions. If no function is provided, the
RetrievalConfig will attempt to import species to use this module, using the name provided as the filter name.
If you are using transmission spectra, your photometric transformation function should model the difference between the clear and occulted stellar spectrum, returning the difference in (planet radius/stellar radius)^2.
[8]:
RunDefinition.add_photometry(path_to_data + 'retrieval/examples/emission/observations/hr8799e_photometry.txt',
emission_model_diseq,
model_resolution = 40)
==============
species v0.7.4
==============
Working folder: /Users/nasedkin/python-packages/petitRADTRANS/docs/content/notebooks
Configuration settings:
- Database: /Users/nasedkin/python-packages/petitRADTRANS/docs/content/notebooks/species_database.hdf5
- Data folder: /Users/nasedkin/python-packages/petitRADTRANS/docs/content/notebooks/data
- Interpolation method: linear
- Magnitude of Vega: 0.03
Adding filter: Keck/NIRC2.Ks... [DONE]
Adding Vega spectrum... [DONE]
Reference: Bohlin et al. 2014, PASP, 126
URL: https://ui.adsabs.harvard.edu/abs/2014PASP..126..711B/abstract
Adding filter: Paranal/NACO.Lp... [DONE]
Adding filter: Paranal/NACO.NB405... [DONE]
Adding filter: Paranal/SPHERE.IRDIS_B_J... [DONE]
Adding filter: Paranal/SPHERE.IRDIS_D_H23_2... [DONE]
Adding filter: Paranal/SPHERE.IRDIS_D_H23_3... [DONE]
Adding filter: Paranal/SPHERE.IRDIS_D_K12_1... [DONE]
Adding filter: Paranal/SPHERE.IRDIS_D_K12_2... [DONE]
Parameters and Priors
Here we add all of the parameters used in the retrieval. Each parameter must be given a name, which matches the name used in the model function. Parameters can be set to fixed or free. Fixed parameters must be given a value, while free parameters are given a function that transforms the unit hypercube generated by multinest into physical prior space. Various prior functions are stored in util, see the API documentation.
[9]:
#################################################
# Add parameters, and priors for free parameters.
#################################################
# This run uses the model of Molliere (2020) for HR8799e
# The lambda function provide uniform priors.
# Distance to the planet in cm
RunDefinition.add_parameter(name = 'D_pl', free = False, value = 41.2925*nc.pc)
# Log of the surface gravity in cgs units.
RunDefinition.add_parameter('log_g',True,
transform_prior_cube_coordinate = \
lambda x : 2.+3.5*x)
# Planet radius in cm
RunDefinition.add_parameter('R_pl', True, \
transform_prior_cube_coordinate = \
lambda x : ( 0.7+1.3*x)*nc.r_jup_mean)
# Temperature in Kelvin
RunDefinition.add_parameter('T_int', True, \
transform_prior_cube_coordinate = \
lambda x : 300.+2000.*x,
value=0.0)
# Spline temperature structure parameters. T1 < T2 < T3
# As these priors depend on each other, they are implemented in the model function
RunDefinition.add_parameter('T3', True, \
transform_prior_cube_coordinate = \
lambda x : x,
value=0.0)
RunDefinition.add_parameter('T2', True, \
transform_prior_cube_coordinate = \
lambda x : x,
value=0.0)
RunDefinition.add_parameter('T1', True, \
transform_prior_cube_coordinate = \
lambda x : x)
# Optical depth model
# power law index in tau = delta * press_cgs**alpha
RunDefinition.add_parameter('alpha', True, \
transform_prior_cube_coordinate = \
lambda x :1.0+x)
# proportionality factor in tau = delta * press_cgs**alpha
RunDefinition.add_parameter('log_delta', True, \
transform_prior_cube_coordinate = \
lambda x : x)
# Chemistry
# A 'free retrieval' would have each line species as a parameter
# Using a (dis)equilibrium model, we only supply bulk parameters.
# Carbon quench pressure
RunDefinition.add_parameter('log_pquench',True,\
transform_prior_cube_coordinate = \
lambda x : -6.0+9.0*x)
# Metallicity
RunDefinition.add_parameter('Fe/H',True,\
transform_prior_cube_coordinate = \
lambda x : -1.5+3.0*x)
# C/O ratio
RunDefinition.add_parameter('C/O',True,\
transform_prior_cube_coordinate = \
lambda x : 0.1+1.5*x)
# Clouds
# Based on an Ackermann-Marley (2001) cloud model
# Width of particle size distribution
RunDefinition.add_parameter('sigma_lnorm', True,
transform_prior_cube_coordinate = \
lambda x : 1.05 + 1.95*x)
# Vertical mixing parameters
RunDefinition.add_parameter('log_kzz',True,\
transform_prior_cube_coordinate = \
lambda x : 5.0+8.*x)
# Sedimentation parameter
RunDefinition.add_parameter('fsed',True,\
transform_prior_cube_coordinate = \
lambda x : 1.0 + 10.0*x)
[10]:
#######################################################
# Define species to be included as absorbers
#######################################################
RunDefinition.set_rayleigh_species(['H2', 'He'])
RunDefinition.set_continuum_opacities(['H2-H2', 'H2-He'])
RunDefinition.set_line_species(['CH4',
'H2O_HITEMP',
'CO2',
'HCN',
'CO_all_iso_HITEMP',
'FeH',
'H2S',
'NH3',
'PH3',
'Na_allard',
'K_allard'],
eq = True)
RunDefinition.add_cloud_species('Fe(c)_cd',eq = True,scaling_factor=(-2.5,4.5))
RunDefinition.add_cloud_species('MgSiO3(c)_cd',eq = True,scaling_factor=(-2.5,4.5))
[11]:
# Before we run the retrieval, let's set up plotting.
##################################################################
# Define what to put into corner plot if run_mode == 'evaluate'
##################################################################
RunDefinition.parameters['R_pl'].plot_in_corner = True
RunDefinition.parameters['R_pl'].corner_label = r'$R_{\rm P}$ ($\rm R_{Jup}$)'
RunDefinition.parameters['R_pl'].corner_transform = lambda x : x/nc.r_jup_mean
RunDefinition.parameters['log_g'].plot_in_corner = True
RunDefinition.parameters['log_g'].corner_ranges = [2., 5.]
RunDefinition.parameters['log_g'].corner_label = "log g"
RunDefinition.parameters['fsed'].plot_in_corner = True
RunDefinition.parameters['log_kzz'].plot_in_corner = True
RunDefinition.parameters['log_kzz'].corner_label = "log Kzz"
RunDefinition.parameters['C/O'].plot_in_corner = True
RunDefinition.parameters['Fe/H'].plot_in_corner = True
RunDefinition.parameters['log_pquench'].plot_in_corner = True
RunDefinition.parameters['log_pquench'].corner_label = "log pquench"
for spec in RunDefinition.cloud_species:
cname = spec.split('_')[0]
RunDefinition.parameters['eq_scaling_'+cname].plot_in_corner = True
RunDefinition.parameters['eq_scaling_'+cname].corner_label = cname
##################################################################
# Define axis properties of spectral plot if run_mode == 'evaluate'
##################################################################
RunDefinition.plot_kwargs["spec_xlabel"] = 'Wavelength [micron]'
RunDefinition.plot_kwargs["spec_ylabel"] = "Flux [W/m2/micron]"
RunDefinition.plot_kwargs["y_axis_scaling"] = 1.0
RunDefinition.plot_kwargs["xscale"] = 'linear'
RunDefinition.plot_kwargs["yscale"] = 'linear'
RunDefinition.plot_kwargs["resolution"] = 1000.# Maximum resolution, will bin any higher resolution data to this resolution
RunDefinition.plot_kwargs["nsample"] = 100 # If we want a plot with many sampled spectra
##################################################################
# Define from which observation object to take P-T
# in evaluation mode (if run_mode == 'evaluate'),
# add PT-envelope plotting options
##################################################################
RunDefinition.plot_kwargs["take_PTs_from"] = 'GRAVITY'
RunDefinition.plot_kwargs["temp_limits"] = [150, 3000]
RunDefinition.plot_kwargs["press_limits"] = [1e2, 1e-5]
Now let’s run the retrieval. Note that for this setup we’re using 4000 live points, and retrieving an emission spectrum on moderate resolution data. This will take some time: probably a day or so running on a cluster, so don’t try to run this on your laptop!
[ ]:
retrieval = Retrieval(RunDefinition,
output_dir = "",
use_MPI=True,
sample_spec = False)
retrieval.run(n_live_points = 4000,
sampling_efficiency=0.05,
const_efficiency_mode=True,
resume = True)
Once the retrieval is complete, the easiest way to generate standard output plots is to use the plot_all function.
[ ]:
retrieval.plot_all(contribution = True)
Contact
If you need any additional help, don’t hesitate to contact Evert Nasedkin.
[ ]: