conspire/math/tensor/
mod.rs

1pub(super) mod list;
2pub(super) mod norm;
3pub(super) mod rank_0;
4pub(super) mod rank_1;
5pub(super) mod rank_2;
6pub(super) mod rank_3;
7pub(super) mod rank_4;
8pub(super) mod tuple;
9pub(super) mod vec;
10
11pub use norm::Norm;
12
13use super::{SquareMatrix, Vector};
14use crate::math::{Style, StyledError, styled_error};
15use rank_0::{
16    TensorRank0,
17    list::{TensorRank0List, vec::TensorRank0ListVec},
18};
19use std::{
20    fmt::{Debug, Display},
21    iter::Sum,
22    ops::{Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Sub, SubAssign},
23};
24
25/// A scalar.
26pub type Scalar = TensorRank0;
27
28/// A vector of scalars.
29pub type Scalars = Vector;
30
31/// A list of scalars.
32pub type ScalarList<const N: usize> = TensorRank0List<N>;
33
34/// A vector of lists of scalars.
35pub type ScalarListVec<const N: usize> = TensorRank0ListVec<N>;
36
37/// Possible errors for tensors.
38#[derive(PartialEq)]
39pub enum TensorError {
40    NotPositiveDefinite,
41    SymmetricMatrixComplexEigenvalues,
42}
43
44impl StyledError for TensorError {
45    fn message(&self, style: &Style) -> String {
46        let h = style.headline;
47        match self {
48            Self::NotPositiveDefinite => format!("{h}Result is not positive definite."),
49            Self::SymmetricMatrixComplexEigenvalues => {
50                format!("{h}Symmetric matrix produced complex eigenvalues")
51            }
52        }
53    }
54}
55
56styled_error!(TensorError);
57
58/// Common methods for solutions.
59pub trait Solution
60where
61    Self: From<Vector> + Tensor,
62{
63    /// Decrements the solution from another vector.
64    fn decrement_from(&mut self, other: &Vector);
65    /// Decrements the solution chained with a vector from another vector.
66    fn decrement_from_chained(&mut self, other: &mut Vector, vector: Vector);
67    /// Decrements the solution from another vector on retained entries.
68    fn decrement_from_retained(&mut self, _retained: &[bool], _other: &Vector) {
69        unimplemented!()
70    }
71}
72
73/// Common methods for Jacobians.
74pub trait Jacobian
75where
76    Self:
77        From<Vector> + Tensor + Sub<Vector, Output = Self> + for<'a> Sub<&'a Vector, Output = Self>,
78{
79    /// Fills the Jacobian into a vector.
80    fn fill_into(self, vector: &mut Vector);
81    /// Fills the Jacobian chained with a vector into another vector.
82    fn fill_into_chained(self, other: Vector, vector: &mut Vector);
83    /// Return only the retained indices.
84    fn retain_from(self, _retained: &[bool]) -> Vector {
85        unimplemented!()
86    }
87    /// Zero out the specified indices.
88    fn zero_out(&mut self, _indices: &[usize]) {
89        unimplemented!()
90    }
91}
92
93/// Common methods for Hessians.
94pub trait Hessian
95where
96    Self: Tensor,
97{
98    /// The entry at the given (row, column) position.
99    fn entry(&self, row: usize, column: usize) -> Scalar;
100    /// Fills the Hessian into a square matrix.
101    fn fill_into(self, square_matrix: &mut SquareMatrix);
102    /// Return only the retained indices.
103    fn retain_from(self, _retained: &[bool]) -> SquareMatrix {
104        unimplemented!()
105    }
106}
107
108/// Accumulates rank-2 blocks into a sparse Hessian-like structure.
109///
110/// Symmetric-safe: the caller guarantees `block` at (a, b) equals the
111/// transpose of the (b, a) contribution, so implementors may store or
112/// mirror as they see fit.
113pub trait HessianAccumulate<const D: usize, const I: usize> {
114    fn accumulate(&mut self, a: usize, b: usize, block: rank_2::TensorRank2<D, I, I>);
115}
116
117/// Accumulates a rank-2 block at an exact (a, b) position.
118///
119/// No symmetry assumed. Only implemented by accumulators that store every
120/// position explicitly.
121pub trait HessianAccumulateGeneral<const D: usize, const I: usize>:
122    HessianAccumulate<D, I>
123{
124    fn accumulate_general(&mut self, a: usize, b: usize, block: rank_2::TensorRank2<D, I, I>);
125}
126
127/// Common methods for rank-2 tensors.
128pub trait Rank2
129where
130    Self: Sized,
131{
132    /// The type that is the transpose of the tensor.
133    type Transpose;
134    /// Returns the deviatoric component of the rank-2 tensor.
135    fn deviatoric(&self) -> Self;
136    /// Returns the deviatoric component and trace of the rank-2 tensor.
137    fn deviatoric_and_trace(&self) -> (Self, TensorRank0);
138    /// Checks whether the tensor is a diagonal tensor.
139    fn is_diagonal(&self) -> bool;
140    /// Checks whether the tensor is the identity tensor.
141    fn is_identity(&self) -> bool;
142    /// Checks whether the tensor is a symmetric tensor.
143    fn is_symmetric(&self) -> bool;
144    /// Returns the second invariant of the rank-2 tensor.
145    fn second_invariant(&self) -> TensorRank0 {
146        0.5 * (self.trace().powi(2) - self.squared_trace())
147    }
148    /// Returns the trace of the rank-2 tensor squared.
149    fn squared_trace(&self) -> TensorRank0;
150    /// Returns the trace of the rank-2 tensor.
151    fn trace(&self) -> TensorRank0;
152    /// Returns the transpose of the rank-2 tensor.
153    fn transpose(&self) -> Self::Transpose;
154}
155
156/// Common methods for tensors.
157#[allow(clippy::len_without_is_empty)]
158pub trait Tensor
159where
160    for<'a> Self: Sized
161        + Add<Self, Output = Self>
162        + Add<&'a Self, Output = Self>
163        + AddAssign
164        + AddAssign<&'a Self>
165        + Clone
166        + Debug
167        + Default
168        + Display
169        + Div<TensorRank0, Output = Self>
170        // + Div<&'a TensorRank0, Output = Self>
171        + DivAssign<TensorRank0>
172        + DivAssign<&'a TensorRank0>
173        + Mul<TensorRank0, Output = Self>
174        // + Mul<&'a TensorRank0, Output = Self>
175        + MulAssign<TensorRank0>
176        + MulAssign<&'a TensorRank0>
177        + Sub<Self, Output = Self>
178        + Sub<&'a Self, Output = Self>
179        + SubAssign
180        + SubAssign<&'a Self>
181        + Sum,
182    Self::Item: Tensor,
183{
184    /// The type of item encountered when iterating over the tensor.
185    type Item;
186    /// Returns number of nonzero entries given absolute and relative tolerances, compared against zero.
187    fn error_count_zero(&self, tol_abs: Scalar, tol_rel: Scalar) -> Option<usize> {
188        let error_count = self
189            .iter()
190            .filter_map(|entry| entry.error_count_zero(tol_abs, tol_rel))
191            .sum();
192        if error_count > 0 {
193            Some(error_count)
194        } else {
195            None
196        }
197    }
198    /// Returns number of different entries given absolute and relative tolerances.
199    fn error_count(&self, other: &Self, tol_abs: Scalar, tol_rel: Scalar) -> Option<usize> {
200        let error_count = self
201            .iter()
202            .zip(other.iter())
203            .filter_map(|(self_entry, other_entry)| {
204                self_entry.error_count(other_entry, tol_abs, tol_rel)
205            })
206            .sum();
207        if error_count > 0 {
208            Some(error_count)
209        } else {
210            None
211        }
212    }
213    /// Returns the full contraction with another tensor.
214    fn full_contraction(&self, tensor: &Self) -> TensorRank0 {
215        self.iter()
216            .zip(tensor.iter())
217            .map(|(self_entry, tensor_entry)| self_entry.full_contraction(tensor_entry))
218            .sum()
219    }
220    /// Checks whether the tensor is the zero tensor.
221    fn is_zero(&self) -> bool {
222        self.iter().filter(|entry| !entry.is_zero()).count() == 0
223    }
224    /// Returns an iterator.
225    ///
226    /// The iterator yields all items from start to end. [Read more](https://doc.rust-lang.org/std/iter/)
227    fn iter(&self) -> impl Iterator<Item = &Self::Item>;
228    /// Returns an iterator that allows modifying each value.
229    ///
230    /// The iterator yields all items from start to end. [Read more](https://doc.rust-lang.org/std/iter/)
231    fn iter_mut(&mut self) -> impl Iterator<Item = &mut Self::Item>;
232    /// Returns the number of elements, also referred to as the ‘length’.
233    fn len(&self) -> usize;
234    /// Returns the tensor norm.
235    fn norm(&self) -> TensorRank0 {
236        self.norm_squared().sqrt()
237    }
238    /// Returns the infinity norm.
239    fn norm_inf(&self) -> TensorRank0 {
240        self.iter()
241            .fold(0.0, |acc, entry| entry.norm_inf().max(acc))
242    }
243    /// Returns the L1 (Manhattan) norm.
244    fn norm_l1(&self) -> TensorRank0 {
245        self.iter().fold(0.0, |acc, entry| acc + entry.norm_l1())
246    }
247    /// Returns the sum of p-th powers of absolute values (used internally by `norm_p`).
248    fn norm_p_sum(&self, p: TensorRank0) -> TensorRank0 {
249        self.iter()
250            .fold(0.0, |acc, entry| acc + entry.norm_p_sum(p))
251    }
252    /// Returns the Minkowski (Lp) norm.
253    fn norm_p(&self, p: TensorRank0) -> TensorRank0 {
254        self.norm_p_sum(p).powf(1.0 / p)
255    }
256    /// Returns the tensor norm squared.
257    fn norm_squared(&self) -> TensorRank0 {
258        self.full_contraction(self)
259    }
260    /// Normalizes the tensor.
261    fn normalize(&mut self) {
262        *self /= self.norm()
263    }
264    /// Returns the tensor normalized.
265    fn normalized(self) -> Self {
266        let norm = self.norm();
267        self / norm
268    }
269    /// Returns the total number of entries.
270    fn size(&self) -> usize;
271    /// Returns the positive difference of the two tensors.
272    fn sub_abs(&self, other: &Self) -> Self {
273        let mut difference = self.clone();
274        difference
275            .iter_mut()
276            .zip(self.iter().zip(other.iter()))
277            .for_each(|(entry, (self_entry, other_entry))| {
278                *entry = self_entry.sub_abs(other_entry)
279            });
280        difference
281    }
282    /// Returns the relative difference of the two tensors.
283    fn sub_rel(&self, other: &Self) -> Self {
284        let mut difference = self.clone();
285        difference
286            .iter_mut()
287            .zip(self.iter().zip(other.iter()))
288            .for_each(|(entry, (self_entry, other_entry))| {
289                *entry = self_entry.sub_rel(other_entry)
290            });
291        difference
292    }
293}
294
295/// Common methods for tensors derived from arrays.
296pub trait TensorArray {
297    /// The type of array corresponding to the tensor.
298    type Array;
299    /// The type of item encountered when iterating over the tensor.
300    type Item;
301    /// Returns the tensor as an array.
302    fn as_array(&self) -> Self::Array;
303    /// Returns the identity tensor.
304    fn identity() -> Self;
305    /// Returns the zero tensor.
306    fn zero() -> Self;
307}
308
309/// Common methods for tensors derived from Vec.
310pub trait TensorVec
311where
312    Self: FromIterator<Self::Item> + Index<usize, Output = Self::Item> + IndexMut<usize>,
313{
314    /// The type of element encountered when iterating over the tensor.
315    type Item;
316    /// Moves all the elements of other into self, leaving other empty.
317    fn append(&mut self, other: &mut Self);
318    /// Returns the total number of elements the vector can hold without reallocating.
319    fn capacity(&self) -> usize;
320    /// Returns `true` if the vector contains no elements.
321    fn is_empty(&self) -> bool;
322    /// Constructs a new, empty Vec, not allocating until elements are pushed onto it.
323    fn new() -> Self;
324    /// Appends an element to the back of the Vec.
325    fn push(&mut self, item: Self::Item);
326    /// Removes an element from the Vec and returns it, shifting elements to the left.
327    fn remove(&mut self, index: usize) -> Self::Item;
328    /// Reserves capacity for at least additional more elements to be inserted in the given Vec.
329    fn reserve(&mut self, additional: usize);
330    /// Retains only the elements specified by the predicate.
331    fn retain<F>(&mut self, f: F)
332    where
333        F: FnMut(&Self::Item) -> bool;
334    /// Removes an element from the Vec and returns it, replacing it with the last element.
335    fn swap_remove(&mut self, index: usize) -> Self::Item;
336    /// Constructs a new, empty vector with at least the specified capacity.
337    fn with_capacity(capacity: usize) -> Self;
338}
conspire/math/tensor/mod.rs

conspire/math/tensor/
mod.rs
