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;
26use rand::rngs::SysRng;
27use rand::{Rng, SeedableRng};
28use rand_pcg::Pcg64Mcg;
29use serde::{Deserialize, Serialize};
30use std::error::Error;
31use std::fs::File;
32use std::ops::{Mul, MulAssign};
33use std::path::Path;
34use std::sync::Arc;
35
36pub mod gauss_markov;
37#[cfg(feature = "premium")]
38pub mod link_specific;
39pub mod white;
40
41#[cfg(feature = "python")]
42use pyo3::prelude::*;
43
44pub use gauss_markov::GaussMarkov;
45pub use white::WhiteNoise;
46
47/// Trait for any kind of stochastic modeling, developing primarily for synthetic orbit determination measurements.
48pub trait Stochastics {
49    /// Return the variance of this stochastic noise model at a given time.
50    fn covariance(&self, epoch: Epoch) -> f64;
51
52    /// Returns a new sample of these stochastics
53    fn sample<R: Rng>(&mut self, epoch: Epoch, rng: &mut R) -> f64;
54}
55
56/// Stochastic noise modeling used primarily for synthetic orbit determination measurements.
57///
58/// This implementation distinguishes between the white noise model and the bias model. It also includes a constant offset.
59#[derive(Copy, Clone, Debug, Default, PartialEq, Serialize, Deserialize)]
60#[cfg_attr(feature = "python", pyclass(get_all, set_all))]
61pub struct StochasticNoise {
62    pub white_noise: Option<WhiteNoise>,
63    pub bias: Option<GaussMarkov>,
64}
65
66impl StochasticNoise {
67    /// Zero noise stochastic process.
68    pub const ZERO: Self = Self {
69        white_noise: None,
70        bias: None,
71    };
72
73    /// The minimum stochastic noise process with a zero mean white noise of 1e-6.
74    pub const MIN: Self = Self {
75        white_noise: Some(WhiteNoise {
76            mean: 0.0,
77            sigma: 1e-6,
78        }),
79        bias: None,
80    };
81
82    /// Default stochastic process of the Deep Space Network, as per DESCANSO Chapter 3, Table 3-3.
83    /// Using the instrument bias as the white noise value, zero constant bias.
84    pub fn default_range_km() -> Self {
85        Self {
86            white_noise: Some(WhiteNoise {
87                sigma: 2.0e-3, // 2 m
88                ..Default::default()
89            }),
90            // Until Nyx can support bias estimation the default bias is None, cf. https://github.com/nyx-space/nyx/issues/326
91            // bias: Some(GaussMarkov::default_range_km()),
92            bias: None,
93        }
94    }
95
96    /// Default stochastic process of the Deep Space Network, using as per DESCANSO Chapter 3, Table 3-3 for the GM process.
97    pub fn default_doppler_km_s() -> Self {
98        Self {
99            white_noise: Some(WhiteNoise {
100                sigma: 3e-6, // 3 mm/s
101                ..Default::default()
102            }),
103            // Until Nyx can support bias estimation the default bias is None, cf. https://github.com/nyx-space/nyx/issues/326
104            // bias: Some(GaussMarkov::default_doppler_km_s()),
105            bias: None,
106        }
107    }
108
109    /// Default stochastic process for an angle measurement (azimuth or elevation)
110    /// Using the instrument bias as the white noise value, zero constant bias.
111    pub fn default_angle_deg() -> Self {
112        Self {
113            white_noise: Some(WhiteNoise {
114                sigma: 1.0e-2, // 0.01 deg
115                ..Default::default()
116            }),
117            // Until Nyx can support bias estimation the default bias is None, cf. https://github.com/nyx-space/nyx/issues/326
118            // bias: Some(GaussMarkov::default_range_km()),
119            bias: None,
120        }
121    }
122
123    /// Sample these stochastics
124    pub fn sample<R: Rng>(&mut self, epoch: Epoch, rng: &mut R) -> f64 {
125        let mut sample = 0.0;
126        if let Some(wn) = &mut self.white_noise {
127            sample += wn.sample(epoch, rng)
128        }
129        if let Some(gm) = &mut self.bias {
130            sample += gm.sample(epoch, rng);
131        }
132        sample
133    }
134
135    /// Return the covariance of these stochastics at a given time.
136    pub fn covariance(&self, epoch: Epoch) -> f64 {
137        let mut variance = 0.0;
138        if let Some(wn) = &self.white_noise {
139            variance += wn.covariance(epoch);
140        }
141        if let Some(gm) = &self.bias {
142            variance += gm.covariance(epoch);
143        }
144        variance
145    }
146
147    /// Simulate the configured stochastic model and store the bias in a parquet file.
148    /// Python: call as `simulate(path, runs=25, unit=None)` where the path is the output Parquet file, runs is the number of runs, and unit is the unit of the bias, reflected only in the headers of the parquet file.
149    ///
150    /// The unit is only used in the headers of the parquet file.
151    ///
152    /// This will simulate the model with "runs" different seeds, sampling the process 500 times for a duration of 5 times the time constant.
153    pub fn simulate<P: AsRef<Path>>(
154        self,
155        path: P,
156        runs: Option<u32>,
157        unit: Option<String>,
158    ) -> Result<Vec<StochasticState>, Box<dyn Error>> {
159        let num_runs = runs.unwrap_or(25);
160
161        let start = Epoch::now().unwrap();
162        let (step, end) = (1.minutes(), start + 1.days());
163
164        let capacity = ((end - start).to_seconds() / step.to_seconds()).ceil() as usize;
165
166        let mut samples = Vec::with_capacity(capacity);
167
168        for run in 0..num_runs {
169            let mut rng = Pcg64Mcg::try_from_rng(&mut SysRng).unwrap();
170
171            let mut mdl = self;
172            for epoch in TimeSeries::inclusive(start, end, step) {
173                if epoch > start + 6.hours() && epoch < start + 12.hours() {
174                    // Skip to see how the variance changes.
175                    continue;
176                }
177                let variance = mdl.covariance(epoch);
178                let sample = mdl.sample(epoch, &mut rng);
179                samples.push(StochasticState {
180                    run,
181                    dt_s: (epoch - start).to_seconds(),
182                    sample,
183                    variance,
184                });
185            }
186        }
187
188        let bias_unit = match unit {
189            Some(unit) => format!("({unit})"),
190            None => "(unitless)".to_string(),
191        };
192
193        // Build the parquet file
194        let hdrs = vec![
195            Field::new("Run", DataType::UInt32, false),
196            Field::new("Delta Time (s)", DataType::Float64, false),
197            Field::new(format!("Bias {bias_unit}"), DataType::Float64, false),
198            Field::new(format!("Variance {bias_unit}"), DataType::Float64, false),
199        ];
200
201        let schema = Arc::new(Schema::new(hdrs));
202        let record = vec![
203            Arc::new(UInt32Array::from(
204                samples.iter().map(|s| s.run).collect::<Vec<u32>>(),
205            )) as ArrayRef,
206            Arc::new(Float64Array::from(
207                samples.iter().map(|s| s.dt_s).collect::<Vec<f64>>(),
208            )) as ArrayRef,
209            Arc::new(Float64Array::from(
210                samples.iter().map(|s| s.sample).collect::<Vec<f64>>(),
211            )) as ArrayRef,
212            Arc::new(Float64Array::from(
213                samples.iter().map(|s| s.variance).collect::<Vec<f64>>(),
214            )) as ArrayRef,
215        ];
216
217        let props = pq_writer(None);
218
219        let file = File::create(path)?;
220        let mut writer = ArrowWriter::try_new(file, schema.clone(), props).unwrap();
221
222        let batch = RecordBatch::try_new(schema, record)?;
223        writer.write(&batch)?;
224        writer.close()?;
225
226        Ok(samples)
227    }
228
229    fn available_data(&self) -> u8 {
230        let mut bits: u8 = 0;
231
232        if self.white_noise.is_some() {
233            bits |= 1 << 0;
234        }
235
236        if self.bias.is_some() {
237            bits |= 1 << 1;
238        }
239
240        bits
241    }
242}
243
244impl Mul<f64> for StochasticNoise {
245    type Output = Self;
246
247    fn mul(mut self, rhs: f64) -> Self::Output {
248        if let Some(mut wn) = &mut self.white_noise {
249            wn *= rhs;
250        }
251        if let Some(mut gm) = &mut self.bias {
252            gm *= rhs;
253        }
254
255        self
256    }
257}
258
259impl MulAssign<f64> for StochasticNoise {
260    fn mul_assign(&mut self, rhs: f64) {
261        *self = *self * rhs;
262    }
263}
264
265impl Encode for StochasticNoise {
266    fn encoded_len(&self) -> der::Result<der::Length> {
267        let flags = self.available_data();
268        flags.encoded_len()? + self.white_noise.encoded_len()? + self.bias.encoded_len()?
269    }
270
271    fn encode(&self, encoder: &mut impl der::Writer) -> der::Result<()> {
272        let flags = self.available_data();
273
274        flags.encode(encoder)?;
275        self.white_noise.encode(encoder)?;
276        self.bias.encode(encoder)
277    }
278}
279
280impl<'a> Decode<'a> for StochasticNoise {
281    fn decode<R: Reader<'a>>(decoder: &mut R) -> der::Result<Self> {
282        let flags: u8 = decoder.decode()?;
283
284        let white_noise = if flags & (1 << 0) != 0 {
285            Some(decoder.decode()?)
286        } else {
287            None
288        };
289
290        let bias = if flags & (1 << 1) != 0 {
291            Some(decoder.decode()?)
292        } else {
293            None
294        };
295
296        Ok(Self { white_noise, bias })
297    }
298}
299
300pub struct StochasticState {
301    pub run: u32,
302    pub dt_s: f64,
303    pub sample: f64,
304    pub variance: f64,
305}
306
307#[cfg(test)]
308mod ut_stochastics {
309    use std::path::PathBuf;
310
311    use super::{white::WhiteNoise, StochasticNoise};
312
313    #[test]
314    fn test_simulate_zero() {
315        let path: PathBuf = [
316            env!("CARGO_MANIFEST_DIR"),
317            "data",
318            "04_output",
319            "stochastics_zero.parquet",
320        ]
321        .iter()
322        .collect();
323
324        let noise = StochasticNoise::default();
325
326        let rslts = noise.simulate(path, None, None).unwrap();
327        assert!(!rslts.is_empty());
328        assert!(rslts.iter().map(|rslt| rslt.sample).sum::<f64>().abs() < f64::EPSILON);
329    }
330
331    #[test]
332    fn test_simulate_constant() {
333        let path: PathBuf = [
334            env!("CARGO_MANIFEST_DIR"),
335            "data",
336            "04_output",
337            "stochastics_constant.parquet",
338        ]
339        .iter()
340        .collect();
341
342        let noise = StochasticNoise {
343            white_noise: Some(WhiteNoise {
344                mean: 15.0,
345                sigma: 2.0,
346            }),
347            ..Default::default()
348        };
349
350        noise.simulate(path, None, None).unwrap();
351    }
352
353    #[test]
354    fn test_simulate_dsn_range() {
355        let path: PathBuf = [
356            env!("CARGO_MANIFEST_DIR"),
357            "data",
358            "04_output",
359            "stochastics_dsn_range.parquet",
360        ]
361        .iter()
362        .collect();
363
364        let noise = StochasticNoise::default_range_km();
365
366        noise
367            .simulate(path, None, Some("kilometer".to_string()))
368            .unwrap();
369    }
370
371    #[test]
372    fn test_simulate_dsn_range_gm_only() {
373        let path: PathBuf = [
374            env!("CARGO_MANIFEST_DIR"),
375            "data",
376            "04_output",
377            "stochastics_dsn_range_gm_only.parquet",
378        ]
379        .iter()
380        .collect();
381
382        let mut noise = StochasticNoise::default_range_km();
383        noise.white_noise = None;
384
385        noise
386            .simulate(path, None, Some("kilometer".to_string()))
387            .unwrap();
388    }
389}
nyx_space/od/noise/mod.rs

nyx_space/od/noise/
mod.rs
