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

1use crate::math::{
2    Hessian, SquareMatrix, Tensor, TensorRank0, TensorRank2, TensorRank2Vec,
3    tensor::vec::TensorVector,
4};
5use std::ops::Mul;
6
7use crate::math::assert::FiniteDifference;
8
9/// A vector of vectors of rank-2 tensors.
10pub type TensorRank2Vec2D<const D: usize, const I: usize, const J: usize> =
11    TensorVector<TensorRank2Vec<D, I, J>>;
12
13impl<const D: usize, const I: usize, const J: usize> TensorRank2Vec2D<D, I, J> {
14    pub fn zero(len: usize) -> Self {
15        (0..len).map(|_| TensorRank2Vec::zero(len)).collect()
16    }
17}
18
19impl<const D: usize, const I: usize, const J: usize> From<TensorRank2Vec2D<D, I, J>>
20    for Vec<TensorRank0>
21{
22    fn from(tensor_rank_2_vec_2d: TensorRank2Vec2D<D, I, J>) -> Self {
23        tensor_rank_2_vec_2d
24            .into_iter()
25            .flat_map(|tensor_rank_2_vec_1d| {
26                tensor_rank_2_vec_1d.into_iter().flat_map(|tensor_rank_2| {
27                    tensor_rank_2
28                        .into_iter()
29                        .flat_map(|tensor_rank_1| tensor_rank_1.into_iter())
30                })
31            })
32            .collect()
33    }
34}
35
36impl<const D: usize, const I: usize, const J: usize> Hessian for TensorRank2Vec2D<D, I, J> {
37    fn entry(&self, row: usize, column: usize) -> TensorRank0 {
38        self[row / D][column / D][row % D][column % D]
39    }
40    fn fill_into(self, square_matrix: &mut SquareMatrix) {
41        self.into_iter().enumerate().for_each(|(a, entry_a)| {
42            entry_a.into_iter().enumerate().for_each(|(b, entry_ab)| {
43                entry_ab
44                    .into_iter()
45                    .enumerate()
46                    .for_each(|(i, entry_ab_i)| {
47                        entry_ab_i
48                            .into_iter()
49                            .enumerate()
50                            .for_each(|(j, entry_ab_ij)| {
51                                square_matrix[D * a + i][D * b + j] = entry_ab_ij
52                            })
53                    })
54            })
55        });
56    }
57    fn retain_from(self, retained: &[bool]) -> SquareMatrix {
58        SquareMatrix::from(self)
59            .into_iter()
60            .zip(retained.iter())
61            .filter(|(_, retained)| **retained)
62            .map(|(self_i, _)| {
63                self_i
64                    .into_iter()
65                    .zip(retained.iter())
66                    .filter(|(_, retained)| **retained)
67                    .map(|(self_ij, _)| self_ij)
68                    .collect()
69            })
70            .collect()
71    }
72}
73
74impl<const D: usize, const I: usize, const J: usize, const K: usize> Mul<TensorRank2<D, J, K>>
75    for TensorRank2Vec2D<D, I, J>
76{
77    type Output = TensorRank2Vec2D<D, I, K>;
78    fn mul(self, tensor_rank_2: TensorRank2<D, J, K>) -> Self::Output {
79        self.iter()
80            .map(|self_entry| {
81                self_entry
82                    .iter()
83                    .map(|self_tensor_rank_2| self_tensor_rank_2 * &tensor_rank_2)
84                    .collect()
85            })
86            .collect()
87    }
88}
89
90impl<const D: usize, const I: usize, const J: usize, const K: usize> Mul<&TensorRank2<D, J, K>>
91    for TensorRank2Vec2D<D, I, J>
92{
93    type Output = TensorRank2Vec2D<D, I, K>;
94    fn mul(self, tensor_rank_2: &TensorRank2<D, J, K>) -> Self::Output {
95        self.iter()
96            .map(|self_entry| {
97                self_entry
98                    .iter()
99                    .map(|self_tensor_rank_2| self_tensor_rank_2 * tensor_rank_2)
100                    .collect()
101            })
102            .collect()
103    }
104}
105
106impl<const D: usize, const I: usize, const J: usize> FiniteDifference
107    for TensorRank2Vec2D<D, I, J>
108{
109    fn error_fd(&self, comparator: &Self, epsilon: TensorRank0) -> Option<(bool, usize)> {
110        let error_count = self
111            .iter()
112            .zip(comparator.iter())
113            .map(|(self_a, comparator_a)| {
114                self_a
115                    .iter()
116                    .zip(comparator_a.iter())
117                    .map(|(self_ab, comparator_ab)| {
118                        self_ab
119                            .iter()
120                            .zip(comparator_ab.iter())
121                            .map(|(self_ab_i, comparator_ab_i)| {
122                                self_ab_i
123                                    .iter()
124                                    .zip(comparator_ab_i.iter())
125                                    .filter(|&(&self_ab_ij, &comparator_ab_ij)| {
126                                        (self_ab_ij / comparator_ab_ij - 1.0).abs() >= epsilon
127                                            && (self_ab_ij.abs() >= epsilon
128                                                || comparator_ab_ij.abs() >= epsilon)
129                                    })
130                                    .count()
131                            })
132                            .sum::<usize>()
133                    })
134                    .sum::<usize>()
135            })
136            .sum();
137        if error_count > 0 {
138            let auxiliary = self
139                .iter()
140                .zip(comparator.iter())
141                .map(|(self_a, comparator_a)| {
142                    self_a
143                        .iter()
144                        .zip(comparator_a.iter())
145                        .map(|(self_ab, comparator_ab)| {
146                            self_ab
147                                .iter()
148                                .zip(comparator_ab.iter())
149                                .map(|(self_ab_i, comparator_ab_i)| {
150                                    self_ab_i
151                                        .iter()
152                                        .zip(comparator_ab_i.iter())
153                                        .filter(|&(&self_ab_ij, &comparator_ab_ij)| {
154                                            (self_ab_ij / comparator_ab_ij - 1.0).abs() >= epsilon
155                                                && (self_ab_ij - comparator_ab_ij).abs() >= epsilon
156                                                && (self_ab_ij.abs() >= epsilon
157                                                    || comparator_ab_ij.abs() >= epsilon)
158                                        })
159                                        .count()
160                                })
161                                .sum::<usize>()
162                        })
163                        .sum::<usize>()
164                })
165                .sum::<usize>()
166                > 0;
167            Some((auxiliary, error_count))
168        } else {
169            None
170        }
171    }
172}
conspire/math/tensor/rank_2/vec_2d/mod.rs

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