Utilities Module (`brutus.utils`)¶

brutus.utils provides the mathematical, photometric, and Monte Carlo sampling routines used across brutus.

Module Organization:

brutus.utils.photometry: Conversions between magnitudes, luptitudes, and linear flux densities in maggies (magnitude, inv_magnitude, luptitude, inv_luptitude, add_mag), plus Gaussian and outlier photometric log-likelihoods (phot_loglike, chisquare_outlier_loglike, uniform_outlier_loglike).
brutus.utils.math: Analytical \(3 \times 3\) matrix inversion and eigenvalue regularization (inverse3, isPSD), coordinate transformations (galactic_to_galactocentric_cyl), truncated normal and \(\chi^2\) densities (truncnorm_pdf, truncnorm_logpdf, chisquare_logpdf), and numerical bounds (LOG_ZERO, AVLIM, RVLIM).
brutus.utils.sampling: Weighted midpoint-CDF quantiles (quantile), antithetic multivariate normal sampling (sample_multivariate_normal), and \((s, A_V, R_V)\) posterior resampling (draw_sar).

Typical Usage:

Photometry utilities:

from brutus.utils.photometry import magnitude, add_mag, phot_loglike
import numpy as np

# Convert linear flux densities (maggies) to magnitudes
flux = np.array([1e-10, 5e-11, 2e-11])
flux_err = np.array([1e-12, 5e-13, 3e-13])
mags, mag_errs = magnitude(flux, flux_err)

# Combine magnitudes of two unresolved binary components
mag1, mag2 = 5.0, 6.0
mag_combined = add_mag(mag1, mag2)

# Evaluate photometric log-likelihood for data of shape (Nobj, Nfilt)
# against a model grid of shape (Nmod, Nfilt), returning (Nobj, Nmod)
obs_flux = np.array([[1.2e-10, 4.8e-11, 2.1e-11]])
obs_err = np.array([[1e-12, 5e-13, 3e-13]])
model_flux = np.array([[1.15e-10, 5.0e-11, 2.0e-11]])
lnl = phot_loglike(obs_flux, obs_err, model_flux, dim_prior=True)

Mathematical utilities:

from brutus.utils.math import inverse3, isPSD, chisquare_logpdf
import numpy as np

cov_matrix = np.array([[1.0, 0.1, 0.0],
                       [0.1, 2.0, 0.2],
                       [0.0, 0.2, 1.5]])
is_valid = isPSD(cov_matrix)
inv_cov = inverse3(cov_matrix, regularize=True)

chi2_values = np.array([1.0, 5.0, 10.0])
log_prob = chisquare_logpdf(chi2_values, dof=3)

Sampling utilities:

from brutus.utils.sampling import sample_multivariate_normal, quantile
import numpy as np

rng = np.random.default_rng(42)
mean = np.array([0.0, 0.0])
cov = np.array([[1.0, 0.5], [0.5, 2.0]])
samples = sample_multivariate_normal(mean, cov, size=1000, rstate=rng)

data = rng.standard_normal(1000)
weights = rng.random(1000)
q16, q50, q84 = quantile(data, [0.16, 0.5, 0.84], weights=weights)

See Also:

Model grids - How utilities are used in fitting algorithm
Understanding results - Sampling and quantile calculations

Photometry Functions¶

brutus.utils.magnitude(phot, err, zeropoints=1.0)[source]¶

Convert linear flux densities and uncertainties to magnitudes.

Parameters:

phot (~numpy.ndarray of shape (Nobs, Nfilt)) – Observed photometric flux densities (in maggies when zeropoints=1.0).
err (~numpy.ndarray of shape (Nobs, Nfilt)) – Observed photometric flux density uncertainties in the same units as phot.
zeropoints (float or ~numpy.ndarray of shape (Nfilt,), optional) – Flux density zero-points. Default is 1.0.

Returns:

