Skip to content

Commit

Permalink
generalize covariance representation
Browse files Browse the repository at this point in the history
  • Loading branch information
dahlend committed Jun 11, 2024
1 parent 944efe8 commit a050cff
Show file tree
Hide file tree
Showing 10 changed files with 140 additions and 71 deletions.
4 changes: 4 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0

### Changed

- Renamed "HorizonsCovariance" to "Covariance" and generalized it to support both
cometary and cartesian representations.
- Made api more consistent for conversion to State objects by renaming all instances of
`.as_state` to `.state`.
- Improved performance of integration and propagation by about 10%.
- Restructured the Rust FOVs to be organized by observatory.
- Renamed `SpiceKernels.cache_kernel_reload` to `kernel_reload` and changed it's args.
Expand Down
2 changes: 1 addition & 1 deletion src/examples/on_sky.py
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@
)

# These can be converted to cartesian coordinates
state = elements.as_state
state = elements.state

# States for objects may also be queried from JPL Horizons
# state = neospy.HorizonsProperties.fetch("Ceres").state
Expand Down
2 changes: 1 addition & 1 deletion src/examples/plot_uncertainty.py
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@
sample = sample + values
params = dict(zip(labels, sample))
states.append(
neospy.CometElements(obj.desig, obj.covariance.epoch, **params).as_state
neospy.CometElements(obj.desig, obj.covariance.epoch, **params).state
)

# Propagate the position of all states to all time steps, recording the V mags
Expand Down
4 changes: 2 additions & 2 deletions src/neospy/horizons.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
import numpy as np

# pylint: disable=import-error
from ._core import HorizonsProperties, HorizonsCovariance # type: ignore
from ._core import HorizonsProperties, Covariance # type: ignore

__all__ = ["HorizonsProperties"]

Expand Down Expand Up @@ -107,7 +107,7 @@ def fetch(name, update_name=True, cache=True, update_cache=False, exact_name=Fal
lab: phys[lookup_rev[lab]] for lab in labels if lab in lookup_rev
}
params = [(lookup_rev.get(x, x), elements.get(x, np.nan)) for x in labels]
phys["covariance"] = HorizonsCovariance(cov_epoch, params, mat)
phys["covariance"] = Covariance(name, cov_epoch, params, mat)
else:
raise ValueError(
"Horizons did not return orbit information for this object:" f"\n{props}"
Expand Down
2 changes: 1 addition & 1 deletion src/neospy/mpc.py
Original file line number Diff line number Diff line change
Expand Up @@ -771,7 +771,7 @@ def table_to_states(orbit_dataframe):
item.peri_arg,
item.peri_time,
item.lon_node,
).as_state
).state
)
return states

Expand Down
111 changes: 111 additions & 0 deletions src/neospy/rust/covariance.rs
Original file line number Diff line number Diff line change
@@ -0,0 +1,111 @@
use std::{collections::HashMap, fmt::Debug};

use crate::{elements::PyCometElements, vector::VectorLike};
use crate::state::PyState;
use neospy_core::{errors::NEOSpyError, io::FileIO};
use pyo3::prelude::*;
use serde::{Deserialize, Serialize};



/// Covariance uncertainty representation of an objects state.
#[pyclass(frozen, get_all, module = "neospy")]
#[derive(Clone, Debug, Deserialize, Serialize)]
pub struct Covariance {
/// Designation of the object
desig: String,

/// Epoch of the covariance matrix fit.
epoch: f64,

/// Name and best estimate of the parameters in the fit.
params: Vec<(String, f64)>,

/// The covariance matrix, where the order of the array corresponds to the parameters.
cov_matrix: Vec<Vec<f64>>,
}

impl neospy_core::io::FileIO for Covariance {}

