gmt_lom/
lib.rs

1//! # Giant Magellan Telescope Linear Optical Model
2//!
3//! Optical linear transformations applied to the rigid body motions of the primary and secondary segmented mirrors of the GMT
4//!
5//! The optical sensitivities can be downloaded from [here](https://s3.us-west-2.amazonaws.com/gmto.modeling/optical_sensitivities.rs.bin),
6//! or they can be recomputed with the [makesens](../makesens/index.html) binary compiled with the `crseo` features
7//! and run on a computer with a NVIDIA GPU.
8//!
9//! # Example
10//! ```
11//! use std::iter::once;
12//! use skyangle::Conversion;
13//! use gmt_lom::LOM;
14//!
15//! let m1_rbm = vec![vec![0f64; 6]; 7];
16//! let mut m2_rbm = vec![vec![0f64; 6]; 7];
17//! m2_rbm
18//!     .iter_mut()
19//!     .step_by(2)
20//!     .for_each(|str| {
21//!        str[3] = 1f64.from_arcsec();
22//!        str[4] = 1f64.from_arcsec();
23//!  });
24//! m2_rbm[6][3] = 1f64.from_arcsec();
25//! m2_rbm[6][4] = 1f64.from_arcsec();
26//! let lom = LOM::builder()
27//!     .into_iter_rigid_body_motions(once((m1_rbm, m2_rbm)))
28//!     .build()
29//!     .unwrap();
30//! let tt = lom.tiptilt();
31//! println!(" Tiptilt: {:.0?}mas", tt);
32//! let stt = lom.segment_tiptilt();
33//! println!("Segment tiptilt:");
34//! println!(" - x: {:.0?}mas", &stt[..7]);
35//! println!(" - y: {:.0?}mas", &stt[7..]);
36//! let sp = lom.segment_piston();
37//! println!("Segment piston:");
38//! sp.chunks(7).enumerate().for_each(|(k, sp)| println!(" - S{}: {:.0?}nm", k + 1, sp) );
39//! ```
40
41use bincode;
42use serde::Serialize;
43use serde_pickle as pickle;
44
45use std::{
46    env,
47    fs::File,
48    marker::PhantomData,
49    ops::{Deref, DerefMut},
50    path::{Path, PathBuf},
51    slice::Chunks,
52};
53
54pub mod lom;
55pub use lom::{LOMBuilder, LOM};
56mod optical_sensitivities;
57pub use optical_sensitivities::{from_opticals, OpticalSensitivities, OpticalSensitivity};
58mod rigid_body_motions;
59pub use rigid_body_motions::RigidBodyMotions;
60#[cfg(feature = "apache")]
61mod table;
62#[cfg(feature = "apache")]
63pub use table::Table;
64
65use crate::rigid_body_motions::RigidBodyMotionsError;
66// pub mod actors_interface;
67
68#[derive(thiserror::Error, Debug)]
69pub enum LinearOpticalModelError {
70    #[error("sensitivities file not found (optical_sensitivities.rs.bin)")]
71    SensitivityFile(#[from] std::io::Error),
72    #[error("parquet file: {1:?} not found")]
73    ParquetFile(#[source] std::io::Error, PathBuf),
74    #[error("sensitivities cannot be loaded from optical_sensitivities.rs.bin")]
75    SensitivityData(#[from] bincode::Error),
76    #[error("segment tip-tilt sensitivity is missing")]
77    SegmentTipTilt,
78    #[error("rigid body motions are missing")]
79    MissingRigidBodyMotions,
80    #[error("failed to write optical metric to pickle file ")]
81    MetricsPickleFile(#[from] pickle::Error),
82    #[error("missing table {0} column ")]
83    Table(String),
84    #[error("PMT read failed")]
85    Read(#[from] csv::Error),
86    #[cfg(feature = "apache")]
87    #[error("failed to read parquet Table")]
88    TableRead(#[from] table::TableError),
89    #[error("failed to process rigid body motions")]
90    RigidBodyMotions(#[from] RigidBodyMotionsError),
91}
92type Result<T> = std::result::Result<T, LinearOpticalModelError>;
93
94#[derive(Debug, Clone)]
95pub enum Formatting {
96    AdHoc,
97    Latex,
98}
99
100/// Sensitivities serialization into a [bincode] file
101pub trait Bin {
102    fn dump<P: AsRef<Path>>(self, path: P) -> Result<Self>
103    where
104        Self: Sized;
105    fn load<P: AsRef<Path>>(path: P) -> Result<Self>
106    where
107        Self: Sized;
108}
109impl<const N: usize> Bin for OpticalSensitivities<N> {
110    /// Saves sensitivities to `path`
111    fn dump<P: AsRef<Path>>(self, path: P) -> Result<Self> {
112        bincode::serialize_into(File::create(path)?, &self)?;
113        Ok(self)
114    }
115    /// Load sensitivities from `path`
116    fn load<P: AsRef<Path>>(path: P) -> Result<Self> {
117        Ok(bincode::deserialize_from(File::open(path)?)?)
118    }
119}
120
121/// Data loader
122pub struct Loader<T> {
123    path: PathBuf,
124    filename: String,
125    phantom: PhantomData<T>,
126}
127/// [Loader] loading interface
128pub trait LoaderTrait<T> {
129    fn load(self) -> Result<T>;
130}
131impl<T> Loader<T> {
132    /// Set the loading path
133    pub fn path<P: AsRef<Path>>(self, path: P) -> Self {
134        Self {
135            path: Path::new(path.as_ref()).to_path_buf(),
136            ..self
137        }
138    }
139    /// Set the loaded file name
140    pub fn filename(self, filename: &str) -> Self {
141        Self {
142            filename: String::from(filename),
143            ..self
144        }
145    }
146}
147impl<const N: usize> Default for Loader<OpticalSensitivities<N>> {
148    /// Default [Loader] for [Vec] of [OpticalSensitivity],
149    /// expecting the file `optical_sensitivities.rs.bin` in the current folder
150    fn default() -> Self {
151        let path = env::var("LOM").unwrap_or_else(|_| ".".to_string());
152        Self {
153            path: Path::new(&path).to_path_buf(),
154            filename: String::from("optical_sensitivities.rs.bin"),
155            phantom: PhantomData,
156        }
157    }
158}
159impl<const N: usize> LoaderTrait<OpticalSensitivities<N>> for Loader<OpticalSensitivities<N>> {
160    /// Loads precomputed optical sensitivities
161    fn load(self) -> Result<OpticalSensitivities<N>> {
162        log::info!("Loading optical sensitivities ...");
163        <OpticalSensitivities<N> as Bin>::load(self.path.join(self.filename))
164    }
165}
166#[cfg(feature = "apache")]
167impl Default for Loader<RigidBodyMotions> {
168    /// Default [Loader] for [RigidBodyMotions] expecting the file `data.parquet` in the current folder
169    fn default() -> Self {
170        Self {
171            path: Path::new(".").to_path_buf(),
172            filename: String::from("data.parquet"),
173            phantom: PhantomData,
174        }
175    }
176}
177#[cfg(feature = "apache")]
178impl LoaderTrait<RigidBodyMotions> for Loader<RigidBodyMotions> {
179    /// Loads M1 and M2 rigid body motions
180    fn load(self) -> Result<RigidBodyMotions> {
181        log::info!("Loading rigid body motions ...");
182        RigidBodyMotions::from_parquet(self.path.join(self.filename), None, None)
183    }
184}
185
186/// Type holding the tip-tilt values
187#[derive(Serialize, Debug, Clone)]
188pub struct TipTilt(Vec<f64>);
189/// Type holding the segment tip-tilt values
190#[derive(Serialize, Debug, Clone)]
191pub struct SegmentTipTilt(Vec<f64>);
192/// Type holding the segment piston values
193#[derive(Serialize, Debug, Clone)]
194pub struct SegmentPiston(Vec<f64>);
195// Dereferencing
196impl Deref for TipTilt {
197    type Target = Vec<f64>;
198    fn deref(&self) -> &Self::Target {
199        &self.0
200    }
201}
202impl DerefMut for TipTilt {
203    fn deref_mut(&mut self) -> &mut Self::Target {
204        &mut self.0
205    }
206}
207impl Deref for SegmentTipTilt {
208    type Target = Vec<f64>;
209    fn deref(&self) -> &Self::Target {
210        &self.0
211    }
212}
213impl DerefMut for SegmentTipTilt {
214    fn deref_mut(&mut self) -> &mut Self::Target {
215        &mut self.0
216    }
217}
218impl Deref for SegmentPiston {
219    type Target = Vec<f64>;
220    fn deref(&self) -> &Self::Target {
221        &self.0
222    }
223}
224impl DerefMut for SegmentPiston {
225    fn deref_mut(&mut self) -> &mut Self::Target {
226        &mut self.0
227    }
228}
229impl From<TipTilt> for Vec<f64> {
230    fn from(value: TipTilt) -> Self {
231        value.0
232    }
233}
234impl From<SegmentPiston> for Vec<f64> {
235    fn from(value: SegmentPiston) -> Self {
236        value.0
237    }
238}
239impl From<SegmentTipTilt> for Vec<f64> {
240    fn from(value: SegmentTipTilt) -> Self {
241        value.0
242    }
243}
244
245pub trait ToPkl {
246    /// Writes optical metrics to a [pickle] file
247    fn to_pkl<P: AsRef<Path>>(&self, path: P) -> Result<()>
248    where
249        Self: Serialize + Sized,
250    {
251        let mut file = File::create(&path)?;
252        pickle::ser::to_writer(&mut file, self, pickle::ser::SerOptions::new())?;
253        Ok(())
254    }
255}
256impl ToPkl for TipTilt {}
257impl ToPkl for SegmentTipTilt {}
258impl ToPkl for SegmentPiston {}
259
260/// Trait for the [LOM] optical metrics
261///
262/// A simple trait looking at the number of items in the [TipTilt], [SegmentTipTilt] and [SegmentPiston] metrics
263pub trait OpticalMetrics {
264    fn n_item(&self) -> usize;
265    /// Returns a [Chunks] iterator with chunks the size of [n_item](OpticalMetrics::n_item)
266    fn items(&self) -> Chunks<'_, f64>
267    where
268        Self: Deref<Target = Vec<f64>>,
269    {
270        self.deref().chunks(self.n_item())
271    }
272    /// Returns the metrics assigning each component in a contiguous time vector
273    fn time_wise(&self, n_sample: Option<usize>) -> Vec<f64>;
274}
275impl OpticalMetrics for TipTilt {
276    /// [TipTilt] `2` x and y items
277    fn n_item(&self) -> usize {
278        2
279    }
280    fn time_wise(&self, n_sample: Option<usize>) -> Vec<f64> {
281        let n_item = self.n_item();
282        let n_total = self.len() / n_item;
283        assert!(n_total >= n_sample.unwrap_or(n_total), "not enough samples");
284        (0..n_item)
285            .flat_map(|i| {
286                self.iter()
287                    .skip(i)
288                    .step_by(n_item)
289                    .skip(n_total - n_sample.unwrap_or(n_total))
290            })
291            .cloned()
292            .collect()
293    }
294}
295impl OpticalMetrics for SegmentTipTilt {
296    /// [SegmentTipTilt] `7x2` x and y items
297    fn n_item(&self) -> usize {
298        14
299    }
300    fn time_wise(&self, n_sample: Option<usize>) -> Vec<f64> {
301        let n_item = self.n_item();
302        let n_total = self.len() / n_item;
303        assert!(n_total >= n_sample.unwrap_or(n_total), "not enough samples");
304        (0..n_item)
305            .flat_map(|i| {
306                self.iter()
307                    .skip(i)
308                    .step_by(n_item)
309                    .skip(n_total - n_sample.unwrap_or(n_total))
310            })
311            .cloned()
312            .collect()
313    }
314}
315impl OpticalMetrics for SegmentPiston {
316    /// [SegmentPiston] `7` items
317    fn n_item(&self) -> usize {
318        7
319    }
320    fn time_wise(&self, n_sample: Option<usize>) -> Vec<f64> {
321        let n_item = self.n_item();
322        let n_total = self.len() / n_item;
323        assert!(n_total >= n_sample.unwrap_or(n_total), "not enough samples");
324        (0..n_item)
325            .flat_map(|i| {
326                self.iter()
327                    .skip(i)
328                    .step_by(n_item)
329                    .skip(n_total - n_sample.unwrap_or(n_total))
330            })
331            .cloned()
332            .collect()
333    }
334}
335
336/// Statistics on [OpticalMetrics]
337pub trait Stats {
338    /// Returns the mean values
339    ///
340    /// Optionally, the statistical moment is evaluated on the last `n_sample`
341    fn mean(&self, n_sample: Option<usize>) -> Vec<f64>
342    where
343        Self: Deref<Target = Vec<f64>> + OpticalMetrics,
344    {
345        let n_item = self.n_item();
346        let n_total = self.len() / n_item;
347        assert!(n_total >= n_sample.unwrap_or(n_total), "not enough samples");
348        (0..n_item)
349            .map(|i| {
350                self.iter()
351                    .skip(i)
352                    .step_by(n_item)
353                    .skip(n_total - n_sample.unwrap_or(n_total))
354                    .sum::<f64>()
355                    / n_sample.unwrap_or(n_total) as f64
356            })
357            .collect()
358    }
359    /// Returns the mean variance values
360    ///
361    /// Optionally, the statistical moment is evaluated on the last `n_sample`
362    fn var(&self, n_sample: Option<usize>) -> Vec<f64>
363    where
364        Self: Deref<Target = Vec<f64>> + OpticalMetrics,
365    {
366        let n_item = self.n_item();
367        let n_total = self.len() / n_item;
368        assert!(n_total >= n_sample.unwrap_or(n_total), "not enough samples");
369        (0..n_item)
370            .zip(self.mean(n_sample))
371            .map(|(i, m)| {
372                self.iter()
373                    .skip(i)
374                    .step_by(n_item)
375                    .skip(n_total - n_sample.unwrap_or(n_total))
376                    .map(|&x| x - m)
377                    .fold(0f64, |a, x| a + x * x)
378                    / n_sample.unwrap_or(n_total) as f64
379            })
380            .collect()
381    }
382    /// Returns the standard deviation values
383    ///
384    /// Optionally, the statistical moment is evaluated on the last `n_sample`
385    fn std(&self, n_sample: Option<usize>) -> Vec<f64>
386    where
387        Self: Deref<Target = Vec<f64>> + OpticalMetrics,
388    {
389        self.var(n_sample).iter().map(|x| x.sqrt()).collect()
390    }
391}
392impl Stats for TipTilt {}
393impl Stats for SegmentTipTilt {}
394impl Stats for SegmentPiston {}
395
396#[cfg(test)]
397mod tests {
398    use super::*;
399    use skyangle::Conversion;
400    use std::iter::once;
401
402    #[test]
403    fn tiptilt() {
404        let mut m1_rbm = vec![vec![0f64; 6]; 7];
405        let mut m2_rbm = vec![vec![0f64; 6]; 7];
406        m1_rbm
407            .iter_mut()
408            .step_by(2)
409            .for_each(|str| str[3] = 1f64.from_arcsec());
410        m2_rbm
411            .iter_mut()
412            .skip(1)
413            .step_by(2)
414            .take(3)
415            .for_each(|str| str[3] = 1f64.from_arcsec());
416        m1_rbm[6][3] = 0.5f64.from_arcsec();
417        m2_rbm[6][3] = (-4f64).from_arcsec();
418        let lom = LOM::builder()
419            .into_iter_rigid_body_motions(once((m1_rbm, m2_rbm)))
420            .build()
421            .unwrap();
422        let stt = lom.segment_tiptilt();
423        let mag: Vec<_> = stt[..7]
424            .iter()
425            .zip(&stt[7..])
426            .map(|(x, y)| x.hypot(*y).to_arcsec())
427            .collect();
428        print!("Segment tiptilt : {:.3?} mas", mag);
429    }
430
431    #[test]
432    fn segment_wfe_rms() {
433        let mut m1_rbm = vec![vec![0f64; 6]; 7];
434        let mut m2_rbm = vec![vec![0f64; 6]; 7];
435        for i in 0..6 {
436            m1_rbm[i][2] = 100e-9 * (i + 1) as f64;
437        }
438        m2_rbm[6][2] = -100e-9;
439        let lom = LOM::builder()
440            .into_iter_rigid_body_motions(once((m1_rbm, m2_rbm)))
441            .build()
442            .unwrap();
443        let swferms = lom.segment_wfe_rms::<-9>();
444        print!("Segment WFE RMS : {:.0?} mas", swferms);
445    }
446}
gmt_lom/lib.rs

gmt_lom/
lib.rs
