Utilities Module (brutus.utils)

Utilities Module (brutus.utils)#

The utils module provides low-level mathematical, photometric, and sampling utilities that support brutus’s core functionality. These functions handle magnitude/flux conversions, likelihood calculations, statistical distributions, and posterior sampling.

Module Organization:

  • ``brutus.utils.photometry``: Magnitude/flux conversions and photometric likelihoods

  • ``brutus.utils.math``: Mathematical utilities (matrix operations, distributions, special functions)

  • ``brutus.utils.sampling``: Posterior sampling and Monte Carlo methods

Typical Usage:

Most users won’t call these utilities directly (they’re used internally by BruteForce and other high-level classes). However, they’re useful for:

  • Custom analysis scripts

  • Implementing new fitting methods

  • Debugging and validation

Photometry utilities:

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

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

# Add magnitudes (combine flux from binaries)
mag1, mag2 = 5.0, 6.0
mag_combined = add_mag(mag1, mag2)  # Brighter than either component

# Compute photometric log-likelihood
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, chisquare_logpdf

# Fast 3x3 matrix inversion (used in covariance calculations)
cov_matrix = np.array([[1.0, 0.1, 0.0],
                       [0.1, 2.0, 0.2],
                       [0.0, 0.2, 1.5]])
inv_cov = inverse3(cov_matrix)

# Chi-square log-PDF (used in outlier models)
chi2_values = np.array([1.0, 5.0, 10.0])
dof = 3
log_prob = chisquare_logpdf(chi2_values, dof)

Sampling utilities:

from brutus.utils.sampling import sample_multivariate_normal, quantile

# Draw samples from multivariate normal
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=10000)

# Compute weighted quantiles
data = np.random.randn(10000)
weights = np.random.rand(10000)
q16, q50, q84 = quantile(data, [0.16, 0.5, 0.84], weights=weights)

Performance Notes:

  • Many functions are JIT-compiled with numba for speed

  • Matrix operations optimized for small dimensions (common in stellar fitting)

  • Photometric likelihoods vectorized for batch processing

See Also:

Photometry Functions#

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

Convert photometry to AB magnitudes.

Parameters:

  • phot (~numpy.ndarray with shape (Nobs, Nfilt)) – Observed photometric flux densities.

  • err (~numpy.ndarray with shape (Nobs, Nfilt)) – Observed photometric flux density errors.

  • zeropoints (float or ~numpy.ndarray with shape (Nfilt,)) – Flux density zero-points. Used as a “location parameter”. Default is 1..

Returns:

  • mag (~numpy.ndarray with shape (Nobs, Nfilt)) – Magnitudes corresponding to input phot.

  • mag_err (~numpy.ndarray with shape (Nobs, Nfilt)) – Magnitudes errors corresponding to input err.

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

Convert AB magnitudes to photometry.

Parameters:

  • mag (~numpy.ndarray with shape (Nobs, Nfilt)) – Magnitudes.

  • err (~numpy.ndarray with shape (Nobs, Nfilt)) – Magnitude errors.

  • 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 mag.

  • phot_err (~numpy.ndarray with shape (Nobs, Nfilt)) – Photometric errors corresponding to input err.

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

Convert photometry to asinh magnitudes (i.e. “Luptitudes”). See Lupton et.

al. (1999) for more details.

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)[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)) – Model fluxes (for each model).

  • mask (~numpy.ndarray with shape (Nobj, Nfilt), optional) – Binary mask indicating whether each observed band can be used (1) or should be skipped (0).

  • 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, -inf 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.

Returns:

lnl – Log-likelihood values.

Return type:

~numpy.ndarray with shape (Nobj, Nmod)

See also

chisquare_outlier_loglike

Outlier model for dim_prior=True

uniform_outlier_loglike

Outlier model for dim_prior=False

brutus.analysis.individual.BruteForce.loglike_grid

Uses this function

Notes

The log-likelihood without dimensionality prior is:

