Analysis Module (`brutus.analysis`)¶

brutus.analysis implements the Bayesian inference engines for fitting individual stars, coeval star clusters, empirical photometric offsets, and 3D line-of-sight dust profiles.

Key Components:

Individual Star Fitting: BruteForce evaluates Bayesian posteriors over pre-computed stellar model grids.
Population Analysis: isochrone_population_loglike and its component functions evaluate cluster likelihoods with binary stars and field-star mixture models.
Photometric Offsets: photometric_offsets and PhotometricOffsetsConfig derive leave-one-out empirical zero-point corrections across survey bands.
Line-of-Sight Dust: los_clouds_priortransform and los_clouds_loglike_samples fit multi-cloud 3D extinction profiles to stellar distance-reddening posteriors.

Inference Workflow:

Grid Likelihood Evaluation: Compute photometric and parallax likelihoods across the stellar model grid.
Continuous Parameter Optimization: Solve for the maximum-likelihood \((s, A_V, R_V)\) and \(3 \times 3\) Fisher matrix per model via Gauss-Newton iteration.
Prior Integration: Integrate Galactic structure, 3D dust, and parallax priors via antithetic Monte Carlo importance sampling in \((\ln d, A_V, R_V)\).
Posterior Resampling: Draw equally weighted posterior samples of stellar models, distances, and extinctions.

Typical Usage Patterns:

For individual field stars with photometry and parallax:

from brutus.analysis import BruteForce
from brutus.core import StarGrid
from brutus.data import load_models
from brutus.utils import inv_magnitude
import numpy as np
import h5py

# Load pre-computed grid
models, labels, label_mask = load_models("grid_mist_v9.h5")
grid = StarGrid(models, labels)

# Initialize fitter
fitter = BruteForce(grid)

# Observed data (Nstars=1 in this example)
# Convert observed magnitudes and errors to linear flux densities (maggies)
mag = np.array([[16.5, 15.2, 14.8, 13.5, 13.1]])  # shape (1, 5)
mag_err = np.array([[0.02, 0.02, 0.02, 0.03, 0.03]])
flux, flux_err = inv_magnitude(mag, mag_err)  # maggies = 10**(-0.4 * mag)
mask = np.ones_like(flux, dtype=bool)
obj_ids = np.array([1])
coords = np.array([[209.1, -19.9]])  # Galactic (l, b) in degrees
parallax = np.array([2.5])  # mas
parallax_err = np.array([0.1])  # mas

# Fit and save results to HDF5
fitter.fit(
    data=flux, data_err=flux_err, data_mask=mask,
    data_labels=obj_ids, save_file="results.h5",
    parallax=parallax, parallax_err=parallax_err,
    data_coords=coords,
    Ndraws=250,
)

# Read posterior samples from HDF5
with h5py.File("results.h5", "r") as f:
    dist = f["samps_dist"][0]  # (Ndraws,) distances in kpc
    av = f["samps_red"][0]     # (Ndraws,) A_V values in mag
print(f"Distance: {np.median(dist) * 1000:.0f} pc")
print(f"Extinction: {np.median(av):.2f} mag")

For stellar clusters with MCMC:

from brutus.analysis import isochrone_population_loglike
from brutus.core import Isochrone, StellarPop
import emcee
import numpy as np

# Initialize population model
iso = Isochrone()
pop = StellarPop(isochrone=iso)

# Observed cluster data (N stars, M filters)
obs_flux = np.array([...])  # shape (N, M)
obs_err = np.array([...])   # shape (N, M)

# Define log-likelihood for MCMC
# theta = [feh, loga, av, rv, dist, field_frac]
# (field_frac is the field contamination fraction, 0 to 1;
# a 5-element theta raises ValueError)
def lnprob(theta):
    return isochrone_population_loglike(
        theta,
        stellarpop=pop,
        obs_phot=obs_flux,
        obs_err=obs_err,
        cluster_prob=0.9  # External membership prior
    )

# Run MCMC
ndim, nwalkers = 6, 32
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob)
sampler.run_mcmc(initial_pos, 5000, progress=True)

# Extract posterior samples
samples = sampler.get_chain(discard=1000, thin=10, flat=True)

Methodological Highlights:

Full-grid evaluation: Evaluates every model in the stellar grid, avoiding MCMC burn-in and multimodal trapping.
Two-stage Gauss-Newton optimization: Profiles distance out in magnitude space across all models before refining high-weight candidates in flux space.
Pointwise cluster mixture modeling: Evaluates the cluster-versus-field mixture probability at each \((m, q)\) grid point prior to marginalizing over mass and binary mass ratio.
Antithetic Monte Carlo integration: Integrates continuous spatial and dust priors in \((\ln d, A_V, R_V)\) using paired draws \((z, -z)\).

See Also:

Model grids - Understanding the brute-force fitting algorithm
Star clusters - Cluster population fitting methodology
Photometric offsets - Empirical calibration procedures
Understanding results - Interpreting posterior distributions

Individual Star Fitting¶

class brutus.analysis.BruteForce(star_grid, verbose=True)[source]¶

Bases: object

Estimate stellar parameters, distances, and extinctions over a StarGrid.

Evaluates the likelihood of every model in star_grid by optimizing (scale, A(V), R(V)) per grid point and integrates the posterior over distance and extinction uncertainty via Monte Carlo sampling.

Parameters:

star_grid (StarGrid) – Pre-loaded stellar model grid for SED generation.
verbose (bool, optional) – Whether to print initialization information. Default is True.

star_grid¶

The underlying stellar model grid.

Type:: StarGrid

models¶

Magnitude coefficients from the grid.

Type:: numpy.ndarray

models_labels¶

Labels for each model in the grid: the fields of star_grid.labels followed by those of star_grid.params.

Type:: structured numpy.ndarray

labels_mask¶

Mask indicating which labels are grid parameters (True) vs predictions (False).

Type:: dict

See also

brutus.core.StarGrid: Stellar model grid infrastructure
brutus.data.load_models: Load pre-computed grids
loglike_grid: Compute log-likelihoods
logpost_grid: Compute log-posteriors

Notes

The BruteForce fitter uses a two-stage approach:

Likelihood computation (loglike_grid): Optimizes distance and extinction for each grid point to find maximum likelihood
Posterior computation (logpost_grid): Integrates over distance and extinction uncertainty using Monte Carlo, applying priors for Galactic structure, dust maps, and astrometry

The fitter automatically handles: - Grid parameter vs. prediction distinction - Age weighting for proper sampling - Grid spacing corrections - Parallax constraints - Galactic structure priors - Dust map priors

Examples

Basic usage with a pre-loaded grid:

>>> from brutus.core import StarGrid
>>> from brutus.analysis.individual import BruteForce
>>> from brutus.data import load_models
>>>
>>> # Load grid
>>> models, labels, label_mask = load_models('grid_mist_v9.h5')
>>> grid = StarGrid(models, labels)
>>>
>>> # Initialize fitter
>>> fitter = BruteForce(grid)
>>>
>>> # Fit data (writes posterior draws and summary metrics to 'results.h5')
>>> results_file = fitter.fit(
...     data, data_err, data_mask,
...     data_labels, save_file='results.h5',
...     parallax=parallax_data, parallax_err=parallax_err,
...     data_coords=coordinates
... )

__init__(star_grid, verbose=True)[source]¶: Initialize BruteForce with a StarGrid instance.

property nmodels¶: Number of models in the grid.

property nfilters¶: Number of filters in the grid.

get_sed_grid(indices=None, av=None, rv=None, return_flux=False)[source]¶

Compute SEDs for multiple grid points simultaneously.

This is the grid-based batch computation method, distinct from StarGrid.get_seds() which handles single star synthesis.

Parameters:

indices (array-like, optional) – Grid indices to compute SEDs for. If None, uses all models.
av (array-like or float, optional) – A(V) values for each model. If None, defaults to 0.
rv (array-like or float, optional) – R(V) values for each model. If None, defaults to 3.3.
return_flux (bool, optional) – If True, return fluxes instead of magnitudes. Default is False.

Returns:

