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

nyx_space/od/noise/
mod.rs