\[\begin{split}\\ln L = -\\frac{1}{2} \\left[ \\chi^2 + N \\ln(2\\pi) + \\ln|\\Sigma| \\right]\end{split}\]

where \(\\chi^2 = \\sum_i (f_i - m_i)^2/\\sigma_i^2\) and \(|\\Sigma|\) is the determinant of the covariance matrix.

With dimensionality prior (dim_prior=True), this becomes the log-PDF of a chi-square distribution, which weights models by goodness-of-fit relative to the number of degrees of freedom.

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} \]

See also

brutus.priors.galactic.logp_galactic_structure

Uses this for coordinate conversion instead of astropy SkyCoord.

Examples

>>> import numpy as np
>>> from brutus.utils.math import galactic_to_galactocentric_cyl
>>> R, Z = galactic_to_galactocentric_cyl(
...     dists=np.array([1.0, 5.0]), ell=90.0, b=0.0
... )
brutus.utils.inverse3(A, regularize=False, min_eigenval_threshold=1e-12)[source]#

Compute the inverse of a series of 3x3 matrices.

When regularize=True, uses diagonal preconditioning: the precision matrix is normalized to a correlation-like matrix (1s on diagonal) before inversion, then transformed back. This reduces condition numbers from O(max_diag/min_diag) to O(1/(1-rho_max)), making the inversion numerically stable even when parameters have very different scales (e.g. scale ~ 10^5 vs R(V) precision ~ 30).

Parameters:

  • A (~numpy.ndarray of shape (…, 3, 3)) – Array of 3x3 matrices.

  • regularize (bool, optional) – Whether to apply diagonal preconditioning and regularization to ensure positive semi-definiteness. Default: False.

  • min_eigenval_threshold (float, optional) – Minimum acceptable eigenvalue for OUTPUT matrices in the normalized space. Default: 1e-12.

Returns:

A_inv – Inverse matrices, guaranteed to be positive semi-definite if regularize=True.

Return type:

~numpy.ndarray of shape (…, 3, 3)

brutus.utils.isPSD(A)[source]#

Check if A is a positive semidefinite matrix.

A matrix is positive semidefinite if all its eigenvalues are non-negative.

Parameters:

A (~numpy.ndarray of shape (N, N)) – Square matrix to test.

Returns:

is_psd – True if the matrix is positive semidefinite, False otherwise.

Return type:

bool

brutus.utils.chisquare_logpdf(x, df, loc=0, scale=1)[source]#

Compute log-PDF of a chi-square distribution.

chisquare_logpdf(x, df, loc, scale) is equal to chisquare_logpdf(y, df) - ln(scale), where y = (x-loc)/scale. NOTE: This function replicates ~scipy.stats.chi2.logpdf.

Parameters:

  • x (~numpy.ndarray of shape (N) or float) – Input values.

  • df (float) – Degrees of freedom. Uses math.lgamma internally, so there is no overflow risk for large df (unlike math.gamma which overflows for df > ~340).

  • loc (float, optional) – Offset of distribution. Default is 0.

  • scale (float, optional) – Scaling of distribution. Default is 1.

Returns:

ans – The natural log of the PDF values.

Return type:

~numpy.ndarray of shape (N) or float

brutus.utils.truncnorm_pdf(x, a, b, loc=0.0, scale=1.0)[source]#

Compute PDF of a truncated normal distribution.

The parent normal distribution has a mean of loc and standard deviation of scale. The distribution is cut off at a and b. NOTE: This function replicates ~scipy.stats.truncnorm.pdf.

Parameters:

  • x (~numpy.ndarray of shape (N) or float) – Input values.

  • a (float) – Lower bound in standardized units. The actual lower cutoff in data space is scale * a + loc.

  • b (float) – Upper bound in standardized units. The actual upper cutoff in data space is scale * b + loc.

  • loc (float, optional) – Mean of normal distribution. Default is 0.0.

  • scale (float, optional) – Standard deviation of normal distribution. Default is 1.0.

