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

conspire/math/tensor/rank_4/
mod.rs
