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

nyx_space/od/noise/
mod.rs