Returns:

ans – The PDF values.

Return type:

~numpy.ndarray of shape (N) or float

brutus.utils.truncnorm_logpdf(x, a, b, loc=0.0, scale=1.0)[source]#

Compute log-PDF of a truncated normal distribution.

The parent normal distribution has a mean of loc and standard deviation of scale. The distribution is cut off at a and b. NOTE: This function replicates ~scipy.stats.truncnorm.logpdf.

Parameters:

  • x (~numpy.ndarray of shape (N) or float) – Input values.

  • a (float) – Lower bound in standardized units. The actual lower cutoff in data space is scale * a + loc.

  • b (float) – Upper bound in standardized units. The actual upper cutoff in data space is scale * b + loc.

  • loc (float, optional) – Mean of normal distribution. Default is 0.0.

  • scale (float, optional) – Standard deviation of normal distribution. Default is 1.0.

Returns:

ans – The natural log PDF values.

Return type:

~numpy.ndarray of shape (N) or float

Sampling Utilities#

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

Compute (weighted) quantiles from an input set of samples.

This function computes quantiles from a set of samples, optionally with weights. For unweighted samples, it uses numpy.percentile. For weighted samples, it computes the cumulative distribution function and interpolates to find the desired quantiles.

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.

Returns:

quantiles – The (weighted) sample quantiles computed at q.

Return type:

~numpy.ndarray with shape (nquantiles,) or float

Raises:

ValueError – If quantiles are outside [0, 1] or if dimensions don’t match.

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([2., 3., 4.])

>>> 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)[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.

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)

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, 1])
>>> cov = np.array([[1, 0.5], [0.5, 1]])
>>> samples = sample_multivariate_normal(mean, cov, size=100)
>>> samples.shape
(2, 100)

>>> # Multiple distributions
>>> means = np.array([[0, 1], [2, 3]])  # 2 distributions, 2D each
>>> covs = np.array([[[1, 0], [0, 1]], [[2, 0.5], [0.5, 2]]])
>>> 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, 6.0), rvlim=(1.0, 8.0), rstate=None)[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 (0., 6.).

  • rvlim (2-tuple, optional) – The R_V limits used to truncate results. Default is (1., 8.).

  • rstate (~numpy.random.RandomState, optional) – ~numpy.random.RandomState instance. If None, uses default numpy random state (or intel-specific version if available).

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 ensure all samples fall within the specified limits for A_V and R_V.

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)

Submodules#

For advanced users who need access to internal implementations:

Mathematical utility functions for brutus.

This module contains mathematical utility functions including matrix operations, statistical distributions, and numerical utilities. Many functions are JIT-compiled with numba for performance.

Functions#

galactic_to_galactocentric_cylCoordinate transform

Convert Galactic (l, b, distance) to galactocentric cylindrical (R, phi, Z)

inverse3Matrix inversion

Fast 3x3 matrix inversion with optional diagonal preconditioning

isPSDMatrix check

Check if matrix is positive semi-definite

chisquare_logpdfChi-square log-PDF

Log-probability density for chi-square distribution

truncnorm_pdfTruncated normal PDF

Probability density for truncated normal

truncnorm_logpdfTruncated normal log-PDF

Log-probability density for truncated normal

See also

brutus.utils.sampling

Sampling utilities

brutus.priors.galactic

Uses truncated normal distributions

Notes

The numba-compiled functions provide significant speedups for tight loops in Bayesian inference. The inverse3 function is specifically optimized for the (scale, A_V, R_V) covariance matrices used throughout brutus.

Matrix regularization in inverse3 prevents numerical issues when matrices are near-singular by adding a small value to the diagonal when eigenvalues are too small.

Examples

