petitRADTRANS.radtrans_core.spectrum_utils

petitRADTRANS.radtrans_core.spectrum_utils#

Attributes#

DEFAULT_BREADS_FILTER_SIZE

_DEFAULT_CONV_HALF_WINDOW

Functions#

compute_bins_edges(→ jax.typing.ArrayLike)

Calculate bin edges for middle bin points.

continuum_subtract_spectrum(→ jax.typing.ArrayLike)

Apply a continuum-removal filter to a spectrum.

convolve_constant_resolution(→ jax.typing.ArrayLike)

convolve_variable_resolution(input_wavelength, ...[, ...])

Convolve a spectrum with a variable-width Gaussian kernel.

convolve(→ jax.typing.ArrayLike)

Convolve a model spectrum with a given instrument resolution.

_convolve_running(input_spectrum, ...)

Convolve a spectrum to a new resolving power.

convolve_and_sample_variable_resolution_breads(...[, ...])

Simulate the observations of a model.

filter_spectrum_with_spline(wavelengths, fluxes[, ...])

Use spline interpolation produce a high-pass filtered version of the input spectrum.

linear_spline_interpolation(x_knots, x_samples[, ...])

From BREADS, BSD 3-Clause License

median_filter_spectrum(spectrum[, size, mode])

High-pass filter spectrum using median filtering.

rebinning_interpolation(input_wavelengths, ...)

_raise_small_size(x)

_raise_incomplete_wavelength_range(args)

rebin_spectrum(input_wavelengths, input_spectrum, ...)

rebin_spectrum_bin(input_wavelengths, input_spectrum, ...)

savgol_filter_spectrum(spectrum[, window_length, ...])

High-pass filter spectrum using Savitzky-Golay filtering.

Module Contents#

petitRADTRANS.radtrans_core.spectrum_utils.DEFAULT_BREADS_FILTER_SIZE = 201#
petitRADTRANS.radtrans_core.spectrum_utils.compute_bins_edges(middle_bin_points: jax.typing.ArrayLike) → jax.typing.ArrayLike#

Calculate bin edges for middle bin points.

Args:

middle_bin_points: array of size N containing the middle bin points to calculate the bin edges of

Returns:

array of size N+1 containing the bins edges

petitRADTRANS.radtrans_core.spectrum_utils.continuum_subtract_spectrum(wavelengths: jax.typing.ArrayLike, spectrum: jax.typing.ArrayLike, parameters: dict, method: str = 'spline', name: str = None) → jax.typing.ArrayLike#

Apply a continuum-removal filter to a spectrum.

The function dispatches to one of three continuum-subtraction methods, selected by method:

  • "spline": fit a smooth spline continuum and return the residuals.

    The spline node locations come from parameters using the first matching key below: - f"{name}_nodes": number of evenly spaced spline nodes across the

    wavelength range.

    • f"{name}_node_array": explicit spline node coordinates.

    • "continuum_subtract_nodes": fallback number of evenly spaced spline

      nodes across the wavelength range.

    At least one of these keys must be present or a ValueError is raised.

  • "medfilt": estimate the continuum with a median filter and subtract

    it from the spectrum. The filter width is taken from the first matching key below, defaulting to 10 if neither is present: - f"{name}_median_filter_width" - "median_filter_width"

  • "savgol": estimate the continuum with a Savitzky-Golay filter and

    subtract it from the spectrum. The filter settings are read from the first matching keys below, with defaults window_length=10 and polyorder=3: - f"{name}_savgol_window_length" or "savgol_window_length" - f"{name}_savgol_polyorder" or "savgol_polyorder"

Values in parameters may be provided directly or as objects exposing a value attribute.

Args:

wavelengths: Wavelength grid of the input spectrum. spectrum: Spectrum to continuum-subtract. parameters: Parameter mapping used to configure the selected method. method: Continuum-subtraction method. Supported values are

"spline", "medfilt", and "savgol".

name: Optional prefix for method-specific parameter keys.

Returns:

Continuum-subtracted spectrum on the original wavelength grid.

petitRADTRANS.radtrans_core.spectrum_utils.convolve_constant_resolution(input_wavelength: jax.typing.ArrayLike, input_flux: jax.typing.ArrayLike, instrument_res: any, max_radius: int = 20) → jax.typing.ArrayLike#
petitRADTRANS.radtrans_core.spectrum_utils._DEFAULT_CONV_HALF_WINDOW = 128#
petitRADTRANS.radtrans_core.spectrum_utils.convolve_variable_resolution(input_wavelength, input_flux, resolutions, half_window=None)#