seds (numpy.ndarray of shape (Nmodels, Nbands)) – Computed SEDs at 1 kpc in apparent magnitudes (or maggies when return_flux=True).
rvecs (numpy.ndarray of shape (Nmodels, Nbands)) – Reddening derivatives \(\partial m / \partial A(V)\) in mag per mag of \(A(V)\) (or \(\partial F / \partial A(V)\) in maggies per mag when return_flux=True).
drvecs (numpy.ndarray of shape (Nmodels, Nbands)) – Scaled differential reddening vectors \(A(V)^{-1} \partial m / \partial R(V)\) in mag per mag of \(A(V)\) (or \(A(V)^{-1} \partial F / \partial R(V)\) in maggies per mag when return_flux=True).

Raises:

ValueError – If av or rv is an array whose shape is not (Nmodels,).

See also

brutus.core.sed_utils.get_seds: Underlying SED computation
StarGrid.get_seds: Single star SED generation

Notes

This method performs batch SED computation for multiple grid points simultaneously, which is more efficient than calling StarGrid.get_seds() repeatedly.

The scale factor relates to distance as \(s = 1/d^2\) where d is distance in kpc. The reddening is applied as:

\[m(\lambda) = m_0(\lambda) + A_V \cdot [r_0(\lambda) + R_V \cdot dr(\lambda)]\]

where \(r_0\) and \(dr\) are the reddening vector components.

loglike_grid(data, data_err, data_mask, avlim=(0.0, 20.0), av_gauss=(0.0, 1000000.0), rvlim=(1.0, 8.0), rv_gauss=(3.32, 0.18), av_init=None, rv_init=None, dim_prior=True, ltol=0.001, init_thresh=1e-06, parallax=None, parallax_err=None, return_vals=False, indices=None)[source]¶

Compute log-likelihood over the stellar model grid.

Optimizes (scale, A(V), R(V)) for each model in the instance’s grid and evaluates the corresponding maximum-likelihood fit to the data.

Parameters:

data (~numpy.ndarray of shape (Nfilt,)) – Observed linear flux densities in maggies (10**(-0.4 * mag)).
data_err (~numpy.ndarray of shape (Nfilt,)) – Flux density measurement errors in maggies.
data_mask (~numpy.ndarray of shape (Nfilt,)) – Boolean mask of valid bands (or 0/1 values, converted to bool).
indices (array-like, optional) – Subset of model indices to use. If None, uses all models.
avlim (tuple, optional) – (min, max) bounds on A(V). Default is AVLIM, (0, 20).
av_gauss (tuple, optional) – (mean, std) for Gaussian prior on A(V). Default is (0.0, 1e6).
rvlim (tuple, optional) – (min, max) bounds on R(V). Default is RVLIM, (1, 8).
rv_gauss (tuple, optional) – (mean, std) for Gaussian prior on R(V). Default is (3.32, 0.18).
av_init (float or ~numpy.ndarray of shape (Nmodels,), optional) – Starting A(V) and R(V) of the magnitude-space fit, one value for every model or one per model. Default is the mean of av_gauss and of rv_gauss.
rv_init (float or ~numpy.ndarray of shape (Nmodels,), optional) – Starting A(V) and R(V) of the magnitude-space fit, one value for every model or one per model. Default is the mean of av_gauss and of rv_gauss.
dim_prior (bool, optional) – Subtract 0.5 * (3 - Ndim) * ln(Ndim) from every log-likelihood, where Ndim is the number of valid bands. It changes only comparisons between objects with different numbers of bands. Default is True.
ltol (float, optional) – Convergence tolerance of both fits: a model stops when the next Gauss-Newton step is predicted to raise its log-posterior in (A(V), R(V)) by less than ltol. Must be positive. Default is 1e-3.
init_thresh (float, optional) – Models whose magnitude-space fit, with the parallax term, is below init_thresh times the best are not refined in flux space and keep their magnitude-space fit. Default is 1e-6.
parallax (float, optional) – Parallax measurement in mas.
parallax_err (float, optional) – Parallax error in mas.
return_vals (bool, optional) – If True, also return the fitted parameters and their precision matrices. Default is False.

Returns:

lnl (~numpy.ndarray of shape (Nmodels,)) – Log-likelihood of each model at its fit, with the Gaussian normalization and, if dim_prior, the dimensional term. It does not include the parallax.
Ndim (int) – Number of valid bands.
chi2 (~numpy.ndarray of shape (Nmodels,)) – Chi-square of each model at its fit, without the parallax.
scale (~numpy.ndarray of shape (Nmodels,)) – Fitted scale, 1 / d**2 with d in kpc. Only with return_vals=True, as are the arrays below.
av, rv (~numpy.ndarray of shape (Nmodels,)) – Fitted A(V) in mag and R(V).
icov_sar (~numpy.ndarray of shape (Nmodels, 3, 3)) – Precision matrix of (scale, A(V), R(V)) at the fit.
Nmodels is len(indices) when indices is given. A model below
init_thresh keeps its magnitude-space fit; the others are refined in
flux space.

Raises:

ValueError – If the shapes of data, data_err and data_mask disagree, if only one of parallax and parallax_err is given, if init_thresh is not in (0, 1), if ltol is not positive, or if no band is valid.

Warns:

RuntimeWarning – If fewer than 4 bands are valid, or if a model’s fit stopped on the step or damping limit before it converged.

See also

logpost_grid: Compute log-posteriors from likelihoods
_optimize_fit_mag: Magnitude-space optimization

Notes

This method optimizes (distance, Av, Rv) for each grid point by:

A magnitude-space fit of every model.
A cut on that fit (init_thresh), with the parallax term if a parallax is given.
A flux-space refinement of the kept models.

Both fits take Gauss-Newton steps in (A(V), R(V)) with the distance profiled out. A step that raises chi2 plus the prior penalties is retried with ten times more Levenberg-Marquardt damping. A model stops when the predicted gain is below ltol, when a step no longer moves it, when the damping passes 1e6, or after 100 steps.

The optimization uses priors on Av and Rv but not on distance/scale (distance priors are applied in logpost_grid).

logpost_grid(results, parallax=None, parallax_err=None, coord=None, Nmc_prior=100, lnprior=None, wt_thresh=0.001, cdf_thresh=0.002, max_models=50000, precision_shrinkage=0.0, lngalprior=None, lndustprior=None, dustfile=None, dlabels=<model grid labels>, avlim=(0.0, 20.0), rvlim=(1.0, 8.0), mem_lim=8000.0, rstate=None, apply_av_prior=True, R_solar=8.2, Z_solar=0.025, return_precision=False)[source]¶

Compute log-posterior over the stellar model grid.

Applies model priors, parallax constraints, Galactic structure and dust-map priors on top of the grid log-likelihoods, integrating over (distance, A(V), R(V)) uncertainty via Monte Carlo sampling.

Parameters:

results (tuple) – The 7-tuple of loglike_grid with return_vals=True. If it was computed with indices, lnprior and dlabels must be subset with the same indices, and the returned sel indexes that subset.
parallax (float, optional) – Parallax and its error, in mas. Both or neither. The selection weight of each model includes the parallax likelihood averaged over the model’s Gaussian in ln(d).
parallax_err (float, optional) – Parallax and its error, in mas. Both or neither. The selection weight of each model includes the parallax likelihood averaged over the model’s Gaussian in ln(d).
coord (tuple of 2 floats, optional) – Galactic (l, b) of the object, in degrees. The Galactic prior applies only when it is given, and the dust-map prior needs it.
Nmc_prior (int, optional) – Number of Monte Carlo draws per model, lowered if needed to fit mem_lim. Default is 100.
lnprior (~numpy.ndarray of shape (Nmodels,), optional) – Log-prior of each model. Default is 0 for every model.
wt_thresh (float or None, optional) – Keep the models whose selection weight is above wt_thresh times the largest. Default is 1e-3.
cdf_thresh (float or None, optional) – Used only if wt_thresh is None: keep the models of highest weight that hold 1 - cdf_thresh of the total. If both are None, every model is kept. Default is 2e-3.
max_models (int or None, optional) – Most models kept after the threshold. Above it, models are kept by systematic sampling with probability proportional to their selection weight, and ln_incl holds the log-probabilities. None sets no limit. Default is 50000.
precision_shrinkage (float, optional) – Fraction by which the off-diagonal terms of each (scale, A(V), R(V)) precision matrix shrink before the draws. Default is 0.0.
lngalprior (callable, optional) – Galactic prior f(dists, coord, labels=...), with dists in kpc. Default is ~brutus.priors.logp_galactic_structure with R_solar and Z_solar.
lndustprior (callable, optional) – Dust-map prior f(avs, dustmap, coord, distance=...). Default is ~brutus.priors.logp_extinction, given avlim unless it is a functools.partial that binds its own.
dustfile (str, path-like or ~brutus.dust.Bayestar, optional) – 3D dust map for the dust-map prior. A path is loaded as a ~brutus.dust.Bayestar map.
dlabels (structured ~numpy.ndarray of shape (Nmodels,) or None, optional) – Model labels passed to the Galactic prior. By default, the labels of the loaded model grid. With None, the prior gets labels=None.
avlim (tuple, optional) – (min, max) bounds on A(V) and R(V). A draw outside them gets zero weight. Defaults are AVLIM, (0, 20), and RVLIM, (1, 8).
rvlim (tuple, optional) – (min, max) bounds on A(V) and R(V). A draw outside them gets zero weight. Defaults are AVLIM, (0, 20), and RVLIM, (1, 8).
mem_lim (float, optional) – Memory limit for the draws, in MB. Default is 8000.0.
rstate (~numpy.random.RandomState or ~numpy.random.Generator, optional) – Random state for the draws and the max_models sampling.
apply_av_prior (bool, optional) – Whether to apply the dust-map extinction prior when dustfile is provided. Default is True.
R_solar (float, optional) – Galactocentric radius and height of the Sun, in kpc, for the default Galactic prior. Defaults are 8.2 and 0.025.
Z_solar (float, optional) – Galactocentric radius and height of the Sun, in kpc, for the default Galactic prior. Defaults are 8.2 and 0.025.
return_precision (bool, optional) – Return the (scale, A(V), R(V)) precision matrices of the selected models in place of their covariances, which skips the inversion. Default is False.

