kshana 0.37.0

// SPDX-License-Identifier: AGPL-3.0-only
//! Single-run orchestration: driving a parsed scenario through the estimator and
//! scoring it into a result.
use crate::estimator::HoldoverEstimator;
use crate::fom::{score, Sample};
use crate::kalman::KalmanClock;
use crate::models::{ClockModel, ErrorModel};
use crate::portable_math::PortableFloat;
use crate::report::{ClockRun, RunResult};
use crate::scenario::{ClockCfg, GnssState, Scenario};
use crate::security::{spoof_detection_score, SPOOF_DETECT_K, SPOOF_MONITOR_S};
use rand::SeedableRng;
use rand_chacha::ChaCha8Rng;

/// GNSS-disciplined phase measurement noise (variance, s^2). Represents the
/// timing noise on the truth observation available while GNSS is nominal
/// (~0.1 ns, 1-sigma); it sets the synchronised covariance floor for the filter.
pub(crate) const PHASE_MEAS_VAR_S2: f64 = 1e-20;

/// 3-sigma protection level for the integrity check.
pub(crate) const PROTECTION_K: f64 = 3.0;

/// Monte-Carlo ensemble size and length for the filter-health (NIS/NEES) check.
/// `seeds × steps` pooled samples (≈ 12 800) give a tight χ² band while staying
/// well under a millisecond per clock.
const HEALTH_STEPS: usize = 200;
const HEALTH_SEEDS: usize = 64;
/// The flicker-aware check runs a 24-state filter, so its ensemble is smaller (960 pooled
/// samples; the NIS band is about +-9 %) to keep an ensemble of runs cheap.
const HEALTH_STEPS_FLICKER: usize = 60;
const HEALTH_SEEDS_FLICKER: usize = 16;
/// Decorrelate the health ensemble's seed stream from the scenario's run seed.
const HEALTH_SEED_SALT: u64 = 0x0F11_7E12_8EA1_7777;

pub(crate) fn run_clock(scn: &Scenario, cfg: &ClockCfg, seed: u64) -> ClockRun {
    run_clock_probed(scn, cfg, seed, &mut |_, _| {})
}

