"""
Magnitude, luptitude, and photometric log-likelihood functions.
See Also
--------
brutus.analysis.individual.BruteForce : Uses :func:`phot_loglike` for stellar fitting
brutus.analysis.populations : Uses outlier log-likelihoods for mixture fitting
brutus.core.sed_utils : SED generation functions
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
from .math import LOG_ZERO
__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 linear flux densities and uncertainties to magnitudes.
Parameters
----------
phot : `~numpy.ndarray` of shape `(Nobs, Nfilt)`
Observed photometric flux densities (in maggies when
``zeropoints=1.0``).
err : `~numpy.ndarray` of shape `(Nobs, Nfilt)`
Observed photometric flux density uncertainties in the same units as
``phot``.
zeropoints : float or `~numpy.ndarray` of shape `(Nfilt,)`, optional
Flux density zero-points. Default is ``1.0``.
Returns
-------
mag : `~numpy.ndarray` of shape `(Nobs, Nfilt)`
Magnitudes (mag) in the survey's native system.
mag_err : `~numpy.ndarray` of shape `(Nobs, Nfilt)`
First-order magnitude uncertainties (mag).
"""
mag = -2.5 * np.log10(phot / zeropoints)
mag_err = 2.5 / np.log(10.0) * err / phot
return mag, mag_err
[docs]
def inv_magnitude(mag, err, zeropoints=1.0):
"""
Convert magnitudes and uncertainties to linear flux densities.
Parameters
----------
mag : `~numpy.ndarray` of shape `(Nobs, Nfilt)`
Magnitudes (mag) in the survey's native system.
err : `~numpy.ndarray` of shape `(Nobs, Nfilt)`
Magnitude uncertainties (mag).
zeropoints : float or `~numpy.ndarray` of shape `(Nfilt,)`, optional
Flux density zero-points. Default is ``1.0``.
Returns
-------
phot : `~numpy.ndarray` of shape `(Nobs, Nfilt)`
Photometric flux densities (in maggies when ``zeropoints=1.0``).
phot_err : `~numpy.ndarray` of shape `(Nobs, Nfilt)`
First-order flux density uncertainties in the same units as ``phot``.
"""
phot = 10 ** (-0.4 * mag) * zeropoints
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 ("Luptitudes", Lupton et al. 1999).
Parameters
----------
phot : `~numpy.ndarray` with shape (Nobs, Nfilt)
Observed photometric flux densities.
err : `~numpy.ndarray` with shape (Nobs, Nfilt)
Observed photometric flux density errors.
skynoise : float or `~numpy.ndarray` with shape (Nfilt,)
Background sky noise. Used as a "softening parameter".
Default is `1.`.
zeropoints : float or `~numpy.ndarray` with shape (Nfilt,)
Flux density zero-points. Used as a "location parameter".
Default is `1.`.
Returns
-------
lupt : `~numpy.ndarray` with shape (Nobs, Nfilt)
Luptitudes corresponding to input `phot`.
lupt_err : `~numpy.ndarray` with shape (Nobs, Nfilt)
Luptitudes errors corresponding to input `err`.
"""
# 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, `LOG_ZERO`
likelihood for DOF≥3. Ensure test data has small but non-zero
residuals.
dof_reduction : int, optional
Number of degrees of freedom to subtract from the effective DOF
when using dim_prior=True. This accounts for parameters being
fitted simultaneously (e.g., scale factors, extinction).
Default is 0.
extra_chi2 : `~numpy.ndarray` with shape (Nobj,) or (Nobj, Nmod), optional
Additional chi-square contributions (e.g., from a parallax
measurement) added to the photometric chi-square *before* any
dimensionality prior is applied. This keeps auxiliary Gaussian
constraints in the same chi-square "measure" as the photometry when
`dim_prior=True` (mirroring how `BruteForce` counts parallax as an
extra band). With `dim_prior=False` the Gaussian normalization of
the extra term is not added in full (see `extra_dims`); callers
needing a proper density should add it themselves.
extra_dims : `~numpy.ndarray` with shape (Nobj,), optional
Number of extra data points behind `extra_chi2`, added to the
degrees of freedom of the dimensionality prior. Use it with
`dim_prior=True`. With `dim_prior=False` it adds only the
``ln(2 pi)`` part of the normalization of the extra term, not the
``ln(var)`` part.
Returns
-------
lnl : `~numpy.ndarray` with shape (Nobj, Nmod)
Log-likelihood values. With `dim_prior=True`, an object with no
degrees of freedom (``Ndim - dof_reduction <= 0``) gets 0 for every
model, so it adds nothing to a sum over objects.
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:
w = mask / var # (Nobj, Nfilt); zero in masked bands
chi2 = np.sum((flux[:, None, :] - mfluxes) ** 2 * w[:, None, :], 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 gives +inf (DOF=1) or LOG_ZERO (DOF>=3); see the docstring
dof = Ndim - dof_reduction # effective degrees of freedom
a = 0.5 * dof[:, None] # shape (Nobj, 1)
mask_valid_dof = dof > 0 # shape (Nobj,)
# An object with no degrees of freedom carries no information: 0.
lnl_dim = np.zeros_like(lnl)
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] = np.maximum(
xlogy(a_valid - 1.0, chi2_valid)
- (chi2_valid / 2.0)
- gammaln(a_valid)
- (np.log(2.0) * a_valid),
LOG_ZERO,
)
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
degrees of freedom (DOF <= 0) get 0, as in `phot_loglike`.
"""
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), 0.0)
return lnl_outlier