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