mag (~numpy.ndarray of shape (Nobs, Nfilt)) – Magnitudes (mag) in the survey’s native system.
mag_err (~numpy.ndarray of shape (Nobs, Nfilt)) – First-order magnitude uncertainties (mag).

brutus.utils.inv_magnitude(mag, err, zeropoints=1.0)[source]¶

Convert magnitudes and uncertainties to linear flux densities.

Parameters:

mag (~numpy.ndarray of shape (Nobs, Nfilt)) – Magnitudes (mag) in the survey’s native system.
err (~numpy.ndarray of shape (Nobs, Nfilt)) – Magnitude uncertainties (mag).
zeropoints (float or ~numpy.ndarray of shape (Nfilt,), optional) – Flux density zero-points. Default is 1.0.

Returns:

phot (~numpy.ndarray of shape (Nobs, Nfilt)) – Photometric flux densities (in maggies when zeropoints=1.0).
phot_err (~numpy.ndarray of shape (Nobs, Nfilt)) – First-order flux density uncertainties in the same units as phot.

brutus.utils.luptitude(phot, err, skynoise=1.0, zeropoints=1.0)[source]¶

Convert photometry to asinh magnitudes (“Luptitudes”, Lupton et al. 1999).

Parameters:

phot (~numpy.ndarray with shape (Nobs, Nfilt)) – Observed photometric flux densities.
err (~numpy.ndarray with shape (Nobs, Nfilt)) – Observed photometric flux density errors.
skynoise (float or ~numpy.ndarray with shape (Nfilt,)) – Background sky noise. Used as a “softening parameter”. Default is 1..
zeropoints (float or ~numpy.ndarray with shape (Nfilt,)) – Flux density zero-points. Used as a “location parameter”. Default is 1..

Returns:

lupt (~numpy.ndarray with shape (Nobs, Nfilt)) – Luptitudes corresponding to input phot.
lupt_err (~numpy.ndarray with shape (Nobs, Nfilt)) – Luptitudes errors corresponding to input err.

brutus.utils.inv_luptitude(lupt, err, skynoise=1.0, zeropoints=1.0)[source]¶

Convert asinh magnitudes (Luptitudes) to photometry.

See Lupton et al. (1999) for more details.

Parameters:

lupt (~numpy.ndarray with shape (Nobs, Nfilt)) – Luptitudes.
err (~numpy.ndarray with shape (Nobs, Nfilt)) – Luptitude errors.
skynoise (float or ~numpy.ndarray with shape (Nfilt,)) – Background sky noise. Used as a “softening parameter”. Default is 1..
zeropoints (float or ~numpy.ndarray with shape (Nfilt,)) – Flux density zero-points. Used as a “location parameter”. Default is 1..

Returns:

phot (~numpy.ndarray with shape (Nobs, Nfilt)) – Photometric flux densities corresponding to input lupt.
phot_err (~numpy.ndarray with shape (Nobs, Nfilt)) – Photometric errors corresponding to input err.

brutus.utils.add_mag(mag1, mag2)[source]¶

Add magnitudes.

Parameters:

mag1 (float or ~numpy.ndarray) – First set of magnitudes.
mag2 (float or ~numpy.ndarray) – Second set of magnitudes.

Returns:

mag_combined – Combined magnitudes corresponding to the combined flux from mag1 and mag2.

Return type:

float or ~numpy.ndarray

brutus.utils.phot_loglike(flux, err, mfluxes, mask=None, dim_prior=False, dof_reduction=0, extra_chi2=None, extra_dims=None)[source]¶

Compute the log-likelihood between observed and model fluxes.

Parameters:

