conspire/math/tensor/
mod.rs

1#[cfg(test)]
2mod test;
3
4pub(super) mod configuration;
5pub(super) mod list;
6pub(super) mod norm;
7pub(super) mod quantity;
8pub(super) mod rank_0;
9pub(super) mod rank_1;
10pub(super) mod rank_2;
11pub(super) mod rank_3;
12pub(super) mod rank_4;
13pub(super) mod tuple;
14pub(super) mod vec;
15
16pub use configuration::{
17    Auxiliary, Configuration, Current, Factor, Flattened, Intermediate, Projection, Reference,
18};
19pub use norm::Norm;
20pub use quantity::{Is, Quantity};
21
22use super::{SquareMatrix, Vector};
23use crate::math::{Style, StyledError, styled_error};
24use crate::units::{Dimensionless, Time, UnitMul};
25use rank_0::{
26    TensorRank0,
27    list::{TensorRank0List, vec::TensorRank0ListVec},
28};
29use std::{
30    fmt::{Debug, Display},
31    iter::Sum,
32    ops::{Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Sub, SubAssign},
33};
34
35/// A scalar.
36pub type Scalar = TensorRank0;
37
38/// A vector of scalars.
39pub type Scalars = Vector;
40
41/// A list of scalars.
42pub type ScalarList<const N: usize> = TensorRank0List<N>;
43
44/// A vector of lists of scalars.
45pub type ScalarListVec<const N: usize> = TensorRank0ListVec<N>;
46
47/// Possible errors for tensors.
48#[derive(PartialEq)]
49pub enum TensorError {
50    NotPositiveDefinite,
51    SquareRootDidNotConverge,
52    SymmetricMatrixComplexEigenvalues,
53}
54
55impl StyledError for TensorError {
56    fn message(&self, style: &Style) -> String {
57        let h = style.headline;
58        match self {
59            Self::NotPositiveDefinite => format!("{h}Result is not positive definite."),
60            Self::SquareRootDidNotConverge => {
61                format!("{h}Matrix square root iteration did not converge.")
62            }
63            Self::SymmetricMatrixComplexEigenvalues => {
64                format!("{h}Symmetric matrix produced complex eigenvalues")
65            }
66        }
67    }
68}
69
70styled_error!(TensorError);
71
72/// A tensor that can be differentiated with respect to a variable of unit `T`.
73///
74/// The derivative is named by the tensor rather than passed in alongside it,
75/// being the same tensor with its unit divided by the variable's. A tuple
76/// computes each half's derivative on its own, so the pair of units its halves
77/// carry never has to be taken apart.
78///
79/// The variable of integration need not be a time — an arclength or a load
80/// parameter is just as ordinary — so it is named rather than assumed, with time
81/// as the default for the common case.
82pub trait Differentiable<T = Time>
83where
84    Self: Tensor,
85{
86    /// The derivative with respect to the variable of integration.
87    type Derivative: Tensor;
88}
89
90/// The derivative of `Y` with respect to a variable of unit `T`.
91///
92/// Spelling the projection out at every use would crowd out the signatures it
93/// appears in, since a tensor names a derivative for each variable it might be
94/// differentiated against.
95pub type Derivative<Y, T = Time> = <Y as Differentiable<T>>::Derivative;
96
97/// The unit a quantity of unit `U` carries once squared.
98///
99/// A tensor squared against itself — its norm squared, the trace of its square,
100/// the invariant built from the two — carries this rather than nothing. Spelling
101/// the projection out at every use would crowd out the signatures, as it would
102/// for a [`Derivative`].
103pub type Square<U> = <U as UnitMul<U>>::Output;
104
105/// The full contraction of two tensors whose units need not agree.
106///
107/// [`Tensor::full_contraction`] contracts a tensor with another of its own type
108/// and gives a number. Contracting tensors of different units gives a quantity
109/// whose unit is the product of theirs — a stress with a rate is a power
110/// density — which is what erased views used to stand in for.
111pub trait ContractWith<Rhs> {
112    /// The quantity the contraction gives.
113    type Output;
114    /// Returns the full contraction with the other tensor.
115    fn contract_with(&self, rhs: &Rhs) -> Self::Output;
116}
117
118/// Views a tensor with its configurations and unit discarded.
119///
120/// Contracting tensors of different units gives a quantity whose unit is the
121/// product, which is not always one this library names. Generic code that only
122/// wants the number contracts the erased views instead.
123pub trait Erase {
124    /// The tensor with its configurations and unit discarded.
125    type Erased: Tensor;
126    /// Views the tensor with its configurations and unit discarded.
127    fn erase(&self) -> &Self::Erased;
128}
129
130impl Erase for TensorRank0 {
131    type Erased = Self;
132    fn erase(&self) -> &Self {
133        self
134    }
135}
136
137/// Common methods for solutions.
138pub trait Solution
139where
140    Self: From<Vector> + Tensor,
141{
142    /// Decrements the solution from another vector.
143    fn decrement_from(&mut self, other: &Vector);
144    /// Decrements the solution chained with a vector from another vector.
145    fn decrement_from_chained(&mut self, other: &mut Vector, vector: &Vector);
146    /// Decrements the solution from another vector on retained entries.
147    fn decrement_from_retained(&mut self, _retained: &[bool], _other: &Vector) {
148        unimplemented!()
149    }
150}
151
152/// Common methods for Jacobians.
153pub trait Jacobian
154where
155    Self:
156        From<Vector> + Tensor + Sub<Vector, Output = Self> + for<'a> Sub<&'a Vector, Output = Self>,
157{
158    /// Fills the Jacobian into a vector.
159    fn fill_into(&self, vector: &mut Vector);
160    /// Fills the Jacobian chained with a vector into another vector.
161    fn fill_into_chained(self, other: Vector, vector: &mut Vector);
162    /// Return only the retained indices.
163    fn retain_from(self, _retained: &[bool]) -> Vector {
164        unimplemented!()
165    }
166    /// Zero out the specified indices.
167    fn zero_out(&mut self, _indices: &[usize]) {
168        unimplemented!()
169    }
170}
171
172/// Common methods for Hessians.
173pub trait Hessian
174where
175    Self: Tensor,
176{
177    /// The entry at the given (row, column) position.
178    fn entry(&self, row: usize, column: usize) -> Scalar;
179    /// Fills the Hessian into a square matrix.
180    fn fill_into(self, square_matrix: &mut SquareMatrix);
181    /// The quadratic form of the Hessian with a vector.
182    ///
183    /// ```math
184    /// \mathbf{v}\cdot\mathbf{H}\cdot\mathbf{v}
185    /// ```
186    fn quadratic_form(&self, _vector: &Vector) -> Scalar {
187        unimplemented!()
188    }
189    /// Return only the retained indices.
190    fn retain_from(self, _retained: &[bool]) -> SquareMatrix {
191        unimplemented!()
192    }
193}
194
195/// Common methods for blocks of a Hessian.
196pub trait HessianBlock {
197    /// The entry of the block at the given row and column within it.
198    fn entry(&self, row: usize, column: usize) -> TensorRank0;
199    /// The number of rows the block occupies.
200    fn height(&self) -> usize;
201    /// The number of columns the block occupies.
202    fn width(&self) -> usize;
203    /// Fills the block into a matrix at the given row and column offsets.
204    fn fill_into_block<M>(&self, matrix: &mut M, row: usize, column: usize)
205    where
206        M: IndexMut<usize, Output = Vector>;
207}
208
209/// A [`HessianBlock`] with its rows and columns swapped.
210pub struct Transposed<H>(pub H);
211
212impl<H: HessianBlock> HessianBlock for Transposed<H> {
213    fn entry(&self, row: usize, column: usize) -> TensorRank0 {
214        self.0.entry(column, row)
215    }
216    fn height(&self) -> usize {
217        self.0.width()
218    }
219    fn width(&self) -> usize {
220        self.0.height()
221    }
222    fn fill_into_block<M>(&self, matrix: &mut M, row: usize, column: usize)
223    where
224        M: IndexMut<usize, Output = Vector>,
225    {
226        (0..self.0.height()).for_each(|i| {
227            (0..self.0.width()).for_each(|j| matrix[row + j][column + i] = self.0.entry(i, j))
228        })
229    }
230}
231
232/// Accumulates rank-2 blocks into a sparse Hessian-like structure.
233///
234/// Symmetric-safe: the caller guarantees `block` at (a, b) equals the
235/// transpose of the (b, a) contribution, so implementors may store or
236/// mirror as they see fit.
237pub trait HessianAccumulate<const D: usize, I, U = Dimensionless> {
238    fn accumulate(&mut self, a: usize, b: usize, block: rank_2::TensorRank2<D, I, I, U>);
239}
240
241/// Common methods for rank-2 tensors.
242pub trait Rank2
243where
244    Self: Sized + Tensor,
245{
246    /// The type that is the transpose of the tensor.
247    type Transpose;
248    /// Returns the deviatoric component of the rank-2 tensor.
249    fn deviatoric(&self) -> Self;
250    /// Returns the deviatoric component and trace of the rank-2 tensor.
251    fn deviatoric_and_trace(&self) -> (Self, Quantity<Self::Unit>);
252    /// Checks whether the tensor is a diagonal tensor.
253    fn is_diagonal(&self) -> bool;
254    /// Checks whether the tensor is the identity tensor.
255    fn is_identity(&self) -> bool;
256    /// Checks whether the tensor is a symmetric tensor.
257    fn is_symmetric(&self) -> bool;
258    /// Returns the second invariant of the rank-2 tensor.
259    fn second_invariant(&self) -> Quantity<Square<Self::Unit>>
260    where
261        Self::Unit: UnitMul<Self::Unit>,
262    {
263        let trace = self.trace();
264        (trace * trace - self.squared_trace()) * 0.5
265    }
266    /// Returns the trace of the rank-2 tensor squared.
267    fn squared_trace(&self) -> Quantity<Square<Self::Unit>>
268    where
269        Self::Unit: UnitMul<Self::Unit>;
270    /// Returns the trace of the rank-2 tensor, which carries its unit.
271    fn trace(&self) -> Quantity<Self::Unit>;
272    /// Returns the transpose of the rank-2 tensor.
273    fn transpose(&self) -> Self::Transpose;
274}
275
276/// Common methods for tensors.
277#[expect(clippy::len_without_is_empty)]
278pub trait Tensor
279where
280    for<'a> Self: Sized
281        + Add<Self, Output = Self>
282        + Add<&'a Self, Output = Self>
283        + AddAssign
284        + AddAssign<&'a Self>
285        + Clone
286        + Debug
287        + Default
288        + Display
289        + Div<TensorRank0, Output = Self>
290        // + Div<&'a TensorRank0, Output = Self>
291        + DivAssign<TensorRank0>
292        + DivAssign<&'a TensorRank0>
293        + Mul<TensorRank0, Output = Self>
294        // + Mul<&'a TensorRank0, Output = Self>
295        + MulAssign<TensorRank0>
296        + MulAssign<&'a TensorRank0>
297        + Sub<Self, Output = Self>
298        + Sub<&'a Self, Output = Self>
299        + SubAssign
300        + SubAssign<&'a Self>
301        + Sum,
302    Self::Item: Tensor,
303{
304    /// The type of item encountered when iterating over the tensor.
305    type Item;
306    /// The physical unit the tensor carries.
307    type Unit;
308    /// Returns number of nonzero entries given absolute and relative tolerances.
309    fn error_count_zero(&self, tol_abs: Scalar, tol_rel: Scalar) -> Option<usize> {
310        let error_count = self
311            .iter()
312            .filter_map(|entry| entry.error_count_zero(tol_abs, tol_rel))
313            .sum();
314        if error_count > 0 {
315            Some(error_count)
316        } else {
317            None
318        }
319    }
320    /// Returns number of different entries given absolute and relative tolerances.
321    fn error_count(&self, other: &Self, tol_abs: Scalar, tol_rel: Scalar) -> Option<usize> {
322        let error_count = self
323            .iter()
324            .zip(other.iter())
325            .filter_map(|(self_entry, other_entry)| {
326                self_entry.error_count(other_entry, tol_abs, tol_rel)
327            })
328            .sum();
329        if error_count > 0 {
330            Some(error_count)
331        } else {
332            None
333        }
334    }
335    /// Returns the full contraction with another tensor.
336    fn full_contraction(&self, tensor: &Self) -> TensorRank0 {
337        self.iter()
338            .zip(tensor.iter())
339            .map(|(self_entry, tensor_entry)| self_entry.full_contraction(tensor_entry))
340            .fold(0.0, f64::algebraic_add)
341    }
342    /// Checks whether the tensor is the zero tensor.
343    fn is_zero(&self) -> bool {
344        self.iter().filter(|entry| !entry.is_zero()).count() == 0
345    }
346    /// Returns an iterator.
347    ///
348    /// The iterator yields all items from start to end. [Read more](https://doc.rust-lang.org/std/iter/)
349    fn iter(&self) -> impl Iterator<Item = &Self::Item>;
350    /// Returns an iterator that allows modifying each value.
351    ///
352    /// The iterator yields all items from start to end. [Read more](https://doc.rust-lang.org/std/iter/)
353    fn iter_mut(&mut self) -> impl Iterator<Item = &mut Self::Item>;
354    /// Returns the number of elements, also referred to as the ‘length’.
355    fn len(&self) -> usize;
356    /// Returns the tensor norm.
357    fn norm(&self) -> Quantity<Self::Unit> {
358        Quantity::new(self.full_contraction(self).sqrt())
359    }
360    /// Returns the infinity norm.
361    fn norm_inf(&self) -> Quantity<Self::Unit> {
362        Quantity::new(
363            self.iter()
364                .fold(0.0, |acc, entry| entry.norm_inf().value().max(acc)),
365        )
366    }
367    /// Returns the L1 (Manhattan) norm.
368    fn norm_l1(&self) -> Quantity<Self::Unit> {
369        Quantity::new(
370            self.iter()
371                .fold(0.0, |acc, entry| acc + entry.norm_l1().value()),
372        )
373    }
374    /// Returns the sum of p-th powers of absolute values (used internally by `norm_p`).
375    fn norm_p_sum(&self, p: TensorRank0) -> TensorRank0 {
376        self.iter()
377            .fold(0.0, |acc, entry| acc + entry.norm_p_sum(p))
378    }
379    /// Returns the Minkowski (Lp) norm.
380    fn norm_p(&self, p: TensorRank0) -> Quantity<Self::Unit> {
381        Quantity::new(self.norm_p_sum(p).powf(1.0 / p))
382    }
383    /// Returns the tensor norm squared, which carries the square of its unit.
384    fn norm_squared(&self) -> Quantity<Square<Self::Unit>>
385    where
386        Self::Unit: UnitMul<Self::Unit>,
387    {
388        Quantity::new(self.full_contraction(self))
389    }
390    /// Normalizes the tensor in place.
391    fn normalize(&mut self) {
392        *self /= self.norm().value()
393    }
394    /// Returns the total number of entries.
395    fn size(&self) -> usize;
396    /// Returns the positive difference of the two tensors.
397    fn sub_abs(&self, other: &Self) -> Self {
398        let mut difference = self.clone();
399        difference
400            .iter_mut()
401            .zip(self.iter().zip(other.iter()))
402            .for_each(|(entry, (self_entry, other_entry))| {
403                *entry = self_entry.sub_abs(other_entry)
404            });
405        difference
406    }
407    /// Returns the relative difference of the two tensors.
408    fn sub_rel(&self, other: &Self) -> Self {
409        let mut difference = self.clone();
410        difference
411            .iter_mut()
412            .zip(self.iter().zip(other.iter()))
413            .for_each(|(entry, (self_entry, other_entry))| {
414                *entry = self_entry.sub_rel(other_entry)
415            });
416        difference
417    }
418}
419
420/// Common methods for tensors derived from arrays.
421pub trait TensorArray {
422    /// The type of array corresponding to the tensor.
423    type Array;
424    /// The type of item encountered when iterating over the tensor.
425    type Item;
426    /// Returns the tensor as an array.
427    fn as_array(&self) -> Self::Array;
428    /// Returns the identity tensor.
429    fn identity() -> Self;
430    /// Returns the zero tensor.
431    fn zero() -> Self;
432}
433
434/// Common methods for tensors derived from Vec.
435pub trait TensorVec
436where
437    Self: FromIterator<Self::Item> + Index<usize, Output = Self::Item> + IndexMut<usize>,
438{
439    /// The type of element encountered when iterating over the tensor.
440    type Item;
441    /// Moves all the elements of other into self, leaving other empty.
442    fn append(&mut self, other: &mut Self);
443    /// Returns the total number of elements the vector can hold without reallocating.
444    fn capacity(&self) -> usize;
445    /// Returns `true` if the vector contains no elements.
446    fn is_empty(&self) -> bool;
447    /// Constructs a new, empty Vec, not allocating until elements are pushed onto it.
448    fn new() -> Self;
449    /// Appends an element to the back of the Vec.
450    fn push(&mut self, item: Self::Item);
451    /// Removes an element from the Vec and returns it, shifting elements to the left.
452    fn remove(&mut self, index: usize) -> Self::Item;
453    /// Reserves capacity for at least additional more elements to be inserted in the given Vec.
454    fn reserve(&mut self, additional: usize);
455    /// Retains only the elements specified by the predicate.
456    fn retain<F>(&mut self, f: F)
457    where
458        F: FnMut(&Self::Item) -> bool;
459    /// Removes an element from the Vec and returns it, replacing it with the last element.
460    fn swap_remove(&mut self, index: usize) -> Self::Item;
461    /// Constructs a new, empty vector with at least the specified capacity.
462    fn with_capacity(capacity: usize) -> Self;
463}
conspire/math/tensor/mod.rs

conspire/math/tensor/
mod.rs