Returns:

sel (~numpy.ndarray of shape (Nsel,)) – Indices of the selected models in results.
cov_sar (~numpy.ndarray of shape (Nsel, 3, 3)) – Covariances of (scale, A(V), R(V)) of the selected models, or their precisions with return_precision=True.
lnp (~numpy.ndarray of shape (Nsel,)) – Log-posterior of each selected model.
dist_mc (~numpy.ndarray of shape (Nsel, Nmc)) – Monte Carlo distance draws, in kpc.
a_mc, r_mc (~numpy.ndarray of shape (Nsel, Nmc)) – Monte Carlo A(V) and R(V) draws.
lnp_mc (~numpy.ndarray of shape (Nsel, Nmc)) – Log-posterior weight of each draw: the model’s log-likelihood and lnprior, the log-distance Jacobian, and the Galactic, dust and parallax terms at the draw. To pool draws across models, subtract ln_incl[:, None].
mc_ess (~numpy.ndarray of shape (Nsel,)) – 1 / sum(w**2) of the normalized weights of each model’s draws. The draws are antithetic pairs, so this measures how concentrated the weights are, not a count of independent draws.
ln_incl (~numpy.ndarray of shape (Nsel,)) – Log-probability that each selected model was kept by the max_models subsampling; 0 when there was none. Use lnp - ln_incl as the log-weights to sum (evidence) or to resample models.

Raises:

ValueError – If only one of parallax and parallax_err is given, if dustfile is given with apply_av_prior but without coord, if max_models is neither None nor a positive integer, or if lnprior or dlabels (with coord) does not have one entry per model in results.

fit(data, data_err, data_mask, data_labels, save_file, phot_offsets=None, parallax=None, parallax_err=None, Nmc_prior=50, avlim=(0.0, 20.0), av_gauss=None, rvlim=(1.0, 8.0), rv_gauss=(3.32, 0.18), lnprior=None, wt_thresh=0.001, cdf_thresh=0.002, Ndraws=250, apply_agewt=True, apply_grad=True, lngalprior=None, lndustprior=None, dustfile=None, apply_dlabels=True, data_coords=None, logl_dim_prior=True, ltol=0.001, logl_initthresh=1e-06, mag_max=50.0, merr_max=0.25, rstate=None, save_dar_draws=True, running_io=True, mem_lim=8000.0, max_models=50000, precision_shrinkage=0.0, verbose=True, R_solar=8.2, Z_solar=0.025)[source]¶

Fit all input models to the input data to compute log-posteriors.

This is the main interface for fitting stellar parameters using grid-based Bayesian inference. Results are saved to an HDF5 file.

Parameters:

data (numpy.ndarray of shape (Ndata, Nfilt)) – Observed linear flux densities in maggies (10**(-0.4 * mag)) for each object.
data_err (numpy.ndarray of shape (Ndata, Nfilt)) – Associated flux density errors in maggies.
data_mask (numpy.ndarray of shape (Ndata, Nfilt)) – Boolean mask of valid measurements (or 0/1 values, converted to bool).
data_labels (numpy.ndarray of shape (Ndata,) or (Ndata, Nlabels)) – Identifier labels for each object stored in the labels dataset of the output file.
save_file (str) – Path to the output HDF5 file (.h5 is appended if absent).
phot_offsets (numpy.ndarray of shape (Nfilt,), optional) – Dimensionless multiplicative photometric offsets applied to data and data_err.
parallax (numpy.ndarray of shape (Ndata,), optional) – Parallax measurements in milliarcseconds (mas) for each object.
parallax_err (numpy.ndarray of shape (Ndata,), optional) – Parallax measurement errors in mas. Required when parallax is provided.
Nmc_prior (int, optional) – Number of Monte Carlo samples for prior integration. Default is 50.
avlim (tuple, optional) – (min, max) bounds on A(V). Default is AVLIM, (0, 20).
av_gauss (tuple, optional) – (mean, std) for a Gaussian prior on A(V) applied during the per-model likelihood optimization. Default is None, which maps to the essentially-flat (0.0, 1e6). Note this is applied in addition to any dust-map prior from dustfile (they are not exclusive), so passing an informative av_gauss together with dustfile counts the A(V) constraint twice.
rvlim (tuple, optional) – (min, max) bounds on R(V). Default is RVLIM, (1, 8).
rv_gauss (tuple, optional) – (mean, std) for Gaussian prior on R(V). Default is (3.32, 0.18) based on Schlafly et al. (2016).
lnprior (numpy.ndarray of shape (Nmodel,), optional) – Log-prior for each model. If not provided, defaults to Kroupa IMF prior on initial mass (for MIST models) or PS1 luminosity function prior (for Bayestar models).
wt_thresh (float, optional) – Threshold wt_thresh * max(weight) for model selection. Default is 1e-3.
cdf_thresh (float, optional) – CDF threshold for model selection (used if wt_thresh is None). Default is 2e-3.
Ndraws (int, optional) – Number of posterior draws to save per object. Default is 250.
apply_agewt (bool, optional) – Whether to add ln(agewt) to the model prior. The grid must then have an agewt label, in labels or params. Default is True.
apply_grad (bool, optional) – Whether to apply grid spacing corrections. Default is True.
lngalprior (callable, optional) – Galactic structure prior function with signature f(dists, coord, labels=None). If not provided, uses ~brutus.priors.logp_galactic_structure. Applied only when data_coords is given.
lndustprior (callable, optional) – Dust prior function with signature f(avs, dustmap, coord, distance=None). If not provided and dustfile is given, uses logp_extinction.
dustfile (str, path-like, or ~brutus.dust.Bayestar, optional) – 3D dust map for extinction priors. Can be a file path (string or os.PathLike, e.g. pathlib.Path) to a Bayestar HDF5 file, which will be loaded automatically, or a pre-loaded Bayestar object.
apply_dlabels (bool, optional) – Whether to pass model labels to the Galactic prior. With False, the prior gets labels=None. Default is True.
data_coords (numpy.ndarray of shape (Ndata, 2), optional) – Galactic (l, b) coordinates for each object in degrees. Required if using default Galactic prior.
logl_dim_prior (bool, optional) – Whether to apply dimensional correction to log-likelihood. Default is True.
ltol (float, optional) – Convergence tolerance of the fits, as in loglike_grid. Default is 1e-3.
logl_initthresh (float, optional) – Models whose magnitude-space fit is below this fraction of the best are not refined (init_thresh of loglike_grid). Default is 1e-6.
mag_max (float, optional) – Maximum allowed magnitude for valid data. Default is 50.0.
merr_max (float, optional) – Maximum allowed magnitude error. Default is 0.25.
rstate (numpy.random.RandomState, optional) – Random state for reproducibility.
save_dar_draws (bool, optional) – Whether to save distance, A(V), and R(V) draws. Default is True.
running_io (bool, optional) – If True, writes results to disk after each object (safer for long runs). If False, accumulates results in memory and writes at end (faster for slow filesystems). Default is True.
mem_lim (float, optional) – Memory limit in MB for Monte Carlo sampling. Default is 8000.0.
verbose (bool, optional) – Whether to print progress to stderr. Default is True.
max_models (int, optional) – When more models than this are selected for an object, about max_models of them are kept at random, each with a probability proportional to its weight (at most 1), to bound memory and runtime. The evidence and the posterior draws divide each kept model’s weight by that probability, so they stay unbiased. None keeps every selected model. Default is 50000.
precision_shrinkage (float, optional) – Fractional shrinkage applied to the off-diagonal terms of the 3x3 (scale, A(V), R(V)) precision matrix before Monte Carlo sampling, which stabilizes strongly-correlated fits. 0 disables it; ~0.03 is a reasonable nonzero value. Default is 0.0.
R_solar (float, optional) – Solar galactocentric radius in kpc, used by the Galactic structure prior. Default is 8.2.
Z_solar (float, optional) – Solar height above the Galactic midplane in kpc, used by the Galactic structure prior. Default is 0.025.

