Results

Output parsing: parameter tables, spectra tables, and FITS HDU lists.

These functions transform raw posterior samples + ModelArgs into user-friendly Table objects and FITS files.

unite.results.count_parameters(model_fn, model_args)[source]

Count the number of free scalar parameters (degrees of freedom) in the model.

Traces the model with a dummy PRNG key and counts every latent (non-observed) sample site, summing the sizes of all their shapes. This gives the total number of unconstrained scalar parameters — i.e. the model degrees of freedom.

Return type:

int

Parameters:
model_fncallable

The numpyro model function returned by build().

model_argsModelArgs

The model arguments returned by build().

Returns:
int

Total number of free scalar parameters.

Examples

>>> model_fn, model_args = builder.build()
>>> print(f'Free parameters: {count_parameters(model_fn, model_args)}')
Free parameters: 14
unite.results.make_parameter_table(samples, args, *, percentiles=None)[source]

Build an Astropy table of posterior parameter samples in physical units.

Return type:

QTable

Parameters:
samplesdict of str to ndarray

Posterior samples in physical space. When using ModelBuilder.fit(), samples are already transformed. When calling mcmc.get_samples() directly, first pass through transform_reparam_samples() to convert any reparameterized (unit-space) parameters back to physical values.

argsModelArgs

Model arguments from ModelBuilder.build().

percentilesndarray of float or None

Array of percentile values in range (0, 1), e.g. [0.16, 0.5, 0.84]. If provided, returns one row per percentile with percentile values. If None (default), returns one row per posterior sample.

Returns:
astropy.table.QTable

If percentiles is None: one row per posterior sample. If percentiles is provided: one row per percentile with a percentile column and one column per parameter. Columns carry physical units where known:

  • Line flux parameters — flux_unit * canonical_wl_unit

  • FWHM parameters — km/s

  • Continuum scale parameters — flux_unit

  • Continuum slope / polynomial coefficients — flux_unit / wl_unit^n

  • Shape / index parameters (beta, temperature, …) — raw values

  • Rest equivalent width columns (rew_{line_label}) — canonical_wl_unit. One column per line whose rest-frame wavelength falls within a continuum region. Appended after all model parameters when a continuum is present.

Parameters:

Notes

For absorption lines, the rest equivalent width is computed by numerically integrating the absorbed flux profile over the spectrum with the finest pixel grid covering the line. Use absorption REW values with caution when the covering spectrum does not fully resolve the absorption profile.

unite.results.make_spectra_tables(samples, args, *, insert_nan=False, percentiles=None, return_hdul=False)[source]
Overloads:

  • samples (dict[str, np.ndarray]), args (ModelArgs), insert_nan (bool), percentiles (np.ndarray | None), return_hdul (Literal[False]) → dict[str, Table]

  • samples (dict[str, np.ndarray]), args (ModelArgs), insert_nan (bool), percentiles (np.ndarray | None), return_hdul (Literal[True]) → fits.HDUList

Parameters:
Return type:

dict[str, Table] | HDUList

Build per-spectrum tables of model decompositions.

Parameters:
samplesdict of str to ndarray

Posterior samples.

argsModelArgs

Model arguments from ModelBuilder.build().

insert_nanbool

If True, insert one NaN row at the midpoint wavelength between each pair of consecutive continuum regions. Default False.

percentilesndarray of float or None

Array of percentile values in range (0, 1), e.g. [0.16, 0.5, 0.84]. If provided, collapses the sample dimension to those percentiles (shape (n_percentiles, n_pixels)). If None (default), returns all samples (shape (n_samples, n_pixels)).

return_hdulbool

If True, wrap the per-spectrum tables in an HDUList and return that instead of a dict. HDU 0 is an empty PrimaryHDU; subsequent HDUs are BinTableHDU entries whose extension names are the spectrum names (upper-cased for FITS compatibility). Default False.

Returns:
dict of str to astropy.table.QTable, or astropy.io.fits.HDUList

When return_hdul=False (default): a dict keyed by spectrum name, one table per spectrum. When return_hdul=True: an HDUList with HDU 0 empty and one BinTableHDU per spectrum. In both cases columns carry physical units where flux_unit was set on the spectrum.

Parameters:
Return type:

dict[str, Table] | HDUList

unite.results.make_hdul(samples, args, *, insert_nan=False, percentiles=None)[source]

Build a FITS HDU list from posterior samples.

Return type:

HDUList

Parameters:
samplesdict of str to ndarray

Posterior samples.

argsModelArgs

Model arguments from ModelBuilder.build().

insert_nanbool

Insert NaN rows between continuum regions. Default False.

percentilesndarray of float or None

Array of percentile values in range (0, 1). If provided, output tables contain percentile rows/columns. If None (default), output tables contain all samples.

Returns:
astropy.io.fits.HDUList

HDU 0: PrimaryHDU (empty, metadata in header). HDU 1: BinTableHDU from parameter table. HDU 2+: BinTableHDU per spectrum.

Parameters:
unite.results.freeze_from_samples(samples, args, *, cenfunc='median')[source]

Convert posterior samples to Fixed priors.