>>> import numpy as np
>>> from brutus.utils.math import inverse3, isPSD
>>>
>>> # Create a 3x3 covariance matrix
>>> cov = np.array([[1.0, 0.1, 0.05],
...                 [0.1, 0.5, 0.02],
...                 [0.05, 0.02, 0.3]])
>>>
>>> # Check if positive semi-definite
>>> is_valid = isPSD(cov)
>>>
>>> # Invert with diagonal preconditioning / regularization
>>> icov = inverse3(cov, regularize=True)
brutus.utils.math.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} \]

See also

brutus.priors.galactic.logp_galactic_structure

Uses this for coordinate conversion instead of astropy SkyCoord.

Examples

>>> import numpy as np
>>> from brutus.utils.math import galactic_to_galactocentric_cyl
>>> R, Z = galactic_to_galactocentric_cyl(
...     dists=np.array([1.0, 5.0]), ell=90.0, b=0.0
... )
brutus.utils.math.inverse3(A, regularize=False, min_eigenval_threshold=1e-12)[source]

Compute the inverse of a series of 3x3 matrices.

When regularize=True, uses diagonal preconditioning: the precision matrix is normalized to a correlation-like matrix (1s on diagonal) before inversion, then transformed back. This reduces condition numbers from O(max_diag/min_diag) to O(1/(1-rho_max)), making the inversion numerically stable even when parameters have very different scales (e.g. scale ~ 10^5 vs R(V) precision ~ 30).

Parameters:

  • A (~numpy.ndarray of shape (…, 3, 3)) – Array of 3x3 matrices.

  • regularize (bool, optional) – Whether to apply diagonal preconditioning and regularization to ensure positive semi-definiteness. Default: False.

  • min_eigenval_threshold (float, optional) – Minimum acceptable eigenvalue for OUTPUT matrices in the normalized space. Default: 1e-12.

Returns:

A_inv – Inverse matrices, guaranteed to be positive semi-definite if regularize=True.

Return type:

~numpy.ndarray of shape (…, 3, 3)

brutus.utils.math.isPSD(A)[source]

Check if A is a positive semidefinite matrix.

A matrix is positive semidefinite if all its eigenvalues are non-negative.

Parameters:

A (~numpy.ndarray of shape (N, N)) – Square matrix to test.

Returns:

is_psd – True if the matrix is positive semidefinite, False otherwise.

Return type:

bool

brutus.utils.math.chisquare_logpdf(x, df, loc=0, scale=1)[source]

Compute log-PDF of a chi-square distribution.

chisquare_logpdf(x, df, loc, scale) is equal to chisquare_logpdf(y, df) - ln(scale), where y = (x-loc)/scale. NOTE: This function replicates ~scipy.stats.chi2.logpdf.

Parameters:

  • x (~numpy.ndarray of shape (N) or float) – Input values.

  • df (float) – Degrees of freedom. Uses math.lgamma internally, so there is no overflow risk for large df (unlike math.gamma which overflows for df > ~340).

  • loc (float, optional) – Offset of distribution. Default is 0.

  • scale (float, optional) – Scaling of distribution. Default is 1.

Returns:

ans – The natural log of the PDF values.

Return type:

~numpy.ndarray of shape (N) or float

brutus.utils.math.truncnorm_pdf(x, a, b, loc=0.0, scale=1.0)[source]

Compute PDF of a truncated normal distribution.

The parent normal distribution has a mean of loc and standard deviation of scale. The distribution is cut off at a and b. NOTE: This function replicates ~scipy.stats.truncnorm.pdf.

Parameters:

  • x (~numpy.ndarray of shape (N) or float) – Input values.

  • a (float) – Lower bound in standardized units. The actual lower cutoff in data space is scale * a + loc.

  • b (float) – Upper bound in standardized units. The actual upper cutoff in data space is scale * b + loc.

  • loc (float, optional) – Mean of normal distribution. Default is 0.0.

  • scale (float, optional) – Standard deviation of normal distribution. Default is 1.0.

Returns:

ans – The PDF values.

Return type:

~numpy.ndarray of shape (N) or float

brutus.utils.math.truncnorm_logpdf(x, a, b, loc=0.0, scale=1.0)[source]

