conspire/math/tensor/rank_4/
mod.rs

1#[cfg(test)]
2mod test;
3
4#[cfg(test)]
5use super::test::ErrorTensor;
6
7use std::{
8    array::from_fn,
9    fmt::{self, Display, Formatter},
10    iter::Sum,
11    mem::transmute,
12    ops::{Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Sub, SubAssign},
13};
14
15use super::{
16    Hessian, Rank2, SquareMatrix, Tensor, TensorArray, Vector,
17    rank_0::TensorRank0,
18    rank_1::TensorRank1,
19    rank_2::TensorRank2,
20    rank_3::{TensorRank3, get_identity_1010_parts},
21};
22
23pub mod list;
24pub mod vec;
25
26/// A *d*-dimensional tensor of rank 4.
27///
28/// `D` is the dimension, `I`, `J`, `K`, `L` are the configurations.
29#[repr(transparent)]
30#[derive(Clone, Debug, PartialEq)]
31pub struct TensorRank4<
32    const D: usize,
33    const I: usize,
34    const J: usize,
35    const K: usize,
36    const L: usize,
37>([TensorRank3<D, J, K, L>; D]);
38
39impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> Default
40    for TensorRank4<D, I, J, K, L>
41{
42    fn default() -> Self {
43        Self::zero()
44    }
45}
46
47impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
48    From<[[[[TensorRank0; D]; D]; D]; D]> for TensorRank4<D, I, J, K, L>
49{
50    fn from(array: [[[[TensorRank0; D]; D]; D]; D]) -> Self {
51        array.into_iter().map(|entry| entry.into()).collect()
52    }
53}
54
55pub const IDENTITY_1010: TensorRank4<3, 1, 0, 1, 0> = TensorRank4(get_identity_1010_parts());
56
57impl<const I: usize, const J: usize, const K: usize> From<TensorRank4<3, I, J, K, 0>>
58    for TensorRank4<3, I, J, K, 2>
59{
60    fn from(tensor_rank_4: TensorRank4<3, I, J, K, 0>) -> Self {
61        unsafe {
62            transmute::<TensorRank4<3, I, J, K, 0>, TensorRank4<3, I, J, K, 2>>(tensor_rank_4)
63        }
64    }
65}
66
67impl<const J: usize, const L: usize> From<TensorRank4<3, 1, J, 1, L>>
68    for TensorRank4<3, 2, J, 2, L>
69{
70    fn from(tensor_rank_4: TensorRank4<3, 1, J, 1, L>) -> Self {
71        unsafe {
72            transmute::<TensorRank4<3, 1, J, 1, L>, TensorRank4<3, 2, J, 2, L>>(tensor_rank_4)
73        }
74    }
75}
76
77impl<const I: usize, const K: usize> From<TensorRank4<3, I, 0, K, 0>>
78    for TensorRank4<3, I, 2, K, 2>
79{
80    fn from(tensor_rank_4: TensorRank4<3, I, 0, K, 0>) -> Self {
81        unsafe {
82            transmute::<TensorRank4<3, I, 0, K, 0>, TensorRank4<3, I, 2, K, 2>>(tensor_rank_4)
83        }
84    }
85}
86
87impl<const K: usize> From<TensorRank4<3, 0, 0, K, 0>> for TensorRank4<3, 2, 2, K, 2> {
88    fn from(tensor_rank_4: TensorRank4<3, 0, 0, K, 0>) -> Self {
89        unsafe {
90            transmute::<TensorRank4<3, 0, 0, K, 0>, TensorRank4<3, 2, 2, K, 2>>(tensor_rank_4)
91        }
92    }
93}
94
95impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
96    From<Vec<Vec<Vec<Vec<TensorRank0>>>>> for TensorRank4<D, I, J, K, L>
97{
98    fn from(vec_rank_4: Vec<Vec<Vec<Vec<TensorRank0>>>>) -> Self {
99        vec_rank_4
100            .into_iter()
101            .map(|vec_rank_3| {
102                vec_rank_3
103                    .into_iter()
104                    .map(|vec_rank_2| {
105                        vec_rank_2
106                            .into_iter()
107                            .map(|vec_rank_1| vec_rank_1.into_iter().collect())
108                            .collect()
109                    })
110                    .collect()
111            })
112            .collect()
113    }
114}
115
116impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
117    From<TensorRank4<D, I, J, K, L>> for Vec<Vec<Vec<Vec<TensorRank0>>>>
118{
119    fn from(tensor_rank_4: TensorRank4<D, I, J, K, L>) -> Self {
120        tensor_rank_4
121            .iter()
122            .map(|tensor_rank_3| {
123                tensor_rank_3
124                    .iter()
125                    .map(|tensor_rank_2| {
126                        tensor_rank_2
127                            .iter()
128                            .map(|tensor_rank_1| tensor_rank_1.iter().copied().collect())
129                            .collect()
130                    })
131                    .collect()
132            })
133            .collect()
134    }
135}
136
137impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
138    From<TensorRank4<D, I, J, K, L>> for Vec<TensorRank0>
139{
140    fn from(tensor_rank_4: TensorRank4<D, I, J, K, L>) -> Self {
141        tensor_rank_4
142            .iter()
143            .flat_map(|tensor_rank_3| {
144                tensor_rank_3.iter().flat_map(|tensor_rank_2| {
145                    tensor_rank_2
146                        .iter()
147                        .flat_map(|tensor_rank_1| tensor_rank_1.iter().copied())
148                })
149            })
150            .collect()
151    }
152}
153
154impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
155    From<TensorRank4<D, I, J, K, L>> for Vector
156{
157    fn from(tensor_rank_4: TensorRank4<D, I, J, K, L>) -> Self {
158        tensor_rank_4
159            .iter()
160            .flat_map(|tensor_rank_3| {
161                tensor_rank_3.iter().flat_map(|tensor_rank_2| {
162                    tensor_rank_2
163                        .iter()
164                        .flat_map(|tensor_rank_1| tensor_rank_1.iter().copied())
165                })
166            })
167            .collect()
168    }
169}
170
171impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> Display
172    for TensorRank4<D, I, J, K, L>
173{
174    fn fmt(&self, f: &mut Formatter) -> fmt::Result {
175        write!(f, "[")?;
176        self.iter()
177            .enumerate()
178            .try_for_each(|(i, entry)| write!(f, "{entry},\n\x1B[u\x1B[{}B\x1B[2D", i + 1))?;
179        write!(f, "\x1B[u\x1B[1A\x1B[{}C]", 16 * D + 2)
180    }
181}
182
183#[cfg(test)]
184impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> ErrorTensor
185    for TensorRank4<D, I, J, K, L>
186{
187    fn error_fd(&self, comparator: &Self, epsilon: TensorRank0) -> Option<(bool, usize)> {
188        let error_count = self
189            .iter()
190            .zip(comparator.iter())
191            .map(|(self_i, comparator_i)| {
192                self_i
193                    .iter()
194                    .zip(comparator_i.iter())
195                    .map(|(self_ij, comparator_ij)| {
196                        self_ij
197                            .iter()
198                            .zip(comparator_ij.iter())
199                            .map(|(self_ijk, comparator_ijk)| {
200                                self_ijk
201                                    .iter()
202                                    .zip(comparator_ijk.iter())
203                                    .filter(|&(&self_ijkl, &comparator_ijkl)| {
204                                        (self_ijkl / comparator_ijkl - 1.0).abs() >= epsilon
205                                            && (self_ijkl.abs() >= epsilon
206                                                || comparator_ijkl.abs() >= epsilon)
207                                    })
208                                    .count()
209                            })
210                            .sum::<usize>()
211                    })
212                    .sum::<usize>()
213            })
214            .sum();
215        if error_count > 0 {
216            let auxiliary = self
217                .iter()
218                .zip(comparator.iter())
219                .map(|(self_i, comparator_i)| {
220                    self_i
221                        .iter()
222                        .zip(comparator_i.iter())
223                        .map(|(self_ij, comparator_ij)| {
224                            self_ij
225                                .iter()
226                                .zip(comparator_ij.iter())
227                                .map(|(self_ijk, comparator_ijk)| {
228                                    self_ijk
229                                        .iter()
230                                        .zip(comparator_ijk.iter())
231                                        .filter(|&(&self_ijkl, &comparator_ijkl)| {
232                                            (self_ijkl / comparator_ijkl - 1.0).abs() >= epsilon
233                                                && (self_ijkl - comparator_ijkl).abs() >= epsilon
234                                                && (self_ijkl.abs() >= epsilon
235                                                    || comparator_ijkl.abs() >= epsilon)
236                                        })
237                                        .count()
238                                })
239                                .sum::<usize>()
240                        })
241                        .sum::<usize>()
242                })
243                .sum::<usize>()
244                > 0;
245            Some((auxiliary, error_count))
246        } else {
247            None
248        }
249    }
250}
251
252impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
253    TensorRank4<D, I, J, K, L>
254{
255    pub fn dyad_ij_kl(
256        tensor_rank_2_a: &TensorRank2<D, I, J>,
257        tensor_rank_2_b: &TensorRank2<D, K, L>,
258    ) -> Self {
259        tensor_rank_2_a
260            .iter()
261            .map(|tensor_rank_2_a_i| {
262                tensor_rank_2_a_i
263                    .iter()
264                    .map(|tensor_rank_2_a_ij| tensor_rank_2_b * tensor_rank_2_a_ij)
265                    .collect()
266            })
267            .collect()
268    }
269    pub fn dyad_ik_jl(
270        tensor_rank_2_a: &TensorRank2<D, I, K>,
271        tensor_rank_2_b: &TensorRank2<D, J, L>,
272    ) -> Self {
273        tensor_rank_2_a
274            .iter()
275            .map(|tensor_rank_2_a_i| {
276                tensor_rank_2_b
277                    .iter()
278                    .map(|tensor_rank_2_b_j| {
279                        tensor_rank_2_a_i
280                            .iter()
281                            .map(|tensor_rank_2_a_ik| tensor_rank_2_b_j * tensor_rank_2_a_ik)
282                            .collect()
283                    })
284                    .collect()
285            })
286            .collect()
287    }
288    pub fn dyad_il_jk(
289        tensor_rank_2_a: &TensorRank2<D, I, L>,
290        tensor_rank_2_b: &TensorRank2<D, J, K>,
291    ) -> Self {
292        tensor_rank_2_a
293            .iter()
294            .map(|tensor_rank_2_a_i| {
295                tensor_rank_2_b
296                    .iter()
297                    .map(|tensor_rank_2_b_j| {
298                        tensor_rank_2_b_j
299                            .iter()
300                            .map(|tensor_rank_2_b_jk| tensor_rank_2_a_i * tensor_rank_2_b_jk)
301                            .collect()
302                    })
303                    .collect()
304            })
305            .collect()
306    }
307    pub fn dyad_il_kj(
308        tensor_rank_2_a: &TensorRank2<D, I, L>,
309        tensor_rank_2_b: &TensorRank2<D, K, J>,
310    ) -> Self {
311        Self::dyad_il_jk(tensor_rank_2_a, &(tensor_rank_2_b.transpose()))
312    }
313}
314
315impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> Hessian
316    for TensorRank4<D, I, J, K, L>
317{
318    fn fill_into(self, square_matrix: &mut SquareMatrix) {
319        self.into_iter().enumerate().for_each(|(i, self_i)| {
320            self_i.into_iter().enumerate().for_each(|(j, self_ij)| {
321                self_ij.into_iter().enumerate().for_each(|(k, self_ijk)| {
322                    self_ijk
323                        .into_iter()
324                        .enumerate()
325                        .for_each(|(l, self_ijkl)| square_matrix[D * i + j][D * k + l] = self_ijkl)
326                })
327            })
328        })
329    }
330}
331
332impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> Tensor
333    for TensorRank4<D, I, J, K, L>
334{
335    type Item = TensorRank3<D, J, K, L>;
336    fn iter(&self) -> impl Iterator<Item = &Self::Item> {
337        self.0.iter()
338    }
339    fn iter_mut(&mut self) -> impl Iterator<Item = &mut Self::Item> {
340        self.0.iter_mut()
341    }
342    fn len(&self) -> usize {
343        D
344    }
345    fn size(&self) -> usize {
346        D * D * D * D
347    }
348}
349
350impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> IntoIterator
351    for TensorRank4<D, I, J, K, L>
352{
353    type Item = TensorRank3<D, J, K, L>;
354    type IntoIter = std::array::IntoIter<Self::Item, D>;
355    fn into_iter(self) -> Self::IntoIter {
356        self.0.into_iter()
357    }
358}
359
360impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> TensorArray
361    for TensorRank4<D, I, J, K, L>
362{
363    type Array = [[[[TensorRank0; D]; D]; D]; D];
364    type Item = TensorRank3<D, J, K, L>;
365    fn as_array(&self) -> Self::Array {
366        let mut array = [[[[0.0; D]; D]; D]; D];
367        array
368            .iter_mut()
369            .zip(self.iter())
370            .for_each(|(entry_rank_3, tensor_rank_3)| *entry_rank_3 = tensor_rank_3.as_array());
371        array
372    }
373    fn identity() -> Self {
374        Self::dyad_ij_kl(&TensorRank2::identity(), &TensorRank2::identity())
375    }
376    fn zero() -> Self {
377        Self(from_fn(|_| Self::Item::zero()))
378    }
379}
380
381impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
382    FromIterator<TensorRank3<D, J, K, L>> for TensorRank4<D, I, J, K, L>
383{
384    fn from_iter<Ii: IntoIterator<Item = TensorRank3<D, J, K, L>>>(into_iterator: Ii) -> Self {
385        let mut tensor_rank_4 = Self::zero();
386        tensor_rank_4
387            .iter_mut()
388            .zip(into_iterator)
389            .for_each(|(tensor_rank_4_i, value_i)| *tensor_rank_4_i = value_i);
390        tensor_rank_4
391    }
392}
393
394impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> Index<usize>
395    for TensorRank4<D, I, J, K, L>
396{
397    type Output = TensorRank3<D, J, K, L>;
398    fn index(&self, index: usize) -> &Self::Output {
399        &self.0[index]
400    }
401}
402
403impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> IndexMut<usize>
404    for TensorRank4<D, I, J, K, L>
405{
406    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
407        &mut self.0[index]
408    }
409}
410
411impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> Sum
412    for TensorRank4<D, I, J, K, L>
413{
414    fn sum<Ii>(iter: Ii) -> Self
415    where
416        Ii: Iterator<Item = Self>,
417    {
418        iter.reduce(|mut acc, item| {
419            acc += item;
420            acc
421        })
422        .unwrap_or_else(Self::default)
423    }
424}
425
426impl<'a, const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
427    Sum<&'a Self> for TensorRank4<D, I, J, K, L>
428{
429    fn sum<Ii>(iter: Ii) -> Self
430    where
431        Ii: Iterator<Item = &'a Self>,
432    {
433        iter.fold(Self::default(), |mut acc, item| {
434            acc += item;
435            acc
436        })
437    }
438}
439
440pub trait ContractAllIndicesWithFirstIndicesOf<TIM, TJN, TKO, TLP> {
441    type Output;
442    fn contract_all_indices_with_first_indices_of(
443        self,
444        tensor_rank_2_a: TIM,
445        tensor_rank_2_b: TJN,
446        tensor_rank_2_c: TKO,
447        tensor_rank_2_d: TLP,
448    ) -> Self::Output;
449}
450
451impl<
452    const D: usize,
453    const I: usize,
454    const J: usize,
455    const K: usize,
456    const L: usize,
457    const M: usize,
458    const N: usize,
459    const O: usize,
460    const P: usize,
461>
462    ContractAllIndicesWithFirstIndicesOf<
463        &TensorRank2<D, I, M>,
464        &TensorRank2<D, J, N>,
465        &TensorRank2<D, K, O>,
466        &TensorRank2<D, L, P>,
467    > for TensorRank4<D, I, J, K, L>
468{
469    type Output = TensorRank4<D, M, N, O, P>;
470    fn contract_all_indices_with_first_indices_of(
471        self,
472        tensor_rank_2_a: &TensorRank2<D, I, M>,
473        tensor_rank_2_b: &TensorRank2<D, J, N>,
474        tensor_rank_2_c: &TensorRank2<D, K, O>,
475        tensor_rank_2_d: &TensorRank2<D, L, P>,
476    ) -> Self::Output {
477        let mut output = TensorRank4::zero();
478        self.into_iter().zip(tensor_rank_2_a.iter()).for_each(|(self_m, tensor_rank_2_a_m)|
479            self_m.into_iter().zip(tensor_rank_2_b.iter()).for_each(|(self_mn, tensor_rank_2_b_n)|
480                self_mn.into_iter().zip(tensor_rank_2_c.iter()).for_each(|(self_mno, tensor_rank_2_c_o)|
481                    self_mno.into_iter().zip(tensor_rank_2_d.iter()).for_each(|(self_mnop, tensor_rank_2_d_p)|
482                        output.iter_mut().zip(tensor_rank_2_a_m.iter()).for_each(|(output_i, tensor_rank_2_a_mi)|
483                            output_i.iter_mut().zip(tensor_rank_2_b_n.iter()).for_each(|(output_ij, tensor_rank_2_b_nj)|
484                                output_ij.iter_mut().zip(tensor_rank_2_c_o.iter()).for_each(|(output_ijk, tensor_rank_2_c_ok)|
485                                    output_ijk.iter_mut().zip(tensor_rank_2_d_p.iter()).for_each(|(output_ijkl, tensor_rank_2_dp)|
486                                        *output_ijkl += self_mnop * tensor_rank_2_a_mi * tensor_rank_2_b_nj * tensor_rank_2_c_ok * tensor_rank_2_dp
487                                    )
488                                )
489                            )
490                        )
491                    )
492                )
493            )
494        );
495        output
496    }
497}
498
499pub trait ContractFirstThirdFourthIndicesWithFirstIndicesOf<TIM, TKO, TLP> {
500    type Output;
501    fn contract_first_third_fourth_indices_with_first_indices_of(
502        self,
503        tensor_rank_2_a: TIM,
504        tensor_rank_2_c: TKO,
505        tensor_rank_2_d: TLP,
506    ) -> Self::Output;
507}
508
509impl<
510    const D: usize,
511    const I: usize,
512    const J: usize,
513    const K: usize,
514    const L: usize,
515    const M: usize,
516    const O: usize,
517    const P: usize,
518>
519    ContractFirstThirdFourthIndicesWithFirstIndicesOf<
520        &TensorRank2<D, I, M>,
521        &TensorRank2<D, K, O>,
522        &TensorRank2<D, L, P>,
523    > for TensorRank4<D, I, J, K, L>
524{
525    type Output = TensorRank4<D, M, J, O, P>;
526    fn contract_first_third_fourth_indices_with_first_indices_of(
527        self,
528        tensor_rank_2_a: &TensorRank2<D, I, M>,
529        tensor_rank_2_b: &TensorRank2<D, K, O>,
530        tensor_rank_2_c: &TensorRank2<D, L, P>,
531    ) -> Self::Output {
532        let mut output = TensorRank4::zero();
533        self.into_iter().zip(tensor_rank_2_a.iter()).for_each(|(self_q, tensor_rank_2_a_q)|
534            output.iter_mut().zip(tensor_rank_2_a_q.iter()).for_each(|(output_i, tensor_rank_2_a_qi)|
535                output_i.iter_mut().zip(self_q.iter()).for_each(|(output_ij, self_qj)|
536                    self_qj.iter().zip(tensor_rank_2_b.iter()).for_each(|(self_qjm, tensor_rank_2_b_m)|
537                        self_qjm.iter().zip(tensor_rank_2_c.iter()).for_each(|(self_qjmn, tensor_rank_2_c_n)|
538                            output_ij.iter_mut().zip(tensor_rank_2_b_m.iter()).for_each(|(output_ijk, tensor_rank_2_b_mk)|
539                                output_ijk.iter_mut().zip(tensor_rank_2_c_n.iter()).for_each(|(output_ijkl, tensor_rank_2_c_nl)|
540                                    *output_ijkl += self_qjmn*tensor_rank_2_a_qi*tensor_rank_2_b_mk*tensor_rank_2_c_nl
541                                )
542                            )
543                        )
544                    )
545                )
546            )
547        );
548        output
549    }
550}
551
552pub trait ContractSecondIndexWithFirstIndexOf<TJN> {
553    type Output;
554    fn contract_second_index_with_first_index_of(self, tensor_rank_2: TJN) -> Self::Output;
555}
556
557impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize, const N: usize>
558    ContractSecondIndexWithFirstIndexOf<&TensorRank2<D, J, N>> for TensorRank4<D, I, J, K, L>
559{
560    type Output = TensorRank4<D, I, N, K, L>;
561    fn contract_second_index_with_first_index_of(
562        self,
563        tensor_rank_2: &TensorRank2<D, J, N>,
564    ) -> Self::Output {
565        let mut output = TensorRank4::zero();
566        output.iter_mut().zip(self).for_each(|(output_i, self_i)| {
567            self_i
568                .into_iter()
569                .zip(tensor_rank_2.iter())
570                .for_each(|(self_is, tensor_rank_2_s)| {
571                    output_i.iter_mut().zip(tensor_rank_2_s.iter()).for_each(
572                        |(output_ij, tensor_rank_2_sj)| *output_ij += &self_is * tensor_rank_2_sj,
573                    )
574                })
575        });
576        output
577    }
578}
579
580pub trait ContractSecondFourthIndicesWithFirstIndicesOf<TJ, TL> {
581    type Output;
582    fn contract_second_fourth_indices_with_first_indices_of(
583        &self,
584        object_a: TJ,
585        object_b: TL,
586    ) -> Self::Output;
587}
588
589impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
590    ContractSecondFourthIndicesWithFirstIndicesOf<&TensorRank1<D, J>, &TensorRank1<D, L>>
591    for TensorRank4<D, I, J, K, L>
592{
593    type Output = TensorRank2<D, I, K>;
594    fn contract_second_fourth_indices_with_first_indices_of(
595        &self,
596        tensor_rank_1_a: &TensorRank1<D, J>,
597        tensor_rank_1_b: &TensorRank1<D, L>,
598    ) -> Self::Output {
599        self.iter()
600            .map(|self_i| {
601                self_i
602                    .iter()
603                    .zip(tensor_rank_1_a.iter())
604                    .map(|(self_ij, tensor_rank_1_a_j)| {
605                        self_ij
606                            .iter()
607                            .map(|self_ijk| self_ijk * (tensor_rank_1_b * tensor_rank_1_a_j))
608                            .collect::<TensorRank1<D, K>>()
609                    })
610                    .sum()
611            })
612            .collect()
613    }
614}
615
616pub trait ContractThirdIndexWithFirstIndexOf<TKL> {
617    type Output;
618    fn contract_third_index_with_first_index_of(&self, tensor: TKL) -> Self::Output;
619}
620
621impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize, const M: usize>
622    ContractThirdIndexWithFirstIndexOf<&TensorRank2<D, M, K>> for TensorRank4<D, I, J, M, L>
623{
624    type Output = TensorRank4<D, I, J, K, L>;
625    fn contract_third_index_with_first_index_of(
626        &self,
627        tensor_rank_2: &TensorRank2<D, M, K>,
628    ) -> Self::Output {
629        self.iter()
630            .map(|self_i| {
631                self_i
632                    .iter()
633                    .map(|self_ij| {
634                        self_ij
635                            .iter()
636                            .zip(tensor_rank_2.iter())
637                            .map(|(self_ijm, tensor_rank_2_m)| {
638                                TensorRank2::<D, K, L>::from((tensor_rank_2_m, self_ijm))
639                            })
640                            .sum()
641                    })
642                    .collect()
643            })
644            .collect()
645    }
646}
647
648pub trait ContractThirdFourthIndicesWithFirstSecondIndicesOf<TKL> {
649    type Output;
650    fn contract_third_fourth_indices_with_first_second_indices_of(
651        self,
652        tensor: TKL,
653    ) -> Self::Output;
654}
655
656impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
657    ContractThirdFourthIndicesWithFirstSecondIndicesOf<&TensorRank2<D, K, L>>
658    for TensorRank4<D, I, J, K, L>
659{
660    type Output = TensorRank2<D, I, J>;
661    fn contract_third_fourth_indices_with_first_second_indices_of(
662        self,
663        tensor_rank_2: &TensorRank2<D, K, L>,
664    ) -> Self::Output {
665        self.into_iter()
666            .map(|self_i| {
667                self_i
668                    .into_iter()
669                    .map(|self_ij| self_ij.full_contraction(tensor_rank_2))
670                    .collect()
671            })
672            .collect()
673    }
674}
675
676impl<
677    const D: usize,
678    const I: usize,
679    const J: usize,
680    const K: usize,
681    const L: usize,
682    const M: usize,
683    const N: usize,
684> ContractThirdFourthIndicesWithFirstSecondIndicesOf<&TensorRank4<D, K, L, M, N>>
685    for TensorRank4<D, I, J, K, L>
686{
687    type Output = TensorRank4<D, I, J, M, N>;
688    fn contract_third_fourth_indices_with_first_second_indices_of(
689        self,
690        tensor: &TensorRank4<D, K, L, M, N>,
691    ) -> Self::Output {
692        self.into_iter()
693            .map(|self_i| {
694                self_i
695                    .into_iter()
696                    .map(|self_ij| {
697                        self_ij
698                            .into_iter()
699                            .zip(tensor.iter())
700                            .map(|(self_ijk, tensor_k)| {
701                                self_ijk
702                                    .into_iter()
703                                    .zip(tensor_k.iter())
704                                    .map(|(self_ijkl, tensor_kl)| tensor_kl * self_ijkl)
705                                    .sum::<TensorRank2<D, M, N>>()
706                            })
707                            .sum()
708                    })
709                    .collect()
710            })
711            .collect()
712    }
713}
714
715pub trait ContractFirstSecondIndicesWithSecondIndicesOf<TI, TJ> {
716    type Output;
717    fn contract_first_second_indices_with_second_indices_of(
718        self,
719        object_a: TI,
720        object_b: TJ,
721    ) -> Self::Output;
722}
723
724impl<
725    const D: usize,
726    const I: usize,
727    const J: usize,
728    const K: usize,
729    const L: usize,
730    const M: usize,
731    const N: usize,
732> ContractFirstSecondIndicesWithSecondIndicesOf<&TensorRank2<D, I, M>, &TensorRank2<D, J, N>>
733    for TensorRank4<D, M, N, K, L>
734{
735    type Output = TensorRank4<D, I, J, K, L>;
736    fn contract_first_second_indices_with_second_indices_of(
737        self,
738        tensor_rank_2_a: &TensorRank2<D, I, M>,
739        tensor_rank_2_b: &TensorRank2<D, J, N>,
740    ) -> Self::Output {
741        let mut output = TensorRank4::zero();
742        output
743            .iter_mut()
744            .zip(tensor_rank_2_a.iter())
745            .for_each(|(output_i, tensor_rank_2_a_i)| {
746                output_i.iter_mut().zip(tensor_rank_2_b.iter()).for_each(
747                    |(output_ij, tensor_rank_2_b_j)| {
748                        self.iter().zip(tensor_rank_2_a_i.iter()).for_each(
749                            |(self_m, tensor_rank_2_a_im)| {
750                                self_m.iter().zip(tensor_rank_2_b_j.iter()).for_each(
751                                    |(self_mn, tensor_rank_2_b_jn)| {
752                                        *output_ij +=
753                                            self_mn * tensor_rank_2_a_im * tensor_rank_2_b_jn
754                                    },
755                                )
756                            },
757                        )
758                    },
759                )
760            });
761        output
762    }
763}
764
765impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
766    Div<TensorRank0> for TensorRank4<D, I, J, K, L>
767{
768    type Output = Self;
769    fn div(mut self, tensor_rank_0: TensorRank0) -> Self::Output {
770        self /= &tensor_rank_0;
771        self
772    }
773}
774
775impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
776    Div<TensorRank0> for &TensorRank4<D, I, J, K, L>
777{
778    type Output = TensorRank4<D, I, J, K, L>;
779    fn div(self, tensor_rank_0: TensorRank0) -> Self::Output {
780        self.iter().map(|self_i| self_i / tensor_rank_0).collect()
781    }
782}
783
784impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
785    Div<&TensorRank0> for TensorRank4<D, I, J, K, L>
786{
787    type Output = Self;
788    fn div(mut self, tensor_rank_0: &TensorRank0) -> Self::Output {
789        self /= tensor_rank_0;
790        self
791    }
792}
793
794impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
795    DivAssign<TensorRank0> for TensorRank4<D, I, J, K, L>
796{
797    fn div_assign(&mut self, tensor_rank_0: TensorRank0) {
798        self.iter_mut().for_each(|self_i| *self_i /= &tensor_rank_0);
799    }
800}
801
802impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
803    DivAssign<&TensorRank0> for TensorRank4<D, I, J, K, L>
804{
805    fn div_assign(&mut self, tensor_rank_0: &TensorRank0) {
806        self.iter_mut().for_each(|self_i| *self_i /= tensor_rank_0);
807    }
808}
809
810impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
811    Mul<TensorRank0> for TensorRank4<D, I, J, K, L>
812{
813    type Output = Self;
814    fn mul(mut self, tensor_rank_0: TensorRank0) -> Self::Output {
815        self *= &tensor_rank_0;
816        self
817    }
818}
819
820impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
821    Mul<&TensorRank0> for TensorRank4<D, I, J, K, L>
822{
823    type Output = Self;
824    fn mul(mut self, tensor_rank_0: &TensorRank0) -> Self::Output {
825        self *= tensor_rank_0;
826        self
827    }
828}
829
830impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
831    MulAssign<TensorRank0> for TensorRank4<D, I, J, K, L>
832{
833    fn mul_assign(&mut self, tensor_rank_0: TensorRank0) {
834        self.iter_mut().for_each(|self_i| *self_i *= &tensor_rank_0);
835    }
836}
837
838impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
839    MulAssign<&TensorRank0> for TensorRank4<D, I, J, K, L>
840{
841    fn mul_assign(&mut self, tensor_rank_0: &TensorRank0) {
842        self.iter_mut().for_each(|self_i| *self_i *= tensor_rank_0);
843    }
844}
845
846impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize, const M: usize>
847    Mul<TensorRank2<D, L, M>> for TensorRank4<D, I, J, K, L>
848{
849    type Output = TensorRank4<D, I, J, K, M>;
850    fn mul(self, tensor_rank_2: TensorRank2<D, L, M>) -> Self::Output {
851        self.into_iter()
852            .map(|self_i| {
853                self_i
854                    .into_iter()
855                    .map(|self_ij| self_ij * &tensor_rank_2)
856                    .collect()
857            })
858            .collect()
859    }
860}
861
862impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize, const M: usize>
863    Mul<&TensorRank2<D, L, M>> for TensorRank4<D, I, J, K, L>
864{
865    type Output = TensorRank4<D, I, J, K, M>;
866    fn mul(self, tensor_rank_2: &TensorRank2<D, L, M>) -> Self::Output {
867        self.into_iter()
868            .map(|self_i| {
869                self_i
870                    .into_iter()
871                    .map(|self_ij| self_ij * tensor_rank_2)
872                    .collect()
873            })
874            .collect()
875    }
876}
877
878impl<const D: usize, const J: usize, const K: usize, const L: usize, const M: usize>
879    Mul<TensorRank4<D, M, J, K, L>> for TensorRank1<D, M>
880{
881    type Output = TensorRank3<D, J, K, L>;
882    fn mul(self, tensor_rank_4: TensorRank4<D, M, J, K, L>) -> Self::Output {
883        self.into_iter()
884            .zip(tensor_rank_4)
885            .map(|(self_m, tensor_rank_4_m)| tensor_rank_4_m * self_m)
886            .sum()
887    }
888}
889
890impl<const D: usize, const J: usize, const K: usize, const L: usize, const M: usize>
891    Mul<&TensorRank4<D, M, J, K, L>> for TensorRank1<D, M>
892{
893    type Output = TensorRank3<D, J, K, L>;
894    fn mul(self, tensor_rank_4: &TensorRank4<D, M, J, K, L>) -> Self::Output {
895        self.into_iter()
896            .zip(tensor_rank_4.iter())
897            .map(|(self_m, tensor_rank_4_m)| tensor_rank_4_m * self_m)
898            .sum()
899    }
900}
901
902impl<const D: usize, const J: usize, const K: usize, const L: usize, const M: usize>
903    Mul<TensorRank4<D, M, J, K, L>> for &TensorRank1<D, M>
904{
905    type Output = TensorRank3<D, J, K, L>;
906    fn mul(self, tensor_rank_4: TensorRank4<D, M, J, K, L>) -> Self::Output {
907        self.iter()
908            .zip(tensor_rank_4)
909            .map(|(self_m, tensor_rank_4_m)| tensor_rank_4_m * self_m)
910            .sum()
911    }
912}
913
914impl<const D: usize, const J: usize, const K: usize, const L: usize, const M: usize>
915    Mul<&TensorRank4<D, M, J, K, L>> for &TensorRank1<D, M>
916{
917    type Output = TensorRank3<D, J, K, L>;
918    fn mul(self, tensor_rank_4: &TensorRank4<D, M, J, K, L>) -> Self::Output {
919        self.iter()
920            .zip(tensor_rank_4.iter())
921            .map(|(self_m, tensor_rank_4_m)| tensor_rank_4_m * self_m)
922            .sum()
923    }
924}
925
926impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize, const M: usize>
927    Mul<TensorRank4<D, M, J, K, L>> for TensorRank2<D, I, M>
928{
929    type Output = TensorRank4<D, I, J, K, L>;
930    fn mul(self, tensor_rank_4: TensorRank4<D, M, J, K, L>) -> Self::Output {
931        self.into_iter()
932            .map(|self_i| self_i * &tensor_rank_4)
933            .collect()
934    }
935}
936
937impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize, const M: usize>
938    Mul<&TensorRank4<D, M, J, K, L>> for TensorRank2<D, I, M>
939{
940    type Output = TensorRank4<D, I, J, K, L>;
941    fn mul(self, tensor_rank_4: &TensorRank4<D, M, J, K, L>) -> Self::Output {
942        self.into_iter()
943            .map(|self_i| self_i * tensor_rank_4)
944            .collect()
945    }
946}
947
948impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize, const M: usize>
949    Mul<TensorRank4<D, M, J, K, L>> for &TensorRank2<D, I, M>
950{
951    type Output = TensorRank4<D, I, J, K, L>;
952    fn mul(self, tensor_rank_4: TensorRank4<D, M, J, K, L>) -> Self::Output {
953        self.iter().map(|self_i| self_i * &tensor_rank_4).collect()
954    }
955}
956
957impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize, const M: usize>
958    Mul<&TensorRank4<D, M, J, K, L>> for &TensorRank2<D, I, M>
959{
960    type Output = TensorRank4<D, I, J, K, L>;
961    fn mul(self, tensor_rank_4: &TensorRank4<D, M, J, K, L>) -> Self::Output {
962        self.iter().map(|self_i| self_i * tensor_rank_4).collect()
963    }
964}
965
966impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> Add
967    for TensorRank4<D, I, J, K, L>
968{
969    type Output = Self;
970    fn add(mut self, tensor_rank_4: Self) -> Self::Output {
971        self += tensor_rank_4;
972        self
973    }
974}
975
976impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> Add<&Self>
977    for TensorRank4<D, I, J, K, L>
978{
979    type Output = Self;
980    fn add(mut self, tensor_rank_4: &Self) -> Self::Output {
981        self += tensor_rank_4;
982        self
983    }
984}
985
986impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
987    Add<TensorRank4<D, I, J, K, L>> for &TensorRank4<D, I, J, K, L>
988{
989    type Output = TensorRank4<D, I, J, K, L>;
990    fn add(self, mut tensor_rank_4: TensorRank4<D, I, J, K, L>) -> Self::Output {
991        tensor_rank_4 += self;
992        tensor_rank_4
993    }
994}
995
996impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> AddAssign
997    for TensorRank4<D, I, J, K, L>
998{
999    fn add_assign(&mut self, tensor_rank_4: Self) {
1000        self.iter_mut()
1001            .zip(tensor_rank_4)
1002            .for_each(|(self_i, tensor_rank_4_i)| *self_i += tensor_rank_4_i);
1003    }
1004}
1005
1006impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
1007    AddAssign<&Self> for TensorRank4<D, I, J, K, L>
1008{
1009    fn add_assign(&mut self, tensor_rank_4: &Self) {
1010        self.iter_mut()
1011            .zip(tensor_rank_4.iter())
1012            .for_each(|(self_i, tensor_rank_4_i)| *self_i += tensor_rank_4_i);
1013    }
1014}
1015
1016impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> Sub
1017    for TensorRank4<D, I, J, K, L>
1018{
1019    type Output = Self;
1020    fn sub(mut self, tensor_rank_4: Self) -> Self::Output {
1021        self -= tensor_rank_4;
1022        self
1023    }
1024}
1025
1026impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> Sub<&Self>
1027    for TensorRank4<D, I, J, K, L>
1028{
1029    type Output = Self;
1030    fn sub(mut self, tensor_rank_4: &Self) -> Self::Output {
1031        self -= tensor_rank_4;
1032        self
1033    }
1034}
1035
1036impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> Sub
1037    for &TensorRank4<D, I, J, K, L>
1038{
1039    type Output = TensorRank4<D, I, J, K, L>;
1040    fn sub(self, tensor_rank_4: Self) -> Self::Output {
1041        tensor_rank_4
1042            .iter()
1043            .zip(self.iter())
1044            .map(|(tensor_rank_4_i, self_i)| self_i - tensor_rank_4_i)
1045            .collect()
1046    }
1047}
1048
1049impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> SubAssign
1050    for TensorRank4<D, I, J, K, L>
1051{
1052    fn sub_assign(&mut self, tensor_rank_4: Self) {
1053        self.iter_mut()
1054            .zip(tensor_rank_4)
1055            .for_each(|(self_i, tensor_rank_4_i)| *self_i -= tensor_rank_4_i);
1056    }
1057}
1058
1059impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
1060    SubAssign<&Self> for TensorRank4<D, I, J, K, L>
1061{
1062    fn sub_assign(&mut self, tensor_rank_4: &Self) {
1063        self.iter_mut()
1064            .zip(tensor_rank_4.iter())
1065            .for_each(|(self_i, tensor_rank_4_i)| *self_i -= tensor_rank_4_i);
1066    }
1067}