Merge pull request #169 from nyx-space/hotfix/python-traj-loading

Fix trajectory loading from Python

Cargo.toml

@@ -1,7 +1,7 @@
 [package]
 name = "nyx-space"
 build = "build.rs"
-version = "2.0.0-alpha.1"
+version = "2.0.0-alpha.2"
 edition = "2021"
 authors = ["Christopher Rabotin <[email protected]>"]
 description = "A high-fidelity space mission toolkit, with orbit propagation, estimation and some systems engineering"
@@ -46,8 +46,8 @@ numpy = {version = "0.17", optional = true}
 indicatif = {version = "0.17", features = ["rayon"]}
 rstats = "1.2.50"
 thiserror = "1.0"
-parquet = {version = "40.0.0", default-features = false, features = ["arrow", "brotli"]}
-arrow = "40.0.0"
+parquet = {version = "41.0.0", default-features = false, features = ["arrow", "brotli"]}
+arrow = "41.0.0"
 shadow-rs = {version = "0.23.0", default-features = false}
 serde_yaml = "0.9.21"
 whoami = "1.3.0"

src/cosmic/mod.rs

@@ -106,14 +106,20 @@ where
     }
 
     /// Return the value of the parameter, returns an error by default
-    fn value(&self, _param: StateParameter) -> Result<f64, NyxError> {
-        Err(NyxError::StateParameterUnavailable)
+    fn value(&self, param: StateParameter) -> Result<f64, NyxError> {
+        Err(NyxError::StateParameterUnavailable(
+            param,
+            "unimplemented in State trait".to_string(),
+        ))
     }
 
     /// Allows setting the value of the given parameter.
     /// NOTE: Most parameters where the `value` is available CANNOT be also set for that parameter (it's a much harder problem!)
-    fn set_value(&mut self, _param: StateParameter, _val: f64) -> Result<(), NyxError> {
-        Err(NyxError::StateParameterUnavailable)
+    fn set_value(&mut self, param: StateParameter, _val: f64) -> Result<(), NyxError> {
+        Err(NyxError::StateParameterUnavailable(
+            param,
+            "unimplemented in State trait".to_string(),
+        ))
     }
 }
 

src/cosmic/orbit.rs

@@ -861,6 +861,41 @@ impl Orbit {
 
 #[cfg_attr(feature = "python", pymethods)]
 impl Orbit {
+    #[cfg(feature = "python")]
+    #[getter]
+    fn get_epoch(&self) -> Epoch {
+        self.epoch
+    }
+    #[cfg(feature = "python")]
+    #[getter]
+    fn x_km(&self) -> f64 {
+        self.x_km
+    }
+    #[cfg(feature = "python")]
+    #[getter]
+    fn y_km(&self) -> f64 {
+        self.x_km
+    }
+    #[cfg(feature = "python")]
+    #[getter]
+    fn z_km(&self) -> f64 {
+        self.x_km
+    }
+    #[cfg(feature = "python")]
+    #[getter]
+    fn vx_km_s(&self) -> f64 {
+        self.vx_km_s
+    }
+    #[cfg(feature = "python")]
+    #[getter]
+    fn vy_km_s(&self) -> f64 {
+        self.vy_km_s
+    }
+    #[cfg(feature = "python")]
+    #[getter]
+    fn vz_km_s(&self) -> f64 {
+        self.vz_km_s
+    }
     /// Returns the magnitude of the radius vector in km
     pub fn rmag_km(&self) -> f64 {
         (self.x_km.powi(2) + self.y_km.powi(2) + self.z_km.powi(2)).sqrt()
@@ -1500,6 +1535,11 @@ impl Orbit {
         format!("{self:x}")
     }
 
+    #[cfg(feature = "python")]
+    fn __eq__(&self, other: &Self) -> bool {
+        self == other
+    }
+
     /// Creates a new Orbit in the provided frame at the provided Epoch given the position and velocity components.
     ///
     /// **Units:** km, km, km, km/s, km/s, km/s
@@ -2061,7 +2101,10 @@ impl State for Orbit {
             StateParameter::VX => Ok(self.vx_km_s),
             StateParameter::VY => Ok(self.vy_km_s),
             StateParameter::VZ => Ok(self.vz_km_s),
-            _ => Err(NyxError::StateParameterUnavailable),
+            _ => Err(NyxError::StateParameterUnavailable(
+                param,
+                "no such parameter for orbit structure".to_string(),
+            )),
         }
     }
 
@@ -2097,7 +2140,12 @@ impl State for Orbit {
                 self.vy_km_s = new_radius.y;
                 self.vz_km_s = new_radius.z;
             }
-            _ => return Err(NyxError::StateParameterUnavailable),
+            _ => {
+                return Err(NyxError::StateParameterUnavailable(
+                    param,
+                    "not settable for orbit structure with set_value".to_string(),
+                ))
+            }
         }
         Ok(())
     }

src/dynamics/guidance/ruggiero.rs

@@ -133,7 +133,10 @@ impl Ruggiero {
                 Ok(num / denom)
             }
             StateParameter::AoP => Ok(1.0),
-            _ => Err(NyxError::StateParameterUnavailable),
+            _ => Err(NyxError::StateParameterUnavailable(
+                *parameter,
+                "not a control variable in Ruggiero".to_string(),
+            )),
         }
     }
 

