conspire/math/tensor/rank_4/
mod.rs

1#[cfg(test)]
2mod test;
3use super::rank_2::relabel as relabel_rank_2;
4use super::rank_3::relabel as relabel_rank_3;
5use crate::math::Quantity;
6use crate::math::{Current, Intermediate, Reference};
7use crate::units::{Dimensionless, UnitDiv, UnitMul};
8
9use crate::math::assert::FiniteDifference;
10
11use std::{
12    array::from_fn,
13    fmt::{self, Debug, Display, Formatter},
14    iter::Sum,
15    marker::PhantomData,
16    mem::transmute,
17    ops::{Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Sub, SubAssign},
18};
19
20use super::{
21    Differentiate, Erase, Hessian, HessianBlock, Rank2, SquareMatrix, Tensor, TensorArray, Vector,
22    rank_0::TensorRank0,
23    rank_1::TensorRank1,
24    rank_2::TensorRank2,
25    rank_3::{TensorRank3, get_identity_1010_parts},
26};
27
28pub(crate) mod list;
29pub(crate) mod vec;
30
31impl<const D: usize, I, J, K, L, U> HessianBlock for TensorRank4<D, I, J, K, L, U> {
32    fn entry(&self, row: usize, column: usize) -> TensorRank0 {
33        self[row / D][row % D][column / D][column % D].value()
34    }
35    fn height(&self) -> usize {
36        D * D
37    }
38    fn width(&self) -> usize {
39        D * D
40    }
41    fn fill_into_block<M>(&self, matrix: &mut M, row: usize, column: usize)
42    where
43        M: IndexMut<usize, Output = Vector>,
44    {
45        self.iter().enumerate().for_each(|(i, self_i)| {
46            self_i.iter().enumerate().for_each(|(j, self_ij)| {
47                let matrix_row = &mut matrix[row + D * i + j];
48                self_ij.iter().enumerate().for_each(|(k, self_ijk)| {
49                    self_ijk.iter().enumerate().for_each(|(l, self_ijkl)| {
50                        matrix_row[column + D * k + l] = self_ijkl.value()
51                    })
52                })
53            })
54        })
55    }
56}
57
58/// A *d*-dimensional tensor of rank 4.
59///
60/// `D` is the dimension, `I`, `J`, `K`, `L` are the configurations.
61#[repr(transparent)]
62pub struct TensorRank4<const D: usize, I, J, K, L, U = Dimensionless>(
63    [TensorRank3<D, J, K, L, U>; D],
64    pub(super) PhantomData<I>,
65);
66
67impl<const D: usize, I, J, K, L, U> Clone for TensorRank4<D, I, J, K, L, U> {
68    fn clone(&self) -> Self {
69        Self(self.0.clone(), PhantomData)
70    }
71}
72
73impl<const D: usize, I, J, K, L, U> Debug for TensorRank4<D, I, J, K, L, U> {
74    fn fmt(&self, f: &mut Formatter) -> fmt::Result {
75        self.0.fmt(f)
76    }
77}
78
79impl<const D: usize, I, J, K, L, U> PartialEq for TensorRank4<D, I, J, K, L, U> {
80    fn eq(&self, other: &Self) -> bool {
81        self.0 == other.0
82    }
83}
84
85impl<const D: usize, I, J, K, L, U> TensorRank4<D, I, J, K, L, U> {
86    pub fn with_unit<V>(self) -> TensorRank4<D, I, J, K, L, V> {
87        relabel(self.into_canonical())
88    }
89    fn canonical(
90        &self,
91    ) -> &TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless> {
92        unsafe {
93            &*(self as *const Self
94                as *const TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless>)
95        }
96    }
97    fn canonical_mut(
98        &mut self,
99    ) -> &mut TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless> {
100        unsafe {
101            &mut *(self as *mut Self
102                as *mut TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless>)
103        }
104    }
105    fn into_canonical(
106        self,
107    ) -> TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless> {
108        unsafe {
109            (&self as *const Self)
110                .cast::<TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless>>()
111                .read()
112        }
113    }
114}
115
116fn relabel<const D: usize, I, J, K, L, U>(
117    tensor: TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless>,
118) -> TensorRank4<D, I, J, K, L, U> {
119    unsafe {
120        (&tensor
121            as *const TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless>)
122            .cast::<TensorRank4<D, I, J, K, L, U>>()
123            .read()
124    }
125}
126
127impl<const D: usize> TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless> {
128    fn as_array_core(&self) -> [[[[TensorRank0; D]; D]; D]; D] {
129        let mut array = [[[[0.0; D]; D]; D]; D];
130        array
131            .iter_mut()
132            .zip(self.iter())
133            .for_each(|(entry_rank_3, tensor_rank_3)| *entry_rank_3 = tensor_rank_3.as_array());
134        array
135    }
136    fn zero_core() -> Self {
137        Self(from_fn(|_| TensorRank3::zero()), PhantomData)
138    }
139    fn add_assign_core(&mut self, tensor: Self) {
140        self.iter_mut()
141            .zip(tensor)
142            .for_each(|(self_i, tensor_i)| *self_i += tensor_i);
143    }
144    fn add_assign_ref_core(&mut self, tensor: &Self) {
145        self.iter_mut()
146            .zip(tensor.iter())
147            .for_each(|(self_i, tensor_i)| *self_i += tensor_i);
148    }
149    fn sub_assign_core(&mut self, tensor: Self) {
150        self.iter_mut()
151            .zip(tensor)
152            .for_each(|(self_i, tensor_i)| *self_i -= tensor_i);
153    }
154    fn sub_assign_ref_core(&mut self, tensor: &Self) {
155        self.iter_mut()
156            .zip(tensor.iter())
157            .for_each(|(self_i, tensor_i)| *self_i -= tensor_i);
158    }
159}
160
161impl<const D: usize, I, J, K, L, U> Default for TensorRank4<D, I, J, K, L, U> {
162    fn default() -> Self {
163        Self::zero()
164    }
165}
166
167impl<const D: usize, I, J, K, L, U> From<[[[[TensorRank0; D]; D]; D]; D]>
168    for TensorRank4<D, I, J, K, L, U>
169{
170    fn from(array: [[[[TensorRank0; D]; D]; D]; D]) -> Self {
171        array.into_iter().map(|entry| entry.into()).collect()
172    }
173}
174
175/// The 3D rank-4 identity, configurations (1, 0, 1, 0).
176pub const IDENTITY_1010: TensorRank4<3, Current, Reference, Current, Reference, Dimensionless> =
177    TensorRank4(get_identity_1010_parts(), PhantomData);
178
179impl<I, J, K, U> From<TensorRank4<3, I, J, K, Reference, U>>
180    for TensorRank4<3, I, J, K, Intermediate, U>
181{
182    fn from(tensor_rank_4: TensorRank4<3, I, J, K, Reference, U>) -> Self {
183        unsafe {
184            transmute::<
185                TensorRank4<3, I, J, K, Reference, U>,
186                TensorRank4<3, I, J, K, Intermediate, U>,
187            >(tensor_rank_4)
188        }
189    }
190}
191
192impl<J, L, U> From<TensorRank4<3, Current, J, Current, L, U>>
193    for TensorRank4<3, Intermediate, J, Intermediate, L, U>
194{
195    fn from(tensor_rank_4: TensorRank4<3, Current, J, Current, L, U>) -> Self {
196        unsafe {
197            transmute::<
198                TensorRank4<3, Current, J, Current, L, U>,
199                TensorRank4<3, Intermediate, J, Intermediate, L, U>,
200            >(tensor_rank_4)
201        }
202    }
203}
204
205impl<I, K, U> From<TensorRank4<3, I, Reference, K, Reference, U>>
206    for TensorRank4<3, I, Intermediate, K, Intermediate, U>
207{
208    fn from(tensor_rank_4: TensorRank4<3, I, Reference, K, Reference, U>) -> Self {
209        unsafe {
210            transmute::<
211                TensorRank4<3, I, Reference, K, Reference, U>,
212                TensorRank4<3, I, Intermediate, K, Intermediate, U>,
213            >(tensor_rank_4)
214        }
215    }
216}
217
218impl<K, U> From<TensorRank4<3, Reference, Reference, K, Reference, U>>
219    for TensorRank4<3, Intermediate, Intermediate, K, Intermediate, U>
220{
221    fn from(tensor_rank_4: TensorRank4<3, Reference, Reference, K, Reference, U>) -> Self {
222        unsafe {
223            transmute::<
224                TensorRank4<3, Reference, Reference, K, Reference, U>,
225                TensorRank4<3, Intermediate, Intermediate, K, Intermediate, U>,
226            >(tensor_rank_4)
227        }
228    }
229}
230
231impl<const D: usize, I, J, K, L, U> From<Vec<Vec<Vec<Vec<TensorRank0>>>>>
232    for TensorRank4<D, I, J, K, L, U>
233{
234    fn from(vec_rank_4: Vec<Vec<Vec<Vec<TensorRank0>>>>) -> Self {
235        vec_rank_4
236            .into_iter()
237            .map(|vec_rank_3| {
238                vec_rank_3
239                    .into_iter()
240                    .map(|vec_rank_2| {
241                        vec_rank_2
242                            .into_iter()
243                            .map(|vec_rank_1| vec_rank_1.into_iter().collect())
244                            .collect()
245                    })
246                    .collect()
247            })
248            .collect()
249    }
250}
251
252impl<const D: usize, I, J, K, L, U> From<TensorRank4<D, I, J, K, L, U>>
253    for Vec<Vec<Vec<Vec<TensorRank0>>>>
254{
255    fn from(tensor_rank_4: TensorRank4<D, I, J, K, L, U>) -> Self {
256        tensor_rank_4
257            .iter()
258            .map(|tensor_rank_3| {
259                tensor_rank_3
260                    .iter()
261                    .map(|tensor_rank_2| {
262                        tensor_rank_2
263                            .iter()
264                            .map(|tensor_rank_1| {
265                                tensor_rank_1.iter().map(|entry| entry.value()).collect()
266                            })
267                            .collect()
268                    })
269                    .collect()
270            })
271            .collect()
272    }
273}
274
275impl<const D: usize, I, J, K, L, U> From<TensorRank4<D, I, J, K, L, U>> for Vec<TensorRank0> {
276    fn from(tensor_rank_4: TensorRank4<D, I, J, K, L, U>) -> Self {
277        tensor_rank_4
278            .iter()
279            .flat_map(|tensor_rank_3| {
280                tensor_rank_3.iter().flat_map(|tensor_rank_2| {
281                    tensor_rank_2
282                        .iter()
283                        .flat_map(|tensor_rank_1| tensor_rank_1.iter().map(|entry| entry.value()))
284                })
285            })
286            .collect()
287    }
288}
289
290impl<const D: usize, I, J, K, L, U> From<TensorRank4<D, I, J, K, L, U>> for Vector {
291    fn from(tensor_rank_4: TensorRank4<D, I, J, K, L, U>) -> Self {
292        tensor_rank_4
293            .iter()
294            .flat_map(|tensor_rank_3| {
295                tensor_rank_3.iter().flat_map(|tensor_rank_2| {
296                    tensor_rank_2
297                        .iter()
298                        .flat_map(|tensor_rank_1| tensor_rank_1.iter().map(|entry| entry.value()))
299                })
300            })
301            .collect()
302    }
303}
304
305impl<const D: usize, I, J, K, L, U> Display for TensorRank4<D, I, J, K, L, U> {
306    fn fmt(&self, f: &mut Formatter) -> fmt::Result {
307        write!(f, "[")?;
308        self.iter()
309            .enumerate()
310            .try_for_each(|(i, entry)| write!(f, "{entry},\n\x1B[u\x1B[{}B\x1B[2D", i + 1))?;
311        write!(f, "\x1B[u\x1B[1A\x1B[{}C]", 16 * D + 2)
312    }
313}
314
315impl<const D: usize, I, J, K, L, U> FiniteDifference for TensorRank4<D, I, J, K, L, U> {
316    fn error_fd(&self, comparator: &Self, epsilon: TensorRank0) -> Option<(bool, usize)> {
317        let error_count = self
318            .iter()
319            .zip(comparator.iter())
320            .map(|(self_i, comparator_i)| {
321                self_i
322                    .iter()
323                    .zip(comparator_i.iter())
324                    .map(|(self_ij, comparator_ij)| {
325                        self_ij
326                            .iter()
327                            .zip(comparator_ij.iter())
328                            .map(|(self_ijk, comparator_ijk)| {
329                                self_ijk
330                                    .iter()
331                                    .zip(comparator_ijk.iter())
332                                    .filter(|&(&self_ijkl, &comparator_ijkl)| {
333                                        self_ijkl.differs(comparator_ijkl, epsilon)
334                                    })
335                                    .count()
336                            })
337                            .sum::<usize>()
338                    })
339                    .sum::<usize>()
340            })
341            .sum();
342        if error_count > 0 {
343            let auxiliary = self
344                .iter()
345                .zip(comparator.iter())
346                .map(|(self_i, comparator_i)| {
347                    self_i
348                        .iter()
349                        .zip(comparator_i.iter())
350                        .map(|(self_ij, comparator_ij)| {
351                            self_ij
352                                .iter()
353                                .zip(comparator_ij.iter())
354                                .map(|(self_ijk, comparator_ijk)| {
355                                    self_ijk
356                                        .iter()
357                                        .zip(comparator_ijk.iter())
358                                        .filter(|&(&self_ijkl, &comparator_ijkl)| {
359                                            self_ijkl.differs_severely(comparator_ijkl, epsilon)
360                                        })
361                                        .count()
362                                })
363                                .sum::<usize>()
364                        })
365                        .sum::<usize>()
366                })
367                .sum::<usize>()
368                > 0;
369            Some((auxiliary, error_count))
370        } else {
371            None
372        }
373    }
374}
375
376impl<const D: usize, I, J, K, L, U> TensorRank4<D, I, J, K, L, U> {
377    pub fn dyad_ij_kl<U1, U2>(
378        tensor_rank_2_a: &TensorRank2<D, I, J, U1>,
379        tensor_rank_2_b: &TensorRank2<D, K, L, U2>,
380    ) -> Self
381    where
382        U1: UnitMul<U2, Output = U>,
383    {
384        relabel(canonical_dyad_ij_kl(
385            tensor_rank_2_a.canonical(),
386            tensor_rank_2_b.canonical(),
387        ))
388    }
389    pub fn dyad_ik_jl<U1, U2>(
390        tensor_rank_2_a: &TensorRank2<D, I, K, U1>,
391        tensor_rank_2_b: &TensorRank2<D, J, L, U2>,
392    ) -> Self
393    where
394        U1: UnitMul<U2, Output = U>,
395    {
396        relabel(canonical_dyad_ik_jl(
397            tensor_rank_2_a.canonical(),
398            tensor_rank_2_b.canonical(),
399        ))
400    }
401    pub fn dyad_il_jk<U1, U2>(
402        tensor_rank_2_a: &TensorRank2<D, I, L, U1>,
403        tensor_rank_2_b: &TensorRank2<D, J, K, U2>,
404    ) -> Self
405    where
406        U1: UnitMul<U2, Output = U>,
407    {
408        relabel(canonical_dyad_il_jk(
409            tensor_rank_2_a.canonical(),
410            tensor_rank_2_b.canonical(),
411        ))
412    }
413    pub fn dyad_il_kj<U1, U2>(
414        tensor_rank_2_a: &TensorRank2<D, I, L, U1>,
415        tensor_rank_2_b: &TensorRank2<D, K, J, U2>,
416    ) -> Self
417    where
418        U1: UnitMul<U2, Output = U>,
419    {
420        Self::dyad_il_jk(tensor_rank_2_a, &(tensor_rank_2_b.transpose()))
421    }
422}
423
424impl<const D: usize, I, J, K, L, U> Hessian for TensorRank4<D, I, J, K, L, U> {
425    fn entry(&self, row: usize, column: usize) -> TensorRank0 {
426        self[row / D][row % D][column / D][column % D].value()
427    }
428    fn quadratic_form(&self, vector: &Vector) -> TensorRank0 {
429        self.iter()
430            .enumerate()
431            .map(|(i, self_i)| {
432                self_i
433                    .iter()
434                    .enumerate()
435                    .map(|(j, self_ij)| {
436                        vector[D * i + j]
437                            * self_ij
438                                .iter()
439                                .enumerate()
440                                .map(|(k, self_ijk)| {
441                                    self_ijk
442                                        .iter()
443                                        .enumerate()
444                                        .map(|(l, self_ijkl)| self_ijkl.value() * vector[D * k + l])
445                                        .sum::<TensorRank0>()
446                                })
447                                .sum::<TensorRank0>()
448                    })
449                    .sum::<TensorRank0>()
450            })
451            .sum()
452    }
453    fn fill_into(self, square_matrix: &mut SquareMatrix) {
454        self.into_iter().enumerate().for_each(|(i, self_i)| {
455            self_i.into_iter().enumerate().for_each(|(j, self_ij)| {
456                self_ij.into_iter().enumerate().for_each(|(k, self_ijk)| {
457                    self_ijk.into_iter().enumerate().for_each(|(l, self_ijkl)| {
458                        square_matrix[D * i + j][D * k + l] = self_ijkl.value()
459                    })
460                })
461            })
462        })
463    }
464    fn retain_from(self, retained: &[bool]) -> SquareMatrix {
465        (0..D * D)
466            .filter(|&row| retained[row])
467            .map(|row| {
468                (0..D * D)
469                    .filter(|&column| retained[column])
470                    .map(|column| self[row / D][row % D][column / D][column % D].value())
471                    .collect()
472            })
473            .collect()
474    }
475}
476
477impl<const D: usize, I, J, K, L, U> Erase for TensorRank4<D, I, J, K, L, U> {
478    type Erased = TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless>;
479    fn erase(&self) -> &Self::Erased {
480        self.canonical()
481    }
482}
483
484impl<const D: usize, I, J, K, L, U> Tensor for TensorRank4<D, I, J, K, L, U> {
485    type Item = TensorRank3<D, J, K, L, U>;
486    type Unit = U;
487    fn iter(&self) -> impl Iterator<Item = &Self::Item> {
488        self.0.iter()
489    }
490    fn iter_mut(&mut self) -> impl Iterator<Item = &mut Self::Item> {
491        self.0.iter_mut()
492    }
493    fn len(&self) -> usize {
494        D
495    }
496    fn size(&self) -> usize {
497        D * D * D * D
498    }
499}
500
501impl<const D: usize, I, J, K, L, U> IntoIterator for TensorRank4<D, I, J, K, L, U> {
502    type Item = TensorRank3<D, J, K, L, U>;
503    type IntoIter = std::array::IntoIter<Self::Item, D>;
504    fn into_iter(self) -> Self::IntoIter {
505        self.0.into_iter()
506    }
507}
508
509impl<const D: usize, I, J, K, L, U> TensorArray for TensorRank4<D, I, J, K, L, U> {
510    type Array = [[[[TensorRank0; D]; D]; D]; D];
511    type Item = TensorRank3<D, J, K, L, U>;
512    fn as_array(&self) -> Self::Array {
513        self.canonical().as_array_core()
514    }
515    fn identity() -> Self {
516        relabel(canonical_dyad_ij_kl(
517            &TensorRank2::identity(),
518            &TensorRank2::identity(),
519        ))
520    }
521    fn zero() -> Self {
522        relabel(TensorRank4::<
523            D,
524            Reference,
525            Reference,
526            Reference,
527            Reference,
528            Dimensionless,
529        >::zero_core())
530    }
531}
532
533impl<const D: usize, I, J, K, L, U> FromIterator<TensorRank3<D, J, K, L, U>>
534    for TensorRank4<D, I, J, K, L, U>
535{
536    fn from_iter<Ii: IntoIterator<Item = TensorRank3<D, J, K, L, U>>>(into_iterator: Ii) -> Self {
537        let mut tensor_rank_4 = Self::zero();
538        tensor_rank_4
539            .iter_mut()
540            .zip(into_iterator)
541            .for_each(|(tensor_rank_4_i, value_i)| *tensor_rank_4_i = value_i);
542        tensor_rank_4
543    }
544}
545
546impl<const D: usize, I, J, K, L, U> Index<usize> for TensorRank4<D, I, J, K, L, U> {
547    type Output = TensorRank3<D, J, K, L, U>;
548    fn index(&self, index: usize) -> &Self::Output {
549        &self.0[index]
550    }
551}
552
553impl<const D: usize, I, J, K, L, U> IndexMut<usize> for TensorRank4<D, I, J, K, L, U> {
554    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
555        &mut self.0[index]
556    }
557}
558
559impl<const D: usize, I, J, K, L, U> Sum for TensorRank4<D, I, J, K, L, U> {
560    fn sum<Ii>(iter: Ii) -> Self
561    where
562        Ii: Iterator<Item = Self>,
563    {
564        iter.reduce(|mut acc, item| {
565            acc += item;
566            acc
567        })
568        .unwrap_or_else(Self::default)
569    }
570}
571
572impl<'a, const D: usize, I, J, K, L, U> Sum<&'a Self> for TensorRank4<D, I, J, K, L, U> {
573    fn sum<Ii>(iter: Ii) -> Self
574    where
575        Ii: Iterator<Item = &'a Self>,
576    {
577        iter.fold(Self::default(), |mut acc, item| {
578            acc += item;
579            acc
580        })
581    }
582}
583
584/// Transforms all four indices of a rank-4 tensor.
585///
586/// Contracts each index of `self` with the first index of a corresponding rank-2 tensor.
587pub trait ContractAllWithFirst<TIM, TJN, TKO, TLP> {
588    type Output;
589    fn contract_all_with_first(
590        self,
591        object_a: TIM,
592        object_b: TJN,
593        object_c: TKO,
594        object_d: TLP,
595    ) -> Self::Output;
596}
597
598impl<const D: usize, I, J, K, L, M, N, O, P, U>
599    ContractAllWithFirst<
600        &TensorRank2<D, I, M, Dimensionless>,
601        &TensorRank2<D, J, N, Dimensionless>,
602        &TensorRank2<D, K, O, Dimensionless>,
603        &TensorRank2<D, L, P, Dimensionless>,
604    > for TensorRank4<D, I, J, K, L, U>
605{
606    type Output = TensorRank4<D, M, N, O, P, U>;
607    fn contract_all_with_first(
608        self,
609        tensor_rank_2_a: &TensorRank2<D, I, M, Dimensionless>,
610        tensor_rank_2_b: &TensorRank2<D, J, N, Dimensionless>,
611        tensor_rank_2_c: &TensorRank2<D, K, O, Dimensionless>,
612        tensor_rank_2_d: &TensorRank2<D, L, P, Dimensionless>,
613    ) -> Self::Output {
614        let first = canonical_transform_first(self.canonical(), tensor_rank_2_a.canonical());
615        let second = canonical_contract_second_with_first(&first, tensor_rank_2_b.canonical());
616        let third = canonical_contract_third_with_first(&second, tensor_rank_2_c.canonical());
617        relabel(canonical_transform_fourth(
618            &third,
619            tensor_rank_2_d.canonical(),
620        ))
621    }
622}
623
624pub trait ContractFirstThirdFourthWithFirst<TIM, TKO, TLP> {
625    type Output;
626    fn contract_first_third_fourth_with_first(
627        self,
628        object_a: TIM,
629        object_b: TKO,
630        object_c: TLP,
631    ) -> Self::Output;
632}
633
634impl<const D: usize, I, J, K, L, M, O, P, U>
635    ContractFirstThirdFourthWithFirst<
636        &TensorRank2<D, I, M, Dimensionless>,
637        &TensorRank2<D, K, O, Dimensionless>,
638        &TensorRank2<D, L, P, Dimensionless>,
639    > for TensorRank4<D, I, J, K, L, U>
640{
641    type Output = TensorRank4<D, M, J, O, P, U>;
642    fn contract_first_third_fourth_with_first(
643        self,
644        tensor_rank_2_a: &TensorRank2<D, I, M, Dimensionless>,
645        tensor_rank_2_b: &TensorRank2<D, K, O, Dimensionless>,
646        tensor_rank_2_c: &TensorRank2<D, L, P, Dimensionless>,
647    ) -> Self::Output {
648        let first = canonical_transform_first(self.canonical(), tensor_rank_2_a.canonical());
649        let third = canonical_contract_third_with_first(&first, tensor_rank_2_b.canonical());
650        relabel(canonical_transform_fourth(
651            &third,
652            tensor_rank_2_c.canonical(),
653        ))
654    }
655}
656
657pub trait ContractSecondWithFirst<TJN> {
658    type Output;
659    fn contract_second_with_first(self, tensor_rank_2: TJN) -> Self::Output;
660}
661
662impl<const D: usize, I, J, K, L, N, U> ContractSecondWithFirst<&TensorRank2<D, J, N, Dimensionless>>
663    for TensorRank4<D, I, J, K, L, U>
664{
665    type Output = TensorRank4<D, I, N, K, L, U>;
666    fn contract_second_with_first(
667        self,
668        tensor_rank_2: &TensorRank2<D, J, N, Dimensionless>,
669    ) -> Self::Output {
670        relabel(canonical_contract_second_with_first(
671            self.canonical(),
672            tensor_rank_2.canonical(),
673        ))
674    }
675}
676
677pub trait ContractSecondFourthWithFirst<TJ, TL> {
678    type Output;
679    fn contract_second_fourth_with_first(&self, object_a: TJ, object_b: TL) -> Self::Output;
680}
681
682impl<const D: usize, I, J, K, L, U, V, W>
683    ContractSecondFourthWithFirst<&TensorRank1<D, J, V>, &TensorRank1<D, L, W>>
684    for TensorRank4<D, I, J, K, L, U>
685where
686    U: UnitMul<V>,
687    <U as UnitMul<V>>::Output: UnitMul<W>,
688{
689    type Output = TensorRank2<D, I, K, <<U as UnitMul<V>>::Output as UnitMul<W>>::Output>;
690    fn contract_second_fourth_with_first(
691        &self,
692        tensor_rank_1_a: &TensorRank1<D, J, V>,
693        tensor_rank_1_b: &TensorRank1<D, L, W>,
694    ) -> Self::Output {
695        relabel_rank_2(canonical_contract_second_fourth_with_first(
696            self.canonical(),
697            tensor_rank_1_a.canonical(),
698            tensor_rank_1_b.canonical(),
699        ))
700    }
701}
702
703fn canonical_contract_second_fourth_with_first<const D: usize>(
704    tensor_rank_4: &TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless>,
705    tensor_rank_1_a: &TensorRank1<D, Reference, Dimensionless>,
706    tensor_rank_1_b: &TensorRank1<D, Reference, Dimensionless>,
707) -> TensorRank2<D, Reference, Reference, Dimensionless> {
708    let mut output = TensorRank2::zero();
709    tensor_rank_4
710        .iter()
711        .zip(output.iter_mut())
712        .for_each(|(tensor_rank_4_i, output_i)| {
713            tensor_rank_4_i.iter().zip(tensor_rank_1_a.iter()).for_each(
714                |(tensor_rank_4_ij, tensor_rank_1_a_j)| {
715                    tensor_rank_4_ij.iter().zip(output_i.iter_mut()).for_each(
716                        |(tensor_rank_4_ijk, output_ik)| {
717                            *output_ik += (tensor_rank_4_ijk * tensor_rank_1_b) * tensor_rank_1_a_j
718                        },
719                    )
720                },
721            )
722        });
723    output
724}
725
726pub trait ContractThirdWithFirst<TKL> {
727    type Output;
728    fn contract_third_with_first(&self, tensor: TKL) -> Self::Output;
729}
730
731impl<const D: usize, I, J, K, L, M, U> ContractThirdWithFirst<&TensorRank2<D, M, K, Dimensionless>>
732    for TensorRank4<D, I, J, M, L, U>
733{
734    type Output = TensorRank4<D, I, J, K, L, U>;
735    fn contract_third_with_first(
736        &self,
737        tensor_rank_2: &TensorRank2<D, M, K, Dimensionless>,
738    ) -> Self::Output {
739        relabel(canonical_contract_third_with_first(
740            self.canonical(),
741            tensor_rank_2.canonical(),
742        ))
743    }
744}
745
746pub trait ContractThirdFourthWithFirstSecond<TKL> {
747    type Output;
748    fn contract_third_fourth_with_first_second(self, tensor: TKL) -> Self::Output;
749}
750
751impl<const D: usize, I, J, K, L, U, V> ContractThirdFourthWithFirstSecond<&TensorRank2<D, K, L, V>>
752    for TensorRank4<D, I, J, K, L, U>
753where
754    U: UnitMul<V>,
755{
756    type Output = TensorRank2<D, I, J, <U as UnitMul<V>>::Output>;
757    fn contract_third_fourth_with_first_second(
758        self,
759        tensor_rank_2: &TensorRank2<D, K, L, V>,
760    ) -> Self::Output {
761        relabel_rank_2(canonical_contract_34_12_rank_2(
762            self.into_canonical(),
763            tensor_rank_2.canonical(),
764        ))
765    }
766}
767
768fn canonical_contract_34_12_rank_2<const D: usize>(
769    tensor_rank_4: TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless>,
770    tensor_rank_2: &TensorRank2<D, Reference, Reference, Dimensionless>,
771) -> TensorRank2<D, Reference, Reference, Dimensionless> {
772    tensor_rank_4
773        .into_iter()
774        .map(|self_i| {
775            self_i
776                .into_iter()
777                .map(|self_ij| self_ij.full_contraction(tensor_rank_2))
778                .collect()
779        })
780        .collect()
781}
782
783impl<const D: usize, I, J, K, L, M, N, U, V>
784    ContractThirdFourthWithFirstSecond<&TensorRank4<D, K, L, M, N, V>>
785    for TensorRank4<D, I, J, K, L, U>
786where
787    U: UnitMul<V>,
788{
789    type Output = TensorRank4<D, I, J, M, N, <U as UnitMul<V>>::Output>;
790    fn contract_third_fourth_with_first_second(
791        self,
792        tensor: &TensorRank4<D, K, L, M, N, V>,
793    ) -> Self::Output {
794        relabel(canonical_contract_34_12_rank_4(
795            self.into_canonical(),
796            tensor.canonical(),
797        ))
798    }
799}
800
801fn canonical_contract_34_12_rank_4<const D: usize>(
802    tensor_rank_4: TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless>,
803    tensor: &TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless>,
804) -> TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless> {
805    tensor_rank_4
806        .into_iter()
807        .map(|self_i| {
808            self_i
809                .into_iter()
810                .map(|self_ij| {
811                    self_ij
812                        .into_iter()
813                        .zip(tensor.iter())
814                        .map(|(self_ijk, tensor_k)| {
815                            self_ijk
816                                .into_iter()
817                                .zip(tensor_k.iter())
818                                .map(|(self_ijkl, tensor_kl)| tensor_kl * self_ijkl)
819                                .sum::<TensorRank2<D, Reference, Reference, Dimensionless>>()
820                        })
821                        .sum()
822                })
823                .collect()
824        })
825        .collect()
826}
827
828pub trait ContractFirstSecondWithSecond<TI, TJ> {
829    type Output;
830    fn contract_first_second_with_second(self, object_a: TI, object_b: TJ) -> Self::Output;
831}
832
833impl<const D: usize, I, J, K, L, M, N, U>
834    ContractFirstSecondWithSecond<
835        &TensorRank2<D, I, M, Dimensionless>,
836        &TensorRank2<D, J, N, Dimensionless>,
837    > for TensorRank4<D, M, N, K, L, U>
838{
839    type Output = TensorRank4<D, I, J, K, L, U>;
840    fn contract_first_second_with_second(
841        self,
842        tensor_rank_2_a: &TensorRank2<D, I, M, Dimensionless>,
843        tensor_rank_2_b: &TensorRank2<D, J, N, Dimensionless>,
844    ) -> Self::Output {
845        let first =
846            canonical_transform_first(self.canonical(), &tensor_rank_2_a.canonical().transpose());
847        relabel(canonical_contract_second_with_first(
848            &first,
849            &tensor_rank_2_b.canonical().transpose(),
850        ))
851    }
852}
853
854impl<const D: usize, I, J, K, L, U> Div<TensorRank0> for TensorRank4<D, I, J, K, L, U> {
855    type Output = Self;
856    fn div(mut self, tensor_rank_0: TensorRank0) -> Self::Output {
857        self /= &tensor_rank_0;
858        self
859    }
860}
861
862impl<const D: usize, I, J, K, L, U> Div<TensorRank0> for &TensorRank4<D, I, J, K, L, U> {
863    type Output = TensorRank4<D, I, J, K, L, U>;
864    fn div(self, tensor_rank_0: TensorRank0) -> Self::Output {
865        self.iter().map(|self_i| self_i / tensor_rank_0).collect()
866    }
867}
868
869impl<const D: usize, I, J, K, L, U> Div<&TensorRank0> for TensorRank4<D, I, J, K, L, U> {
870    type Output = Self;
871    fn div(mut self, tensor_rank_0: &TensorRank0) -> Self::Output {
872        self /= tensor_rank_0;
873        self
874    }
875}
876
877impl<const D: usize, I, J, K, L, U> DivAssign<TensorRank0> for TensorRank4<D, I, J, K, L, U> {
878    fn div_assign(&mut self, tensor_rank_0: TensorRank0) {
879        self.iter_mut().for_each(|self_i| *self_i /= &tensor_rank_0);
880    }
881}
882
883impl<const D: usize, I, J, K, L, U> DivAssign<&TensorRank0> for TensorRank4<D, I, J, K, L, U> {
884    fn div_assign(&mut self, tensor_rank_0: &TensorRank0) {
885        self.iter_mut().for_each(|self_i| *self_i /= tensor_rank_0);
886    }
887}
888
889impl<const D: usize, I, J, K, L, U> Mul<TensorRank0> for TensorRank4<D, I, J, K, L, U> {
890    type Output = Self;
891    fn mul(mut self, tensor_rank_0: TensorRank0) -> Self::Output {
892        self *= &tensor_rank_0;
893        self
894    }
895}
896
897impl<const D: usize, I, J, K, L, U> Mul<&TensorRank0> for TensorRank4<D, I, J, K, L, U> {
898    type Output = Self;
899    fn mul(mut self, tensor_rank_0: &TensorRank0) -> Self::Output {
900        self *= tensor_rank_0;
901        self
902    }
903}
904
905impl<const D: usize, I, J, K, L, U> MulAssign<TensorRank0> for TensorRank4<D, I, J, K, L, U> {
906    fn mul_assign(&mut self, tensor_rank_0: TensorRank0) {
907        self.iter_mut().for_each(|self_i| *self_i *= &tensor_rank_0);
908    }
909}
910
911impl<const D: usize, I, J, K, L, U> MulAssign<&TensorRank0> for TensorRank4<D, I, J, K, L, U> {
912    fn mul_assign(&mut self, tensor_rank_0: &TensorRank0) {
913        self.iter_mut().for_each(|self_i| *self_i *= tensor_rank_0);
914    }
915}
916
917impl<const D: usize, I, J, K, L, M, U, V> Mul<TensorRank2<D, L, M, V>>
918    for TensorRank4<D, I, J, K, L, U>
919where
920    U: UnitMul<V>,
921{
922    type Output = TensorRank4<D, I, J, K, M, <U as UnitMul<V>>::Output>;
923    fn mul(self, tensor_rank_2: TensorRank2<D, L, M, V>) -> Self::Output {
924        self.into_iter()
925            .map(|self_i| {
926                self_i
927                    .into_iter()
928                    .map(|self_ij| self_ij * &tensor_rank_2)
929                    .collect()
930            })
931            .collect()
932    }
933}
934
935impl<const D: usize, I, J, K, L, M, U, V> Mul<&TensorRank2<D, L, M, V>>
936    for TensorRank4<D, I, J, K, L, U>
937where
938    U: UnitMul<V>,
939{
940    type Output = TensorRank4<D, I, J, K, M, <U as UnitMul<V>>::Output>;
941    fn mul(self, tensor_rank_2: &TensorRank2<D, L, M, V>) -> Self::Output {
942        self.into_iter()
943            .map(|self_i| {
944                self_i
945                    .into_iter()
946                    .map(|self_ij| self_ij * tensor_rank_2)
947                    .collect()
948            })
949            .collect()
950    }
951}
952
953impl<const D: usize, J, K, L, M, U, V> Mul<TensorRank4<D, M, J, K, L, V>> for TensorRank1<D, M, U>
954where
955    U: UnitMul<V>,
956{
957    type Output = TensorRank3<D, J, K, L, <U as UnitMul<V>>::Output>;
958    fn mul(self, tensor_rank_4: TensorRank4<D, M, J, K, L, V>) -> Self::Output {
959        relabel_rank_3(canonical_rank_1_times_rank_4(
960            self.canonical(),
961            tensor_rank_4.canonical(),
962        ))
963    }
964}
965
966impl<const D: usize, J, K, L, M, U, V> Mul<&TensorRank4<D, M, J, K, L, V>> for TensorRank1<D, M, U>
967where
968    U: UnitMul<V>,
969{
970    type Output = TensorRank3<D, J, K, L, <U as UnitMul<V>>::Output>;
971    fn mul(self, tensor_rank_4: &TensorRank4<D, M, J, K, L, V>) -> Self::Output {
972        relabel_rank_3(canonical_rank_1_times_rank_4(
973            self.canonical(),
974            tensor_rank_4.canonical(),
975        ))
976    }
977}
978
979impl<const D: usize, J, K, L, M, U, V> Mul<TensorRank4<D, M, J, K, L, V>> for &TensorRank1<D, M, U>
980where
981    U: UnitMul<V>,
982{
983    type Output = TensorRank3<D, J, K, L, <U as UnitMul<V>>::Output>;
984    fn mul(self, tensor_rank_4: TensorRank4<D, M, J, K, L, V>) -> Self::Output {
985        relabel_rank_3(canonical_rank_1_times_rank_4(
986            self.canonical(),
987            tensor_rank_4.canonical(),
988        ))
989    }
990}
991
992impl<const D: usize, J, K, L, M, U, V> Mul<&TensorRank4<D, M, J, K, L, V>> for &TensorRank1<D, M, U>
993where
994    U: UnitMul<V>,
995{
996    type Output = TensorRank3<D, J, K, L, <U as UnitMul<V>>::Output>;
997    fn mul(self, tensor_rank_4: &TensorRank4<D, M, J, K, L, V>) -> Self::Output {
998        relabel_rank_3(canonical_rank_1_times_rank_4(
999            self.canonical(),
1000            tensor_rank_4.canonical(),
1001        ))
1002    }
1003}
1004
1005impl<const D: usize, I, J, K, L, M, U, V> Mul<TensorRank4<D, M, J, K, L, V>>
1006    for TensorRank2<D, I, M, U>
1007where
1008    U: UnitMul<V>,
1009{
1010    type Output = TensorRank4<D, I, J, K, L, <U as UnitMul<V>>::Output>;
1011    fn mul(self, tensor_rank_4: TensorRank4<D, M, J, K, L, V>) -> Self::Output {
1012        relabel(canonical_rank_2_times_rank_4(
1013            self.canonical(),
1014            tensor_rank_4.canonical(),
1015        ))
1016    }
1017}
1018
1019impl<const D: usize, I, J, K, L, M, U, V> Mul<&TensorRank4<D, M, J, K, L, V>>
1020    for TensorRank2<D, I, M, U>
1021where
1022    U: UnitMul<V>,
1023{
1024    type Output = TensorRank4<D, I, J, K, L, <U as UnitMul<V>>::Output>;
1025    fn mul(self, tensor_rank_4: &TensorRank4<D, M, J, K, L, V>) -> Self::Output {
1026        relabel(canonical_rank_2_times_rank_4(
1027            self.canonical(),
1028            tensor_rank_4.canonical(),
1029        ))
1030    }
1031}
1032
1033impl<const D: usize, I, J, K, L, M, U, V> Mul<TensorRank4<D, M, J, K, L, V>>
1034    for &TensorRank2<D, I, M, U>
1035where
1036    U: UnitMul<V>,
1037{
1038    type Output = TensorRank4<D, I, J, K, L, <U as UnitMul<V>>::Output>;
1039    fn mul(self, tensor_rank_4: TensorRank4<D, M, J, K, L, V>) -> Self::Output {
1040        relabel(canonical_rank_2_times_rank_4(
1041            self.canonical(),
1042            tensor_rank_4.canonical(),
1043        ))
1044    }
1045}
1046
1047impl<const D: usize, I, J, K, L, M, U, V> Mul<&TensorRank4<D, M, J, K, L, V>>
1048    for &TensorRank2<D, I, M, U>
1049where
1050    U: UnitMul<V>,
1051{
1052    type Output = TensorRank4<D, I, J, K, L, <U as UnitMul<V>>::Output>;
1053    fn mul(self, tensor_rank_4: &TensorRank4<D, M, J, K, L, V>) -> Self::Output {
1054        relabel(canonical_rank_2_times_rank_4(
1055            self.canonical(),
1056            tensor_rank_4.canonical(),
1057        ))
1058    }
1059}
1060
1061impl<const D: usize, I, J, K, L, U> Add for TensorRank4<D, I, J, K, L, U> {
1062    type Output = Self;
1063    fn add(mut self, tensor_rank_4: Self) -> Self::Output {
1064        self += tensor_rank_4;
1065        self
1066    }
1067}
1068
1069impl<const D: usize, I, J, K, L, U> Add<&Self> for TensorRank4<D, I, J, K, L, U> {
1070    type Output = Self;
1071    fn add(mut self, tensor_rank_4: &Self) -> Self::Output {
1072        self += tensor_rank_4;
1073        self
1074    }
1075}
1076
1077impl<const D: usize, I, J, K, L, U> Add<TensorRank4<D, I, J, K, L, U>>
1078    for &TensorRank4<D, I, J, K, L, U>
1079{
1080    type Output = TensorRank4<D, I, J, K, L, U>;
1081    fn add(self, mut tensor_rank_4: TensorRank4<D, I, J, K, L, U>) -> Self::Output {
1082        tensor_rank_4 += self;
1083        tensor_rank_4
1084    }
1085}
1086
1087impl<const D: usize, I, J, K, L, U> AddAssign for TensorRank4<D, I, J, K, L, U> {
1088    fn add_assign(&mut self, tensor_rank_4: Self) {
1089        self.canonical_mut()
1090            .add_assign_core(tensor_rank_4.into_canonical());
1091    }
1092}
1093
1094impl<const D: usize, I, J, K, L, U> AddAssign<&Self> for TensorRank4<D, I, J, K, L, U> {
1095    fn add_assign(&mut self, tensor_rank_4: &Self) {
1096        self.canonical_mut()
1097            .add_assign_ref_core(tensor_rank_4.canonical());
1098    }
1099}
1100
1101impl<const D: usize, I, J, K, L, U> Sub for TensorRank4<D, I, J, K, L, U> {
1102    type Output = Self;
1103    fn sub(mut self, tensor_rank_4: Self) -> Self::Output {
1104        self -= tensor_rank_4;
1105        self
1106    }
1107}
1108
1109impl<const D: usize, I, J, K, L, U> Sub<&Self> for TensorRank4<D, I, J, K, L, U> {
1110    type Output = Self;
1111    fn sub(mut self, tensor_rank_4: &Self) -> Self::Output {
1112        self -= tensor_rank_4;
1113        self
1114    }
1115}
1116
1117impl<const D: usize, I, J, K, L, U> Sub for &TensorRank4<D, I, J, K, L, U> {
1118    type Output = TensorRank4<D, I, J, K, L, U>;
1119    fn sub(self, tensor_rank_4: Self) -> Self::Output {
1120        tensor_rank_4
1121            .iter()
1122            .zip(self.iter())
1123            .map(|(tensor_rank_4_i, self_i)| self_i - tensor_rank_4_i)
1124            .collect()
1125    }
1126}
1127
1128impl<const D: usize, I, J, K, L, U> SubAssign for TensorRank4<D, I, J, K, L, U> {
1129    fn sub_assign(&mut self, tensor_rank_4: Self) {
1130        self.canonical_mut()
1131            .sub_assign_core(tensor_rank_4.into_canonical());
1132    }
1133}
1134
1135impl<const D: usize, I, J, K, L, U> SubAssign<&Self> for TensorRank4<D, I, J, K, L, U> {
1136    fn sub_assign(&mut self, tensor_rank_4: &Self) {
1137        self.canonical_mut()
1138            .sub_assign_ref_core(tensor_rank_4.canonical());
1139    }
1140}
1141
1142impl<const D: usize, I, J, K, L, U, V> Mul<Quantity<V>> for TensorRank4<D, I, J, K, L, U>
1143where
1144    U: UnitMul<V>,
1145{
1146    type Output = TensorRank4<D, I, J, K, L, <U as UnitMul<V>>::Output>;
1147    fn mul(self, quantity: Quantity<V>) -> Self::Output {
1148        relabel(self.into_canonical() * quantity.value())
1149    }
1150}
1151
1152impl<const D: usize, I, J, K, L, U, V> Mul<&Quantity<V>> for TensorRank4<D, I, J, K, L, U>
1153where
1154    U: UnitMul<V>,
1155{
1156    type Output = TensorRank4<D, I, J, K, L, <U as UnitMul<V>>::Output>;
1157    fn mul(self, quantity: &Quantity<V>) -> Self::Output {
1158        self * *quantity
1159    }
1160}
1161
1162impl<const D: usize, I, J, K, L, U, V> Mul<Quantity<V>> for &TensorRank4<D, I, J, K, L, U>
1163where
1164    U: UnitMul<V>,
1165{
1166    type Output = TensorRank4<D, I, J, K, L, <U as UnitMul<V>>::Output>;
1167    fn mul(self, quantity: Quantity<V>) -> Self::Output {
1168        relabel(self.canonical().clone() * quantity.value())
1169    }
1170}
1171
1172impl<const D: usize, I, J, K, L, U, V> Mul<&Quantity<V>> for &TensorRank4<D, I, J, K, L, U>
1173where
1174    U: UnitMul<V>,
1175{
1176    type Output = TensorRank4<D, I, J, K, L, <U as UnitMul<V>>::Output>;
1177    fn mul(self, quantity: &Quantity<V>) -> Self::Output {
1178        self * *quantity
1179    }
1180}
1181
1182impl<const D: usize, I, J, K, L, U, V> Div<Quantity<V>> for TensorRank4<D, I, J, K, L, U>
1183where
1184    U: UnitDiv<V>,
1185{
1186    type Output = TensorRank4<D, I, J, K, L, <U as UnitDiv<V>>::Output>;
1187    fn div(self, quantity: Quantity<V>) -> Self::Output {
1188        relabel(self.into_canonical() / quantity.value())
1189    }
1190}
1191
1192fn canonical_dyad_ij_kl<const D: usize>(
1193    tensor_rank_2_a: &TensorRank2<D, Reference, Reference, Dimensionless>,
1194    tensor_rank_2_b: &TensorRank2<D, Reference, Reference, Dimensionless>,
1195) -> TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless> {
1196    tensor_rank_2_a
1197        .iter()
1198        .map(|tensor_rank_2_a_i| {
1199            tensor_rank_2_a_i
1200                .iter()
1201                .map(|tensor_rank_2_a_ij| tensor_rank_2_b * tensor_rank_2_a_ij)
1202                .collect()
1203        })
1204        .collect()
1205}
1206
1207fn canonical_dyad_ik_jl<const D: usize>(
1208    tensor_rank_2_a: &TensorRank2<D, Reference, Reference, Dimensionless>,
1209    tensor_rank_2_b: &TensorRank2<D, Reference, Reference, Dimensionless>,
1210) -> TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless> {
1211    tensor_rank_2_a
1212        .iter()
1213        .map(|tensor_rank_2_a_i| {
1214            tensor_rank_2_b
1215                .iter()
1216                .map(|tensor_rank_2_b_j| {
1217                    tensor_rank_2_a_i
1218                        .iter()
1219                        .map(|tensor_rank_2_a_ik| tensor_rank_2_b_j * tensor_rank_2_a_ik)
1220                        .collect()
1221                })
1222                .collect()
1223        })
1224        .collect()
1225}
1226
1227fn canonical_dyad_il_jk<const D: usize>(
1228    tensor_rank_2_a: &TensorRank2<D, Reference, Reference, Dimensionless>,
1229    tensor_rank_2_b: &TensorRank2<D, Reference, Reference, Dimensionless>,
1230) -> TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless> {
1231    tensor_rank_2_a
1232        .iter()
1233        .map(|tensor_rank_2_a_i| {
1234            tensor_rank_2_b
1235                .iter()
1236                .map(|tensor_rank_2_b_j| {
1237                    tensor_rank_2_b_j
1238                        .iter()
1239                        .map(|tensor_rank_2_b_jk| tensor_rank_2_a_i * tensor_rank_2_b_jk)
1240                        .collect()
1241                })
1242                .collect()
1243        })
1244        .collect()
1245}
1246
1247fn canonical_rank_1_times_rank_4<const D: usize>(
1248    tensor_rank_1: &TensorRank1<D, Reference, Dimensionless>,
1249    tensor_rank_4: &TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless>,
1250) -> TensorRank3<D, Reference, Reference, Reference, Dimensionless> {
1251    tensor_rank_1
1252        .iter()
1253        .zip(tensor_rank_4.iter())
1254        .map(|(&tensor_rank_1_m, tensor_rank_4_m)| tensor_rank_4_m * tensor_rank_1_m.value())
1255        .sum()
1256}
1257
1258fn canonical_rank_2_times_rank_4<const D: usize>(
1259    tensor_rank_2: &TensorRank2<D, Reference, Reference, Dimensionless>,
1260    tensor_rank_4: &TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless>,
1261) -> TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless> {
1262    tensor_rank_2
1263        .iter()
1264        .map(|tensor_rank_2_i| tensor_rank_2_i * tensor_rank_4)
1265        .collect()
1266}
1267
1268fn canonical_contract_second_with_first<const D: usize>(
1269    tensor_rank_4: &TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless>,
1270    tensor_rank_2: &TensorRank2<D, Reference, Reference, Dimensionless>,
1271) -> TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless> {
1272    let mut output = TensorRank4::zero();
1273    tensor_rank_4
1274        .iter()
1275        .zip(output.iter_mut())
1276        .for_each(|(tensor_rank_4_i, output_i)| {
1277            tensor_rank_4_i.iter().zip(tensor_rank_2.iter()).for_each(
1278                |(tensor_rank_4_is, tensor_rank_2_s)| {
1279                    tensor_rank_2_s.iter().zip(output_i.iter_mut()).for_each(
1280                        |(tensor_rank_2_sj, output_ij)| {
1281                            output_ij.iter_mut().zip(tensor_rank_4_is.iter()).for_each(
1282                                |(output_ijk, tensor_rank_4_isk)| {
1283                                    output_ijk
1284                                        .iter_mut()
1285                                        .zip(tensor_rank_4_isk.iter())
1286                                        .for_each(|(output_ijkl, tensor_rank_4_iskl)| {
1287                                            *output_ijkl += tensor_rank_4_iskl * tensor_rank_2_sj
1288                                        })
1289                                },
1290                            )
1291                        },
1292                    )
1293                },
1294            )
1295        });
1296    output
1297}
1298
1299fn canonical_contract_third_with_first<const D: usize>(
1300    tensor_rank_4: &TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless>,
1301    tensor_rank_2: &TensorRank2<D, Reference, Reference, Dimensionless>,
1302) -> TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless> {
1303    let mut output = TensorRank4::zero();
1304    tensor_rank_4
1305        .iter()
1306        .zip(output.iter_mut())
1307        .for_each(|(tensor_rank_4_i, output_i)| {
1308            tensor_rank_4_i.iter().zip(output_i.iter_mut()).for_each(
1309                |(tensor_rank_4_ij, output_ij)| {
1310                    tensor_rank_4_ij.iter().zip(tensor_rank_2.iter()).for_each(
1311                        |(tensor_rank_4_ijm, tensor_rank_2_m)| {
1312                            tensor_rank_2_m.iter().zip(output_ij.iter_mut()).for_each(
1313                                |(tensor_rank_2_mk, output_ijk)| {
1314                                    output_ijk
1315                                        .iter_mut()
1316                                        .zip(tensor_rank_4_ijm.iter())
1317                                        .for_each(|(output_ijkl, tensor_rank_4_ijml)| {
1318                                            *output_ijkl += tensor_rank_2_mk * tensor_rank_4_ijml
1319                                        })
1320                                },
1321                            )
1322                        },
1323                    )
1324                },
1325            )
1326        });
1327    output
1328}
1329
1330fn canonical_transform_first<const D: usize>(
1331    tensor_rank_4: &TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless>,
1332    tensor_rank_2: &TensorRank2<D, Reference, Reference, Dimensionless>,
1333) -> TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless> {
1334    let mut output = TensorRank4::zero();
1335    tensor_rank_4.iter().zip(tensor_rank_2.iter()).for_each(
1336        |(tensor_rank_4_m, tensor_rank_2_m)| {
1337            tensor_rank_2_m.iter().zip(output.iter_mut()).for_each(
1338                |(tensor_rank_2_mi, output_i)| {
1339                    output_i.iter_mut().zip(tensor_rank_4_m.iter()).for_each(
1340                        |(output_ij, tensor_rank_4_mj)| {
1341                            output_ij.iter_mut().zip(tensor_rank_4_mj.iter()).for_each(
1342                                |(output_ijk, tensor_rank_4_mjk)| {
1343                                    output_ijk
1344                                        .iter_mut()
1345                                        .zip(tensor_rank_4_mjk.iter())
1346                                        .for_each(|(output_ijkl, tensor_rank_4_mjkl)| {
1347                                            *output_ijkl += tensor_rank_4_mjkl * tensor_rank_2_mi
1348                                        })
1349                                },
1350                            )
1351                        },
1352                    )
1353                },
1354            )
1355        },
1356    );
1357    output
1358}
1359
1360fn canonical_transform_fourth<const D: usize>(
1361    tensor_rank_4: &TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless>,
1362    tensor_rank_2: &TensorRank2<D, Reference, Reference, Dimensionless>,
1363) -> TensorRank4<D, Reference, Reference, Reference, Reference, Dimensionless> {
1364    let mut output = TensorRank4::zero();
1365    tensor_rank_4
1366        .iter()
1367        .zip(output.iter_mut())
1368        .for_each(|(tensor_rank_4_i, output_i)| {
1369            tensor_rank_4_i.iter().zip(output_i.iter_mut()).for_each(
1370                |(tensor_rank_4_ij, output_ij)| {
1371                    tensor_rank_4_ij.iter().zip(output_ij.iter_mut()).for_each(
1372                        |(tensor_rank_4_ijk, output_ijk)| {
1373                            tensor_rank_4_ijk.iter().zip(tensor_rank_2.iter()).for_each(
1374                                |(tensor_rank_4_ijkm, tensor_rank_2_m)| {
1375                                    output_ijk.iter_mut().zip(tensor_rank_2_m.iter()).for_each(
1376                                        |(output_ijkl, tensor_rank_2_ml)| {
1377                                            *output_ijkl += tensor_rank_4_ijkm * tensor_rank_2_ml
1378                                        },
1379                                    )
1380                                },
1381                            )
1382                        },
1383                    )
1384                },
1385            )
1386        });
1387    output
1388}
1389
1390impl<const D: usize, I, J, K, L, U, T> Differentiate<T> for TensorRank4<D, I, J, K, L, U>
1391where
1392    U: UnitDiv<T>,
1393{
1394    type Derivative = TensorRank4<D, I, J, K, L, <U as UnitDiv<T>>::Output>;
1395}
conspire/math/tensor/rank_4/mod.rs

conspire/math/tensor/rank_4/
mod.rs
