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

nyx_space/od/noise/
mod.rs
