hoomd_vector/
cartesian.rs

1// Copyright (c) 2024-2026 The Regents of the University of Michigan.
2// Part of hoomd-rs, released under the BSD 3-Clause License.
3
4//! Implement canonical vector types.
5
6use serde::{Deserialize, Serialize};
7use serde_with::serde_as;
8use std::{
9    array, fmt,
10    iter::{Sum, zip},
11    ops::{Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign},
12};
13
14use approxim::approx_derive::RelativeEq;
15use rand::{
16    Rng,
17    distr::{Distribution, StandardUniform, Uniform},
18};
19
20use crate::{Cross, Error, InnerProduct, Metric, Rotate, Unit, Vector};
21use hoomd_linear_algebra::{MatMul, matrix::Matrix};
22
23/// A [`Vector`] represented by `N` `f64` coordinates.
24///
25/// [`Cartesian`] is the canonical implementation of [`InnerProduct`].
26///
27/// ## Constructing vectors
28///
29/// The default is the 0 vector:
30/// ```
31/// use hoomd_vector::Cartesian;
32///
33/// let v = Cartesian::<3>::default();
34/// assert_eq!(v, [0.0; 3].into())
35/// ```
36///
37/// Create a vector with an array of coordinates:
38/// ```
39/// use hoomd_vector::Cartesian;
40///
41/// let v = Cartesian::from([1.0, 2.0, 3.0, 4.0, 5.0]);
42/// ```
43///
44/// 2D and 3D vectors can also be initialized from tuples:
45/// ```
46/// use hoomd_vector::Cartesian;
47///
48/// let a = Cartesian::from((1.0, 2.0, 3.0));
49/// let b = Cartesian::from((4.0, 5.0));
50/// ```
51///
52/// Construct a random vector in the [-1, 1] hypercube:
53///
54/// ```
55/// use hoomd_vector::Cartesian;
56/// use rand::{RngExt, SeedableRng, rngs::StdRng};
57///
58/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
59/// let mut rng = StdRng::seed_from_u64(1);
60/// let v: Cartesian<3> = rng.random();
61/// # Ok(())
62/// # }
63/// ```
64///
65/// ## Operating on vectors
66///
67/// Use vector math operations when you can:
68/// ```
69/// use hoomd_vector::{Cartesian, InnerProduct};
70///
71/// let a = Cartesian::from([1.0, 2.0]);
72/// let b = Cartesian::from([4.0, 8.0]);
73/// let c = (a + b).dot(&a);
74/// ```
75///
76/// Access the coordinates directly when needed:
77/// ```
78/// use hoomd_vector::Cartesian;
79///
80/// let a = Cartesian::from((1.0, 2.0));
81/// let b = Cartesian::from((a[1], 0.0));
82/// ```
83///
84/// Compute the sum of an iterator over vectors:
85/// ```
86/// use hoomd_vector::Cartesian;
87///
88/// let total: Cartesian<2> =
89///     [Cartesian::from((1.0, 2.0)), Cartesian::from((3.0, 4.0))]
90///         .into_iter()
91///         .sum();
92/// ```
93#[serde_as]
94#[derive(Clone, Copy, Debug, PartialEq, RelativeEq, Serialize, Deserialize)]
95#[approx(epsilon_type = f64)]
96pub struct Cartesian<const N: usize> {
97    /// The vector's coordinates.
98    #[serde_as(as = "[_; N]")]
99    #[approx(into_iter)]
100    pub coordinates: [f64; N],
101}
102
103impl<const N: usize> Default for Cartesian<N> {
104    /// Create a 0 vector.
105    ///
106    /// # Example
107    /// ```
108    /// use hoomd_vector::Cartesian;
109    ///
110    /// let v = Cartesian::<3>::default();
111    /// assert_eq!(v, [0.0; 3].into())
112    /// ```
113    #[inline]
114    fn default() -> Self {
115        Cartesian::from([0.0; N])
116    }
117}
118
119impl<const N: usize> From<[f64; N]> for Cartesian<N> {
120    /// Create a Cartesian vector with the given coordinates.
121    ///
122    /// # Example
123    /// ```
124    /// use hoomd_vector::Cartesian;
125    ///
126    /// let v = Cartesian::from([4.0, 3.0]);
127    /// ```
128    #[inline]
129    fn from(coordinates: [f64; N]) -> Self {
130        Self { coordinates }
131    }
132}
133
134impl<const N: usize> IntoIterator for Cartesian<N> {
135    type Item = f64;
136    type IntoIter = <[f64; N] as IntoIterator>::IntoIter;
137
138    /// Iterate over the components of the vector.
139    ///
140    /// # Example
141    /// ```
142    /// use hoomd_vector::Cartesian;
143    ///
144    /// let a = Cartesian::from([1.0, 2.0]);
145    /// let mut iter = a.into_iter();
146    ///
147    /// assert_eq!(iter.next(), Some(1.0));
148    /// assert_eq!(iter.next(), Some(2.0));
149    /// assert_eq!(iter.next(), None);
150    /// ```
151    #[inline]
152    fn into_iter(self) -> Self::IntoIter {
153        self.coordinates.into_iter()
154    }
155}
156
157impl From<(f64, f64)> for Cartesian<2> {
158    #[inline]
159    fn from(coordinates: (f64, f64)) -> Self {
160        Self {
161            coordinates: coordinates.into(),
162        }
163    }
164}
165
166impl From<(f64, f64, f64)> for Cartesian<3> {
167    #[inline]
168    fn from(coordinates: (f64, f64, f64)) -> Self {
169        Self {
170            coordinates: coordinates.into(),
171        }
172    }
173}
174
175impl<const N: usize> TryFrom<Vec<f64>> for Cartesian<N> {
176    type Error = Error;
177
178    /// Create a Cartesian vector with coordinates given by a [`Vec<f64>`]
179    ///
180    /// # Example
181    /// ```
182    /// use hoomd_vector::Cartesian;
183    ///
184    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
185    /// let v = Cartesian::<3>::try_from(vec![3.0, 4.0, 5.0])?;
186    /// assert_eq!(v, [3.0, 4.0, 5.0].into());
187    /// # Ok(())
188    /// # }
189    /// ```
190    /// <div class="warning">
191    ///
192    /// This method deallocates the Vec after copying it.
193    /// Use `Cartesian::From<[f64; N]>` in performance critical code.
194    ///
195    /// </div>
196    #[inline]
197    fn try_from(value: Vec<f64>) -> Result<Self, Self::Error> {
198        let coordinates = value.try_into().map_err(|_| Error::InvalidVectorLength)?;
199        Ok(Self { coordinates })
200    }
201}
202
203impl<const N: usize> TryFrom<std::ops::Range<usize>> for Cartesian<N> {
204    type Error = Error;
205
206    /// Create a Cartesian vector with coordinates given by a range.
207    ///
208    /// # Example
209    /// ```
210    /// use hoomd_vector::Cartesian;
211    ///
212    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
213    /// let v = Cartesian::<3>::try_from(3..6)?;
214    /// assert_eq!(v, [3.0, 4.0, 5.0].into());
215    /// # Ok(())
216    /// # }
217    /// ```
218    #[inline]
219    fn try_from(value: std::ops::Range<usize>) -> Result<Self, Self::Error> {
220        if value.len() != N {
221            return Err(Error::InvalidVectorLength);
222        }
223
224        let coordinates = array::from_fn(|i| (value.start + i) as f64);
225        Ok(Self { coordinates })
226    }
227}
228
229impl<const N: usize> TryFrom<[f64; N]> for Unit<Cartesian<N>> {
230    type Error = Error;
231
232    /// Create a unit Cartesian vector by normalizing the given coordinates.
233    ///
234    /// # Example
235    /// ```
236    /// use hoomd_vector::{Cartesian, Unit};
237    ///
238    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
239    /// let n = Unit::<Cartesian<3>>::try_from([2.0, 0.0, 0.0])?;
240    /// assert_eq!(*n.get(), [1.0, 0.0, 0.0].into());
241    /// # Ok(())
242    /// # }
243    /// ```
244    #[inline]
245    fn try_from(value: [f64; N]) -> Result<Self, Self::Error> {
246        Cartesian::from(value).to_unit().map(|t| t.0)
247    }
248}
249
250impl<const N: usize> Vector for Cartesian<N> {}
251
252impl<const N: usize> InnerProduct for Cartesian<N> {
253    #[inline]
254    fn dot(&self, other: &Self) -> f64 {
255        zip(self.coordinates.iter(), other.coordinates.iter())
256            .fold(0.0, |product, (&x, &y)| x.mul_add(y, product))
257    }
258
259    /// Create a unit vector in the space.
260    ///
261    /// The default unit vector in Cartesian space is `[0.0, 0.0, ...., 1.0]`.
262    ///
263    /// # Example
264    /// ```
265    /// use hoomd_vector::{Cartesian, InnerProduct};
266    ///
267    /// let u = Cartesian::<2>::default_unit();
268    /// assert_eq!(*u.get(), [0.0, 1.0].into());
269    ///
270    /// let u = Cartesian::<3>::default_unit();
271    /// assert_eq!(*u.get(), [0.0, 0.0, 1.0].into());
272    /// ```
273    #[inline]
274    fn default_unit() -> Unit<Self> {
275        assert!(N >= 1);
276        let mut coordinates = [0.0; N];
277        coordinates[N - 1] = 1.0;
278        Unit(Self { coordinates })
279    }
280}
281
282impl<const N: usize> Metric for Cartesian<N> {
283    /// Computes the squared distance between two points in Euclidean space.
284    /// ```math
285    /// d^2(\vec{x},\vec{y}) = \sum_{i=1}^{N} (x_i - y_i)^2
286    /// ```
287    ///
288    /// # Example
289    /// ```
290    /// use hoomd_vector::{Cartesian, Metric};
291    ///
292    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
293    /// let x = Cartesian::from([0.0, 1.0, 1.0]);
294    /// let y = Cartesian::from([1.0, 0.0, 0.0]);
295    /// assert_eq!(3.0, x.distance_squared(&y));
296    /// # Ok(())
297    /// # }
298    /// ```
299    #[inline]
300    fn distance_squared(&self, other: &Self) -> f64 {
301        zip(self.coordinates.iter(), other.coordinates.iter())
302            .fold(0.0, |product, x| product + (x.0 - x.1).powi(2))
303    }
304    #[inline]
305    fn distance(&self, other: &Self) -> f64 {
306        (self.distance_squared(other)).sqrt()
307    }
308    /// Return the number of dimensions in this Cartesian vector space.
309    ///
310    /// # Example
311    /// ```
312    /// use hoomd_vector::{Cartesian, Metric};
313    ///
314    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
315    /// let vec2 = Cartesian::<2>::default();
316    /// let vec3 = Cartesian::<3>::default();
317    /// assert_eq!(2, vec2.n_dimensions());
318    /// assert_eq!(3, vec3.n_dimensions());
319    /// # Ok(())
320    /// # }
321    /// ```
322    #[inline]
323    fn n_dimensions(&self) -> usize {
324        N
325    }
326}
327
328impl<const N: usize> Add for Cartesian<N> {
329    type Output = Self;
330
331    #[inline]
332    fn add(self, rhs: Self) -> Self {
333        let mut coordinates = [0.0; N];
334
335        for (result, (a, b)) in coordinates
336            .iter_mut()
337            .zip(self.coordinates.iter().zip(rhs.coordinates.iter()))
338        {
339            *result = a + b;
340        }
341        Self { coordinates }
342    }
343}
344
345impl<const N: usize> AddAssign for Cartesian<N> {
346    #[inline]
347    fn add_assign(&mut self, rhs: Self) {
348        for (result, a) in self.coordinates.iter_mut().zip(rhs.coordinates.iter()) {
349            *result += a;
350        }
351    }
352}
353
354impl<const N: usize> Div<f64> for Cartesian<N> {
355    type Output = Self;
356
357    #[inline]
358    fn div(self, rhs: f64) -> Self {
359        Self {
360            coordinates: self.coordinates.map(|x| x / rhs),
361        }
362    }
363}
364
365impl<const N: usize> DivAssign<f64> for Cartesian<N> {
366    #[inline]
367    fn div_assign(&mut self, rhs: f64) {
368        self.coordinates = self.coordinates.map(|x| x / rhs);
369    }
370}
371
372impl<const N: usize> Mul<f64> for Cartesian<N> {
373    type Output = Self;
374
375    #[inline]
376    fn mul(self, rhs: f64) -> Self {
377        Self {
378            coordinates: self.coordinates.map(|x| x * rhs),
379        }
380    }
381}
382
383impl<const N: usize> MulAssign<f64> for Cartesian<N> {
384    #[inline]
385    fn mul_assign(&mut self, rhs: f64) {
386        self.coordinates = self.coordinates.map(|x| x * rhs);
387    }
388}
389
390impl<const N: usize> Sub for Cartesian<N> {
391    type Output = Self;
392
393    #[inline]
394    fn sub(self, rhs: Self) -> Self {
395        let mut coordinates = [0.0; N];
396        for (result, (a, b)) in coordinates
397            .iter_mut()
398            .zip(self.coordinates.iter().zip(rhs.coordinates.iter()))
399        {
400            *result = a - b;
401        }
402        Self { coordinates }
403    }
404}
405
406impl<const N: usize> SubAssign for Cartesian<N> {
407    #[inline]
408    fn sub_assign(&mut self, rhs: Self) {
409        for (result, a) in self.coordinates.iter_mut().zip(rhs.coordinates) {
410            *result -= a;
411        }
412    }
413}
414
415impl<const N: usize> fmt::Display for Cartesian<N> {
416    #[inline]
417    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
418        write!(
419            f,
420            "[{}]",
421            self.coordinates
422                .iter()
423                .map(f64::to_string)
424                .collect::<Vec<String>>()
425                .join(", ")
426        )
427    }
428}
429
430impl<const N: usize> Neg for Cartesian<N> {
431    type Output = Self;
432    #[inline]
433    fn neg(self) -> Self::Output {
434        Self {
435            coordinates: self.coordinates.map(|x| -x),
436        }
437    }
438}
439
440impl Cross for Cartesian<3> {
441    #[inline]
442    fn cross(&self, other: &Self) -> Self {
443        Cartesian::from((
444            self[1] * other[2] - self[2] * other[1],
445            self[2] * other[0] - self[0] * other[2],
446            self[0] * other[1] - self[1] * other[0],
447        ))
448    }
449}
450
451impl<const N: usize> Distribution<Cartesian<N>> for StandardUniform {
452    /// Sample a Cartesian vector from the uniform [-1, 1] hypercube.
453    ///
454    /// Each coordinate in the vector is in the closed range [-1, 1].
455    ///
456    /// # Example
457    ///
458    /// ```
459    /// use hoomd_vector::Cartesian;
460    /// use rand::{RngExt, SeedableRng, rngs::StdRng};
461    ///
462    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
463    /// let mut rng = StdRng::seed_from_u64(1);
464    /// let v: Cartesian<3> = rng.random();
465    /// # Ok(())
466    /// # }
467    /// ```
468    #[inline]
469    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> Cartesian<N> {
470        let uniform = Uniform::new_inclusive(-1.0, 1.0)
471            .expect("hard-coded range should form a valid distribution");
472        Cartesian {
473            coordinates: array::from_fn(|_| uniform.sample(rng)),
474        }
475    }
476}
477
478impl<const N: usize, T> Index<T> for Cartesian<N>
479where
480    T: Into<usize> + std::slice::SliceIndex<[f64], Output = f64>,
481{
482    type Output = f64;
483    /// Get the value of the vector at coordinate i.
484    ///
485    /// # Example
486    /// ```
487    /// use hoomd_vector::Cartesian;
488    ///
489    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
490    /// let v = Cartesian::<3>::try_from(3..6)?;
491    /// assert_eq!((v[0], v[1], v[2]), (3.0, 4.0, 5.0));
492    /// # Ok(())
493    /// # }
494    /// ```
495    #[inline]
496    fn index(&self, index: T) -> &Self::Output {
497        &self.coordinates[index]
498    }
499}
500
501impl Cartesian<2> {
502    /// Construct a 2-vector perpendicular to self.
503    ///
504    /// Given a vector $`(v_x, v_y)`$ `perpendicular` returns the vector
505    /// rotated by $`\pi/2`$:
506    /// ```math
507    /// (-v_y, v_x)
508    /// ```
509    ///
510    /// # Example
511    /// ```
512    /// use hoomd_vector::Cartesian;
513    ///
514    /// let v = Cartesian::from([1.0, -4.5]);
515    /// assert_eq!(v.perpendicular(), [4.5, 1.0].into());
516    /// ```
517    #[inline]
518    #[must_use]
519    pub fn perpendicular(self) -> Self {
520        Cartesian::from([-self[1], self[0]])
521    }
522}
523
524impl<const N: usize, T> IndexMut<T> for Cartesian<N>
525where
526    T: Into<usize> + std::slice::SliceIndex<[f64], Output = f64>,
527{
528    /// Get a mutable reference to the value of the vector at coordinate i.
529    ///
530    /// # Example
531    /// ```
532    /// use hoomd_vector::Cartesian;
533    ///
534    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
535    /// let mut v = Cartesian::<3>::try_from(3..6)?;
536    /// assert_eq!((v[0], v[1], v[2]), (3.0, 4.0, 5.0));
537    /// v[0] += 1.0;
538    /// assert_eq!(v[0], 4.0);
539    /// # Ok(())
540    /// # }
541    /// ```
542    #[inline]
543    fn index_mut(&mut self, index: T) -> &mut Self::Output {
544        &mut self.coordinates[index]
545    }
546}
547
548impl<const N: usize> Sum for Cartesian<N> {
549    #[inline]
550    fn sum<I>(iter: I) -> Self
551    where
552        I: Iterator<Item = Self>,
553    {
554        iter.fold(Cartesian::default(), |acc, x| acc + x)
555    }
556}
557
558/// Rotate vectors efficiently.
559///
560/// Construct a [`RotationMatrix`] to efficiently rotate many vectors by the same rotation.
561///
562/// See:
563/// * [`RotationMatrix::from<Angle>`]
564/// * [`RotationMatrix::from<Versor>`]
565///
566/// [`RotationMatrix`] _intentionally_ does not implement [`Rotation`](crate::Rotation).
567/// [`Angle`](crate::Angle) and [`Versor`](crate::Versor) are representations of
568/// rotations that are often the most effective and numerically stable to
569/// manipulate.
570#[serde_as]
571#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
572pub struct RotationMatrix<const N: usize> {
573    /// Rows of the rotation matrix.
574    #[serde_as(as = "[_; N]")]
575    pub(crate) rows: [Cartesian<N>; N],
576}
577
578impl<const N: usize> From<RotationMatrix<N>> for Matrix<N, N> {
579    #[inline]
580    fn from(value: RotationMatrix<N>) -> Self {
581        Self {
582            rows: value.rows().map(|arr| arr.coordinates),
583        }
584    }
585}
586
587impl<const N: usize> RotationMatrix<N> {
588    /// Get the rows of the rotation matrix.
589    ///
590    /// # Example
591    /// ```
592    /// use hoomd_vector::{Angle, InnerProduct, RotationMatrix, Vector};
593    /// use std::f64::consts::PI;
594    ///
595    /// let a = Angle::from(PI / 2.0);
596    ///
597    /// let matrix = RotationMatrix::from(a);
598    /// assert!(matrix.rows()[0].dot(&[0.0, -1.0].into()) > 0.99);
599    /// assert!(matrix.rows()[1].dot(&[1.0, 0.0].into()) > 0.99);
600    /// ```
601    #[inline]
602    #[must_use]
603    pub fn rows(&self) -> [Cartesian<N>; N] {
604        self.rows
605    }
606
607    /// Create a matrix that performs the inverse rotation.
608    ///
609    /// Matrix inversion is cheaper than [`Angle`](crate::Angle) ->
610    /// [`RotationMatrix`] and [`Versor`](crate::Versor) -> [`RotationMatrix`]
611    /// conversions. When you need both rotations, convert once and then invert.
612    ///
613    /// # Example
614    /// ```
615    /// use hoomd_vector::{Angle, InnerProduct, RotationMatrix, Vector};
616    /// use std::f64::consts::PI;
617    ///
618    /// let a = Angle::from(PI / 2.0);
619    ///
620    /// let matrix = RotationMatrix::from(a);
621    /// let inverted_matrix = matrix.inverted();
622    /// assert!(inverted_matrix.rows()[0].dot(&[0.0, 1.0].into()) > 0.99);
623    /// assert!(inverted_matrix.rows()[1].dot(&[-1.0, 0.0].into()) > 0.99);
624    /// ```
625    #[inline]
626    #[must_use]
627    pub fn inverted(self) -> Self {
628        Self {
629            rows: array::from_fn(|i| array::from_fn::<_, N, _>(|j| self.rows()[j][i]).into()),
630        }
631    }
632}
633
634impl<const N: usize> Default for RotationMatrix<N> {
635    /// Create an identity matrix.
636    ///
637    /// ```math
638    /// \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}
639    /// ```
640    /// ,
641    /// ```math
642    /// \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}
643    /// ```
644    /// , and so on.
645    ///
646    /// # Example
647    ///
648    /// ```
649    /// use hoomd_vector::RotationMatrix;
650    ///
651    /// let identity = RotationMatrix::<3>::default();
652    /// ```
653    #[inline]
654    fn default() -> RotationMatrix<N> {
655        RotationMatrix {
656            rows: array::from_fn(|i| array::from_fn(|j| if i == j { 1.0 } else { 0.0 }).into()),
657        }
658    }
659}
660
661impl<const N: usize> fmt::Display for RotationMatrix<N> {
662    #[inline]
663    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
664        write!(
665            f,
666            "[{}]",
667            self.rows
668                .iter()
669                .map(Cartesian::<N>::to_string)
670                .collect::<Vec<String>>()
671                .join("\n ")
672        )
673    }
674}
675
676impl<const N: usize> Rotate<Cartesian<N>> for RotationMatrix<N> {
677    type Matrix = RotationMatrix<N>;
678
679    #[inline]
680    /// Rotate a [`Cartesian<N>`] by a [`RotationMatrix`]
681    ///
682    /// # Examples
683    /// ```
684    /// use approxim::assert_relative_eq;
685    /// use hoomd_vector::{Angle, Cartesian, Rotate, RotationMatrix};
686    /// use std::f64::consts::PI;
687    ///
688    /// let v = Cartesian::from([-1.0, 0.0]);
689    /// let a = Angle::from(PI / 2.0);
690    ///
691    /// let matrix = RotationMatrix::from(a);
692    /// let rotated = matrix.rotate(&v);
693    /// assert_relative_eq!(rotated, [0.0, -1.0].into());
694    /// ```
695    ///
696    /// ```
697    /// use approxim::assert_relative_eq;
698    /// use hoomd_vector::{Cartesian, Rotate, RotationMatrix, Versor};
699    /// use std::f64::consts::PI;
700    ///
701    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
702    /// let a = Cartesian::from([-1.0, 0.0, 0.0]);
703    /// let v = Versor::from_axis_angle([0.0, 0.0, 1.0].try_into()?, PI / 2.0);
704    ///
705    /// let matrix = RotationMatrix::from(v);
706    /// let b = matrix.rotate(&a);
707    /// assert_relative_eq!(b, [0.0, -1.0, 0.0].into());
708    /// # Ok(())
709    /// # }
710    /// ```
711    fn rotate(&self, vector: &Cartesian<N>) -> Cartesian<N> {
712        let mut coordinates = [0.0; N];
713
714        for (result, row) in coordinates.iter_mut().zip(self.rows.iter()) {
715            *result = row.dot(vector);
716        }
717
718        Cartesian { coordinates }
719    }
720}
721
722impl<const N: usize, const K: usize> MatMul<Matrix<N, K>> for RotationMatrix<N> {
723    type Output = Matrix<N, K>;
724
725    #[inline]
726    fn matmul(&self, rhs: &Matrix<N, K>) -> Self::Output {
727        Matrix::from(*self).matmul(rhs)
728    }
729}
730
731impl<const N: usize> Cartesian<N> {
732    /// Convert a [`Cartesian<N>`] into a row matrix [`Matrix<1, N>`].
733    ///
734    /// # Example
735    /// ```
736    /// use hoomd_linear_algebra::matrix::Matrix;
737    /// use hoomd_vector::Cartesian;
738    ///
739    /// let a = Cartesian::from([1.0, -2.0, 3.0]);
740    ///
741    /// let b = a.to_row_matrix();
742    /// assert_eq!(b.rows, [[1.0, -2.0, 3.0]]);
743    /// ```
744    #[inline]
745    #[must_use]
746    pub fn to_row_matrix(self) -> Matrix<1, N> {
747        Matrix {
748            rows: [self.coordinates],
749        }
750    }
751    /// Convert a [`Cartesian<N>`] into a column matrix [`Matrix<N, 1>`].
752    ///
753    /// # Example
754    /// ```
755    /// use hoomd_linear_algebra::matrix::Matrix;
756    /// use hoomd_vector::Cartesian;
757    ///
758    /// let a = Cartesian::from([1.0, -2.0, 3.0]);
759    ///
760    /// let b = a.to_column_matrix();
761    /// assert_eq!(b.rows, [[1.0], [-2.0], [3.0]]);
762    /// ```
763    #[inline]
764    #[must_use]
765    pub fn to_column_matrix(self) -> Matrix<N, 1> {
766        Matrix {
767            rows: std::array::from_fn(|i| [self[i]]),
768        }
769    }
770}
771
772impl<const N: usize> From<Matrix<1, N>> for Cartesian<N> {
773    #[inline]
774    fn from(value: Matrix<1, N>) -> Self {
775        value.rows[0].into()
776    }
777}
778
779#[cfg(test)]
780mod tests {
781    use crate::{Angle, Rotation, Versor};
782
783    use super::*;
784    use approxim::assert_relative_eq;
785    use paste::paste;
786    use rand::{RngExt, SeedableRng, rngs::StdRng};
787    use rstest::rstest;
788
789    // Parameterize a test function over an array of vector lengths
790    macro_rules! parameterize_vector_length {
791        // macro with name as above that takes an identifier (fn) and an expression
792        // $(...),* matches 0 or more expressions (values) separated by commas
793        ($test_body:ident, [$($dim:expr),*]) => {
794
795            // Now, we repeat the test block 0 or more times, one for each $dim
796            $(
797                // paste package combines values in [< >] to form a new ident
798                paste! {
799                    #[test]
800                    fn [< $test_body "_" $dim>]() {
801                        const DIM: usize = $dim;
802                        $test_body::<DIM>();
803                    }
804                }
805            )*
806        };
807    }
808
809    /// Generate a pair of length N vectors.
810    /// The first vector ranges from [0, N-1] and the second ranges from [N, 2*N-1]
811    fn generate_vector_pair<const N: usize>() -> (Cartesian<N>, Cartesian<N>) {
812        (
813            Cartesian::try_from(0..N).unwrap(),
814            Cartesian::try_from(N..N * 2).unwrap(),
815        )
816    }
817
818    fn index<const N: usize>() {
819        let (_, b) = generate_vector_pair::<N>();
820        assert!(zip(0..N, b.coordinates.iter()).all(|(i, &x)| b[i] == x));
821    }
822    parameterize_vector_length!(index, [2, 3, 4, 8, 16, 32]);
823
824    fn index_mut<const N: usize>() {
825        let (a, mut b) = generate_vector_pair::<N>();
826        zip(0..N, b.coordinates.iter_mut()).for_each(|(i, x)| *x = a[i]);
827        assert_eq!(a, b);
828    }
829    parameterize_vector_length!(index_mut, [2, 3, 4, 8, 16, 32]);
830
831    fn add_explicit<const N: usize>() {
832        let (a, b) = generate_vector_pair::<N>();
833        let c = a.add(b);
834
835        let addition_answer: Vec<f64> = (0..(2 * N))
836            .step_by(2)
837            .map(|x| (x + N) as f64)
838            .collect::<Vec<_>>();
839
840        assert_eq!(c, Cartesian::try_from(addition_answer).unwrap());
841    }
842    parameterize_vector_length!(add_explicit, [2, 3, 4, 8, 16, 32]);
843
844    fn add_operator<const N: usize>() {
845        let (a, b) = generate_vector_pair::<N>();
846        let c = a + b;
847
848        let addition_answer: Vec<f64> = (0..(2 * N))
849            .step_by(2)
850            .map(|x| (x + N) as f64)
851            .collect::<Vec<_>>();
852
853        assert_eq!(c, Cartesian::try_from(addition_answer).unwrap());
854    }
855    parameterize_vector_length!(add_operator, [2, 3, 4, 8, 16, 32]);
856
857    fn add_assign<const N: usize>() {
858        let (a, b) = generate_vector_pair::<N>();
859        let mut c = a;
860        c += b;
861
862        let addition_answer: Vec<f64> = (0..(2 * N))
863            .step_by(2)
864            .map(|x| (x + N) as f64)
865            .collect::<Vec<_>>();
866
867        assert_eq!(c, Cartesian::try_from(addition_answer).unwrap());
868    }
869    parameterize_vector_length!(add_assign, [2, 3, 4, 8, 16, 32]);
870
871    fn sub_operator<const N: usize>() {
872        let (a, b) = generate_vector_pair::<N>();
873        let c = a - b;
874
875        let subtraction_answer = [-(N as f64); N];
876
877        assert_eq!(c, subtraction_answer.into());
878    }
879    parameterize_vector_length!(sub_operator, [2, 3, 4, 8, 16, 32]);
880
881    fn sub_assign<const N: usize>() {
882        let (a, b) = generate_vector_pair::<N>();
883        let mut c = a;
884        c -= b;
885
886        let subtraction_answer = [-(N as f64); N];
887
888        assert_eq!(c, subtraction_answer.into());
889    }
890
891    parameterize_vector_length!(sub_assign, [2, 3, 4, 8, 16, 32]);
892
893    fn mul_operator<const N: usize>() {
894        let (a, _) = generate_vector_pair::<N>();
895        let b = 12.0;
896        let c = a * b;
897
898        let multiplication_answer: Vec<f64> = (0..N).map(|x| (x as f64) * b).collect::<Vec<_>>();
899
900        assert_eq!(c, Cartesian::try_from(multiplication_answer).unwrap());
901    }
902    parameterize_vector_length!(mul_operator, [2, 3, 4, 8, 16, 32]);
903
904    fn mul_assign<const N: usize>() {
905        let (mut a, _) = generate_vector_pair::<N>();
906        let b = 12.0;
907        a *= b;
908
909        let multiplication_answer: Vec<f64> = (0..N).map(|x| (x as f64) * b).collect::<Vec<_>>();
910
911        assert_eq!(a, Cartesian::try_from(multiplication_answer).unwrap());
912    }
913
914    parameterize_vector_length!(mul_assign, [2, 3, 4, 8, 16, 32]);
915
916    fn div_operator<const N: usize>() {
917        let (a, _) = generate_vector_pair::<N>();
918        let b = 12.0;
919        let c = a / b;
920
921        let division_answer: Vec<f64> = (0..N).map(|x| (x as f64) / b).collect::<Vec<_>>();
922
923        assert_eq!(c, Cartesian::try_from(division_answer).unwrap());
924    }
925    parameterize_vector_length!(div_operator, [2, 3, 4, 8, 16, 32]);
926
927    fn div_assign<const N: usize>() {
928        let (mut a, _) = generate_vector_pair::<N>();
929        let b = 12.0;
930        a /= b;
931
932        let division_answer: Vec<f64> = (0..N).map(|x| (x as f64) / b).collect::<Vec<_>>();
933
934        assert_eq!(a, Cartesian::try_from(division_answer).unwrap());
935    }
936
937    parameterize_vector_length!(div_assign, [2, 3, 4, 8, 16, 32]);
938
939    fn compute_add_ref_ref<const N: usize>(a: &Cartesian<N>, b: &Cartesian<N>) -> Cartesian<N> {
940        *a + *b
941    }
942
943    fn compute_add_ref_type<const N: usize>(a: &Cartesian<N>, b: Cartesian<N>) -> Cartesian<N> {
944        *a + b
945    }
946
947    fn compute_add_type_ref<const N: usize>(a: Cartesian<N>, b: &Cartesian<N>) -> Cartesian<N> {
948        a + *b
949    }
950
951    fn add_with_refs<const N: usize>() {
952        let (a, b) = generate_vector_pair::<N>();
953
954        let addition_answer = Cartesian::try_from(
955            (0..(2 * N))
956                .step_by(2)
957                .map(|x| (x + N) as f64)
958                .collect::<Vec<_>>(),
959        )
960        .unwrap();
961
962        let c = compute_add_ref_ref(&a, &b);
963        assert_eq!(c, addition_answer);
964
965        let c = compute_add_ref_type(&a, b);
966        assert_eq!(c, addition_answer);
967
968        let c = compute_add_type_ref(a, &b);
969        assert_eq!(c, addition_answer);
970    }
971    parameterize_vector_length!(add_with_refs, [2, 3, 4, 8, 16, 32]);
972
973    fn dot<const N: usize>() {
974        let (a, b) = generate_vector_pair::<N>();
975        let c = a.dot(&b);
976
977        let n = N as f64;
978
979        // Analytical solution to sum of k * (n+k), k = 0 to n-1
980        let dot_ans = (5.0 * n.powi(3) - 6.0 * n.powi(2) + n) / 6.0;
981
982        assert_eq!(c, dot_ans);
983    }
984    parameterize_vector_length!(dot, [2, 3, 4, 8, 16, 32]);
985
986    fn neg<const N: usize>() {
987        let a: Cartesian<N> = Cartesian::try_from(0..N).unwrap();
988        let b = -a;
989        for (i, j) in zip(a.coordinates.iter(), b.coordinates.iter()) {
990            assert_eq!(*i, -j);
991        }
992    }
993    parameterize_vector_length!(neg, [2, 3, 4, 8, 16, 32]);
994
995    #[test]
996    fn cross() {
997        let (a, b) = generate_vector_pair::<3>();
998        let c = a.cross(&b);
999
1000        // Analytical solution
1001        let cross_ans = Cartesian::from((-3.0, 6.0, -3.0));
1002
1003        assert_eq!(c, cross_ans);
1004
1005        let a = Cartesian::from([1.0, 0.0, 0.0]);
1006        let b = Cartesian::from([0.0, 1.0, 0.0]);
1007        assert_eq!(a.cross(&b), [0.0, 0.0, 1.0].into());
1008
1009        let a = Cartesian::from([0.0, 3.0, 0.0]);
1010        let b = Cartesian::from([0.0, 0.0, 2.0]);
1011        assert_eq!(a.cross(&b), [6.0, 0.0, 0.0].into());
1012
1013        let a = Cartesian::from([2.0, 0.0, 0.0]);
1014        let b = Cartesian::from([0.0, 0.0, 4.0]);
1015        assert_eq!(a.cross(&b), [0.0, -8.0, 0.0].into());
1016    }
1017
1018    #[test]
1019    fn display() {
1020        let a = Cartesian::from([1.5, 2.125, -3.875]);
1021        let s = format!("{a}");
1022        assert_eq!(s, "[1.5, 2.125, -3.875]");
1023
1024        let a = Cartesian::from([10.0, 20.0, 30.0, 40.0]);
1025        let s = format!("{a}");
1026        assert_eq!(s, "[10, 20, 30, 40]");
1027    }
1028
1029    #[test]
1030    fn from_2_tuple() {
1031        let a = Cartesian::from((3.0, 0.125));
1032        assert_eq!(a.coordinates, [3.0, 0.125]);
1033    }
1034
1035    #[test]
1036    fn from_3_tuple() {
1037        let a = Cartesian::from((-0.5, 2.0, 18.125));
1038        assert_eq!(a.coordinates, [-0.5, 2.0, 18.125]);
1039    }
1040
1041    fn from_vec<const N: usize>() {
1042        let mut vec = Vec::with_capacity(N);
1043
1044        assert_eq!(
1045            Cartesian::<N>::try_from(vec.clone()),
1046            Err(Error::InvalidVectorLength)
1047        );
1048
1049        for i in 0..N {
1050            vec.push(i as f64 * 0.5);
1051        }
1052        let a = Cartesian::<N>::try_from(vec.clone()).unwrap();
1053
1054        assert_eq!(vec, Vec::from(a.coordinates));
1055
1056        vec.push(1.0);
1057        assert_eq!(
1058            Cartesian::<N>::try_from(vec.clone()),
1059            Err(Error::InvalidVectorLength)
1060        );
1061    }
1062    parameterize_vector_length!(from_vec, [2, 3, 4, 8, 16, 32]);
1063
1064    fn random_in_range<const N: usize>() {
1065        // Loosely verify we are drawing from the correct distribution
1066        let mut rng = StdRng::seed_from_u64(1);
1067        let a: Cartesian<N> = rng.random();
1068
1069        assert!(a.coordinates.iter().all(|&x| -1.0 < x && x < 1.0));
1070
1071        // This test will fail ~1e-3008 percent of the time - it's probably fine
1072        if N == 10_000 {
1073            assert!(a.coordinates.iter().any(|&x| x < 0.0));
1074        }
1075    }
1076
1077    parameterize_vector_length!(random_in_range, [2, 3, 4, 8, 16, 32, 10_000]);
1078
1079    #[test]
1080    fn to_unit() {
1081        let a = Cartesian::from((2.0, 0.0, 0.0));
1082        let (Unit(unit_a), _) = a
1083            .to_unit()
1084            .expect("hard-coded vector should have non-zero length");
1085        assert_eq!(unit_a, [1.0, 0.0, 0.0].into());
1086
1087        let (Unit(unit_a), _) = a.to_unit_unchecked();
1088        assert_eq!(unit_a, [1.0, 0.0, 0.0].into());
1089
1090        let Unit(unit_a) = Unit::<Cartesian<3>>::try_from([1.0, 0.0, 0.0])
1091            .expect("hard-coded vector should have non-zero length");
1092        assert_eq!(unit_a, [1.0, 0.0, 0.0].into());
1093
1094        let a = Cartesian::from((3.0, 0.0, 4.0));
1095        let (Unit(unit_a), _) = a
1096            .to_unit()
1097            .expect("hard-coded vector should have non-zero length");
1098        assert_eq!(unit_a, [3.0 / 5.0, 0.0, 4.0 / 5.0].into());
1099
1100        let (Unit(unit_a), _) = a.to_unit_unchecked();
1101        assert_eq!(unit_a, [3.0 / 5.0, 0.0, 4.0 / 5.0].into());
1102
1103        let Unit(unit_a) = Unit::<Cartesian<3>>::try_from([3.0, 0.0, 4.0])
1104            .expect("hard-coded vector should have non-zero length");
1105        assert_eq!(unit_a, [3.0 / 5.0, 0.0, 4.0 / 5.0].into());
1106
1107        let a = Cartesian::from((0.0, 0.0, 0.0));
1108        assert!(matches!(a.to_unit(), Err(Error::InvalidVectorMagnitude)));
1109        assert!(matches!(
1110            Unit::<Cartesian<3>>::try_from([0.0, 0.0, 0.0]),
1111            Err(Error::InvalidVectorMagnitude)
1112        ));
1113    }
1114
1115    #[test]
1116    fn sum() {
1117        let total: Cartesian<2> = [Cartesian::from((1.0, 2.0)), Cartesian::from((-2.0, -1.0))]
1118            .into_iter()
1119            .sum();
1120
1121        assert_eq!(total, [-1.0, 1.0].into());
1122    }
1123
1124    #[test]
1125    fn perpendicular() {
1126        let v = Cartesian::from([1.0, -4.5]);
1127        assert_eq!(v.perpendicular(), [4.5, 1.0].into());
1128    }
1129
1130    #[test]
1131    fn test_rotationmatrix_display() {
1132        let m = RotationMatrix::<3>::default();
1133        let s = format!("{m}");
1134        assert_eq!(
1135            s,
1136            "[[1, 0, 0]
1137 [0, 1, 0]
1138 [0, 0, 1]]"
1139        );
1140
1141        let m = RotationMatrix::<2>::default();
1142        let s = format!("{m}");
1143        assert_eq!(
1144            s,
1145            "[[1, 0]
1146 [0, 1]]"
1147        );
1148    }
1149
1150    #[rstest(
1151        angle => [
1152            Angle::default(),
1153            Angle::from(std::f64::consts::FRAC_PI_3),
1154            Angle::from(999.9 * std::f64::consts::PI / 1.234)
1155        ]
1156    )]
1157    fn test_rotationmatrix_invert_2d(angle: Angle) {
1158        let transposed = RotationMatrix::from(angle).inverted();
1159        let inverted = RotationMatrix::from(angle.inverted());
1160        for (a, b) in transposed.rows().iter().zip(inverted.rows().iter()) {
1161            assert_relative_eq!(a, b);
1162        }
1163    }
1164    #[rstest(
1165        versor => [
1166            Versor::default(),
1167            Versor::from_axis_angle(
1168                Unit::try_from([1.0, 0.0, 0.0]).unwrap(), std::f64::consts::PI / 1.234
1169            )
1170        ]
1171    )]
1172    fn test_rotationmatrix_invert_3d(versor: Versor) {
1173        let transposed = RotationMatrix::from(versor).inverted();
1174        let inverted = RotationMatrix::from(versor.inverted());
1175        for (a, b) in transposed.rows().iter().zip(inverted.rows().iter()) {
1176            assert_relative_eq!(a, b);
1177        }
1178    }
1179}
hoomd_vector/cartesian.rs

hoomd_vector/
cartesian.rs
