Skip to content

Commit

Permalink
Added normal_prior parameter to FitModel and deprecated the prior par…
Browse files Browse the repository at this point in the history 
…ameter in run_multinest and run_ultranest
  • Loading branch information
@tomasstolker
tomasstolker committed Feb 4, 2024
1 parent b887861 commit 936dd32
Showing 1 changed file with 82 additions and 64 deletions.
146 changes: 82 additions & 64 deletions species/fit/fit_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,9 +10,6 @@
import numpy as np
import spectres

from PyAstronomy.pyasl import fastRotBroad
from scipy import interpolate, stats

try:
import ultranest
except:
Expand All @@ -31,6 +28,8 @@
"(Linux) or DYLD_LIBRARY_PATH (Mac)?"
)

from PyAstronomy.pyasl import fastRotBroad
from scipy import interpolate, stats
from typeguard import typechecked

from species.core import constants
Expand Down Expand Up @@ -92,6 +91,7 @@ def __init__(
fit_corr: Optional[List[str]] = None,
apply_weights: Union[bool, Dict[str, Union[float, np.ndarray]]] = False,
ext_filter: Optional[str] = None,
normal_prior: Optional[Dict[str, Tuple[float, float]]] = None,
) -> None:
"""
Parameters
Expand Down Expand Up @@ -403,6 +403,16 @@ def __init__(
svo/theory/fps/>`_ as argument then the extinction
``ism_ext`` is fitted in that filter instead of the
$V$ band.
normal_prior : dict(str, tuple(float, float)), None
Dictionary with normal priors for one or multiple
parameters. The prior can be set for any of the
atmosphere or calibration parameters, e.g.
``prior={'teff': (1200., 100.)}``, for a prior distribution
with a mean of 1200 K and a standard deviation of 100 K.
Additionally, a prior can be set for the mass, e.g.
``prior={'mass': (13., 3.)}`` for an expected mass of 13
Mjup with an uncertainty of 3 Mjup. The parameter is not
used if the argument is set to ``None``.
Returns
-------
Expand All @@ -418,6 +428,17 @@ def __init__(
"The 'bounds' dictionary should contain 'teff' and 'radius'."
)

if model == "bt-settl":
warnings.warn(
"It is recommended to use the CIFIST "
"grid of the BT-Settl, because it is "
"a newer version. In that case, set "
"model='bt-settl-cifist' when using "
"add_model of Database."
)

# Set attributes

self.object = ReadObject(object_name)
self.parallax = self.object.get_parallax()
self.binary = False
Expand All @@ -431,14 +452,12 @@ def __init__(
self.model = model
self.bounds = bounds

if self.model == "bt-settl":
warnings.warn(
"It is recommended to use the CIFIST "
"grid of the BT-Settl, because it is "
"a newer version. In that case, set "
"model='bt-settl-cifist' when using "
"add_model of Database."
)
if normal_prior is None:
self.normal_prior = {}
else:
self.normal_prior = normal_prior

# Set model parameters and boundaries

if self.model == "planck":
# Fitting blackbody radiation
Expand Down Expand Up @@ -761,17 +780,17 @@ def __init__(
self.modelspec = []

if self.model != "planck":
for key, value in self.spectrum.items():
print(f"\rInterpolating {key}...", end="", flush=True)
for spec_key, spec_value in self.spectrum.items():
print(f"\rInterpolating {spec_key}...", end="", flush=True)

wavel_range = (0.9 * value[0][0, 0], 1.1 * value[0][-1, 0])
wavel_range = (0.9 * spec_value[0][0, 0], 1.1 * spec_value[0][-1, 0])

readmodel = ReadModel(self.model, wavel_range=wavel_range)

readmodel.interpolate_grid(
wavel_resample=self.spectrum[key][0][:, 0],
wavel_resample=spec_value[0][:, 0],
smooth=True,
spec_res=self.spectrum[key][3],
spec_res=spec_value[3],
)

self.modelspec.append(readmodel)
Expand Down Expand Up @@ -819,17 +838,17 @@ def __init__(
self.diskphot.append(readmodel)
print(" [DONE]")

for key, value in self.spectrum.items():
print(f"\rInterpolating {key}...", end="", flush=True)
for spec_key, spec_value in self.spectrum.items():
print(f"\rInterpolating {spec_key}...", end="", flush=True)

wavel_range = (0.9 * value[0][0, 0], 1.1 * value[0][-1, 0])
wavel_range = (0.9 * spec_value[0][0, 0], 1.1 * spec_value[0][-1, 0])

readmodel = ReadModel("blackbody", wavel_range=wavel_range)

readmodel.interpolate_grid(
wavel_resample=self.spectrum[key][0][:, 0],
wavel_resample=spec_value[0][:, 0],
smooth=True,
spec_res=self.spectrum[key][3],
spec_res=spec_value[3],
)

self.diskspec.append(readmodel)
Expand Down Expand Up @@ -1087,26 +1106,25 @@ def _prior_transform(
Cube with the sampled model parameters.
"""

n_modelpar = len(self.modelpar)

if isinstance(cube, np.ndarray):
# Create new array for UltraNest
param_out = cube.copy()
else:
# Convert from ctypes.c_double to np.ndarray
# Only required with MultiNest
# n_modelpar = len(self.modelpar)
# cube = np.ctypeslib.as_array(cube, shape=(n_modelpar,))

# Use the same object with MultiNest
param_out = cube