#[pymethods]
impl Covariance {
#[new]
#[allow(clippy::too_many_arguments)]
pub fn new(desig:String, epoch: f64, params: Vec<(String, f64)>, cov_matrix: Vec<Vec<f64>>) -> Self {
Self {
desig,
epoch,
params,
cov_matrix,
}
}

fn __repr__(&self) -> String {
format!(
"Covariance(desig={:?}, epoch={:?}, params={:?}, cov_matrix={:?})",
self.desig, self.epoch, self.params, self.cov_matrix,
)
}

/// Create a State object from the fit values.
///
/// This looks for either cometary elements or cartesian elements in the parameters.
/// Cometary Elements must contain the following complete set of keys:
/// ["eccentricity", "peri_dist", "peri_time", "lon_of_ascending", "peri_arg", "inclination"]
/// Cartesian must contain the following complete set of keys:
/// ['x', 'y', 'z', 'vx', 'vy', 'vz']
///
/// All units must be in degrees, AU, or AU/Day as appropriate.
#[getter]
pub fn state(&self) -> PyResult<PyState>{
let epoch = self.epoch;
let desig = self.desig.clone();
let mut hash = HashMap::new();
for (key, val) in self.params.iter(){
let lower_key = key.to_lowercase();
if hash.insert(lower_key, *val).is_some(){
return Err(NEOSpyError::IOError(format!("Repeat parameter {:?} present in covariance", &key)))?;
}
}

if hash.contains_key("eccentricity"){
let eccentricity = *hash.get("eccentricity").ok_or(NEOSpyError::ValueError("Covariance missing 'eccentricity'".into()))?;
let peri_dist = *hash.get("peri_dist").ok_or(NEOSpyError::ValueError("Covariance missing 'peri_dist'".into()))?;
let peri_time = *hash.get("peri_time").ok_or(NEOSpyError::ValueError("Covariance missing 'peri_time'".into()))?;
let lon_of_ascending = *hash.get("lon_of_ascending").ok_or(NEOSpyError::ValueError("Covariance missing 'lon_of_ascending'".into()))?;
let peri_arg = *hash.get("peri_arg").ok_or(NEOSpyError::ValueError("Covariance missing 'peri_arg'".into()))?;
let inclination = *hash.get("inclination").ok_or(NEOSpyError::ValueError("Covariance missing 'inclination'".into()))?;
let elem = PyCometElements::new(desig, epoch, eccentricity, inclination, peri_dist, peri_arg, peri_time, lon_of_ascending);
elem.state()

} else if hash.contains_key("x"){
let x = *hash.get("x").ok_or(NEOSpyError::ValueError("Covariance missing 'x'".into()))?;
let y = *hash.get("y").ok_or(NEOSpyError::ValueError("Covariance missing 'y'".into()))?;
let z = *hash.get("z").ok_or(NEOSpyError::ValueError("Covariance missing 'z'".into()))?;
let vx = *hash.get("vx").ok_or(NEOSpyError::ValueError("Covariance missing 'vx'".into()))?;
let vy = *hash.get("vy").ok_or(NEOSpyError::ValueError("Covariance missing 'vy'".into()))?;
let vz = *hash.get("vz").ok_or(NEOSpyError::ValueError("Covariance missing 'vz'".into()))?;
let pos = VectorLike::Arr([x, y, z]);
let vel = VectorLike::Arr([vx, vy, vz]);
Ok(PyState::new(Some(desig), epoch, pos, vel, None, None))
} else{
Err(NEOSpyError::ValueError("Covariance cannot be converted to a state, \
the covariance parameters must either contain: \
['x', 'y', 'z', 'vx', 'vy', 'vz'] or \
['eccentricity', 'peri_dist', 'peri_time', 'lon_of_ascending', 'peri_arg', 'inclination']".into()))?
}
}

/// Save the covariance matrix to a file.
#[pyo3(name = "save")]
pub fn py_save(&self, filename: String) -> PyResult<usize> {
Ok(self.save(filename)?)
}

/// Load a covariance matrix from a file.
#[staticmethod]
#[pyo3(name = "load")]
pub fn py_load(filename: String) -> PyResult<Self> {
Ok(Self::load(filename)?)
}
}
2 changes: 1 addition & 1 deletion src/neospy/rust/elements.rs
Original file line number Diff line number Diff line change
Expand Up @@ -166,7 +166,7 @@ impl PyCometElements {

/// Convert the orbital elements into a cartesian State.
#[getter]
pub fn as_state(&self) -> PyResult<PyState> {
pub fn state(&self) -> PyResult<PyState> {
Ok(self.0.try_to_state()?.into())
}

Expand Down
76 changes: 14 additions & 62 deletions src/neospy/rust/horizons.rs
Original file line number Diff line number Diff line change
Expand Up @@ -2,59 +2,11 @@ use std::fmt::Debug;

use crate::elements::PyCometElements;
use crate::state::PyState;
use neospy_core::{io::FileIO, prelude};
use crate::covariance::Covariance;
use neospy_core::{ io::FileIO, prelude};
use pyo3::prelude::*;
use serde::{Deserialize, Serialize};

/// Horizons object properties
/// Physical, orbital, and observational properties of a solar system object as recorded in JPL Horizons.
#[pyclass(frozen, get_all, module = "neospy")]
#[derive(Clone, Debug, Deserialize, Serialize)]
pub struct HorizonsCovariance {
/// Epoch of the covariance matrix fit.
epoch: f64,

/// Name and best estimate of the parameters in the fit.
params: Vec<(String, f64)>,

/// The covariance matrix, where the order of the array corresponds to the parameters.
cov_matrix: Vec<Vec<f64>>,
}

impl neospy_core::io::FileIO for HorizonsCovariance {}

#[pymethods]
impl HorizonsCovariance {
#[new]
#[allow(clippy::too_many_arguments)]
pub fn new(epoch: f64, params: Vec<(String, f64)>, cov_matrix: Vec<Vec<f64>>) -> Self {
Self {
epoch,
params,
cov_matrix,
}
}

fn __repr__(&self) -> String {
format!(
"HorizonsCovariance(epoch={:?}, params={:?}, cov_matrix={:?})",
self.epoch, self.params, self.cov_matrix,
)
}

/// Save the covariance matrix to a file.
#[pyo3(name = "save")]
pub fn py_save(&self, filename: String) -> PyResult<usize> {
Ok(self.save(filename)?)
}

/// Load a covariance matrix from a file.
#[staticmethod]
#[pyo3(name = "load")]
pub fn py_load(filename: String) -> PyResult<Self> {
Ok(Self::load(filename)?)
}
}

/// Horizons object properties
/// Physical, orbital, and observational properties of a solar system object as recorded in JPL Horizons.
Expand Down Expand Up @@ -108,7 +60,7 @@ pub struct HorizonsProperties {
arc_len: Option<f64>,

/// Covariance values in the orbit fit.
covariance: Option<HorizonsCovariance>,
covariance: Option<Covariance>,
}

impl neospy_core::io::FileIO for HorizonsProperties {}
Expand All @@ -133,7 +85,7 @@ impl HorizonsProperties {
moid: Option<f64>,
g_phase: Option<f64>,
arc_len: Option<f64>,
covariance: Option<HorizonsCovariance>,
covariance: Option<Covariance>,
) -> Self {
Self {
desig,
Expand Down Expand Up @@ -162,31 +114,31 @@ impl HorizonsProperties {
desig: prelude::Desig::Name(self.desig.clone()),
epoch: self
.epoch
.ok_or_else(|| prelude::NEOSpyError::ValueError("No Epoch defined".into()))?,
eccentricity: self.eccentricity.ok_or_else(|| {
.ok_or(prelude::NEOSpyError::ValueError("No Epoch defined".into()))?,
eccentricity: self.eccentricity.ok_or(
prelude::NEOSpyError::ValueError("No Eccentricity defined".into())
})?,
)?,
inclination: self
.inclination
.ok_or_else(|| prelude::NEOSpyError::ValueError("No Inclination defined".into()))?
.ok_or(prelude::NEOSpyError::ValueError("No Inclination defined".into()))?
.to_radians(),
peri_arg: self
.peri_arg
.ok_or_else(|| prelude::NEOSpyError::ValueError("No peri_arg defined".into()))?
.ok_or(prelude::NEOSpyError::ValueError("No peri_arg defined".into()))?
.to_radians(),
peri_dist: self
.peri_dist
.ok_or_else(|| prelude::NEOSpyError::ValueError("No peri_dist defined".into()))?,
.ok_or(prelude::NEOSpyError::ValueError("No peri_dist defined".into()))?,
peri_time: self
.peri_time
.ok_or_else(|| prelude::NEOSpyError::ValueError("No peri_time defined".into()))?,
.ok_or(prelude::NEOSpyError::ValueError("No peri_time defined".into()))?,
lon_of_ascending: self
.lon_of_ascending
.ok_or_else(|| {
.ok_or(
prelude::NEOSpyError::ValueError(
"No longitude of ascending node defined".into(),
)
})?
)?
.to_radians(),
frame: prelude::Frame::Ecliptic,
}))
Expand All @@ -195,7 +147,7 @@ impl HorizonsProperties {
/// Convert the orbital elements of the object to a State.
#[getter]
pub fn state(&self) -> PyResult<PyState> {
self.elements()?.as_state()
self.elements()?.state()
}

fn __repr__(&self) -> String {
Expand Down
4 changes: 3 additions & 1 deletion src/neospy/rust/lib.rs
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
use pyo3::prelude::{pymodule, wrap_pyfunction, Bound, PyModule, PyResult, Python};

mod covariance;
mod elements;
mod fitting;
mod flux;
Expand Down Expand Up @@ -46,7 +47,8 @@ fn _core(_py: Python, m: &Bound<'_, PyModule>) -> PyResult<()> {
m.add_class::<flux::PyModelResults>()?;

m.add_class::<horizons::HorizonsProperties>()?;
m.add_class::<horizons::HorizonsCovariance>()?;

m.add_class::<covariance::Covariance>()?;

m.add_function(wrap_pyfunction!(frame::frame_change_py, m)?)?;
m.add_function(wrap_pyfunction!(frame::wgs_lat_lon_to_ecef, m)?)?;
Expand Down
4 changes: 2 additions & 2 deletions src/tests/test_state.py
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ def test_from_orbital_elements_specific(self):
peri_arg=0,
lon_of_ascending=0,
peri_time=123456,
).as_state
).state
assert vs.jd == 123456
elem = vs.elements
assert np.isclose(elem.eccentricity, 0.1)
Expand All @@ -53,7 +53,7 @@ def test_from_orbital_elements(self, ecc, inc, peri_dist, peri_arg, lon, time):
peri_arg=peri_arg,
lon_of_ascending=lon,
peri_time=123456 + time,
).as_state
).state
assert state.jd == 123456
elements = state.elements
assert np.isclose(elements.eccentricity, ecc)
Expand Down

0 comments on commit a050cff

Please sign in to comment.