Skip to main content

sidereon_core/estimation/
recipe.rs

1//! Named operation-order recipes for the GNSS estimation substrate.
2//!
3//! Phase-2 collapses the three thick estimator stacks (`spp`, `rtk`/`rtk_filter`,
4//! `precise_positioning`) onto one shared substrate plus thin, runtime-selectable
5//! strategies. The single hard constraint is that each external reference's
6//! bit-exactness (Skyfield for SPP, RTKLIB for RTK, the PPP oracle for PPP) must
7//! be preserved to 0 ULP. Different references need different floating-point
8//! operation orders for the *same* physical quantity, so the substrate never
9//! "simplifies" a parity-sensitive formula into one shared form. Instead every
10//! such choice is a NAMED variant: a strategy selects the op-order it needs by
11//! enum value rather than by owning a copy of the helper.
12//!
13//! This module *names* the recipes; the substrate and strategies route every
14//! caller through them. Each reference-faithful strategy resolves to the single
15//! op-order it was already using, so threading the recipe through the shared
16//! spine reproduces the prior code path bit-for-bit and leaves every existing
17//! 0-ULP golden unchanged.
18//!
19//! The `Canonical*` variants belong to the single consistent IERS-rigorous
20//! model (the bounded-tolerance canonical strategy, P6). They are NOT used by
21//! any reference-faithful strategy; canonical is an additional selectable
22//! strategy that changes nothing about the references. The SPP canonical range
23//! and frame variants ([`RangeRecipe::CanonicalLightTimeClosedFormSagnac`],
24//! [`FrameRecipe::CanonicalWgs84`]) are implemented and driven by
25//! [`EstimationRecipe::canonical_spp`]; the RTK and PPP canonical square-root
26//! solve ([`NormalRecipe::CanonicalSquareRoot`] on
27//! [`SolverRecipe::OwnedDeterministicCholesky`]) by
28//! [`EstimationRecipe::canonical_rtk`] and [`EstimationRecipe::canonical_ppp`].
29//! Canonical SPP, RTK, and PPP are all wired.
30
31/// Estimation technique: which physical observation model and parameter set a
32/// strategy estimates.
33#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Default)]
34pub enum Technique {
35    /// Single-point positioning: undifferenced pseudorange PVT.
36    #[default]
37    Spp,
38    /// Real-time kinematic: double-differenced code/phase baseline.
39    Rtk,
40    /// Precise point positioning: undifferenced ionosphere-free code/phase.
41    Ppp,
42}
43
44/// The reference a reference-faithful strategy is bit-exact against. The
45/// external oracles (Skyfield, RTKLIB, the PPP oracle) are CI validation targets
46/// whose goldens stay 0-ULP unchanged through P0-P5; [`Self::OwnedDeterministic`]
47/// is instead pinned to the owned solver's own frozen-bits golden (P5).
48#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Default)]
49pub enum ReferenceTarget {
50    /// Skyfield (the SPP geometry/clock/Sagnac reference).
51    #[default]
52    Skyfield,
53    /// RTKLIB (the RTK double-difference baseline reference).
54    Rtklib,
55    /// scipy least-squares host solve (the SPP solver-agreement reference).
56    /// Named for the `trust-region-least-squares` host-LAPACK fingerprint study;
57    /// not a runtime
58    /// estimation strategy (it is not wired into the SPP solve path), so it is
59    /// not a valid [`StrategyId`] target.
60    Scipy,
61    /// The PPP float/fixed oracle arc.
62    PppOracle,
63    /// The SPP owned deterministic trust-region solver
64    /// ([`SolverRecipe::OwnedDeterministicTrf`]): a fixed-reduction-order dense
65    /// subproblem factorization with no nalgebra LU and no black-box BLAS in that
66    /// solve. It is pinned to its own frozen-bits golden rather than to an
67    /// external library, and is selectable only for [`Technique::Spp`]. The owned
68    /// kernel owns only the subproblem factorization (the surrounding
69    /// normal-matrix / gradient / norm reductions stay on nalgebra), so its
70    /// cross-platform bit guarantee is scoped to the factorization; see
71    /// [`SolverRecipe::OwnedDeterministicTrf`] for the precise scope.
72    OwnedDeterministic,
73}
74
75/// Runtime-selectable strategy identity. `Reference` strategies are 0-ULP
76/// bit-exact to an external reference and remain the validation oracles;
77/// `Canonical` is the single bounded-tolerance "best" model (P6). Canonical SPP,
78/// RTK, and PPP are all wired.
79#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
80pub enum StrategyId {
81    /// A reference-faithful strategy: 0-ULP to `target` for `technique`.
82    Reference {
83        technique: Technique,
84        target: ReferenceTarget,
85    },
86    /// The canonical strategy for `technique` (bounded-tolerance, truth-gated).
87    Canonical { technique: Technique },
88}
89
90impl Default for StrategyId {
91    fn default() -> Self {
92        Self::Reference {
93            technique: Technique::Spp,
94            target: ReferenceTarget::Skyfield,
95        }
96    }
97}
98
99impl StrategyId {
100    /// SPP, bit-exact to Skyfield (`spp::solve`).
101    pub const fn spp_reference() -> Self {
102        Self::Reference {
103            technique: Technique::Spp,
104            target: ReferenceTarget::Skyfield,
105        }
106    }
107
108    /// RTK, bit-exact to RTKLIB (`rtk` / `rtk_filter`).
109    pub const fn rtk_reference() -> Self {
110        Self::Reference {
111            technique: Technique::Rtk,
112            target: ReferenceTarget::Rtklib,
113        }
114    }
115
116    /// PPP, bit-exact to the PPP oracle arc (`precise_positioning`).
117    pub const fn ppp_reference() -> Self {
118        Self::Reference {
119            technique: Technique::Ppp,
120            target: ReferenceTarget::PppOracle,
121        }
122    }
123
124    /// SPP via the owned deterministic trust-region solver
125    /// ([`SolverRecipe::OwnedDeterministicTrf`]): the owned dense subproblem
126    /// factorization, pinned to its own frozen-bits golden (its cross-platform
127    /// bit guarantee is scoped to the factorization). Selecting this through
128    /// [`crate::estimation::strategies::estimate`] drives the owned solver
129    /// rather than the legacy nalgebra LU path.
130    pub const fn spp_owned_deterministic() -> Self {
131        Self::Reference {
132            technique: Technique::Spp,
133            target: ReferenceTarget::OwnedDeterministic,
134        }
135    }
136}
137
138/// Geometric range / light-time / transmit-time operation order. Each variant
139/// names an existing range model; the substrate selects the op-order rather
140/// than copying the helper.
141#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Default)]
142pub enum RangeRecipe {
143    /// SPP closed-form light-time with a fixed transmit-time iteration count and
144    /// a measured-pseudorange seed (`spp/mod.rs` `sat_model`).
145    #[default]
146    SppMeasuredPseudorangeFixedIter,
147    /// `observables::predict` rounded-microsecond transmit time with a fixed
148    /// light-time iteration count (PPP / forward-prediction model).
149    ObservableRoundedMicrosecondFixedIter,
150    /// RTK provided-transmit-position range with the RTKLIB first-order Sagnac
151    /// scalar (`rtk_filter::model` line-of-sight / geometric range).
152    RtkProvidedTxFirstOrderSagnac,
153    /// Canonical: full iterative light-time (iterated to convergence, not a
154    /// fixed truncation) with the closed-form Sagnac Z-rotation, never a
155    /// first-order scalar Sagnac. Driven by [`EstimationRecipe::canonical_spp`]
156    /// in the SPP measurement model; not used by any reference strategy.
157    CanonicalLightTimeClosedFormSagnac,
158}
159
160/// Earth-rotation (Sagnac) correction operation order.
161#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Default)]
162pub enum SagnacRecipe {
163    /// Closed-form z-axis rotation of the satellite ECEF position by
164    /// `OMEGA_E_DOT * tau` (SPP / observables).
165    #[default]
166    ClosedFormZRotation,
167    /// RTKLIB first-order scalar Sagnac term added to the range
168    /// (`rtk_filter::model`).
169    RtklibFirstOrderScalar,
170    /// No Sagnac correction (synthetic / ECI-consistent inputs).
171    Off,
172}
173
174/// Local-frame / ENU / az-el basis construction operation order.
175#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Default)]
176pub enum FrameRecipe {
177    /// SPP Skyfield-parity ECEF->geodetic with the three-iteration AU-scaled
178    /// latitude solve (`spp` geodetic conversion).
179    #[default]
180    SppSkyfieldAuThreeIter,
181    /// Geocentric-up local frame used by the RTK elevation reference
182    /// (`rtk_filter` elevation mask / antenna projection).
183    GeocentricUpRtkReference,
184    /// Geodetic NEU basis built from the cross-product convention
185    /// (`precise_positioning::model` troposphere geometry).
186    GeodeticNeuCrossProduct,
187    /// DOP ENU rotation basis (`dop`).
188    DopEnuRotation,
189    /// Canonical: one consistent meters-native WGS84/ITRF geodetic basis under
190    /// IERS conventions (the core PROJ-pinned closed-form solve), never a
191    /// reference-specific AU-scaled path. Driven by
192    /// [`EstimationRecipe::canonical_spp`]; not used by any reference strategy.
193    CanonicalWgs84,
194}
195
196/// Normal-equation assembly tie-breaking / fold order. The tie order is the
197/// pivot/elimination convention that fixes the bit pattern of the reduced
198/// system.
199#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Default)]
200pub enum NormalRecipe {
201    /// SPP weighted-residual rows with a finite-difference design matrix
202    /// (`spp` least-squares).
203    #[default]
204    SppWeightedResidualFiniteDifference,
205    /// RTK double-difference block assembly with the first-tie covariance fold
206    /// (`rtk_filter::normal` first-tie block).
207    RtkDoubleDifferenceBlockFirstTie,
208    /// PPP weighted normal equations with epoch-local receiver clocks eliminated
209    /// and the reduced static system solved last-tie.
210    PppDenseLastTie,
211    /// Canonical square-root-information solve, shared by canonical RTK and
212    /// canonical PPP: the SPD normal system is solved by the owned deterministic
213    /// Cholesky factorization `Λ = L Lᵀ` plus forward/back substitution, where
214    /// `L` is the information-matrix square root. For RTK this is the
215    /// double-difference information system `Λ x = η` assembled by the same shared
216    /// block fold the RTK reference uses; for PPP it is the weighted normal
217    /// system assembled from the same undifferenced rows the PPP reference uses,
218    /// with epoch-local receiver clocks eliminated before factorization. This is
219    /// the numerically rigorous op-order for an SPD normal
220    /// matrix (no pivoting; exploits symmetry), distinct from the reference RTK
221    /// general first-tie Gaussian elimination
222    /// ([`Self::RtkDoubleDifferenceBlockFirstTie`]) and the reference PPP last-tie
223    /// Gaussian elimination ([`Self::PppDenseLastTie`]). Driven by
224    /// [`EstimationRecipe::canonical_rtk`] and [`EstimationRecipe::canonical_ppp`]
225    /// on the owned [`SolverRecipe::OwnedDeterministicCholesky`] kernel; not used
226    /// by any reference strategy.
227    CanonicalSquareRoot,
228}
229
230/// Linear-solve / factorization operation order. Determinism note: the legacy
231/// SPP path is nalgebra LU (not bit-portable end-to-end), preserved as a named
232/// variant; the owned deterministic kernel (P5) owns the dense subproblem
233/// factorization with its own goldens. Its determinism scope is the
234/// factorization, not the surrounding nalgebra reductions that build the
235/// subproblem -- see [`Self::OwnedDeterministicTrf`].
236#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Default)]
237pub enum SolverRecipe {
238    /// nalgebra trust-region least squares, the current SPP solver
239    /// (`spp` / `crate::astro::math::least_squares`). Existing SPP goldens use
240    /// this; kept unchanged.
241    #[default]
242    NalgebraTrfLegacy,
243    /// Flat first-tie Gaussian elimination (RTK baseline/filter solve).
244    FlatGaussianFirstTie,
245    /// Dense last-tie Gaussian elimination (PPP solve,
246    /// `crate::astro::math::linear::solve_linear_last_tie`).
247    DenseGaussianLastTie,
248    /// scipy host LAPACK reference solve (machine-dependent; only as a
249    /// fingerprinted CI reference, never canonical).
250    ScipyHostLapackReference,
251    /// Owned deterministic Cholesky (square-root) linear solve, the canonical RTK
252    /// (P6 increment 2) and canonical PPP (P6 increment 3) solver: the SPD normal
253    /// system is factored `Λ = L Lᵀ` and solved by forward/back substitution
254    /// through the owned
255    /// [`crate::astro::math::linear::solve_flat_normal_square_root_into`] kernel,
256    /// with no nalgebra LU and no black-box BLAS. Paired with
257    /// [`NormalRecipe::CanonicalSquareRoot`]. Both the RTK and PPP canonical paths
258    /// are owned scalar arithmetic and f64 sqrt is IEEE-754 correctly rounded, so
259    /// unlike [`Self::OwnedDeterministicTrf`] (whose surrounding reductions ride
260    /// nalgebra) its bit guarantee covers the full solve and is portable across
261    /// platforms.
262    OwnedDeterministicCholesky,
263    /// Owned deterministic trust-region subproblem solve added in P5: a
264    /// fixed-reduction-order dense Gaussian elimination (the
265    /// `OwnedGaussianFirstTie` kernel) with no nalgebra LU and no black-box BLAS
266    /// in the factorization, pinned to its OWN frozen-bits goldens. Scope: it
267    /// owns ONLY the subproblem factorization; the normal-matrix / gradient /
268    /// norm reductions that build the subproblem still flow through nalgebra's
269    /// CPU-dispatched dense algebra, so the cross-platform bit guarantee is
270    /// scoped to the factorization, not the full solve.
271    OwnedDeterministicTrf,
272}
273
274/// The full operation-order recipe a strategy composes: one variant per stage.
275/// `Default` and the named constructors reproduce the CURRENT behavior of each
276/// existing strategy, so selecting a recipe never changes a reference golden
277/// (PPP goldens were re-frozen once when static PPP moved to clock-eliminated
278/// reduced normals; see `estimation::strategies`).
279#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Default)]
280pub struct EstimationRecipe {
281    pub range: RangeRecipe,
282    pub sagnac: SagnacRecipe,
283    pub frame: FrameRecipe,
284    pub normal: NormalRecipe,
285    pub solver: SolverRecipe,
286}
287
288impl EstimationRecipe {
289    /// The current SPP reference recipe (`spp::solve`, Skyfield-parity).
290    pub const fn spp() -> Self {
291        Self {
292            range: RangeRecipe::SppMeasuredPseudorangeFixedIter,
293            sagnac: SagnacRecipe::ClosedFormZRotation,
294            frame: FrameRecipe::SppSkyfieldAuThreeIter,
295            normal: NormalRecipe::SppWeightedResidualFiniteDifference,
296            solver: SolverRecipe::NalgebraTrfLegacy,
297        }
298    }
299
300    /// The current RTK reference recipe (`rtk` / `rtk_filter`, RTKLIB-parity).
301    pub const fn rtk() -> Self {
302        Self {
303            range: RangeRecipe::RtkProvidedTxFirstOrderSagnac,
304            sagnac: SagnacRecipe::RtklibFirstOrderScalar,
305            frame: FrameRecipe::GeocentricUpRtkReference,
306            normal: NormalRecipe::RtkDoubleDifferenceBlockFirstTie,
307            solver: SolverRecipe::FlatGaussianFirstTie,
308        }
309    }
310
311    /// The current PPP reference recipe (`precise_positioning`, oracle-parity).
312    pub const fn ppp() -> Self {
313        Self {
314            range: RangeRecipe::ObservableRoundedMicrosecondFixedIter,
315            sagnac: SagnacRecipe::ClosedFormZRotation,
316            frame: FrameRecipe::GeodeticNeuCrossProduct,
317            normal: NormalRecipe::PppDenseLastTie,
318            solver: SolverRecipe::DenseGaussianLastTie,
319        }
320    }
321
322    /// The SPP recipe driving the owned deterministic trust-region solver: the
323    /// SPP reference model with [`SolverRecipe::OwnedDeterministicTrf`] swapped
324    /// in for the legacy nalgebra LU linear-solve stage. Every other stage is the
325    /// SPP reference op-order, so only the factorization changes.
326    pub const fn spp_owned_deterministic() -> Self {
327        let mut recipe = Self::spp();
328        recipe.solver = SolverRecipe::OwnedDeterministicTrf;
329        recipe
330    }
331
332    /// The canonical SPP recipe: the single consistent IERS-rigorous SPP
333    /// measurement model. It diverges from [`Self::spp`] (the Skyfield-faithful
334    /// reference) only where the physics says to:
335    /// - range: [`RangeRecipe::CanonicalLightTimeClosedFormSagnac`] iterates the
336    ///   light-time loop to convergence (vs the reference's fixed
337    ///   transmit-time truncation), with the closed-form Sagnac Z-rotation (never
338    ///   a first-order scalar Sagnac).
339    /// - frame: [`FrameRecipe::CanonicalWgs84`] solves ECEF->geodetic directly in
340    ///   meters on the WGS84 ellipsoid (vs the reference's Skyfield AU-scaled
341    ///   three-iteration latitude loop).
342    /// - solver: [`SolverRecipe::OwnedDeterministicTrf`] owns the trust-region
343    ///   subproblem factorization so canonical is deterministic run-to-run on a
344    ///   pinned build (its cross-platform bit guarantee is scoped to the
345    ///   factorization; the surrounding reductions ride nalgebra).
346    ///
347    /// The Sagnac stage is the closed-form Z-rotation the SPP reference already
348    /// uses (the rigorous form), and the normal stage is the SPP
349    /// weighted-residual finite-difference assembly the trust-region solver
350    /// consumes; neither needs a separate canonical variant for SPP.
351    pub const fn canonical_spp() -> Self {
352        Self {
353            range: RangeRecipe::CanonicalLightTimeClosedFormSagnac,
354            sagnac: SagnacRecipe::ClosedFormZRotation,
355            frame: FrameRecipe::CanonicalWgs84,
356            normal: NormalRecipe::SppWeightedResidualFiniteDifference,
357            solver: SolverRecipe::OwnedDeterministicTrf,
358        }
359    }
360
361    /// The canonical RTK recipe: the double-difference baseline under the
362    /// numerically rigorous square-root-information solve. It keeps the RTK
363    /// reference's double-difference measurement physics (the provided-transmit
364    /// range with the RTKLIB first-order Sagnac scalar, the geocentric-up
365    /// elevation frame), because the canonical RTK divergence the physics calls
366    /// for is in the linear algebra, not the observation model: the same SPD
367    /// information system the reference assembles is solved by the owned
368    /// deterministic Cholesky square-root factorization
369    /// ([`NormalRecipe::CanonicalSquareRoot`] on
370    /// [`SolverRecipe::OwnedDeterministicCholesky`]) instead of the reference's
371    /// general first-tie Gaussian elimination. The square-root solve needs no
372    /// pivoting, exploits the symmetry of the SPD normal matrix, and is entirely
373    /// owned scalar arithmetic (no nalgebra, no BLAS), so canonical RTK is
374    /// well-conditioned and bit-reproducible across platforms.
375    pub const fn canonical_rtk() -> Self {
376        Self {
377            range: RangeRecipe::RtkProvidedTxFirstOrderSagnac,
378            sagnac: SagnacRecipe::RtklibFirstOrderScalar,
379            frame: FrameRecipe::GeocentricUpRtkReference,
380            normal: NormalRecipe::CanonicalSquareRoot,
381            solver: SolverRecipe::OwnedDeterministicCholesky,
382        }
383    }
384
385    /// The canonical PPP recipe: the undifferenced ionosphere-free PPP arc under
386    /// the numerically rigorous square-root-information solve. Like
387    /// [`Self::canonical_rtk`] it keeps the PPP reference's measurement physics
388    /// (the rounded-microsecond fixed-iteration light-time with the rigorous
389    /// closed-form Sagnac Z-rotation, and the geodetic NEU antenna frame), because
390    /// the canonical PPP divergence the physics calls for is in the linear
391    /// algebra, not the observation model: the same weighted normal equations
392    /// the reference assembles from the undifferenced rows are reduced by
393    /// eliminating epoch-local receiver clocks, then solved by the owned
394    /// deterministic Cholesky square-root factorization
395    /// ([`NormalRecipe::CanonicalSquareRoot`] on
396    /// [`SolverRecipe::OwnedDeterministicCholesky`]) instead of the reference's
397    /// dense last-tie Gaussian elimination ([`SolverRecipe::DenseGaussianLastTie`]).
398    /// The square-root solve needs no pivoting, exploits the symmetry of the SPD
399    /// normal matrix, and is entirely owned scalar arithmetic (no nalgebra, no
400    /// BLAS), so it is well-conditioned and the solve itself is bit-portable.
401    /// Determinism scope (calibrated, not overstated): unlike canonical RTK, the PPP
402    /// measurement model that builds the rows evaluates troposphere mapping,
403    /// antenna, and geodetic-frame transcendentals through the platform math
404    /// library, so canonical PPP's overall output is bit-reproducible run-to-run on
405    /// a pinned build but is not claimed bit-portable across platforms; only the
406    /// owned Cholesky solve carries the cross-platform guarantee.
407    pub const fn canonical_ppp() -> Self {
408        Self {
409            range: RangeRecipe::ObservableRoundedMicrosecondFixedIter,
410            sagnac: SagnacRecipe::ClosedFormZRotation,
411            frame: FrameRecipe::GeodeticNeuCrossProduct,
412            normal: NormalRecipe::CanonicalSquareRoot,
413            solver: SolverRecipe::OwnedDeterministicCholesky,
414        }
415    }
416
417    /// The canonical recipe for a `technique`. Canonical SPP (P6 increment 1),
418    /// canonical RTK (P6 increment 2), and canonical PPP (P6 increment 3) are all
419    /// wired, so every technique has a canonical strategy. Returns `Option` to keep
420    /// the resolver's "not yet implemented" surface stable.
421    pub const fn for_canonical(technique: Technique) -> Option<Self> {
422        match technique {
423            Technique::Spp => Some(Self::canonical_spp()),
424            Technique::Rtk => Some(Self::canonical_rtk()),
425            Technique::Ppp => Some(Self::canonical_ppp()),
426        }
427    }
428
429    /// The reference recipe for an explicit `(technique, target)` pair, or `None`
430    /// if the pair is not a supported reference strategy. This is the single
431    /// source of truth for which targets each technique can run: only the wired
432    /// reference oracles (Skyfield for SPP, RTKLIB for RTK, the PPP oracle for
433    /// PPP) and the SPP owned deterministic solver are valid. Every other pair
434    /// (a cross-technique oracle, or the unwired scipy host-LAPACK reference) is
435    /// rejected so an impossible strategy can never silently run a mismatched
436    /// recipe.
437    pub const fn for_reference(technique: Technique, target: ReferenceTarget) -> Option<Self> {
438        match (technique, target) {
439            (Technique::Spp, ReferenceTarget::Skyfield) => Some(Self::spp()),
440            (Technique::Spp, ReferenceTarget::OwnedDeterministic) => {
441                Some(Self::spp_owned_deterministic())
442            }
443            (Technique::Rtk, ReferenceTarget::Rtklib) => Some(Self::rtk()),
444            (Technique::Ppp, ReferenceTarget::PppOracle) => Some(Self::ppp()),
445            _ => None,
446        }
447    }
448}
449
450/// How a strategy forms its integer-ambiguity identifiers, and against what they
451/// are referenced. Naming this lets the RTK and PPP fixed solvers share one
452/// LAMBDA resolution kernel
453/// ([`crate::estimation::substrate::ambiguity::resolve_integer_lattice`]) and
454/// differ only in DATA rather than in separate algorithm trees.
455#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
456pub enum DifferencingMode {
457    /// Double-differenced ambiguities, one reference satellite per constellation
458    /// (the RTK baseline / sequential-filter convention: each non-reference
459    /// satellite is differenced against its own system's reference).
460    DoubleDifferencePerSystemReference,
461    /// Undifferenced ambiguities, one per satellite per receiver (the PPP
462    /// convention: no reference satellite, all satellites carry their own
463    /// ionosphere-free ambiguity).
464    Undifferenced,
465}
466
467/// Whether partial ambiguity resolution is attempted when the full-set integer
468/// fix fails its ratio test, and with what floor on the retained subset size.
469#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
470pub enum PartialResolution {
471    /// Full-set only: a failed ratio test means "not fixed" (PPP, and the RTK
472    /// sequential filter, both take the full set or nothing).
473    Disabled,
474    /// Confidence-ranked then exhaustive subset fallback down to
475    /// `min_ambiguities` retained (the RTK static fixed solver,
476    /// `rtk_filter::search::search_partial_fixed_ambiguities`).
477    Exhaustive { min_ambiguities: usize },
478}
479
480/// The integer-ambiguity identity/eligibility policy a strategy resolves under:
481/// the strategy DATA that replaces the RTK-vs-PPP algorithm-tree split. The
482/// LAMBDA resolution kernel is common; only these fields differ between the
483/// reference strategies. Named in P3; consumed by the runtime selector in P4.
484#[derive(Debug, Clone, Copy, PartialEq)]
485pub struct AmbiguityIdPolicy {
486    pub differencing: DifferencingMode,
487    /// Exclude float-only constellations from the integer search set.
488    pub float_only_gating: bool,
489    pub partial: PartialResolution,
490    /// Ratio-test acceptance threshold passed to the LAMBDA kernel.
491    pub ratio_threshold: f64,
492}
493
494impl AmbiguityIdPolicy {
495    /// The static RTK fixed-baseline policy (`rtk_filter::fixed`): per-system
496    /// double differences, float-only constellations excluded from the search,
497    /// partial resolution down to `partial_min_ambiguities`.
498    pub const fn rtk_static(ratio_threshold: f64, partial_min_ambiguities: usize) -> Self {
499        Self {
500            differencing: DifferencingMode::DoubleDifferencePerSystemReference,
501            float_only_gating: true,
502            partial: PartialResolution::Exhaustive {
503                min_ambiguities: partial_min_ambiguities,
504            },
505            ratio_threshold,
506        }
507    }
508
509    /// The sequential RTK filter policy (`rtk_filter::update`): per-system double
510    /// differences, float-only constellations excluded, full set or nothing.
511    pub const fn rtk_sequential(ratio_threshold: f64) -> Self {
512        Self {
513            differencing: DifferencingMode::DoubleDifferencePerSystemReference,
514            float_only_gating: true,
515            partial: PartialResolution::Disabled,
516            ratio_threshold,
517        }
518    }
519
520    /// The static PPP fixed policy (`precise_positioning::fixed`): undifferenced
521    /// per-satellite ambiguities, no constellation gating, full set or nothing.
522    pub const fn ppp(ratio_threshold: f64) -> Self {
523        Self {
524            differencing: DifferencingMode::Undifferenced,
525            float_only_gating: false,
526            partial: PartialResolution::Disabled,
527            ratio_threshold,
528        }
529    }
530}
531
532/// The operation order used to normalize one residual against its weight before
533/// the sigma comparison in a per-residual screen. Naming the order keeps each
534/// screen bit-identical while the formula lives in exactly one place
535/// ([`crate::estimation::substrate::qc::normalized_residual`]).
536#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
537pub enum ResidualNormRecipe {
538    /// `value · weight` where `weight` is an inverse double-difference *variance*
539    /// (`1/(sigma_sat^2 + sigma_ref^2)`), so the normalized innovation is
540    /// `value / sigma^2`. The RTK sequential information filter, whose DD rows
541    /// weight by inverse variance, screens its predicted innovations this way.
542    RtkInverseVarianceInnovation,
543    /// `value · weight` where `weight` is an inverse *sigma*
544    /// (`1/sqrt(sigma_sat^2 + sigma_ref^2)`), so the normalized residual is
545    /// `value / sigma`. The RTK static float/fixed least-squares baselines, whose
546    /// DD rows weight by inverse sigma, screen their post-fit residuals this way.
547    RtkInverseSigmaResidual,
548    /// `|value| · sqrt(weight)` where `weight` is an inverse *sigma* (`1/sigma`):
549    /// the residual magnitude times the square root of the inverse-sigma weight.
550    /// The PPP float leave-one-out screen (PPP rows weight by inverse sigma, as
551    /// `MeasurementWeights` documents).
552    PppInverseSigmaMagnitude,
553}
554
555/// The residual-screen family a strategy applies after (or, for the filter,
556/// before) a solve. Strategy DATA for the P4 selector; the chi-square variant is
557/// the SPP RAIM aggregate test, the rest are per-residual sigma screens that
558/// share [`crate::estimation::substrate::qc::normalized_residual`].
559#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
560pub enum ScreenKind {
561    /// SPP RAIM: aggregate chi-square on the weighted residual sum, then FDE
562    /// leave-one-out exclusion (`quality::raim` / `quality::fde`).
563    RaimChiSquare,
564    /// RTK static fixed: worst information-weighted residual vs a sigma gate,
565    /// excluding the worst satellite within a budget
566    /// (`rtk_filter::fixed::solve_fixed_baseline_validated`).
567    RtkFixedResidualValidation,
568    /// RTK sequential filter: information-weighted innovation screen on predicted
569    /// DD rows, masking rejected rows and coasting (`rtk_filter::update`).
570    RtkSequentialInnovation,
571    /// PPP float: worst studentized residual vs a sigma gate, leave-one-out prune
572    /// and re-solve while WRMS improves (`precise_positioning::float`).
573    PppFloatLeaveOneOut,
574}
575
576impl ScreenKind {
577    /// The per-residual normalization op-order this screen uses, or `None` for
578    /// the aggregate chi-square RAIM screen (which scores the weighted residual
579    /// sum, not individual residuals).
580    pub const fn residual_norm(self) -> Option<ResidualNormRecipe> {
581        match self {
582            Self::RaimChiSquare => None,
583            Self::RtkFixedResidualValidation => Some(ResidualNormRecipe::RtkInverseSigmaResidual),
584            Self::RtkSequentialInnovation => Some(ResidualNormRecipe::RtkInverseVarianceInnovation),
585            Self::PppFloatLeaveOneOut => Some(ResidualNormRecipe::PppInverseSigmaMagnitude),
586        }
587    }
588}
589
590#[cfg(test)]
591mod tests {
592    use super::*;
593
594    #[test]
595    fn defaults_name_current_spp_behavior() {
596        // The per-stage defaults are the SPP reference op-orders, so an
597        // unspecified recipe reproduces the current SPP path.
598        assert_eq!(EstimationRecipe::default(), EstimationRecipe::spp());
599        assert_eq!(
600            RangeRecipe::default(),
601            RangeRecipe::SppMeasuredPseudorangeFixedIter
602        );
603        assert_eq!(SagnacRecipe::default(), SagnacRecipe::ClosedFormZRotation);
604        assert_eq!(FrameRecipe::default(), FrameRecipe::SppSkyfieldAuThreeIter);
605        assert_eq!(
606            NormalRecipe::default(),
607            NormalRecipe::SppWeightedResidualFiniteDifference
608        );
609        assert_eq!(SolverRecipe::default(), SolverRecipe::NalgebraTrfLegacy);
610        assert_eq!(StrategyId::default(), StrategyId::spp_reference());
611    }
612
613    #[test]
614    fn strategy_constructors_match_reference_targets() {
615        assert_eq!(
616            StrategyId::spp_reference(),
617            StrategyId::Reference {
618                technique: Technique::Spp,
619                target: ReferenceTarget::Skyfield,
620            }
621        );
622        assert_eq!(
623            StrategyId::rtk_reference(),
624            StrategyId::Reference {
625                technique: Technique::Rtk,
626                target: ReferenceTarget::Rtklib,
627            }
628        );
629        assert_eq!(
630            StrategyId::ppp_reference(),
631            StrategyId::Reference {
632                technique: Technique::Ppp,
633                target: ReferenceTarget::PppOracle,
634            }
635        );
636    }
637
638    #[test]
639    fn for_reference_selects_each_supported_pairs_recipe() {
640        assert_eq!(
641            EstimationRecipe::for_reference(Technique::Spp, ReferenceTarget::Skyfield),
642            Some(EstimationRecipe::spp())
643        );
644        assert_eq!(
645            EstimationRecipe::for_reference(Technique::Rtk, ReferenceTarget::Rtklib),
646            Some(EstimationRecipe::rtk())
647        );
648        assert_eq!(
649            EstimationRecipe::for_reference(Technique::Ppp, ReferenceTarget::PppOracle),
650            Some(EstimationRecipe::ppp())
651        );
652    }
653
654    #[test]
655    fn owned_deterministic_recipe_swaps_only_the_solver() {
656        let owned = EstimationRecipe::spp_owned_deterministic();
657        assert_eq!(owned.solver, SolverRecipe::OwnedDeterministicTrf);
658        // Every non-solver stage is the SPP reference op-order.
659        assert_eq!(
660            EstimationRecipe {
661                solver: SolverRecipe::NalgebraTrfLegacy,
662                ..owned
663            },
664            EstimationRecipe::spp()
665        );
666        assert_eq!(
667            EstimationRecipe::for_reference(Technique::Spp, ReferenceTarget::OwnedDeterministic),
668            Some(owned)
669        );
670    }
671
672    #[test]
673    fn canonical_spp_recipe_uses_the_rigorous_op_orders() {
674        let canonical = EstimationRecipe::canonical_spp();
675        // Range: full iterative light-time with closed-form Sagnac, not the SPP
676        // reference's fixed-iteration measured-pseudorange recipe.
677        assert_eq!(
678            canonical.range,
679            RangeRecipe::CanonicalLightTimeClosedFormSagnac
680        );
681        assert_ne!(canonical.range, EstimationRecipe::spp().range);
682        // Frame: one consistent meters-native WGS84 basis, not the Skyfield AU
683        // path.
684        assert_eq!(canonical.frame, FrameRecipe::CanonicalWgs84);
685        assert_ne!(canonical.frame, EstimationRecipe::spp().frame);
686        // Sagnac stays the closed-form Z-rotation (the rigorous form the SPP
687        // reference already uses); the canonical divergence is never a
688        // first-order scalar Sagnac.
689        assert_eq!(canonical.sagnac, SagnacRecipe::ClosedFormZRotation);
690        assert_ne!(canonical.sagnac, SagnacRecipe::RtklibFirstOrderScalar);
691        // Solver: the owned deterministic factorization, for run-to-run
692        // determinism on a pinned build.
693        assert_eq!(canonical.solver, SolverRecipe::OwnedDeterministicTrf);
694        assert_eq!(
695            EstimationRecipe::for_canonical(Technique::Spp),
696            Some(canonical)
697        );
698    }
699
700    #[test]
701    fn canonical_rtk_recipe_uses_the_square_root_solve() {
702        let canonical = EstimationRecipe::canonical_rtk();
703        // Normal + solver: the owned Cholesky square-root information solve, not
704        // the reference RTK first-tie Gaussian elimination.
705        assert_eq!(canonical.normal, NormalRecipe::CanonicalSquareRoot);
706        assert_eq!(canonical.solver, SolverRecipe::OwnedDeterministicCholesky);
707        assert_ne!(canonical.normal, EstimationRecipe::rtk().normal);
708        assert_ne!(canonical.solver, EstimationRecipe::rtk().solver);
709        // Measurement physics stays the RTK reference double-difference model: the
710        // canonical RTK divergence is in the linear algebra, not the observation
711        // model, so range/sagnac/frame match the reference.
712        assert_eq!(canonical.range, EstimationRecipe::rtk().range);
713        assert_eq!(canonical.sagnac, EstimationRecipe::rtk().sagnac);
714        assert_eq!(canonical.frame, EstimationRecipe::rtk().frame);
715        assert_eq!(
716            EstimationRecipe::for_canonical(Technique::Rtk),
717            Some(canonical)
718        );
719    }
720
721    #[test]
722    fn canonical_ppp_recipe_uses_the_square_root_solve() {
723        let canonical = EstimationRecipe::canonical_ppp();
724        // Normal + solver: the owned Cholesky square-root information solve, not
725        // the reference PPP dense last-tie Gaussian elimination.
726        assert_eq!(canonical.normal, NormalRecipe::CanonicalSquareRoot);
727        assert_eq!(canonical.solver, SolverRecipe::OwnedDeterministicCholesky);
728        assert_ne!(canonical.normal, EstimationRecipe::ppp().normal);
729        assert_ne!(canonical.solver, EstimationRecipe::ppp().solver);
730        // Measurement physics stays the PPP reference undifferenced model: the
731        // canonical PPP divergence is in the linear algebra, not the observation
732        // model, so range/sagnac/frame match the reference.
733        assert_eq!(canonical.range, EstimationRecipe::ppp().range);
734        assert_eq!(canonical.sagnac, EstimationRecipe::ppp().sagnac);
735        assert_eq!(canonical.frame, EstimationRecipe::ppp().frame);
736        // Canonical RTK and PPP share the square-root normal + owned Cholesky
737        // solver (the same numerically rigorous SPD op-order).
738        assert_eq!(canonical.normal, EstimationRecipe::canonical_rtk().normal);
739        assert_eq!(canonical.solver, EstimationRecipe::canonical_rtk().solver);
740        assert_eq!(
741            EstimationRecipe::for_canonical(Technique::Ppp),
742            Some(canonical)
743        );
744    }
745
746    #[test]
747    fn for_canonical_wires_all_three_techniques() {
748        assert_eq!(
749            EstimationRecipe::for_canonical(Technique::Spp),
750            Some(EstimationRecipe::canonical_spp())
751        );
752        assert_eq!(
753            EstimationRecipe::for_canonical(Technique::Rtk),
754            Some(EstimationRecipe::canonical_rtk())
755        );
756        assert_eq!(
757            EstimationRecipe::for_canonical(Technique::Ppp),
758            Some(EstimationRecipe::canonical_ppp())
759        );
760    }
761
762    #[test]
763    fn for_reference_rejects_impossible_pairs() {
764        // Cross-technique oracles and the unwired scipy reference are not
765        // supported reference strategies.
766        for (technique, target) in [
767            (Technique::Spp, ReferenceTarget::Rtklib),
768            (Technique::Spp, ReferenceTarget::PppOracle),
769            (Technique::Spp, ReferenceTarget::Scipy),
770            (Technique::Rtk, ReferenceTarget::Skyfield),
771            (Technique::Rtk, ReferenceTarget::OwnedDeterministic),
772            (Technique::Rtk, ReferenceTarget::PppOracle),
773            (Technique::Ppp, ReferenceTarget::Skyfield),
774            (Technique::Ppp, ReferenceTarget::OwnedDeterministic),
775        ] {
776            assert_eq!(
777                EstimationRecipe::for_reference(technique, target),
778                None,
779                "{technique:?} + {target:?} must be rejected"
780            );
781        }
782    }
783
784    #[test]
785    fn reference_ambiguity_policies_name_current_behavior() {
786        let rtk_static = AmbiguityIdPolicy::rtk_static(3.0, 4);
787        assert_eq!(
788            rtk_static.differencing,
789            DifferencingMode::DoubleDifferencePerSystemReference
790        );
791        assert!(rtk_static.float_only_gating);
792        assert_eq!(
793            rtk_static.partial,
794            PartialResolution::Exhaustive { min_ambiguities: 4 }
795        );
796
797        let rtk_seq = AmbiguityIdPolicy::rtk_sequential(3.0);
798        assert_eq!(
799            rtk_seq.differencing,
800            DifferencingMode::DoubleDifferencePerSystemReference
801        );
802        assert!(rtk_seq.float_only_gating);
803        assert_eq!(rtk_seq.partial, PartialResolution::Disabled);
804
805        let ppp = AmbiguityIdPolicy::ppp(2.5);
806        assert_eq!(ppp.differencing, DifferencingMode::Undifferenced);
807        assert!(!ppp.float_only_gating);
808        assert_eq!(ppp.partial, PartialResolution::Disabled);
809    }
810
811    #[test]
812    fn rtk_and_ppp_id_policies_differ_only_in_data() {
813        // Same LAMBDA kernel, different identity/eligibility data: the two stacks
814        // are no longer separate algorithm trees, only different policy values.
815        let rtk = AmbiguityIdPolicy::rtk_static(3.0, 1);
816        let ppp = AmbiguityIdPolicy::ppp(3.0);
817        assert_ne!(rtk.differencing, ppp.differencing);
818        assert_ne!(rtk.float_only_gating, ppp.float_only_gating);
819        assert_ne!(rtk.partial, ppp.partial);
820    }
821
822    #[test]
823    fn screen_kinds_select_their_normalization_order() {
824        assert_eq!(ScreenKind::RaimChiSquare.residual_norm(), None);
825        assert_eq!(
826            ScreenKind::RtkFixedResidualValidation.residual_norm(),
827            Some(ResidualNormRecipe::RtkInverseSigmaResidual)
828        );
829        assert_eq!(
830            ScreenKind::RtkSequentialInnovation.residual_norm(),
831            Some(ResidualNormRecipe::RtkInverseVarianceInnovation)
832        );
833        assert_eq!(
834            ScreenKind::PppFloatLeaveOneOut.residual_norm(),
835            Some(ResidualNormRecipe::PppInverseSigmaMagnitude)
836        );
837    }
838
839    #[test]
840    fn each_strategy_selects_a_distinct_solver_order() {
841        // The three reference strategies must not collapse onto one solver
842        // op-order; that distinction is what preserves their separate goldens.
843        assert_ne!(
844            EstimationRecipe::spp().solver,
845            EstimationRecipe::rtk().solver
846        );
847        assert_ne!(
848            EstimationRecipe::rtk().solver,
849            EstimationRecipe::ppp().solver
850        );
851        assert_ne!(
852            EstimationRecipe::spp().solver,
853            EstimationRecipe::ppp().solver
854        );
855    }
856}