For each parameter in args, wraps its central value in a Fixed prior. Free parameters are summarised by cenfunc applied to their posterior samples. Parameters that were already Fixed in the original fit (such as norm_wav_a) are passed through with their constant value, making the returned dict a complete set of frozen priors – no separate lookup of region centres or other fixed quantities is needed.

Intended for re-fitting the same observation with a different model (e.g. adding or removing a line component while holding kinematics fixed). Always reuse the same Spectra object between fits so that line_scale and continuum_scale are identical and the frozen amplitude values remain physically consistent.

Return type:

dict[str, Fixed]

Parameters:
samplesdict

Posterior samples, e.g. from mcmc.get_samples().

argsModelArgs

Model arguments from build() for the original fit.

cenfuncstr or callable, optional

Central-value function. Default 'median'.

If a callable, it is called with a 1-D numpy.ndarray of posterior samples and must return a scalar.

If a string, recognised presets are:

'median'

Coordinate-wise posterior median.

'mean'

Coordinate-wise posterior mean.

'map'

Maximum-a-posteriori sample – the single sample with the highest log-joint samples['log_prob']. Requires 'log_prob' in samples (produced automatically by ModelBuilder.fit()).

'mle'

Maximum-likelihood sample – the single sample with the highest total log-likelihood samples['log_likelihood']. Requires 'log_likelihood' in samples (produced automatically by ModelBuilder.fit()).

Returns:
dict[str, Fixed]

Mapping from numpyro site name to Fixed prior. Keys cover all parameters in args.dependency_order: free parameters frozen at the chosen central value, already-Fixed parameters passed through at their constant value.

Raises:
ValueError

If cenfunc='map' is requested but 'log_prob' is not in samples.

ValueError

If cenfunc='mle' is requested but 'log_likelihood' is not in samples.

ValueError

If cenfunc is an unrecognised string.

Parameters:

Examples

Freeze all parameters at their posterior median for a re-fit on the same spectra (e.g. to add a component while holding kinematics fixed). norm_wav_a is automatically included – no manual lookup needed:

>>> frozen = freeze_from_samples(samples, args)
>>> fwhm_narrow.prior = frozen['fwhm_gauss_narrow']
>>> norm_wav_token.prior = frozen['norm_wav_a']  # was already Fixed

Use the MAP sample (highest log posterior) – more robust when scale and slope parameters are correlated. Requires samples to contain a 'log_prob' key (produced automatically by ModelBuilder.fit()):

>>> frozen_map = freeze_from_samples(samples, args, cenfunc='map')

Use the MLE sample (highest total log-likelihood):

>>> frozen_mle = freeze_from_samples(samples, args, cenfunc='mle')

Use the mean instead of the median:

>>> frozen = freeze_from_samples(samples, args, cenfunc='mean')

Pass a callable for custom logic (e.g. a high percentile):

>>> import numpy as np
>>> frozen = freeze_from_samples(samples, args, cenfunc=lambda x: np.percentile(x, 75))

Model evaluator: decompose posterior predictions into per-line and per-region contributions.

Given posterior samples and ModelArgs, this module reconstructs the full model prediction for each spectrum, broken down into individual line and continuum-region contributions in original flux units.

class unite.compute.SpectrumPrediction(wavelength, total, lines, continuum_regions, tau_profiles)[source]

Bases: object

Decomposed model prediction for a single spectrum.

All arrays are in original (un-normalized) flux units.

For emission lines, each entry in lines is the intrinsic (un-attenuated) flux profile: flux * profile. Summing all line contributions and continuum regions always reconstructs total regardless of zorder configuration.

For absorption lines, each entry in lines is the flux removed by that absorber (negative): total - total_without_j.

Parameters:
wavelength: ndarray

Pixel-center wavelengths in the disperser’s unit. Shape (n_pixels,).

total: ndarray

Total model flux (lines + continuum). Shape (n_samples, n_pixels).

lines: dict[str, ndarray]

Per-line contributions keyed by informative line labels (e.g. 'Ha', '[NII]_6585'). For emission lines: intrinsic (un-attenuated) flux profile (positive). For absorption lines: flux removed by the absorber (negative). Shape (n_samples, n_pixels) each.

continuum_regions: dict[str, ndarray]

Per-continuum-region contributions keyed by informative region labels (e.g. 'linear_6400_6700', 'powerlaw_0.95_2.5'). Shape (n_samples, n_pixels) each.

tau_profiles: dict[str, ndarray]

LSF-convolved optical depth profiles tau_j * phi_j(λ) for absorption lines, evaluated at pixel midpoints. Dimensionless and non-negative. Keyed by line label; empty dict for emission-only models. Shape (n_samples, n_pixels) each.

unite.compute.evaluate_model(samples, args)[source]

Evaluate the model for each posterior sample and decompose contributions.

Uses jax.vmap() to evaluate all samples in a single vectorised XLA kernel launch rather than a Python loop, giving a large speed-up when the number of posterior samples is large.

Return type:

list[SpectrumPrediction]

Parameters:
samplesdict of str to ndarray

Posterior samples as returned by mcmc.get_samples() or Predictive. Each value has shape (n_samples,) or (n_samples, ...).

argsModelArgs

Pre-built data bundle from ModelBuilder.build().

Returns:
list of SpectrumPrediction

One prediction per spectrum in args.spectra.

Parameters: