Source code for pyretlife.retrieval.priors
"""
PRIORS.py
Here, the unity cube is converted to prior cube by using either an
uninformative, uniform prior, or a Gaussian prior.
The variable cube has length equal to the number of parameters to be
retrieved; the indexing follows the order of the params global list.
"""
# -----------------------------------------------------------------------------
# IMPORTS
# -----------------------------------------------------------------------------
import scipy.stats as stat
import numpy as np
from typing import Union
from pathlib import Path
from numpy import ndarray
from scipy.interpolate import PchipInterpolator
from astropy import constants as const
# -----------------------------------------------------------------------------
# DEFINITIONS
# -----------------------------------------------------------------------------
[docs]
def assign_priors(dictionary: dict) -> dict:
"""
The assign_priors function takes a dictionary of parameters and assigns the appropriate prior function to each
parameter. The prior functions are defined in the priors module. The assign_priors function is called by the
read_parameters function, which reads in the parameters from a .yaml file.
:param dictionary: Pass the dictionary containing all the parameters and their values
:return: A dictionary with the prior function assigned to each parameter
"""
for parameter in dictionary.keys():
prior_kind = dictionary[parameter]["prior"]["kind"]
if prior_kind == "uniform":
dictionary[parameter]["prior"]["function"] = uniform_prior
elif prior_kind == "log-uniform":
dictionary[parameter]["prior"]["function"] = log_uniform_prior
elif prior_kind == "upper-log-uniform":
dictionary[parameter]["prior"]["function"] = upper_log_uniform_prior
elif prior_kind == "gaussian":
dictionary[parameter]["prior"]["function"] = gaussian_prior
elif prior_kind == "log-gaussian":
dictionary[parameter]["prior"]["function"] = log_gaussian_prior
elif prior_kind == "fourth-uniform":
dictionary[parameter]["prior"]["function"] = fourth_power_uniform_prior
elif prior_kind == "third-uniform":
dictionary[parameter]["prior"]["function"] = third_power_uniform_prior
elif prior_kind == "second-uniform":
dictionary[parameter]["prior"]["function"] = second_power_uniform_prior
elif prior_kind == "2d-uniform":
dictionary[parameter]["prior"]["function"] = twod_uniform_prior
elif prior_kind == "custom":
raise ValueError(
f"Custom prior selected for {parameter} not currently implemented. Please choose another valid prior! "
f"Exiting the run..."
)
# dictionary[parameter]["prior"]["function"] = custom_prior
# dictionary[parameter]["prior"]["prior_specs"]["prior_data"] = read_custom_prior(
# dictionary[parameter]["prior"]["prior_specs"]["prior_path"]
# )
else:
raise ValueError(
f"{parameter} does not have a valid prior. Please! choose a valid prior! "
f"Exiting the run..."
)
return dictionary
#
# def read_custom_prior(path: Union[str, Path]) -> ndarray:
# """
# The read_custom_prior function reads in a custom prior from the path specified by the user.
# The function takes one argument, which is a string or Path object specifying where to find the file containing
# the custom prior. The function returns a ndarray containing the data from the custom prior distribution.
#
# :param path: Specify that the path parameter can be either a string or a path object
# :return: A dictionary with the prior function assigned to each parameter
# """
# return np.loadtxt(path)
[docs]
def uniform_prior(r: float, prior_specs: dict) -> float:
"""
Scales a random number generated in a uniform prior between 0 and 1 to the respective value corresponding to a
uniform prior ranged between x1 and x2. This function is called by the assign_priors function.
A random number generated from a uniform prior between [x1, x2].
:param r: A random float generated from the uniform prior between [0, 1].
:param prior_specs: A dictionary of prior_specs containing the keywords "lower" and "upper"
:return: A random number generated from a uniform prior between [x1, x2].
"""
x1 = prior_specs["lower"]
x2 = prior_specs["upper"]
return x1 + r * (x2 - x1)
[docs]
def gaussian_prior(r: float, prior_specs: dict) -> float:
"""
Scales a random number generated in a uniform prior between 0 and 1 to the respective value corresponding to a
gaussian prior centered at mu and of standard deviation sigma. This function is called by the assign_priors
function.
:param r: A random float generated from the uniform prior between [0, 1].
:param prior_specs: A dictionary of prior_specs containing the keywords "mean" and "sigma"
:return: A random float generated from a gaussian prior G(mu,sigma).
"""
mu = prior_specs["mean"]
sigma = prior_specs["sigma"]
return stat.norm.ppf(r) * sigma + mu
[docs]
def log_uniform_prior(r: float, prior_specs: dict) -> float:
"""
Returns a number like 10^x, where x is generated by scaling a random number generated in a uniform prior between
0 and 1 to the respective value corresponding to a uniform prior ranged between x1 and x2. This function is
called by the assign_priors function.
:param r: A random float generated from the uniform prior between [0, 1].
:param prior_specs: A dictionary of prior_specs containing the keywords "log_lower" and "log_upper"
:return: A random number whose exponent in base 10 is generated from a uniform prior between [x1, x2].
"""
prior_logspace = {
"lower": prior_specs["log_lower"],
"upper": prior_specs["log_upper"],
}
return np.power(10, uniform_prior(r, prior_logspace))
[docs]
def upper_log_uniform_prior(r: float, prior_specs: dict) -> float:
"""
Returns a number like (1-10^x), where x is generated by scaling a random number generated in a uniform prior
between 0 and 1 to the respective value corresponding to a uniform prior ranged between x1 and x2. This function
is called by the assign_priors function.
:param r: A random float generated from the uniform prior between [0, 1].
:param prior_specs: A dictionary of prior_specs containing the keywords "log_lower" and "log_upper"
:return: A random number like (1-10^x) whose exponent in base 10 is generated from a uniform prior between [x1, x2].
"""
prior_logspace = {
"lower": prior_specs["log_lower"],
"upper": prior_specs["log_upper"],
}
return 1 - np.power(10, uniform_prior(r, prior_logspace))
[docs]
def log_gaussian_prior(r: float, prior_specs: dict) -> float:
"""
Returns a number like 10^x, where x is generated by scaling a random number generated in a uniform prior between
0 and 1 to the respective value corresponding to a gaussian prior ranged between x1 and x2. This function is
called by the assign_priors function.
:param r: A random float generated from the uniform prior between [0, 1].
:param prior_specs: A dictionary of prior_specs containing the keywords "log_mean" and "log_sigma"
:return: A random number whose exponent in base 10 is generated from a gaussian prior G(mu,sigma).
"""
prior_logspace = {
"mean": prior_specs["log_mean"],
"sigma": prior_specs["log_sigma"],
}
return np.power(10, gaussian_prior(r, prior_logspace))
[docs]
def fourth_power_uniform_prior(r: float, prior_specs: dict) -> float:
"""
The fourth_power_uniform_prior function takes in a random number and prior specifications. It then returns the
fourth power of a uniform distribution with lower and upper bounds specified by the prior_specs.
:param r: A random float generated from the uniform prior between [0, 1].
:param prior_specs: A dictionary of prior_specs containing the keywords "fourth_lower" and "fourth_upper"
:return: A random number like x^4 where x is generated from a uniform prior between [x1, x2]
"""
prior_fourth = {
"lower": prior_specs["fourth_lower"],
"upper": prior_specs["fourth_upper"],
}
return np.sign(uniform_prior(r, prior_fourth))*np.power(uniform_prior(r, prior_fourth), 4)
[docs]
def third_power_uniform_prior(r: float, prior_specs: dict) -> float:
"""
The fourth_power_uniform_prior function takes in a random number and prior specifications. It then returns the
fourth power of a uniform distribution with lower and upper bounds specified by the prior_specs.
:param r: A random float generated from the uniform prior between [0, 1].
:param prior_specs: A dictionary of prior_specs containing the keywords "fourth_lower" and "fourth_upper"
:return: A random number like x^4 where x is generated from a uniform prior between [x1, x2]
"""
prior_fourth = {
"lower": prior_specs["third_lower"],
"upper": prior_specs["third_upper"],
}
return np.power(uniform_prior(r, prior_fourth), 3)
[docs]
def second_power_uniform_prior(r: float, prior_specs: dict) -> float:
"""
The fourth_power_uniform_prior function takes in a random number and prior specifications. It then returns the
fourth power of a uniform distribution with lower and upper bounds specified by the prior_specs.
:param r: A random float generated from the uniform prior between [0, 1].
:param prior_specs: A dictionary of prior_specs containing the keywords "fourth_lower" and "fourth_upper"
:return: A random number like x^4 where x is generated from a uniform prior between [x1, x2]
"""
prior_fourth = {
"lower": prior_specs["second_lower"],
"upper": prior_specs["second_upper"],
}
return np.sign(uniform_prior(r, prior_fourth))*np.power(uniform_prior(r, prior_fourth), 2)
[docs]
def twod_uniform_prior(r: float, R_pl: int) -> float:
"""
Scales a random number generated in a uniform prior between 0 and 1 to the respective value corresponding to a
uniform prior ranged between the allowed mass values allowed by the current radius assigned prior.
This function is only used for M_pl. Formula taken from Zeng et al. (2016) doi.org/10.3847/0004-637X/819/2/127.
:param r: A random float generated from the uniform prior between [0, 1].
:param R_pl: The radius of the Planet in Earth radii.
:return: A random number generated from a uniform prior between [x1, x2].
"""
# defining the lower and upper bounds for the mass-radius relation
# values from Zeng et al. (2016), table 2
masses = np.array([0.125,0.25,0.5,1.0,2.0,4.0,8.0,16.0,32.0])
R100Fe = np.array([0.445,0.55,0.676,0.823,0.99,1.176,1.38,1.59,1.82])
R100H2O = np.array([0.776,0.952,1.163,1.41,1.71,2.05,2.45,2.9,3.36])
lowerf = PchipInterpolator(R100Fe, masses)
upperf = PchipInterpolator(R100H2O, masses)
x1 = lowerf(R_pl/const.R_earth.cgs.value)
x2 = upperf(R_pl/const.R_earth.cgs.value)
return (x1 + r * (x2 - x1)) * const.M_earth.cgs.value
[docs]
def custom_prior(r: float, prior_specs: dict) -> float:
"""
The custom_prior function takes in a random number r and the prior_specs dictionary.
It then returns the quantile of the data corresponding to that random number.
This is useful for when you want to use your own data as a prior.
:param r: A random float generated from the uniform prior between [0, 1].
:param prior_specs: A dictionary of prior_specs containing the keywords "data"
:return: The number corresponding to the r-th quantile of the custom prior
"""
return np.quantile(prior_specs["data"], r, axis=0)
[docs]
def invalid_prior(par: str):
"""
The invalid_prior function is used to raise an error if the user provides
an invalid prior. It takes a single argument, par, which is the name of
the parameter for which an invalid prior was provided. The function then
raises a ValueError with a message explaining that the prior provided for
par is invalid.
:param par: A string containing the name of the parameter
:raise ValueError:
"""
# Note: If you are exiting the run with an error, you should not use
# `sys.exit(0)` because that will produce a return code of 0, which
# usually means "success" or "no errors". In any case, raising a
# `ValueError` is probably the best way to go here, because it will
# give the user a clear error message and traceback, and will also
# exit the run with a non-zero return code.
# TODO: Add functionality for invalid priors
raise ValueError("The prior provided for " + str(par) + "is invalid.")