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

1#[cfg(test)]
2mod test;
3
4use super::TensorRank2;
5use crate::math::{
6    Hessian, HessianAccumulate, Rank2, Scalar, SquareMatrix, Tensor, Vector,
7    tensor::vec::TensorVector,
8};
9use crate::units::{Dimensionless, UnitMul};
10use std::ops::Mul;
11
12use super::sparse_vec::TensorRank2SparseVec;
13
14use crate::math::{TensorArray, TensorRank0, assert::FiniteDifference};
15
16/// A vector of sparse vectors of rank-2 tensors, storing only inserted entries.
17pub type TensorRank2SparseVec2D<const D: usize, I, J, U = Dimensionless> =
18    TensorVector<TensorRank2SparseVec<D, I, J, U>>;
19
20impl<const D: usize, I, J, U> TensorRank2SparseVec2D<D, I, J, U> {
21    pub fn zero(len: usize) -> Self {
22        (0..len).map(|_| TensorRank2SparseVec::default()).collect()
23    }
24}
25
26impl<const D: usize, I, U> HessianAccumulate<D, I, U> for TensorRank2SparseVec2D<D, I, I, U> {
27    fn accumulate(&mut self, a: usize, b: usize, block: TensorRank2<D, I, I, U>) {
28        if a == b {
29            self[a][b] += block;
30        } else {
31            self[b][a] += block.transpose();
32            self[a][b] += block;
33        }
34    }
35}
36
37impl<const D: usize, I, J, K, U, V> Mul<TensorRank2SparseVec2D<D, J, K, V>>
38    for TensorRank2<D, I, J, U>
39where
40    U: UnitMul<V>,
41{
42    type Output = TensorRank2SparseVec2D<D, I, K, <U as UnitMul<V>>::Output>;
43    fn mul(self, tensor_rank_2_sparse_vec_2d: TensorRank2SparseVec2D<D, J, K, V>) -> Self::Output {
44        tensor_rank_2_sparse_vec_2d
45            .into_iter()
46            .map(|row| {
47                TensorRank2SparseVec(
48                    row.0
49                        .into_iter()
50                        .map(|(column, block)| (column, &self * block))
51                        .collect(),
52                )
53            })
54            .collect()
55    }
56}
57
58impl<const D: usize, I, J, K, U, V> Mul<TensorRank2<D, J, K, V>>
59    for TensorRank2SparseVec2D<D, I, J, U>
60where
61    U: UnitMul<V>,
62{
63    type Output = TensorRank2SparseVec2D<D, I, K, <U as UnitMul<V>>::Output>;
64    fn mul(self, tensor_rank_2: TensorRank2<D, J, K, V>) -> Self::Output {
65        self.into_iter()
66            .map(|row| {
67                TensorRank2SparseVec(
68                    row.0
69                        .into_iter()
70                        .map(|(column, block)| (column, block * &tensor_rank_2))
71                        .collect(),
72                )
73            })
74            .collect()
75    }
76}
77
78impl<const D: usize, I, J, U> Hessian for TensorRank2SparseVec2D<D, I, J, U> {
79    fn quadratic_form(&self, vector: &Vector) -> Scalar {
80        self.iter()
81            .enumerate()
82            .map(|(a, row)| {
83                row.entries()
84                    .map(|(b, block)| {
85                        block
86                            .iter()
87                            .enumerate()
88                            .map(|(i, block_i)| {
89                                block_i
90                                    .iter()
91                                    .enumerate()
92                                    .map(|(j, block_ij)| {
93                                        block_ij.value() * vector[D * a + i] * vector[D * b + j]
94                                    })
95                                    .sum::<Scalar>()
96                            })
97                            .sum::<Scalar>()
98                    })
99                    .sum::<Scalar>()
100            })
101            .sum()
102    }
103    fn entry(&self, row: usize, column: usize) -> Scalar {
104        match self[row / D]
105            .0
106            .binary_search_by_key(&(column / D), |&(b, _)| b)
107        {
108            Ok(k) => self[row / D].0[k].1[row % D][column % D].value(),
109            Err(_) => 0.0,
110        }
111    }
112    fn fill_into(self, square_matrix: &mut SquareMatrix) {
113        self.iter().enumerate().for_each(|(a, row)| {
114            row.entries().for_each(|(b, block)| {
115                block.iter().enumerate().for_each(|(i, block_i)| {
116                    block_i.iter().enumerate().for_each(|(j, block_ij)| {
117                        square_matrix[D * a + i][D * b + j] = block_ij.value()
118                    })
119                })
120            })
121        });
122    }
123    fn retain_from(self, retained: &[bool]) -> SquareMatrix {
124        let mut remap = vec![0; retained.len()];
125        let mut count = 0;
126        retained.iter().enumerate().for_each(|(p, &keep)| {
127            if keep {
128                remap[p] = count;
129                count += 1;
130            }
131        });
132        let mut square_matrix = SquareMatrix::zero(count);
133        self.iter().enumerate().for_each(|(a, row)| {
134            row.entries().for_each(|(b, block)| {
135                block.iter().enumerate().for_each(|(i, block_i)| {
136                    block_i.iter().enumerate().for_each(|(j, block_ij)| {
137                        if retained[D * a + i] && retained[D * b + j] {
138                            square_matrix[remap[D * a + i]][remap[D * b + j]] = block_ij.value()
139                        }
140                    })
141                })
142            })
143        });
144        square_matrix
145    }
146}
147
148impl<const D: usize, I, J, U> FiniteDifference for TensorRank2SparseVec2D<D, I, J, U> {
149    fn error_fd(&self, comparator: &Self, epsilon: TensorRank0) -> Option<(bool, usize)> {
150        let zero = TensorRank2::zero();
151        let block_errors =
152            |self_ab: &TensorRank2<D, I, J, U>, comparator_ab: &TensorRank2<D, I, J, U>| {
153                let mut errors = (0, 0);
154                self_ab.iter().zip(comparator_ab.iter()).for_each(
155                    |(self_ab_i, comparator_ab_i)| {
156                        self_ab_i.iter().zip(comparator_ab_i.iter()).for_each(
157                            |(&self_ab_ij, &comparator_ab_ij)| {
158                                if self_ab_ij.differs(comparator_ab_ij, epsilon) {
159                                    errors.0 += 1;
160                                    if self_ab_ij.differs_severely(comparator_ab_ij, epsilon) {
161                                        errors.1 += 1;
162                                    }
163                                }
164                            },
165                        )
166                    },
167                );
168                errors
169            };
170        let (error_count, severe_count) = self
171            .iter()
172            .zip(comparator.iter())
173            .map(|(self_a, comparator_a)| {
174                let mut errors = (0, 0);
175                let (mut p, mut q) = (0, 0);
176                while p < self_a.0.len() || q < comparator_a.0.len() {
177                    let b = self_a.0.get(p).map(|&(b, _)| b);
178                    let c = comparator_a.0.get(q).map(|&(c, _)| c);
179                    let block = match (b, c) {
180                        (Some(b), Some(c)) if b == c => {
181                            p += 1;
182                            q += 1;
183                            block_errors(&self_a.0[p - 1].1, &comparator_a.0[q - 1].1)
184                        }
185                        (Some(b), Some(c)) if b < c => {
186                            p += 1;
187                            block_errors(&self_a.0[p - 1].1, &zero)
188                        }
189                        (Some(_), None) => {
190                            p += 1;
191                            block_errors(&self_a.0[p - 1].1, &zero)
192                        }
193                        _ => {
194                            q += 1;
195                            block_errors(&zero, &comparator_a.0[q - 1].1)
196                        }
197                    };
198                    errors.0 += block.0;
199                    errors.1 += block.1;
200                }
201                errors
202            })
203            .fold((0, 0), |sum, errors| (sum.0 + errors.0, sum.1 + errors.1));
204        if error_count > 0 {
205            Some((severe_count > 0, error_count))
206        } else {
207            None
208        }
209    }
210}
conspire/math/tensor/rank_2/sparse_vec_2d/mod.rs

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