src/errors.rs

@@ -19,6 +19,7 @@
 use super::thiserror::Error;
 use crate::io::ConfigError;
 use crate::md::trajectory::TrajError;
+use crate::md::StateParameter;
 pub use crate::md::TargetingError;
 pub use crate::time::Errors as TimeErrors;
 use crate::Spacecraft;
@@ -75,8 +76,8 @@ pub enum NyxError {
     #[error("Partials for this model are not defined")]
     PartialsUndefined,
     /// State parameter cannot be used in this function
-    #[error("State parameter cannot be used in this function")]
-    StateParameterUnavailable,
+    #[error("Unavailable parameter {0:?}: {1}")]
+    StateParameterUnavailable(StateParameter, String),
     /// Could not load file
     #[error("Could not load file: {0}")]
     LoadingError(String),

src/io/trajectory_data.rs

@@ -100,7 +100,7 @@ impl TrajectoryLoader {
         let mut has_epoch = false; // Required
         let mut frame = None;
 
-        let mut found_fields = [
+        let mut found_fields = vec![
             (StateParameter::X, false),
             (StateParameter::Y, false),
             (StateParameter::Z, false),
@@ -166,12 +166,21 @@ impl TrajectoryLoader {
 
         let expected_type = std::any::type_name::<S>().split("::").last().unwrap();
 
-        if expected_type == "Spacecraft" && !sc_compat {
-            // Oops, missing fuel data (using the full call in case the field name changes in the future)
-            return Err(Box::new(NyxError::FileUnreadable(format!(
-                "Missing `{}` field",
-                found_fields.last().unwrap().0.to_field(None).name()
-            ))));
+        if expected_type == "Spacecraft" {
+            if !sc_compat {
+                // Oops, missing fuel data (using the full call in case the field name changes in the future)
+                return Err(Box::new(NyxError::FileUnreadable(format!(
+                    "Missing `{}` field",
+                    found_fields.last().unwrap().0.to_field(None).name()
+                ))));
+            }
+        } else if sc_compat {
+            // Not a spacecraft, remove the fuel mass
+            if let Some(last_field) = found_fields.last_mut() {
+                if last_field.0 == StateParameter::FuelMass && last_field.1 {
+                    last_field.1 = false;
+                }
+            }
         }
 
         // At this stage, we know that the measurement is valid and the conversion is supported.
@@ -201,8 +210,8 @@ impl TrajectoryLoader {
                 );
             }
 
-            if sc_compat {
-                // Read the fuel
+            if expected_type == "Spacecraft" {
+                // Read the fuel only if this is a spacecraft we're building
                 shared_data.push(
                     batch
                         .column_by_name("fuel_mass (kg)")

src/mc/results.rs

@@ -225,7 +225,10 @@ where
                 }
             }
             // Oh, this parameter was not found!
-            return Err(NyxError::StateParameterUnavailable);
+            return Err(NyxError::StateParameterUnavailable(
+                param,
+                "not among dispersions of Monte Carlo setup".to_string(),
+            ));
         }
         Ok(report)
     }

src/md/trajectory/orbit_traj.rs

@@ -192,6 +192,18 @@ impl Traj<Orbit> {
                             stm: None,
                         };
 
+                        // Check that we don't have duplicate epochs
+                        if let Some(prev) = traj.states.last() {
+                            if prev.epoch == orbit.epoch {
+                                warn!(
+                                    "[line: {}] duplicate epoch {} in OEM -- overwriting",
+                                    lno + 1,
+                                    orbit.epoch
+                                );
+                                traj.states.pop();
+                            }
+                        }
+
                         traj.states.push(orbit);
                     }
                     Err(_) => {

src/polyfit/hermite.rs

@@ -235,14 +235,14 @@ pub fn hermite_eval(
     x_eval: f64,
 ) -> Result<(f64, f64), NyxError> {
     if xs.len() != ys.len() || xs.len() != ydots.len() {
-        error!("Abscissas (xs), ordinates (ys), and first derivatives (ydots) must contain the same number of items, but they are of lengths {}, {}, and {}", xs.len(), ys.len(), ydots.len());
-        return Err(NyxError::MathDomain("=(".to_string()));
+        let msg = format!("Abscissas (xs), ordinates (ys), and first derivatives (ydots) must contain the same number of items, but they are of lengths {}, {}, and {}", xs.len(), ys.len(), ydots.len());
+        return Err(NyxError::MathDomain(msg));
     } else if xs.is_empty() {
-        error!("No interpolation data provided");
-        return Err(NyxError::MathDomain("=(".to_string()));
+        let msg = "No interpolation data provided".to_string();
+        return Err(NyxError::MathDomain(msg));
     } else if xs.len() > 3 * INTERPOLATION_SAMPLES {
-        error!("More than {} samples provided, which is the maximum number of items allowed for a Hermite interpolation", 3 * INTERPOLATION_SAMPLES);
-        return Err(NyxError::MathDomain("=(".to_string()));
+        let msg = format!("More than {} samples provided, which is the maximum number of items allowed for a Hermite interpolation", 3 * INTERPOLATION_SAMPLES);
+        return Err(NyxError::MathDomain(msg));
     }
 
     // At this point, we know that the lengths of items is correct, so we can directly address them without worry for overflowing the array.
@@ -274,7 +274,8 @@ pub fn hermite_eval(
         let c2 = x_eval - xs[i - 1];
         let denom = xs[i] - xs[i - 1];
         if denom.abs() < f64::EPSILON {
-            return Err(NyxError::MathDomain("=(".to_string()));
+            let msg = format!("duplicate abscissa data: denominator near zero ({denom:e})");
+            return Err(NyxError::MathDomain(msg));
         }
 
         /*        The second column of WORK contains interpolated derivative */
@@ -329,7 +330,10 @@ pub fn hermite_eval(
             let c2 = x_eval - xs[xi - 1];
             let denom = xs[xij - 1] - xs[xi - 1];
             if denom.abs() < f64::EPSILON {
-                return Err(NyxError::MathDomain("=(".to_string()));
+                let msg = format!(
+                    "duplicate abscissa data in interp table: denominator near zero ({denom:e})"
+                );
+                return Err(NyxError::MathDomain(msg));
             }
 
             /*           Compute the interpolated derivative at X for the Ith */

src/python/mission_design/spacecraft.rs

@@ -77,8 +77,8 @@ impl Spacecraft {
         format!("{self}\n{self:x}")
     }
 
-    fn __eq__(&self, other: Self) -> bool {
-        *self == other
+    fn __eq__(&self, other: &Self) -> bool {
+        self == other
     }
 
     /// Note: this returns a COPY of the orbit, not a mutable reference to it!

tests/python/test_mission_design.py

@@ -10,7 +10,7 @@
     SpacecraftDynamics,
     TrajectoryLoader,
 )
-from nyx_space.time import Duration, Unit, Epoch
+from nyx_space.time import Duration, Unit, Epoch, TimeSeries
 
 
 def test_propagate():
@@ -147,6 +147,25 @@ def test_build_spacecraft():
     traj_pkl = pickle.dumps(traj)
     traj_unpkl = pickle.loads(traj_pkl)
     assert traj_unpkl.__eq__(traj)
+    # Check that we can convert this to a spacecraft trajectory
+    traj_sc = traj.to_spacecraft_traj()
+    traj_orbit = traj.to_orbit_traj()
+    traj_orbit_dc = traj_sc.downcast()
+    # Check that we can query it (will raise an exception if we can't, thereby failing the test)
+    ts = TimeSeries(Epoch("2020-06-01T12:00:00.000000"), Epoch("2020-06-01T13:00:00.000000"), step=Unit.Minute*17 + Unit.Second*13.8, inclusive=True)
+    for epoch in ts:
+        orbit = traj_orbit.at(epoch)
+        dc_orbit = traj_orbit_dc.at(epoch)
+        sc_orbit = traj_sc.at(epoch).orbit
+        # Check params individually
+        assert orbit.epoch == sc_orbit.epoch
+        assert orbit.x_km == sc_orbit.x_km
+        assert orbit.vx_km_s == sc_orbit.vx_km_s
+        # Check the downcasted version
+        assert dc_orbit.x_km == sc_orbit.x_km
+        assert dc_orbit.vx_km_s == sc_orbit.vx_km_s
+        assert dc_orbit.__eq__(sc_orbit)
+
 
 
 if __name__ == "__main__":