Skip to main content

conspire/math/tensor/
mod.rs

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