use std::collections::HashMap;
use pleiades_backend::{Angle, CelestialBody, CustomBodyId};
use pleiades_compression::{cartesian_state_to_spherical, CartesianState, CompressedArtifact};
use pleiades_jpl::{production_holdout_corpus, reference_snapshot, SnapshotEntry};
use crate::regenerate::{build_packaged_artifact, coordinates, normalize_lookup_instant};
use crate::AU_IN_KM;
#[derive(Clone, Debug)]
pub struct BodyChannelError {
pub body: CelestialBody,
pub comparison_count: usize,
pub max_longitude_arcsec: f64,
pub rms_longitude_arcsec: f64,
pub max_latitude_arcsec: f64,
pub rms_latitude_arcsec: f64,
pub max_distance_km: f64,
pub rms_distance_km: f64,
pub max_lon_speed_arcsec_per_day: f64,
pub max_lat_speed_arcsec_per_day: f64,
pub max_radial_speed_au_per_day: f64,
}
impl BodyChannelError {
fn label(&self) -> String {
format!("{:?}", self.body)
}
pub fn summary_line(&self) -> String {
format!(
"{}: n={} max_lon={:.4} arcsec rms_lon={:.4} arcsec max_lat={:.4} arcsec rms_lat={:.4} arcsec max_dist={:.3} km rms_dist={:.3} km max_lon_speed={:.4} arcsec/day max_lat_speed={:.4} arcsec/day max_radial_speed={:.6} AU/day",
self.label(),
self.comparison_count,
self.max_longitude_arcsec,
self.rms_longitude_arcsec,
self.max_latitude_arcsec,
self.rms_latitude_arcsec,
self.max_distance_km,
self.rms_distance_km,
self.max_lon_speed_arcsec_per_day,
self.max_lat_speed_arcsec_per_day,
self.max_radial_speed_au_per_day,
)
}
}
struct BodyAccumulator {
body: CelestialBody,
count: usize,
max_lon_arcsec: f64,
sum_sq_lon_arcsec: f64,
max_lat_arcsec: f64,
sum_sq_lat_arcsec: f64,
max_dist_km: f64,
sum_sq_dist_km: f64,
max_lon_speed_arcsec_per_day: f64,
max_lat_speed_arcsec_per_day: f64,
max_radial_speed_au_per_day: f64,
}
impl BodyAccumulator {
fn new(body: CelestialBody) -> Self {
Self {
body,
count: 0,
max_lon_arcsec: 0.0,
sum_sq_lon_arcsec: 0.0,
max_lat_arcsec: 0.0,
sum_sq_lat_arcsec: 0.0,
max_dist_km: 0.0,
sum_sq_dist_km: 0.0,
max_lon_speed_arcsec_per_day: 0.0,
max_lat_speed_arcsec_per_day: 0.0,
max_radial_speed_au_per_day: 0.0,
}
}
fn accumulate(&mut self, lon_arcsec: f64, lat_arcsec: f64, dist_km: Option<f64>) {
self.count += 1;
self.max_lon_arcsec = self.max_lon_arcsec.max(lon_arcsec);
self.sum_sq_lon_arcsec += lon_arcsec * lon_arcsec;
self.max_lat_arcsec = self.max_lat_arcsec.max(lat_arcsec);
self.sum_sq_lat_arcsec += lat_arcsec * lat_arcsec;
if let Some(d) = dist_km {
self.max_dist_km = self.max_dist_km.max(d);
self.sum_sq_dist_km += d * d;
}
}
fn accumulate_speed(
&mut self,
lon_speed_arcsec_per_day: f64,
lat_speed_arcsec_per_day: f64,
radial_speed_au_per_day: f64,
) {
self.max_lon_speed_arcsec_per_day = self
.max_lon_speed_arcsec_per_day
.max(lon_speed_arcsec_per_day);
self.max_lat_speed_arcsec_per_day = self
.max_lat_speed_arcsec_per_day
.max(lat_speed_arcsec_per_day);
self.max_radial_speed_au_per_day = self
.max_radial_speed_au_per_day
.max(radial_speed_au_per_day);
}
fn finish(self) -> BodyChannelError {
let n = self.count as f64;
let rms = |sum_sq: f64| if n > 0.0 { (sum_sq / n).sqrt() } else { 0.0 };
BodyChannelError {
body: self.body,
comparison_count: self.count,
max_longitude_arcsec: self.max_lon_arcsec,
rms_longitude_arcsec: rms(self.sum_sq_lon_arcsec),
max_latitude_arcsec: self.max_lat_arcsec,
rms_latitude_arcsec: rms(self.sum_sq_lat_arcsec),
max_distance_km: self.max_dist_km,
rms_distance_km: rms(self.sum_sq_dist_km),
max_lon_speed_arcsec_per_day: self.max_lon_speed_arcsec_per_day,
max_lat_speed_arcsec_per_day: self.max_lat_speed_arcsec_per_day,
max_radial_speed_au_per_day: self.max_radial_speed_au_per_day,
}
}
}
pub fn accuracy_baseline_against(
holdout: &[SnapshotEntry],
artifact: &CompressedArtifact,
) -> Vec<BodyChannelError> {
let mut accumulators: HashMap<String, BodyAccumulator> = HashMap::new();
for entry in holdout {
let key = format!("{:?}", entry.body);
let acc = accumulators
.entry(key)
.or_insert_with(|| BodyAccumulator::new(entry.body.clone()));
let lookup_instant = normalize_lookup_instant(entry.epoch);
let artifact_result = artifact.lookup_ecliptic(&entry.body, lookup_instant);
let artifact_coords = match artifact_result {
Ok(coords) => coords,
Err(_) => continue, };
let holdout_coords = coordinates(entry);
let lon_diff_deg = Angle::from_degrees(
artifact_coords.longitude.degrees() - holdout_coords.longitude.degrees(),
)
.normalized_signed()
.degrees();
let lon_arcsec = lon_diff_deg.abs() * 3600.0;
let lat_arcsec =
(artifact_coords.latitude.degrees() - holdout_coords.latitude.degrees()).abs() * 3600.0;
let dist_km = match (artifact_coords.distance_au, holdout_coords.distance_au) {
(Some(a), Some(h)) => Some((a - h).abs() * AU_IN_KM),
_ => None,
};
acc.accumulate(lon_arcsec, lat_arcsec, dist_km);
if let (Some(vx), Some(vy), Some(vz)) = (entry.vx_km_s, entry.vy_km_s, entry.vz_km_s) {
let pos_au = [
entry.x_km / AU_IN_KM,
entry.y_km / AU_IN_KM,
entry.z_km / AU_IN_KM,
];
let vel_au_per_day = [
vx * 86400.0 / AU_IN_KM,
vy * 86400.0 / AU_IN_KM,
vz * 86400.0 / AU_IN_KM,
];
let truth_spherical = cartesian_state_to_spherical(CartesianState {
pos_au,
vel_au_per_day,
});
let truth_lon_deg_per_day = truth_spherical.lon_rate_rad_per_day.to_degrees();
let truth_lat_deg_per_day = truth_spherical.lat_rate_rad_per_day.to_degrees();
let truth_radial_au_per_day = truth_spherical.dist_rate_au_per_day;
if let Ok(motion) = artifact.lookup_motion(&entry.body, lookup_instant) {
let art_lon = motion.longitude_deg_per_day.unwrap_or(0.0);
let art_lat = motion.latitude_deg_per_day.unwrap_or(0.0);
let art_radial = motion.distance_au_per_day.unwrap_or(0.0);
let lon_speed_arcsec = (art_lon - truth_lon_deg_per_day).abs() * 3600.0;
let lat_speed_arcsec = (art_lat - truth_lat_deg_per_day).abs() * 3600.0;
let radial_speed_au = (art_radial - truth_radial_au_per_day).abs();
acc.accumulate_speed(lon_speed_arcsec, lat_speed_arcsec, radial_speed_au);
}
}
}
let mut results: Vec<BodyChannelError> = accumulators
.into_values()
.filter(|acc| acc.count > 0)
.map(|acc| acc.finish())
.collect();
results.sort_by_key(|a| a.label());
results
}
pub fn packaged_artifact_accuracy_baseline() -> Vec<BodyChannelError> {
let holdout = production_holdout_corpus();
let artifact = build_packaged_artifact();
accuracy_baseline_against(holdout, &artifact)
}
#[allow(dead_code)]
pub(crate) fn eros_self_consistency_max_longitude_arcsec() -> f64 {
let eros_body = CelestialBody::Custom(CustomBodyId::new("asteroid", "433-Eros"));
let artifact = build_packaged_artifact();
let mut max_lon_arcsec: f64 = 0.0;
for entry in reference_snapshot() {
if entry.body != eros_body {
continue;
}
let lookup_instant = normalize_lookup_instant(entry.epoch);
let artifact_coords = match artifact.lookup_ecliptic(&eros_body, lookup_instant) {
Ok(c) => c,
Err(_) => continue,
};
let snapshot_coords = coordinates(entry);
let lon_diff_deg = Angle::from_degrees(
artifact_coords.longitude.degrees() - snapshot_coords.longitude.degrees(),
)
.normalized_signed()
.degrees();
let lon_arcsec = lon_diff_deg.abs() * 3600.0;
max_lon_arcsec = max_lon_arcsec.max(lon_arcsec);
}
max_lon_arcsec
}
#[cfg(test)]
mod tests {
use pleiades_backend::{Instant, JulianDay, TimeScale};
use pleiades_compression::{
ArtifactHeader, BodyArtifact, ChannelKind, CompressedArtifact, PolynomialChannel, Segment,
};
use super::*;
fn synthetic_holdout_with_velocity() -> Vec<SnapshotEntry> {
let au = AU_IN_KM;
let vy_km_s = 0.01 * au / 86400.0;
vec![SnapshotEntry {
body: CelestialBody::Sun,
epoch: Instant::new(JulianDay::from_days(2_451_545.0), TimeScale::Tt),
x_km: au,
y_km: 0.0,
z_km: 0.0,
vx_km_s: Some(0.0),
vy_km_s: Some(vy_km_s),
vz_km_s: Some(0.0),
}]
}
fn synthetic_artifact_linear() -> CompressedArtifact {
let jd0 = 2_451_545.0_f64;
let span_days = 10.0_f64;
let start = Instant::new(JulianDay::from_days(jd0), TimeScale::Tt);
let end = Instant::new(JulianDay::from_days(jd0 + span_days), TimeScale::Tt);
let lon_rate_deg_per_day = 0.01_f64.to_degrees();
let lon_end = lon_rate_deg_per_day * span_days;
let segment = Segment::new(
start,
end,
vec![
PolynomialChannel::linear(ChannelKind::Longitude, 9, 0.0, lon_end),
PolynomialChannel::linear(ChannelKind::Latitude, 9, 0.0, 0.0),
PolynomialChannel::linear(ChannelKind::DistanceAu, 10, 1.0, 1.0),
],
);
CompressedArtifact::new(
ArtifactHeader::new(
"synthetic-linear-test",
"synthetic linear velocity test source",
),
vec![BodyArtifact::new(CelestialBody::Sun, vec![segment])],
)
}
fn synthetic_holdout() -> Vec<SnapshotEntry> {
let au = AU_IN_KM;
vec![SnapshotEntry {
body: CelestialBody::Sun,
epoch: Instant::new(JulianDay::from_days(2_451_545.0), TimeScale::Tt),
x_km: au,
y_km: 0.0,
z_km: 0.0,
vx_km_s: None,
vy_km_s: None,
vz_km_s: None,
}]
}
fn synthetic_artifact() -> CompressedArtifact {
let jd = JulianDay::from_days(2_451_545.0);
let instant = Instant::new(jd, TimeScale::Tt);
let segment = Segment::new(
instant,
instant,
vec![
PolynomialChannel::linear(ChannelKind::Longitude, 9, 0.0, 0.0),
PolynomialChannel::linear(ChannelKind::Latitude, 9, 0.0, 0.0),
PolynomialChannel::linear(ChannelKind::DistanceAu, 10, 1.0, 1.0),
],
);
CompressedArtifact::new(
ArtifactHeader::new("synthetic-test", "synthetic test source"),
vec![BodyArtifact::new(CelestialBody::Sun, vec![segment])],
)
}
#[test]
fn baseline_reports_speed_error_fields() {
let errors = accuracy_baseline_against(
&synthetic_holdout_with_velocity(),
&synthetic_artifact_linear(),
);
assert_eq!(
errors.len(),
1,
"expected exactly 1 body in synthetic speed baseline"
);
let sun = &errors[0];
assert_eq!(sun.comparison_count, 1, "Sun should have 1 comparison");
assert!(
sun.max_lon_speed_arcsec_per_day < 1e-3,
"Sun max longitude speed error too large: {} arcsec/day",
sun.max_lon_speed_arcsec_per_day
);
assert!(
sun.max_lat_speed_arcsec_per_day < 1e-3,
"Sun max latitude speed error too large: {} arcsec/day",
sun.max_lat_speed_arcsec_per_day
);
assert!(
sun.max_radial_speed_au_per_day < 1e-6,
"Sun max radial speed error too large: {} AU/day",
sun.max_radial_speed_au_per_day
);
}
#[test]
fn baseline_reports_zero_error_for_an_artifact_that_matches_holdout() {
let errors = accuracy_baseline_against(&synthetic_holdout(), &synthetic_artifact());
assert_eq!(
errors.len(),
1,
"expected exactly 1 body in synthetic baseline"
);
let sun = &errors[0];
assert_eq!(sun.comparison_count, 1, "Sun should have 1 comparison");
assert!(
sun.max_longitude_arcsec < 1e-3,
"Sun max longitude error too large: {} arcsec",
sun.max_longitude_arcsec
);
assert!(
sun.max_latitude_arcsec < 1e-3,
"Sun max latitude error too large: {} arcsec",
sun.max_latitude_arcsec
);
assert!(
sun.max_distance_km < 1.0,
"Sun max distance error too large: {} km",
sun.max_distance_km
);
}
#[test]
fn baseline_excludes_body_when_all_lookups_fail() {
let holdout = synthetic_holdout(); let empty_artifact = CompressedArtifact::new(
ArtifactHeader::new("empty-test", "empty test source"),
vec![],
);
let errors = accuracy_baseline_against(&holdout, &empty_artifact);
assert!(
errors.is_empty(),
"vacuity guard must exclude bodies with zero successful comparisons; got {} entries",
errors.len()
);
}
#[test]
fn packaged_artifact_baseline_is_non_vacuous() {
let errors = packaged_artifact_accuracy_baseline();
let expected_bodies = [
CelestialBody::Sun,
CelestialBody::Moon,
CelestialBody::Mercury,
CelestialBody::Venus,
CelestialBody::Mars,
CelestialBody::Jupiter,
CelestialBody::Saturn,
CelestialBody::Uranus,
CelestialBody::Neptune,
CelestialBody::Pluto,
];
assert_eq!(
errors.len(),
10,
"expected 10 base bodies in packaged baseline; got {}: {:?}",
errors.len(),
errors
.iter()
.map(|e| format!("{:?}", e.body))
.collect::<Vec<_>>()
);
for body in &expected_bodies {
let entry = errors
.iter()
.find(|e| &e.body == body)
.unwrap_or_else(|| panic!("{body:?} must appear in the packaged baseline"));
assert!(
entry.comparison_count > 0,
"{body:?} must have at least one successful comparison (got 0 — vacuous baseline)"
);
}
let sub_arcsec_bodies = [
CelestialBody::Sun,
CelestialBody::Moon,
CelestialBody::Mercury,
CelestialBody::Venus,
CelestialBody::Mars,
];
for body in &sub_arcsec_bodies {
let entry = errors.iter().find(|e| &e.body == body).unwrap();
assert!(
entry.max_longitude_arcsec < 1.0,
"{body:?} max longitude error must be <1 arcsec (got {:.6}\")",
entry.max_longitude_arcsec
);
}
let outer_bodies = [
CelestialBody::Jupiter,
CelestialBody::Saturn,
CelestialBody::Uranus,
CelestialBody::Neptune,
CelestialBody::Pluto,
];
for body in &outer_bodies {
let entry = errors.iter().find(|e| &e.body == body).unwrap();
assert!(
entry.max_longitude_arcsec < 1.0,
"{body:?} max longitude error must be <1 arcsec after SP2 heliocentric reframe (got {:.6}\")",
entry.max_longitude_arcsec
);
}
let uranus = errors
.iter()
.find(|e| e.body == CelestialBody::Uranus)
.expect("Uranus must appear in the packaged baseline");
assert!(
uranus.max_longitude_arcsec > 0.0001,
"Uranus max longitude error must be >0.0001\" (got {:.6}\" — baseline may be vacuous)",
uranus.max_longitude_arcsec
);
assert!(
uranus.max_longitude_arcsec < 1.0,
"Uranus max longitude error must be <1\" after SP2 heliocentric reframe (got {:.4}\" — reframe may be broken)",
uranus.max_longitude_arcsec
);
}
#[test]
fn outer_planet_longitude_meets_astrology_grade_envelope() {
let baseline = crate::accuracy_baseline::packaged_artifact_accuracy_baseline();
for body_error in &baseline {
let c = crate::thresholds::accuracy_ceiling(&body_error.body).lon_arcsec;
assert!(
body_error.max_longitude_arcsec <= c,
"{:?} longitude {:.3}\" exceeds ceiling {:.1}\"",
body_error.body,
body_error.max_longitude_arcsec,
c
);
}
}
#[test]
fn all_channels_within_published_ceilings_for_major_bodies() {
let baseline = packaged_artifact_accuracy_baseline();
for e in &baseline {
let c = crate::thresholds::accuracy_ceiling(&e.body);
assert!(
e.max_longitude_arcsec <= c.lon_arcsec,
"{:?} lon {:.4}\" exceeds ceiling {:.1}\"",
e.body,
e.max_longitude_arcsec,
c.lon_arcsec
);
assert!(
e.max_latitude_arcsec <= c.lat_arcsec,
"{:?} lat {:.4}\" exceeds ceiling {:.1}\"",
e.body,
e.max_latitude_arcsec,
c.lat_arcsec
);
assert!(
e.max_distance_km <= c.dist_km,
"{:?} dist {:.3} km exceeds ceiling {:.0} km",
e.body,
e.max_distance_km,
c.dist_km
);
assert!(
e.max_lon_speed_arcsec_per_day <= c.lon_speed_arcsec_per_day,
"{:?} lon speed {:.4} arcsec/day exceeds ceiling {:.2} arcsec/day",
e.body,
e.max_lon_speed_arcsec_per_day,
c.lon_speed_arcsec_per_day
);
assert!(
e.max_lat_speed_arcsec_per_day <= c.lat_speed_arcsec_per_day,
"{:?} lat speed {:.4} arcsec/day exceeds ceiling {:.2} arcsec/day",
e.body,
e.max_lat_speed_arcsec_per_day,
c.lat_speed_arcsec_per_day
);
assert!(
e.max_radial_speed_au_per_day <= c.radial_speed_au_per_day,
"{:?} radial speed {:.6} AU/day exceeds ceiling {:.2e} AU/day",
e.body,
e.max_radial_speed_au_per_day,
c.radial_speed_au_per_day
);
}
}
#[test]
fn encoded_artifact_within_size_budget() {
let bytes_len = crate::data::packaged_artifact_bytes().len();
assert!(
bytes_len <= crate::thresholds::PACKAGED_BUDGETS.max_encoded_bytes,
"encoded artifact {} bytes exceeds budget {} bytes",
bytes_len,
crate::thresholds::PACKAGED_BUDGETS.max_encoded_bytes
);
}
#[test]
fn speed_channels_are_non_vacuous_for_major_bodies() {
let errors = packaged_artifact_accuracy_baseline();
let inner_bodies = [
CelestialBody::Sun,
CelestialBody::Mercury,
CelestialBody::Venus,
CelestialBody::Mars,
];
for body in &inner_bodies {
let e = errors
.iter()
.find(|e| &e.body == body)
.unwrap_or_else(|| panic!("{body:?} must appear in the packaged baseline"));
assert!(
e.max_lon_speed_arcsec_per_day > 0.0001,
"{body:?} max_lon_speed {:.6} arcsec/day is not strictly positive (speed channel may be vacuous — all velocity rows skipped?)",
e.max_lon_speed_arcsec_per_day
);
}
let moon = errors
.iter()
.find(|e| e.body == CelestialBody::Moon)
.expect("Moon must appear in the packaged baseline");
assert!(
moon.max_lon_speed_arcsec_per_day > 0.005,
"Moon max_lon_speed {:.6} arcsec/day is not strictly positive (speed channel may be vacuous — all velocity rows skipped?)",
moon.max_lon_speed_arcsec_per_day
);
let outer_bodies = [
CelestialBody::Jupiter,
CelestialBody::Saturn,
CelestialBody::Uranus,
CelestialBody::Neptune,
CelestialBody::Pluto,
];
for body in &outer_bodies {
let e = errors
.iter()
.find(|e| &e.body == body)
.unwrap_or_else(|| panic!("{body:?} must appear in the packaged baseline"));
assert!(
e.max_lon_speed_arcsec_per_day > 0.0001,
"{body:?} max_lon_speed {:.6} arcsec/day is not strictly positive (speed channel may be vacuous — all velocity rows skipped?)",
e.max_lon_speed_arcsec_per_day
);
}
}
#[test]
#[ignore = "maintainer helper: prints the Eros self-consistency max longitude error"]
fn print_eros_self_consistency_max_longitude_arcsec() {
let v = crate::accuracy_baseline::eros_self_consistency_max_longitude_arcsec();
eprintln!("EROS_SELF_CONSISTENCY_MAX_LON_ARCSEC = {v:.6}\"");
}
#[test]
fn eros_round_trips_against_its_reference_snapshot_within_documented_target() {
let eros_body = pleiades_backend::CelestialBody::Custom(
pleiades_backend::CustomBodyId::new("asteroid", "433-Eros"),
);
let eros_row_count = pleiades_jpl::reference_snapshot()
.iter()
.filter(|e| e.body == eros_body)
.count();
assert!(
eros_row_count > 0,
"reference snapshot contains no Eros rows — self-consistency check would be vacuous (iterate zero rows → max=0.0 → trivially passes)"
);
let ceiling = crate::thresholds::accuracy_ceiling(&eros_body);
let max_lon_arcsec = crate::accuracy_baseline::eros_self_consistency_max_longitude_arcsec();
assert!(
max_lon_arcsec <= ceiling.lon_arcsec,
"Eros self-consistency {max_lon_arcsec:.4}\" > {:.1}\" (artifact does not reproduce the reference snapshot it was fit from within the Asteroid-class ceiling)",
ceiling.lon_arcsec
);
}
}