Compute log-PDF of a truncated normal distribution.

The parent normal distribution has a mean of loc and standard deviation of scale. The distribution is cut off at a and b. NOTE: This function replicates ~scipy.stats.truncnorm.logpdf.

Parameters:

  • x (~numpy.ndarray of shape (N) or float) – Input values.

  • a (float) – Lower bound in standardized units. The actual lower cutoff in data space is scale * a + loc.

  • b (float) – Upper bound in standardized units. The actual upper cutoff in data space is scale * b + loc.

  • loc (float, optional) – Mean of normal distribution. Default is 0.0.

  • scale (float, optional) – Standard deviation of normal distribution. Default is 1.0.

Returns:

ans – The natural log PDF values.

Return type:

~numpy.ndarray of shape (N) or float

Photometric utility functions for brutus.

This module contains functions for converting between magnitudes and fluxes, handling asinh magnitudes (“luptitudes”), and computing photometric likelihoods. These utilities are fundamental for all photometric analysis in brutus.

Functions#

magnitudeConvert flux to magnitudes

AB magnitude system conversion

inv_magnitudeConvert magnitudes to flux

Inverse of magnitude conversion

luptitudeConvert flux to asinh magnitudes

Asinh magnitude system (Lupton et al. 1999)

inv_luptitudeConvert asinh magnitudes to flux

Inverse of luptitude conversion

add_magCombine magnitudes

Add fluxes in magnitude space

phot_loglikePhotometric log-likelihood

Core likelihood function for stellar fitting

chisquare_outlier_loglikeChi-square outlier model

Outlier likelihood for mixture models

uniform_outlier_loglikeUniform outlier model

Alternative outlier likelihood

See also

brutus.analysis.individual.BruteForce

Uses phot_loglike for fitting

brutus.analysis.populations

Uses outlier models for mixture fitting

brutus.core.sed_utils

SED generation functions

Notes

The photometric likelihood function phot_loglike is the core of brutus’s Bayesian inference machinery. It supports:

  • Multi-filter photometry with flexible masking

  • Optional dimensionality prior (chi-square distribution)

  • Degrees of freedom reduction for fitted parameters

The asinh magnitude system (luptitudes) provides better behavior than standard magnitudes for faint sources with high noise, preventing negative flux issues.

Examples

Basic magnitude conversion:

>>> import numpy as np
>>> from brutus.utils.photometry import magnitude, inv_magnitude
>>>
>>> # Convert flux to magnitude
>>> flux = np.array([[100, 200]])  # arbitrary flux units
>>> flux_err = np.array([[10, 20]])
>>> mags, mag_errs = magnitude(flux, flux_err)
>>>
>>> # Convert back
>>> flux_recovered, flux_err_recovered = inv_magnitude(mags, mag_errs)

Photometric likelihood:

>>> from brutus.utils.photometry import phot_loglike
>>>
>>> # Observed photometry
>>> obs_flux = np.array([[1.0, 2.0, 3.0]])  # 1 object, 3 filters
>>> obs_err = np.array([[0.1, 0.2, 0.3]])
>>>
>>> # Model predictions for 10 models
>>> model_flux = np.random.uniform(0.5, 4.0, (1, 10, 3))
>>>
>>> # Compute likelihoods
>>> lnl = phot_loglike(obs_flux, obs_err, model_flux)
>>> # lnl.shape is (1, 10) - likelihood for each model
brutus.utils.photometry.magnitude(phot, err, zeropoints=1.0)[source]

Convert photometry to AB magnitudes.

Parameters:

  • phot (~numpy.ndarray with shape (Nobs, Nfilt)) – Observed photometric flux densities.

  • err (~numpy.ndarray with shape (Nobs, Nfilt)) – Observed photometric flux density errors.

  • zeropoints (float or ~numpy.ndarray with shape (Nfilt,)) – Flux density zero-points. Used as a “location parameter”. Default is 1..

