nyx_space/od/noise/
mod.rs

1/*
2    Nyx, blazing fast astrodynamics
3    Copyright (C) 2018-onwards Christopher Rabotin <christopher.rabotin@gmail.com>
4
5    This program is free software: you can redistribute it and/or modify
6    it under the terms of the GNU Affero General Public License as published
7    by the Free Software Foundation, either version 3 of the License, or
8    (at your option) any later version.
9
10    This program is distributed in the hope that it will be useful,
11    but WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13    GNU Affero General Public License for more details.
14
15    You should have received a copy of the GNU Affero General Public License
16    along with this program.  If not, see <https://www.gnu.org/licenses/>.
17*/
18
19use crate::io::watermark::pq_writer;
20use arrow::array::{ArrayRef, Float64Array, UInt32Array};
21use arrow::datatypes::{DataType, Field, Schema};
22use arrow::record_batch::RecordBatch;
23use der::{Decode, Encode, Reader};
24use hifitime::{Epoch, TimeSeries, TimeUnits};
25use parquet::arrow::ArrowWriter;
26
27use rand::rngs::SysRng;
28use rand::{Rng, SeedableRng};
29use rand_pcg::Pcg64Mcg;
30use serde::{Deserialize, Serialize};
31use std::error::Error;
32use std::fmt::Display;
33use std::fs::File;
34use std::ops::{Mul, MulAssign};
35use std::path::Path;
36use std::sync::Arc;
37
38pub mod gauss_markov;
39pub mod link_specific;
40pub mod white;
41
42#[cfg(feature = "python")]
43use hifitime::Duration;
44#[cfg(feature = "python")]
45use pyo3::exceptions::PyValueError;
46#[cfg(feature = "python")]
47use pyo3::prelude::*;
48#[cfg(feature = "python")]
49use pyo3::types::PyType;
50
51pub use gauss_markov::GaussMarkov;
52pub use white::WhiteNoise;
53
54/// Trait for any kind of stochastic modeling, developing primarily for synthetic orbit determination measurements.
55pub trait Stochastics {
56    /// Return the variance of this stochastic noise model at a given time.
57    fn covariance(&self, epoch: Epoch) -> f64;
58
59    /// Returns a new sample of these stochastics
60    fn sample<R: Rng>(&mut self, epoch: Epoch, rng: &mut R) -> f64;
61}
62
63/// Stochastic noise modeling used primarily for synthetic orbit determination measurements.
64///
65/// This implementation distinguishes between the white noise model and the bias model. It also includes a constant offset.
66#[derive(Copy, Clone, Debug, Default, PartialEq, Serialize, Deserialize)]
67#[cfg_attr(feature = "python", pyclass(from_py_object, get_all, set_all))]
68pub struct StochasticNoise {
69    pub white_noise: Option<WhiteNoise>,
70    pub bias: Option<GaussMarkov>,
71}
72
73impl StochasticNoise {
74    /// Zero noise stochastic process.
75    pub const ZERO: Self = Self {
76        white_noise: None,
77        bias: None,
78    };
79
80    /// The minimum stochastic noise process with a zero mean white noise of 1e-6.
81    pub const MIN: Self = Self {
82        white_noise: Some(WhiteNoise {
83            mean: 0.0,
84            sigma: 1e-6,
85        }),
86        bias: None,
87    };
88
89    /// Default stochastic process of the Deep Space Network, as per DESCANSO Chapter 3, Table 3-3.
90    /// Using the instrument bias as the white noise value, zero constant bias.
91    pub fn default_range_km() -> Self {
92        Self {
93            white_noise: Some(WhiteNoise {
94                sigma: 2.0e-3, // 2 m
95                ..Default::default()
96            }),
97            // Until Nyx can support bias estimation the default bias is None, cf. https://github.com/nyx-space/nyx/issues/326
98            // bias: Some(GaussMarkov::default_range_km()),
99            bias: None,
100        }
101    }
102
103    /// Default stochastic process of the Deep Space Network, using as per DESCANSO Chapter 3, Table 3-3 for the GM process.
104    pub fn default_doppler_km_s() -> Self {
105        Self {
106            white_noise: Some(WhiteNoise {
107                sigma: 3e-6, // 3 mm/s
108                ..Default::default()
109            }),
110            // Until Nyx can support bias estimation the default bias is None, cf. https://github.com/nyx-space/nyx/issues/326
111            // bias: Some(GaussMarkov::default_doppler_km_s()),
112            bias: None,
113        }
114    }
115
116    /// Default stochastic process for an angle measurement (azimuth or elevation)
117    /// Using the instrument bias as the white noise value, zero constant bias.
118    pub fn default_angle_deg() -> Self {
119        Self {
120            white_noise: Some(WhiteNoise {
121                sigma: 1.0e-2, // 0.01 deg
122                ..Default::default()
123            }),
124            // Until Nyx can support bias estimation the default bias is None, cf. https://github.com/nyx-space/nyx/issues/326
125            // bias: Some(GaussMarkov::default_range_km()),
126            bias: None,
127        }
128    }
129
130    /// Sample these stochastics
131    pub fn sample<R: Rng>(&mut self, epoch: Epoch, rng: &mut R) -> f64 {
132        let mut sample = 0.0;
133        if let Some(wn) = &mut self.white_noise {
134            sample += wn.sample(epoch, rng)
135        }
136        if let Some(gm) = &mut self.bias {
137            sample += gm.sample(epoch, rng);
138        }
139        sample
140    }
141
142    /// Executes a hardcoded 24-hour Monte Carlo simulation of the stochastic model, exporting the time history to a Parquet file.
143    ///
144    /// # Warning: Hardcoded Time Series & Diagnostic Data Gaps
145    /// This method does *not* accept a user-defined tracking schedule or time series. It inherently evaluates the stochastic process
146    /// over a strict 24-hour period, beginning at the exact system clock moment of method execution, utilizing a 1-minute step size.
147    ///
148    /// Furthermore, users will observe exactly 1,082 samples per simulation run, rather than the 1,441 samples expected from a
149    /// continuous 24-hour 1-minute cadence. The simulation intentionally drops all epochs strictly greater than +6 hours and
150    /// strictly less than +12 hours from the start time. This hardcoded artifact is designed to demonstrate variance bounds
151    /// expansion in the absence of measurements (e.g., simulating a tracking dropout for a Gauss-Markov bias).
152    ///
153    /// # Algorithm
154    /// 1. Establish `start` as the system clock time at invocation.
155    /// 2. Construct an inclusive time series from `start` to `start + 24 hours` at 1-minute intervals.
156    /// 3. For each configured run, seed a PRNG (`Pcg64Mcg`) using system entropy.
157    /// 4. Evaluate the process covariance and sample the stochastic noise at each epoch.
158    /// 5. Discard all epochs inside the `(start + 6h, start + 12h)` open interval.
159    /// 6. Export the remaining 1,082 samples per run to an Apache Arrow RecordBatch and write to disk via Parquet.
160    pub fn simulate<P: AsRef<Path>>(
161        self,
162        path: P,
163        runs: Option<u32>,
164        unit: Option<String>,
165    ) -> Result<Vec<StochasticState>, Box<dyn Error>> {
166        let num_runs = runs.unwrap_or(25);
167
168        let start = Epoch::now().unwrap();
169        let (step, end) = (1.minutes(), start + 1.days());
170
171        let capacity = ((end - start).to_seconds() / step.to_seconds()).ceil() as usize;
172
173        let mut samples = Vec::with_capacity(capacity);
174
175        for run in 0..num_runs {
176            let mut rng = Pcg64Mcg::try_from_rng(&mut SysRng).unwrap();
177
178            let mut mdl = self;
179            for epoch in TimeSeries::inclusive(start, end, step) {
180                if epoch > start + 6.hours() && epoch < start + 12.hours() {
181                    // Skip to see how the variance changes.
182                    continue;
183                }
184                let variance = mdl.covariance(epoch);
185                let sample = mdl.sample(epoch, &mut rng);
186                samples.push(StochasticState {
187                    run,
188                    dt_s: (epoch - start).to_seconds(),
189                    sample,
190                    variance,
191                });
192            }
193        }
194
195        let bias_unit = match unit {
196            Some(unit) => format!("({unit})"),
197            None => "(unitless)".to_string(),
198        };
199
200        // Build the parquet file
201        let hdrs = vec![
202            Field::new("Run", DataType::UInt32, false),
203            Field::new("Delta Time (s)", DataType::Float64, false),
204            Field::new(format!("Bias {bias_unit}"), DataType::Float64, false),
205            Field::new(format!("Variance {bias_unit}"), DataType::Float64, false),
206        ];
207
208        let schema = Arc::new(Schema::new(hdrs));
209        let record = vec![
210            Arc::new(UInt32Array::from(
211                samples.iter().map(|s| s.run).collect::<Vec<u32>>(),
212            )) as ArrayRef,
213            Arc::new(Float64Array::from(
214                samples.iter().map(|s| s.dt_s).collect::<Vec<f64>>(),
215            )) as ArrayRef,
216            Arc::new(Float64Array::from(
217                samples.iter().map(|s| s.sample).collect::<Vec<f64>>(),
218            )) as ArrayRef,
219            Arc::new(Float64Array::from(
220                samples.iter().map(|s| s.variance).collect::<Vec<f64>>(),
221            )) as ArrayRef,
222        ];
223
224        let props = pq_writer(None);
225
226        let file = File::create(path)?;
227        let mut writer = ArrowWriter::try_new(file, schema.clone(), props).unwrap();
228
229        let batch = RecordBatch::try_new(schema, record)?;
230        writer.write(&batch)?;
231        writer.close()?;
232
233        Ok(samples)
234    }
235
236    fn available_data(&self) -> u8 {
237        let mut bits: u8 = 0;
238
239        if self.white_noise.is_some() {
240            bits |= 1 << 0;
241        }
242
243        if self.bias.is_some() {
244            bits |= 1 << 1;
245        }
246
247        bits
248    }
249}
250
251#[cfg_attr(feature = "python", pymethods)]
252impl StochasticNoise {
253    #[cfg(feature = "python")]
254    #[pyo3(signature=(white_noise=None, bias=None, name=None))]
255    #[new]
256    fn py_new(
257        white_noise: Option<WhiteNoise>,
258        bias: Option<GaussMarkov>,
259        name: Option<String>,
260    ) -> PyResult<Self> {
261        if let Some(name) = name {
262            match name.to_ascii_lowercase().as_str() {
263                "range" => Ok(Self::default_range_km()),
264                "doppler" => Ok(Self::default_doppler_km_s()),
265                "angles" => Ok(Self::default_angle_deg()),
266                _ => Err(PyValueError::new_err(format!(
267                    "name must be `range`, `doppler`, or `angles` (received `{name}`)"
268                ))),
269            }
270        } else {
271            Ok(Self { white_noise, bias })
272        }
273    }
274
275    /// Return the covariance of these stochastics at a given time.
276    pub fn covariance(&self, epoch: Epoch) -> f64 {
277        let mut variance = 0.0;
278        if let Some(wn) = &self.white_noise {
279            variance += wn.covariance(epoch);
280        }
281        if let Some(gm) = &self.bias {
282            variance += gm.covariance(epoch);
283        }
284        variance
285    }
286
287    /// Executes a hardcoded 24-hour Monte Carlo simulation of the stochastic model, exporting the time history to a Parquet file.
288    ///
289    /// # Warning: Hardcoded Time Series & Diagnostic Data Gaps
290    /// This method does *not* accept a user-defined tracking schedule or time series. It inherently evaluates the stochastic process
291    /// over a strict 24-hour period, beginning at the exact system clock moment of method execution, utilizing a 1-minute step size.
292    ///
293    /// Furthermore, users will observe exactly 1,082 samples per simulation run, rather than the 1,441 samples expected from a
294    /// continuous 24-hour 1-minute cadence. The simulation intentionally drops all epochs strictly greater than +6 hours and
295    /// strictly less than +12 hours from the start time. This hardcoded artifact is designed to demonstrate variance bounds
296    /// expansion in the absence of measurements (e.g., simulating a tracking dropout for a Gauss-Markov bias).
297    ///
298    /// # Algorithm
299    /// 1. Establish `start` as the system clock time at invocation.
300    /// 2. Construct an inclusive time series from `start` to `start + 24 hours` at 1-minute intervals.
301    /// 3. For each configured run, seed a PRNG (`Pcg64Mcg`) using system entropy.
302    /// 4. Evaluate the process covariance and sample the stochastic noise at each epoch.
303    /// 5. Discard all epochs inside the `(start + 6h, start + 12h)` open interval.
304    /// 6. Export the remaining 1,082 samples per run to an Apache Arrow RecordBatch and write to disk via Parquet.
305    ///
306    /// :param path: The filesystem path for the output Parquet file.
307    /// :type path: str
308    /// :param runs: The number of Monte Carlo runs. Defaults to 25 if not provided.
309    /// :type runs: int | None
310    /// :param unit: An optional string appended to the Parquet column headers for plotting clarity.
311    /// :type unit: str | None
312    /// :rtype: list[StochasticState]
313    /// :raises Exception: If the underlying Apache Arrow RecordBatch fails to allocate or write to the specified filesystem path.
314    #[cfg(feature = "python")]
315    #[pyo3(name = "simulate")]
316    fn py_simulate(
317        &self,
318        path: &str,
319        runs: Option<u32>,
320        unit: Option<String>,
321    ) -> PyResult<Vec<StochasticState>> {
322        self.simulate(path, runs, unit)
323            .map_err(|e| PyValueError::new_err(e.to_string()))
324    }
325
326    /// Constructs a high precision zero-mean range noise model (accounting for clock error and thermal error) from
327    /// the Allan deviation of the clock, integration time, chip rate (depends on the ranging code), and
328    /// signal-power-to-noise-density ratio (S/N₀).
329    ///
330    /// NOTE: The Allan Deviation should be provided given the integration time. For example, if the integration time
331    /// is one second, the Allan Deviation should be the deviation over one second.
332    ///
333    /// IMPORTANT: These do NOT include atmospheric noises, which add up to ~10 cm one-sigma.
334    #[cfg(feature = "python")]
335    #[pyo3(name = "from_hardware_range_km")]
336    #[classmethod]
337    fn py_from_hardware_range_km(
338        _cls: &Bound<'_, PyType>,
339
340        allan_deviation: f64,
341        integration_time: Duration,
342        chip_rate: link_specific::ChipRate,
343        s_n0: link_specific::SN0,
344    ) -> Self {
345        Self::from_hardware_range_km(allan_deviation, integration_time, chip_rate, s_n0)
346    }
347
348    #[cfg(feature = "python")]
349    #[pyo3(name = "from_hardware_doppler_km_s")]
350    #[classmethod]
351    fn py_from_hardware_doppler_km_s(
352        _cls: &Bound<'_, PyType>,
353        allan_deviation: f64,
354        integration_time: Duration,
355        carrier: link_specific::CarrierFreq,
356        c_n0: link_specific::CN0,
357    ) -> Self {
358        Self::from_hardware_doppler_km_s(allan_deviation, integration_time, carrier, c_n0)
359    }
360
361    fn __str__(&self) -> String {
362        format!("{self}")
363    }
364
365    fn __repr__(&self) -> String {
366        format!("{self} @ {self:p}")
367    }
368}
369
370impl Display for StochasticNoise {
371    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
372        match (self.white_noise, self.bias) {
373            (Some(wn), None) => write!(f, "Stochastics with {wn:?}"),
374            (None, Some(bias)) => write!(f, "Stochastics with bias {bias}"),
375            (None, None) => write!(f, "Noiseless stochastics"),
376            (Some(wn), Some(bias)) => write!(f, "Stochastics with {wn:?} and bias {bias}"),
377        }
378    }
379}
380
381impl Mul<f64> for StochasticNoise {
382    type Output = Self;
383
384    fn mul(mut self, rhs: f64) -> Self::Output {
385        if let Some(wn) = &mut self.white_noise {
386            *wn *= rhs;
387        }
388        if let Some(gm) = &mut self.bias {
389            *gm *= rhs;
390        }
391
392        self
393    }
394}
395
396impl MulAssign<f64> for StochasticNoise {
397    fn mul_assign(&mut self, rhs: f64) {
398        *self = *self * rhs;
399    }
400}
401
402impl Encode for StochasticNoise {
403    fn encoded_len(&self) -> der::Result<der::Length> {
404        let flags = self.available_data();
405        flags.encoded_len()? + self.white_noise.encoded_len()? + self.bias.encoded_len()?
406    }
407
408    fn encode(&self, encoder: &mut impl der::Writer) -> der::Result<()> {
409        let flags = self.available_data();
410
411        flags.encode(encoder)?;
412        self.white_noise.encode(encoder)?;
413        self.bias.encode(encoder)
414    }
415}
416
417impl<'a> Decode<'a> for StochasticNoise {
418    fn decode<R: Reader<'a>>(decoder: &mut R) -> der::Result<Self> {
419        let flags: u8 = decoder.decode()?;
420
421        let white_noise = if flags & (1 << 0) != 0 {
422            Some(decoder.decode()?)
423        } else {
424            None
425        };
426
427        let bias = if flags & (1 << 1) != 0 {
428            Some(decoder.decode()?)
429        } else {
430            None
431        };
432
433        Ok(Self { white_noise, bias })
434    }
435}
436
437#[derive(Copy, Clone, Debug)]
438#[cfg_attr(feature = "python", pyclass(from_py_object, get_all))]
439pub struct StochasticState {
440    pub run: u32,
441    pub dt_s: f64,
442    pub sample: f64,
443    pub variance: f64,
444}
445
446#[cfg(feature = "python")]
447#[cfg_attr(feature = "python", pymethods)]
448impl StochasticState {
449    fn __str__(&self) -> String {
450        format!("{self:?}")
451    }
452    fn __repr__(&self) -> String {
453        format!("{self:?} @ {self:p}")
454    }
455}
456
457#[cfg(test)]
458mod ut_stochastics {
459    use std::path::PathBuf;
460
461    use super::{StochasticNoise, white::WhiteNoise};
462
463    #[test]
464    fn test_simulate_zero() {
465        let path: PathBuf = [
466            env!("CARGO_MANIFEST_DIR"),
467            "../data",
468            "04_output",
469            "stochastics_zero.parquet",
470        ]
471        .iter()
472        .collect();
473
474        let noise = StochasticNoise::default();
475
476        let rslts = noise.simulate(path, None, None).unwrap();
477        assert!(!rslts.is_empty());
478        assert!(rslts.iter().map(|rslt| rslt.sample).sum::<f64>().abs() < f64::EPSILON);
479    }
480
481    #[test]
482    fn test_simulate_constant() {
483        let path: PathBuf = [
484            env!("CARGO_MANIFEST_DIR"),
485            "../data",
486            "04_output",
487            "stochastics_constant.parquet",
488        ]
489        .iter()
490        .collect();
491
492        let noise = StochasticNoise {
493            white_noise: Some(WhiteNoise {
494                mean: 15.0,
495                sigma: 2.0,
496            }),
497            ..Default::default()
498        };
499
500        noise.simulate(path, None, None).unwrap();
501    }
502
503    #[test]
504    fn test_simulate_dsn_range() {
505        let path: PathBuf = [
506            env!("CARGO_MANIFEST_DIR"),
507            "../data",
508            "04_output",
509            "stochastics_dsn_range.parquet",
510        ]
511        .iter()
512        .collect();
513
514        let noise = StochasticNoise::default_range_km();
515
516        noise
517            .simulate(path, None, Some("kilometer".to_string()))
518            .unwrap();
519    }
520
521    #[test]
522    fn test_simulate_dsn_range_gm_only() {
523        let path: PathBuf = [
524            env!("CARGO_MANIFEST_DIR"),
525            "../data",
526            "04_output",
527            "stochastics_dsn_range_gm_only.parquet",
528        ]
529        .iter()
530        .collect();
531
532        let mut noise = StochasticNoise::default_range_km();
533        noise.white_noise = None;
534
535        noise
536            .simulate(path, None, Some("kilometer".to_string()))
537            .unwrap();
538    }
539}
nyx_space/od/noise/mod.rs

nyx_space/od/noise/
mod.rs