/// [`run_clock`] with a probe called for every outage sample with the predictor's timing
/// error (s) and the 1-sigma bound (s) the integrity check scales by [`PROTECTION_K`]. The
/// probe only observes; `run_clock` passes a no-op, so the result is the same.
pub(crate) fn run_clock_probed(
    scn: &Scenario,
    cfg: &ClockCfg,
    seed: u64,
    probe: &mut dyn FnMut(f64, f64),
) -> ClockRun {
    let mut rng = ChaCha8Rng::seed_from_u64(seed);
    let mut clock = ClockModel::new(&cfg.id, &cfg.provenance, cfg.y0, cfg.q_wf, cfg.q_rw)
        .with_drift(cfg.drift)
        .with_flicker(cfg.flicker_floor);
    let mut est = HoldoverEstimator::new();
    // Kalman estimator running alongside the analytic predictor: it is disciplined
    // to the truth while GNSS is nominal and coasts open-loop during the outage,
    // where its 1-sigma phase uncertainty is the protection bound. Integrity is the
    // fraction of outage samples whose actual error stays inside the k-sigma bound.
    let mut kf = KalmanClock::new(cfg.q_wf, cfg.q_rw, PHASE_MEAS_VAR_S2);
    let dt = scn.time.step_s;
    let n = (scn.time.duration_s / dt).round() as usize;
    let mut series = Vec::with_capacity(n + 1);
    // Raw clock phase over the whole run, for the Allan-deviation curve.
    let mut phase = Vec::with_capacity(n + 1);
    let (mut outage_samples, mut contained) = (0u64, 0u64);
    // Steps since the predictor was last anchored to the observed phase.
    let mut steps_since_sync = 0usize;
    for i in 0..=n {
        let t = i as f64 * dt;
        if i > 0 {
            clock.step(dt, &mut rng);
            kf.predict(dt);
            steps_since_sync += 1;
        }
        let gnss = scn.gnss.state_at(t);
        let ph = clock.phase();
        phase.push(ph);
        let err_s = est.timing_error(t, ph, clock.det_freq(), clock.drift_rate(), gnss);
        if gnss == GnssState::Nominal {
            // Truth is observed: the timing error is zero and the filter re-syncs.
            kf.update(0.0);
            steps_since_sync = 0;
        } else {
            outage_samples += 1;
            // The series is the deterministic predictor's error, which knows nothing of the
            // flicker frequency, so the bound adds the flicker phase variance accumulated
            // since the last anchor (unconditional: the 2-state filter's own variance plus
            // the truth bank's). Zero without a flicker floor.
            let var = kf.phase_sigma().ppowi(2)
                + clock.flicker_phase_variance_after_steps(steps_since_sync, dt);
            probe(err_s, var.sqrt());
            if err_s.abs() <= PROTECTION_K * var.sqrt() {
                contained += 1;
            }
        }
        series.push(Sample {
            t,
            error_ns: err_s * 1e9,
            gnss,
        });
    }
    let mut fom = score(&series, scn.threshold_ns);
    // Integrity: how reliably the filter's protection bound contains the truth.
    if outage_samples > 0 {
        fom.integrity = Some(contained as f64 / outage_samples as f64);
    }
    // Security: clock-aided spoof-detection score relative to the timing spec.
    fom.security = Some(spoof_detection_score(
        cfg.q_wf,
        cfg.q_rw,
        PHASE_MEAS_VAR_S2,
        scn.threshold_ns,
        SPOOF_MONITOR_S,
        dt,
        SPOOF_DETECT_K,
    ));
    // Filter-consistency health: a Monte-Carlo NIS/NEES check that the deployed
    // Kalman tuning (Q matched to the truth model, q_factor = 1) is self-consistent.
    let health_cfg = crate::filter_health::HealthConfig {
        q_wf: cfg.q_wf,
        q_rw: cfg.q_rw,
        r: PHASE_MEAS_VAR_S2,
        dt,
        steps: HEALTH_STEPS,
        seeds: HEALTH_SEEDS,
        q_factor: 1.0,
        base_seed: seed ^ HEALTH_SEED_SALT,
    };
    // With a flicker floor the truth the check draws includes the same flicker bank as the
    // clock above, and the filter is the matched extended filter; without one, the
    // two-state check as before.
    let filter_health = Some(match clock.flicker_bank() {
        Some(bank) => crate::filter_health::assess_with_flicker(
            crate::filter_health::HealthConfig {
                steps: HEALTH_STEPS_FLICKER,
                seeds: HEALTH_SEEDS_FLICKER,
                ..health_cfg
            },
            &bank,
        ),
        None => crate::filter_health::assess(health_cfg),
    });
    ClockRun {
        spec: clock.spec(),
        series,
        fom,
        adev_curve: crate::allan::overlapping_adev_curve(&phase, dt),
        filter_health,
    }
}

/// Run a clock-holdover scenario whose GNSS availability is derived from orbital
/// geometry. The user orbit, constellation, and elevation mask produce a
/// visibility timeline, which then drives the standard clock-holdover run.
pub fn run_orbit_clock(scn: &crate::orbit::OrbitClockScenario) -> Result<RunResult, String> {
    let sats = scn.all_satellites()?;
    let (user, sky) = scn.user_and_sky_for(&sats)?;
    let timeline = crate::orbit::build_timeline_at(
        &user,
        &sats,
        scn.time.step_s,
        scn.time.duration_s,
        scn.mask_deg,
        sky.as_ref(),
    );
    let inner = Scenario {
        seed: scn.seed,
        threshold_ns: scn.threshold_ns,
        runs: 1,
        time: scn.time.clone(),
        gnss: timeline,
        clock_quantum: scn.clock_quantum.clone(),
        clock_classical: scn.clock_classical.clone(),
        integrity: None,
    };
    let mut result = run(&inner);
    // Emit the propagated user track (km) so the playground can draw the 3D orbit.
    // This is the orbit pack's only extra output; non-orbit runs leave it `None`.
    result.eci_track = Some(sample_eci_track_km(
        &user,
        scn.time.step_s,
        scn.time.duration_s,
    ));
    Ok(result)
}

