matrix/
vector.rs

1use std::{
2    fmt::Debug,
3    ops::{Add, AddAssign, Index, Mul, MulAssign, Sub, SubAssign},
4};
5
6use crate::scalar::{MulAdd, Scalar, Sqrt};
7
8#[derive(Clone, Default)]
9pub struct Vector<K> {
10    pub _d: Vec<K>,
11}
12
13pub trait Dot<K> {
14    fn dot(&self, v: &Vector<K>) -> K;
15}
16
17pub trait Angle {
18    type Output;
19    fn angle_cos(u: &Self, v: &Self) -> Self::Output;
20}
21
22pub fn angle_cos<V: Angle>(u: &V, v: &V) -> V::Output {
23    V::angle_cos(u, v)
24}
25
26#[macro_export]
27macro_rules! V {
28    () => {
29        Vector::default()
30    };
31
32    ($values:expr) => {
33        Vector::from($values)
34    };
35}
36
37impl<K: Scalar> Debug for Vector<K> {
38    fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
39        f.debug_list().entries(self._d.iter()).finish()
40    }
41}
42
43impl<K: Clone, const N: usize> From<[K; N]> for Vector<K> {
44    fn from(data: [K; N]) -> Self {
45        Vector {
46            _d: Vec::from(data),
47        }
48    }
49}
50
51impl<K> From<Vec<K>> for Vector<K> {
52    fn from(data: Vec<K>) -> Self {
53        Vector { _d: data }
54    }
55}
56
57impl<K> Index<usize> for Vector<K> {
58    type Output = K;
59
60    fn index(&self, index: usize) -> &Self::Output {
61        &self._d[index]
62    }
63}
64
65impl<K: Scalar> Add for Vector<K> {
66    type Output = Self;
67
68    fn add(self, other: Self) -> Self::Output {
69        assert_eq!(self.size(), other.size(), "vectors must be the same size");
70
71        let mut vec = Vec::with_capacity(self.size());
72        for i in 0..self.size() {
73            vec.push(self[i] + other[i]);
74        }
75
76        V!(vec)
77    }
78}
79
80impl<K: Scalar> AddAssign<&Vector<K>> for Vector<K> {
81    fn add_assign(&mut self, rhs: &Vector<K>) {
82        assert_eq!(self.size(), rhs.size(), "vectors must be the same size");
83
84        for i in 0..self.size() {
85            self._d[i] += rhs[i];
86        }
87    }
88}
89
90impl<K: Scalar> Sub for Vector<K> {
91    type Output = Self;
92
93    fn sub(self, rhs: Self) -> Self::Output {
94        assert_eq!(self.size(), rhs.size(), "vectors must be the same size");
95
96        let mut vec = Vec::with_capacity(self.size());
97        for i in 0..self.size() {
98            vec.push(self[i] - rhs[i]);
99        }
100
101        V!(vec)
102    }
103}
104
105impl<K: Scalar> SubAssign<&Vector<K>> for Vector<K> {
106    fn sub_assign(&mut self, rhs: &Vector<K>) {
107        assert_eq!(self.size(), rhs.size(), "vectors must be the same size");
108
109        for i in 0..self.size() {
110            self._d[i] -= rhs[i];
111        }
112    }
113}
114
115impl<K: Scalar + Mul<U, Output = K>, U: Scalar> Mul<U> for Vector<K> {
116    type Output = Self;
117
118    fn mul(self, a: U) -> Self::Output {
119        &self * &a
120    }
121}
122
123impl<K: Scalar + Mul<U, Output = K>, U: Scalar> Mul<&U> for &Vector<K> {
124    type Output = Vector<K>;
125
126    fn mul(self, &a: &U) -> Self::Output {
127        let mut vec = Vec::with_capacity(self.size());
128
129        for i in 0..self.size() {
130            vec.push(self[i] * a);
131        }
132
133        V!(vec)
134    }
135}
136
137impl<K: Scalar + MulAssign<U>, U: Scalar> MulAssign<&U> for Vector<K> {
138    fn mul_assign(&mut self, a: &U) {
139        for i in &mut self._d {
140            *i *= *a;
141        }
142    }
143}
144
145pub fn linear_combination<K: Scalar>(
146    u: &[&Vector<K>],
147    coefs: &[K],
148) -> Vector<K> {
149    assert_eq!(
150        u.len(),
151        coefs.len(),
152        "vectors and scalers must be the same size"
153    );
154
155    assert!(
156        u.iter().all(|v| v.size() == u[0].size()),
157        "vectors must be have the same dimention"
158    );
159
160    let mut iter = u.iter().zip(coefs);
161
162    if let Some(mut sum) = iter.next().map(|(&v, &k)| v.clone() * k) {
163        for (v, k) in iter {
164            for i in 0..sum.size() {
165                sum._d[i] = v[i].mul_add(k, &sum[i]);
166            }
167        }
168        sum
169    } else {
170        V!()
171    }
172}
173
174pub fn cross_product<K: Scalar>(u: &Vector<K>, v: &Vector<K>) -> Vector<K> {
175    assert!(
176        u.size() == 3 && v.size() == u.size(),
177        "vectors must have be of size 3"
178    );
179
180    V!([
181        (u[1] * v[2]) - (u[2] * v[1]),
182        (u[2] * v[0]) - (u[0] * v[2]),
183        (u[0] * v[1]) - (u[1] * v[0]),
184    ])
185}
186
187impl<K: Scalar> Vector<K> {
188    pub fn zero(size: usize) -> Self {
189        V!(vec![K::default(); size])
190    }
191
192    pub fn size(&self) -> usize {
193        self._d.len()
194    }
195
196    pub fn first(&self) -> Option<&K> {
197        self._d.first()
198    }
199
200    pub fn is_empty(&self) -> bool {
201        self._d.is_empty()
202    }
203
204    pub fn add(&mut self, v: &Vector<K>) {
205        *self += v;
206    }
207
208    pub fn sub(&mut self, vec: &Vector<K>) {
209        *self -= vec;
210    }
211
212    pub fn scl(&mut self, a: K) {
213        *self *= &a;
214    }
215
216    pub fn norm_1(&self) -> K::AbsOutput {
217        let mut sum = K::AbsOutput::default();
218        for x in &self._d {
219            sum += x.abs();
220        }
221        sum
222    }
223
224    pub fn norm(&self) -> K::AbsOutput {
225        let mut sum = K::AbsOutput::default();
226        for x in &self._d {
227            let a = x.abs().clone();
228            sum = a.mul_add(&a, &sum);
229        }
230        sum.sqrt()
231    }
232
233    pub fn norm_inf(&self) -> K::AbsOutput {
234        let mut max = K::AbsOutput::default();
235        for x in &self._d {
236            let abs_x = x.abs();
237            if max < abs_x {
238                max = abs_x;
239            }
240        }
241        max
242    }
243}