"""Orbit determination result types and configuration.
This module mirrors the Rust wrapper's ``empyrean::ODConfig`` and
``empyrean::DetermineResult`` field-for-field, with the same nested
structure (no flattening). ``ODConfig()`` defaults are identical to
``ODConfig::default()`` on the Rust side so no surprises round-
tripping through the C ABI.
"""
import math
from collections.abc import Iterator
from dataclasses import dataclass, field
from enum import Enum
from typing import TypeAlias
import numpy as np
from empyrean.coordinates.enums import Frame, Origin
from empyrean.od.disposition import ParamDisposition
from empyrean.od.residuals import (
AcceptabilityReport,
FitSummary,
ObservationResults,
ResidualSummary,
StationBiases,
)
from empyrean.orbits.orbits import (
CartesianOrbits,
CometaryOrbits,
KeplerianOrbits,
SphericalOrbits,
)
# JSON-like value type for the nested wire dicts marshaled across the
# C ABI boundary (str / numeric / bool / None leaves, plus nested dicts
# and lists of the same).
WireValue: TypeAlias = str | int | float | bool | list["WireValue"] | dict[str, "WireValue"] | None
# ── Enums ────────────────────────────────────────────────────
class ForceModelTier(str, Enum):
"""Force-model tier for OD propagation."""
APPROXIMATE = "approximate"
BASIC = "basic"
STANDARD = "standard"
[docs]
class SolveForParams(str, Enum):
"""Parameters to solve for in differential correction.
Mirrors ``scott::od::SolveForParams``.
"""
STATE_ONLY = "state_only"
"""Solve only for the 6-element state vector."""
STATE_AND_NONGRAV = "state_and_nongrav"
"""Solve for state + (A1, A2, A3) non-grav coefficients (9 params)."""
AUTO = "auto"
"""Start with state-only, escalate to 9-param on poor fit. Tuned
via :class:`AutoEscalationPolicy`."""
EXPLICIT = "explicit"
"""An explicit per-axis solve requested via ``solve_for_flags``
(:class:`SolveFor`) — e.g. Marsden + DT, or state + AMRAT. Reported
on :attr:`DetermineResult.solve_for_used` when the fit used the
explicit flag surface rather than one of the coarse presets above."""
[docs]
class CovarianceRepresentation(str, Enum):
"""Coordinate basis the OD output covariance is reported in."""
CARTESIAN = "cartesian"
KEPLERIAN = "keplerian"
COMETARY = "cometary"
SPHERICAL = "spherical"
[docs]
class PhotometryModel(str, Enum):
"""Photometric model for the post-OD phase-function fit.
Mirrors ``empyrean::PhotometryModel``. In ``AUTO`` the fit climbs a
model ladder -- H-only -> HG12 -> HG1G2 -- admitting the richest
model the arc's phase-angle coverage and magnitude count support,
and a :class:`PhotometryResult` reports the model it actually fitted
on ``model_used`` (never ``AUTO``). An explicit value pins a
specific model. HG12 / HG1G2 follow Muinonen et al. (2010); H-only
holds the slope fixed.
"""
AUTO = "auto"
"""Auto-select up the ladder (H-only -> HG12 -> HG1G2) by data richness."""
HONLY = "honly"
"""Fit only the absolute magnitude H (fixed slope)."""
HG = "hg"
"""Two-parameter H, G."""
HG12 = "hg12"
"""Two-parameter H, G12 (Muinonen et al. 2010)."""
HG1G2 = "hg1g2"
"""Three-parameter H, G1, G2 (Muinonen et al. 2010)."""
# ── Output epoch (tagged-union dataclass) ────────────────────
[docs]
class OutputEpochMode(str, Enum):
"""How :class:`OutputEpoch` selects the fitted-orbit epoch.
Mirrors the discriminant on ``scott::od::OutputEpoch``.
"""
MID_ARC = "mid_arc"
"""Midpoint of the observation arc (default). Resolved against
the active observation set (not the full input arc) so multi-year
arcs whose mid-arc target lies in a chaotic interval keep the
integrator anchor inside the IOD opposition window."""
LAST_OBSERVATION = "last_observation"
"""Epoch of the last observation, resolved against the active set."""
IOD_EPOCH = "iod_epoch"
"""Anchor at the IOD-derived epoch — the state stays where the
initial-orbit determination produced it. Matches OrbFit's
``epoch.eq0`` and find_orb's "anchor at most recent good fit"
pattern."""
EXPLICIT = "explicit"
""":attr:`OutputEpoch.mjd_tdb` is honored."""
[docs]
@dataclass
class OutputEpoch:
"""Output epoch for the fitted orbit. Mirrors
``scott::od::OutputEpoch``.
"""
mode: OutputEpochMode = OutputEpochMode.MID_ARC
mjd_tdb: float | None = None
"""Required when ``mode == OutputEpochMode.EXPLICIT``."""
# ── Origin policy (tagged-union dataclass) ───────────────────
[docs]
class OriginPolicyMode(str, Enum):
"""How :class:`OriginPolicy` selects the central body for IOD + DC.
Mirrors the discriminant on ``scott::od::OriginPolicy``.
"""
AUTO = "auto"
"""Selects the central body (heliocentric vs Earth-centric)
automatically. Default."""
EXPLICIT = "explicit"
"""Pin IOD + DC to a specific central body — set
:attr:`OriginPolicy.origin` to the desired :class:`Origin`. Skips
the cascade."""
[docs]
@dataclass
class OriginPolicy:
"""Origin-policy selector for the OD pipeline. Mirrors
``scott::od::OriginPolicy``.
Auto handles TCOs / minimoons / geocentric impactors / chaotic-
capture interiors without per-object regime classification by the
caller. Explicit is required for cataloged satellites where
heliocentric Gauss is unphysical, and recommended for pipelines
that already know the regime.
"""
mode: OriginPolicyMode = OriginPolicyMode.AUTO
origin: "Origin | str | None" = None
"""The central body to pin to. Required when
``mode == OriginPolicyMode.EXPLICIT``. Pass an :class:`Origin`
instance or a canonical name string."""
# ── Nested config bundles ────────────────────────────────────
[docs]
@dataclass
class IODConfig:
"""IOD ranging tuning. Mirrors the IOD section of
``scott::od::ODConfig``. Defaults match ``ODConfig::default()``."""
max_triplet_attempts: int = 10
max_triplet_span_days: float = 30.0
opposition_gap_days: float = 90.0
"""Set to a negative value to disable opposition splitting."""
max_iod_arc_days: float = 30.0
"""Maximum arc length (days) used for IOD."""
curvature_snr_threshold: float = 3.0
max_iod_fractional_sigma_a: float = 1.0
[docs]
@dataclass
class AutoEscalationPolicy:
"""Trigger thresholds for :attr:`SolveForParams.AUTO` escalation.
Mirrors ``scott::od::AutoEscalationPolicy``."""
reduced_chi2: float = 10.0
at_ct_ratio: float = 3.0
min_arc_days: float = 30.0
min_n_obs: int = 50
[docs]
@dataclass
class AcceptabilityThresholds:
"""Thresholds for the post-DC fit-acceptability sub-checks. Mirrors
``scott::od::AcceptabilityThresholds``.
Defaults are tuned for production NEO survey work; tighten for
Sentry-grade impact-monitoring orbits (e.g. ``fractional_sigma_a =
1e-4``), loosen for short-arc discovery fits.
"""
reduced_chi2: float = 3.0
rms_arcsec: float = 1.0
at_ct_ratio: float = 3.0
min_arc_days: float = 7.0
fractional_sigma_a: float = 0.1
# ── Weighting (mirrors empyrean::WeightingConfig) ─────────────
[docs]
class WeightingPreset(str, Enum):
"""Preset selector for :class:`WeightingConfig`.
Picking a preset seeds the layer chain with scott's curated
layers; entries in :attr:`WeightingConfig.additional_layers` are
placed ahead of the preset's rules (sigma resolution is
first-match-wins, so they override the preset for their stations
and the preset is the fallback).
"""
NONE = "none"
"""No preset — only ``additional_layers`` apply."""
VFCC2017 = "vfcc2017"
"""Vereš, Farnocchia, Chesley & Chamberlin 2017 station floors +
nightly de-weighting at floor-σ policy. Production default."""
NEODYS = "neodys"
"""NEODyS production preset."""
[docs]
class SigmaPolicy(str, Enum):
"""How a weighting layer's σ combines with the per-observation
reported σ. Mirrors ``scott::weighting::SigmaPolicy``."""
DEFAULT_ONLY = "default_only"
"""``σ = reported`` if present, else ``σ = rule``. Default."""
FLOOR = "floor"
"""``σ = max(reported, rule)``. VFCC2017 / NEODyS production policy."""
[docs]
class WeightingLayerKind(str, Enum):
"""Discriminator for :class:`WeightingLayer` variants."""
OBSERVATORY_RULE = "observatory_rule"
NIGHTLY_DEWEIGHTING = "nightly_deweighting"
[docs]
@dataclass
class WeightingLayer:
"""One element of the weighting pipeline. Tagged-union shape —
the active fields depend on :attr:`kind`, and fields belonging to
the *other* kind must be left at their defaults: a
``NIGHTLY_DEWEIGHTING`` layer reads only :attr:`max_gap_days`
(nightly de-weighting cannot be scoped by station or time range),
and an ``OBSERVATORY_RULE`` layer never reads
:attr:`max_gap_days`. Setting an inapplicable field raises
``ValueError`` at construction rather than being silently
ignored.
Mirrors ``scott::weighting::WeightingLayer``.
"""
kind: WeightingLayerKind = WeightingLayerKind.OBSERVATORY_RULE
# ── ObservatoryRule fields ─────────────────────────────────
obs_code: str = ""
"""MPC observatory code (e.g. ``"F51"``). Matched exactly and
case-sensitively."""
sigma: tuple[float, float] = (1.0, 1.0)
"""1σ (RA·cos(δ), Dec) in arcseconds."""
start_epoch_mjd_tdb: float | None = None
"""Start of applicable time range (MJD TDB). ``None`` =
unbounded."""
end_epoch_mjd_tdb: float | None = None
"""End of applicable time range (MJD TDB). ``None`` =
unbounded."""
scale: float = 1.0
"""Scale factor on the final weight."""
# ── NightlyDeweighting fields ──────────────────────────────
max_gap_days: float = 0.5
"""Max gap (days) between observations to count as the same
night (NightlyDeweighting only)."""
def __post_init__(self) -> None:
if self.kind == WeightingLayerKind.NIGHTLY_DEWEIGHTING:
# NightlyDeweighting reads only max_gap_days; a scoping
# field here would be silently inert engine-side.
inert = []
if self.obs_code != "":
inert.append(f"obs_code={self.obs_code!r}")
if tuple(self.sigma) != (1.0, 1.0):
inert.append(f"sigma={tuple(self.sigma)!r}")
if self.start_epoch_mjd_tdb is not None:
inert.append(f"start_epoch_mjd_tdb={self.start_epoch_mjd_tdb!r}")
if self.end_epoch_mjd_tdb is not None:
inert.append(f"end_epoch_mjd_tdb={self.end_epoch_mjd_tdb!r}")
if self.scale != 1.0:
inert.append(f"scale={self.scale!r}")
if inert:
raise ValueError(
f"NIGHTLY_DEWEIGHTING layer reads only max_gap_days, but "
f"{', '.join(inert)} was set. Nightly de-weighting groups "
f"per station internally and always applies to every "
f"station; use an OBSERVATORY_RULE layer for per-station "
f"sigmas."
)
if not math.isfinite(self.max_gap_days) or self.max_gap_days <= 0.0:
raise ValueError(
f"NIGHTLY_DEWEIGHTING layer: max_gap_days must be finite "
f"and > 0 (days), got {self.max_gap_days!r}; the "
f"production default is 0.5."
)
elif self.kind == WeightingLayerKind.OBSERVATORY_RULE:
if self.max_gap_days != 0.5:
raise ValueError(
f"OBSERVATORY_RULE layer (obs_code={self.obs_code!r}) does "
f"not read max_gap_days (got {self.max_gap_days!r}); use a "
f"NIGHTLY_DEWEIGHTING layer instead."
)
# obs_code must survive the 4-byte C-ABI field byte-exactly:
# matching is exact and case-sensitive engine-side, so a
# truncated / trimmed / repaired code would match the wrong
# station, or none.
if self.obs_code == "":
raise ValueError("OBSERVATORY_RULE layer: obs_code must be a non-empty MPC code.")
if len(self.obs_code.encode("utf-8", errors="surrogatepass")) > 4:
raise ValueError(
f"OBSERVATORY_RULE layer: obs_code {self.obs_code!r} is "
f"longer than the 4-byte MPC station field; matching is "
f"exact — a truncated code would match the wrong station."
)
if not all(
ch.isascii() and ch.isprintable() and not ch.isspace() for ch in self.obs_code
):
raise ValueError(
f"OBSERVATORY_RULE layer: obs_code {self.obs_code!r} must "
f"be printable ASCII with no whitespace — MPC station "
f"matching is exact and case-sensitive."
)
# Defense-in-depth value checks (the engine-side layer-sigma
# validation ships separately): a non-finite or non-positive
# sigma / scale produces NaN or infinite weights downstream.
sigma = tuple(self.sigma)
if len(sigma) != 2:
raise ValueError(
f"OBSERVATORY_RULE layer (obs_code={self.obs_code!r}): "
f"sigma must be a (ra, dec) pair, got {sigma!r}."
)
if not all(math.isfinite(s) and s > 0.0 for s in sigma):
raise ValueError(
f"OBSERVATORY_RULE layer (obs_code={self.obs_code!r}): "
f"sigma must be finite and > 0 arcsec, got {sigma!r}."
)
if not math.isfinite(self.scale) or self.scale <= 0.0:
raise ValueError(
f"OBSERVATORY_RULE layer (obs_code={self.obs_code!r}): "
f"scale must be finite and > 0, got {self.scale!r}; use "
f"1.0 for no scaling."
)
for name, value in (
("start_epoch_mjd_tdb", self.start_epoch_mjd_tdb),
("end_epoch_mjd_tdb", self.end_epoch_mjd_tdb),
):
if value is not None and not math.isfinite(value):
raise ValueError(
f"OBSERVATORY_RULE layer (obs_code={self.obs_code!r}): "
f"{name} must be a finite MJD TDB or None for "
f"unbounded, got {value!r}."
)
def _default_weighting_layers() -> list["WeightingLayer"]:
# Mirrors `scott::od::ODConfig::default()` — the production preset is
# VFCC2017 station floors WITH NightlyDeweighting (1/√N within 0.5 days)
# appended. Without the nightly layer, OD on objects with clustered
# same-station-same-night observations diverges from `validate-core`'s
# direct-scott path because the rejection layer treats high-weight
# cluster residuals as outliers.
return [
WeightingLayer(
kind=WeightingLayerKind.NIGHTLY_DEWEIGHTING,
max_gap_days=0.5,
)
]
[docs]
@dataclass
class WeightingConfig:
"""Observation weighting pipeline. Mirrors
``empyrean::WeightingConfig``.
Default = enabled with the VFCC2017 preset + a NightlyDeweighting layer
(production hot path; matches ``scott::od::ODConfig::default()``).
Set ``enabled=False`` for uniform 1″ weighting; pick a different
preset or replace ``additional_layers`` for custom pipelines.
.. warning::
Passing ``additional_layers`` **replaces** the default layer
list — it does not append to it. The default list is
``[WeightingLayer(kind=WeightingLayerKind.NIGHTLY_DEWEIGHTING,
max_gap_days=0.5)]``, the production per-night 1/√N
de-weighting; supplying e.g. a single observatory rule
therefore drops nightly de-weighting from the production
default. To add a rule while keeping production behavior,
include the nightly layer explicitly::
WeightingConfig(
additional_layers=[
WeightingLayer(
kind=WeightingLayerKind.OBSERVATORY_RULE,
obs_code="F51",
sigma=(0.1, 0.1),
),
WeightingLayer(kind=WeightingLayerKind.NIGHTLY_DEWEIGHTING),
],
)
At most one NIGHTLY_DEWEIGHTING layer is accepted per chain
(duplicates would compound the 1/√N factor and are rejected).
"""
enabled: bool = True
preset: WeightingPreset = WeightingPreset.VFCC2017
default_sigma_arcsec: float = 1.0
"""Default 1σ when no rule applies (arcsec). Used only when
``preset = NONE``."""
sigma_policy: SigmaPolicy | None = None
"""Sigma combination policy override. ``None`` = use the
preset's policy."""
additional_layers: list[WeightingLayer] = field(default_factory=_default_weighting_layers)
"""Layers placed ahead of the preset's chain — first-match-wins,
so they override preset rules for their stations; the preset is
the fallback."""
# ── Debiasing (mirrors empyrean::DebiasingConfig) ─────────────
[docs]
class DebiasingResolution(str, Enum):
"""Healpix resolution of a debiasing table."""
STANDARD = "standard"
"""NSIDE = 64, ~35 MB. Production default."""
HIRES = "hires"
"""NSIDE = 256, ~567 MB."""
[docs]
@dataclass
class DebiasingConfig:
"""Catalog-bias-correction configuration. Mirrors scott's
``Option<Arc<DebiasingTable>>`` field on ``ODConfig``.
Default = enabled at standard resolution with no explicit path
(uses the DataManager default lookup at
``~/.empyrean/data/bias.dat``). Set ``enabled=False`` to disable
catalog debiasing entirely.
"""
enabled: bool = True
resolution: DebiasingResolution = DebiasingResolution.STANDARD
bias_dat_path: str | None = None
[docs]
class RejectionKind(str, Enum):
"""Outlier-rejection strategy selector.
Pick the variant that matches the reference pipeline you're
interoperating with — `ADAPTIVE` is the production default
(information-loss-weighted, Layer 3); `CMC2003` matches the
OrbFit / NEODyS χ²-with-hysteresis scheme of Carpino, Milani
& Chesley (2003).
"""
ADAPTIVE = "adaptive"
CMC2003 = "cmc2003"
[docs]
@dataclass
class RejectionConfig:
"""Outlier-rejection configuration. The active fields are
determined by :attr:`kind`.
Mirrors ``scott::rejection::RejectionStrategy`` plus the upstream
``max_rejection_passes`` knob. Set ``enabled=False`` to disable
the rejection pass entirely.
"""
enabled: bool = True
kind: RejectionKind = RejectionKind.ADAPTIVE
"""Strategy selector. Default :attr:`RejectionKind.ADAPTIVE`."""
# ── Adaptive (kind = ADAPTIVE) ──────────────────────────
chi2_base: float = 9.21
"""χ²(2 dof, p = 0.01) — Carpino, Milani & Chesley 2003.
Adaptive rejection only."""
lambda_: float = 1.0
"""Adaptation strength. ``0`` reduces to standard χ² rejection;
higher values protect informative observations more.
Adaptive rejection only."""
max_threshold: float = 100.0
"""Effective-threshold cap for adaptive rejection."""
# ── CMC2003 (kind = CMC2003) ────────────────────────────
chi2_rej: float = 8.0
"""χ²-with-hysteresis upper threshold — reject when χ² > chi2_rej.
CMC2003 only. Default 8.0 (≈ 98.2% confidence at 2 DOF)."""
chi2_rec: float = 7.0
"""χ²-with-hysteresis lower threshold — recover a previously-
rejected observation when χ² < chi2_rec. CMC2003 only. Must
satisfy ``chi2_rec < chi2_rej`` for hysteresis to break cycles.
Default 7.0 (≈ 96.9% confidence at 2 DOF)."""
# ── Both ────────────────────────────────────────────────
max_passes: int = 4
[docs]
@dataclass
class StationRaDecConfig:
"""Per-station RA/Dec bias-fit configuration.
Schur-eliminated nuisance parameters that absorb per-station
pointing offsets, fit alongside the orbit. Default thresholds
target modern survey arcs.
Attributes
----------
sigma_prior_arcsec : float
1-σ Gaussian prior on the per-station offset, in arcseconds.
Default 0.3.
min_obs_per_station : int
Minimum observations per station required to allocate a
bias parameter for that station. Default 5.
"""
sigma_prior_arcsec: float = 0.3
min_obs_per_station: int = 5
# ── Wide solve-for + photometry request (mirrors empyrean::SolveFor /
# empyrean::PhotometryConfig) ─────────────────────────────────
MAX_THRUST_SEGMENTS = 3
"""Largest number of thrust Δv segments one fit can declare."""
[docs]
@dataclass
class SolveFor:
"""What the fit does with each parameter axis. Mirrors ``empyrean::SolveFor``.
Set on :attr:`ODConfig.solve_for_flags` to request an explicit
multi-axis fit that the coarse :class:`SolveForParams` variants
can't name. Each axis carries a
:class:`~empyrean.od.disposition.ParamDisposition` rather than a
flag, subject to its own precondition (a declared prior on the
orbit) enforced by the engine.
Why not booleans
----------------
"Not solved" is two different answers. A **fixed** axis is
marginalized out and contributes nothing; a **considered** axis is
not estimated but its prior uncertainty still reaches the posterior
through its measurement partials, so the reported σ accounts for an
error source the fit did not absorb. Both produce well-formed
covariances, and ``False`` cannot say which was meant — so the axes
take tags, and a bool is refused rather than coerced.
"""
marsden: ParamDisposition = ParamDisposition.FIXED
"""Disposition of the Marsden A1/A2/A3 block (requires a non-grav covariance)."""
dt: ParamDisposition = ParamDisposition.FIXED
"""Disposition of the non-grav time delay DT (requires ``marsden`` solved + a DT prior)."""
amrat: ParamDisposition = ParamDisposition.FIXED
"""Disposition of the SRP AMRAT (requires an SRP AMRAT prior)."""
thrust: list[ParamDisposition] = field(default_factory=list)
"""Disposition of each **declared** thrust Δv segment, positional with
the orbit's ``correction_covariances``.
Positional rather than a count: a considered or fixed burn sits
between solved ones as readily as after them, and a count cannot say
which burn is which.
"""
def __post_init__(self) -> None:
# Parse and validate eagerly, so a bad tag fails where it was
# written rather than deep in the marshal.
self.marsden = ParamDisposition.parse(self.marsden)
self.dt = ParamDisposition.parse(self.dt)
self.amrat = ParamDisposition.parse(self.amrat)
if len(self.thrust) > MAX_THRUST_SEGMENTS:
raise ValueError(
f"thrust has {len(self.thrust)} entries but the engine's maximum is "
f"{MAX_THRUST_SEGMENTS} — the 17-column wide budget is shared across "
"the state, Marsden, DT, AMRAT and thrust axes"
)
self.thrust = [ParamDisposition.parse(d) for d in self.thrust]
[docs]
@dataclass
class PhotometryConfig:
"""Post-OD photometric-fit configuration. Mirrors
``empyrean::PhotometryConfig``.
Attach via :attr:`ODConfig.photometry`. The fit runs after the
orbit is solved and never touches the state. Sentinel rule:
``0`` / ``0.0`` on a tuning field requests the engine default.
"""
model: PhotometryModel = PhotometryModel.AUTO
"""Model to fit. Default :attr:`PhotometryModel.AUTO`."""
sigma_lightcurve: float = 0.0
"""1σ lightcurve scatter floor (mag). ``0.0`` → engine default (0.2)."""
include_rejected: bool = False
"""Include astrometrically-rejected observations' magnitudes."""
max_irls_iterations: int = 0
"""Max Huber-IRLS iterations. ``0`` → engine default (30)."""
huber_k: float = 0.0
"""Huber tuning constant. ``0.0`` → engine default (1.5)."""
# ── Top-level config ─────────────────────────────────────────
[docs]
@dataclass
class ODConfig:
"""Unified orbit-determination configuration.
Sensible production defaults out of the box:
- VFCC2017 station weighting + nightly de-weighting
(:attr:`WeightingConfig.preset`)
- EFCC2020 catalog debiasing enabled
(:attr:`DebiasingConfig.enabled`)
- :attr:`SolveForParams.AUTO` (escalates 6→9 parameters on poor fit)
- Adaptive outlier rejection enabled, ``max_passes = 4``
"""
# ── Shared (all OD entry points) ────────────────────────
force_model: ForceModelTier = ForceModelTier.STANDARD
epsilon: float = 1e-9
"""Adaptive integrator truncation-error tolerance."""
max_light_time_iterations: int = 3
num_threads: int = 0
"""``0`` = use all available cores."""
frame: Frame = Frame.ICRF
weighting: "WeightingConfig" = field(default_factory=lambda: WeightingConfig())
"""Observation weighting pipeline. Default = enabled + VFCC2017
preset. See :class:`WeightingConfig` for full layered control."""
debiasing: "DebiasingConfig" = field(default_factory=lambda: DebiasingConfig())
"""Catalog-bias-correction configuration. Default = EFCC2020
standard resolution loaded from the engine's default data location.
See :class:`DebiasingConfig`."""
excluded_perturbers: list[Origin | str] = field(default_factory=list)
"""Bodies to omit from the perturber set. Pass :class:`Origin`
instances (or canonical names). Useful when fitting an asteroid
that the force model would otherwise include as a perturber —
e.g. fitting Eros while excluding ``Origin.asteroid(433)``."""
origin: OriginPolicy = field(default_factory=OriginPolicy)
"""Origin-policy selector. Default :attr:`OriginPolicyMode.AUTO`
(heliocentric → geocentric Earth cascade). Set
``origin=OriginPolicy(mode=OriginPolicyMode.EXPLICIT, origin=Origin.EARTH)``
to pin the pipeline to a specific central body for catalog
satellites or regime-classified workflows."""
# ── IOD (determine only) ────────────────────────────────
iod: IODConfig = field(default_factory=IODConfig)
# ── Differential correction ─────────────────────────────
output_epoch: OutputEpoch = field(default_factory=OutputEpoch)
max_iterations: int = 100
convergence_tol: float = 1e-5
"""DC convergence tolerance, tested on the **undamped** Gauss–Newton
step's quadratic form — the number published back as
:attr:`DetermineResult.gn_step_qnorm`, NOT the μ-damped
:attr:`DetermineResult.update_norm`, which this tolerance does not
bound. The engine default is 1e-5."""
allow_arc_truncation: bool = True
"""Allow the outward-expansion pipeline to truncate a sub-arc it
cannot fit as one piece.
Setting ``False`` makes an arc spanning a dynamical discontinuity
FAIL loudly instead of delivering a fit of the reconcilable sub-arc
with the rest tagged ``outside_arc``. Per-observation rejection is
orthogonal and still runs, and under an ``AUTO`` origin policy the
refusal feeds the origin cascade rather than surfacing — pin the
origin to get a pure loud failure."""
coorbital_enabled: bool = True
"""Master switch for the co-orbital IOD lane. Enabling it does not
route ordinary objects through the lane — it still fires only when
every co-orbitality gate passes."""
solve_for: SolveForParams = SolveForParams.AUTO
auto_escalation: AutoEscalationPolicy = field(default_factory=AutoEscalationPolicy)
acceptability: AcceptabilityThresholds = field(default_factory=AcceptabilityThresholds)
fit_station_biases: bool = False
"""Enable Schur-eliminated per-station RA/Dec bias fitting."""
station_radec: StationRaDecConfig = field(default_factory=StationRaDecConfig)
use_span_grouping: bool = False
# ── Rejection ──────────────────────────────────────────
rejection: RejectionConfig = field(default_factory=RejectionConfig)
auto_force_model: bool = False
"""Auto-select force-model tier from IOD orbital elements."""
output_representation: CovarianceRepresentation = CovarianceRepresentation.CARTESIAN
solve_for_flags: SolveFor | None = None
"""Explicit per-axis wide solve request. When set, overrides the
coarse :attr:`solve_for` and asks the engine for an ``Explicit``
fit over the requested axes (Marsden / DT / AMRAT / thrust).
``None`` = use :attr:`solve_for`."""
allow_unbracketed_maneuvers: bool = False
"""Permit solving a thrust Δv segment whose burn window is not
bracketed by observations (the state absorbs it otherwise). Default
``False`` — refuse loudly."""
photometry: PhotometryConfig | None = None
"""Post-OD photometric fit. ``None`` (default) disables it; the fit
runs after the orbit is solved and never touches the state."""
def _to_wire_dict(self) -> dict[str, WireValue]:
"""Serialize to the nested dict shape the binding consumes.
Internal — called by :func:`empyrean.determine` /
:func:`empyrean.evaluate` / :func:`empyrean.refine` to marshal
the config across the FFI boundary. For user-facing
serialization (saving config to JSON, displaying it in a
notebook, etc.), use :func:`dataclasses.asdict`.
"""
wire: dict[str, WireValue] = {
"force_model": _enum_value(self.force_model),
"epsilon": self.epsilon,
"max_light_time_iterations": self.max_light_time_iterations,
"num_threads": self.num_threads,
"frame": _enum_value(self.frame),
"weighting": _weighting_to_dict(self.weighting),
"debiasing": _debiasing_to_dict(self.debiasing),
"excluded_perturbers_naif": [_origin_to_naif(o) for o in self.excluded_perturbers],
"origin": {
"mode": _enum_value(self.origin.mode),
"naif_id": (
_origin_to_naif(self.origin.origin) if self.origin.origin is not None else None
),
},
"iod": {
"max_triplet_attempts": self.iod.max_triplet_attempts,
"max_triplet_span_days": self.iod.max_triplet_span_days,
"opposition_gap_days": self.iod.opposition_gap_days,
"max_iod_arc_days": self.iod.max_iod_arc_days,
"curvature_snr_threshold": self.iod.curvature_snr_threshold,
"max_iod_fractional_sigma_a": self.iod.max_iod_fractional_sigma_a,
},
"output_epoch": {
"mode": self.output_epoch.mode,
"mjd_tdb": self.output_epoch.mjd_tdb,
},
"max_iterations": self.max_iterations,
"convergence_tol": self.convergence_tol,
"allow_arc_truncation": self.allow_arc_truncation,
"coorbital_enabled": self.coorbital_enabled,
"solve_for": _enum_value(self.solve_for),
"auto_escalation": {
"reduced_chi2": self.auto_escalation.reduced_chi2,
"at_ct_ratio": self.auto_escalation.at_ct_ratio,
"min_arc_days": self.auto_escalation.min_arc_days,
"min_n_obs": self.auto_escalation.min_n_obs,
},
"acceptability": {
"reduced_chi2": self.acceptability.reduced_chi2,
"rms_arcsec": self.acceptability.rms_arcsec,
"at_ct_ratio": self.acceptability.at_ct_ratio,
"min_arc_days": self.acceptability.min_arc_days,
"fractional_sigma_a": self.acceptability.fractional_sigma_a,
},
"fit_station_biases": self.fit_station_biases,
"station_radec": {
"sigma_prior_arcsec": self.station_radec.sigma_prior_arcsec,
"min_obs_per_station": self.station_radec.min_obs_per_station,
},
"use_span_grouping": self.use_span_grouping,
"rejection": {
"enabled": self.rejection.enabled,
"kind": _enum_value(self.rejection.kind),
"chi2_base": self.rejection.chi2_base,
# Python alias `lambda_` keeps the keyword from
# collising with the language; wire format uses bare
# `lambda` so it round-trips through Rust unchanged.
"lambda": self.rejection.lambda_,
"max_threshold": self.rejection.max_threshold,
"chi2_rej": self.rejection.chi2_rej,
"chi2_rec": self.rejection.chi2_rec,
"max_passes": self.rejection.max_passes,
},
"auto_force_model": self.auto_force_model,
"output_representation": _enum_value(self.output_representation),
"allow_unbracketed_maneuvers": self.allow_unbracketed_maneuvers,
}
# Explicit per-axis solve request and photometry config are only
# emitted when set — the Rust parser reads these keys only when
# present, so their absence leaves the coarse `solve_for` /
# photometry-off defaults untouched.
if self.solve_for_flags is not None:
wire["solve_for_flags"] = {
"marsden": self.solve_for_flags.marsden.value,
"dt": self.solve_for_flags.dt.value,
"amrat": self.solve_for_flags.amrat.value,
"thrust": [d.value for d in self.solve_for_flags.thrust],
}
if self.photometry is not None:
wire["photometry"] = {
"model": _enum_value(self.photometry.model),
"sigma_lightcurve": self.photometry.sigma_lightcurve,
"include_rejected": self.photometry.include_rejected,
"max_irls_iterations": self.photometry.max_irls_iterations,
"huber_k": self.photometry.huber_k,
}
return wire
def _weighting_to_dict(w: WeightingConfig) -> dict[str, WireValue]:
"""Serialize a :class:`WeightingConfig` to the wire dict the
PyO3 bridge expects."""
nightly_indices = [
i
for i, layer in enumerate(w.additional_layers)
if layer.kind == WeightingLayerKind.NIGHTLY_DEWEIGHTING
]
if len(nightly_indices) > 1:
raise ValueError(
f"WeightingConfig.additional_layers contains "
f"{len(nightly_indices)} NIGHTLY_DEWEIGHTING layers (at indices "
f"{nightly_indices}); each additional pass compounds the "
f"per-night 1/sqrt(N) de-weighting multiplicatively — include "
f"exactly one."
)
return {
"enabled": w.enabled,
"preset": _enum_value(w.preset),
"default_sigma_arcsec": w.default_sigma_arcsec,
"sigma_policy": _enum_value(w.sigma_policy) if w.sigma_policy is not None else None,
"additional_layers": [_weighting_layer_to_dict(layer) for layer in w.additional_layers],
}
def _weighting_layer_to_dict(layer: WeightingLayer) -> dict[str, WireValue]:
# Emit only the fields the layer's kind reads: the binding rejects
# inapplicable fields loudly (strict per-kind validation), and a
# nightly layer must not smuggle inert ObservatoryRule scoping
# across the wire.
if layer.kind == WeightingLayerKind.NIGHTLY_DEWEIGHTING:
return {
"kind": _enum_value(layer.kind),
"max_gap_days": layer.max_gap_days,
}
return {
"kind": _enum_value(layer.kind),
"obs_code": layer.obs_code,
"sigma": list(layer.sigma),
"start_epoch_mjd_tdb": layer.start_epoch_mjd_tdb,
"end_epoch_mjd_tdb": layer.end_epoch_mjd_tdb,
"scale": layer.scale,
}
def _debiasing_to_dict(d: DebiasingConfig) -> dict[str, WireValue]:
return {
"enabled": d.enabled,
"resolution": _enum_value(d.resolution),
"bias_dat_path": d.bias_dat_path,
}
def _enum_value(v: Enum | str) -> str:
"""Accept either an Enum or a bare string; return a string."""
return str(v.value) if isinstance(v, Enum) else str(v)
def _origin_to_naif(o: Origin | str) -> int:
"""Internal — resolve an :class:`Origin` (or canonical name) to the
integer body code the binding wire format uses."""
from empyrean._convert import origin_to_naif
return origin_to_naif(o)
# ── Result types ─────────────────────────────────────────────
# Re-export StationBiases at the result module so callers can import
# it alongside the other OD types.
__all__ = []
# Any of the four orbit flavors that can come back from a determine /
# refine, depending on `ODConfig.output_representation`.
OrbitsTable = CartesianOrbits | KeplerianOrbits | CometaryOrbits | SphericalOrbits
[docs]
@dataclass
class EvaluateResult:
"""Result of orbit evaluation (residuals only, no fitting)."""
observations: ObservationResults
summary: ResidualSummary
[docs]
@dataclass
class SolvedCovariance:
"""Full tagged solved-parameter covariance from a wide OD fit.
Mirrors ``empyrean::SolvedCovariance``.
:attr:`matrix` is the real solved covariance, sized
``width × width``; parameters are located by the slot fields, never
by width (width 9 is Marsden-only OR one thrust segment). Canonical
order is ``[state 6 | Marsden 3 | DT 1 | AMRAT 1 | thrust 3×k]``. The
Δv axes are integration-frame components (see
:attr:`DetermineResult.dv_frame`).
"""
matrix: np.ndarray
"""The solved covariance, shaped ``(width, width)``."""
width: int
"""Solved width (6..=17 under the current engine)."""
marsden_slot: int | None
"""Column of the first Marsden coefficient, when Marsden was solved."""
dt_slot: int | None
"""Column of the DT scalar, when DT was solved."""
amrat_slot: int | None
"""Column of the AMRAT scalar, when AMRAT was solved."""
thrust_slots: list[tuple[int, int, int]]
"""Column triples of each fitted thrust Δv segment (one
``(i, i+1, i+2)`` per solved segment). Empty when no thrust was
solved."""
[docs]
@dataclass
class BandStat:
"""Per-band photometric fit statistics. Mirrors ``empyrean::BandStat``."""
band: str
"""Photometric band tag."""
n: int
"""Number of observations in this band."""
offset_applied: float
"""Band→V offset applied (mag)."""
mean_residual: float
"""Mean residual in V (mag)."""
rms: float
"""RMS residual in V (mag)."""
[docs]
@dataclass
class GateRecord:
"""One model-ladder gate decision from the photometric fit. Mirrors
``empyrean::GateRecord``."""
model: PhotometryModel
"""Model the gate evaluated."""
passed: bool
"""Whether the model was admitted."""
reason: str
"""Human-readable gate reason."""
[docs]
@dataclass
class PhotometryResult:
"""Post-OD photometric solution — an H/G fit over the arc's
magnitudes, run after the orbit is solved. Mirrors
``empyrean::PhotometryResult``.
Photometry has no astrometric partials, so it never touches the
state. H carries honest σ via :attr:`covariance`.
"""
h: float
"""Fitted absolute magnitude H (mag)."""
slope1: float
"""First slope parameter (G / G12 / G1 by model)."""
slope2: float
"""Second slope parameter (G2 for HG1G2; unused otherwise)."""
covariance: np.ndarray | None
"""Parameter covariance over (H, slope1, slope2), shaped ``(3, 3)``
when available. ``None`` otherwise."""
model_used: PhotometryModel
"""Model actually fitted (never :attr:`PhotometryModel.AUTO`)."""
reduced_chi2: float
"""Reduced χ² of the photometric fit over its used magnitudes."""
constraint_active: bool
"""Whether a simplex constraint was active on the fitted slopes."""
n_mags_used: int
"""Magnitudes used in the fit."""
n_mags_rejected_photometric: int
"""Magnitudes rejected by the photometric outlier pass."""
n_obs_without_mags: int
"""Observations carrying no magnitude."""
n_mags_from_astrometric_selected: int
"""Magnitudes drawn from astrometrically-selected observations."""
n_mags_from_astrometric_rejected: int
"""Magnitudes drawn from astrometrically-rejected observations."""
alpha_min_deg: float
"""Minimum phase angle of the fitted magnitudes (deg)."""
alpha_max_deg: float
"""Maximum phase angle of the fitted magnitudes (deg)."""
alpha_span_deg: float
"""Phase-angle span of the fitted magnitudes (deg)."""
per_band: list[BandStat]
"""Per-band statistics."""
gates: list[GateRecord]
"""Model-ladder gate records."""
n_mags_dropped_unconvertible: int
"""Magnitudes excluded from the fit because their photometric band
has no adopted V-band conversion (unknown/unspecified band codes,
comet total/nuclear magnitudes). Never silent: each exclusion is
counted here and the distinct offending band codes are listed in
:attr:`dropped_bands`. The observations' astrometry is
unaffected."""
dropped_bands: list[str]
"""Distinct band codes that were dropped, sorted."""
[docs]
@dataclass
class TrustGateEvent:
"""The intervening event named by an ``encounter_intervenes``
covariance-trust verdict."""
kind: str
"""``"close_approach"`` or ``"high_nonlinearity"``."""
epoch_mjd_tdb: float
"""Epoch of the event (MJD TDB)."""
body: str | None = None
"""Name of the approached body (close-approach events only)."""
distance_au: float | None = None
"""Approach distance at the signal (AU; close-approach only)."""
nonlinearity: float | None = None
"""Nonlinearity ratio at the crossing (high-nonlinearity only)."""
threshold: float | None = None
"""Threshold the nonlinearity exceeded (high-nonlinearity only)."""
[docs]
@dataclass
class CovarianceTrust:
"""Event-aware trust verdict on the delivered covariance, evaluated
over its validity window on the converged orbit.
``trusted``: no intervening close approach and a 6-state solve — the
linear covariance may be used as delivered. ``encounter_intervenes``:
a close approach (or high-nonlinearity crossing) lies inside the
window; do not extrapolate the linear covariance across it —
escalate to nonlinear uncertainty propagation (second-order when
:attr:`second_order_recoverable`, otherwise sampling).
``weakly_determined_high_n``: the fit solved more than the 6-state,
so the delivered 6×6 is a marginal of a wider fit (conservative
flag). A ``DetermineResult.covariance_trust`` of ``None`` means the
call path ran no gate — absence of a verdict is not trust."""
verdict: str
"""``"trusted"`` / ``"encounter_intervenes"`` /
``"weakly_determined_high_n"``."""
solved_width: int | None = None
"""Solved-for width of the fit the verdict refers to (absent for
``trusted``)."""
second_order_recoverable: bool | None = None
"""Whether a second-order state-only correction can recover the
encounter (``encounter_intervenes`` only)."""
event: TrustGateEvent | None = None
"""The earliest intervening event (``encounter_intervenes``
only)."""
[docs]
@dataclass
class StallDelivery:
"""The numbers behind a fit delivered by the stable-stall acceptance
— one whose final solve latched no convergence criterion but was
stationary to a small fraction of its own formal σ, with χ² clearing
the fit-quality bars.
Present as :attr:`DetermineResult.stall_delivery` exactly when
:attr:`DetermineResult.termination` is ``"stalled_delivered"``.
``None`` is the ordinary case; a consumer that wants the strict
convergence contract checks for it."""
underlying_stop: str
"""What the SOLVER reported before the acceptance overrode the
verdict — always ``"damping_exhausted"`` or
``"inner_trials_exhausted"``, the only two stops the acceptance
considers. Kept beside :attr:`DetermineResult.termination` so the
delivery's verdict and the solver's own are both readable and
neither is inferred from the other."""
gn_qnorm: float
"""q of the undamped Gauss-Newton step at the delivered iterate — the
quantity ``convergence_tol`` bounds and did not bound here. ``√q`` is
the remaining step in the fit's own formal σ."""
convergence_tol: float
"""The **resolved** tolerance the stall was judged against.
Delivered because :attr:`ODConfig.convergence_tol` carries a "use the
engine default" sentinel, so the value actually enforced is not
recoverable from the config that was passed in."""
optical_reduced_chi2: float
"""Reduced χ² of the OPTICAL family alone at the delivered iterate.
The full objective's reduced χ² and the bar both were held to are
already on :attr:`DetermineResult.acceptability`, and the delivery's
μ and iteration count are already
:attr:`DetermineResult.mu_final` and
:attr:`DetermineResult.final_solve_iterations` — a stall-delivered
fit's final solve IS the stalled solve."""
@property
def step_sigmas(self) -> float:
"""The remaining undamped Gauss-Newton step in units of the fit's
own formal 1σ, ``√q`` — the number to quote when saying how far
from stationary the delivered iterate is."""
return float(np.sqrt(self.gn_qnorm))
[docs]
@dataclass
class DetermineResult:
"""Result of orbit determination — returned by both
:func:`~empyrean.od.determine.determine` (full IOD + DC pipeline)
and :func:`~empyrean.od.determine.refine` (Bayesian-prior fit
against an existing orbit + covariance).
Mirrors the Rust wrapper's ``empyrean::DetermineResult``.
"""
orbit: OrbitsTable
"""Fitted orbit. Coordinate flavor matches
:attr:`ODConfig.output_representation`."""
observations: ObservationResults
"""Per-observation residuals + rejection / influence diagnostics."""
summary: ResidualSummary
iterations: int
update_norm: float
"""**NOT the quantity :attr:`ODConfig.convergence_tol` bounds.** The
μ-DAMPED last ACCEPTED step's q-norm — the size of the step the
solver actually took, after Levenberg-Marquardt damping.
The tolerance is tested on the *undamped* Gauss-Newton step, reported
as :attr:`gn_step_qnorm`. The two are incomparable in BOTH
directions: a converged fit routinely reports an ``update_norm``
orders of magnitude ABOVE ``convergence_tol`` (damping is light near
the optimum, so the step taken is large), while a solve thrashing to
``"damping_exhausted"`` reports one orders of magnitude BELOW it (μ
has crushed every step to nothing). Read :attr:`termination` and
:attr:`gn_step_qnorm` for convergence, or :attr:`acceptability` for
the fit-quality verdict.
On the Schur path (station-bias / nuisance fits) this carries that
loop's own step norm under the Schur complement instead — the
quantity IT tests against ``convergence_tol``. :attr:`termination`
says which loop ran.
Exactly ``0.0`` with :attr:`accepted_steps` ``== 0`` means the solver
latched a criterion at its STARTING point and never accepted a step:
the incoming iterate was already stationary. It does not mean "a step
of size zero was taken"."""
converged: bool
covariance: np.ndarray
"""Fitted 6×6 state covariance, in :attr:`covariance_representation`."""
covariance_representation: CovarianceRepresentation
covariance_9x9: np.ndarray | None
"""Full 9×9 covariance over (state, A1, A2, A3) when solving for non-grav."""
non_grav_delta: np.ndarray | None
"""Cumulative non-grav corrections (ΔA1, ΔA2, ΔA3) when present."""
rejection_passes: int
num_oppositions_fit: int
force_model_used: ForceModelTier
solve_for_used: SolveForParams
acceptability: AcceptabilityReport
station_biases: StationBiases
"""Per-station fitted nuisance biases when
:attr:`ODConfig.fit_station_biases` was active. Empty quivr table
otherwise."""
solved_covariance: SolvedCovariance | None
"""Full tagged solved-parameter covariance when the fit solved any
wide axis (Marsden / DT / AMRAT / thrust). ``None`` for a state-only
fit — read it, not the width, to locate solved parameters."""
dt_delta: float | None
"""Cumulative non-grav time-delay correction ΔDT (days), when DT was
solved. ``None`` otherwise."""
amrat_delta: float | None
"""Cumulative SRP AMRAT correction (m²/kg), when AMRAT was solved.
``None`` otherwise."""
thrust_delta_m_per_s: np.ndarray | None
"""Per-segment fitted thrust Δv (m/s), shaped ``(k, 3)`` and
expressed in :attr:`dv_frame`. ``None`` when no thrust was solved."""
dv_frame: str | None
"""Integration frame the thrust Δv components are expressed in
(``"icrf"`` / ``"eclipticj2000"`` / ``"itrf93"``). ``None`` when no
thrust was solved."""
photometry: PhotometryResult | None
"""Post-OD photometric solution when photometry was requested and
ran. ``None`` otherwise."""
covariance_trust: CovarianceTrust | None
"""Event-aware trust verdict on the delivered covariance. ``None``
when the call path ran no trust gate — absence of a verdict is not
trust."""
dispositions: SolveFor
"""What the fit actually did with each parameter axis.
The partition the engine ran, not the one that was requested: an
axis can be requested and then not opened. Read this rather than
:attr:`solve_for_used` to learn whether an axis was *considered* —
not estimated, but contributing its prior uncertainty to the
posterior through its measurement partials — because a solved
covariance's slot tags record only what occupied a column, and a
considered axis occupies none.
It is also what tells you whether re-attaching a prior to an axis
would double-count it."""
warnings: list[str]
"""Covariance the fit was given and deliberately did not use.
Delivered as payload rather than written to a log, because a dropped
prior cross term changes how the σ for that slot should be read.
Empty when the fit used everything it was given."""
termination: str | None
"""Which criterion ended the solve that produced the published state,
as a stable lowercase tag — ``"gradient_tolerance"``,
``"step_tolerance"``, ``"cost_tolerance"``, ``"max_iterations"``,
``"damping_exhausted"``, ``"inner_trials_exhausted"``,
``"stalled_delivered"``, ``"schur_step_tolerance"``, or
``"unrecognized"`` for a stop this build cannot name.
``None`` only on a result that did not come from a solver run; every
fit this package delivers names its stop. An unrecognized stop is
``"unrecognized"`` rather than ``None``: a stop that fired but cannot
be named is not the same as no stop at all."""
gn_step_qnorm: float | None
"""Quadratic form of the **undamped** Gauss-Newton step — the only
number here comparable to :attr:`ODConfig.convergence_tol`, and the
one the solver's step test is decided on. ``√q`` is the remaining
step measured in the fit's own formal σ.
``q <= convergence_tol`` is guaranteed only under
``"step_tolerance"`` (or ``"schur_step_tolerance"``). Under the
gradient, cost and stall verdicts the fit is delivered as converged
with ``q`` legitimately above the tolerance.
``None`` when the undamped system was singular at the returned point,
when the Schur loop ran at solved width > 6 (its only step norm there
is λ-shrunken, which is what this field promises it is not), or on a
non-solver result."""
mu_final: float | None
"""Levenberg-Marquardt damping at termination. Large μ beside a small
:attr:`update_norm` is the signature of a stalled solve rather than a
converged one.
``None`` on the Schur path, which runs its own λ schedule and forms
no comparable quantity — so its absence, alongside
``termination == "schur_step_tolerance"``, is what identifies that
path — and on non-solver results."""
accepted_steps: int
"""Steps the solve actually ACCEPTED — trials that passed the
gain-ratio test and moved the iterate.
``0`` says the solver latched a criterion at its starting point and
never moved: the incoming orbit was already stationary for this
observation set. That is what disambiguates :attr:`update_norm`'s
``0.0`` sentinel. Counted per SOLVE, so it describes the same solve
as :attr:`final_solve_iterations`, never a total across the
rejection–refit passes."""
final_solve_iterations: int | None
"""Iterations spent by the solve that produced the published state —
the one :attr:`update_norm`, :attr:`termination`,
:attr:`gn_step_qnorm` and :attr:`accepted_steps` all describe.
Additive disambiguator for :attr:`iterations`, whose composition
varies by entry point; this one never does. ``None`` on a result that
did not come from a solver run."""
stall_delivery: StallDelivery | None
"""Set when this fit was delivered by the stable-stall acceptance
rather than by a convergence criterion — see :class:`StallDelivery`.
``None`` is the ordinary case. A consumer that wants the strict
convergence contract checks ``stall_delivery is None``."""
[docs]
@dataclass
class DetermineFailure:
"""Why one object in a batch produced no orbit."""
object_id: str
"""ADES object identifier whose fit failed."""
message: str
"""The engine's message."""
kind: str
"""Classified cause, as a stable snake_case name — ``"iod"``,
``"od"``, ``"radar_only"``, ``"observer_construction"``,
``"earth_orientation_coverage"``, ``"observation_conversion"``,
``"unsupported_coordinate_system"``, ``"duplicate_obs_ids"``,
``"non_grav_not_recovered"``, or ``"unknown_<code>"`` for a cause
this build of the bindings does not name. Branch on this rather than
on :attr:`message`."""
[docs]
@dataclass
class DetermineResults:
"""The result of a batch orbit determination — every object the
observations grouped into, delivered or not.
:func:`~empyrean.od.determine.determine` fits **every** object in the
input. Three tables describe the run:
- :attr:`orbits` — one row per object that produced an orbit.
- :attr:`summary` — one row per **input** object, so an object that
failed is a row saying so rather than an absence.
- :attr:`residuals` — every delivered fit's per-observation rows,
each tagged with the ``object_id`` it belongs to.
Index by object identifier for the full single-object view::
fits = determine(observations)
print(len(fits), "object(s);", len(fits.delivered), "delivered")
yr4 = fits["2024 YR4"] # DetermineResult
print(yr4.acceptability.extrapolation_acceptable)
Fitting one object is the one-entry case, not a separate call:
:meth:`single` unwraps it and refuses (loudly) if the batch turned
out to hold more than one.
"""
orbits: OrbitsTable
"""Fitted orbits — one row per **delivered** object, carrying
``object_id``, covariance, and the fitted non-grav / SRP slots. Feed
straight back into propagation or ephemeris generation."""
summary: FitSummary
"""One row per **input** object, delivered or failed."""
residuals: ObservationResults
"""Per-observation residuals of every delivered fit, each row tagged
with its ``object_id``."""
failures: dict[str, DetermineFailure]
"""Objects that produced no orbit, keyed by identifier."""
unmatched_orbit_ids: list[str]
"""Seed orbits that matched no observation group and therefore
constrained nothing. Reported rather than dropped."""
_results: dict[str, DetermineResult] = field(default_factory=dict, repr=False)
"""Per-object single-fit views, keyed by identifier. Read through
:meth:`__getitem__`."""
def __len__(self) -> int:
"""Number of objects the observations grouped into — delivered
and failed alike."""
return len(self.summary)
def __iter__(self) -> Iterator[str]:
"""Iterate the object identifiers, in table order."""
return iter(self.object_ids)
def __contains__(self, object_id: str) -> bool:
return object_id in self._results or object_id in self.failures
def __getitem__(self, object_id: str) -> DetermineResult:
"""The single-object view of one delivered fit.
Raises ``KeyError`` for an object the batch never saw, and
``ValueError`` — naming the reason — for one that failed. A
failed object never returns an empty result.
"""
if object_id in self._results:
return self._results[object_id]
failure = self.failures.get(object_id)
if failure is not None:
raise ValueError(f"orbit determination failed for {object_id}: {failure.message}")
raise KeyError(f"no object {object_id!r} in this batch; it holds {self.object_ids}")
@property
def object_ids(self) -> list[str]:
"""Every object identifier, in table order."""
ids: list[str] = self.summary.object_id.to_pylist()
return ids
@property
def delivered(self) -> list[str]:
"""Identifiers of the objects that produced an orbit."""
return list(self._results)
@property
def all_failed(self) -> bool:
"""True when the batch ran but no object produced an orbit. The
per-object :attr:`failures` are the diagnosis."""
return len(self.summary) > 0 and not self._results
[docs]
def single(self) -> DetermineResult:
"""The one fit, for the common single-object call.
Raises — loudly — when the batch holds anything other than
exactly one delivered object: zero objects, more than one (naming
them, so the caller can index instead), or one whose fit failed.
It never chooses among several fits on the caller's behalf.
"""
ids = self.object_ids
if len(ids) == 1:
return self[ids[0]]
if not ids:
raise ValueError(
"orbit determination produced no objects: the observations carried no rows to group"
)
raise ValueError(
f"orbit determination fitted {len(ids)} objects "
f"({', '.join(ids)}); single() refuses to choose one — "
f"iterate the batch or index it by object_id"
)