Returns:

Path to the output HDF5 file.

Return type:

str

Raises:

ValueError – If the shapes of data, data_err, data_mask, data_labels, parallax, parallax_err, data_coords or lnprior disagree, if only one of parallax and parallax_err is given, if logl_initthresh is not in (0, 1), if ltol is not positive, if max_models is neither None nor a positive integer, or if apply_agewt is True and the grid has no agewt label. These checks run before the output file is created.

Notes

The output HDF5 file contains the following datasets:

labels: Object labels (Ndata, Nlabels)
model_idx: Resampled model indices (Ndata, Ndraws)
ml_scale: Maximum likelihood scale factors (Ndata, Ndraws)
ml_av: Maximum likelihood A(V) values (Ndata, Ndraws)
ml_rv: Maximum likelihood R(V) values (Ndata, Ndraws)
ml_cov_sar: Covariance matrices (Ndata, Ndraws, 3, 3)
obj_log_post: Log-posteriors (Ndata, Ndraws)
obj_log_evid: Log-evidence per object (Ndata,)
obj_chi2min: Minimum chi-squared per object (Ndata,)
obj_Nbands: Number of bands used per object (Ndata,)
mc_ess: Monte Carlo weight-concentration index (Ndata, Ndraws). For each posterior draw, 1 / sum(w_i^2) over the MC integration weights for that draw’s source model. Low values indicate poor overlap between the Gaussian proposal and the target posterior. Because the MC draws are antithetic (correlated in pairs) this is a concentration diagnostic rather than a strict count of independent samples, and may read up to ~2x optimistic for well-mixed models.

An object with no valid bands is not fitted. Its row has model_idx = -99, obj_Nbands = 0 and NaN in every float field, and fit warns once with the number of such objects.

If save_dar_draws=True, also includes:

samps_dist: Distance draws in kpc (Ndata, Ndraws)
samps_red: A(V) draws (Ndata, Ndraws)
samps_dred: R(V) draws (Ndata, Ndraws)
samps_logp: Log-weights for draws (Ndata, Ndraws)

Examples

>>> from brutus.analysis import BruteForce
>>> from brutus.core import StarGrid
>>> from brutus.data import load_models
>>>
>>> # Load model grid
>>> models, labels, label_mask = load_models('grid.h5')
>>> grid = StarGrid(models, labels)
>>> fitter = BruteForce(grid)
>>>
>>> # Fit photometry
>>> results_file = fitter.fit(
...     phot, phot_err, phot_mask, obj_labels,
...     save_file='results.h5',
...     parallax=parallax, parallax_err=parallax_err,
...     data_coords=coords
... )

Population Analysis¶

brutus.analysis.isochrone_population_loglike(theta, stellarpop, obs_phot, obs_err, parallax=None, parallax_err=None, cluster_prob=0.95, dim_prior=True, outlier_model_func=None, smf_grid=None, eep_grid=None, mini_bound=0.08, eep_binary_max=480.0, corr_params=None, return_components=False, mask=None, **outlier_kwargs)[source]¶

Compute log-likelihood for a coeval stellar population using isochrone fitting.

Evaluates the population log-likelihood in five steps: 1. Generate isochrone grid over (mass, SMF) 2. Compute cluster and outlier likelihoods at each grid point 3. Apply mixture model at each grid point 4. Marginalize over (mass, SMF) with proper geometric factors 5. Sum over all objects

Parameters:

theta (array-like, shape (6,)) – Population parameters: [feh, loga, av, rv, dist, field_frac], with dist in pc and field_frac the field contamination fraction (0 to 1).
stellarpop (StellarPop object) – StellarPop model from core.populations module with get_seds method
obs_phot (array-like, shape (N_objects, N_filters)) – Observed flux densities in units of 10**(-0.4 * mag)
obs_err (array-like, shape (N_objects, N_filters)) – Flux density errors in same units
parallax (array-like, shape (N_objects,), optional) – Parallax measurements (mas)
parallax_err (array-like, shape (N_objects,), optional) – Parallax errors (mas)
cluster_prob (float or array-like of shape (N_objects,), optional) – Prior probability of cluster membership. Default 0.95
dim_prior (bool, optional) – Use chi-square (True) or normal (False) likelihood. Default True
outlier_model_func (callable, optional) – Custom outlier model function
smf_grid (array-like, optional) – Secondary mass fraction grid for binary modeling
eep_grid (array-like, optional) – EEP grid for isochrone evaluation
mini_bound (float, optional) – Minimum initial mass for isochrone. Default 0.08
eep_binary_max (float, optional) – Maximum EEP for binary modeling. Default 480.0
corr_params (array-like, optional) – Empirical correction parameters [dtdm, drdm, msto_smooth, feh_scale] forwarded to the isochrone grid generation (and from there to StellarPop.get_seds). Default None (model defaults).
return_components (bool, optional) – Return intermediate results for debugging. Default False
mask (array-like, shape (N_objects, N_filters), optional) – Data validity mask
**outlier_kwargs (dict) – Additional arguments passed to outlier model

Returns:

lnl_total (float) – Total log-likelihood summed over all objects; LOG_ZERO (-1e300) when it is not finite or no model in the isochrone grid is valid.
components (dict, optional) – Intermediate results (if return_components=True) containing: - ‘lnl_total’: same as primary return value - ‘lnl_per_object’: array of per-object likelihoods - ‘isochrone_grid’: the generated grid dictionary - ‘lnl_cluster’: cluster likelihood array - ‘lnl_outlier’: outlier likelihood array - ‘lnl_mixture’: mixed likelihood array

Raises:

ValueError – If theta does not have 6 entries; if obs_phot, obs_err and mask are not 2-D arrays of the same shape; if only one of parallax and parallax_err is given, or either does not have shape (N_objects,); or if smf_grid or eep_grid is invalid (see generate_isochrone_population_grid).

See also

generate_isochrone_population_grid: Step 1 - Grid generation
compute_isochrone_cluster_loglike: Step 2 - Cluster likelihood
compute_isochrone_outlier_loglike: Step 3 - Outlier likelihood
apply_isochrone_mixture_model: Step 4 - Mixture model
marginalize_isochrone_grid: Step 5 - Marginalization
brutus.core.populations.StellarPop: Stellar population model

Notes

This function implements the complete mixture-before-marginalization workflow:

\[\ln L(\theta) = \sum_i \ln \left[ \int \int \left( w_c P_c(d_i|m, s, \theta) + w_o P_o(d_i) \right) dm \, ds \right]\]

where: - \(\theta = [{\rm Fe/H}, \log {\rm age}, A_V, R_V, d, f_{\rm field}]\) are population parameters - \(m\) is stellar mass - \(s\) is secondary mass fraction (SMF) - \(w_c, w_o\) are mixture weights - \(P_c, P_o\) are cluster and outlier likelihoods - \(d_i\) is data for object i

The function is designed for use with MCMC samplers like emcee or nested sampling codes like dynesty. It returns LOG_ZERO (-1e300), with a warning, when no model in the isochrone grid is valid at theta.

Examples

Use with emcee for MCMC sampling:

>>> def lnprob(theta):
...     if not in_prior_bounds(theta):
...         return -np.inf
...     return isochrone_population_loglike(theta, pop, flux, err)
>>>
>>> sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob)
>>> sampler.run_mcmc(p0, nsteps)

Extract intermediate results for diagnostics:

