Merge pull request #398 from nyx-space/feat/gh-396-moar-serde

Orbit determination constant biases + measurement ambiguity and moduli

Author: ChristopherRabotin

Cargo.toml

@@ -1,7 +1,7 @@
 [package]
 name = "nyx-space"
 build = "build.rs"
-version = "2.0.0"
+version = "2.0.1"
 edition = "2021"
 authors = ["Christopher Rabotin <[email protected]>"]
 description = "A high-fidelity space mission toolkit, with orbit propagation, estimation and some systems engineering"

README.md

@@ -62,6 +62,12 @@ For Python projects, get started by installing the library via `pip`: `pip insta
 
 Nyx is provided under the [AGPLv3 License](./LICENSE). By using this software, you assume responsibility for adhering to the license. Refer to [the pricing page](https://nyxspace.com/pricing/?utm_source=readme-price) for an FAQ on the AGPLv3 license. Notably, any software that incorporates, links to, or depends on Nyx must also be released under the AGPLv3 license, even if you distribute an unmodified version of Nyx.
 
+# Versioning
+
+Nyx mostly adheres to SemVer. New patch versions should be rare. Updated dependencies trigger a new minor version. _However_ new fields in structures and new behavior may also be added with minor releases, but the public facing initializers and functions should not significantly change (but may still change).
+
+Major releases are for dramatic changes.
+
 
 [cratesio-image]: https://img.shields.io/crates/v/nyx-space.svg
 [cratesio]: https://crates.io/crates/nyx-space

examples/04_lro_od/plot_od_rslt.py

@@ -9,7 +9,8 @@
 
 @click.command
 @click.option("-p", "--path", type=str, default="./04_lro_od_results.parquet")
-def main(path: str):
+@click.option("-f", "--full", type=bool, default=True)
+def main(path: str, full: bool):
     df = pl.read_parquet(path)
 
     df = (
@@ -123,6 +124,9 @@ def main(path: str):
         y=["Sigma Vx (RIC) (km/s)", "Sigma Vy (RIC) (km/s)", "Sigma Vz (RIC) (km/s)"],
     ).show()
 
+    if not full:
+        return
+
     # Load the RIC diff.
     for fname, errname in [
         ("04_lro_od_truth_error", "OD vs Flown"),

src/cosmic/orbitdual.rs

@@ -476,7 +476,7 @@ impl OrbitDual {
                 param: StateParameter::MeanAnomaly,
             })
         } else if self.ecc()?.real() > 1.0 {
-            info!("computing the hyperbolic anomaly");
+            debug!("computing the hyperbolic anomaly");
             // From GMAT's TrueToHyperbolicAnomaly
             Ok(OrbitPartial {
                 dual: ((self.ta_deg()?.dual.to_radians().sin()

src/dynamics/guidance/ruggiero.rs

@@ -17,6 +17,7 @@
 */
 
 use anise::prelude::Almanac;
+use serde::{Deserialize, Serialize};
 use snafu::ResultExt;
 
 use super::{
@@ -32,7 +33,7 @@ use std::fmt;
 use std::sync::Arc;
 
 /// Ruggiero defines the closed loop guidance law from IEPC 2011-102
-#[derive(Copy, Clone, Default)]
+#[derive(Copy, Clone, Default, Serialize, Deserialize)]
 pub struct Ruggiero {
     /// Stores the objectives
     pub objectives: [Option<Objective>; 5],

src/lib.rs

@@ -41,7 +41,7 @@ pub mod cosmic;
 /// Utility functions shared by different modules, and which may be useful to engineers.
 pub mod utils;
 
-mod errors;
+pub mod errors;
 /// Nyx will (almost) never panic and functions which may fail will return an error.
 pub use self::errors::NyxError;
 

src/mc/dispersion.rs

@@ -17,10 +17,11 @@
 */
 
 use crate::md::StateParameter;
+use serde::{Deserialize, Serialize};
 use typed_builder::TypedBuilder;
 
 /// A dispersions configuration, allows specifying min/max bounds (by default, they are not set)
-#[derive(Copy, Clone, TypedBuilder)]
+#[derive(Copy, Clone, TypedBuilder, Serialize, Deserialize)]
 pub struct StateDispersion {
     pub param: StateParameter,
     #[builder(default, setter(strip_option))]

src/md/events/details.rs

@@ -25,12 +25,13 @@ use crate::md::prelude::{Interpolatable, Traj};
 use crate::md::EventEvaluator;
 use crate::time::Duration;
 use core::fmt;
+use serde::{Deserialize, Serialize};
 use std::sync::Arc;
 
 /// Enumerates the possible edges of an event in a trajectory.
 ///
 /// `EventEdge` is used to describe the nature of a trajectory event, particularly in terms of its temporal dynamics relative to a specified condition or threshold. This enum helps in distinguishing whether the event is occurring at a rising edge, a falling edge, or if the edge is unclear due to insufficient data or ambiguous conditions.
-#[derive(Copy, Clone, Debug, PartialEq)]
+#[derive(Copy, Clone, Debug, PartialEq, Serialize, Deserialize)]
 pub enum EventEdge {
     /// Represents a rising edge of the event. This indicates that the event is transitioning from a lower to a higher evaluation of the event. For example, in the context of elevation, a rising edge would indicate an increase in elevation from a lower angle.
     Rising,
@@ -73,14 +74,14 @@ where
 {
     /// Generates detailed information about an event at a specific epoch in a trajectory.
     ///
-    /// This takes an `Epoch` as an input and returns a `Result<Self, NyxError>`.
+    /// This takes an `Epoch` as an input and returns a `Result<Self, EventError>`.
     /// It is designed to determine the state of a trajectory at a given epoch, evaluate the specific event at that state, and ascertain the nature of the event (rising, falling, or unclear).
     /// The initialization intelligently determines the edge type of the event by comparing the event's value at the current, previous, and next epochs.
     /// It ensures robust event characterization in trajectories.
     ///
     /// # Returns
     /// - `Ok(EventDetails<S>)` if the state at the given epoch can be determined and the event details are successfully evaluated.
-    /// - `Err(NyxError)` if there is an error in retrieving the state at the specified epoch.
+    /// - `Err(EventError)` if there is an error in retrieving the state at the specified epoch.
     ///
     pub fn new<E: EventEvaluator<S>>(
         state: S,

src/md/events/mod.rs

@@ -27,6 +27,7 @@ use crate::time::{Duration, Unit};
 use crate::State;
 use anise::prelude::{Almanac, Frame};
 use anise::structure::planetocentric::ellipsoid::Ellipsoid;
+use serde::{Deserialize, Serialize};
 
 use std::default::Default;
 use std::fmt;
@@ -59,15 +60,14 @@ where
 }
 
 /// Defines a state parameter event finder
-#[derive(Clone, Debug)]
-
+#[derive(Copy, Clone, Debug, Serialize, Deserialize)]
 pub struct Event {
     /// The state parameter
     pub parameter: StateParameter,
     /// The desired self.desired_value, must be in the same units as the state parameter
     pub desired_value: f64,
-    /// The time precision after which the solver will report that it cannot find any more precise
-    pub epoch_precision: Unit,
+    /// The duration precision after which the solver will report that it cannot find any more precise
+    pub epoch_precision: Duration,
     /// The precision on the desired value
     pub value_precision: f64,
     /// An optional frame in which to search this -- it IS recommended to convert the whole trajectory instead of searching in a given frame!
@@ -133,12 +133,12 @@ impl Event {
         parameter: StateParameter,
         desired_value: f64,
         value_precision: f64,
-        epoch_precision: Unit,
+        unit_precision: Unit,
     ) -> Self {
         Self {
             parameter,
             desired_value,
-            epoch_precision,
+            epoch_precision: 1 * unit_precision,
             value_precision,
             obs_frame: None,
         }
@@ -166,7 +166,7 @@ impl Event {
         Self {
             parameter,
             desired_value,
-            epoch_precision: Unit::Millisecond,
+            epoch_precision: Unit::Millisecond * 1,
             value_precision: 1e-3,
             obs_frame: Some(target_frame),
         }
@@ -179,7 +179,7 @@ impl Default for Event {
             parameter: StateParameter::Periapsis,
             desired_value: 0.0,
             value_precision: 1e-3,
-            epoch_precision: Unit::Second,
+            epoch_precision: Unit::Second * 1,
             obs_frame: None,
         }
     }

src/md/objective.rs

@@ -18,10 +18,11 @@
 
 use super::StateParameter;
 use crate::{errors::StateError, Spacecraft, State};
+use serde::{Deserialize, Serialize};
 use std::fmt;
 
 /// Defines a state parameter event finder
-#[derive(Copy, Clone, Debug)]
+#[derive(Copy, Clone, Debug, Serialize, Deserialize)]
 pub struct Objective {
     /// The state parameter to target
     pub parameter: StateParameter,

src/od/filter/kalman.rs

@@ -300,7 +300,14 @@ where
         // Compute the prefit ratio for the automatic rejection.
         // The measurement covariance is the square of the measurement itself.
         // So we compute its Cholesky decomposition to return to the non squared values.
-        let r_k_chol = s_k.clone().cholesky().ok_or(ODError::SingularNoiseRk)?.l();
+        let r_k_chol = match s_k.clone().cholesky() {
+            Some(r_k_clone) => r_k_clone.l(),
+            None => {
+                // In very rare case, when there isn't enough noise in the measurements,
+                // the inverting of S_k fails. If so, we revert back to the nominal Kalman derivation.
+                r_k.clone().cholesky().ok_or(ODError::SingularNoiseRk)?.l()
+            }
+        };
 
         // Compute the ratio as the average of each component of the prefit over the square root of the measurement
         // matrix r_k. Refer to ODTK MathSpec equation 4.10.

src/od/ground_station/mod.rs

@@ -24,7 +24,7 @@ use indexmap::{IndexMap, IndexSet};
 use snafu::ensure;
 
 use super::msr::MeasurementType;
-use super::noise::StochasticNoise;
+use super::noise::{GaussMarkov, StochasticNoise};
 use super::{ODAlmanacSnafu, ODError, ODTrajSnafu, TrackingDevice};
 use crate::io::ConfigRepr;
 use crate::od::NoiseNotConfiguredSnafu;
@@ -119,6 +119,34 @@ impl GroundStation {
         self
     }
 
+    /// Returns a copy of this ground station with the measurement type noises' constant bias set to the provided value.
+    pub fn with_msr_bias_constant(
+        mut self,
+        msr_type: MeasurementType,
+        bias_constant: f64,
+    ) -> Result<Self, ODError> {
+        if self.stochastic_noises.is_none() {
+            self.stochastic_noises = Some(IndexMap::new());
+        }
+
+        let stochastics = self.stochastic_noises.as_mut().unwrap();
+
+        let this_noise = stochastics
+            .get_mut(&msr_type)
+            .ok_or(ODError::NoiseNotConfigured {
+                kind: format!("{msr_type:?}"),
+            })
+            .unwrap();
+
+        if this_noise.bias.is_none() {
+            this_noise.bias = Some(GaussMarkov::ZERO);
+        }
+
+        this_noise.bias.unwrap().constant = Some(bias_constant);
+
+        Ok(self)
+    }
+
     /// Computes the azimuth and elevation of the provided object seen from this ground station, both in degrees.
     /// This is a shortcut to almanac.azimuth_elevation_range_sez.
     pub fn azimuth_elevation_of(

src/od/ground_station/trk_device.rs

@@ -174,4 +174,20 @@ impl TrackingDevice<Spacecraft> for GroundStation {
             })?
             .covariance(epoch))
     }
+
+    fn measurement_bias(&self, msr_type: MeasurementType, _epoch: Epoch) -> Result<f64, ODError> {
+        let stochastics = self.stochastic_noises.as_ref().unwrap();
+
+        if let Some(gm) = stochastics
+            .get(&msr_type)
+            .ok_or(ODError::NoiseNotConfigured {
+                kind: format!("{msr_type:?}"),
+            })?
+            .bias
+        {
+            Ok(gm.constant.unwrap_or(0.0))
+        } else {
+            Ok(0.0)
+        }
+    }
 }

src/od/msr/measurement.rs

@@ -23,7 +23,10 @@ use nalgebra::{allocator::Allocator, DefaultAllocator, DimName, OVector};
 use std::fmt;
 
 /// A type-agnostic simultaneous measurement storage structure. Allows storing any number of simultaneous measurement of a given taker.
-#[derive(Clone, Debug, PartialEq)]
+///
+/// Note that two measurements are considered equal if the tracker and epoch match exactly, and if both have the same measurement types,
+/// and those measurements are equal to within 1e-10 (this allows for some leeway in TDM producers).
+#[derive(Clone, Debug)]
 pub struct Measurement {
     /// Tracker alias which made this measurement
     pub tracker: String,
@@ -58,6 +61,7 @@ impl Measurement {
     where
         DefaultAllocator: Allocator<S>,
     {
+        // Consider adding a modulo modifier here, any bias should be configured by each ground station.
         let mut obs = OVector::zeros();
         for (i, t) in types.iter().enumerate() {
             if let Some(msr_value) = self.data.get(t) {
@@ -91,3 +95,17 @@ impl fmt::Display for Measurement {
         write!(f, "{} measured {} on {}", self.tracker, msrs, self.epoch)
     }
 }
+
+impl PartialEq for Measurement {
+    fn eq(&self, other: &Self) -> bool {
+        self.tracker == other.tracker
+            && self.epoch == other.epoch
+            && self.data.iter().all(|(key, &value)| {
+                if let Some(&other_value) = other.data.get(key) {
+                    (value - other_value).abs() < 1e-10
+                } else {
+                    false
+                }
+            })
+    }
+}