conspire/math/matrix/vector/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::math::assert::FiniteDifference;
5use crate::units::Dimensionless;
6
7use crate::math::{
8    Erase, Jacobian, Matrix, Quantity, QuantityVector, Scalar, Solution, SquareMatrix, Tensor,
9    TensorRank1Vec, TensorRank2, TensorTuple, TensorVec, write_tensor_rank_0,
10};
11use std::{
12    fmt::{Display, Formatter, Result},
13    iter::Sum,
14    mem::forget,
15    ops::{
16        Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, RangeFrom, RangeTo, Sub,
17        SubAssign,
18    },
19    slice, vec,
20};
21
22/// A vector.
23#[derive(Clone, Debug, PartialEq)]
24pub struct Vector(Vec<Scalar>);
25
26impl Vector {
27    /// Returns a raw pointer to the vector’s buffer, or a dangling raw pointer valid for zero sized reads if the vector didn’t allocate.
28    pub const fn as_ptr(&self) -> *const Scalar {
29        self.0.as_ptr()
30    }
31    pub fn as_slice(&self) -> &[Scalar] {
32        self.0.as_slice()
33    }
34    pub fn as_mut_slice(&mut self) -> &mut [Scalar] {
35        self.0.as_mut_slice()
36    }
37    pub fn ones(len: usize) -> Self {
38        Self(vec![1.0; len])
39    }
40    pub fn zero(len: usize) -> Self {
41        Self(vec![0.0; len])
42    }
43}
44
45impl Default for Vector {
46    fn default() -> Self {
47        Self::new()
48    }
49}
50
51impl FiniteDifference for Vector {
52    fn error_fd(&self, comparator: &Self, epsilon: Scalar) -> Option<(bool, usize)> {
53        let error_count = self
54            .iter()
55            .zip(comparator.iter())
56            .map(|(entry, comparator_entry)| {
57                entry
58                    .iter()
59                    .zip(comparator_entry.iter())
60                    .filter(|&(&entry_i, &comparator_entry_i)| {
61                        (entry_i / comparator_entry_i - 1.0).abs() >= epsilon
62                            && (entry_i.abs() >= epsilon || comparator_entry_i.abs() >= epsilon)
63                    })
64                    .count()
65            })
66            .sum();
67        if error_count > 0 {
68            let auxiliary = self
69                .iter()
70                .zip(comparator.iter())
71                .map(|(entry, comparator_entry)| {
72                    entry
73                        .iter()
74                        .zip(comparator_entry.iter())
75                        .filter(|&(&entry_i, &comparator_entry_i)| {
76                            (entry_i / comparator_entry_i - 1.0).abs() >= epsilon
77                                && (entry_i - comparator_entry_i).abs() >= epsilon
78                                && (entry_i.abs() >= epsilon || comparator_entry_i.abs() >= epsilon)
79                        })
80                        .count()
81                })
82                .sum::<usize>()
83                > 0;
84            Some((auxiliary, error_count))
85        } else {
86            None
87        }
88    }
89}
90
91impl Display for Vector {
92    fn fmt(&self, f: &mut Formatter) -> Result {
93        write!(f, "\x1B[s")?;
94        write!(f, "[")?;
95        self.0.chunks(5).enumerate().try_for_each(|(i, chunk)| {
96            chunk
97                .iter()
98                .try_for_each(|entry| write_tensor_rank_0(f, entry))?;
99            if (i + 1) * 5 < self.len() {
100                writeln!(f, "\x1B[2D,")?;
101                write!(f, "\x1B[u")?;
102                write!(f, "\x1B[{}B ", i + 1)?;
103            }
104            Ok(())
105        })?;
106        write!(f, "\x1B[2D]")?;
107        Ok(())
108    }
109}
110
111impl<const N: usize> From<[Scalar; N]> for Vector {
112    fn from(array: [Scalar; N]) -> Self {
113        Self(array.to_vec())
114    }
115}
116
117impl From<&[Scalar]> for Vector {
118    fn from(slice: &[Scalar]) -> Self {
119        Self(slice.to_vec())
120    }
121}
122
123impl From<Scalar> for Vector {
124    fn from(scalar: Scalar) -> Self {
125        Vector(vec![scalar])
126    }
127}
128
129impl From<Vec<Scalar>> for Vector {
130    fn from(vec: Vec<Scalar>) -> Self {
131        Self(vec)
132    }
133}
134
135impl From<Vector> for Vec<Scalar> {
136    fn from(vector: Vector) -> Self {
137        vector.0
138    }
139}
140
141impl<const D: usize, I> From<TensorRank1Vec<D, I>> for Vector {
142    fn from(tensor_rank_1_vec: TensorRank1Vec<D, I>) -> Self {
143        let length = tensor_rank_1_vec.len() * D;
144        let capacity = tensor_rank_1_vec.capacity() * D;
145        let pointer = tensor_rank_1_vec.as_ptr() as *mut Scalar;
146        forget(tensor_rank_1_vec);
147        unsafe { Self(Vec::from_raw_parts(pointer, length, capacity)) }
148    }
149}
150
151impl<const D: usize, I, J> From<TensorRank2<D, I, J>> for Vector {
152    fn from(tensor_rank_2: TensorRank2<D, I, J>) -> Self {
153        let length = D * D;
154        let capacity = length;
155        let pointer = tensor_rank_2.as_ptr() as *mut Scalar;
156        unsafe { Self(Vec::from_raw_parts(pointer, length, capacity)) }
157    }
158}
159
160impl FromIterator<Scalar> for Vector {
161    fn from_iter<Ii: IntoIterator<Item = Scalar>>(into_iterator: Ii) -> Self {
162        Self(Vec::from_iter(into_iterator))
163    }
164}
165
166impl Index<usize> for Vector {
167    type Output = Scalar;
168    fn index(&self, index: usize) -> &Self::Output {
169        &self.0[index]
170    }
171}
172
173impl Index<RangeTo<usize>> for Vector {
174    type Output = [Scalar];
175    fn index(&self, indices: RangeTo<usize>) -> &Self::Output {
176        &self.0[indices]
177    }
178}
179
180impl Index<RangeFrom<usize>> for Vector {
181    type Output = [Scalar];
182    fn index(&self, indices: RangeFrom<usize>) -> &Self::Output {
183        &self.0[indices]
184    }
185}
186
187impl IndexMut<usize> for Vector {
188    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
189        &mut self.0[index]
190    }
191}
192
193impl Tensor for Vector {
194    type Item = Scalar;
195    type Unit = Dimensionless;
196    fn iter(&self) -> impl Iterator<Item = &Self::Item> {
197        self.0.iter()
198    }
199    fn iter_mut(&mut self) -> impl Iterator<Item = &mut Self::Item> {
200        self.0.iter_mut()
201    }
202    fn len(&self) -> usize {
203        self.0.len()
204    }
205    fn norm_inf(&self) -> Quantity<Dimensionless> {
206        Quantity::new(self.iter().fold(0.0, |acc, entry| entry.abs().max(acc)))
207    }
208    fn size(&self) -> usize {
209        self.len()
210    }
211}
212
213impl Solution for Vector {
214    fn decrement_from(&mut self, other: &Vector) {
215        self.iter_mut()
216            .zip(other.iter())
217            .for_each(|(self_i, vector_i)| *self_i -= vector_i)
218    }
219    fn decrement_from_chained(&mut self, other: &mut Self, vector: &Vector) {
220        self.iter_mut()
221            .chain(other.iter_mut())
222            .zip(vector.iter())
223            .for_each(|(entry_i, vector_i)| *entry_i -= vector_i)
224    }
225    fn decrement_from_retained(&mut self, retained: &[bool], other: &Vector) {
226        self.iter_mut()
227            .zip(retained.iter())
228            .filter(|(_, retained_i)| **retained_i)
229            .zip(other.iter())
230            .for_each(|((self_i, _), vector_i)| *self_i -= vector_i)
231    }
232}
233
234impl Jacobian for Vector {
235    fn fill_into(&self, vector: &mut Vector) {
236        self.iter()
237            .zip(vector.iter_mut())
238            .for_each(|(self_i, vector_i)| *vector_i = *self_i)
239    }
240    fn fill_into_chained(self, other: Self, vector: &mut Self) {
241        self.into_iter()
242            .chain(other)
243            .zip(vector.iter_mut())
244            .for_each(|(entry_i, vector_i)| *vector_i = entry_i)
245    }
246    fn retain_from(self, retained: &[bool]) -> Vector {
247        self.into_iter()
248            .zip(retained.iter())
249            .filter(|(_, retained_i)| **retained_i)
250            .map(|(entry, _)| entry)
251            .collect()
252    }
253    fn zero_out(&mut self, indices: &[usize]) {
254        indices.iter().for_each(|&index| self[index] = 0.0)
255    }
256}
257
258impl IntoIterator for Vector {
259    type Item = Scalar;
260    type IntoIter = vec::IntoIter<Self::Item>;
261    fn into_iter(self) -> Self::IntoIter {
262        self.0.into_iter()
263    }
264}
265
266impl<'a> IntoIterator for &'a Vector {
267    type Item = &'a Scalar;
268    type IntoIter = slice::Iter<'a, Scalar>;
269    fn into_iter(self) -> Self::IntoIter {
270        self.0.iter()
271    }
272}
273
274impl Extend<Scalar> for Vector {
275    fn extend<I>(&mut self, iter: I)
276    where
277        I: IntoIterator<Item = Scalar>,
278    {
279        self.0.extend(iter)
280    }
281}
282
283impl TensorVec for Vector {
284    type Item = Scalar;
285    fn append(&mut self, other: &mut Self) {
286        self.0.append(&mut other.0)
287    }
288    fn capacity(&self) -> usize {
289        self.0.capacity()
290    }
291    fn is_empty(&self) -> bool {
292        self.0.is_empty()
293    }
294    fn new() -> Self {
295        Self(Vec::new())
296    }
297    fn push(&mut self, item: Self::Item) {
298        self.0.push(item)
299    }
300    fn remove(&mut self, index: usize) -> Self::Item {
301        self.0.remove(index)
302    }
303    fn reserve(&mut self, additional: usize) {
304        self.0.reserve(additional)
305    }
306    fn retain<F>(&mut self, f: F)
307    where
308        F: FnMut(&Self::Item) -> bool,
309    {
310        self.0.retain(f)
311    }
312    fn swap_remove(&mut self, index: usize) -> Self::Item {
313        self.0.swap_remove(index)
314    }
315    fn with_capacity(capacity: usize) -> Self {
316        Self(Vec::with_capacity(capacity))
317    }
318}
319
320impl Sum for Vector {
321    fn sum<Ii>(iter: Ii) -> Self
322    where
323        Ii: Iterator<Item = Self>,
324    {
325        iter.reduce(|mut acc, item| {
326            acc += item;
327            acc
328        })
329        .unwrap_or_else(Self::default)
330    }
331}
332
333impl Div<Scalar> for Vector {
334    type Output = Self;
335    fn div(mut self, scalar: Scalar) -> Self::Output {
336        self /= &scalar;
337        self
338    }
339}
340
341impl Div<&Scalar> for Vector {
342    type Output = Self;
343    fn div(mut self, scalar: &Scalar) -> Self::Output {
344        self /= scalar;
345        self
346    }
347}
348
349impl DivAssign<Scalar> for Vector {
350    fn div_assign(&mut self, scalar: Scalar) {
351        self.iter_mut().for_each(|entry| *entry /= &scalar);
352    }
353}
354
355impl DivAssign<&Scalar> for Vector {
356    fn div_assign(&mut self, scalar: &Scalar) {
357        self.iter_mut().for_each(|entry| *entry /= scalar);
358    }
359}
360
361impl Erase for Vector {
362    type Erased = Self;
363    fn erase(&self) -> &Self {
364        self
365    }
366}
367
368impl Mul<Quantity<Dimensionless>> for Vector {
369    type Output = Self;
370    fn mul(self, quantity: Quantity<Dimensionless>) -> Self::Output {
371        self * quantity.value()
372    }
373}
374
375impl Mul<Quantity<Dimensionless>> for &Vector {
376    type Output = Vector;
377    fn mul(self, quantity: Quantity<Dimensionless>) -> Self::Output {
378        self * quantity.value()
379    }
380}
381
382impl Mul<Scalar> for Vector {
383    type Output = Self;
384    fn mul(mut self, scalar: Scalar) -> Self::Output {
385        self *= &scalar;
386        self
387    }
388}
389
390impl Mul<&Scalar> for Vector {
391    type Output = Self;
392    fn mul(mut self, scalar: &Scalar) -> Self::Output {
393        self *= scalar;
394        self
395    }
396}
397
398impl Mul<Scalar> for &Vector {
399    type Output = Vector;
400    fn mul(self, scalar: Scalar) -> Self::Output {
401        self.iter().map(|self_i| self_i * scalar).collect()
402    }
403}
404
405impl Mul<&Scalar> for &Vector {
406    type Output = Vector;
407    fn mul(self, scalar: &Scalar) -> Self::Output {
408        self.iter().map(|self_i| self_i * scalar).collect()
409    }
410}
411
412impl MulAssign<Scalar> for Vector {
413    fn mul_assign(&mut self, scalar: Scalar) {
414        self.iter_mut().for_each(|entry| *entry *= &scalar);
415    }
416}
417
418impl MulAssign<&Scalar> for Vector {
419    fn mul_assign(&mut self, scalar: &Scalar) {
420        self.iter_mut().for_each(|entry| *entry *= scalar);
421    }
422}
423
424impl Add for Vector {
425    type Output = Self;
426    fn add(mut self, vector: Self) -> Self::Output {
427        self += vector;
428        self
429    }
430}
431
432impl Add<&Self> for Vector {
433    type Output = Self;
434    fn add(mut self, vector: &Self) -> Self::Output {
435        self += vector;
436        self
437    }
438}
439
440impl Add<Vector> for &Vector {
441    type Output = Vector;
442    fn add(self, mut vector: Vector) -> Self::Output {
443        vector += self;
444        vector
445    }
446}
447
448impl Add for &Vector {
449    type Output = Vector;
450    fn add(self, vector: Self) -> Self::Output {
451        vector
452            .iter()
453            .zip(self.iter())
454            .map(|(vector_i, self_i)| self_i + vector_i)
455            .collect()
456    }
457}
458
459impl AddAssign for Vector {
460    fn add_assign(&mut self, vector: Self) {
461        self.iter_mut()
462            .zip(vector.iter())
463            .for_each(|(self_entry, scalar)| *self_entry += scalar);
464    }
465}
466
467impl AddAssign<&Self> for Vector {
468    fn add_assign(&mut self, vector: &Self) {
469        self.iter_mut()
470            .zip(vector.iter())
471            .for_each(|(self_entry, scalar)| *self_entry += scalar);
472    }
473}
474
475impl Mul for Vector {
476    type Output = Scalar;
477    fn mul(self, vector: Self) -> Self::Output {
478        self.iter()
479            .zip(vector.iter())
480            .map(|(self_i, vector_i)| self_i.algebraic_mul(*vector_i))
481            .fold(0.0, f64::algebraic_add)
482    }
483}
484
485impl Mul<&Self> for Vector {
486    type Output = Scalar;
487    fn mul(self, vector: &Self) -> Self::Output {
488        self.iter()
489            .zip(vector.iter())
490            .map(|(self_i, vector_i)| self_i.algebraic_mul(*vector_i))
491            .fold(0.0, f64::algebraic_add)
492    }
493}
494
495impl Mul<Vector> for &Vector {
496    type Output = Scalar;
497    fn mul(self, vector: Vector) -> Self::Output {
498        self.iter()
499            .zip(vector.iter())
500            .map(|(self_i, vector_i)| self_i.algebraic_mul(*vector_i))
501            .fold(0.0, f64::algebraic_add)
502    }
503}
504
505impl Mul for &Vector {
506    type Output = Scalar;
507    fn mul(self, vector: Self) -> Self::Output {
508        self.iter()
509            .zip(vector.iter())
510            .map(|(self_i, vector_i)| self_i.algebraic_mul(*vector_i))
511            .fold(0.0, f64::algebraic_add)
512    }
513}
514
515impl Sub for Vector {
516    type Output = Self;
517    fn sub(mut self, vector: Self) -> Self::Output {
518        self -= vector;
519        self
520    }
521}
522
523impl Sub<&Self> for Vector {
524    type Output = Self;
525    fn sub(mut self, vector: &Self) -> Self::Output {
526        self -= vector;
527        self
528    }
529}
530
531impl Sub<Vector> for &Vector {
532    type Output = Vector;
533    fn sub(self, mut vector: Vector) -> Self::Output {
534        vector
535            .iter_mut()
536            .zip(self.iter())
537            .for_each(|(vector_i, self_i)| *vector_i = self_i - *vector_i);
538        vector
539    }
540}
541
542impl Sub for &Vector {
543    type Output = Vector;
544    fn sub(self, vector: Self) -> Self::Output {
545        vector
546            .iter()
547            .zip(self.iter())
548            .map(|(vector_i, self_i)| self_i - vector_i)
549            .collect()
550    }
551}
552
553impl SubAssign for Vector {
554    fn sub_assign(&mut self, vector: Self) {
555        self.iter_mut()
556            .zip(vector.iter())
557            .for_each(|(self_entry, tensor_rank_1)| *self_entry -= tensor_rank_1);
558    }
559}
560
561impl SubAssign<&Self> for Vector {
562    fn sub_assign(&mut self, vector: &Self) {
563        self.iter_mut()
564            .zip(vector.iter())
565            .for_each(|(self_entry, tensor_rank_1)| *self_entry -= tensor_rank_1);
566    }
567}
568
569impl SubAssign<&[Scalar]> for Vector {
570    fn sub_assign(&mut self, slice: &[Scalar]) {
571        self.iter_mut()
572            .zip(slice.iter())
573            .for_each(|(self_entry, tensor_rank_1)| *self_entry -= tensor_rank_1);
574    }
575}
576
577impl Mul<&Matrix> for &Vector {
578    type Output = Vector;
579    fn mul(self, matrix: &Matrix) -> Self::Output {
580        let mut output = Vector::zero(matrix.width());
581        self.iter()
582            .zip(matrix.iter())
583            .for_each(|(self_i, matrix_i)| {
584                output
585                    .iter_mut()
586                    .zip(matrix_i.iter())
587                    .for_each(|(output_j, matrix_ij)| *output_j += self_i * matrix_ij)
588            });
589        output
590    }
591}
592
593impl<const D: usize, I, U> Mul<&TensorRank1Vec<D, I, U>> for &Vector {
594    type Output = Scalar;
595    fn mul(self, tensor_rank_1_vec: &TensorRank1Vec<D, I, U>) -> Self::Output {
596        tensor_rank_1_vec
597            .iter()
598            .enumerate()
599            .map(|(a, entry_a)| {
600                entry_a
601                    .iter()
602                    .enumerate()
603                    .map(|(i, entry_a_i)| self[D * a + i].algebraic_mul(entry_a_i.value()))
604                    .fold(0.0, f64::algebraic_add)
605            })
606            .fold(0.0, f64::algebraic_add)
607    }
608}
609
610impl<U> Mul<&QuantityVector<U>> for &Vector {
611    type Output = Scalar;
612    fn mul(self, quantity_vector: &QuantityVector<U>) -> Self::Output {
613        quantity_vector
614            .iter()
615            .enumerate()
616            .map(|(a, entry_a)| self[a].algebraic_mul(entry_a.value()))
617            .fold(0.0, f64::algebraic_add)
618    }
619}
620
621impl<const D: usize, I, J, U> Mul<&TensorRank2<D, I, J, U>> for &Vector {
622    type Output = Scalar;
623    fn mul(self, tensor_rank_2: &TensorRank2<D, I, J, U>) -> Self::Output {
624        tensor_rank_2
625            .iter()
626            .enumerate()
627            .map(|(i, entry_i)| {
628                entry_i
629                    .iter()
630                    .enumerate()
631                    .map(|(j, entry_ij)| self[D * i + j].algebraic_mul(entry_ij.value()))
632                    .fold(0.0, f64::algebraic_add)
633            })
634            .fold(0.0, f64::algebraic_add)
635    }
636}
637
638impl<const D: usize, I, J, K, L, U, V>
639    Mul<&TensorTuple<TensorRank2<D, I, J, U>, TensorRank2<D, K, L, V>>> for &Vector
640{
641    type Output = Scalar;
642    fn mul(
643        self,
644        tensor_tuple: &TensorTuple<TensorRank2<D, I, J, U>, TensorRank2<D, K, L, V>>,
645    ) -> Self::Output {
646        let (tensor_rank_2_a, tensor_rank_2_b) = tensor_tuple.into();
647        &self.iter().take(D * D).copied().collect::<Vector>() * tensor_rank_2_a
648            + &self.iter().skip(D * D).copied().collect::<Vector>() * tensor_rank_2_b
649    }
650}
651
652impl Div<SquareMatrix> for &Vector {
653    type Output = Vector;
654    fn div(self, square_matrix: SquareMatrix) -> Self::Output {
655        match square_matrix.solve_lu(self) {
656            Ok(solution) => solution,
657            Err(error) => panic!("{error:?}"),
658        }
659    }
660}
conspire/math/matrix/vector/mod.rs

conspire/math/matrix/vector/
mod.rs