flux (~numpy.ndarray with shape (Nobj, Nfilt)) – Observed flux values.
err (~numpy.ndarray with shape (Nobj, Nfilt)) – Associated flux errors.
mfluxes (~numpy.ndarray with shape (Nobj, Nmod, Nfilt) or (Nmod, Nfilt)) – Model fluxes (for each model). A 2-D array means every object is compared against the same shared model grid; this triggers a much faster matrix-multiplication evaluation of the chi-square that avoids materializing any (Nobj, Nmod, Nfilt) intermediates.
mask (~numpy.ndarray with shape (Nobj, Nfilt), optional) – Binary mask indicating whether each observed band can be used (1) or should be skipped (0). Bands with non-finite flux/error or non-positive error are always excluded (they contribute neither to the chi-square nor to the effective dimensionality).
dim_prior (bool, optional) –
Whether to apply a dimensionality prior from the chi-squared distribution. Default is False.

Warning

When dim_prior=True, perfect model matches (chi2≈0) can cause problematic behavior: +inf likelihood for DOF=1, LOG_ZERO likelihood for DOF≥3. Ensure test data has small but non-zero residuals.
dof_reduction (int, optional) – Number of degrees of freedom to subtract from the effective DOF when using dim_prior=True. This accounts for parameters being fitted simultaneously (e.g., scale factors, extinction). Default is 0.
extra_chi2 (~numpy.ndarray with shape (Nobj,) or (Nobj, Nmod), optional) – Additional chi-square contributions (e.g., from a parallax measurement) added to the photometric chi-square before any dimensionality prior is applied. This keeps auxiliary Gaussian constraints in the same chi-square “measure” as the photometry when dim_prior=True (mirroring how BruteForce counts parallax as an extra band). With dim_prior=False the Gaussian normalization of the extra term is not added in full (see extra_dims); callers needing a proper density should add it themselves.
extra_dims (~numpy.ndarray with shape (Nobj,), optional) – Number of extra data points behind extra_chi2, added to the degrees of freedom of the dimensionality prior. Use it with dim_prior=True. With dim_prior=False it adds only the ln(2 pi) part of the normalization of the extra term, not the ln(var) part.

Returns:

lnl – Log-likelihood values. With dim_prior=True, an object with no degrees of freedom (Ndim - dof_reduction <= 0) gets 0 for every model, so it adds nothing to a sum over objects.

Return type:

~numpy.ndarray with shape (Nobj, Nmod)

Constants¶

brutus.utils.math.LOG_ZERO = -1e+300¶: Proxy for log(0) in log-probabilities. It stays finite under arithmetic, where -inf would give NaN (for example -inf - -inf).

brutus.utils.math.AVLIM = (0.0, 20.0)¶: Default fit limits on A(V), in mag.

brutus.utils.math.RVLIM = (1.0, 8.0)¶: Default fit limits on R(V).

Mathematical Functions¶

brutus.utils.galactic_to_galactocentric_cyl(dists, ell, b, R_solar=8.2, Z_solar=0.025)[source]¶

Convert Galactic coordinates to Galactocentric cylindrical coordinates.

Converts heliocentric Galactic coordinates (l, b, distance) to Galactocentric cylindrical coordinates (R, Z) using a simple rotation and translation. This is a fast NumPy-based replacement for astropy SkyCoord coordinate transformations.

Parameters:

dists (array_like) – Heliocentric distances in kpc.
ell (float or array_like) – Galactic longitude in degrees.
b (float or array_like) – Galactic latitude in degrees.
R_solar (float, optional) – Solar Galactocentric radius in kpc. Default is 8.2.
Z_solar (float, optional) – Solar height above the Galactic midplane in kpc. Default is 0.025.

Returns:

R (ndarray) – Galactocentric cylindrical radius in kpc.
Z (ndarray) – Height above/below the Galactic midplane in kpc.

Notes

The conversion assumes:

The Sun is located at (x, y, z) = (R_solar, 0, Z_solar) in Galactocentric Cartesian coordinates.
Galactic longitude l=0 points toward the Galactic center.
The Galactic midplane is at Z=0.

The Cartesian positions relative to the Sun are:

\[ \begin{align}\begin{aligned}x = d \cos(b) \cos(l)\\y = d \cos(b) \sin(l)\end{aligned}\end{align} \]

