Skip to main content

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

1use super::{Quantity, sparse_vec::QuantitySparseVec};
2use crate::math::{
3    Hessian, Scalar, SquareMatrix, Tensor, TensorRank0, Vector, assert::FiniteDifference,
4    tensor::vec::TensorVector,
5};
6use crate::units::Dimensionless;
7
8/// A vector of sparse vectors of quantities, storing only inserted entries.
9pub type QuantitySparseVec2D<U = Dimensionless> = TensorVector<QuantitySparseVec<U>>;
10
11impl<U> QuantitySparseVec2D<U> {
12    pub fn zero(len: usize) -> Self {
13        (0..len).map(|_| QuantitySparseVec::default()).collect()
14    }
15}
16
17impl<U> Hessian for QuantitySparseVec2D<U> {
18    fn quadratic_form(&self, vector: &Vector) -> Scalar {
19        self.iter()
20            .enumerate()
21            .map(|(a, row)| {
22                row.entries()
23                    .map(|(b, entry)| entry.value() * vector[a] * vector[b])
24                    .sum::<Scalar>()
25            })
26            .sum()
27    }
28    fn entry(&self, row: usize, column: usize) -> Scalar {
29        match self[row].0.binary_search_by_key(&column, |&(b, _)| b) {
30            Ok(k) => self[row].0[k].1.value(),
31            Err(_) => 0.0,
32        }
33    }
34    fn fill_into(self, square_matrix: &mut SquareMatrix) {
35        self.iter().enumerate().for_each(|(a, row)| {
36            row.entries()
37                .for_each(|(b, entry)| square_matrix[a][b] = entry.value())
38        });
39    }
40    fn retain_from(self, retained: &[bool]) -> SquareMatrix {
41        let mut remap = vec![0; retained.len()];
42        let mut count = 0;
43        retained.iter().enumerate().for_each(|(p, &keep)| {
44            if keep {
45                remap[p] = count;
46                count += 1;
47            }
48        });
49        let mut square_matrix = SquareMatrix::zero(count);
50        self.iter().enumerate().for_each(|(a, row)| {
51            row.entries().for_each(|(b, entry)| {
52                if retained[a] && retained[b] {
53                    square_matrix[remap[a]][remap[b]] = entry.value()
54                }
55            })
56        });
57        square_matrix
58    }
59}
60
61impl<U> FiniteDifference for QuantitySparseVec2D<U> {
62    fn error_fd(&self, comparator: &Self, epsilon: TensorRank0) -> Option<(bool, usize)> {
63        let zero = Quantity::new(0.0);
64        let entry_errors = |self_ab: &Quantity<U>, comparator_ab: &Quantity<U>| {
65            if self_ab.differs(*comparator_ab, epsilon) {
66                (
67                    1,
68                    self_ab.differs_severely(*comparator_ab, epsilon) as usize,
69                )
70            } else {
71                (0, 0)
72            }
73        };
74        let (error_count, severe_count) = self
75            .iter()
76            .zip(comparator.iter())
77            .map(|(self_a, comparator_a)| {
78                let mut errors = (0, 0);
79                let (mut p, mut q) = (0, 0);
80                while p < self_a.0.len() || q < comparator_a.0.len() {
81                    let b = self_a.0.get(p).map(|&(b, _)| b);
82                    let c = comparator_a.0.get(q).map(|&(c, _)| c);
83                    let entry = match (b, c) {
84                        (Some(b), Some(c)) if b == c => {
85                            p += 1;
86                            q += 1;
87                            entry_errors(&self_a.0[p - 1].1, &comparator_a.0[q - 1].1)
88                        }
89                        (Some(b), Some(c)) if b < c => {
90                            p += 1;
91                            entry_errors(&self_a.0[p - 1].1, &zero)
92                        }
93                        (Some(_), None) => {
94                            p += 1;
95                            entry_errors(&self_a.0[p - 1].1, &zero)
96                        }
97                        _ => {
98                            q += 1;
99                            entry_errors(&zero, &comparator_a.0[q - 1].1)
100                        }
101                    };
102                    errors.0 += entry.0;
103                    errors.1 += entry.1;
104                }
105                errors
106            })
107            .fold((0, 0), |sum, errors| (sum.0 + errors.0, sum.1 + errors.1));
108        if error_count > 0 {
109            Some((severe_count > 0, error_count))
110        } else {
111            None
112        }
113    }
114}