conspire/math/tensor/quantity/sparse_vec_2d/
mod.rs1use 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
8pub 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}