maria_linalg/
vector.rs

1//! Implements necessary methods on vectors.
2
3use std::{
4    fmt,
5    ops::{
6        Add,
7        Sub,
8        Index,
9        IndexMut,
10    },
11};
12
13use rand::{
14    random,
15    thread_rng,
16    prelude::SliceRandom,
17};
18
19use rand_distr::{
20    Distribution,
21    Normal,
22};
23
24use super::Matrix;
25
26#[derive(Clone, Copy, PartialEq, Debug)]
27/// Abstracts over a vector of arbitary dimension.
28pub struct Vector<const N: usize> {
29    /// Contains the values of this vector.
30    pub values: [f64; N],
31}
32
33/// Implements necessary behaviors of a vector.
34impl<const N: usize> Vector<N> {
35    /// Constructs a zero vector.
36    pub fn zero() -> Self {
37        Self {
38            values: [0.0; N],
39        }
40    }
41
42    /// Constructs a vector of provided values.
43    pub fn new(values: [f64; N]) -> Self {
44        Self {
45            values,
46        }
47    }
48
49    /// Constructs an orthogonal basis of vectors.
50    pub fn basis() -> [Self; N] {
51        let i = Matrix::identity();
52
53        i.decompose()
54    }
55
56    /// Scales a vector by a provided scalar, returning the new vector.
57    pub fn scale(&self, scalar: f64) -> Self {
58        let mut newvalues = [0.0; N];
59        for i in 0..N {
60            newvalues[i] = scalar * self[i];
61        }
62
63        Self {
64            values: newvalues,
65        }
66    }
67    
68    /// Dots this vector with another vector.
69    pub fn dot(&self, other: Self) -> f64 {
70        let mut output = 0.0;
71
72        for i in 0..N {
73            output += self[i] * other[i];
74        }
75
76        output
77    }
78
79    /// Takes the norm of a vector.
80    pub fn norm(&self) -> f64 {
81        let mut output = 0.0;
82
83        for i in 0..N {
84            output += self[i].powf(2.0);
85        }
86
87        output.sqrt()
88    }
89
90    /// Returns the unit vector parallel to this vector.
91    pub fn normalize(&self) -> Self {
92        self.scale(1.0 / self.norm())
93    }
94
95    /// Left-multiplies the provided matrix by the transpose of this vector, returning the result.
96    pub fn mult(&self, matrix: Matrix<N>) -> Self {
97        let mut output = Self::zero();
98
99        for i in 0..N {
100            for j in 0..N {
101                output[i] += matrix[(i, j)] * self[j];
102            }
103        }
104
105        output
106    }
107
108    /// Given two vectors, generate a "child" vector.
109    /// This function is useful for genetic optimization algorithms.
110    pub fn child(mother: &Self, father: &Self, stdev: f64) -> Self {
111        let mut child = Self::zero();
112
113        for i in 0..N {
114            // Select gene for child
115            child[i] = if random::<f64>() < 0.5 {
116                mother[i]
117            } else {
118                father[i]
119            };
120
121            // Mutate this gene
122            // NOTE: it's ok to use `unwrap` here because we
123            // know that we will always be able to create a normal
124            // distribution of type N(0, `stdev`)
125            let normal = Normal::new(0.0, stdev).unwrap();
126            let v = normal.sample(&mut thread_rng());
127            child[i] += v;
128        }
129
130        child
131    }
132
133    /// Given two discrete vectors, generate a "child" discrete vector.
134    /// This function is useful for *discrete* genetic optimization algorithms.
135    pub fn child_discrete(mother: &Self, father: &Self, permitted: &[f64]) -> Self {
136        let mut child = Self::zero();
137
138        for i in 0..N {
139            // Select gene for child
140            child[i] = if random::<f64>() < 0.5 {
141                mother[i]
142            } else {
143                father[i]
144            };
145
146            if random::<f64>() < 1.0 / (N as f64) {
147                let mut rng = thread_rng();
148                let v = permitted.choose(&mut rng).unwrap();
149                child[i] = *v;
150            }
151        }
152
153        child
154    }
155
156    /// Determines if this vector is within the element-wise contraints.
157    pub fn check(&self, lower: [Option<f64>; N], upper: [Option<f64>; N]) -> bool {
158        for i in 0..N {
159            if let Some (l) = lower[i] {
160                if self[i] < l {
161                    return false;
162                }
163            }
164            
165            if let Some (u) = upper[i] {
166                if self[i] > u {
167                    return false;
168                }
169            }
170        }
171
172        true
173    }
174}
175
176impl Vector<3> {
177    /// Rotates this vector by the provided axis by the provided angle
178    ///     using Rodrigues' rotation formula.
179    /// 
180    /// *Note*: the provided angle is in radians.
181    pub fn rotate(&self, axis: Vector<3>, angle: f64) -> Self {
182        let k = axis.normalize();
183
184        self.scale(angle.cos())
185            + k.scale(self.dot(k) * (1.0 - angle.cos()))
186            + k.cross(*self).scale(angle.sin())
187    }
188
189    /// Crosses this vector with another vector.
190    pub fn cross(&self, other: Self) -> Self {
191        let x = self[1] * other[2] - self[2] * other[1];
192        let y = self[2] * other[0] - self[0] * other[2];
193        let z = self[0] * other[1] - self[1] * other[0];
194
195        [x, y, z].into()
196    }
197}
198
199impl<const N: usize> From<[f64; N]> for Vector<N> {
200    fn from(values: [f64; N]) -> Self {
201        Self::new(values)
202    }
203}
204
205impl<const N: usize> From<Vector<N>> for [f64; N] {
206    fn from(vector: Vector<N>) -> Self {
207        vector.values
208    }
209}
210
211impl<const N: usize> Index<usize> for Vector<N> {
212    type Output = f64;
213
214    fn index(&self, idx: usize) -> &Self::Output {
215        &self.values[idx]
216    }
217}
218
219impl<const N: usize> IndexMut<usize> for Vector<N> {
220    fn index_mut(&mut self, idx: usize) -> &mut Self::Output {
221        &mut self.values[idx]
222    }
223}
224
225impl<const N: usize> Add for Vector<N> {
226    type Output = Self;
227
228    fn add(self, other: Self) -> Self {
229        let mut new = Self::zero();
230        
231        for i in 0..self.values.len() {
232            new[i] = self[i] + other[i];
233        }
234
235        new
236    }
237}
238
239impl<const N: usize> Sub for Vector<N> {
240    type Output = Self;
241
242    fn sub(self, other: Self) -> Self {
243        let mut new = Self::zero();
244        
245        for i in 0..N {
246            new[i] = self[i] - other[i];
247        }
248
249        new
250    }
251}
252
253impl<const N: usize> fmt::Display for Vector<N> {
254    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
255        let mut rows = Vec::new();
256        let mut maxlen = 0;
257        for i in 0..N {
258            let value = self[i];
259            let row = if value >= 0.0 {
260                format!(" {:.8}", value)
261            } else {
262                format!("{:.8}", value)
263            };
264            let l = row.len();
265            rows.push(row);
266            if l > maxlen {
267                maxlen = l;
268            }
269        }
270
271        let mut output = String::new();
272        for i in 0..N {
273            output.push_str("[");
274            output.push_str(
275                &format!("{:^i$}", rows[i], i = maxlen + 2)
276            );
277            output.push_str("]\n");
278        }
279
280        write!(f, "{}", output)
281    }
282}
283
284#[test]
285fn display_vector() {
286    let vec: Vector<3> = [-0.15, 10.0, 1000.0].into();
287    println!("{}", vec);
288}
maria_linalg/vector.rs

maria_linalg/
vector.rs
