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

nyx_space/od/noise/
mod.rs