>>> lnl, components = isochrone_population_loglike(
...     theta, pop, flux, err, return_components=True
... )
>>> print(f"Per-object likelihoods: {components['lnl_per_object']}")

brutus.analysis.generate_isochrone_population_grid(stellarpop, feh, loga, av, rv, dist, smf_grid=None, eep_grid=None, mini_bound=0.08, eep_binary_max=480.0, corr_params=None)[source]¶

Generate isochrone population grid over (mass, SMF) parameter space.

Parameters:

stellarpop (StellarPop object) – StellarPop model from core.populations module with get_seds method
feh (float) – Stellar population parameters: metallicity [Fe/H] (dex), log10(age/yr), extinction A(V) (mag), R(V) (dimensionless) and distance (pc)
loga (float) – Stellar population parameters: metallicity [Fe/H] (dex), log10(age/yr), extinction A(V) (mag), R(V) (dimensionless) and distance (pc)
av (float) – Stellar population parameters: metallicity [Fe/H] (dex), log10(age/yr), extinction A(V) (mag), R(V) (dimensionless) and distance (pc)
rv (float) – Stellar population parameters: metallicity [Fe/H] (dex), log10(age/yr), extinction A(V) (mag), R(V) (dimensionless) and distance (pc)
dist (float) – Stellar population parameters: metallicity [Fe/H] (dex), log10(age/yr), extinction A(V) (mag), R(V) (dimensionless) and distance (pc)
smf_grid (array-like, optional) – Secondary mass fraction grid. Default is 21 uniform points from 0.0 to 1.0
eep_grid (array-like, optional) – EEP grid for isochrone evaluation. Default is 1000 points from 202 to 808
mini_bound (float, optional) – Minimum initial mass for evaluation. Default 0.08 solar masses
eep_binary_max (float, optional) – Maximum EEP for binary modeling. Default 480.0
corr_params (array-like, optional) – Empirical correction parameters [dtdm, drdm, msto_smooth, feh_scale]

Returns:

grid – Dictionary containing:

'photometry': numpy.ndarray of shape (N_total_points, N_filters), model flux densities in maggies (10**(-0.4 * mag)).
'masses': numpy.ndarray of shape (N_total_points,), initial stellar masses in \(M_\odot\).
'smf_values': numpy.ndarray of shape (N_total_points,), secondary mass fractions in [0, 1].
'mass_jacobians': numpy.ndarray of shape (N_total_points,), mass grid spacings \(\Delta M\) in \(M_\odot\).
'smf_jacobians': numpy.ndarray of shape (N_total_points,), SMF grid spacings \(\Delta({\rm SMF})\).
'grid_info': dict with 'smf_grid', 'eep_grid', and 'n_total_points'.

Return type:

dict

Raises:

ValueError – If smf_grid is not 1-D, finite, strictly increasing and within [0, 1], if eep_grid is not 1-D, finite and strictly increasing, or if no model in the grid is valid (for example, theta outside the isochrone grid).

See also

compute_isochrone_cluster_loglike: Use this grid for likelihood computation
StellarPop.get_seds: Underlying SED generation

Notes

The grid is constructed by:

Looping over SMF values (binary mass ratios, passed to StellarPop.get_seds as binary_fraction)
For each SMF, computing isochrone along EEP dimension
Extracting masses from the isochrone
Computing geometric jacobians (grid spacings) for proper integration
Keeping the models with positive mass spacing, storing the post-main-sequence block once with the full SMF measure (see below), and leaving non-finite SEDs to the likelihood

The jacobians are critical for proper marginalization - they represent the geometric factors \(dm\) and \(d({\rm SMF})\) in the integral:

\[P({\rm data}) = \int \int P({\rm data}|m, {\rm SMF}) \, dm \, d({\rm SMF})\]

Binary models are computed only for eep <= eep_binary_max. Models above eep_binary_max are independent of SMF, so they are stored once (from the first SMF slice) with the total SMF measure (the sum of all SMF grid spacings) so that SMF marginalization weights them consistently with main-sequence models that accumulate over every slice. Their recorded smf_values entry is the slice they were generated from and is not meaningful for these SMF-independent models.

brutus.analysis.compute_isochrone_cluster_loglike(obs_flux, obs_err, isochrone_grid, parallax=None, parallax_err=None, distance=None, dim_prior=True, mask=None)[source]¶

Compute cluster membership likelihood using existing photometry infrastructure.

Parameters:

obs_flux (array-like, shape (N_objects, N_filters)) – Observed linear flux densities in maggies (10**(-0.4 * mag)).
obs_err (array-like, shape (N_objects, N_filters)) – Flux density errors in maggies.
isochrone_grid (dict) – Isochrone grid from generate_isochrone_population_grid()
parallax (array-like, shape (N_objects,), optional) – Parallax measurements (mas)
parallax_err (array-like, shape (N_objects,), optional) – Parallax errors (mas)
distance (float, optional) – Population distance (pc). Required if parallax provided.
dim_prior (bool, optional) – Whether to use chi-square (True) or normal (False) likelihood
mask (array-like, shape (N_objects, N_filters), optional) – Data mask (1=use, 0=skip)

Returns:

lnl_cluster – Cluster membership log-likelihood for each grid point and object. Invalid models (NaN photometry) are assigned NaN likelihood.

Return type:

array-like, shape (N_grid_points, N_objects)

Raises:

ValueError – If only one of parallax and parallax_err is given, or if parallax is given without distance.

See also

brutus.utils.photometry.phot_loglike: Photometric likelihood function
compute_isochrone_outlier_loglike: Complementary outlier likelihood
generate_isochrone_population_grid: Creates the isochrone_grid input

Notes

The likelihood includes both photometric and parallax components:

\[\ln L_{\rm cluster} = \ln L_{\rm phot} + \ln L_{\rm parallax}\]

where the photometric likelihood uses either chi-square (dim_prior=True) or Gaussian (dim_prior=False) formulation.

Invalid models with NaN photometry are preserved as NaN in the output to be properly handled during marginalization (they contribute zero probability via logsumexp).

brutus.analysis.compute_isochrone_outlier_loglike(obs_flux, obs_err, isochrone_grid=None, parallax=None, parallax_err=None, dim_prior=True, outlier_model_func=None, mask=None, **outlier_kwargs)[source]¶

Compute outlier likelihood with stellar-parameter-aware interface.

Parameters:

obs_flux (array-like, shape (N_objects, N_filters)) – Observed linear flux densities in maggies (10**(-0.4 * mag)).
obs_err (array-like, shape (N_objects, N_filters)) – Flux density errors in maggies.
isochrone_grid (dict, optional) – Isochrone grid containing stellar parameters for potential dependence
parallax (array-like, shape (N_objects,), optional) – Parallax measurements (mas)
parallax_err (array-like, shape (N_objects,), optional) – Parallax errors (mas)
dim_prior (bool, optional) – Use chi-square (True) or uniform (False) outlier model
outlier_model_func (callable, optional) – Custom outlier model function. Must accept a mask keyword (the same band mask seen by the cluster likelihood) in addition to stellar_params, parallax, and parallax_err.
mask (array-like, shape (N_objects, N_filters), optional) – Data mask (1=use, 0=skip). Should be the same mask used for the cluster likelihood so both mixture components describe the same data.
**outlier_kwargs (dict) – Additional arguments for outlier model

Returns:

lnl_outlier – Outlier likelihood for each grid point and object

Return type:

array-like, shape (N_grid_points, N_objects)

brutus.analysis.apply_isochrone_mixture_model(lnl_cluster, lnl_outlier, cluster_prob, field_fraction)[source]¶

Combine cluster and outlier log-likelihoods at each isochrone grid point.

For each (grid_point, object) pair: P(data|mass,SMF) = P_cluster * P(data|cluster) + P_outlier * P(data|outlier)

Parameters:

lnl_cluster (array-like, shape (N_grid_points, N_objects)) – Cluster membership likelihoods
lnl_outlier (array-like, shape (N_grid_points, N_objects)) – Outlier model likelihoods
cluster_prob (float or array-like of shape (N_objects,)) – Prior probability of cluster membership (external)
field_fraction (float) – Fraction of cluster stars that are field contaminants (fitted parameter)

Returns:

lnl_mixture – Mixed log-likelihoods for each grid point and object

Return type:

array-like, shape (N_grid_points, N_objects)

See also

