#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
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
---------
magnitude : Convert flux to magnitudes
AB magnitude system conversion
inv_magnitude : Convert magnitudes to flux
Inverse of magnitude conversion
luptitude : Convert flux to asinh magnitudes
Asinh magnitude system (Lupton et al. 1999)
inv_luptitude : Convert asinh magnitudes to flux
Inverse of luptitude conversion
add_mag : Combine magnitudes
Add fluxes in magnitude space
phot_loglike : Photometric log-likelihood
Core likelihood function for stellar fitting
chisquare_outlier_loglike : Chi-square outlier model
Outlier likelihood for mixture models
uniform_outlier_loglike : Uniform 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
"""
import warnings
import numpy as np
from scipy.special import gammaln, xlogy
__all__ = [
"magnitude",
"inv_magnitude",
"luptitude",
"inv_luptitude",
"add_mag",
"phot_loglike",
"chisquare_outlier_loglike",
"uniform_outlier_loglike",
]
[docs]
def magnitude(phot, err, zeropoints=1.0):
"""
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`.
"""
# Compute magnitudes.
mag = -2.5 * np.log10(phot / zeropoints)
# Compute errors.
mag_err = 2.5 / np.log(10.0) * err / phot
return mag, mag_err
[docs]
def inv_magnitude(mag, err, zeropoints=1.0):
"""
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`.
"""
# Compute magnitudes.
phot = 10 ** (-0.4 * mag) * zeropoints
# Compute errors.
phot_err = err * 0.4 * np.log(10.0) * phot
return phot, phot_err
[docs]
def luptitude(phot, err, skynoise=1.0, zeropoints=1.0):
"""
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`.
"""
# Normalize photometry.
f = phot / zeropoints
df = err / zeropoints
b = skynoise / zeropoints
# Compute luptitudes.
lupt = -2.5 / np.log(10.0) * (np.arcsinh(0.5 * f / b) + np.log(b))
# Compute errors.
lupt_err = 2.5 / np.log(10.0) * df / (2.0 * b * np.sqrt(1.0 + (0.5 * f / b) ** 2))
return lupt, lupt_err
[docs]
def inv_luptitude(lupt, err, skynoise=1.0, zeropoints=1.0):
"""
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`.
"""
# Normalize sky noise.
b = skynoise / zeropoints
# Compute flux.
f = 2.0 * b * np.sinh(-0.4 * np.log(10.0) * lupt - np.log(b))
# Compute errors.
df = (
err
* 0.4
* np.log(10.0)
* 2.0
* b
* np.cosh(-0.4 * np.log(10.0) * lupt - np.log(b))
)
# Convert back to original units.
phot = f * zeropoints
phot_err = df * zeropoints
return phot, phot_err
[docs]
def add_mag(mag1, mag2):
"""
Add magnitudes.
Parameters
----------
mag1 : float or `~numpy.ndarray`
First set of magnitudes.
mag2 : float or `~numpy.ndarray`
Second set of magnitudes.
Returns
-------
mag_combined : float or `~numpy.ndarray`
Combined magnitudes corresponding to the combined flux from
`mag1` and `mag2`.
"""
# Compute combined flux.
flux_combined = 10 ** (-0.4 * mag1) + 10 ** (-0.4 * mag2)
# Convert back to magnitudes.
mag_combined = -2.5 * np.log10(flux_combined)
return mag_combined
[docs]
def phot_loglike(
flux,
err,
mfluxes,
mask=None,
dim_prior=False,
dof_reduction=0,
extra_chi2=None,
extra_dims=None,
):
r"""
Compute the log-likelihood between observed and model fluxes.
Parameters
----------
flux : `~numpy.ndarray` with shape (Nobj, Nfilt)
Observed flux values.
err : `~numpy.ndarray` with shape (Nobj, Nfilt)
Associated flux errors.
mfluxes : `~numpy.ndarray` with shape (Nobj, Nmod, Nfilt) or (Nmod, Nfilt)
Model fluxes (for each model). A 2-D array means every object is
compared against the same shared model grid; this triggers a much
faster matrix-multiplication evaluation of the chi-square that avoids
materializing any (Nobj, Nmod, Nfilt) intermediates.
mask : `~numpy.ndarray` with shape (Nobj, Nfilt), optional
Binary mask indicating whether each observed band can be
used (1) or should be skipped (0). Bands with non-finite flux/error
or non-positive error are always excluded (they contribute neither
to the chi-square nor to the effective dimensionality).
dim_prior : bool, optional
Whether to apply a dimensionality prior from the
chi-squared distribution. Default is `False`.
.. warning::
When dim_prior=True, perfect model matches (chi2≈0) can cause
problematic behavior: +inf likelihood for DOF=1, -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.
extra_chi2 : `~numpy.ndarray` with shape (Nobj,) or (Nobj, Nmod), optional
Additional chi-square contributions (e.g., from a parallax
measurement) added to the photometric chi-square *before* any
dimensionality prior is applied. This keeps auxiliary Gaussian
constraints in the same chi-square "measure" as the photometry when
`dim_prior=True` (mirroring how `BruteForce` counts parallax as an
extra band). Note that for `dim_prior=False` the Gaussian
normalization constant of the extra term is *not* added; callers
needing a proper density should add it themselves.
extra_dims : `~numpy.ndarray` with shape (Nobj,), optional
Additional dimensionality (number of extra data points) associated
with `extra_chi2`, added to the per-object effective DOF.
Returns
-------
lnl : `~numpy.ndarray` with shape (Nobj, Nmod)
Log-likelihood values.
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:
.. math::
\\ln L = -\\frac{1}{2} \\left[ \\chi^2 + N \\ln(2\\pi) + \\ln|\\Sigma| \\right]
where :math:`\\chi^2 = \\sum_i (f_i - m_i)^2/\\sigma_i^2` and
:math:`|\\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.
For the shared-model (2-D `mfluxes`) fast path, the chi-square is
expanded as :math:`\\sum w f^2 - 2 (w f) M^T + w (M^2)^T` and evaluated
with matrix products. This is algebraically identical to the direct
residual sum but reassociates the floating-point additions, so results
can differ from the 3-D path at machine precision.
"""
flux = np.asarray(flux)
err = np.asarray(err)
mfluxes = np.asarray(mfluxes)
# Initialize values.
Nobj, Nfilt = flux.shape[:2]
shared_models = mfluxes.ndim == 2
# Ensure proper dimensions.
if shared_models:
Nmod = mfluxes.shape[0]
if mfluxes.shape != (Nmod, Nfilt):
raise ValueError("Inconsistent dimensions between flux and mfluxes")
else:
Nmod = mfluxes.shape[1]
if mfluxes.shape != (Nobj, Nmod, Nfilt):
raise ValueError("Inconsistent dimensions between flux and mfluxes")
if mask is None:
mask = np.ones_like(flux)
# Fold data validity into the mask: non-finite fluxes/errors and
# non-positive errors cannot contribute a meaningful chi-square term.
# Without this, a single unmasked NaN flux poisons every model for that
# object, and a zero error yields +inf log-likelihood via log(var).
with np.errstate(invalid="ignore"):
valid = np.isfinite(flux) & np.isfinite(err) & (err > 0)
mask = np.asarray(mask) * valid
# Sanitize NaN/invalid values in masked bands to prevent NaN propagation.
# In numpy, NaN * 0 = NaN, so masked bands with NaN flux or error would
# corrupt the entire computation. Replace with safe placeholders.
flux = np.where(mask > 0, flux, 0.0)
err = np.where(mask > 0, err, 1.0)
# Apply mask to get effective dimensionality per object.
Ndim = np.sum(mask, axis=1)
# Compute variance (including errors).
var = err**2
if shared_models:
# Shared-model fast path: expand the chi-square into three matrix
# products instead of materializing (Nobj, Nmod, Nfilt) temporaries.
w = mask / var # (Nobj, Nfilt); zero in masked bands
chi2 = (
np.sum(flux**2 * w, axis=1)[:, None]
- 2.0 * (flux * w) @ mfluxes.T
+ w @ (mfluxes.T**2)
)
# Guard against tiny negative values from floating-point cancellation.
np.maximum(chi2, 0.0, out=chi2)
else:
# Mask fluxes and model fluxes.
flux_masked = flux[:, None, :] * mask[:, None, :] # (Nobj, 1, Nfilt)
mfluxes_masked = mfluxes * mask[:, None, :] # (Nobj, Nmod, Nfilt)
var_masked = var[:, None, :] * mask[:, None, :] # (Nobj, 1, Nfilt)
# Compute residuals.
resid = flux_masked - mfluxes_masked # (Nobj, Nmod, Nfilt)
# Compute chi-squared.
chi2 = np.sum(resid**2 / np.where(var_masked > 0, var_masked, np.inf), axis=2)
# Fold in any auxiliary chi-square contributions (e.g. parallax) so they
# share the same measure as the photometry under the dimensionality prior.
if extra_chi2 is not None:
extra_chi2 = np.asarray(extra_chi2)
if extra_chi2.ndim == 1:
chi2 = chi2 + extra_chi2[:, None]
else:
chi2 = chi2 + extra_chi2
if extra_dims is not None:
Ndim = Ndim + np.asarray(extra_dims)
# Compute log-likelihood.
lnl = -0.5 * chi2
# Add normalization term.
log_det_term = np.sum(np.log(var) * mask, axis=1) # (Nobj,)
lnl += -0.5 * (Ndim[:, None] * np.log(2.0 * np.pi) + log_det_term[:, None])
# Apply dimensionality prior if requested.
if dim_prior:
# Compute log-pdf of chi2 distribution in a vectorized way.
# NOTE: chi2=0 causes +inf (DOF=1) or -inf (DOF>=3) - see docstring warning
dof = Ndim - dof_reduction # effective degrees of freedom
a = 0.5 * dof[:, None] # shape (Nobj, 1)
mask_valid_dof = dof > 0 # shape (Nobj,)
# Broadcast chi2 to (Nobj, Nmod)
lnl_dim = np.full_like(lnl, -np.inf)
valid_idx = np.where(mask_valid_dof)[0]
if valid_idx.size > 0:
chi2_valid = chi2[valid_idx]
a_valid = a[valid_idx]
lnl_dim[valid_idx] = (
xlogy(a_valid - 1.0, chi2_valid)
- (chi2_valid / 2.0)
- gammaln(a_valid)
- (np.log(2.0) * a_valid)
)
lnl = lnl_dim
return lnl
[docs]
def chisquare_outlier_loglike(
flux,
err,
stellar_params=None,
parallax=None,
parallax_err=None,
p_value_cut=1e-5,
mask=None,
dof_reduction=0,
):
"""
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:
an outlier is "as likely as" an inlier model sitting at the chi-square
value corresponding to the given p-value tail, evaluated with the same
degrees of freedom the inlier likelihood uses.
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.
mask : array-like, shape (Nobj, Nfilt), optional
Binary mask (1=use, 0=skip) of the bands used by the *inlier*
likelihood. Combined with the internal finiteness check so that the
outlier model counts exactly the same bands as the inlier model it
is mixed with.
dof_reduction : int, optional
Number of degrees of freedom to subtract from the effective DOF,
matching the `dof_reduction` used by the inlier `phot_loglike`.
Default is 0.
Returns
-------
lnl_outlier : array-like, shape (Nobj,)
Log-likelihood for outlier model for each object. Objects with no
usable bands (DOF <= 0) get -inf.
"""
from scipy.stats import chi2 as chisquare
flux = np.asarray(flux)
err = np.asarray(err)
# Get effective dimensionality per object
with np.errstate(invalid="ignore"):
valid = np.isfinite(flux) & np.isfinite(err) & (err > 0)
if mask is not None:
valid = valid & (np.asarray(mask) > 0)
ndim = np.sum(valid, axis=1) # shape (Nobj,)
# Add parallax contribution to dimensionality
if parallax is not None and parallax_err is not None:
parallax = np.asarray(parallax)
parallax_err = np.asarray(parallax_err)
parallax_mask = (
np.isfinite(parallax) & np.isfinite(parallax_err) & (parallax_err > 0)
)
ndim = ndim + parallax_mask.astype(int)
dof = ndim - dof_reduction
# Compute chi-square threshold and log-probability (guarding DOF <= 0,
# where the chi-square distribution is undefined and scipy returns NaN).
dof_safe = np.where(dof > 0, dof, 1)
chi2_threshold = chisquare.ppf(1.0 - p_value_cut, dof_safe)
lnl_outlier = np.where(dof > 0, chisquare.logpdf(chi2_threshold, dof_safe), -np.inf)
return lnl_outlier