And the Galactocentric cylindrical coordinates are:

\[ \begin{align}\begin{aligned}R = \sqrt{(x - R_\odot)^2 + y^2}\\Z = d \sin(b) + Z_\odot\end{aligned}\end{align} \]

Sampling Utilities¶

brutus.utils.quantile(x, q, weights=None)[source]¶

Compute weighted or unweighted quantiles using the midpoint-CDF convention.

Evaluates the empirical cumulative distribution function at the sample midpoints (cumsum(w) - 0.5 * w) / sum(w) and interpolates linearly at q. Omitting weights is equivalent to passing uniform weights.

Parameters:

x (~numpy.ndarray with shape (nsamps,)) – Input samples.
q (~numpy.ndarray with shape (nquantiles,) or float) – The list of quantiles to compute from [0., 1.].
weights (~numpy.ndarray with shape (nsamps,), optional) – The associated weight from each sample. If None, all samples are weighted equally. Samples with zero weight are dropped before the CDF is built, so they never move a quantile.

Returns:

quantiles – The (weighted) sample quantiles computed at q. A scalar q gives shape (1,).

Return type:

~numpy.ndarray with shape (nquantiles,)

Raises:

ValueError – If quantiles are outside [0, 1], if dimensions don’t match, or if a weight is negative or not finite, or the weights sum to 0.

Examples

>>> import numpy as np
>>> x = np.array([1, 2, 3, 4, 5])
>>> q = np.array([0.25, 0.5, 0.75])
>>> quantile(x, q)
array([1.75, 3.  , 4.25])

>>> weights = np.array([1, 1, 1, 1, 10])  # Last sample heavily weighted
>>> quantile(x, q, weights=weights)
array([4.        , 4.63636364, 5.        ])

brutus.utils.sample_multivariate_normal(mean, cov, size=1, eps=1e-30, rstate=None, antithetic=False)[source]¶

Draw samples from many multivariate normal distributions.

Returns samples from an arbitrary number of multivariate distributions. The multivariate distributions must all have the same dimension. This function is optimized for drawing from many distributions simultaneously using Cholesky decomposition.

Parameters:

mean (~numpy.ndarray of shape (Ndist, dim) or (dim,)) – Means of the various multivariate distributions, where Ndist is the number of desired distributions and dim is the dimension of the distributions.
cov (~numpy.ndarray of shape (Ndist, dim, dim) or (dim, dim)) – Covariances of the various multivariate distributions, where Ndist is the number of desired distributions and dim is the dimension of the distributions.
size (int, optional) – Number of samples to draw from each distribution. Default is 1.
eps (float, optional) – Small factor added to covariances prior to Cholesky decomposition. Helps ensure numerical stability and should have no effect on the outcome. Default is 1e-30.
rstate (~numpy.random.RandomState, optional) – ~numpy.random.RandomState instance. If None, uses default numpy random state.
antithetic (bool, optional) – If True (only effective for the 3D fast path), draw the underlying standard normals in antithetic pairs (z, -z) along the sample axis. This leaves every Monte-Carlo estimate unbiased but reduces its variance (and halves the number of Gaussian draws generated). Default is False.

Returns:

samples – Sampled values. For a single distribution, returns (dim, size). For multiple distributions, returns (dim, size, Ndist).

Return type:

~numpy.ndarray of shape (dim, size, Ndist) or (dim, size)

Raises:

ValueError – If mean is not 1-D or 2-D, or if cov does not have shape (Ndist, dim, dim) for a 2-D mean.

Notes

Provided covariances must be positive semi-definite. Use the isPSD function from brutus.utils.math to check individual matrices if unsure.

For a single distribution, this function simply calls numpy’s multivariate_normal. For multiple distributions, it uses Cholesky decomposition for efficiency.

Examples