Returns:

  • mag (~numpy.ndarray with shape (Nobs, Nfilt)) – Magnitudes corresponding to input phot.

  • mag_err (~numpy.ndarray with shape (Nobs, Nfilt)) – Magnitudes errors corresponding to input err.

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

Convert AB magnitudes to photometry.

Parameters:

  • mag (~numpy.ndarray with shape (Nobs, Nfilt)) – Magnitudes.

  • err (~numpy.ndarray with shape (Nobs, Nfilt)) – Magnitude errors.

  • 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 mag.

  • phot_err (~numpy.ndarray with shape (Nobs, Nfilt)) – Photometric errors corresponding to input err.

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

Convert photometry to asinh magnitudes (i.e. “Luptitudes”). See Lupton et.

al. (1999) for more details.

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.photometry.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.photometry.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.photometry.phot_loglike(flux, err, mfluxes, mask=None, dim_prior=False, dof_reduction=0)[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)) – Model fluxes (for each model).

  • mask (~numpy.ndarray with shape (Nobj, Nfilt), optional) – Binary mask indicating whether each observed band can be used (1) or should be skipped (0).

  • 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, -inf 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.

Returns:

lnl – Log-likelihood values.

Return type:

~numpy.ndarray with shape (Nobj, Nmod)

See also

chisquare_outlier_loglike

Outlier model for dim_prior=True

uniform_outlier_loglike

Outlier model for dim_prior=False

brutus.analysis.individual.BruteForce.loglike_grid

Uses this function

Notes

The log-likelihood without dimensionality prior is:

\[\begin{split}\\ln L = -\\frac{1}{2} \\left[ \\chi^2 + N \\ln(2\\pi) + \\ln|\\Sigma| \\right]\end{split}\]

where \(\\chi^2 = \\sum_i (f_i - m_i)^2/\\sigma_i^2\) and \(|\\Sigma|\) is the determinant of the covariance matrix.

With dimensionality prior (dim_prior=True), this becomes the log-PDF of a chi-square distribution, which weights models by goodness-of-fit relative to the number of degrees of freedom.

brutus.utils.photometry.chisquare_outlier_loglike(flux, err, stellar_params=None, parallax=None, parallax_err=None, p_value_cut=1e-05)[source]

Compute chi-square based outlier model log-likelihood.

Uses a chi-square distribution with a p-value cut to model outlier probabilities. This is the default outlier model for dim_prior=True.

Parameters:

  • flux (array-like, shape (Nobj, Nfilt)) – Observed flux values

  • err (array-like, shape (Nobj, Nfilt)) – Flux errors

  • stellar_params (dict, optional) – Stellar parameters (masses, colors, etc.) - not used in current implementation but provided for future stellar-dependent outlier models

  • parallax (array-like, shape (Nobj,), optional) – Parallax measurements (mas)

  • parallax_err (array-like, shape (Nobj,), optional) – Parallax errors (mas)

  • p_value_cut (float, optional) – P-value threshold for outlier definition. Default 1e-5.

Returns:

lnl_outlier – Log-likelihood for outlier model for each object

Return type:

array-like, shape (Nobj,)

brutus.utils.photometry.uniform_outlier_loglike(flux, err, stellar_params=None, parallax=None, parallax_err=None, sigma_clip=3.0)[source]

Compute quasi-uniform outlier model log-likelihood.

Assumes uniform distribution within +/- sigma_clip * error bounds around the data. This is the default outlier model for dim_prior=False.

Parameters:

  • flux (array-like, shape (Nobj, Nfilt)) – Observed flux values

  • err (array-like, shape (Nobj, Nfilt)) – Flux errors

  • stellar_params (dict, optional) – Stellar parameters - not used in current implementation

  • parallax (array-like, shape (Nobj,), optional) – Parallax measurements (mas)

  • parallax_err (array-like, shape (Nobj,), optional) – Parallax errors (mas)

  • sigma_clip (float, optional) – Number of sigma for uniform bounds. Default 3.0.

