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

1#[cfg(test)]
2mod test;
3
4use std::f64::consts::TAU;
5
6use super::{
7    super::{
8        super::assert_eq_within_tols,
9        Rank2, Tensor, TensorArray, TensorError,
10        rank_0::{TensorRank0, list::TensorRank0List},
11        rank_1::{CrossProduct, TensorRank1},
12        rank_4::TensorRank4,
13    },
14    TensorRank2,
15};
16use crate::ABS_TOL;
17
18impl<const I: usize> TensorRank2<3, I, I> {
19    /// Returns the matrix logarithm of the 3x3 symmetric tensor.
20    pub fn logm(&self) -> Result<Self, TensorError> {
21        if self.is_diagonal() {
22            let mut logm = TensorRank2::zero();
23            logm.iter_mut()
24                .enumerate()
25                .zip(self.iter())
26                .for_each(|((i, logm_i), self_i)| logm_i[i] = self_i[i].ln());
27            Ok(logm)
28        } else {
29            let tensor = self - &TensorRank2::identity();
30            let norm = tensor.norm();
31            if norm < 1e-2 {
32                let num_terms = if norm < 1e-4 {
33                    2
34                } else if norm < 1e-3 {
35                    3
36                } else {
37                    5
38                };
39                let mut logm = tensor.clone();
40                let mut power = tensor.clone();
41                (2..=num_terms).for_each(|k| {
42                    power *= &tensor;
43                    logm += &power * (if k % 2 == 0 { -1.0 } else { 1.0 } / k as f64);
44                });
45                Ok(logm)
46            } else if self.is_symmetric() {
47                let mut eigenvalues = solve_cubic_symmetric(self.invariants())?;
48                if eigenvalues.iter().any(|eigenvalue| eigenvalue <= &0.0) {
49                    panic!("Symmetric matrix has a non-positive eigenvalue")
50                }
51                let eigenvectors = find_orthonormal_eigenvectors(&eigenvalues, self);
52                eigenvalues
53                    .iter_mut()
54                    .for_each(|eigenvalue| *eigenvalue = eigenvalue.ln());
55                Ok(reconstruct_symmetric(eigenvalues, eigenvectors))
56            } else {
57                panic!("Matrix logarithm only implemented for symmetric cases")
58            }
59        }
60    }
61    /// Returns the derivative of the matrix logarithm of the 3x3 symmetric tensor.
62    pub fn dlogm(&self) -> Result<TensorRank4<3, I, I, I, I>, TensorError> {
63        if self.is_diagonal() {
64            let mut dlogm = TensorRank4::zero();
65            dlogm.iter_mut().enumerate().for_each(|(i, dlogm_i)| {
66                dlogm_i.iter_mut().enumerate().for_each(|(j, dlogm_ij)| {
67                    dlogm_ij.iter_mut().enumerate().for_each(|(k, dlogm_ijk)| {
68                        dlogm_ijk
69                            .iter_mut()
70                            .enumerate()
71                            .filter(|(l, _)| i == k && &j == l)
72                            .for_each(|(_, dlogm_ijkl)| {
73                                *dlogm_ijkl = if assert_eq_within_tols(&self[i][i], &self[j][j])
74                                    .is_ok()
75                                {
76                                    1.0 / self[j][j]
77                                } else {
78                                    (self[i][i].ln() - self[j][j].ln()) / (self[i][i] - self[j][j])
79                                }
80                            })
81                    })
82                })
83            });
84            Ok(dlogm)
85        } else if self.is_symmetric() {
86            let eigenvalues = solve_cubic_symmetric(self.invariants())?;
87            if eigenvalues.iter().any(|eigenvalue| eigenvalue <= &0.0) {
88                panic!("Symmetric matrix has a non-positive eigenvalue")
89            }
90            let divided_difference: Self = eigenvalues
91                .iter()
92                .map(|eigenvalue_i| {
93                    eigenvalues
94                        .iter()
95                        .map(|eigenvalue_j| {
96                            if assert_eq_within_tols(eigenvalue_i, eigenvalue_j).is_ok() {
97                                1.0 / eigenvalue_j
98                            } else {
99                                (eigenvalue_i.ln() - eigenvalue_j.ln())
100                                    / (eigenvalue_i - eigenvalue_j)
101                            }
102                        })
103                        .collect()
104                })
105                .collect();
106            let eigenvectors = find_orthonormal_eigenvectors(&eigenvalues, self).transpose();
107            Ok(eigenvectors.iter().map(|eigenvector_i|
108                eigenvectors.iter().map(|eigenvector_j|
109                    eigenvectors.iter().map(|eigenvector_k|
110                        eigenvectors.iter().map(|eigenvector_l|
111                            eigenvector_i.iter().zip(eigenvector_k.iter().zip(divided_difference.iter())).map(|(eigenvector_ip, (eigenvector_kp, divided_difference_p))|
112                                eigenvector_j.iter().zip(eigenvector_l.iter().zip(divided_difference_p.iter())).map(|(eigenvector_jq, (eigenvector_lq, divided_difference_pq))|
113                                    eigenvector_ip * eigenvector_kp * divided_difference_pq * eigenvector_jq * eigenvector_lq
114                                ).sum::<TensorRank0>()
115                            ).sum()
116                        ).collect()
117                    ).collect()
118                ).collect()
119            ).collect())
120        } else {
121            panic!("Matrix logarithm only implemented for symmetric cases")
122        }
123    }
124    /// Returns the invariants of the 3x3 symmetric tensor.
125    pub fn invariants(&self) -> TensorRank0List<3> {
126        let trace = self.trace();
127        TensorRank0List::from([
128            trace,
129            0.5 * (trace.powi(2) - self.squared_trace()),
130            self.determinant(),
131        ])
132    }
133}
134
135fn solve_cubic_symmetric(
136    coefficients: TensorRank0List<3>,
137) -> Result<TensorRank0List<3>, TensorError> {
138    let c2 = coefficients[0];
139    let c1 = coefficients[1];
140    let c0 = coefficients[2];
141    let p = c1 - c2 * c2 / 3.0;
142    let q = -(2.0 * c2.powi(3) - 9.0 * c2 * c1 + 27.0 * c0) / 27.0;
143    if p.abs() < ABS_TOL {
144        let t = (-q).cbrt();
145        let lambda = t + c2 / 3.0;
146        return Ok(TensorRank0List::from([lambda; _]));
147    }
148    let discriminant = -4.0 * p * p * p - 27.0 * q * q;
149    let scale = (4.0 * p * p * p).abs().max(27.0 * q * q);
150    if discriminant.abs() <= 1e-13 * scale {
151        let r = (q / 2.0).cbrt();
152        let lambda_double = r + c2 / 3.0;
153        let lambda_simple = -2.0 * r + c2 / 3.0;
154        let lambdas = if lambda_double >= lambda_simple {
155            [lambda_double, lambda_double, lambda_simple]
156        } else {
157            [lambda_simple, lambda_double, lambda_double]
158        };
159        Ok(TensorRank0List::from(lambdas))
160    } else if discriminant > 0.0 {
161        let sqrt_term = (-p / 3.0).sqrt();
162        let cos_arg = 3.0 * q / (2.0 * p * (-p / 3.0).sqrt());
163        let cos_arg = cos_arg.clamp(-1.0, 1.0);
164        let theta = cos_arg.acos();
165        let mut lambdas = [
166            2.0 * sqrt_term * (theta / 3.0).cos() + c2 / 3.0,
167            2.0 * sqrt_term * ((theta + TAU) / 3.0).cos() + c2 / 3.0,
168            2.0 * sqrt_term * ((theta + 2.0 * TAU) / 3.0).cos() + c2 / 3.0,
169        ];
170        lambdas.iter_mut().for_each(|lambda| {
171            for _ in 0..2 {
172                let x = *lambda;
173                let f = x * x * x - c2 * x * x + c1 * x - c0;
174                let f_prime = 3.0 * x * x - 2.0 * c2 * x + c1;
175                if f_prime.abs() < ABS_TOL {
176                    break;
177                }
178                *lambda -= f / f_prime;
179            }
180        });
181        lambdas.sort_by(|a, b| b.partial_cmp(a).unwrap());
182        Ok(TensorRank0List::from(lambdas))
183    } else {
184        Err(TensorError::SymmetricMatrixComplexEigenvalues)
185    }
186}
187
188fn find_orthonormal_eigenvectors<const I: usize>(
189    eigenvalues: &TensorRank0List<3>,
190    tensor: &TensorRank2<3, I, I>,
191) -> TensorRank2<3, I, I> {
192    if assert_eq_within_tols(&eigenvalues[0], &eigenvalues[1]).is_ok() {
193        let mut eigenvectors = TensorRank2::zero();
194        eigenvectors[2] = eigenvector_symmetric(eigenvalues[2], tensor);
195        eigenvectors[0] = orthogonal_unit_vector(&eigenvectors[2]);
196        eigenvectors[1] = eigenvectors[2].cross(&eigenvectors[0]);
197        eigenvectors
198    } else if assert_eq_within_tols(&eigenvalues[1], &eigenvalues[2]).is_ok() {
199        let mut eigenvectors = TensorRank2::zero();
200        eigenvectors[0] = eigenvector_symmetric(eigenvalues[0], tensor);
201        eigenvectors[1] = orthogonal_unit_vector(&eigenvectors[0]);
202        eigenvectors[2] = eigenvectors[0].cross(&eigenvectors[1]);
203        eigenvectors
204    } else {
205        let mut eigenvectors = eigenvalues
206            .iter()
207            .map(|&eigenvalue| eigenvector_symmetric(eigenvalue, tensor))
208            .collect::<TensorRank2<3, I, I>>();
209        eigenvectors[0].normalize();
210        let proj1 = &eigenvectors[1] * &eigenvectors[0];
211        for i in 0..3 {
212            eigenvectors[1][i] -= proj1 * eigenvectors[0][i];
213        }
214        eigenvectors[1].normalize();
215        eigenvectors[2] = eigenvectors[0].cross(&eigenvectors[1]);
216        eigenvectors
217    }
218}
219
220fn orthogonal_unit_vector<const I: usize>(vector: &TensorRank1<3, I>) -> TensorRank1<3, I> {
221    let axis = vector
222        .iter()
223        .enumerate()
224        .min_by(|(_, a), (_, b)| a.abs().partial_cmp(&b.abs()).unwrap())
225        .map(|(i, _)| i)
226        .unwrap();
227    let mut other = TensorRank1::<3, I>::zero();
228    other[axis] = 1.0;
229    vector.cross(&other).normalized()
230}
231
232fn eigenvector_symmetric<const I: usize>(
233    eigenvalue: TensorRank0,
234    tensor: &TensorRank2<3, I, I>,
235) -> TensorRank1<3, I> {
236    let m = tensor - TensorRank2::identity() * eigenvalue;
237    [m[1].cross(&m[2]), m[0].cross(&m[2]), m[0].cross(&m[1])]
238        .into_iter()
239        .max_by(|a, b| a.norm().partial_cmp(&b.norm()).unwrap())
240        .unwrap()
241        .normalized()
242}
243
244fn reconstruct_symmetric<const I: usize>(
245    eigenvalues: TensorRank0List<3>,
246    eigenvectors: TensorRank2<3, I, I>,
247) -> TensorRank2<3, I, I> {
248    let mut tensor = TensorRank2::zero();
249    eigenvalues
250        .iter()
251        .zip(eigenvectors.iter())
252        .for_each(|(eigenvalue, eigenvector)| {
253            tensor
254                .iter_mut()
255                .zip(eigenvector.iter())
256                .for_each(|(tensor_i, eigenvector_i)| {
257                    tensor_i.iter_mut().zip(eigenvector.iter()).for_each(
258                        |(tensor_ij, eigenvector_j)| {
259                            *tensor_ij += eigenvalue * eigenvector_i * eigenvector_j
260                        },
261                    )
262                })
263        });
264    tensor
265}
conspire/math/tensor/rank_2/logarithm/mod.rs

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