marginalize_isochrone_grid: Next step after mixture application
compute_isochrone_cluster_loglike: Cluster likelihood component
compute_isochrone_outlier_loglike: Outlier likelihood component

Notes

The mixture model is applied as:

\[P({\rm data}|m, {\rm SMF}) = w_c \cdot P_c + w_o \cdot P_o\]

where: - \(w_c = P_{\rm cluster} \cdot (1 - f_{\rm field})\) - \(w_o = 1 - w_c\) - \(P_{\rm cluster}\) is the prior probability (cluster_prob) - \(f_{\rm field}\) is the field contamination fraction

This is computed in log-space using logsumexp for numerical stability:

\[\ln L_{\rm mix} = \ln(\exp(\ln L_c + \ln w_c) + \exp(\ln L_o + \ln w_o))\]

Applying the mixture at each (mass, SMF) grid point before marginalizing evaluates the cluster and outlier hypotheses under the same per-star integral.

brutus.analysis.marginalize_isochrone_grid(lnl_mixture, mass_jacobians, smf_jacobians)[source]¶

Marginalize mixed likelihoods over (mass, SMF) grid with geometric jacobians.

Performs: P(data|population_params) = ∫∫ P(data|mass,SMF) dm d(SMF)

Parameters:

lnl_mixture (array-like, shape (N_grid_points, N_objects)) – Mixed likelihoods at each grid point
mass_jacobians (array-like, shape (N_grid_points,)) – Mass grid spacing (geometric factors for integration)
smf_jacobians (array-like, shape (N_grid_points,)) – SMF grid spacing (geometric factors for integration)

Returns:

lnl_marginalized – Marginalized log-likelihoods for each object; LOG_ZERO (-1e300) for a zero likelihood.

Return type:

array-like, shape (N_objects,)

See also

apply_isochrone_mixture_model: Previous step before marginalization
generate_isochrone_population_grid: Provides the jacobians

Notes

Performs the integration:

\[P({\rm data}|\theta) = \int \int P({\rm data}|m, {\rm SMF}, \theta) \, dm \, d({\rm SMF})\]

numerically using a grid-based approach:

\[\ln P \approx \ln \frac{\sum_{i,j} \exp(\ln L_{i,j}) \cdot \Delta m_i \cdot \Delta({\rm SMF})_j}{\sum_{i,j} \Delta m_i \cdot \Delta({\rm SMF})_j}\]

where: - \(\ln L_{i,j}\) is the mixed likelihood at grid point (i,j) - \(\Delta m_i\) is the mass grid spacing (mass_jacobians) - \(\Delta({\rm SMF})_j\) is the SMF grid spacing (smf_jacobians)

Invalid models (with NaN likelihood) are set to -np.inf before logsumexp so that they contribute zero probability. Normalizing the measure over the valid grid points (the denominator above) makes \(p(m, {\rm SMF} \mid \theta)\) a proper uniform prior because an unnormalized \(\theta\)-dependent grid volume would multiply every mixture component, including the \(\theta\)-independent outlier model, and bias the inferred population parameters and field fraction.

The jacobians represent geometric integration weights and are crucial for obtaining unbiased parameter estimates.

Photometric Offsets¶

class brutus.analysis.PhotometricOffsetsConfig(min_bands_used: int = 4, min_bands_unused: int = 3, n_bootstrap: int = 300, uncertainty_method: str = 'bootstrap_iqr')[source]¶

Bases: object

Configure band-count thresholds and bootstrap settings for photometric_offsets.

Parameters:

min_bands_used (int, optional) – Minimum number of bands, excluding the band being calibrated, required for objects where the current band was used in fitting. Default is 4 (i.e. at least 5 observed bands in total, matching the legacy > 3 + 1 cut).
min_bands_unused (int, optional) – Minimum number of observed bands required for objects where the current band was not used in fitting. Default is 3. (The legacy v0.8.3 pipeline effectively required 4; pass min_bands_unused=4 to reproduce it.)
n_bootstrap (int, optional) – Number of bootstrap realizations for uncertainty estimation. Default is 300.
uncertainty_method (str, optional) – Method for uncertainty estimation. Options: - ‘bootstrap_std’: Standard deviation of bootstrap medians - ‘bootstrap_iqr’: Scaled interquartile range of bootstrap medians Default is ‘bootstrap_iqr’.

See also

photometric_offsets: Main function using this configuration

Notes

The configuration defaults are chosen to balance statistical robustness with computational efficiency:

min_bands_used=4: Keeps at least 4 bands without the band being calibrated, enough for a fit of 3 parameters with 1 degree of freedom
min_bands_unused=3: Minimum for meaningful photometric constraints
n_bootstrap=300: Sufficient for stable uncertainty estimates
bootstrap_iqr: More robust to outliers than standard deviation

Examples

>>> config = PhotometricOffsetsConfig(min_bands_used=5, n_bootstrap=500)
>>> offsets, errors, n_used = photometric_offsets(
...     phot, err, mask, models, idxs, avs, rvs, dists,
...     config=config, rng=np.random.default_rng(42)
... )

__init__(min_bands_used: int = 4, min_bands_unused: int = 3, n_bootstrap: int = 300, uncertainty_method: str = 'bootstrap_iqr')[source]¶

brutus.analysis.photometric_offsets(phot: ndarray, err: ndarray, mask: ndarray, models: ndarray, idxs: ndarray, reds: ndarray, dreds: ndarray, dists: ndarray, sel: ndarray | None = None, weights: ndarray | None = None, mask_fit: ndarray | None = None, old_offsets: ndarray | None = None, prior_mean: ndarray | None = None, prior_std: ndarray | None = None, verbose: bool = True, config: PhotometricOffsetsConfig | None = None, rng: Generator | None = None) → Tuple[ndarray, ndarray, ndarray][source]¶

Compute multiplicative photometric offsets between data and models.

This function computes photometric offsets that account for systematic differences between observed photometry and model predictions. The offsets are computed by comparing model/data flux ratios across a sample of objects, with proper uncertainty estimation and optional prior constraints.

Parameters:

phot (np.ndarray of shape (n_objects, n_filters)) – Observed linear flux densities in maggies (10**(-0.4 * mag)) for all objects. Entries in masked bands may be non-finite placeholders; observed entries must be finite, and non-positive observed fluxes are excluded per band.
err (np.ndarray of shape (n_objects, n_filters)) – Associated flux density errors in maggies for all objects. Must be positive in observed bands; masked entries are ignored.
mask (np.ndarray of shape (n_objects, n_filters)) – Binary mask (0/1, bool or numeric) indicating observed bands for each object.
models (np.ndarray of shape (n_models, n_filters, 3)) – Magnitude polynomial coefficients (m_0, R_0, dR/dR_V) at 1 kpc.
idxs (np.ndarray of shape (n_objects, n_samples)) – Posterior-resampled grid model indices for each object.
reds (np.ndarray of shape (n_objects, n_samples)) – Posterior \(A(V)\) extinction draws in magnitudes for each object.
dreds (np.ndarray of shape (n_objects, n_samples)) – Posterior \(R(V)\) reddening parameter draws (dimensionless) for each object.
dists (np.ndarray of shape (n_objects, n_samples)) – Posterior distance draws in kiloparsecs (kpc) for each object.
sel (np.ndarray of bool or 0/1, shape (n_objects), optional) – Selection of objects to use. Default uses all objects.
weights (np.ndarray of shape (n_objects, n_samples), optional) – Sample weights for each object, finite and >= 0. Default uses uniform weights.
mask_fit (np.ndarray of shape (n_filters), optional) – Boolean mask indicating which filters were used in fitting. Default assumes all filters were used.
old_offsets (np.ndarray of shape (n_filters), optional) – Offsets the fit applied, finite and > 0. They multiply the fluxes and errors in the leave-one-out weights, so the weights use the fluxes of the fit. Default is no previous offsets.
prior_mean (np.ndarray of shape (n_filters), optional) – Gaussian prior means for offsets, finite. Must be provided with prior_std.
prior_std (np.ndarray of shape (n_filters), optional) – Gaussian prior standard deviations for offsets, finite and > 0.
verbose (bool, optional) – Whether to print progress information. Default is True.
config (PhotometricOffsetsConfig, optional) – Configuration object with analysis parameters. Default uses standard configuration.
rng (np.random.Generator, optional) – Random number generator, the only source of randomness. Pass a seeded one for reproducible results. Default is an unseeded generator.

Returns:

