treams_core/tmatrix/sphere/
saved.rs

1//! Fixed derivative state: shape, physical radii and complete per-degree Mie states.
2//! treams-rs extension.
3
4use super::SphereResidual;
5use crate::{
6    MAX_DEGREE, Result,
7    coeffs::MieResidual,
8    saved::{Reader, SavedState, Writer, invalid},
9};
10
11impl SphereResidual {
12    pub(crate) fn lmax(&self) -> usize {
13        self.degrees.len()
14    }
15
16    /// Byte count of saved derivative state for this degree and layer count.
17    pub fn state_size(lmax: usize, boundaries: usize) -> Result<usize> {
18        if lmax == 0 || lmax > MAX_DEGREE.unsigned_abs() as usize {
19            return Err(invalid());
20        }
21        let degrees = MieResidual::state_size(boundaries)?
22            .checked_mul(lmax)
23            .ok_or_else(invalid)?;
24        boundaries
25            .checked_mul(8)
26            .and_then(|radii| radii.checked_add(24))
27            .and_then(|header| header.checked_add(degrees))
28            .ok_or_else(invalid)
29    }
30
31    pub(crate) fn write_state(&self, writer: &mut Writer) {
32        writer.usize(self.degrees.len());
33        writer.usize(self.radii.len());
34        writer.f64(self.k0);
35        for &radius in &self.radii {
36            writer.f64(radius);
37        }
38        for degree in &self.degrees {
39            degree.write_state(writer);
40        }
41    }
42
43    pub(crate) fn read_state(reader: &mut Reader<'_>) -> Result<Self> {
44        let available = reader.remaining_len();
45        let lmax = reader.usize()?;
46        let boundaries = reader.usize()?;
47        if Self::state_size(lmax, boundaries)? > available {
48            return Err(invalid());
49        }
50        let k0 = reader.f64()?;
51        let radii = (0..boundaries)
52            .map(|_| reader.f64())
53            .collect::<Result<Vec<_>>>()?;
54        let degrees = (0..lmax)
55            .map(|_| {
56                let degree = MieResidual::read_state(reader)?;
57                if degree.boundaries() != boundaries {
58                    return Err(invalid());
59                }
60                Ok(degree)
61            })
62            .collect::<Result<Vec<_>>>()?;
63        Ok(Self { k0, radii, degrees })
64    }
65}
66
67impl SavedState for SphereResidual {
68    fn save_state(&self) -> Result<Vec<u8>> {
69        let mut writer = Writer::new(Self::state_size(self.degrees.len(), self.radii.len())?);
70        self.write_state(&mut writer);
71        Ok(writer.finish())
72    }
73
74    fn from_state(bytes: &[u8]) -> Result<Self> {
75        let mut reader = Reader::new(bytes);
76        let residual = Self::read_state(&mut reader)?;
77        reader.finish()?;
78        Ok(residual)
79    }
80}
81
82#[cfg(test)]
83mod tests {
84    #![allow(clippy::float_cmp)] // Serialization preserves numerical bits exactly.
85
86    use super::*;
87    use crate::{Complex, coeffs::Material, tmatrix::sphere};
88    use nalgebra::DMatrix;
89
90    #[test]
91    fn saved_sphere_reuses_identical_derivatives() -> Result<()> {
92        let materials = [
93            Material {
94                epsilon: Complex::new(2.2, 0.1),
95                kappa: Complex::new(0.06, 0.0),
96                ..Material::default()
97            },
98            Material::default(),
99        ];
100        let (_, original) = sphere(2, 1.3, &[0.7], &materials)?;
101        let bytes = original.save_state()?;
102        assert_eq!(bytes.len(), SphereResidual::state_size(2, 1)?);
103        let restored = SphereResidual::from_state(&bytes)?;
104        let direction = [Material {
105            epsilon: Complex::new(0.2, 0.1),
106            mu: Complex::new(-0.1, 0.0),
107            kappa: Complex::new(0.03, 0.0),
108        }; 2];
109        let tangent = original.pushforward(0.2, &[-0.1], &direction)?;
110        let (rows, cols) = original.shape();
111        let cotangent = DMatrix::from_element(rows, cols, Complex::new(0.3, -0.2));
112        let expected = original.pullback(&cotangent)?;
113        for _ in 0..2 {
114            assert_eq!(restored.pushforward(0.2, &[-0.1], &direction)?, tangent);
115            let actual = restored.pullback(&cotangent)?;
116            assert_eq!(actual.k0, expected.k0);
117            assert_eq!(actual.layers.radii, expected.layers.radii);
118            assert_eq!(actual.layers.epsilon, expected.layers.epsilon);
119            assert_eq!(actual.layers.mu, expected.layers.mu);
120            assert_eq!(actual.layers.kappa, expected.layers.kappa);
121        }
122        assert_eq!(restored.save_state()?, bytes);
123        Ok(())
124    }
125
126    #[test]
127    fn saved_sphere_checks_lengths_before_allocation() -> Result<()> {
128        let (_, residual) = sphere(1, 1.0, &[0.5], &[Material::default(); 2])?;
129        let bytes = residual.save_state()?;
130        assert!(SphereResidual::from_state(&bytes[..bytes.len() - 1]).is_err());
131        let mut extended = bytes.clone();
132        extended.push(0);
133        assert!(SphereResidual::from_state(&extended).is_err());
134        for offset in [0, 8, 32] {
135            let mut corrupt = bytes.clone();
136            corrupt[offset..offset + 8].copy_from_slice(&u64::MAX.to_le_bytes());
137            assert!(SphereResidual::from_state(&corrupt).is_err());
138        }
139        assert!(SphereResidual::state_size(1, usize::MAX).is_err());
140        Ok(())
141    }
142}
treams_core/tmatrix/sphere/saved.rs

treams_core/tmatrix/sphere/
saved.rs
