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
426pub trait ContractAllIndicesWithFirstIndicesOf<TIM, TJN, TKO, TLP> {
427    type Output;
428    fn contract_all_indices_with_first_indices_of(
429        self,
430        tensor_rank_2_a: TIM,
431        tensor_rank_2_b: TJN,
432        tensor_rank_2_c: TKO,
433        tensor_rank_2_d: TLP,
434    ) -> Self::Output;
435}
436
437impl<
438    const D: usize,
439    const I: usize,
440    const J: usize,
441    const K: usize,
442    const L: usize,
443    const M: usize,
444    const N: usize,
445    const O: usize,
446    const P: usize,
447>
448    ContractAllIndicesWithFirstIndicesOf<
449        &TensorRank2<D, I, M>,
450        &TensorRank2<D, J, N>,
451        &TensorRank2<D, K, O>,
452        &TensorRank2<D, L, P>,
453    > for TensorRank4<D, I, J, K, L>
454{
455    type Output = TensorRank4<D, M, N, O, P>;
456    fn contract_all_indices_with_first_indices_of(
457        self,
458        tensor_rank_2_a: &TensorRank2<D, I, M>,
459        tensor_rank_2_b: &TensorRank2<D, J, N>,
460        tensor_rank_2_c: &TensorRank2<D, K, O>,
461        tensor_rank_2_d: &TensorRank2<D, L, P>,
462    ) -> Self::Output {
463        let mut output = TensorRank4::zero();
464        self.into_iter().zip(tensor_rank_2_a.iter()).for_each(|(self_m, tensor_rank_2_a_m)|
465            self_m.into_iter().zip(tensor_rank_2_b.iter()).for_each(|(self_mn, tensor_rank_2_b_n)|
466                self_mn.into_iter().zip(tensor_rank_2_c.iter()).for_each(|(self_mno, tensor_rank_2_c_o)|
467                    self_mno.into_iter().zip(tensor_rank_2_d.iter()).for_each(|(self_mnop, tensor_rank_2_d_p)|
468                        output.iter_mut().zip(tensor_rank_2_a_m.iter()).for_each(|(output_i, tensor_rank_2_a_mi)|
469                            output_i.iter_mut().zip(tensor_rank_2_b_n.iter()).for_each(|(output_ij, tensor_rank_2_b_nj)|
470                                output_ij.iter_mut().zip(tensor_rank_2_c_o.iter()).for_each(|(output_ijk, tensor_rank_2_c_ok)|
471                                    output_ijk.iter_mut().zip(tensor_rank_2_d_p.iter()).for_each(|(output_ijkl, tensor_rank_2_dp)|
472                                        *output_ijkl += self_mnop * tensor_rank_2_a_mi * tensor_rank_2_b_nj * tensor_rank_2_c_ok * tensor_rank_2_dp
473                                    )
474                                )
475                            )
476                        )
477                    )
478                )
479            )
480        );
481        output
482    }
483}
484
485pub trait ContractFirstThirdFourthIndicesWithFirstIndicesOf<TIM, TKO, TLP> {
486    type Output;
487    fn contract_first_third_fourth_indices_with_first_indices_of(
488        self,
489        tensor_rank_2_a: TIM,
490        tensor_rank_2_c: TKO,
491        tensor_rank_2_d: TLP,
492    ) -> Self::Output;
493}
494
495impl<
496    const D: usize,
497    const I: usize,
498    const J: usize,
499    const K: usize,
500    const L: usize,
501    const M: usize,
502    const O: usize,
503    const P: usize,
504>
505    ContractFirstThirdFourthIndicesWithFirstIndicesOf<
506        &TensorRank2<D, I, M>,
507        &TensorRank2<D, K, O>,
508        &TensorRank2<D, L, P>,
509    > for TensorRank4<D, I, J, K, L>
510{
511    type Output = TensorRank4<D, M, J, O, P>;
512    fn contract_first_third_fourth_indices_with_first_indices_of(
513        self,
514        tensor_rank_2_a: &TensorRank2<D, I, M>,
515        tensor_rank_2_b: &TensorRank2<D, K, O>,
516        tensor_rank_2_c: &TensorRank2<D, L, P>,
517    ) -> Self::Output {
518        let mut output = TensorRank4::zero();
519        self.into_iter().zip(tensor_rank_2_a.iter()).for_each(|(self_q, tensor_rank_2_a_q)|
520            output.iter_mut().zip(tensor_rank_2_a_q.iter()).for_each(|(output_i, tensor_rank_2_a_qi)|
521                output_i.iter_mut().zip(self_q.iter()).for_each(|(output_ij, self_qj)|
522                    self_qj.iter().zip(tensor_rank_2_b.iter()).for_each(|(self_qjm, tensor_rank_2_b_m)|
523                        self_qjm.iter().zip(tensor_rank_2_c.iter()).for_each(|(self_qjmn, tensor_rank_2_c_n)|
524                            output_ij.iter_mut().zip(tensor_rank_2_b_m.iter()).for_each(|(output_ijk, tensor_rank_2_b_mk)|
525                                output_ijk.iter_mut().zip(tensor_rank_2_c_n.iter()).for_each(|(output_ijkl, tensor_rank_2_c_nl)|
526                                    *output_ijkl += self_qjmn*tensor_rank_2_a_qi*tensor_rank_2_b_mk*tensor_rank_2_c_nl
527                                )
528                            )
529                        )
530                    )
531                )
532            )
533        );
534        output
535    }
536}
537
538pub trait ContractSecondIndexWithFirstIndexOf<TJN> {
539    type Output;
540    fn contract_second_index_with_first_index_of(self, tensor_rank_2: TJN) -> Self::Output;
541}
542
543impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize, const N: usize>
544    ContractSecondIndexWithFirstIndexOf<&TensorRank2<D, J, N>> for TensorRank4<D, I, J, K, L>
545{
546    type Output = TensorRank4<D, I, N, K, L>;
547    fn contract_second_index_with_first_index_of(
548        self,
549        tensor_rank_2: &TensorRank2<D, J, N>,
550    ) -> Self::Output {
551        let mut output = TensorRank4::zero();
552        output.iter_mut().zip(self).for_each(|(output_i, self_i)| {
553            self_i
554                .into_iter()
555                .zip(tensor_rank_2.iter())
556                .for_each(|(self_is, tensor_rank_2_s)| {
557                    output_i.iter_mut().zip(tensor_rank_2_s.iter()).for_each(
558                        |(output_ij, tensor_rank_2_sj)| *output_ij += &self_is * tensor_rank_2_sj,
559                    )
560                })
561        });
562        output
563    }
564}
565
566pub trait ContractSecondFourthIndicesWithFirstIndicesOf<TJ, TL> {
567    type Output;
568    fn contract_second_fourth_indices_with_first_indices_of(
569        &self,
570        object_a: TJ,
571        object_b: TL,
572    ) -> Self::Output;
573}
574
575impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
576    ContractSecondFourthIndicesWithFirstIndicesOf<&TensorRank1<D, J>, &TensorRank1<D, L>>
577    for TensorRank4<D, I, J, K, L>
578{
579    type Output = TensorRank2<D, I, K>;
580    fn contract_second_fourth_indices_with_first_indices_of(
581        &self,
582        tensor_rank_1_a: &TensorRank1<D, J>,
583        tensor_rank_1_b: &TensorRank1<D, L>,
584    ) -> Self::Output {
585        self.iter()
586            .map(|self_i| {
587                self_i
588                    .iter()
589                    .zip(tensor_rank_1_a.iter())
590                    .map(|(self_ij, tensor_rank_1_a_j)| {
591                        self_ij
592                            .iter()
593                            .map(|self_ijk| self_ijk * (tensor_rank_1_b * tensor_rank_1_a_j))
594                            .collect()
595                    })
596                    .sum()
597            })
598            .collect()
599    }
600}
601
602pub trait ContractThirdIndexWithFirstIndexOf<TKL> {
603    type Output;
604    fn contract_third_index_with_first_index_of(&self, tensor: TKL) -> Self::Output;
605}
606
607impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize, const M: usize>
608    ContractThirdIndexWithFirstIndexOf<&TensorRank2<D, M, K>> for TensorRank4<D, I, J, M, L>
609{
610    type Output = TensorRank4<D, I, J, K, L>;
611    fn contract_third_index_with_first_index_of(
612        &self,
613        tensor_rank_2: &TensorRank2<D, M, K>,
614    ) -> Self::Output {
615        self.iter()
616            .map(|self_i| {
617                self_i
618                    .iter()
619                    .map(|self_ij| {
620                        self_ij
621                            .iter()
622                            .zip(tensor_rank_2.iter())
623                            .map(|(self_ijm, tensor_rank_2_m)| (tensor_rank_2_m, self_ijm).into())
624                            .sum()
625                    })
626                    .collect()
627            })
628            .collect()
629    }
630}
631
632pub trait ContractThirdFourthIndicesWithFirstSecondIndicesOf<TKL> {
633    type Output;
634    fn contract_third_fourth_indices_with_first_second_indices_of(
635        self,
636        tensor: TKL,
637    ) -> Self::Output;
638}
639
640impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
641    ContractThirdFourthIndicesWithFirstSecondIndicesOf<&TensorRank2<D, K, L>>
642    for TensorRank4<D, I, J, K, L>
643{
644    type Output = TensorRank2<D, I, J>;
645    fn contract_third_fourth_indices_with_first_second_indices_of(
646        self,
647        tensor_rank_2: &TensorRank2<D, K, L>,
648    ) -> Self::Output {
649        self.into_iter()
650            .map(|self_i| {
651                self_i
652                    .into_iter()
653                    .map(|self_ij| self_ij.full_contraction(tensor_rank_2))
654                    .collect()
655            })
656            .collect()
657    }
658}
659
660impl<
661    const D: usize,
662    const I: usize,
663    const J: usize,
664    const K: usize,
665    const L: usize,
666    const M: usize,
667    const N: usize,
668> ContractThirdFourthIndicesWithFirstSecondIndicesOf<&TensorRank4<D, K, L, M, N>>
669    for TensorRank4<D, I, J, K, L>
670{
671    type Output = TensorRank4<D, I, J, M, N>;
672    fn contract_third_fourth_indices_with_first_second_indices_of(
673        self,
674        tensor: &TensorRank4<D, K, L, M, N>,
675    ) -> Self::Output {
676        self.into_iter()
677            .map(|self_i| {
678                self_i
679                    .into_iter()
680                    .map(|self_ij| {
681                        self_ij
682                            .into_iter()
683                            .zip(tensor.iter())
684                            .map(|(self_ijk, tensor_k)| {
685                                self_ijk
686                                    .into_iter()
687                                    .zip(tensor_k.iter())
688                                    .map(|(self_ijkl, tensor_kl)| tensor_kl * self_ijkl)
689                                    .sum()
690                            })
691                            .sum()
692                    })
693                    .collect()
694            })
695            .collect()
696    }
697}
698
699pub trait ContractFirstSecondIndicesWithSecondIndicesOf<TI, TJ> {
700    type Output;
701    fn contract_first_second_indices_with_second_indices_of(
702        self,
703        object_a: TI,
704        object_b: TJ,
705    ) -> Self::Output;
706}
707
708impl<
709    const D: usize,
710    const I: usize,
711    const J: usize,
712    const K: usize,
713    const L: usize,
714    const M: usize,
715    const N: usize,
716> ContractFirstSecondIndicesWithSecondIndicesOf<&TensorRank2<D, I, M>, &TensorRank2<D, J, N>>
717    for TensorRank4<D, M, N, K, L>
718{
719    type Output = TensorRank4<D, I, J, K, L>;
720    fn contract_first_second_indices_with_second_indices_of(
721        self,
722        tensor_rank_2_a: &TensorRank2<D, I, M>,
723        tensor_rank_2_b: &TensorRank2<D, J, N>,
724    ) -> Self::Output {
725        let mut output = TensorRank4::zero();
726        output
727            .iter_mut()
728            .zip(tensor_rank_2_a.iter())
729            .for_each(|(output_i, tensor_rank_2_a_i)| {
730                output_i.iter_mut().zip(tensor_rank_2_b.iter()).for_each(
731                    |(output_ij, tensor_rank_2_b_j)| {
732                        self.iter().zip(tensor_rank_2_a_i.iter()).for_each(
733                            |(self_m, tensor_rank_2_a_im)| {
734                                self_m.iter().zip(tensor_rank_2_b_j.iter()).for_each(
735                                    |(self_mn, tensor_rank_2_b_jn)| {
736                                        *output_ij +=
737                                            self_mn * tensor_rank_2_a_im * tensor_rank_2_b_jn
738                                    },
739                                )
740                            },
741                        )
742                    },
743                )
744            });
745        output
746    }
747}
748
749impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
750    Div<TensorRank0> for TensorRank4<D, I, J, K, L>
751{
752    type Output = Self;
753    fn div(mut self, tensor_rank_0: TensorRank0) -> Self::Output {
754        self /= &tensor_rank_0;
755        self
756    }
757}
758
759impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
760    Div<TensorRank0> for &TensorRank4<D, I, J, K, L>
761{
762    type Output = TensorRank4<D, I, J, K, L>;
763    fn div(self, tensor_rank_0: TensorRank0) -> Self::Output {
764        self.iter().map(|self_i| self_i / tensor_rank_0).collect()
765    }
766}
767
768impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
769    Div<&TensorRank0> for TensorRank4<D, I, J, K, L>
770{
771    type Output = Self;
772    fn div(mut self, tensor_rank_0: &TensorRank0) -> Self::Output {
773        self /= tensor_rank_0;
774        self
775    }
776}
777
778impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
779    DivAssign<TensorRank0> for TensorRank4<D, I, J, K, L>
780{
781    fn div_assign(&mut self, tensor_rank_0: TensorRank0) {
782        self.iter_mut().for_each(|self_i| *self_i /= &tensor_rank_0);
783    }
784}
785
786impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
787    DivAssign<&TensorRank0> for TensorRank4<D, I, J, K, L>
788{
789    fn div_assign(&mut self, tensor_rank_0: &TensorRank0) {
790        self.iter_mut().for_each(|self_i| *self_i /= tensor_rank_0);
791    }
792}
793
794impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
795    Mul<TensorRank0> for TensorRank4<D, I, J, K, L>
796{
797    type Output = Self;
798    fn mul(mut self, tensor_rank_0: TensorRank0) -> Self::Output {
799        self *= &tensor_rank_0;
800        self
801    }
802}
803
804impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
805    Mul<&TensorRank0> for TensorRank4<D, I, J, K, L>
806{
807    type Output = Self;
808    fn mul(mut self, tensor_rank_0: &TensorRank0) -> Self::Output {
809        self *= tensor_rank_0;
810        self
811    }
812}
813
814impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
815    MulAssign<TensorRank0> for TensorRank4<D, I, J, K, L>
816{
817    fn mul_assign(&mut self, tensor_rank_0: TensorRank0) {
818        self.iter_mut().for_each(|self_i| *self_i *= &tensor_rank_0);
819    }
820}
821
822impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
823    MulAssign<&TensorRank0> for TensorRank4<D, I, J, K, L>
824{
825    fn mul_assign(&mut self, tensor_rank_0: &TensorRank0) {
826        self.iter_mut().for_each(|self_i| *self_i *= tensor_rank_0);
827    }
828}
829
830impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize, const M: usize>
831    Mul<TensorRank2<D, L, M>> for TensorRank4<D, I, J, K, L>
832{
833    type Output = TensorRank4<D, I, J, K, M>;
834    fn mul(self, tensor_rank_2: TensorRank2<D, L, M>) -> Self::Output {
835        self.into_iter()
836            .map(|self_i| {
837                self_i
838                    .into_iter()
839                    .map(|self_ij| self_ij * &tensor_rank_2)
840                    .collect()
841            })
842            .collect()
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 J: usize, const K: usize, const L: usize, const M: usize>
863    Mul<TensorRank4<D, M, J, K, L>> for TensorRank1<D, M>
864{
865    type Output = TensorRank3<D, J, K, L>;
866    fn mul(self, tensor_rank_4: TensorRank4<D, M, J, K, L>) -> Self::Output {
867        self.into_iter()
868            .zip(tensor_rank_4)
869            .map(|(self_m, tensor_rank_4_m)| tensor_rank_4_m * self_m)
870            .sum()
871    }
872}
873
874impl<const D: usize, const J: usize, const K: usize, const L: usize, const M: usize>
875    Mul<&TensorRank4<D, M, J, K, L>> for TensorRank1<D, M>
876{
877    type Output = TensorRank3<D, J, K, L>;
878    fn mul(self, tensor_rank_4: &TensorRank4<D, M, J, K, L>) -> Self::Output {
879        self.into_iter()
880            .zip(tensor_rank_4.iter())
881            .map(|(self_m, tensor_rank_4_m)| tensor_rank_4_m * self_m)
882            .sum()
883    }
884}
885
886impl<const D: usize, const J: usize, const K: usize, const L: usize, const M: usize>
887    Mul<TensorRank4<D, M, J, K, L>> for &TensorRank1<D, M>
888{
889    type Output = TensorRank3<D, J, K, L>;
890    fn mul(self, tensor_rank_4: TensorRank4<D, M, J, K, L>) -> Self::Output {
891        self.iter()
892            .zip(tensor_rank_4)
893            .map(|(self_m, tensor_rank_4_m)| tensor_rank_4_m * self_m)
894            .sum()
895    }
896}
897
898impl<const D: usize, const J: usize, const K: usize, const L: usize, const M: usize>
899    Mul<&TensorRank4<D, M, J, K, L>> for &TensorRank1<D, M>
900{
901    type Output = TensorRank3<D, J, K, L>;
902    fn mul(self, tensor_rank_4: &TensorRank4<D, M, J, K, L>) -> Self::Output {
903        self.iter()
904            .zip(tensor_rank_4.iter())
905            .map(|(self_m, tensor_rank_4_m)| tensor_rank_4_m * self_m)
906            .sum()
907    }
908}
909
910impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize, const M: usize>
911    Mul<TensorRank4<D, M, J, K, L>> for TensorRank2<D, I, M>
912{
913    type Output = TensorRank4<D, I, J, K, L>;
914    fn mul(self, tensor_rank_4: TensorRank4<D, M, J, K, L>) -> Self::Output {
915        self.into_iter()
916            .map(|self_i| self_i * &tensor_rank_4)
917            .collect()
918    }
919}
920
921impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize, const M: usize>
922    Mul<&TensorRank4<D, M, J, K, L>> for TensorRank2<D, I, M>
923{
924    type Output = TensorRank4<D, I, J, K, L>;
925    fn mul(self, tensor_rank_4: &TensorRank4<D, M, J, K, L>) -> Self::Output {
926        self.into_iter()
927            .map(|self_i| self_i * tensor_rank_4)
928            .collect()
929    }
930}
931
932impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize, const M: usize>
933    Mul<TensorRank4<D, M, J, K, L>> for &TensorRank2<D, I, M>
934{
935    type Output = TensorRank4<D, I, J, K, L>;
936    fn mul(self, tensor_rank_4: TensorRank4<D, M, J, K, L>) -> Self::Output {
937        self.iter().map(|self_i| self_i * &tensor_rank_4).collect()
938    }
939}
940
941impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize, const M: usize>
942    Mul<&TensorRank4<D, M, J, K, L>> for &TensorRank2<D, I, M>
943{
944    type Output = TensorRank4<D, I, J, K, L>;
945    fn mul(self, tensor_rank_4: &TensorRank4<D, M, J, K, L>) -> Self::Output {
946        self.iter().map(|self_i| self_i * tensor_rank_4).collect()
947    }
948}
949
950impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> Add
951    for TensorRank4<D, I, J, K, L>
952{
953    type Output = Self;
954    fn add(mut self, tensor_rank_4: Self) -> Self::Output {
955        self += tensor_rank_4;
956        self
957    }
958}
959
960impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> Add<&Self>
961    for TensorRank4<D, I, J, K, L>
962{
963    type Output = Self;
964    fn add(mut self, tensor_rank_4: &Self) -> Self::Output {
965        self += tensor_rank_4;
966        self
967    }
968}
969
970impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
971    Add<TensorRank4<D, I, J, K, L>> for &TensorRank4<D, I, J, K, L>
972{
973    type Output = TensorRank4<D, I, J, K, L>;
974    fn add(self, mut tensor_rank_4: TensorRank4<D, I, J, K, L>) -> Self::Output {
975        tensor_rank_4 += self;
976        tensor_rank_4
977    }
978}
979
980impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> AddAssign
981    for TensorRank4<D, I, J, K, L>
982{
983    fn add_assign(&mut self, tensor_rank_4: Self) {
984        self.iter_mut()
985            .zip(tensor_rank_4)
986            .for_each(|(self_i, tensor_rank_4_i)| *self_i += tensor_rank_4_i);
987    }
988}
989
990impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
991    AddAssign<&Self> for TensorRank4<D, I, J, K, L>
992{
993    fn add_assign(&mut self, tensor_rank_4: &Self) {
994        self.iter_mut()
995            .zip(tensor_rank_4.iter())
996            .for_each(|(self_i, tensor_rank_4_i)| *self_i += tensor_rank_4_i);
997    }
998}
999
1000impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> Sub
1001    for TensorRank4<D, I, J, K, L>
1002{
1003    type Output = Self;
1004    fn sub(mut self, tensor_rank_4: Self) -> Self::Output {
1005        self -= tensor_rank_4;
1006        self
1007    }
1008}
1009
1010impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> Sub<&Self>
1011    for TensorRank4<D, I, J, K, L>
1012{
1013    type Output = Self;
1014    fn sub(mut self, tensor_rank_4: &Self) -> Self::Output {
1015        self -= tensor_rank_4;
1016        self
1017    }
1018}
1019
1020impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> Sub
1021    for &TensorRank4<D, I, J, K, L>
1022{
1023    type Output = TensorRank4<D, I, J, K, L>;
1024    fn sub(self, tensor_rank_4: Self) -> Self::Output {
1025        tensor_rank_4
1026            .iter()
1027            .zip(self.iter())
1028            .map(|(tensor_rank_4_i, self_i)| self_i - tensor_rank_4_i)
1029            .collect()
1030    }
1031}
1032
1033impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize> SubAssign
1034    for TensorRank4<D, I, J, K, L>
1035{
1036    fn sub_assign(&mut self, tensor_rank_4: Self) {
1037        self.iter_mut()
1038            .zip(tensor_rank_4)
1039            .for_each(|(self_i, tensor_rank_4_i)| *self_i -= tensor_rank_4_i);
1040    }
1041}
1042
1043impl<const D: usize, const I: usize, const J: usize, const K: usize, const L: usize>
1044    SubAssign<&Self> for TensorRank4<D, I, J, K, L>
1045{
1046    fn sub_assign(&mut self, tensor_rank_4: &Self) {
1047        self.iter_mut()
1048            .zip(tensor_rank_4.iter())
1049            .for_each(|(self_i, tensor_rank_4_i)| *self_i -= tensor_rank_4_i);
1050    }
1051}