for item in cube_index:
if item in self.prior:
if item in self.normal_prior:
# Gaussian prior
param_out[cube_index[item]] = stats.norm.ppf(
param_out[cube_index[item]],
loc=self.prior[item][0],
scale=self.prior[item][1],
loc=self.normal_prior[item][0],
scale=self.normal_prior[item][1],
)

else:
Expand Down Expand Up @@ -1341,7 +1359,7 @@ def _lnlike_func(

ln_like = 0.0

for key, value in self.prior.items():
for key, value in self.normal_prior.items():
if key == "mass":
if "logg" in self.modelpar:
mass = logg_to_mass(
Expand Down Expand Up @@ -1827,7 +1845,7 @@ def run_multinest(
tag: str,
n_live_points: int = 1000,
output: str = "multinest/",
prior: Optional[Dict[str, Tuple[float, float]]] = None,
**kwargs,
) -> None:
"""
Function to run the ``PyMultiNest`` wrapper of the
Expand Down Expand Up @@ -1856,14 +1874,6 @@ def run_multinest(
Number of live points.
output : str
Path that is used for the output files from MultiNest.
prior : dict(str, tuple(float, float)), None
Dictionary with Gaussian priors for one or multiple
parameters. The prior can be set for any of the
atmosphere or calibration parameters, e.g.
``prior={'teff': (1200., 100.)}``. Additionally, a prior
can be set for the mass, e.g. ``prior={'mass': (13., 3.)}``
for an expected mass of 13 Mjup with an uncertainty of
3 Mjup. The parameter is not used if set to ``None``.
Returns
-------
Expand All @@ -1875,10 +1885,19 @@ def run_multinest(

# Set attributes

if prior is None:
self.prior = {}
else:
self.prior = prior
if "prior" in kwargs:
warnings.warn(
"The 'prior' parameter has been deprecated "
"and will be removed in a future release. "
"Please use the 'normal_prior' of FitModel "
"instead.",
DeprecationWarning,
)

if kwargs["prior"] is None:
self.normal_prior = {}
else:
self.normal_prior = kwargs["prior"]

# Get the MPI rank of the process

Expand All @@ -1898,7 +1917,7 @@ def run_multinest(
# Add parallax to dictionary with Gaussian priors

if "parallax" in self.modelpar:
self.prior["parallax"] = (self.parallax[0], self.parallax[1])
self.normal_prior["parallax"] = (self.parallax[0], self.parallax[1])

@typechecked
def _lnprior_multinest(cube, n_dim: int, n_param: int) -> None:
Expand Down Expand Up @@ -2040,7 +2059,8 @@ def _lnlike_multinest(
# Add samples to the database

if mpi_rank == 0:
# Writing the samples to the database is only possible when using a single process
# Writing the samples to the database is only
# possible when using a single process
from species.data.database import Database

species_db = Database()
Expand All @@ -2057,16 +2077,12 @@ def _lnlike_multinest(

@typechecked
def run_ultranest(
self,
tag: str,
min_num_live_points=400,
output: str = "ultranest/",
prior: Optional[Dict[str, Tuple[float, float]]] = None,
self, tag: str, min_num_live_points=400, output: str = "ultranest/", **kwargs
) -> None:
"""
Function to run ``UltraNest`` for constructing the posterior
probability distributions on model parameters and computing
the marginal likelihood (i.e. "evidence").
Function to run ``UltraNest`` for estimating the posterior
distributions of model parameters and computing the
marginalized likelihood (i.e. "evidence").
Parameters
----------
Expand All @@ -2083,14 +2099,6 @@ def run_ultranest(
Therefore, a value above 100 is typically useful.
output : str
Path that is used for the output files from ``UltraNest``.
prior : dict(str, tuple(float, float)), None
Dictionary with Gaussian priors for one or multiple
parameters. The prior can be set for any of the
atmosphere or calibration parameters, e.g.
``prior={'teff': (1200., 100.)}``. Additionally, a prior
can be set for the mass, e.g. ``prior={'mass': (13., 3.)}``
for an expected mass of 13 Mjup with an uncertainty of
3 Mjup. The parameter is not used if set to ``None``.
Returns
-------
Expand All @@ -2102,10 +2110,19 @@ def run_ultranest(

# Set attributes

if prior is None:
self.prior = {}
else:
self.prior = prior
if "prior" in kwargs:
warnings.warn(
"The 'prior' parameter has been deprecated "
"and will be removed in a future release. "
"Please use the 'normal_prior' of FitModel "
"instead.",
DeprecationWarning,
)

if kwargs["prior"] is None:
self.normal_prior = {}
else:
self.normal_prior = kwargs["prior"]

# Get the MPI rank of the process

Expand All @@ -2125,7 +2142,7 @@ def run_ultranest(
# Add parallax to dictionary with Gaussian priors

if "parallax" in self.modelpar:
self.prior["parallax"] = (self.parallax[0], self.parallax[1])
self.normal_prior["parallax"] = (self.parallax[0], self.parallax[1])

@typechecked
def _lnprior_ultranest(cube: np.ndarray) -> np.ndarray:
Expand All @@ -2150,7 +2167,8 @@ def _lnprior_ultranest(cube: np.ndarray) -> np.ndarray:
@typechecked
def _lnlike_ultranest(params: np.ndarray) -> Union[float, np.float64]:
"""
Function for returning the log-likelihood for the sampled parameter cube.
Function for returning the log-likelihood for the
sampled parameter cube.
Parameters
----------
Expand Down

0 comments on commit 936dd32

Please sign in to comment.