/// Run the clock-holdover scenario for both clocks and assemble the result.
pub fn run(scn: &Scenario) -> RunResult {
    RunResult {
        schema_version: crate::interchange::SCHEMA_VERSION.into(),
        engine_version: env!("CARGO_PKG_VERSION").into(),
        scenario_hash: crate::report::hash_scenario(scn),
        seed: scn.seed,
        threshold_ns: scn.threshold_ns,
        quantum: run_clock(scn, &scn.clock_quantum, scn.seed),
        classical: run_clock(
            scn,
            &scn.clock_classical,
            scn.seed.wrapping_add(0x9e3779b97f4a7c15),
        ),
        eci_track: None,
        meta: None,
        integrity_fom: None,
    }
}

/// Cap on the propagated `eci_track` length, so the orbit-pack JSON stays small
/// (one sample per ~2 min over a day ≈ 720 points). The track is decimated to at
/// most this many points by striding the time grid.
const MAX_ECI_TRACK: usize = 720;

/// Sample the user spacecraft's propagated ECI position over the scenario's time
/// grid (`step_s`, `duration_s`), in kilometres, decimated to at most
/// [`MAX_ECI_TRACK`] points. The magnitude of each sample is an independent
/// physical fact (e.g. a GPS orbit is ~26,560 km — IS-GPS-200), so it is a
/// non-circular oracle for the orbit visualisation.
fn sample_eci_track_km(
    user: &crate::orbit::UserPath,
    step_s: f64,
    duration_s: f64,
) -> Vec<[f64; 3]> {
    let n = (duration_s / step_s).round() as usize;
    // The full grid has n+1 samples (i = 0..=n). Stride by ceil((n+1)/cap) so the
    // decimated track keeps at most MAX_ECI_TRACK points while still spanning the
    // whole grid (the final sample is included explicitly below).
    let stride = (n + 1).div_ceil(MAX_ECI_TRACK).max(1);
    let sample = |i: usize| {
        let p = user.position_eci(i as f64 * step_s);
        [p[0] / 1000.0, p[1] / 1000.0, p[2] / 1000.0]
    };
    let mut track = Vec::new();
    let mut last_i = 0;
    let mut i = 0;
    while i <= n {
        track.push(sample(i));
        last_i = i;
        i += stride;
    }
    // Always include the final grid sample so the track spans the whole run, even
    // when n is not a multiple of the stride.
    if last_i != n {
        track.push(sample(n));
    }
    track
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::scenario::*;

    fn demo() -> Scenario {
        Scenario {
            seed: 7,
            threshold_ns: 100.0,
            runs: 1,
            time: TimeCfg {
                step_s: 10.0,
                duration_s: 3600.0,
            },
            gnss: GnssTimeline {
                windows: vec![
                    GnssWindow {
                        t0: 0.0,
                        t1: 600.0,
                        state: GnssState::Nominal,
                    },
                    GnssWindow {
                        t0: 600.0,
                        t1: 3600.0,
                        state: GnssState::Denied,
                    },
                ],
            },
            clock_quantum: ClockCfg {
                id: "optical".into(),
                provenance: "demo".into(),
                y0: 1e-13,
                q_wf: 1e-26,
                q_rw: 1e-34,
                drift: 0.0,
                flicker_floor: 0.0,
            },
            clock_classical: ClockCfg {
                id: "csac".into(),
                provenance: "demo".into(),
                y0: 1e-11,
                q_wf: 1e-24,
                q_rw: 1e-32,
                drift: 0.0,
                flicker_floor: 0.0,
            },
            integrity: None,
        }
    }

    #[test]
    fn nominal_window_has_zero_error() {
        let r = run(&demo());
        assert_eq!(r.quantum.series[0].error_ns, 0.0);
    }

    #[test]
    fn quantum_outperforms_classical() {
        let r = run(&demo());
        assert!(r.quantum.fom.timing_p95_ns < r.classical.fom.timing_p95_ns);
        assert!(r.quantum.fom.availability >= r.classical.fom.availability);
    }

    #[test]
    fn integrity_is_populated_and_bound_is_reliable() {
        // The Kalman protection bound, whose process noise matches the truth model,
        // should contain the actual error on the overwhelming majority of outage
        // samples (3-sigma => ~99.7% for a well-matched filter).
        let r = run(&demo());
        let qi = r
            .quantum
            .fom
            .integrity
            .expect("quantum integrity populated");
        let ci = r
            .classical
            .fom
            .integrity
            .expect("classical integrity populated");
        assert!(qi >= 0.95, "quantum integrity too low: {qi}");
        assert!(ci >= 0.95, "classical integrity too low: {ci}");
        assert!(qi <= 1.0 && ci <= 1.0);
    }

    #[test]
    fn security_is_populated_and_quantum_leads() {
        // Both clocks report a spoof-detection score; the quieter quantum clock
        // has a tighter detection floor and so scores at least as high.
        let r = run(&demo());
        let qs = r.quantum.fom.security.expect("quantum security populated");
        let cs = r
            .classical
            .fom
            .security
            .expect("classical security populated");
        assert!((0.0..=1.0).contains(&qs) && (0.0..=1.0).contains(&cs));
        assert!(qs >= cs, "quantum security {qs} < classical {cs}");
    }

    // An orbit scenario whose user spacecraft sits at the GPS orbital radius:
    // semi-major axis 26,559.7 km (IS-GPS-200 nominal) → altitude above the mean
    // Earth radius (6371 km) of 20,188.7 km. The constellation is irrelevant to
    // the user track, so a small synthetic Walker pattern keeps the test fast.
    const ORBIT_GPS_USER: &str = r#"
kind = "orbit"
seed = 7
threshold_ns = 10.0
mask_deg = 5.0

[time]
step_s = 120.0
duration_s = 43200.0

[user]
altitude_km = 20188.7
inclination_deg = 55.0
u0_deg = 0.0

[constellation]
altitude_km = 20180.0
inclination_deg = 55.0
planes = 3
sats_per_plane = 3
phasing_f = 1.0

[clock_quantum]
id = "optical"
provenance = "test"
y0 = 1.0e-15
q_wf = 1.0e-30
q_rw = 1.0e-40

[clock_classical]
id = "csac"
provenance = "test"
y0 = 1.0e-11
q_wf = 9.0e-20
q_rw = 1.0e-28
"#;

    #[test]
    fn orbit_run_emits_eci_track_with_grid_length_and_consistent_first_radius() {
        let scn: crate::orbit::OrbitClockScenario = toml::from_str(ORBIT_GPS_USER).unwrap();
        let r = run_orbit_clock(&scn).expect("orbit run");
        let track = r.eci_track.expect("orbit run emits eci_track");
        assert!(!track.is_empty(), "track is non-empty");

        // Length: the full grid has n+1 = duration/step + 1 = 361 samples; with a
        // ≤720-point cap and stride 1 it is emitted in full.
        let n = (scn.time.duration_s / scn.time.step_s).round() as usize;
        assert_eq!(track.len(), n + 1, "track length matches grid samples");

        // Engine-internal consistency: the first track sample's magnitude equals
        // |user.position_eci(0)| (in km).
        let p0 = scn.user.as_ref().unwrap().to_orbit().position_eci(0.0);
        let r0_km = (p0[0].ppowi(2) + p0[1].ppowi(2) + p0[2].ppowi(2)).sqrt() / 1000.0;
        let first = track[0];
        let track_r0 = (first[0].ppowi(2) + first[1].ppowi(2) + first[2].ppowi(2)).sqrt();
        assert!(
            (track_r0 - r0_km).abs() < 1e-6,
            "first track radius {track_r0} != user radius {r0_km}"
        );

        // External oracle (non-circular): the GPS nominal semi-major axis is
        // 26,559.7 km (IS-GPS-200 / published constellation parameter). The user
        // here is configured at that radius, so |r| must match within tolerance.
        assert!(
            (track_r0 - 26_559.7).abs() < 2_000.0,
            "GPS-altitude user radius {track_r0} km not ≈ 26,559.7 km"
        );
    }

    #[test]
    fn eci_track_is_decimated_below_the_cap() {
        // A day at a 60 s step would be 1441 grid samples; the cap (720) decimates.
        let mut scn: crate::orbit::OrbitClockScenario = toml::from_str(ORBIT_GPS_USER).unwrap();
        scn.time.step_s = 60.0;
        scn.time.duration_s = 86400.0;
        let r = run_orbit_clock(&scn).expect("orbit run");
        let track = r.eci_track.expect("eci_track");
        assert!(
            track.len() <= MAX_ECI_TRACK,
            "decimated track {} exceeds cap {MAX_ECI_TRACK}",
            track.len()
        );
        assert!(track.len() > 1, "decimated track keeps more than one point");
    }

    #[test]
    fn non_orbit_run_has_no_eci_track() {
        let r = run(&demo());
        assert!(r.eci_track.is_none(), "clock run carries no eci_track");
    }

    /// D11 / B4' (replacement for the withdrawn B4, fixed by the coordinator BEFORE any
    /// computation): the pooled RMS of `error / 1-sigma bound` over all outage samples and
    /// all runs lies in `[0.8, 1.2]` for EACH clock, on BOTH the shipped scenario (seed 42)
    /// and a fresh seed (7), 200 runs each.
    ///
    /// Basis, a priori: the flicker error within a run is nearly a constant random
    /// frequency, so each run contributes about one chi-square(1) draw; over 200 runs the mean
    /// square has a standard deviation of about sqrt(2/200) = 0.10 and the RMS about 0.05, so
    /// [0.8, 1.2] is about +-4 sigma. It still catches a bound inflated by more than 1.25x.
    /// B4 itself (coverage <= 0.9995) is withdrawn: it sits inside the metric's own sampling
    /// noise (std about 0.0037 against a 0.0022 margin below 1) and the exploratory prototype's
    /// 0.99968 would have failed it. If this test fails, report it; do not adjust the band.
    #[test]
    fn the_bound_is_not_inflated_pooled_rms_of_error_over_sigma() {
        let src = include_str!("../scenarios/clock-ensemble.toml");
        for seed in [42u64, 7] {
            let mut scn: Scenario = toml::from_str(src).expect("shipped scenario parses");
            scn.seed = seed;
            assert_eq!(scn.runs, 200);
            for (name, quantum) in [("quantum", true), ("classical", false)] {
                let cfg = if quantum {
                    &scn.clock_quantum
                } else {
                    &scn.clock_classical
                };
                let (mut sum_sq, mut n) = (0.0f64, 0u64);
                for k in 0..scn.runs as u64 {
                    let mut s = scn.seed.wrapping_add(k);
                    if !quantum {
                        s = s.wrapping_add(crate::ensemble::GOLDEN);
                    }
                    run_clock_probed(&scn, cfg, s, &mut |err, sigma| {
                        let z = err / sigma;
                        sum_sq += z * z;
                        n += 1;
                    });
                }
                assert!(n > 0);
                let rms = (sum_sq / n as f64).sqrt();
                eprintln!("B4' seed {seed} {name}: pooled RMS {rms:.4} over {n} samples");
                assert!(
                    (0.8..=1.2).contains(&rms),
                    "seed {seed}, {name}: pooled RMS(error / 1-sigma bound) = {rms:.4}, band [0.8, 1.2]"
                );
            }
        }
    }
}