"""
Low-resolution likelihood calculations for atmospheric retrieval.
Manages spectral binning (emission and transmission), continuum estimation,
and the full low-resolution likelihood evaluation. Extracted from ModelData.
"""
import numpy as np
from petitRADTRANS import physical_constants as phys_const
from scipy.integrate import trapezoid
from scipy.interpolate import splrep, splev
import pyratbay.spectrum as ps
from GUIBRUSHR.General_Constants.FunctionsAndConstants.Constant_Variables import ConstantVariables
from GUIBRUSHR.General_Constants.FunctionsAndConstants.general_functions import (
convolve_resolution, estimate_continuum,
)
[docs]
class LikelihoodLR:
"""
Low-resolution likelihood evaluator for atmospheric retrieval.
Handles spectral binning to instrument resolution, emission/transmission
mode handling, offset corrections, and chi-square likelihood evaluation.
Attributes:
atmosphere: Atmosphere object with species, resolution, stellar data.
bestpars_data: Bestpars object with nfit.
model_type: Model type string ("Retrieval", "Model", etc.).
clight: Speed of light constant.
"""
[docs]
def __init__(self, atmosphere, bestpars_data, model_type, clight):
self.atmosphere = atmosphere
self.bestpars_data = bestpars_data
self.model_type = model_type
self.clight = clight
# Reason of the most recent low_resolution_lhood failure (None on
# success). Read by ModelData.lh_function_gib: "lr_model_nan".
self.failure_reason = None
def _bin_emission_lr(self, instrument, ind_planet, model_LR, wl_full_resolution_LR, i, rp):
"""
Compute the emission spectrum value for a single wavelength bin.
Calculates the planet-to-star flux ratio for one instrument bin, then
scales by the (Rp/R*)^2 geometric factor to obtain the observed eclipse
depth. The result is returned as a relative value (1 + eclipse_depth),
where 1.0 represents the stellar continuum.
Two methods are supported depending on ``stellar_spectrum_type_lr_mode``:
- ``"Integral"``: the precomputed stellar flux is a per-bin integrated
value; the planet flux is integrated via the trapezoidal rule and
divided by the stellar integral.
- ``"Planck"`` / ``"Phoenix"``: the stellar spectrum is sampled
point-by-point; the eclipse ratio is the mean of the point-wise
planet/star flux ratios within the bin.
Args:
instrument: Instrument object with lr_data attributes
(``lbottom``, ``ltop`` bin edges in nm).
ind_planet (np.ndarray): Indices into ``model_LR`` and
``wl_full_resolution_LR`` that fall inside this bin.
model_LR (np.ndarray): Full-resolution planet flux array [W/m^2/nm].
wl_full_resolution_LR (np.ndarray): Full-resolution wavelength array [nm].
i (int): Index of the current wavelength bin.
rp (float): Planet radius [cm].
Returns:
float: Relative emission spectrum value for the bin:
``1 + eclipse_depth``, where values > 1 indicate planetary
emission above the stellar continuum.
"""
# Get stellar flux for this instrument and bin
# Note: stellar_spectrum_LR[instrument.name] can be either:
# - A single integrated flux value (if stellar_spectrum_type_lr == "Integral")
# - An array of point values (if stellar_spectrum_type_lr == "Planck" or "Phoenix")
flux_star_instrument_in_bin = self.atmosphere.stellar_spectrum[instrument.name]
# Extract planet flux values within this wavelength bin
flux_planet_in_bin = model_LR[ind_planet]
# INTEGRAL METHOD: Integrate planet flux over bin and normalize by stellar flux
if self.atmosphere.stellar_spectrum_type_lr_mode == "Integral":
if len(ind_planet) == 0:
flux_planet_full_in_bin = 0
else:
if len(ind_planet) == 1:
# Create simple 2-point integration using bin boundaries
x = np.array([instrument.lr_data.lbottom[i], instrument.lr_data.ltop[i]])
y_temp = model_LR[ind_planet]
# Duplicate flux values at boundaries for integration
y = np.append(y_temp, y_temp)
else:
y = flux_planet_in_bin
x = wl_full_resolution_LR[ind_planet]
# Perform trapezoidal integration of planet flux over wavelength bin
flux_planet_full_in_bin = trapezoid(y, x)
# Calculate eclipse depth: (integrated planet flux) / (stellar flux)
eclipse_ratio_bin = flux_planet_full_in_bin / flux_star_instrument_in_bin[i]
# POINT-SAMPLED METHOD: Use flux ratios at individual wavelength points
else: # Planck or Phoenix stellar spectrum
# Direct ratio of planet to stellar flux
# stellar_spectrum_LR contains point-by-point stellar flux values
eclipse_ratio_bin = np.mean(flux_planet_in_bin / flux_star_instrument_in_bin[i])
# Scale eclipse depth by (Rp/R*)^2 to get observed contrast
# This converts from flux ratio to observable transit/eclipse depth
Fscaled = (
eclipse_ratio_bin *
(rp / (self.atmosphere.target.stellar_radius * phys_const.r_sun)) ** 2
)
# Store as relative spectrum: 1 + eclipse_depth
# (1.0 = stellar continuum, values > 1 indicate planetary emission)
return 1 + Fscaled
@staticmethod
def _bin_transmission_lr(instrument, model_after_preprocessing, depth_full_resolution_LR, ind_planet):
"""
Compute the transmission spectrum value for a single wavelength bin.
Returns the mean spectral value within the bin. Two sources are used
depending on whether the instrument expects flux-normalised spectra:
- ``flux_norm=True``: averages the continuum-normalised model
(``model_after_preprocessing``) produced by
``_preprocess_model_lr_transmission``.
- ``flux_norm=False``: averages the raw transit depth array
(``depth_full_resolution_LR``), which already encodes
``1 - transit_depth``.
Args:
instrument: Instrument object; ``lr_data.flux_norm`` selects the
averaging source.
model_after_preprocessing (np.ndarray): Continuum-normalised
spectrum after resolution convolution and continuum estimation.
depth_full_resolution_LR (np.ndarray): Full-resolution transit
depth array (``1 - transit_depth``).
ind_planet (np.ndarray): Indices of full-resolution points that
fall inside the current wavelength bin.
Returns:
float: Mean spectral value for the bin.
"""
if instrument.lr_data.flux_norm:
# Flux-normalised mode: average the continuum-divided residuals
# produced by _preprocess_model_lr_transmission
value = np.mean(model_after_preprocessing[ind_planet])
else:
# Raw transit depth mode: average depth values directly within the bin
# depth_full_resolution_LR already encodes (1 - transit_depth)
value = np.mean(depth_full_resolution_LR[ind_planet])
return value
def _preprocess_model_lr_transmission(self, instrument, wl_full_resolution_LR, depth_full_resolution_LR,
num_windows=40, top_window=10):
"""
Preprocess the transmission spectrum before binning to instrument resolution.
When the instrument expects flux-normalised spectra (``flux_norm=True``),
this method:
1. Convolves the model to the instrument's spectral resolution if the
model resolution is higher than the dataset resolution; otherwise
resamples via cubic spline onto the instrument wavelength grid.
2. Estimates and divides out the spectral continuum using a sliding-window
approach (``estimate_continuum``), returning the normalised residuals.
When ``flux_norm=False`` (emission or simple transit depth), no
processing is applied and the inputs are returned unchanged.
Args:
instrument: Instrument object with lr_data attributes:
``flux_norm`` (bool), ``resolution_dataset`` (float),
``ldata_LR`` (np.ndarray, wavelength grid in nm).
wl_full_resolution_LR (np.ndarray): Full-resolution wavelength
array [nm].
depth_full_resolution_LR (np.ndarray): Full-resolution transit
depth array.
num_windows (int): Number of sliding windows for continuum
estimation. Default is 40.
top_window (int): Number of top points per window used to anchor
the continuum. Default is 10.
Returns:
tuple:
- **wl_after_processing** (np.ndarray): Wavelength array after
preprocessing [nm].
- **model_after_processing** (np.ndarray): Spectrum after
continuum normalisation, or the raw depth array if
``flux_norm=False``.
"""
if instrument.lr_data.flux_norm:
if self.atmosphere.resolution_obj.LR_resolution > instrument.lr_data.resolution_dataset:
# Model resolution is higher than the dataset: convolve down to
# the instrument resolution using a Gaussian kernel defined by
# the half-width at half-maximum (hfhm) in velocity space
hfhm = ConstantVariables.CLIGHT / (2 * instrument.lr_data.resolution_dataset)
wl_new, spectrum_new = convolve_resolution(
wl_full_resolution_LR, depth_full_resolution_LR, hfhm,
self.atmosphere.rv_sampling
)
else:
# Model resolution is lower than or equal to the dataset:
# resample onto the instrument wavelength grid via cubic spline
wl_new = instrument.lr_data.ldata_LR
spline_model = splrep(wl_full_resolution_LR, depth_full_resolution_LR)
spectrum_new = splev(wl_new, spline_model)
wl_after_processing = wl_new
# Divide out the spectral continuum so that only relative absorption
# features remain (continuum-normalised spectrum)
_, model_after_processing = estimate_continuum(
spectrum_new, num_windows=num_windows, top_window=top_window, check_gaussian=False
)
else:
# Emission or raw transit depth: no resolution matching or continuum
# normalisation needed — pass arrays through unchanged
wl_after_processing = wl_full_resolution_LR
model_after_processing = depth_full_resolution_LR
return wl_after_processing, model_after_processing
def _bin_spectrum_to_instrument_resolution(
self,
instrument,
wl_full_resolution_LR,
model_LR,
depth_full_resolution_LR,
rp,
):
"""
Bin a full-resolution spectrum to the instrument's spectral resolution.
This method performs spectral binning by integrating the high-resolution
atmospheric model spectrum over the wavelength bins defined by the
instrument's resolution. It handles both emission and transmission modes,
with special treatment for emission spectroscopy including stellar flux
normalization and eclipse depth calculations.
The binning process:
1. For each wavelength bin defined by the instrument
2. Find all full-resolution wavelength points within the bin
3. For emission: calculate eclipse depth using planet/star flux ratio
4. For transmission: average the spectrum values in the bin
5. Handle edge cases (e.g., zero flux) with appropriate fallbacks
Args:
instrument: Instrument object containing wavelength bin definitions
(lr_data.lbottom, lr_data.ltop, lr_data.npix)
wl_full_resolution_LR: Full-resolution wavelength array [microns]
model_LR: Full-resolution model flux array (planet flux for emission,
transit depth for transmission)
depth_full_resolution_LR: Full-resolution spectrum array
rp: Planet radius in cm
Returns:
final_spectrum_LR: 1D array of binned spectrum values with length
equal to instrument.lr_data.npix. Values are set
to unity (1.0) where NaN occurs, representing a
neutral (featureless) spectrum.
Notes:
- For emission mode with integral stellar spectrum type, uses
trapezoidal integration over each bin
- For emission mode with point-sampled stellar spectrum (Planck/Phoenix),
uses point-by-point ratios
- NaN values in the output are replaced with 1.0 (neutral spectrum)
- The returned spectrum is in relative units (1.0 = continuum level)
"""
# --- Step 1: mode-dependent preprocessing ---
if self.atmosphere.rad_mode == "Emission":
# Emission: no continuum normalisation needed; use arrays as-is
wl_after_processing = wl_full_resolution_LR
model_after_processing = model_LR
else:
# Transmission: convolve/resample to instrument resolution and
# normalise to the spectral continuum when flux_norm is enabled
wl_after_processing, model_after_processing = self._preprocess_model_lr_transmission(
instrument, wl_full_resolution_LR, depth_full_resolution_LR
)
# --- Step 2: apply systemic velocity shift ---
# Doppler-shift the wavelength axis by the target's systemic velocity
# so that the model is aligned with the observed rest frame
vsys = self.atmosphere.target.systemic_velocity
wl_after_processing = wl_after_processing * (1 + vsys / self.clight)
# --- Step 3: binning strategy selection ---
# use_hr_linelists_for_lr / temporary_force_new_binning activate the
# faster vectorised path (spline or top-hat via pysynphot); the default
# path loops bin-by-bin using the dedicated helper methods
if self.atmosphere.use_hr_linelists_for_lr or self.atmosphere.temporary_force_new_binning:
if instrument.lr_data.use_spline_instead_of_tophat:
# SPLINE PATH: fit a cubic spline to the full-resolution model
# and evaluate it at the instrument wavelength grid points
spline_model = splrep(wl_after_processing, model_after_processing)
final_spectrum_LR = splev(instrument.lr_data.ldata_LR, spline_model)
# TODO: see if emission is necessary. If yes, adapt the star flux to the spline
# if self.atmosphere.rad_mode == "Emission":
else:
# TOP-HAT PATH: bin the spectrum using rectangular top-hat
# kernels centred on each instrument channel (ps.bin_spectrum
# expects wavelengths in microns, hence the /1e3 conversion)
final_spectrum_LR = ps.bin_spectrum(
instrument.lr_data.ldata_LR_micron, # channel centres [µm]
wl_after_processing / 1e3, # model wavelengths [µm]
model_after_processing, # model flux / depth
half_widths=instrument.lr_data.step_LR_micron, # half bin widths [µm]
gaps="ignore" # skip gaps between channels
)
if self.atmosphere.rad_mode == "Emission":
# Normalise binned planet flux by the stellar flux to get
# the planet-to-star flux ratio per channel
eclipse_ratio_bin = final_spectrum_LR / self.atmosphere.stellar_spectrum[instrument.name]
# Scale by (Rp/R*)^2 to convert to observable eclipse depth
spectrum_binned_LR_scaled = (
eclipse_ratio_bin *
(rp / (self.atmosphere.target.stellar_radius * phys_const.r_sun)) ** 2
)
# Express as relative spectrum: 1 = stellar continuum
final_spectrum_LR = 1 + spectrum_binned_LR_scaled
else:
# BIN-BY-BIN PATH: explicit loop over each instrument wavelength bin
# Initialise output array with unity (featureless continuum baseline)
final_spectrum_LR = np.ones(instrument.lr_data.npix)
# Loop over bin edges defined by the instrument (lbottom, ltop)
for i, (bottom, top) in enumerate(zip(instrument.lr_data.lbottom, instrument.lr_data.ltop)):
# Select full-resolution points inside the open-closed interval (bottom, top]
ind_planet = np.where(
(wl_after_processing > bottom) & (wl_after_processing <= top)
)[0]
if self.atmosphere.rad_mode == "Emission":
# Emission: compute planet/star flux ratio and scale by (Rp/R*)^2
final_spectrum_LR[i] = self._bin_emission_lr(
instrument, ind_planet, model_after_processing, wl_after_processing, i, rp
)
else:
# Transmission: average the (normalised or raw) depth values in the bin
final_spectrum_LR[i] = self._bin_transmission_lr(
instrument, model_after_processing, depth_full_resolution_LR, ind_planet
)
# --- Step 4: sanitise output ---
# Replace NaN entries (e.g. from empty bins or 0/0 divisions) with 1.0,
# which represents a neutral, featureless spectrum at the continuum level
final_spectrum_LR[np.isnan(final_spectrum_LR)] = 1
return final_spectrum_LR
[docs]
def low_resolution_lhood(
self,
temperature,
mass_fraction,
vmr,
MMW,
dict_calc_model,
rp,
):
"""
Compute likelihood for low-resolution observations.
This method calculates the likelihood by comparing the binned theoretical
spectrum with low-resolution observational data. It handles spectral
binning, offset corrections, and chi-square calculation for each instrument.
Args:
temperature: Temperature profile
mass_fraction: Mass fraction profiles for species
vmr: Volume mixing ratios
MMW: Mean molecular weight profile
dict_calc_model: Dictionary containing model calculation parameters
rp: Planet radius in cm
Returns:
Tuple containing:
- lhood: Log-likelihood value
- wl_full_resolution_LR: Full resolution wavelength array
- depth_full_resolution_LR: Full resolution spectrum
- wl_binned_LR_complete: Binned wavelength arrays for all instruments
- spectrum_binned_LR_complete: Binned spectra for all instruments
- opacity_contribution_LR: Opacity contributions (if model plotting)
- lh_low_resolution: Dictionary of per-instrument likelihood details
"""
lhood = 0
opacity_contribution_LR = None
wl_binned_LR_complete = []
spectrum_binned_LR_complete = []
self.failure_reason = None # cleared on every evaluation
# Calculate full resolution model spectrum
wl_full_resolution_LR, model_LR, depth_full_resolution_LR = self.atmosphere.calc_model(
temperature, mass_fraction, vmr, MMW, dict_calc_model, False,
)
# MODEL FAILURE GUARD: stop here if the LR model spectrum is NaN/Inf.
# This is a model-evaluation failure upstream of the binning/likelihood;
# returning now (before the binning step, which would silently mask NaN
# to 1.0) makes the cause explicit and distinct from any HR/PCA issue.
# depth_full_resolution_LR is populated for both engines; model_LR is
# None under PyratBay, so it is guarded conditionally.
if (not np.all(np.isfinite(depth_full_resolution_LR))
or (model_LR is not None and not np.all(np.isfinite(model_LR)))):
self.failure_reason = "lr_model_nan"
print("Low-resolution model failed: NaN/Inf in the model spectrum "
"(model evaluation, not PCA).")
return (
-np.inf, wl_full_resolution_LR, depth_full_resolution_LR, [], [], None, {}
)
# Calculate opacity contributions for model plotting
if "Opacity" in self.model_type:
opacity_contribution_LR = self.atmosphere.opacity_contribution(
temperature, mass_fraction, MMW, dict_calc_model, False
)
lh_low_resolution = {}
# Process each low-resolution instrument
for index_LR, instrument in enumerate(self.atmosphere.resolution_obj.instruments_LR):
# Bin the full-resolution spectrum to instrument resolution
final_spectrum_LR = self._bin_spectrum_to_instrument_resolution(
instrument=instrument,
wl_full_resolution_LR=wl_full_resolution_LR,
model_LR=model_LR,
depth_full_resolution_LR=depth_full_resolution_LR,
rp=rp
)
# Apply instrument-specific offset correction
offsetLR_arr = dict_calc_model["offsetLR_arr"]
offsetLR = 0 if (index_LR == 0 or offsetLR_arr is None) else offsetLR_arr[index_LR - 1]
final_spectrum_LR += offsetLR
# Store results for output
wl_binned_LR_complete.append(instrument.lr_data.ldata_LR)
spectrum_binned_LR_complete.append(final_spectrum_LR)
# Calculate chi-square and likelihood contribution
# Compute residuals
residuals = instrument.lr_data.rdata_LR - final_spectrum_LR
# Select error array element-wise based on the sign of each residual:
# - positive residual → upper error (sigma_plus)
# - negative residual → lower error (sigma_minus)
sigma_used = np.where(
residuals >= 0,
instrument.lr_data.errdata_LR_plus,
instrument.lr_data.errdata_LR_minus
)
# Per-instrument noise scale (beta_LR). One value per LR instrument; the
# error bars are rescaled multiplicatively (sigma -> beta_LR * sigma), with
# the matching -N ln(beta_LR) normalization term below. Defaults to 1 when
# beta_LR is not fitted. The sigma sign-selection above uses the raw
# residuals, so the positive beta_LR scaling does not change which branch
# (sigma_plus / sigma_minus) is picked.
beta_LR_arr = dict_calc_model.get("beta_LR_arr")
beta_LR = (
beta_LR_arr[index_LR]
if beta_LR_arr is not None and np.size(beta_LR_arr) > 0
else 1
)
chiquadro_LR = np.sum((residuals / (beta_LR * sigma_used)) ** 2)
dof = instrument.lr_data.npix
chi2_reduced = chiquadro_LR / dof
# Gaussian log-likelihood (same form as the HR likelihood, Gibson et al. 2021):
# ln L = -(N/2) ln(2pi) - N ln(beta) - sum(ln sigma_i) - chi2 / 2
# chiquadro_LR is already computed with beta-scaled errors.
lh_single_instrument = (
-0.5 * dof * np.log(2.0 * np.pi)
- dof * np.log(beta_LR)
- np.sum(np.log(sigma_used))
- 0.5 * chiquadro_LR
)
lhood += lh_single_instrument
lh_low_resolution[
instrument.name] = f"{chiquadro_LR};{lh_single_instrument};{chi2_reduced};{instrument.lr_data.npix};{self.bestpars_data.nfit}"
return (
lhood, wl_full_resolution_LR, depth_full_resolution_LR, wl_binned_LR_complete, spectrum_binned_LR_complete,
opacity_contribution_LR, lh_low_resolution
)