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

gmt_lom/
lib.rs
