diff --git a/pyproject.toml b/pyproject.toml index e7c558ee7..46a51a54a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -63,6 +63,9 @@ plots = [ oem = [ "oem>=0.4.5", ] +mpc = [ + "mpcq>=0.4.0; python_version < '3.13'", +] docs = [ "sphinx>=8.2.0,<9.0.0", "furo>=2024.8.6", diff --git a/src/adam_core/photometry/__init__.py b/src/adam_core/photometry/__init__.py index 490230539..825a286f2 100644 --- a/src/adam_core/photometry/__init__.py +++ b/src/adam_core/photometry/__init__.py @@ -4,6 +4,7 @@ estimate_absolute_magnitude_v_from_detections, estimate_absolute_magnitude_v_from_detections_grouped, ) +from .color_determination import ColorFit, estimate_colors from .lightcurve import reduced_magnitude from .magnitude import ( calculate_apparent_magnitude_v, @@ -47,4 +48,7 @@ "RotationPeriodObservations", "RotationPeriodResult", "GroupedRotationPeriodResults", + # Color determination + "estimate_colors", + "ColorFit", ] diff --git a/src/adam_core/photometry/color_determination.py b/src/adam_core/photometry/color_determination.py new file mode 100644 index 000000000..6c08bd9a7 --- /dev/null +++ b/src/adam_core/photometry/color_determination.py @@ -0,0 +1,622 @@ +# Obtaining different colors for asteroids. + +from __future__ import annotations + +import logging +from typing import TYPE_CHECKING, Literal + +import numpy as np +import pyarrow.compute as pc +import quivr as qv +import scipy.optimize + +from ..dynamics.propagation import propagate_2body +from ..observers.observers import Observers +from .bandpasses.api import bandpass_delta_mag, map_to_canonical_filter_bands +from .hg12star import hg12star_correction +from .lightcurve import reduced_magnitude +from .magnitude import observing_geometry +from .magnitude_common import hg_phase_correction + +if TYPE_CHECKING: + # mpcq is an optional dependency (install via `adam_core[mpc]`); it is only + # referenced in type annotations here, which `from __future__ import + # annotations` keeps from being evaluated at runtime. This keeps + # `import adam_core.photometry` working without mpcq installed. + from mpcq import MPCObservations + from mpcq.orbits import MPCOrbits + +logger = logging.getLogger(__name__) + +# Color channels we fit an absolute magnitude for. These are the base band +# letters of the canonical vendored filter IDs (e.g. SDSS_g, LSST_g, PS1_g all +# reduce to the "g" channel); see `_resolve_channels`. +_BANDS = ("g", "i", "r", "u") +_PHI_TYPES = ("HG12star", "HG", "c1c2") +# If fewer than this fraction of an object's observations survive validity +# filtering, band-recognition filtering, and outlier rejection, the fit is +# not trustworthy enough to report silently. +_MIN_RETAINED_FRACTION = 0.5 +# Physically meaningful range of the H-G / HG12* slope parameter. Both G +# (Bowell et al. 1989) and G12* (Penttilä 2016) are defined on [0, 1]. +_G_BOUNDS = (0.0, 1.0) +_G_LABELS = {"HG12star": "G12*", "HG": "G"} + + +def _validate_g_bounds( + G: float, + phi_type: str, + obj_id: str, + force_g_bounds: bool, +) -> None: + """ + Check that a fitted slope parameter lies within its physical range. + + "c1c2" has no slope parameter (``G`` is NaN) and is skipped. When ``G`` is + out of range: raise ``ValueError`` if ``force_g_bounds`` is True, otherwise + log a warning and keep the fit. + """ + if phi_type == "c1c2" or not np.isfinite(G): + return + lo, hi = _G_BOUNDS + if lo <= G <= hi: + return + label = _G_LABELS[phi_type] + msg = ( + f"Fitted {label} = {G:.4f} for {obj_id} is outside the physical " + f"[{lo:g}, {hi:g}] range" + ) + if force_g_bounds: + raise ValueError(msg) + logger.warning("%s; keeping it because force_g_bounds=False", msg) + + +class ColorFit(qv.Table): + object_id = qv.LargeStringColumn() + g_mag = qv.Float64Column(nullable=True) + i_mag = qv.Float64Column(nullable=True) + r_mag = qv.Float64Column(nullable=True) + u_mag = qv.Float64Column(nullable=True) + # 1-sigma formal uncertainties on the per-band absolute magnitudes, rescaled + # to the observed scatter (see `_fit_per_band_h`). NaN for unobserved bands. + g_mag_sigma = qv.Float64Column(nullable=True) + i_mag_sigma = qv.Float64Column(nullable=True) + r_mag_sigma = qv.Float64Column(nullable=True) + u_mag_sigma = qv.Float64Column(nullable=True) + g_r = qv.Float64Column(nullable=True) + g_i = qv.Float64Column(nullable=True) + r_i = qv.Float64Column(nullable=True) + # Color uncertainties, propagated from the full parameter covariance (so the + # H_x/H_y correlation through the shared phase parameter is accounted for). + g_r_sigma = qv.Float64Column(nullable=True) + g_i_sigma = qv.Float64Column(nullable=True) + r_i_sigma = qv.Float64Column(nullable=True) + # Fitted phase slope parameter (G for "HG", G12* for "HG12star"; NaN for + # "c1c2") and its 1-sigma uncertainty. + phase_param = qv.Float64Column(nullable=True) + phase_param_sigma = qv.Float64Column(nullable=True) + # Fit-quality diagnostics over the finally-included observations. + chi2 = qv.Float64Column(nullable=True) + reduced_chi2 = qv.Float64Column(nullable=True) + dof = qv.Int64Column(nullable=True) + rank = qv.Int64Column(nullable=True) + converged = qv.BooleanColumn(nullable=True) + num_obs = qv.Int64Column(nullable=True) + num_outliers = qv.Int64Column(nullable=True) + + +def _resolve_channels( + stn: np.ndarray, bands: np.ndarray +) -> tuple[np.ndarray, np.ndarray]: + """ + Map raw (observatory_code, reported_band) pairs to g/i/r/u color channels. + + Rather than matching MPC band strings literally, this routes each observation + through the shared `map_to_canonical_filter_bands` utility, which resolves + `(observatory_code, band)` to a canonical vendored filter ID (handling MPC/ADES + label quirks, e.g. G96 "G" -> SDSS_g, ATLAS "o" -> ATLAS_o, LSST encodings, etc). + The canonical filter's base band letter is then taken as the color channel, so + every g-like filter (SDSS_g, LSST_g, PS1_g, DECam_g, so on) contributes to the "g" + fit, and filters outside the g/i/r/u set (V, PS1_w, SkyMapper_v) or rows with + no resolvable filter are returned as ``None`` (excluded downstream). + + Returns ``(channels, filter_ids)``, each an object array of length + ``len(bands)``. ``channels`` entries are one of ``"g"``, ``"i"``, ``"r"``, + ``"u"`` or ``None``; ``filter_ids`` holds the canonical vendored filter ID (or + ``None``) and is kept so callers can apply inter-system color-term corrections + (see `_apply_color_terms`). + """ + # on_unknown="skip" leaves unresolvable rows as None instead of raising, so + # unfiltered reports and unmapped bands are simply dropped from the fit. + filter_ids = map_to_canonical_filter_bands(stn, bands, on_unknown="skip") + channels = np.empty(len(filter_ids), dtype=object) + for i, fid in enumerate(filter_ids.tolist()): + if fid is None: + channels[i] = None + continue + base = str(fid).rsplit("_", 1)[-1].lower() + channels[i] = base if base in _BANDS else None + + unresolved = channels == None # noqa: E711 + if np.any(unresolved): + dropped = sorted( + {f"{s}|{b}" for s, b in zip(stn[unresolved], bands[unresolved])} + ) + logger.warning( + "Excluding %d observation(s) whose (station, band) does not resolve to a " + "g/i/r/u color channel: %s", + int(np.sum(unresolved)), + dropped, + ) + return channels, filter_ids + + +def _apply_color_terms( + m_red: np.ndarray, + filter_ids: np.ndarray, + channels: np.ndarray, + composition: str | tuple[float, float], +) -> np.ndarray: + """ + Reconcile reduced magnitudes onto a single reference filter per color channel. + + A channel may pool observations taken through different but same-letter filters + (e.g. SDSS_g and LSST_g both feed the "g" channel). Merging them directly biases + the per-channel H by the inter-system color term. This converts every row onto + the channel's reference filter using `bandpass_delta_mag`: + + m_red_ref = m_red + Δm(composition, filter_id -> reference_filter) + + Because the reference is the dominant filter, rows already in it are unchanged, + and a channel containing a single filter system is a no-op. ``composition`` (a + template id "C"/"S"/"NEO"/"MBA" or a ``(weight_C, weight_S)`` tuple) therefore + only influences channels that genuinely mix filter systems. + """ + out = np.asarray(m_red, dtype=np.float64).copy() + fid_str = np.array([str(f) for f in filter_ids.tolist()], dtype=object) + for ch in _BANDS: + in_ch = channels == ch + if not np.any(in_ch): + continue + uniq, counts = np.unique(fid_str[in_ch], return_counts=True) + if len(uniq) < 2: + continue # single filter system in this channel: nothing to reconcile + ref = str(uniq[int(np.argmax(counts))]) + for src in uniq.tolist(): + if src == ref: + continue + delta = bandpass_delta_mag(composition, src, ref) + out[in_ch & (fid_str == src)] += delta + logger.debug( + f"Color-term correction {src} -> {ref} ({ch} channel): {delta:.4f} mag" + ) + return out + + +def _prepare_geometry( + obs: MPCObservations, + object_coords, +) -> tuple[ + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, +]: + """ + Extract geometry and photometry arrays needed for per-band H fitting. + + Returns (mag, rmsmag, channels, filter_ids, r, delta, alpha_deg, valid_mask). + ``channels`` holds the resolved g/i/r/u color channel (or ``None``) for each + row and ``filter_ids`` the canonical vendored filter ID (or ``None``); see + `_resolve_channels`. valid_mask selects rows with finite mag, finite positive + rmsmag. + """ + stn = np.asarray(obs.stn.to_numpy(zero_copy_only=False), dtype=object).astype(str) + observers = Observers.from_codes(stn, obs.obstime) + + mag = obs.mag.to_numpy(zero_copy_only=False).astype(np.float64) + rmsmag = obs.rmsmag.to_numpy(zero_copy_only=False).astype(np.float64) + bands = np.asarray(obs.band.to_numpy(zero_copy_only=False), dtype=object).astype( + str + ) + channels, filter_ids = _resolve_channels(stn, bands) + + r, delta, alpha_deg = observing_geometry(object_coords, observers) + valid = np.isfinite(mag) & np.isfinite(rmsmag) & (rmsmag > 0) + return mag, rmsmag, channels, filter_ids, r, delta, alpha_deg, valid + + +def _band_selector_matrix(channels: np.ndarray) -> np.ndarray: + """N×4 selector matrix for H_g, H_i, H_r, H_u columns.""" + return np.column_stack([(channels == b).astype(float) for b in _BANDS]) + + +def _fit_per_band_h( + m_red: np.ndarray, + alpha_deg: np.ndarray, + channels: np.ndarray, + root_weights: np.ndarray, + phi_type: Literal["HG12star", "HG", "c1c2"], +) -> dict[str, float]: + """ + Fit per-band absolute magnitudes (H_g, H_i, H_r, H_u) with g(t)=0 (no rotation + term), using one of three phase-function models: + + - "HG12star": Penttilä (2016) HG12* phase function; G12* fit jointly (nonlinear). + - "HG": standard Bowell et al. H-G phase function; G fit jointly (nonlinear). + - "c1c2": polynomial phase correction c1*alpha + c2*alpha^2 (alpha in radians); + purely linear. + + ``channels`` is the resolved g/i/r/u color channel for each row (or ``None``); + see `_resolve_channels`. Observations whose channel is not one of "g", "i", "r", + "u" are excluded up front and counted as outliers. A channel with zero surviving + observations is reported as NaN. If, after all exclusions and outlier rejection, + fewer than `_MIN_RETAINED_FRACTION` of the input rows remain, the fit is + considered unreliable and raises. + + In all cases the fit is solved with iterative 3-sigma outlier rejection. + + Returns a dict of fit results and diagnostics: + + - "H_g"/"H_i"/"H_r"/"H_u" and their "_sigma": per-band absolute magnitudes and + 1-sigma uncertainties (NaN for an unobserved band). + - "g_r_sigma"/"g_i_sigma"/"r_i_sigma": color uncertainties, propagated from the + full parameter covariance so the H_x/H_y correlation is included. + - "G"/"G_sigma": fitted slope parameter (G for "HG", G12* for "HG12star"; NaN + for "c1c2") and its uncertainty. + - "chi2"/"reduced_chi2"/"dof"/"rank": goodness-of-fit over the finally-included + rows and the design-matrix rank. + - "converged": whether the (nonlinear) optimizer reported success; always True + for the linear "c1c2" solve. + - "num_obs"/"num_outliers". + + Uncertainties come from the (J'*W*J)^-1 covariance rescaled by the reduced + chi-square, i.e. errors are matched to the observed scatter rather than trusting + the absolute rmsmag calibration. + """ + n = len(m_red) + H_sel = _band_selector_matrix(channels) + full_weights = root_weights**2 + + # Rows whose (station, band) did not resolve to a color channel are already + # logged in `_resolve_channels`; here they are simply excluded from the fit. + known_band_mask = np.isin(channels, _BANDS) + included = known_band_mask.copy() + + if phi_type == "c1c2": + alpha_rad = np.deg2rad(alpha_deg) + A = np.column_stack([alpha_rad, alpha_rad**2, H_sel]) + H_idx = 2 + else: + A = H_sel + correction_fn = ( + hg12star_correction if phi_type == "HG12star" else hg_phase_correction + ) + H_init = np.array( + [ + float(np.mean(m_red[channels == b])) if np.any(channels == b) else 0.0 + for b in _BANDS + ] + ) + params0 = np.concatenate([[0.15], H_init]) + H_idx = 1 + # Loop-invariant: A, root_weights, and m_red never change across + # outlier-rejection iterations, only the `included` mask does. + Aw = A * root_weights[:, None] + Bw = m_red * root_weights + + num_params = A.shape[1] + (0 if phi_type == "c1c2" else 1) + values = np.zeros(num_params) + converged = False + while not converged: + if phi_type == "c1c2": + Aw = A[included] * root_weights[included, None] + Bw = m_red[included] * root_weights[included] + values, _, _, _ = np.linalg.lstsq(Aw, Bw, rcond=None) + res = (A @ values - m_red) ** 2 * full_weights + else: + + def func(par): + corr = correction_fn(alpha_deg, par[0]) + return ( + Aw[included] @ par[1:] + + (corr * root_weights)[included] + - Bw[included] + ) + + result = scipy.optimize.least_squares(func, params0, verbose=0) + values = result.x + corr = correction_fn(alpha_deg, values[0]) + res = (A @ values[1:] + corr - m_red) ** 2 * full_weights + params0 = values + + n_incl = int(np.sum(included)) + sigma2 = ( + np.dot(res, included) / (n_incl - num_params) + if n_incl > num_params + else np.inf + ) + outliers = res > 9 * sigma2 + new_outliers = outliers & included + converged = not np.any(new_outliers) + included &= ~outliers + + num_outliers = int(np.sum(~included)) + if n - num_outliers < _MIN_RETAINED_FRACTION * n: + raise ValueError( + f"Outlier/band rejection removed {num_outliers}/{n} observations " + f"(more than {1 - _MIN_RETAINED_FRACTION:.0%} of the data); fit is unreliable." + ) + + # Fit diagnostics: goodness of fit and parameter covariance. + # + # chi2 is the weighted sum of squared residuals over the finally-included + # rows; dof = n_incl - num_params. `J` is the weighted design matrix (c1c2) + # or the optimizer's residual Jacobian (nonlinear), both equal + # d(weighted residual)/d(params), so cov = (J'J)^-1, rescaled by the reduced + # chi-square to match the observed scatter. pinv keeps this well-defined when + # an unobserved band leaves its H column at zero (rank-deficient normal + # matrix); those bands are then masked out to NaN below. + n_incl = int(np.sum(included)) + dof = n_incl - num_params + chi2 = float(np.dot(res, included)) + reduced_chi2 = chi2 / dof if dof > 0 else float("nan") + + if phi_type == "c1c2": + J = A[included] * root_weights[included, None] + optimizer_converged = True + else: + J = np.asarray(result.jac, dtype=np.float64) + optimizer_converged = bool(result.success) + rank = int(np.linalg.matrix_rank(J)) if J.size else 0 + + if dof > 0: + cov = np.linalg.pinv(J.T @ J) * reduced_chi2 + param_sigma = np.sqrt(np.clip(np.diag(cov), 0.0, np.inf)) + else: + cov = np.full((num_params, num_params), np.nan) + param_sigma = np.full(num_params, np.nan) + + band_present = [bool(np.any(channels[known_band_mask] == b)) for b in _BANDS] + H_values = [ + float(values[H_idx + i]) if band_present[i] else float("nan") + for i in range(len(_BANDS)) + ] + H_sigma = [ + float(param_sigma[H_idx + i]) if band_present[i] else float("nan") + for i in range(len(_BANDS)) + ] + + def _color_sigma(i: int, j: int) -> float: + if not (band_present[i] and band_present[j]): + return float("nan") + a, b = H_idx + i, H_idx + j + var = float(cov[a, a] + cov[b, b] - 2.0 * cov[a, b]) + return float(np.sqrt(var)) if var > 0 else float("nan") + + # _BANDS order is (g, i, r, u) -> indices g=0, i=1, r=2, u=3. + g_r_sigma = _color_sigma(0, 2) + g_i_sigma = _color_sigma(0, 1) + r_i_sigma = _color_sigma(2, 1) + + G_fit = float(values[0]) if phi_type != "c1c2" else float("nan") + G_sigma = float(param_sigma[0]) if phi_type != "c1c2" else float("nan") + + return { + "H_g": H_values[0], + "H_i": H_values[1], + "H_r": H_values[2], + "H_u": H_values[3], + "H_g_sigma": H_sigma[0], + "H_i_sigma": H_sigma[1], + "H_r_sigma": H_sigma[2], + "H_u_sigma": H_sigma[3], + "g_r_sigma": g_r_sigma, + "g_i_sigma": g_i_sigma, + "r_i_sigma": r_i_sigma, + "G": G_fit, + "G_sigma": G_sigma, + "chi2": chi2, + "reduced_chi2": reduced_chi2, + "dof": dof, + "rank": rank, + "converged": optimizer_converged, + "num_obs": n, + "num_outliers": num_outliers, + } + + +def estimate_colors( + observations: MPCObservations, + orbits: MPCOrbits, + phi_type: Literal["HG12star", "HG", "c1c2"], + force_g_bounds: bool = True, + color_term_composition: str | tuple[float, float] | None = None, +) -> ColorFit: + """ + Estimate per-band absolute magnitudes and colors for each object. + + Inputs can contain data for multiple objects, multiple observers, and + multiple color bands. + + Parameters + ---------- + observations + MPC astrometric/photometric observations. Must have valid ``requested_provid``, + ``obstime``, ``mag``, ``band``, and ``stn`` columns. + orbits + MPC fitted orbits for the same objects. Used to propagate positions + to each observation epoch. + phi_type + Phase function type: "HG12star" (Penttilä 2016), "HG" (standard H-G), + or "c1c2" (polynomial). + force_g_bounds + Whether to enforce the physical [0, 1] range on the fitted slope + parameter (G for "HG", G12* for "HG12star"; ignored for "c1c2"). If + True (default), an out-of-range fit raises ``ValueError``. If False, it + is logged as a warning and the out-of-range value is kept -- some + analyses (e.g. Greenstreet et al.) only reproduce when values outside + [0, 1] are allowed. + color_term_composition + If set, reconcile observations from different filter systems within a + color channel (e.g. SDSS_g and LSST_g both feeding "g") onto the channel's + most-observed filter using `bandpass_delta_mag`, assuming this reflectance + spectrum: a template id ("C", "S", "NEO", "MBA") or a ``(weight_C, + weight_S)`` tuple. Channels observed through a single filter system are + unaffected, so this is a no-op unless a channel actually mixes systems. + If ``None`` (default), no color-term correction is applied and same-letter + filters are pooled directly (griz inter-system terms are ~0.01 mag; see + `_apply_color_terms`). + + Returns + ------- + ColorFit + One row per unique object found in both ``observations`` and ``orbits``. + """ + if phi_type not in _PHI_TYPES: + raise ValueError( + f"Unsupported phi_type {phi_type!r}; expected one of {_PHI_TYPES}" + ) + + len_before = len(observations) + observations = observations.apply_mask(pc.is_valid(observations.band)) + observations = observations.apply_mask(pc.is_valid(observations.mag)) + if len(observations) != len_before: + logger.info("Removed %d null bands", len_before - len(observations)) + unique_ids = [ + x for x in pc.unique(observations.requested_provid).to_pylist() if x is not None + ] + + rows: list[dict[str, object]] = [] + + for obj_id in unique_ids: + obs_mask = pc.equal(observations.requested_provid, obj_id) + obs = observations.apply_mask(obs_mask) + + orb_mask = pc.equal(orbits.requested_provid, obj_id) + orb = orbits.apply_mask(orb_mask) + if len(orb) == 0: + continue + if len(orb) > 1: + raise ValueError(f"Expected exactly one orbit for {obj_id}, got {len(orb)}") + + adam_orbits = orb.orbits() + propagated = propagate_2body(adam_orbits, obs.obstime) + object_coords = propagated.coordinates + + # Per-object outputs default to None (no fit produced) and are overwritten + # when a fit runs. G_fit stays NaN so `_validate_g_bounds` skips objects + # without a slope parameter (no valid data, or phi_type="c1c2"). + row: dict[str, object] = { + "object_id": obj_id, + "g_mag": None, + "i_mag": None, + "r_mag": None, + "u_mag": None, + "g_mag_sigma": None, + "i_mag_sigma": None, + "r_mag_sigma": None, + "u_mag_sigma": None, + "g_r": None, + "g_i": None, + "r_i": None, + "g_r_sigma": None, + "g_i_sigma": None, + "r_i_sigma": None, + "phase_param": None, + "phase_param_sigma": None, + "chi2": None, + "reduced_chi2": None, + "dof": None, + "rank": None, + "converged": None, + "num_obs": len(obs), + "num_outliers": None, + } + G_fit: float = float("nan") + + try: + mag, rmsmag, channels, filter_ids, r, delta, alpha_deg, valid = ( + _prepare_geometry(obs, object_coords) + ) + n_invalid = len(obs) - int(np.sum(valid)) + if np.any(valid): + m_red = reduced_magnitude(mag[valid], r[valid], delta[valid]) + if color_term_composition is not None: + m_red = _apply_color_terms( + m_red, + filter_ids[valid], + channels[valid], + color_term_composition, + ) + root_weights = 1.0 / rmsmag[valid] + fit = _fit_per_band_h( + m_red, alpha_deg[valid], channels[valid], root_weights, phi_type + ) + G_fit = fit["G"] + row.update( + g_mag=fit["H_g"], + i_mag=fit["H_i"], + r_mag=fit["H_r"], + u_mag=fit["H_u"], + g_mag_sigma=fit["H_g_sigma"], + i_mag_sigma=fit["H_i_sigma"], + r_mag_sigma=fit["H_r_sigma"], + u_mag_sigma=fit["H_u_sigma"], + g_r=fit["H_g"] - fit["H_r"], + g_i=fit["H_g"] - fit["H_i"], + r_i=fit["H_r"] - fit["H_i"], + g_r_sigma=fit["g_r_sigma"], + g_i_sigma=fit["g_i_sigma"], + r_i_sigma=fit["r_i_sigma"], + phase_param=G_fit, + phase_param_sigma=fit["G_sigma"], + chi2=fit["chi2"], + reduced_chi2=fit["reduced_chi2"], + dof=fit["dof"], + rank=fit["rank"], + converged=fit["converged"], + num_outliers=n_invalid + int(fit["num_outliers"]), + ) + else: + row["num_outliers"] = n_invalid + except Exception: + logger.exception("Problem when fitting colors for %s", obj_id) + raise + + _validate_g_bounds(G_fit, phi_type, obj_id, force_g_bounds) + rows.append(row) + + def _col(name: str) -> list[object]: + return [row[name] for row in rows] + + return ColorFit.from_kwargs( + object_id=_col("object_id"), + g_mag=_col("g_mag"), + i_mag=_col("i_mag"), + r_mag=_col("r_mag"), + u_mag=_col("u_mag"), + g_mag_sigma=_col("g_mag_sigma"), + i_mag_sigma=_col("i_mag_sigma"), + r_mag_sigma=_col("r_mag_sigma"), + u_mag_sigma=_col("u_mag_sigma"), + g_r=_col("g_r"), + g_i=_col("g_i"), + r_i=_col("r_i"), + g_r_sigma=_col("g_r_sigma"), + g_i_sigma=_col("g_i_sigma"), + r_i_sigma=_col("r_i_sigma"), + phase_param=_col("phase_param"), + phase_param_sigma=_col("phase_param_sigma"), + chi2=_col("chi2"), + reduced_chi2=_col("reduced_chi2"), + dof=_col("dof"), + rank=_col("rank"), + converged=_col("converged"), + num_obs=_col("num_obs"), + num_outliers=_col("num_outliers"), + ) diff --git a/src/adam_core/photometry/hg12star.py b/src/adam_core/photometry/hg12star.py new file mode 100644 index 000000000..5bf25929c --- /dev/null +++ b/src/adam_core/photometry/hg12star.py @@ -0,0 +1,131 @@ +import numpy as np +import numpy.typing as npt + + +# Basis functions Phi1, Phi2, Phi3 from Penttila et al. (2016) +# Hermite cubic spline (Appendix A, Eq. A.1). +# xs in degrees; derivatives ds are in d(y)/d(alpha_rad) as given in Penttila Table A.2/A.3. +def _hermite_spline( + x_deg: npt.NDArray[np.float64] | float, + xs_deg: npt.NDArray[np.float64], + ys: npt.NDArray[np.float64], + ds_per_rad: npt.NDArray[np.float64], +) -> npt.NDArray[np.float64] | float: + x = np.deg2rad(np.asarray(x_deg, dtype=float)) + xs = np.deg2rad(xs_deg) + scalar = x.ndim == 0 + x = np.atleast_1d(x) + j = np.clip(np.searchsorted(xs, x, side="right") - 1, 0, len(xs) - 2) + dx = xs[j + 1] - xs[j] + dy = ys[j + 1] - ys[j] + t = (x - xs[j]) / dx + a = ds_per_rad[j] * dx - dy + b = -ds_per_rad[j + 1] * dx + dy + result = (1 - t) * ys[j] + t * ys[j + 1] + t * (1 - t) * ((1 - t) * a + t * b) + return float(result[0]) if scalar else result + + +# Spline knots Table A.2: xi1 +_XI1_X = np.array([7.5, 30.0, 60.0, 90.0, 120.0, 150.0]) +_XI1_Y = np.array( + [7.5e-1, 3.3486016e-1, 1.3410560e-1, 5.1104756e-2, 2.1465687e-2, 3.6396989e-3] +) +_XI1_D = np.array( + [ + -1.9098593, + -5.5463432e-1, + -2.4404599e-1, + -9.4980438e-2, + -2.1411424e-2, + -9.1328612e-2, + ] +) + +# xi2 +_XI2_X = np.array([7.5, 30.0, 60.0, 90.0, 120.0, 150.0]) +_XI2_Y = np.array( + [9.25e-1, 6.2884169e-1, 3.1755495e-1, 1.2716367e-1, 2.2373903e-2, 1.6505689e-4] +) +_XI2_D = np.array( + [ + -5.7295780e-1, + -7.6705367e-1, + -4.5665789e-1, + -2.8071809e-1, + -1.1173257e-1, + -8.6573138e-8, + ] +) + +# xi3 Table A.3 +_XI3_X = np.array([0.0, 0.3, 1.0, 2.0, 4.0, 8.0, 12.0, 20.0, 30.0]) +_XI3_Y = np.array( + [ + 1.0, + 8.3381185e-1, + 5.7735424e-1, + 4.2144772e-1, + 2.3174230e-1, + 1.0348178e-1, + 6.1733473e-2, + 1.6107006e-2, + 0.0, + ] +) +_XI3_D = np.array( + [ + -1.0630097e-1, + -4.1180439e1, + -1.0366915e1, + -7.5784615, + -3.6960950, + -7.8605652e-1, + -4.6527012e-1, + -2.0459545e-1, + 0.0, + ] +) + + +def _phi1(alpha_deg: npt.NDArray[np.float64] | float) -> npt.NDArray[np.float64]: + a = np.asarray(alpha_deg, dtype=float) + lin = 1.0 - (6.0 / np.pi) * np.deg2rad(a) + spl = _hermite_spline(a, _XI1_X, _XI1_Y, _XI1_D) + return np.where(a <= 7.5, lin, spl) + + +def _phi2(alpha_deg: npt.NDArray[np.float64] | float) -> npt.NDArray[np.float64]: + a = np.asarray(alpha_deg, dtype=float) + lin = 1.0 - (9.0 / (5.0 * np.pi)) * np.deg2rad(a) + spl = _hermite_spline(a, _XI2_X, _XI2_Y, _XI2_D) + return np.where(a <= 7.5, lin, spl) + + +def _phi3(alpha_deg: npt.NDArray[np.float64] | float) -> npt.NDArray[np.float64]: + a = np.asarray(alpha_deg, dtype=float) + spl = _hermite_spline(a, _XI3_X, _XI3_Y, _XI3_D) + return np.where(a <= 30.0, spl, 0.0) + + +def hg12star_correction( + alpha_deg: npt.NDArray[np.float64] | float, g12star: float +) -> npt.NDArray[np.float64]: + """Compute alpha correction using H,G12* approximation. + + Parameters: + ----------- + alpha_deg: np.ndarray + angle Sun-object-observer in degrees + g12star: float + value of G12* parameter to use for computing G1 and G2 + + Returns: + -------- + Magnitude correction for the given alphas. + """ + G1 = 0.84293649 * g12star + G2 = 0.53513350 * (1.0 - g12star) + G3 = 1.0 - G1 - G2 + combined = G1 * _phi1(alpha_deg) + G2 * _phi2(alpha_deg) + G3 * _phi3(alpha_deg) + combined = np.maximum(combined, 1e-10) + return -2.5 * np.log10(combined) diff --git a/src/adam_core/photometry/magnitude.py b/src/adam_core/photometry/magnitude.py index 2a5e3410e..cbadf3bad 100644 --- a/src/adam_core/photometry/magnitude.py +++ b/src/adam_core/photometry/magnitude.py @@ -174,6 +174,39 @@ def calculate_phase_angle( return out[:n_obj] +def observing_geometry( + object_coords: CartesianCoordinates, + observers: Observers, +) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64]]: + """ + Per-observation observing geometry: heliocentric distance, observer distance, + and solar phase angle. + + Both inputs must be heliocentric (origin = SUN) and paired row-for-row; + the phase angle is delegated to `calculate_phase_angle`, which validates the + geometry (finite, r > 0, delta > 0). + + Parameters + ---------- + object_coords + Object Cartesian coordinates in AU. + observers + Observer states, aligned with ``object_coords``. + + Returns + ------- + r_au, delta_au, phase_angle_deg + Heliocentric distance (AU), observer distance (AU), and solar phase angle + (degrees) for each paired row. + """ + object_pos = np.asarray(object_coords.r, dtype=np.float64) + observer_pos = np.asarray(observers.coordinates.r, dtype=np.float64) + r_au = np.linalg.norm(object_pos, axis=1) + delta_au = np.linalg.norm(object_pos - observer_pos, axis=1) + phase_angle_deg = calculate_phase_angle(object_coords, observers) + return r_au, delta_au, phase_angle_deg + + def convert_magnitude( magnitude: npt.NDArray[np.float64], source_filter_id: npt.NDArray[np.object_], diff --git a/src/adam_core/photometry/rotation/wrappers.py b/src/adam_core/photometry/rotation/wrappers.py index 497cf7fdf..009bbc242 100644 --- a/src/adam_core/photometry/rotation/wrappers.py +++ b/src/adam_core/photometry/rotation/wrappers.py @@ -15,7 +15,7 @@ from ...observations.exposures import Exposures from ...observers.observers import Observers from ...observers.utils import calculate_observing_night -from ..magnitude import calculate_phase_angle +from ..magnitude import observing_geometry from .core import ( GroupedRotationPeriodResults, RotationPeriodObservations, @@ -125,13 +125,8 @@ def build_rotation_period_observations_from_detections( time = object_coords_helio.time.rescale("tdb") mag = _as_float64_nan(detections.mag) mag_sigma = _as_float64_nan(detections.mag_sigma) - r_au = np.linalg.norm(np.asarray(object_coords_helio.r, dtype=np.float64), axis=1) - delta_vec = np.asarray(object_coords_helio.r, dtype=np.float64) - np.asarray( - observers_helio.coordinates.r, dtype=np.float64 - ) - delta_au = np.linalg.norm(delta_vec, axis=1) - phase_angle_deg = np.asarray( - calculate_phase_angle(object_coords_helio, observers_helio), dtype=np.float64 + r_au, delta_au, phase_angle_deg = observing_geometry( + object_coords_helio, observers_helio ) if np.any(~np.isfinite(mag)): diff --git a/src/adam_core/photometry/tests/data/color_fixture_2025_MF76.npz b/src/adam_core/photometry/tests/data/color_fixture_2025_MF76.npz new file mode 100644 index 000000000..74ce17a64 Binary files /dev/null and b/src/adam_core/photometry/tests/data/color_fixture_2025_MF76.npz differ diff --git a/src/adam_core/photometry/tests/data/color_fixture_2025_MN25.npz b/src/adam_core/photometry/tests/data/color_fixture_2025_MN25.npz new file mode 100644 index 000000000..f381efcd0 Binary files /dev/null and b/src/adam_core/photometry/tests/data/color_fixture_2025_MN25.npz differ diff --git a/src/adam_core/photometry/tests/data/color_fixture_2025_MO35.npz b/src/adam_core/photometry/tests/data/color_fixture_2025_MO35.npz new file mode 100644 index 000000000..755d28906 Binary files /dev/null and b/src/adam_core/photometry/tests/data/color_fixture_2025_MO35.npz differ diff --git a/src/adam_core/photometry/tests/data/color_fixture_2025_MS34.npz b/src/adam_core/photometry/tests/data/color_fixture_2025_MS34.npz new file mode 100644 index 000000000..49a7e077e Binary files /dev/null and b/src/adam_core/photometry/tests/data/color_fixture_2025_MS34.npz differ diff --git a/src/adam_core/photometry/tests/data/color_fixture_2025_MU8.npz b/src/adam_core/photometry/tests/data/color_fixture_2025_MU8.npz new file mode 100644 index 000000000..44f89eee8 Binary files /dev/null and b/src/adam_core/photometry/tests/data/color_fixture_2025_MU8.npz differ diff --git a/src/adam_core/photometry/tests/data/color_fixture_2025_MV71.npz b/src/adam_core/photometry/tests/data/color_fixture_2025_MV71.npz new file mode 100644 index 000000000..bf64e2a7c Binary files /dev/null and b/src/adam_core/photometry/tests/data/color_fixture_2025_MV71.npz differ diff --git a/src/adam_core/photometry/tests/data/generate_color_fixtures.py b/src/adam_core/photometry/tests/data/generate_color_fixtures.py new file mode 100644 index 000000000..e4642a63d --- /dev/null +++ b/src/adam_core/photometry/tests/data/generate_color_fixtures.py @@ -0,0 +1,179 @@ +""" +Generate color determination test fixtures. + +Run from the repo root: + pdm run python src/adam_core/photometry/tests/data/generate_color_fixtures.py + +Requires a BigQuery MPC replica connection (set MPCQ_PROJECT_ID). The generated +.npz files ARE committed, so this script only needs to be re-run when adding or +refreshing fixtures. +""" + +from __future__ import annotations + +import os +from pathlib import Path + +import numpy as np +import pyarrow.compute as pc +from mpcq import MPCObservations +from mpcq.client import BigQueryMPCClient +from mpcq.orbits import MPCOrbits + +from adam_core.dynamics.propagation import propagate_2body +from adam_core.observers.observers import Observers +from adam_core.photometry.bandpasses.api import map_to_canonical_filter_bands + +OUT_DIR = Path(__file__).parent + +# Paper values from Greenstreet et al. 2026 (ApJL 996 L33) Table 2, for the +# subset of objects with a good Fourier/LSM color match. +PAPER_COLORS = { + "2025 MF76": { + "g_r_fourier": 0.58, + "g_i_fourier": 0.75, + "r_i_fourier": 0.17, + "g_r_lsm": 0.56, + "g_i_lsm": 0.75, + "r_i_lsm": 0.19, + }, + "2025 MN25": { + "g_r_fourier": 0.41, + "g_i_fourier": 0.54, + "r_i_fourier": 0.13, + "g_r_lsm": 0.42, + "g_i_lsm": 0.51, + "r_i_lsm": 0.08, + }, + "2025 MO35": { + "g_r_fourier": 0.56, + "g_i_fourier": 0.71, + "r_i_fourier": 0.15, + "g_r_lsm": 0.57, + "g_i_lsm": 0.71, + "r_i_lsm": 0.14, + }, + "2025 MS34": { + "g_r_fourier": 0.58, + "g_i_fourier": 0.75, + "r_i_fourier": 0.17, + "g_r_lsm": 0.59, + "g_i_lsm": 0.74, + "r_i_lsm": 0.15, + }, + "2025 MU8": { + "g_r_fourier": 0.46, + "g_i_fourier": 0.58, + "r_i_fourier": 0.12, + "g_r_lsm": 0.47, + "g_i_lsm": 0.58, + "r_i_lsm": 0.12, + }, + "2025 MV71": { + "g_r_fourier": 0.43, + "g_i_fourier": 0.57, + "r_i_fourier": 0.13, + "g_r_lsm": 0.44, + "g_i_lsm": 0.57, + "r_i_lsm": 0.13, + }, +} + + +def build_fixture( + obj_id: str, + mpc_obs: MPCObservations, + mpc_orb: MPCOrbits, + out_path: Path, +) -> None: + mask_obs = pc.equal(mpc_obs.provid, obj_id) + obs = mpc_obs.apply_mask(mask_obs) + mask_orb = pc.equal(mpc_orb.provid, obj_id) + orb = mpc_orb.apply_mask(mask_orb) + if len(orb) == 0: + raise ValueError(f"No orbit found for {obj_id}") + + times = obs.obstime + stns = obs.stn.to_numpy(zero_copy_only=False).astype(str) + bands_raw = obs.band.to_numpy(zero_copy_only=False).astype(str) + + orbits = orb.orbits() + prop = propagate_2body(orbits, times) + observers = Observers.from_codes(stns, times) + filter_ids = map_to_canonical_filter_bands(stns, bands_raw, on_unknown="skip") + + obj_pos = prop.coordinates.r + obs_pos = observers.coordinates.r + + paper = PAPER_COLORS[obj_id] + + orb0 = orb[0] + epoch_days = int(orb0.epoch.days[0].as_py()) + epoch_nanos = int(orb0.epoch.nanos[0].as_py()) + + np.savez_compressed( + out_path, + object_id=np.array([obj_id], dtype=object), + station=stns, + # Orbit parameters + H_v_mpc=np.array([float(orb0.h[0].as_py())], dtype=np.float64), + G_mpc=np.array([float(orb0.g[0].as_py())], dtype=np.float64), + epoch_days=np.array([epoch_days], dtype=np.int64), + epoch_nanos=np.array([epoch_nanos], dtype=np.int64), + q=np.array([float(orb0.q[0].as_py())], dtype=np.float64), + e=np.array([float(orb0.e[0].as_py())], dtype=np.float64), + inc=np.array([float(orb0.i[0].as_py())], dtype=np.float64), + node=np.array([float(orb0.node[0].as_py())], dtype=np.float64), + argperi=np.array([float(orb0.argperi[0].as_py())], dtype=np.float64), + peri_time=np.array([float(orb0.peri_time[0].as_py())], dtype=np.float64), + # Per-observation data + obsid=np.array(obs.obsid.to_pylist(), dtype=object), + obstime_days=np.array(times.days.to_pylist(), dtype=np.int64), + obstime_nanos=np.array(times.nanos.to_pylist(), dtype=np.int64), + band=bands_raw, + filter_id=np.array( + [f if f is not None else "" for f in filter_ids.tolist()], dtype=object + ), + mag_obs=np.asarray(obs.mag.to_numpy(zero_copy_only=False), dtype=np.float64), + rmsmag=np.asarray(obs.rmsmag.to_numpy(zero_copy_only=False), dtype=np.float64), + ra=np.asarray(obs.ra.to_numpy(zero_copy_only=False), dtype=np.float64), + dec=np.asarray(obs.dec.to_numpy(zero_copy_only=False), dtype=np.float64), + object_pos=np.asarray(obj_pos, dtype=np.float64), + observer_pos=np.asarray(obs_pos, dtype=np.float64), + # Paper reference colors + paper_g_r_fourier=np.array([paper["g_r_fourier"]], dtype=np.float64), + paper_g_i_fourier=np.array([paper["g_i_fourier"]], dtype=np.float64), + paper_r_i_fourier=np.array([paper["r_i_fourier"]], dtype=np.float64), + paper_g_r_lsm=np.array([paper["g_r_lsm"]], dtype=np.float64), + paper_g_i_lsm=np.array([paper["g_i_lsm"]], dtype=np.float64), + paper_r_i_lsm=np.array([paper["r_i_lsm"]], dtype=np.float64), + ) + print( + f"Wrote {out_path.name} " + f"(n={len(obs)}, bands={dict(zip(*np.unique(bands_raw, return_counts=True)))})" + ) + + +def _query_mpc_data(object_ids: list[str]) -> tuple[MPCObservations, MPCOrbits]: + client = BigQueryMPCClient( + dataset_id="asteroid_institute_mpc_replica", + views_dataset_id="asteroid_institute___mpc_replica_views", + project=os.environ["MPCQ_PROJECT_ID"], + ) + observations = client.query_observations(object_ids) + orbits = client.query_orbits(object_ids) + return observations, orbits + + +def main() -> None: + object_ids = list(PAPER_COLORS) + mpc_obs, mpc_orb = _query_mpc_data(object_ids) + + for obj_id in object_ids: + slug = obj_id.replace(" ", "_") + out_path = OUT_DIR / f"color_fixture_{slug}.npz" + build_fixture(obj_id, mpc_obs, mpc_orb, out_path) + + +if __name__ == "__main__": + main() diff --git a/src/adam_core/photometry/tests/test_color_determination.py b/src/adam_core/photometry/tests/test_color_determination.py new file mode 100644 index 000000000..f5edab963 --- /dev/null +++ b/src/adam_core/photometry/tests/test_color_determination.py @@ -0,0 +1,509 @@ +from __future__ import annotations + +from pathlib import Path +from typing import Literal + +import numpy as np +import pyarrow as pa +import pyarrow.compute as pc +import pytest +import quivr as qv +from mpcq import MPCObservations +from mpcq.orbits import MPCOrbits + +from adam_core.time import Timestamp + +from ..bandpasses.api import bandpass_delta_mag +from ..color_determination import ColorFit, _apply_color_terms, estimate_colors + +DATA_DIR = Path(__file__).parent / "data" + +COLOR_FIXTURES: list[str] = sorted(p.name for p in DATA_DIR.glob("color_fixture_*.npz")) +if not COLOR_FIXTURES: + COLOR_FIXTURES = ["__NO_FIXTURES__"] + +# Tolerances: compare within the given margin of the paper (Fourier) values. +# Greenstreet et al. 2026 reports colors from Fourier fits that include a +# rotational period term. Our implementation sets g(t)=0 (no rotation), so +# per-band H values can be biased when multi-band observations sample +# different rotational phases. 2025 MO35 will have a separate larger tolerance. +HG12STAR_TOLERANCE = 0.06 +HG_TOLERANCE = 0.06 +C1C2_TOLERANCE = 0.06 + + +def _load_fixture_observations(fx: np.lib.npyio.NpzFile) -> MPCObservations: + n = int(fx["mag_obs"].shape[0]) + obstime = Timestamp.from_kwargs( + days=pa.array(fx["obstime_days"].tolist(), type=pa.int64()), + nanos=pa.array(fx["obstime_nanos"].tolist(), type=pa.int64()), + scale="utc", + ) + return MPCObservations.from_kwargs( + requested_provid=[str(fx["object_id"][0])] * n, + primary_designation=[None] * n, + obsid=fx["obsid"].astype(str).tolist(), + trksub=[None] * n, + provid=[str(fx["object_id"][0])] * n, + permid=[None] * n, + submission_id=[None] * n, + obssubid=[None] * n, + obstime=obstime, + ra=fx["ra"].tolist(), + dec=fx["dec"].tolist(), + rmsra=[None] * n, + rmsdec=[None] * n, + rmscorr=[None] * n, + mag=fx["mag_obs"].tolist(), + rmsmag=fx["rmsmag"].tolist(), + band=fx["band"].astype(str).tolist(), + stn=fx["station"].astype(str).tolist(), + updated_at=None, + created_at=None, + status=[None] * n, + astcat=[None] * n, + mode=[None] * n, + ) + + +def _load_fixture_orbits(fx: np.lib.npyio.NpzFile) -> MPCOrbits: + obj_id = str(fx["object_id"][0]) + epoch = Timestamp.from_kwargs( + days=pa.array([int(fx["epoch_days"][0])], type=pa.int64()), + nanos=pa.array([int(fx["epoch_nanos"][0])], type=pa.int64()), + scale="tdb", + ) + return MPCOrbits.from_kwargs( + requested_provid=[obj_id], + primary_designation=[None], + id=[None], + provid=[obj_id], + epoch=epoch, + q=fx["q"].tolist(), + e=fx["e"].tolist(), + i=fx["inc"].tolist(), + node=fx["node"].tolist(), + argperi=fx["argperi"].tolist(), + peri_time=fx["peri_time"].tolist(), + q_unc=[None], + e_unc=[None], + i_unc=[None], + node_unc=[None], + argperi_unc=[None], + peri_time_unc=[None], + a1=[None], + a2=[None], + a3=[None], + h=fx["H_v_mpc"].tolist(), + g=fx["G_mpc"].tolist(), + created_at=None, + updated_at=None, + ) + + +def _paper_colors(fx: np.lib.npyio.NpzFile) -> dict[str, float]: + return { + "g_r": float(fx["paper_g_r_fourier"][0]), + "g_i": float(fx["paper_g_i_fourier"][0]), + "r_i": float(fx["paper_r_i_fourier"][0]), + } + + +def _assert_colors_close( + result: ColorFit, object_id: str, paper: dict[str, float], tolerance: float +) -> None: + # Keep tighter tolerances for all but this one + if object_id == "2025 MO35": + tolerance = max(tolerance, 0.11) + row = result.apply_mask(pc.equal(result.object_id, object_id)) + assert len(row) == 1, f"Expected 1 result row for {object_id}, got {len(row)}" + + g_r = row.g_r[0].as_py() + g_i = row.g_i[0].as_py() + r_i = row.r_i[0].as_py() + + assert np.isfinite(g_r), f"g-r not finite for {object_id}" + assert np.isfinite(g_i), f"g-i not finite for {object_id}" + assert np.isfinite(r_i), f"r-i not finite for {object_id}" + + assert abs(g_r - paper["g_r"]) <= tolerance, ( + f"{object_id} g-r: got {g_r:.3f}, paper Fourier {paper['g_r']:.3f}, " + f"diff={g_r - paper['g_r']:+.3f} > tol={tolerance}" + ) + assert abs(g_i - paper["g_i"]) <= tolerance, ( + f"{object_id} g-i: got {g_i:.3f}, paper Fourier {paper['g_i']:.3f}, " + f"diff={g_i - paper['g_i']:+.3f} > tol={tolerance}" + ) + assert abs(r_i - paper["r_i"]) <= tolerance, ( + f"{object_id} r-i: got {r_i:.3f}, paper Fourier {paper['r_i']:.3f}, " + f"diff={r_i - paper['r_i']:+.3f} > tol={tolerance}" + ) + + +@pytest.mark.parametrize( + "phi_type,tolerance", + [("HG12star", HG12STAR_TOLERANCE), ("HG", HG_TOLERANCE), ("c1c2", C1C2_TOLERANCE)], +) +@pytest.mark.parametrize("fixture_name", COLOR_FIXTURES) +def test_estimate_colors_from_fixture( + fixture_name: str, phi_type: Literal["HG", "c1c2"], tolerance: float +) -> None: + if fixture_name == "__NO_FIXTURES__": + pytest.skip("No color fixtures found on disk.") + + fixture_path = DATA_DIR / fixture_name + if not fixture_path.exists(): + pytest.skip(f"Missing fixture {fixture_name}") + + fx = np.load(fixture_path, allow_pickle=True) + object_id = str(fx["object_id"][0]) + + observations = _load_fixture_observations(fx) + orbits = _load_fixture_orbits(fx) + + # Greenstreet et al. is only reproduced when slope parameters outside the + # physical [0, 1] range are allowed, so relax the bound here. + result = estimate_colors(observations, orbits, phi_type, force_g_bounds=False) + + assert isinstance(result, ColorFit) + assert len(result) >= 1 + + _assert_colors_close(result, object_id, _paper_colors(fx), tolerance) + + +_BAND_MAG_FIELD = {"g": "g_mag", "i": "i_mag", "r": "r_mag", "u": "u_mag"} + + +def _channels_present(fx: np.lib.npyio.NpzFile) -> set[str]: + """ + Color channels (g/i/r/u) present in a fixture, derived from the canonical + ``filter_id`` the fixture generator resolved via `map_to_canonical_filter_bands` + (e.g. ``SDSS_g``/``LSST_g`` -> "g"). This mirrors how `estimate_colors` groups + observations into channels, rather than matching raw MPC band strings. + """ + present: set[str] = set() + for fid in fx["filter_id"].astype(str).tolist(): + if not fid: + continue + base = fid.rsplit("_", 1)[-1].lower() + if base in _BAND_MAG_FIELD: + present.add(base) + return present + + +@pytest.mark.parametrize("fixture_name", COLOR_FIXTURES) +def test_estimate_colors_missing_band_is_nan(fixture_name: str) -> None: + """ + A band with zero recognized-band observations for an object must be + reported as NaN, not a spuriously finite value from an unconstrained fit. + """ + if fixture_name == "__NO_FIXTURES__": + pytest.skip("No color fixtures found on disk.") + + fixture_path = DATA_DIR / fixture_name + fx = np.load(fixture_path, allow_pickle=True) + object_id = str(fx["object_id"][0]) + bands_present = _channels_present(fx) + missing_bands = set(_BAND_MAG_FIELD) - bands_present + if not missing_bands: + print(f"{fixture_name} has observations in every band; nothing to check.") + # Declare this test passing instead of skipped, to avoid making people wonder + return + + observations = _load_fixture_observations(fx) + orbits = _load_fixture_orbits(fx) + result = estimate_colors(observations, orbits, "HG12star", force_g_bounds=False) + row = result.apply_mask(pc.equal(result.object_id, object_id)) + assert len(row) == 1 + + for band, field in _BAND_MAG_FIELD.items(): + value = getattr(row, field)[0].as_py() + sigma = getattr(row, f"{field}_sigma")[0].as_py() + if band in missing_bands: + assert value is not None and np.isnan( + value + ), f"{object_id} {field}: expected NaN for unobserved band {band!r}, got {value}" + assert sigma is not None and np.isnan( + sigma + ), f"{object_id} {field}_sigma: expected NaN for unobserved band {band!r}, got {sigma}" + else: + assert value is not None and np.isfinite( + value + ), f"{object_id} {field}: expected a finite value for observed band {band!r}, got {value}" + assert ( + sigma is not None and np.isfinite(sigma) and sigma > 0 + ), f"{object_id} {field}_sigma: expected a positive finite value for band {band!r}, got {sigma}" + + +def test_estimate_colors_multi_object() -> None: + """ + estimate_colors should produce identical per-object results whether + objects are passed in one at a time or batched together. + """ + fixture_paths = [ + DATA_DIR / name for name in COLOR_FIXTURES if name != "__NO_FIXTURES__" + ] + if len(fixture_paths) < 2: + pytest.skip("Need at least two color fixtures to test multi-object batching.") + + fixtures = [np.load(p, allow_pickle=True) for p in fixture_paths] + object_ids = [str(fx["object_id"][0]) for fx in fixtures] + + observations = qv.concatenate([_load_fixture_observations(fx) for fx in fixtures]) + orbits = qv.concatenate([_load_fixture_orbits(fx) for fx in fixtures]) + + result = estimate_colors(observations, orbits, "HG12star", force_g_bounds=False) + + assert isinstance(result, ColorFit) + assert len(result) == len(object_ids) + assert set(result.object_id.to_pylist()) == set(object_ids) + + for fx, object_id in zip(fixtures, object_ids): + _assert_colors_close(result, object_id, _paper_colors(fx), HG12STAR_TOLERANCE) + + +# 2025 MN25 fits G12* well below 0 with the HG12* model, so it exercises the +# out-of-range slope-parameter handling. +_OUT_OF_BOUNDS_FIXTURE = "color_fixture_2025_MN25.npz" + + +def _load_out_of_bounds_case() -> tuple[MPCObservations, MPCOrbits, str]: + fixture_path = DATA_DIR / _OUT_OF_BOUNDS_FIXTURE + if not fixture_path.exists(): + pytest.skip(f"Missing fixture {_OUT_OF_BOUNDS_FIXTURE}") + fx = np.load(fixture_path, allow_pickle=True) + return ( + _load_fixture_observations(fx), + _load_fixture_orbits(fx), + str(fx["object_id"][0]), + ) + + +def test_force_g_bounds_true_raises_on_out_of_range() -> None: + """With force_g_bounds=True (default), an out-of-[0,1] slope fit raises.""" + observations, orbits, _ = _load_out_of_bounds_case() + + with pytest.raises(ValueError, match=r"G12\*.*outside the physical"): + estimate_colors(observations, orbits, "HG12star") + + +def test_force_g_bounds_false_warns_and_returns( + caplog: pytest.LogCaptureFixture, +) -> None: + """With force_g_bounds=False, the out-of-range fit is kept and a warning logged.""" + observations, orbits, object_id = _load_out_of_bounds_case() + + with caplog.at_level("WARNING", logger="adam_core.photometry.color_determination"): + result = estimate_colors(observations, orbits, "HG12star", force_g_bounds=False) + + assert isinstance(result, ColorFit) + row = result.apply_mask(pc.equal(result.object_id, object_id)) + assert len(row) == 1 + assert any( + "outside the physical" in record.message + and "force_g_bounds=False" in record.message + for record in caplog.records + ), "Expected an out-of-range warning mentioning force_g_bounds=False" + + +# --------------------------------------------------------------------------- +# Inter-system color-term correction (opt-in color_term_composition) +# --------------------------------------------------------------------------- + + +def test_apply_color_terms_reconciles_mixed_channel() -> None: + """ + In a channel that mixes filter systems, minority-filter rows are shifted onto + the most-observed filter by exactly `bandpass_delta_mag`, while the majority + filter and any single-system channel are left untouched. + """ + # g channel: 3x LSST_g (majority) + 1x SDSS_g (minority); r channel single-system. + filter_ids = np.array( + ["LSST_g", "LSST_g", "LSST_g", "SDSS_g", "LSST_r", "LSST_r"], dtype=object + ) + channels = np.array(["g", "g", "g", "g", "r", "r"], dtype=object) + m_red = np.array([20.0, 20.0, 20.0, 20.0, 19.0, 19.0], dtype=np.float64) + + out = _apply_color_terms(m_red, filter_ids, channels, "S") + delta = bandpass_delta_mag("S", "SDSS_g", "LSST_g") + + assert delta != 0.0 + # Majority LSST_g rows unchanged (reference filter). + np.testing.assert_array_equal(out[:3], 20.0) + # Minority SDSS_g row converted onto the LSST_g reference. + assert np.isclose(out[3], 20.0 + delta) + # Single-system r channel untouched. + np.testing.assert_array_equal(out[4:], 19.0) + # Input is not mutated. + np.testing.assert_array_equal(m_red[3], 20.0) + + +def test_apply_color_terms_single_system_is_noop() -> None: + """A channel observed through a single filter system is never corrected.""" + filter_ids = np.array(["LSST_g", "LSST_g", "LSST_r"], dtype=object) + channels = np.array(["g", "g", "r"], dtype=object) + m_red = np.array([20.0, 21.0, 19.0], dtype=np.float64) + out = _apply_color_terms(m_red, filter_ids, channels, "C") + np.testing.assert_array_equal(out, m_red) + + +def _has_mixed_channel(fx: np.lib.npyio.NpzFile) -> bool: + """True if any g/i/r/u channel draws on more than one canonical filter.""" + by_channel: dict[str, set[str]] = {} + for fid in fx["filter_id"].astype(str).tolist(): + if not fid: + continue + base = fid.rsplit("_", 1)[-1].lower() + if base in _BAND_MAG_FIELD: + by_channel.setdefault(base, set()).add(fid) + return any(len(fids) > 1 for fids in by_channel.values()) + + +@pytest.mark.parametrize("fixture_name", COLOR_FIXTURES) +def test_color_term_composition_noop_on_single_system_fixtures( + fixture_name: str, +) -> None: + """ + On a fixture whose channels are each single-system, the color-term correction + changes nothing: colors are bit-for-bit identical with and without it. + """ + if fixture_name == "__NO_FIXTURES__": + pytest.skip("No color fixtures found on disk.") + fx = np.load(DATA_DIR / fixture_name, allow_pickle=True) + if _has_mixed_channel(fx): + pytest.skip(f"{fixture_name} mixes filter systems; covered elsewhere.") + + observations = _load_fixture_observations(fx) + orbits = _load_fixture_orbits(fx) + + base = estimate_colors(observations, orbits, "HG12star", force_g_bounds=False) + corrected = estimate_colors( + observations, + orbits, + "HG12star", + force_g_bounds=False, + color_term_composition="S", + ) + for field in ("g_r", "g_i", "r_i", "g_mag", "r_mag", "i_mag"): + assert getattr(base, field).to_pylist() == getattr(corrected, field).to_pylist() + + +@pytest.mark.parametrize("fixture_name", COLOR_FIXTURES) +def test_color_term_composition_preserves_paper_on_mixed_fixtures( + fixture_name: str, +) -> None: + """ + On a fixture that mixes filter systems within a channel, applying the + correction keeps colors finite and still within tolerance of the paper (the + griz inter-system terms are small, so reproduction is preserved). + """ + if fixture_name == "__NO_FIXTURES__": + pytest.skip("No color fixtures found on disk.") + fx = np.load(DATA_DIR / fixture_name, allow_pickle=True) + if not _has_mixed_channel(fx): + pytest.skip(f"{fixture_name} has no mixed-system channel.") + + object_id = str(fx["object_id"][0]) + observations = _load_fixture_observations(fx) + orbits = _load_fixture_orbits(fx) + + result = estimate_colors( + observations, + orbits, + "HG12star", + force_g_bounds=False, + color_term_composition="S", + ) + _assert_colors_close(result, object_id, _paper_colors(fx), HG12STAR_TOLERANCE) + + +# --------------------------------------------------------------------------- +# Fit diagnostics (uncertainties, chi-square, DOF/rank, convergence) +# --------------------------------------------------------------------------- + +_DIAGNOSTICS_FIXTURE = "color_fixture_2025_MF76.npz" + + +def _load_diagnostics_fixture() -> np.lib.npyio.NpzFile: + path = DATA_DIR / _DIAGNOSTICS_FIXTURE + if not path.exists(): + pytest.skip(f"Missing fixture {_DIAGNOSTICS_FIXTURE}") + fx: np.lib.npyio.NpzFile = np.load(path, allow_pickle=True) + return fx + + +@pytest.mark.parametrize("phi_type", ["HG12star", "HG", "c1c2"]) +def test_fit_diagnostics_are_populated( + phi_type: Literal["HG12star", "HG", "c1c2"], +) -> None: + """Every fit reports goodness-of-fit, covariance-based errors, and status.""" + fx = _load_diagnostics_fixture() + observations = _load_fixture_observations(fx) + orbits = _load_fixture_orbits(fx) + + row = estimate_colors(observations, orbits, phi_type, force_g_bounds=False) + assert len(row) == 1 + + chi2 = row.chi2[0].as_py() + dof = row.dof[0].as_py() + reduced_chi2 = row.reduced_chi2[0].as_py() + assert chi2 > 0 + assert dof > 0 + assert np.isfinite(reduced_chi2) and reduced_chi2 > 0 + assert np.isclose(reduced_chi2, chi2 / dof) + assert row.converged[0].as_py() is True + + # DOF invariant: included observations minus the number of fitted parameters. + num_params = 6 if phi_type == "c1c2" else 5 + assert dof == row.num_obs[0].as_py() - row.num_outliers[0].as_py() - num_params + + # Design-matrix rank = one column per observed band, plus the phase columns + # (G for HG/HG12star; c1*alpha + c2*alpha^2 for c1c2). + present = _channels_present(fx) + phase_cols = 2 if phi_type == "c1c2" else 1 + assert row.rank[0].as_py() == len(present) + phase_cols + + # Phase slope parameter: fitted (with an uncertainty) for HG/HG12star, and + # NaN for c1c2 which has no such parameter. + phase_param = row.phase_param[0].as_py() + phase_param_sigma = row.phase_param_sigma[0].as_py() + if phi_type == "c1c2": + assert np.isnan(phase_param) and np.isnan(phase_param_sigma) + else: + assert np.isfinite(phase_param) + assert np.isfinite(phase_param_sigma) and phase_param_sigma > 0 + + # Per-band uncertainties are positive-finite for observed bands (NaN handling + # for unobserved bands is covered by test_estimate_colors_missing_band_is_nan). + for band, field in _BAND_MAG_FIELD.items(): + if band not in present: + continue + sigma = getattr(row, f"{field}_sigma")[0].as_py() + assert np.isfinite(sigma) and sigma > 0 + + for field in ("g_r_sigma", "g_i_sigma", "r_i_sigma"): + value = getattr(row, field)[0].as_py() + assert np.isfinite(value) and value > 0 + + +def test_color_sigma_propagates_covariance_not_quadrature() -> None: + """ + Color uncertainties use the full parameter covariance, so the H_x/H_y + correlation through the shared phase parameter reduces them below a naive + quadrature sum. The HG model makes G and H strongly degenerate, so the effect + is pronounced there. + """ + fx = _load_diagnostics_fixture() + observations = _load_fixture_observations(fx) + orbits = _load_fixture_orbits(fx) + + row = estimate_colors(observations, orbits, "HG", force_g_bounds=False) + g_sigma = row.g_mag_sigma[0].as_py() + r_sigma = row.r_mag_sigma[0].as_py() + g_r_sigma = row.g_r_sigma[0].as_py() + + quadrature = float(np.hypot(g_sigma, r_sigma)) + assert g_r_sigma < quadrature + # The per-band magnitudes share the G degeneracy, so each is far more + # uncertain than the color itself. + assert g_r_sigma < g_sigma diff --git a/src/adam_core/photometry/tests/test_color_fit_synthetic.py b/src/adam_core/photometry/tests/test_color_fit_synthetic.py new file mode 100644 index 000000000..984b4be3b --- /dev/null +++ b/src/adam_core/photometry/tests/test_color_fit_synthetic.py @@ -0,0 +1,144 @@ +""" +Synthetic, ground-truth unit tests for the per-band color fit. + +The fixture tests exercise the full pipeline against real MPC data and paper +values, but cannot pin down exact behaviour. Here we build reduced magnitudes +directly from known per-band absolute magnitudes and a known phase function, then +check `_fit_per_band_h` recovers the injected colors, phase parameter, outlier +count, missing-band handling, and error scaling. +""" + +from __future__ import annotations + +from typing import Literal + +import numpy as np +import pytest + +from ..color_determination import _fit_per_band_h +from ..hg12star import hg12star_correction +from ..magnitude_common import hg_phase_correction + +PhiType = Literal["HG12star", "HG", "c1c2"] + +# Injected truth: g-r = 0.6, g-i = 0.8, r-i = 0.2. +_H_TRUE = {"g": 18.0, "r": 17.4, "i": 17.2} + + +def _synthesize( + phi_type: PhiType, + phase_param: float | tuple[float, float], + H_true: dict[str, float] = _H_TRUE, + n_per_band: int = 60, + noise: float = 0.0, + seed: int = 0, +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """ + Build (m_red, alpha_deg, channels, root_weights) for a known model. + + For "HG12star"/"HG", ``phase_param`` is the scalar slope (G12*/G); for "c1c2" + it is a ``(c1, c2)`` pair with alpha in radians. + """ + rng = np.random.default_rng(seed) + bands = list(H_true) + channels = np.array([b for b in bands for _ in range(n_per_band)], dtype=object) + alpha = rng.uniform(1.0, 40.0, size=len(channels)) + base = np.array([H_true[c] for c in channels], dtype=np.float64) + + if phi_type == "c1c2": + assert isinstance(phase_param, tuple) + c1, c2 = phase_param + alpha_rad = np.deg2rad(alpha) + m_red = base + c1 * alpha_rad + c2 * alpha_rad**2 + else: + assert not isinstance(phase_param, tuple) + correction = ( + hg12star_correction(alpha, phase_param) + if phi_type == "HG12star" + else hg_phase_correction(alpha, phase_param) + ) + m_red = base + np.asarray(correction) + + sigma = noise if noise > 0 else 1.0 + if noise > 0: + m_red = m_red + rng.normal(0.0, noise, size=len(m_red)) + root_weights = np.full(len(m_red), 1.0 / sigma) + return m_red, alpha, channels, root_weights + + +@pytest.mark.parametrize( + "phi_type, phase_param", + [("HG12star", 0.4), ("HG", 0.15)], +) +def test_fit_recovers_known_colors_and_phase( + phi_type: PhiType, phase_param: float +) -> None: + """With noiseless data the fit recovers the injected colors and slope exactly.""" + m_red, alpha, channels, rw = _synthesize(phi_type, phase_param) + fit = _fit_per_band_h(m_red, alpha, channels, rw, phi_type) + + assert fit["H_g"] - fit["H_r"] == pytest.approx(0.6, abs=1e-4) + assert fit["H_g"] - fit["H_i"] == pytest.approx(0.8, abs=1e-4) + assert fit["H_r"] - fit["H_i"] == pytest.approx(0.2, abs=1e-4) + assert fit["G"] == pytest.approx(phase_param, abs=1e-4) + assert fit["converged"] is True + assert fit["num_outliers"] == 0 + + +def test_fit_recovers_known_colors_c1c2() -> None: + """The linear c1c2 model recovers colors to machine precision and G is NaN.""" + m_red, alpha, channels, rw = _synthesize("c1c2", (0.03, -5e-4)) + fit = _fit_per_band_h(m_red, alpha, channels, rw, "c1c2") + + assert fit["H_g"] - fit["H_r"] == pytest.approx(0.6, abs=1e-6) + assert fit["H_g"] - fit["H_i"] == pytest.approx(0.8, abs=1e-6) + assert fit["H_r"] - fit["H_i"] == pytest.approx(0.2, abs=1e-6) + assert np.isnan(fit["G"]) + + +def test_fit_rejects_injected_outlier() -> None: + """A single gross outlier is flagged and does not corrupt the recovered color.""" + m_red, alpha, channels, rw = _synthesize("HG12star", 0.4) + m_red = m_red.copy() + m_red[0] += 2.0 # 2-magnitude blunder on a g-band point + + fit = _fit_per_band_h(m_red, alpha, channels, rw, "HG12star") + assert fit["num_outliers"] >= 1 + assert fit["H_g"] - fit["H_r"] == pytest.approx(0.6, abs=1e-3) + + +def test_fit_reports_nan_for_absent_band() -> None: + """A band with no observations yields NaN magnitude and NaN uncertainty.""" + m_red, alpha, channels, rw = _synthesize( + "HG12star", 0.4, H_true={"g": 18.0, "r": 17.4} + ) + fit = _fit_per_band_h(m_red, alpha, channels, rw, "HG12star") + + assert np.isnan(fit["H_i"]) and np.isnan(fit["H_i_sigma"]) + assert np.isnan(fit["H_u"]) and np.isnan(fit["H_u_sigma"]) + assert np.isfinite(fit["H_g"]) and np.isfinite(fit["H_r"]) + assert fit["H_g"] - fit["H_r"] == pytest.approx(0.6, abs=1e-4) + + +def test_fit_uncertainties_scale_with_injected_noise() -> None: + """ + When the weights match the true noise, the reduced chi-square is ~1 and the + reported errors scale linearly with the noise level. Using the same seed makes + the doubled-noise realization exactly twice the smaller one, so the reported + color sigma must double. + """ + + def run(sigma: float) -> dict[str, float]: + m_red, alpha, channels, rw = _synthesize( + "HG12star", 0.4, n_per_band=300, noise=sigma, seed=7 + ) + return _fit_per_band_h(m_red, alpha, channels, rw, "HG12star") + + small = run(0.03) + large = run(0.06) + + assert small["reduced_chi2"] == pytest.approx(1.0, abs=0.2) + assert large["reduced_chi2"] == pytest.approx(1.0, abs=0.2) + assert large["g_r_sigma"] == pytest.approx(2.0 * small["g_r_sigma"], rel=1e-6) + # ~sigma / sqrt(N) per band scatter, loosely (covariance with G inflates it a bit). + assert 0.0 < small["g_r_sigma"] < 0.03 diff --git a/src/adam_core/photometry/tests/test_hg12star.py b/src/adam_core/photometry/tests/test_hg12star.py new file mode 100644 index 000000000..18d5d00dc --- /dev/null +++ b/src/adam_core/photometry/tests/test_hg12star.py @@ -0,0 +1,126 @@ +""" +Direct unit tests for the HG12* phase function against Penttila et al. (2016). + +The basis functions Phi1, Phi2, Phi3 (Appendix A, Eq. A.1) are Hermite cubic +splines. The reference tables A.2/A.3 give, at each knot, the value of the basis +function and its derivative d(Phi)/d(alpha_rad). Evaluating a basis function at +a knot must return the tabulated value exactly, and a central difference at an +interior knot must reproduce the tabulated derivative, so the tables are used +as the ground truth here. +""" + +from __future__ import annotations + +from typing import Callable + +import numpy as np +import numpy.typing as npt +import pytest + +from ..hg12star import ( + _XI1_D, + _XI1_X, + _XI1_Y, + _XI2_D, + _XI2_X, + _XI2_Y, + _XI3_X, + _XI3_Y, + _phi1, + _phi2, + _phi3, + hg12star_correction, +) + + +def test_basis_functions_are_unity_at_opposition() -> None: + """All three basis functions are normalized to 1 at zero phase angle.""" + assert _phi1(0.0) == pytest.approx(1.0) + assert _phi2(0.0) == pytest.approx(1.0) + assert _phi3(0.0) == pytest.approx(1.0) + + +@pytest.mark.parametrize( + "phi, xs, ys", + [ + (_phi1, _XI1_X, _XI1_Y), + (_phi2, _XI2_X, _XI2_Y), + (_phi3, _XI3_X, _XI3_Y), + ], +) +def test_basis_functions_interpolate_reference_knots( + phi: Callable[[float], npt.NDArray[np.float64]], + xs: np.ndarray, + ys: np.ndarray, +) -> None: + """Each basis function reproduces its tabulated knot values exactly.""" + got = np.array([phi(float(x)) for x in xs]) + np.testing.assert_allclose(got, ys, atol=1e-9) + + +@pytest.mark.parametrize( + "phi, xs, ds", + [ + (_phi1, _XI1_X, _XI1_D), + (_phi2, _XI2_X, _XI2_D), + ], +) +def test_basis_function_slopes_match_reference_table( + phi: Callable[[float], npt.NDArray[np.float64]], + xs: np.ndarray, + ds: np.ndarray, +) -> None: + """ + A central difference at each interior knot reproduces the tabulated + derivative d(Phi)/d(alpha in radians), confirming the derivative tables feed + the spline correctly. + """ + h = 1e-5 # radians + for k in range(1, len(xs) - 1): + x_rad = np.deg2rad(xs[k]) + num = (phi(np.rad2deg(x_rad + h)) - phi(np.rad2deg(x_rad - h))) / (2 * h) + assert num == pytest.approx(ds[k], rel=1e-3, abs=1e-3) + + +@pytest.mark.parametrize("alpha", [0.0, 1.5, 3.0, 7.5]) +def test_phi1_phi2_follow_closed_form_below_7p5_deg(alpha: float) -> None: + """Below 7.5 deg the first two basis functions are exactly linear in alpha.""" + assert _phi1(alpha) == pytest.approx(1.0 - (6.0 / np.pi) * np.deg2rad(alpha)) + assert _phi2(alpha) == pytest.approx( + 1.0 - (9.0 / (5.0 * np.pi)) * np.deg2rad(alpha) + ) + + +@pytest.mark.parametrize("alpha", [30.0, 45.0, 90.0, 150.0]) +def test_phi3_is_zero_beyond_30_deg(alpha: float) -> None: + """Phi3 is clamped to 0 for phase angles at or beyond 30 deg.""" + assert float(np.atleast_1d(_phi3(alpha))[0]) == 0.0 + + +@pytest.mark.parametrize("g12star", [0.0, 0.2, 0.5, 0.8, 1.0]) +def test_correction_is_zero_at_opposition(g12star: float) -> None: + """At zero phase the combined phase function is 1, so the correction is 0.""" + value = float(np.atleast_1d(hg12star_correction(np.array([0.0]), g12star))[0]) + assert value == pytest.approx(0.0, abs=1e-9) + + +def test_correction_is_monotonic_in_phase() -> None: + """ + For a physical G12* the magnitude correction grows monotonically with phase + (the object dims as it moves away from opposition). + """ + alpha = np.linspace(0.0, 60.0, 61) + corr = hg12star_correction(alpha, 0.5) + assert corr[0] == pytest.approx(0.0, abs=1e-9) + assert np.all(np.diff(corr) >= -1e-9) + assert corr[-1] > 0.0 + + +def test_correction_scalar_and_vector_agree() -> None: + """Elementwise scalar calls match a single vectorized call.""" + alpha = np.array([2.0, 10.0, 25.0, 55.0]) + vector = np.asarray(hg12star_correction(alpha, 0.4)) + elementwise = np.array( + [float(np.atleast_1d(hg12star_correction(float(a), 0.4))[0]) for a in alpha] + ) + np.testing.assert_allclose(vector, elementwise, atol=1e-12)