offsets (np.ndarray of shape (n_filters)) – Multiplicative photometric offsets (model/data ratios).
offset_errors (np.ndarray of shape (n_filters)) – Uncertainties on the photometric offsets.
n_objects_used (np.ndarray of shape (n_filters)) – Number of objects that actually informed each offset (objects with non-positive fluxes or fully degenerate weights are excluded).

Raises:

TypeError – If one of the required input arrays is not a numpy array.
ValueError – If an input has the wrong shape, dtype or values, or if only one of prior_mean and prior_std is given. All checks run before the model SEDs are generated.

Examples

>>> import numpy as np
>>> from brutus.analysis.offsets import photometric_offsets
>>>
>>> # Mock data for demonstration
>>> n_obj, n_filt, n_samp = 100, 5, 50
>>> phot = np.random.uniform(0.1, 10, (n_obj, n_filt))
>>> err = 0.1 * phot
>>> mask = np.random.choice([0, 1], (n_obj, n_filt), p=[0.1, 0.9])
>>>
>>> # Mock fitted parameters
>>> models = np.random.random((1000, n_filt, 3))
>>> idxs = np.random.randint(0, 1000, (n_obj, n_samp))
>>> reds = np.random.uniform(0, 2, (n_obj, n_samp))
>>> dreds = np.random.uniform(2.5, 4.5, (n_obj, n_samp))
>>> dists = np.random.uniform(0.1, 10, (n_obj, n_samp))
>>>
>>> # Compute offsets
>>> offsets, errors, n_used = photometric_offsets(
...     phot, err, mask, models, idxs, reds, dreds, dists
... )
>>> print(f"Computed offsets: {offsets}")

See also

PhotometricOffsetsConfig: Configuration options
brutus.core.sed_utils.get_seds: Model SED generation
brutus.utils.photometry.phot_loglike: Likelihood computation
brutus.analysis.individual.BruteForce: Source of fitted parameters

Notes

The photometric offset for each band is computed as:

Generate model SEDs for all fitted objects and posterior samples
Scale by distance: \(F_{\rm model} = F_0 / d^2\)
Compute flux ratios: \(r = F_{\rm model} / F_{\rm obs}\)
Reweight samples: For bands used in fitting, importance-reweight the full-posterior draws to \(P(M|D_{-i})\) (band i excluded) with weights \(w \propto L_{-i}/L_{\rm full}\) to avoid circularity
Bootstrap: Resample objects and models with weights, compute median
Uncertainty: From bootstrap distribution (IQR or std)
Apply priors (optional): Bayesian combination with prior

Reweighting in step 4 prevents circularity when calibrating a band that was included in the original fit. Because the input draws target the full posterior \(P(M \mid D)\), weighting each draw by the likelihood ratio \(L_{-i} / L_{\rm full}\) (the reciprocal of band \(i\)’s likelihood up to normalization) targets the leave-one-out posterior \(P(M \mid D_{-i})\), whereas weighting by \(L_{-i}\) alone would only partially remove the dependence on band \(i\). With Gaussian likelihoods, \(\ln w = \chi^2_i / 2\) up to a per-object constant, where \(\chi^2_i\) is band i’s chi-square. This holds for fits with and without logl_dim_prior, whose term is constant for one object. The parallax does not enter the weights.

The offsets should be applied as:

\[F_{\rm corrected} = F_{\rm observed} \times {\rm offset}\]

For iterative refinement, provide old_offsets from previous iteration.

References

The bootstrap methodology follows standard non-parametric uncertainty estimation. The likelihood-ratio importance reweighting ensures unbiased estimates for bands included in the original fit.

Line-of-Sight Dust¶

brutus.analysis.los_clouds_priortransform(u, rlims=(0.0, 6.0), dlims=(4.0, 19.0), pb_params=(-3.0, 0.7, -inf, 0.0), s_params=(-3.0, 0.3, -inf, 0.0), dust_template=False, nlims=(0.2, 2))[source]¶

Transform unit cube samples to physical parameters for LOS dust fitting.

The “prior transform” for the LOS fit that converts from draws on the N-dimensional unit cube to samples from the prior. Used in nested sampling methods. Assumes uniform priors for distance and reddening and a (truncated) log-normal in outlier fraction.

Parameters:

u (array_like, shape (Nparams,)) – The Nparams values drawn from the unit cube. Contains the portion of outliers P_b, followed by the foreground smoothing sfore and background smoothing sback, followed by the foreground reddening fred, followed by a series of (dist, red) pairs for each “cloud” along the LOS.
rlims (tuple of float, optional) – The reddening bounds within which we’d like to sample. Default is (0., 6.), which assumes reddening is in units of A_V.
dlims (tuple of float, optional) – The distance bounds within which we’d like to sample. Default is (4., 19.), which assumes distance is in units of distance modulus.
pb_params (tuple of float, optional) – Mean, standard deviation, lower bound, and upper bound for a truncated log-normal distribution used as a prior for the outlier model. The default is (-3., 0.7, -np.inf, 0.), which corresponds to a mean of 0.05, a standard deviation of a factor of 2, a lower bound of 0, and an upper bound of 1.
s_params (tuple of float, optional) – Mean, standard deviation, lower bound, and upper bound for a truncated log-normal distribution used as a prior for the smoothing along the reddening axis, as a fraction of the rlims range (the kernel width is s * (rlims[1] - rlims[0])). The default is (-3., 0.3, -np.inf, 0.), which corresponds to a mean of 0.05, a standard deviation of a factor of 1.35, a lower bound of 0, and an upper bound of 1.
dust_template (bool, optional) – Whether or not to use a spatial distribution for the dust based on a particular template. If true, dust along the line of sight will be in terms of rescalings of the template rather than A_V. Default is False.
nlims (tuple of float, optional) – Lower and upper bounds for the uniform prior for the rescaling applied to the Planck spatial reddening template. Default is (0.2, 2.).

Returns:

x – The transformed physical parameters in order: [pb, s_fore, s_back, A_V_fore, d1, A_V1, d2, A_V2, …] where di are sorted in increasing order.

Return type:

ndarray, shape (Nparams,)

Raises:

ValueError – If u is not 1-D with 4 + 2 * n entries (n >= 1) of finite values in [0, 1], or if a limit or hyperparameter is degenerate.

See also

los_clouds_loglike_samples: Likelihood function for these parameters
kernel_gauss: Gaussian smoothing kernel
kernel_lorentz: Lorentzian smoothing kernel

Notes

The prior distributions are:

Outlier fraction (P_b): Truncated log-normal with median ~0.05
Smoothing (s_fore, s_back): Truncated log-normal with median ~0.05
Foreground extinction: Uniform over rlims
Cloud distances: Uniform over dlims, sorted in increasing order
Cloud extinctions: Uniform over rlims, ordered by distance

The sorting ensures clouds are ordered by increasing distance, which is required by the likelihood’s per-segment extinction model.

The log-normal priors on outlier fraction and smoothing concentrate probability near zero while allowing occasional larger values.

The reddening amplitudes are drawn as independent uniforms (only the distances are sorted). When the likelihood is evaluated with monotonic=True (the default in los_clouds_loglike_samples), non-monotone amplitude orderings are rejected with LOG_ZERO. In the default mode the foreground and the n cloud amplitudes form one ordered chain, so only a fraction 1/(n+1)! of the sampled prior volume is retained (the nested-sampling evidence carries a -ln((n+1)!) offset); with a dust template or an additive foreground the foreground is not part of the ordering, so the retained fraction is 1/n! and the offset is -ln(n!). Add the corresponding ln((n+1)!) or ln(n!) back to each lnZ before comparing models with different numbers of clouds.

Examples

>>> import numpy as np
>>> # For 1-cloud model: pb, s0, s, fred, dist1, red1 (6 parameters)
>>> u = np.random.uniform(0, 1, 6)
>>> params = los_clouds_priortransform(u)
>>> print(f"Outlier fraction: {params[0]:.3f}")
>>> print(f"Cloud distance (DM): {params[4]:.1f}")
>>> print(f"Cloud extinction (A_V): {params[5]:.2f}")

brutus.analysis.los_clouds_loglike_samples(theta, dsamps, rsamps, kernel='gauss', rlims=(0.0, 6.0), template_reds=None, Ndraws=25, additive_foreground=False, monotonic=True)[source]¶

Compute log-likelihood for multi-cloud extinction model along line-of-sight.