Convolve a spectrum with a variable-width Gaussian kernel. Based on implementation from Ben Burningham for Brewster. fwang23nh/brewster_v2

Memory-bounded sliding-window implementation: the original built the full dense (N x N) Gaussian-LSF matrix over the entire input grid, which is O(N_model^2) in memory and made reverse-mode AD (NUTS gradients) allocate hundreds of GB at line-by-line model resolution. Because a Gaussian LSF has finite support, each output channel only needs the model points within a few sigma of it, so we evaluate the kernel on a fixed +/-half_window slice of the grid centred on each channel. This is O(N_model * (2*half_window + 1)) and returns the same result as the full-matrix version to within the negligible far-wing contribution (window edges carry weight ~exp(-half_window^2/…) ~ 0). Output grid == input grid, so the surrounding RV-shift / rebin pipeline is unchanged.

Args:

input_wavelength - units of micron (the model grid; also the output grid) input_flux - arbitrary units resolutions - spectral resolving power, R, per input wavelength half_window - number of grid points on each side of a channel to include

in the kernel. Defaults to _DEFAULT_CONV_HALF_WINDOW.

Returns:

convolved_spectrum - same units / grid as input_flux

petitRADTRANS.radtrans_core.spectrum_utils.convolve(input_wavelength: jax.typing.ArrayLike, input_flux: jax.typing.ArrayLike, instrument_res: any) → jax.typing.ArrayLike#

Convolve a model spectrum with a given instrument resolution. This method replicates the logic from petitRADTRANS.retrieval.data.Data.convolve().

Args:
input_wavelengthArrayLike

The wavelength array of the model spectrum.

input_fluxArrayLike

The flux array of the model spectrum.

instrument_resfloat or ArrayLike

The instrument resolution. Can be a single float (constant resolution) or an array (variable resolution per wavelength point).

Returns:
flux_lsfArrayLike

The convolved flux array.

petitRADTRANS.radtrans_core.spectrum_utils._convolve_running(input_spectrum, sigma_lsf_gauss_filter, **kwargs)#

Convolve a spectrum to a new resolving power. The spectrum is convolved using Gaussian filters with a standard deviation

std_dev = R_in(lambda) / R_new(lambda) * input_wavelengths_bins.

Both the input resolving power and output resolving power can vary with wavelength. The input resolving power is given by:

lambda / Delta_lambda

where lambda is the center of a wavelength bin and Delta_lambda is the difference between the edges of the bin.

The weights of the convolution are stored in a (N, M) matrix, with N being the size of the input, and M the size of the convolution kernels. To speed-up calculations, a matrix A of shape (N, M) is built from the inputs such as:

A[i, :] = s[i - M/2], s[i - M/2 + 1], …, s[i - M/2 + M],

with s the input spectrum. The definition of the convolution C of s by constant weights with wavelength is:

C[i] = sum_{j=0}^{j=M-1} s[i - M/2 + j] * weights[j].

Thus, the convolution of s by weights at index i is:

C[i] = sum_{j=0}^{j=M-1} A[i, j] * weights[i, j].

Args:

input_wavelengths: (cm) wavelengths of the input spectrum input_spectrum: input spectrum convolve_resolving_power: resolving power of output spectrum input_resolving_power: if not None, skip its calculation using input_wavelengths

Returns:

convolved_spectrum: the convolved spectrum at the new resolving power

petitRADTRANS.radtrans_core.spectrum_utils.convolve_and_sample_variable_resolution_breads(wavelengths, resolutions, model_wavelengths, model_fluxes, num_sigma=3, filter_size=None)#

Simulate the observations of a model. Convolves the model with a variable Gaussian LSF, sampled at each desired spectral channel.

From Jerry Xuan circa 2024, jruffio/breads

Args:
wavelengths: jnp.ndarray

the wavelengths desired (length of N_output)

resolutions: np.ndarray

the spectral resolving power at each wavelengths (length of N_output)

model_wavelengths: np.ndarray

the wavelengths of the model (length of N_model)

model_fluxes: np.ndarray

the fluxes of the model (length of N_model)

num_sigma: float

number of +/- sigmas to evaluate the LSF to.

