conspire/math/matrix/vector/
mod.rs

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

conspire/math/matrix/vector/
mod.rs