Returns:

lnl_outlier – Log-likelihood for outlier model for each object

Return type:

array-like, shape (Nobj,)

Sampling utility functions for brutus.

This module contains functions for statistical sampling, quantile computation, and random number generation used in Bayesian inference workflows. These utilities are essential for posterior sampling and uncertainty quantification.

Functions#

quantileWeighted quantiles

Compute (weighted) quantiles from samples

draw_sarPosterior sampling

Draw from (scale, A_V, R_V) posterior

sample_multivariate_normalGaussian sampling

Sample from multivariate normal with bounds

See also

brutus.analysis.individual.BruteForce

Uses these for posterior sampling

brutus.utils.math

Mathematical utilities

Notes

The quantile function supports weighted samples, which is crucial for computing credible intervals from posterior samples with non-uniform weights (e.g., from importance sampling or nested sampling).

The draw_sar function is specifically designed for sampling the joint posterior of distance scale, extinction, and reddening curve shape from BruteForce fitting results.

Examples

Weighted quantile computation:

>>> import numpy as np
>>> from brutus.utils.sampling import quantile
>>>
>>> # Samples with different weights
>>> samples = np.array([1, 2, 3, 4, 5])
>>> weights = np.array([1, 1, 1, 1, 10])  # Last sample heavily weighted
>>>
>>> # Compute median and 68% credible interval
>>> q = np.array([0.16, 0.5, 0.84])
>>> intervals = quantile(samples, q, weights=weights)
>>> print(f"Median: {intervals[1]:.2f}")

Drawing from posterior:

>>> from brutus.utils.sampling import draw_sar
>>>
>>> # Posterior means and covariances from fitting
>>> # scales, avs, rvs, covs_sar = ... (from BruteForce)
>>>
>>> # Generate posterior samples
>>> # samples = draw_sar(scales, avs, rvs, covs_sar, ndraws=1000)
brutus.utils.sampling.quantile(x, q, weights=None)[source]

Compute (weighted) quantiles from an input set of samples.

This function computes quantiles from a set of samples, optionally with weights. For unweighted samples, it uses numpy.percentile. For weighted samples, it computes the cumulative distribution function and interpolates to find the desired quantiles.

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.

Returns:

quantiles – The (weighted) sample quantiles computed at q.

Return type:

~numpy.ndarray with shape (nquantiles,) or float

Raises:

ValueError – If quantiles are outside [0, 1] or if dimensions don’t match.

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([2., 3., 4.])

>>> weights = np.array([1, 1, 1, 1, 10])  # Last sample heavily weighted
>>> quantile(x, q, weights=weights)
array([4.        , 4.63636364, 5.        ])
brutus.utils.sampling.draw_sar(scales, avs, rvs, covs_sar, ndraws=500, avlim=(0.0, 6.0), rvlim=(1.0, 8.0), rstate=None)[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 (0., 6.).

  • rvlim (2-tuple, optional) – The R_V limits used to truncate results. Default is (1., 8.).

  • rstate (~numpy.random.RandomState, optional) – ~numpy.random.RandomState instance. If None, uses default numpy random state (or intel-specific version if available).

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 ensure all samples fall within the specified limits for A_V and R_V.

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)
brutus.utils.sampling.sample_multivariate_normal(mean, cov, size=1, eps=1e-30, rstate=None)[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.

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)

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, 1])
>>> cov = np.array([[1, 0.5], [0.5, 1]])
>>> samples = sample_multivariate_normal(mean, cov, size=100)
>>> samples.shape
(2, 100)

>>> # Multiple distributions
>>> means = np.array([[0, 1], [2, 3]])  # 2 distributions, 2D each
>>> covs = np.array([[[1, 0], [0, 1]], [[2, 0.5], [0.5, 2]]])
>>> samples = sample_multivariate_normal(means, covs, size=50)
>>> samples.shape
(2, 50, 2)