>>> import numpy as np
>>> # Single distribution
>>> mean = np.array([0.0, 1.0])
>>> cov = np.array([[1.0, 0.5], [0.5, 1.0]])
>>> samples = sample_multivariate_normal(mean, cov, size=100)
>>> samples.shape
(2, 100)

>>> # Multiple distributions
>>> means = np.array([[0.0, 1.0], [2.0, 3.0]])  # 2 distributions, 2D each
>>> covs = np.array([[[1.0, 0.0], [0.0, 1.0]], [[2.0, 0.5], [0.5, 2.0]]])
>>> samples = sample_multivariate_normal(means, covs, size=50)
>>> samples.shape
(2, 50, 2)

brutus.utils.draw_sar(scales, avs, rvs, covs_sar, ndraws=500, avlim=(0.0, 20.0), rvlim=(1.0, 8.0), rstate=None, max_attempts=10000)[source]¶

Generate random draws from the joint scale-A_V-R_V posterior for a given object.

This function generates Monte Carlo samples from the joint posterior of scale factors, reddening (A_V), and reddening curve shape (R_V) for stellar fitting applications.

Parameters:

scales (~numpy.ndarray of shape (Nsamps)) – An array of scale factors s derived between the models and the data.
avs (~numpy.ndarray of shape (Nsamps)) – An array of reddenings A(V) derived for the models.
rvs (~numpy.ndarray of shape (Nsamps)) – An array of reddening shapes R(V) derived for the models.
covs_sar (~numpy.ndarray of shape (Nsamps, 3, 3)) – An array of covariance matrices corresponding to (scales, avs, rvs).
ndraws (int, optional) – The number of desired random draws. Default is 500.
avlim (2-tuple, optional) – The A_V limits used to truncate results. Default is AVLIM, (0, 20), the BruteForce.fit default.
rvlim (2-tuple, optional) – The R_V limits used to truncate results. Default is RVLIM, (1, 8).
rstate (~numpy.random.RandomState, ~numpy.random.Generator, or) – module, optional Random state used for the draws (anything exposing a normal(loc, scale, size) method, including the numpy.random module itself). If None, uses the default numpy random state.
max_attempts (int, optional) – Maximum number of rejection-sampling passes per posterior sample before the remaining slots are padded with the mean values (with a warning). A sample with no accepted draw after 100 passes is padded at once. Default is 10000.

Returns:

sdraws (~numpy.ndarray of shape (Nsamps, Ndraws)) – Scale-factor samples.
adraws (~numpy.ndarray of shape (Nsamps, Ndraws)) – Reddening (A_V) samples.
rdraws (~numpy.ndarray of shape (Nsamps, Ndraws)) – Reddening shape (R_V) samples.

Notes

The function samples from multivariate normal distributions defined by the means (scales, avs, rvs) and covariances (covs_sar), then applies rejection sampling to keep the samples within the specified limits for A_V and R_V. A sample that cannot fill ndraws this way is padded with its mean, which can lie outside the limits.

All Nsamps distributions are drawn in a single batched sample_multivariate_normal call per rejection pass (rather than one numpy.random.multivariate_normal call per distribution), which avoids per-call dispatch and per-matrix SVD overhead.

Examples

>>> import numpy as np
>>> scales = np.array([1.0, 1.1])
>>> avs = np.array([0.1, 0.2])
>>> rvs = np.array([3.1, 3.3])
>>> covs_sar = np.array([[[0.01, 0, 0], [0, 0.01, 0], [0, 0, 0.1]],
...                      [[0.01, 0, 0], [0, 0.01, 0], [0, 0, 0.1]]])
>>> sdraws, adraws, rdraws = draw_sar(scales, avs, rvs, covs_sar, ndraws=100)
>>> sdraws.shape
(2, 100)

Utilities Module (brutus.utils)¶

Photometry Functions¶

Constants¶

Mathematical Functions¶

Sampling Utilities¶

Utilities Module (`brutus.utils`)¶