filter_size: int or None

Number of points in the LSF filter. Must be provided as a static Python int when this function is JIT-compiled. If None, a conservative static default is used so the function remains JIT-compatible.

Returns:
output_model: np.ndarray

the fluxes in each of the wavelength channels (length of N_output)

petitRADTRANS.radtrans_core.spectrum_utils.filter_spectrum_with_spline(wavelengths, fluxes, uncertainties=None, x_nodes=None, m_spline=None)#

Use spline interpolation produce a high-pass filtered version of the input spectrum.

From BREADS, BSD 3-Clause License Copyright (c) 2024, jruffio jruffio/breads

Args:
wavelengths: np.ndarray

Wavelengths of the input spectrum

fluxes: np.ndarray

Fluxes of the input spectrum

uncertainties: np.ndarray

Uncertainties of the input spectrum

x_nodes: np.ndarray

Nodes for the spline interpolation as np.ndarray in the same units as wavelengths. x_nodes can also be a list of ndarrays/list to model discontinous functions.

m_spline: np.ndarray

Pre-computed spline matrix. If None, it will be computed from x_nodes and wavelengths.

Returns:
high_pass_filtered_fluxes: np.ndarray

High-pass filtered fluxes of the input spectrum

petitRADTRANS.radtrans_core.spectrum_utils.linear_spline_interpolation(x_knots, x_samples, spline_degree=3)#

From BREADS, BSD 3-Clause License

Copyright (c) 2024, jruffio

Compute a spline based linear model. If Y=[y1,y2,..] are the values of the function at the location of the node [x1,x2,…]. jnp.dot(M,Y) is the interpolated spline corresponding to the sampling of the x-axis (x_samples)

Uses interpax (JAX-native) for spline basis evaluation instead of SciPy’s InterpolatedUnivariateSpline. This enables JIT compilation and GPU acceleration. Note: interpax uses natural cubic splines while the original SciPy implementation used B-splines, resulting in different basis functions. However, the practical impact on spectrum fitting is negligible (~0.01 mean difference in fitted spectra).

Args:
x_knots: List of nodes for the spline interpolation as jnp.ndarray in the same units as x_samples.

x_knots can also be a list of ndarrays/list to model discontinous functions.

x_samples: Vector of x values. ie, the sampling of the data. spline_degree: Degree of the spline interpolation (default: 3).

if jnp.size(x_knots) <= spline_degree, then spline_degree = jnp.size(x_knots)-1

Returns:

M: Matrix of size (D,N) with D the size of x_samples and N the total number of nodes.

petitRADTRANS.radtrans_core.spectrum_utils.median_filter_spectrum(spectrum, size=10, mode='nearest')#

High-pass filter spectrum using median filtering. Returns the residuals after subtracting the median-filtered (low-pass) version.

Args:

spectrum: Input spectrum array size: Median filter window size mode: Padding mode for edges (‘nearest’, ‘reflect’, ‘constant’, etc.)

Returns:

High-pass filtered spectrum (original - median_filtered)

petitRADTRANS.radtrans_core.spectrum_utils.rebinning_interpolation(input_wavelengths, input_spectrum, rebin_bin_low, rebin_bin_high, n_rebin)#
petitRADTRANS.radtrans_core.spectrum_utils._raise_small_size(x)#
petitRADTRANS.radtrans_core.spectrum_utils._raise_incomplete_wavelength_range(args)#
petitRADTRANS.radtrans_core.spectrum_utils.rebin_spectrum(input_wavelengths, input_spectrum, rebinned_wavelengths)#
petitRADTRANS.radtrans_core.spectrum_utils.rebin_spectrum_bin(input_wavelengths, input_spectrum, rebinned_wavelengths, bin_widths)#
petitRADTRANS.radtrans_core.spectrum_utils.savgol_filter_spectrum(spectrum, window_length: int = 11, polyorder: int = 3, mode: str = 'nearest')#

High-pass filter spectrum using Savitzky-Golay filtering. Returns the residuals after subtracting the smoothed (low-pass) version.

Args:

spectrum: Input spectrum array window_length: Window size for polynomial fitting (will be made odd) polyorder: Order of polynomial to fit (must be < window_length) mode: Padding mode for edges (‘nearest’, ‘reflect’, ‘constant’, etc.)

Returns:

High-pass filtered spectrum (original - smoothed)