conspire/math/tensor/rank_2/sparse_symmetric_vec_2d/
mod.rs

1#[cfg(test)]
2mod test;
3
4use super::TensorRank2;
5use crate::math::{
6    Hessian, HessianAccumulate, Rank2, Scalar, SquareMatrix, Tensor, TensorRank0, Vector,
7};
8use crate::units::Dimensionless;
9use std::{
10    fmt::{self, Debug, Display, Formatter},
11    iter::Sum,
12    ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Sub, SubAssign},
13};
14
15use super::sparse_vec::TensorRank2SparseVec;
16use super::sparse_vec_2d::TensorRank2SparseVec2D;
17
18use crate::math::{TensorArray, assert::FiniteDifference};
19
20/// A vector of sparse vectors of rank-2 tensors, storing only the symmetric half.
21///
22/// The underlying block matrix is known to be symmetric under index-pair
23/// transpose, meaning block(a, b) == block(b, a)ᵀ for every pair of block
24/// indices. Only the canonical (row <= column) half of the blocks is stored;
25/// entries on the other side are reconstructed by transposing on lookup
26/// instead of being duplicated in memory.
27pub struct TensorRank2SparseVec2DSymmetric<const D: usize, I, J, U = Dimensionless>(
28    TensorRank2SparseVec2D<D, I, J, U>,
29);
30
31impl<const D: usize, I, J, U> Clone for TensorRank2SparseVec2DSymmetric<D, I, J, U> {
32    fn clone(&self) -> Self {
33        Self(self.0.clone())
34    }
35}
36
37impl<const D: usize, I, J, U> Debug for TensorRank2SparseVec2DSymmetric<D, I, J, U> {
38    fn fmt(&self, f: &mut Formatter) -> fmt::Result {
39        Debug::fmt(&self.0, f)
40    }
41}
42
43impl<const D: usize, I, J, U> Default for TensorRank2SparseVec2DSymmetric<D, I, J, U> {
44    fn default() -> Self {
45        Self(Default::default())
46    }
47}
48
49impl<const D: usize, I, J, U> PartialEq for TensorRank2SparseVec2DSymmetric<D, I, J, U> {
50    fn eq(&self, other: &Self) -> bool {
51        self.0 == other.0
52    }
53}
54
55impl<const D: usize, I, J, U> TensorRank2SparseVec2DSymmetric<D, I, J, U> {
56    pub fn zero(len: usize) -> Self {
57        Self(TensorRank2SparseVec2D::zero(len))
58    }
59}
60
61impl<const D: usize, I, J, U> Display for TensorRank2SparseVec2DSymmetric<D, I, J, U> {
62    fn fmt(&self, f: &mut Formatter) -> fmt::Result {
63        write!(f, "Need to implement Display")
64    }
65}
66
67impl<const D: usize, I, J, U> Tensor for TensorRank2SparseVec2DSymmetric<D, I, J, U> {
68    type Item = TensorRank2SparseVec<D, I, J, U>;
69    type Unit = U;
70    fn iter(&self) -> impl Iterator<Item = &Self::Item> {
71        self.0.iter()
72    }
73    fn iter_mut(&mut self) -> impl Iterator<Item = &mut Self::Item> {
74        self.0.iter_mut()
75    }
76    fn len(&self) -> usize {
77        self.0.len()
78    }
79    fn size(&self) -> usize {
80        self.0.size()
81    }
82}
83
84impl<const D: usize, I, J, U> Add for TensorRank2SparseVec2DSymmetric<D, I, J, U> {
85    type Output = Self;
86    fn add(self, other: Self) -> Self {
87        Self(self.0 + other.0)
88    }
89}
90
91impl<const D: usize, I, J, U> Add<&Self> for TensorRank2SparseVec2DSymmetric<D, I, J, U> {
92    type Output = Self;
93    fn add(self, other: &Self) -> Self {
94        Self(self.0 + &other.0)
95    }
96}
97
98impl<const D: usize, I, J, U> AddAssign for TensorRank2SparseVec2DSymmetric<D, I, J, U> {
99    fn add_assign(&mut self, other: Self) {
100        self.0 += other.0;
101    }
102}
103
104impl<const D: usize, I, J, U> AddAssign<&Self> for TensorRank2SparseVec2DSymmetric<D, I, J, U> {
105    fn add_assign(&mut self, other: &Self) {
106        self.0 += &other.0;
107    }
108}
109
110impl<const D: usize, I, J, U> Sub for TensorRank2SparseVec2DSymmetric<D, I, J, U> {
111    type Output = Self;
112    fn sub(self, other: Self) -> Self {
113        Self(self.0 - other.0)
114    }
115}
116
117impl<const D: usize, I, J, U> Sub<&Self> for TensorRank2SparseVec2DSymmetric<D, I, J, U> {
118    type Output = Self;
119    fn sub(self, other: &Self) -> Self {
120        Self(self.0 - &other.0)
121    }
122}
123
124impl<const D: usize, I, J, U> SubAssign for TensorRank2SparseVec2DSymmetric<D, I, J, U> {
125    fn sub_assign(&mut self, other: Self) {
126        self.0 -= other.0;
127    }
128}
129
130impl<const D: usize, I, J, U> SubAssign<&Self> for TensorRank2SparseVec2DSymmetric<D, I, J, U> {
131    fn sub_assign(&mut self, other: &Self) {
132        self.0 -= &other.0;
133    }
134}
135
136impl<const D: usize, I, J, U> Mul<TensorRank0> for TensorRank2SparseVec2DSymmetric<D, I, J, U> {
137    type Output = Self;
138    fn mul(self, scalar: TensorRank0) -> Self {
139        Self(self.0 * scalar)
140    }
141}
142
143impl<const D: usize, I, J, U> MulAssign<TensorRank0>
144    for TensorRank2SparseVec2DSymmetric<D, I, J, U>
145{
146    fn mul_assign(&mut self, scalar: TensorRank0) {
147        self.0 *= scalar;
148    }
149}
150
151impl<const D: usize, I, J, U> MulAssign<&TensorRank0>
152    for TensorRank2SparseVec2DSymmetric<D, I, J, U>
153{
154    fn mul_assign(&mut self, scalar: &TensorRank0) {
155        self.0 *= scalar;
156    }
157}
158
159impl<const D: usize, I, J, U> Div<TensorRank0> for TensorRank2SparseVec2DSymmetric<D, I, J, U> {
160    type Output = Self;
161    fn div(self, scalar: TensorRank0) -> Self {
162        Self(self.0 / scalar)
163    }
164}
165
166impl<const D: usize, I, J, U> DivAssign<TensorRank0>
167    for TensorRank2SparseVec2DSymmetric<D, I, J, U>
168{
169    fn div_assign(&mut self, scalar: TensorRank0) {
170        self.0 /= scalar;
171    }
172}
173
174impl<const D: usize, I, J, U> DivAssign<&TensorRank0>
175    for TensorRank2SparseVec2DSymmetric<D, I, J, U>
176{
177    fn div_assign(&mut self, scalar: &TensorRank0) {
178        self.0 /= scalar;
179    }
180}
181
182impl<const D: usize, I, J, U> Sum for TensorRank2SparseVec2DSymmetric<D, I, J, U> {
183    fn sum<T>(iter: T) -> Self
184    where
185        T: Iterator<Item = Self>,
186    {
187        iter.fold(Self::default(), |sum, entry| sum + entry)
188    }
189}
190
191impl<const D: usize, I, U> HessianAccumulate<D, I, U>
192    for TensorRank2SparseVec2DSymmetric<D, I, I, U>
193{
194    fn accumulate(&mut self, a: usize, b: usize, block: TensorRank2<D, I, I, U>) {
195        if a <= b {
196            self.0[a][b] += block;
197        } else {
198            self.0[b][a] += block.transpose();
199        }
200    }
201}
202
203impl<const D: usize, I, J, U> Hessian for TensorRank2SparseVec2DSymmetric<D, I, J, U> {
204    fn entry(&self, row: usize, column: usize) -> Scalar {
205        let (a, b, i, j) = if row / D <= column / D {
206            (row / D, column / D, row % D, column % D)
207        } else {
208            (column / D, row / D, column % D, row % D)
209        };
210        match self.0[a].0.binary_search_by_key(&b, |&(c, _)| c) {
211            Ok(k) => self.0[a].0[k].1[i][j].value(),
212            Err(_) => 0.0,
213        }
214    }
215    fn quadratic_form(&self, vector: &Vector) -> Scalar {
216        //
217        // Only one triangle is stored, and the entry mirroring a stored one
218        // contributes the very same product, so an off-diagonal block counts
219        // twice. A diagonal block is stored whole and counts once.
220        //
221        self.0
222            .iter()
223            .enumerate()
224            .map(|(a, row)| {
225                row.entries()
226                    .map(|(b, block)| {
227                        block
228                            .iter()
229                            .enumerate()
230                            .map(|(i, block_i)| {
231                                block_i
232                                    .iter()
233                                    .enumerate()
234                                    .map(|(j, block_ij)| {
235                                        block_ij.value() * vector[D * a + i] * vector[D * b + j]
236                                    })
237                                    .sum::<Scalar>()
238                            })
239                            .sum::<Scalar>()
240                            * if a == b { 1.0 } else { 2.0 }
241                    })
242                    .sum::<Scalar>()
243            })
244            .sum()
245    }
246    fn fill_into(self, square_matrix: &mut SquareMatrix) {
247        self.0.iter().enumerate().for_each(|(a, row)| {
248            row.entries().for_each(|(b, block)| {
249                block.iter().enumerate().for_each(|(i, block_i)| {
250                    block_i.iter().enumerate().for_each(|(j, block_ij)| {
251                        square_matrix[D * a + i][D * b + j] = block_ij.value();
252                        if a != b {
253                            square_matrix[D * b + j][D * a + i] = block_ij.value();
254                        }
255                    })
256                })
257            })
258        });
259    }
260    fn retain_from(self, retained: &[bool]) -> SquareMatrix {
261        let mut remap = vec![0; retained.len()];
262        let mut count = 0;
263        retained.iter().enumerate().for_each(|(p, &keep)| {
264            if keep {
265                remap[p] = count;
266                count += 1;
267            }
268        });
269        let mut square_matrix = SquareMatrix::zero(count);
270        self.0.iter().enumerate().for_each(|(a, row)| {
271            row.entries().for_each(|(b, block)| {
272                block.iter().enumerate().for_each(|(i, block_i)| {
273                    block_i.iter().enumerate().for_each(|(j, block_ij)| {
274                        if retained[D * a + i] && retained[D * b + j] {
275                            square_matrix[remap[D * a + i]][remap[D * b + j]] = block_ij.value();
276                            if a != b {
277                                square_matrix[remap[D * b + j]][remap[D * a + i]] =
278                                    block_ij.value();
279                            }
280                        }
281                    })
282                })
283            })
284        });
285        square_matrix
286    }
287}
288
289impl<const D: usize, I, J, U> FiniteDifference for TensorRank2SparseVec2DSymmetric<D, I, J, U> {
290    fn error_fd(&self, comparator: &Self, epsilon: TensorRank0) -> Option<(bool, usize)> {
291        let zero = TensorRank2::zero();
292        let block_errors =
293            |self_ab: &TensorRank2<D, I, J, U>, comparator_ab: &TensorRank2<D, I, J, U>| {
294                let mut errors = (0, 0);
295                self_ab.iter().zip(comparator_ab.iter()).for_each(
296                    |(self_ab_i, comparator_ab_i)| {
297                        self_ab_i.iter().zip(comparator_ab_i.iter()).for_each(
298                            |(&self_ab_ij, &comparator_ab_ij)| {
299                                if self_ab_ij.differs(comparator_ab_ij, epsilon) {
300                                    errors.0 += 1;
301                                    if self_ab_ij.differs_severely(comparator_ab_ij, epsilon) {
302                                        errors.1 += 1;
303                                    }
304                                }
305                            },
306                        )
307                    },
308                );
309                errors
310            };
311        let (error_count, severe_count) = self
312            .0
313            .iter()
314            .zip(comparator.0.iter())
315            .map(|(self_a, comparator_a)| {
316                let mut errors = (0, 0);
317                let (mut p, mut q) = (0, 0);
318                while p < self_a.0.len() || q < comparator_a.0.len() {
319                    let b = self_a.0.get(p).map(|&(b, _)| b);
320                    let c = comparator_a.0.get(q).map(|&(c, _)| c);
321                    let block = match (b, c) {
322                        (Some(b), Some(c)) if b == c => {
323                            p += 1;
324                            q += 1;
325                            block_errors(&self_a.0[p - 1].1, &comparator_a.0[q - 1].1)
326                        }
327                        (Some(b), Some(c)) if b < c => {
328                            p += 1;
329                            block_errors(&self_a.0[p - 1].1, &zero)
330                        }
331                        (Some(_), None) => {
332                            p += 1;
333                            block_errors(&self_a.0[p - 1].1, &zero)
334                        }
335                        _ => {
336                            q += 1;
337                            block_errors(&zero, &comparator_a.0[q - 1].1)
338                        }
339                    };
340                    errors.0 += block.0;
341                    errors.1 += block.1;
342                }
343                errors
344            })
345            .fold((0, 0), |sum, errors| (sum.0 + errors.0, sum.1 + errors.1));
346        if error_count > 0 {
347            Some((severe_count > 0, error_count))
348        } else {
349            None
350        }
351    }
352}
conspire/math/tensor/rank_2/sparse_symmetric_vec_2d/mod.rs

conspire/math/tensor/rank_2/sparse_symmetric_vec_2d/
mod.rs