Compute the log-likelihood for the cumulative reddening along the line of sight (LOS) parameterized by theta, given a set of input reddening and distance draws. Assumes an outlier model that is uniform in reddening over rlims, independent of distance.

Parameters:

theta (array_like, shape (Nparams,)) – A collection of parameters that characterizes the cumulative reddening along the LOS. Contains the fraction of outliers P_b followed by the fractional reddening smoothing for the foreground s0 and background s followed by the foreground reddening fred followed by a series of (dist, red) pairs for each “cloud” along the LOS. Each red is the absolute mean reddening level for stars between dist and the next cloud (values are not summed across clouds); if template_reds is provided the cloud values are instead dimensionless rescalings of each star’s template reddening, and if additive_foreground=True they are increments above the foreground fred.
dsamps (array_like, shape (Nobj, Nsamps)) – Distance samples for each object. Follows the units used in theta. Non-finite draws are handled as described for rsamps.
rsamps (array_like, shape (Nobj, Nsamps)) – Reddening samples for each object. Follows the units in theta. Each star’s likelihood is the average of the kernel over its finite draws (finite distance and reddening); a draw with a NaN or infinite value is left out. A finite distance outside [0, 1e10) counts as a draw with zero likelihood. A star with no finite draw is left out of the sum.
kernel (str or callable, optional) – The kernel used to weight the samples along the LOS. If a string is passed, a pre-specified kernel will be used. Options include ‘lorentz’, ‘gauss’, and ‘tophat’. The built-in kernels are renormalized for truncation to rlims (see the kernels’ rlims parameter), since the stellar reddening samples cannot scatter outside that range. A callable must accept (reds, (kmean, ksig)) where kmean/ksig are arrays broadcastable against reds; it is responsible for its own normalization. Default is ‘gauss’.
rlims (tuple of float, optional) – The reddening bounds within which we’d like to sample, used to normalize both the outlier model and the built-in kernels. Default is (0., 6.), which assumes reddening is in units of A_V.
template_reds (array_like, shape (Nobj,), optional) – Reddenings for each star based on a spatial dust template. If not provided, the same reddening value in a given distance bin will be fit to all stars. If provided, a rescaled version of the individual reddenings will be fit instead.
Ndraws (int, optional) – The number of draws to use for each star, at least 1; the first Ndraws columns of dsamps and rsamps are used. Default is 25.
additive_foreground (bool, optional) – Whether the foreground is treated as just another value or added to all background values. Default is False.
monotonic (bool, optional) –
Whether to enforce monotonicity in the fits so that the effective extinction values must get larger with distance. The constraint is applied after the template_reds/additive_foreground semantics are accounted for: with a template the (dimensionless) rescalings must be non-decreasing (the foreground, which is in A_V, is not compared against them); with an additive foreground the increments must be non-negative and non-decreasing. Default is True.

Note that with monotonic=True the sampled prior (los_clouds_priortransform) is not renormalized to the monotone region, so nested-sampling evidences carry a mode-dependent offset for an n-cloud model: -ln((n+1)!) in the default mode (foreground and clouds form one ordered chain), but -ln(n!) with a dust template or additive foreground (the foreground is not part of the ordering). Add the corresponding ln((n+1)!) or ln(n!) back before comparing evidences across different numbers of clouds.

Returns:

loglike – The computed log-likelihood; LOG_ZERO (-1e300) for a zero likelihood, or when it is NaN.

Return type:

float

Raises:

ValueError – If the shapes of theta, dsamps, rsamps or template_reds do not match, if pb is outside [0, 1], if s0 or s is not positive, if Ndraws is not an integer >= 1, if the cloud distances decrease, or if kernel is unknown.

Examples

>>> import numpy as np
>>> # Generate synthetic stellar samples
>>> nstars = 10
>>> nsamps = 50
>>> dsamps = np.random.uniform(6, 12, (nstars, nsamps))  # distance modulus
>>> rsamps = np.random.uniform(0, 2, (nstars, nsamps))   # A_V extinction
>>>
>>> # 1-cloud model parameters: [pb, s0, s, fred, dist1, red1]
>>> theta = [0.1, 0.05, 0.05, 0.2, 8.0, 0.5]
>>>
>>> # Compute likelihood
>>> loglike = los_clouds_loglike_samples(theta, dsamps, rsamps)
>>> print(f"Log-likelihood: {loglike:.2f}")

brutus.analysis.kernel_tophat(reds, kp, rlims=None)[source]¶

Compute weighted log-probabilities using a Top-Hat kernel.

Assigns uniform probability density to extinction samples within mean +/- half_bin_width and LOG_ZERO outside.

Parameters:

reds (array_like of shape (Nsamps,) or (Nobj, Nsamps)) – Reddening samples in magnitudes (\(A(V)\)).
kp (2-tuple of float or ndarray) – Kernel parameters (mean, half_bin_width) in magnitudes, broadcastable against reds.
rlims (tuple of float, optional) – If provided, the kernel support is clipped to [rlims[0], rlims[1]] and renormalized so the kernel integrates to 1 over the admissible reddening range. Without this, kernels centered near a boundary (e.g. A_V ~ 0) lose probability mass to values the reddening samples can never take, biasing boundary fits. Default is None (no truncation).

Returns:

logw – Log-weights for each sample.

Return type:

ndarray, shape (Nsamps,)

Examples

>>> import numpy as np
>>> reds = np.array([0.1, 0.3, 0.5, 0.7, 0.9])
>>> kp = (0.5, 0.2)  # mean=0.5, half-width=0.2, so range [0.3, 0.7]
>>> logw = kernel_tophat(reds, kp)
>>> # Samples within [0.3, 0.7] should have higher weights

brutus.analysis.kernel_gauss(reds, kp, rlims=None)[source]¶

Compute weighted log-probabilities using a Gaussian kernel.

Evaluates the log-density of a normal distribution centered on mean with standard deviation standard_deviation at each reddening sample.

Parameters:

reds (array_like of shape (Nsamps,) or (Nobj, Nsamps)) – Reddening samples in magnitudes (\(A(V)\)).
kp (2-tuple of float or ndarray) – Kernel parameters (mean, standard_deviation) in magnitudes, broadcastable against reds.
rlims (tuple of float, optional) – If provided, the kernel is renormalized for truncation to [rlims[0], rlims[1]] (divided by its probability mass over that range) so it is a proper density on the admissible reddening range. Without this, kernels centered near a boundary (e.g. A_V ~ 0) lose probability mass to values the reddening samples can never take, biasing boundary fits. Default is None (no truncation).

Returns:

logw – Log-weights for each sample.

Return type:

ndarray, shape (Nsamps,)

Examples

>>> import numpy as np
>>> reds = np.array([0.1, 0.3, 0.5, 0.7, 0.9])
>>> kp = (0.5, 0.1)  # mean=0.5, std=0.1
>>> logw = kernel_gauss(reds, kp)
>>> # Sample at 0.5 should have highest weight

brutus.analysis.kernel_lorentz(reds, kp, rlims=None)[source]¶

Compute weighted log-probabilities using a Lorentzian kernel.

Evaluates the log-density of a Cauchy (Lorentzian) distribution centered on mean with half-width at half-maximum HWHM at each reddening sample.

Parameters:

reds (array_like of shape (Nsamps,) or (Nobj, Nsamps)) – Reddening samples in magnitudes (\(A(V)\)).
kp (2-tuple of float or ndarray) – Kernel parameters (mean, HWHM) in magnitudes, broadcastable against reds.
rlims (tuple of float, optional) – If provided, the kernel is renormalized for truncation to [rlims[0], rlims[1]] (divided by its probability mass over that range) so it is a proper density on the admissible reddening range. Without this, kernels centered near a boundary (e.g. A_V ~ 0) lose probability mass to values the reddening samples can never take, biasing boundary fits. Default is None (no truncation).

Returns:

logw – Log-weights for each sample.

Return type:

ndarray, shape (Nsamps,)

Examples

>>> import numpy as np
>>> reds = np.array([0.1, 0.3, 0.5, 0.7, 0.9])
>>> kp = (0.5, 0.1)  # mean=0.5, HWHM=0.1
>>> logw = kernel_lorentz(reds, kp)
>>> # Sample at 0.5 should have highest weight, with heavy tails

Analysis Module (brutus.analysis)¶

Individual Star Fitting¶

Population Analysis¶

Photometric Offsets¶

Line-of-Sight Dust¶

Analysis Module (`brutus.analysis`)¶