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, Outer, Rotate, Unit, Vector, Wedge};
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    /// assert_eq!(2, Cartesian::<2>::n_dimensions());
316    /// assert_eq!(3, Cartesian::<3>::n_dimensions());
317    /// # Ok(())
318    /// # }
319    /// ```
320    #[inline]
321    fn n_dimensions() -> usize {
322        N
323    }
324}
325
326impl<const N: usize> Add for Cartesian<N> {
327    type Output = Self;
328
329    #[inline]
330    fn add(self, rhs: Self) -> Self {
331        let mut coordinates = [0.0; N];
332
333        for (result, (a, b)) in coordinates
334            .iter_mut()
335            .zip(self.coordinates.iter().zip(rhs.coordinates.iter()))
336        {
337            *result = a + b;
338        }
339        Self { coordinates }
340    }
341}
342
343impl<const N: usize> AddAssign for Cartesian<N> {
344    #[inline]
345    fn add_assign(&mut self, rhs: Self) {
346        for (result, a) in self.coordinates.iter_mut().zip(rhs.coordinates.iter()) {
347            *result += a;
348        }
349    }
350}
351
352impl<const N: usize> Div<f64> for Cartesian<N> {
353    type Output = Self;
354
355    #[inline]
356    fn div(self, rhs: f64) -> Self {
357        Self {
358            coordinates: self.coordinates.map(|x| x / rhs),
359        }
360    }
361}
362
363impl<const N: usize> DivAssign<f64> for Cartesian<N> {
364    #[inline]
365    fn div_assign(&mut self, rhs: f64) {
366        self.coordinates = self.coordinates.map(|x| x / rhs);
367    }
368}
369
370impl<const N: usize> Mul<f64> for Cartesian<N> {
371    type Output = Self;
372
373    #[inline]
374    fn mul(self, rhs: f64) -> Self {
375        Self {
376            coordinates: self.coordinates.map(|x| x * rhs),
377        }
378    }
379}
380
381impl<const N: usize> MulAssign<f64> for Cartesian<N> {
382    #[inline]
383    fn mul_assign(&mut self, rhs: f64) {
384        self.coordinates = self.coordinates.map(|x| x * rhs);
385    }
386}
387
388impl<const N: usize> Sub for Cartesian<N> {
389    type Output = Self;
390
391    #[inline]
392    fn sub(self, rhs: Self) -> Self {
393        let mut coordinates = [0.0; N];
394        for (result, (a, b)) in coordinates
395            .iter_mut()
396            .zip(self.coordinates.iter().zip(rhs.coordinates.iter()))
397        {
398            *result = a - b;
399        }
400        Self { coordinates }
401    }
402}
403
404impl<const N: usize> SubAssign for Cartesian<N> {
405    #[inline]
406    fn sub_assign(&mut self, rhs: Self) {
407        for (result, a) in self.coordinates.iter_mut().zip(rhs.coordinates) {
408            *result -= a;
409        }
410    }
411}
412
413impl<const N: usize> fmt::Display for Cartesian<N> {
414    #[inline]
415    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
416        write!(
417            f,
418            "[{}]",
419            self.coordinates
420                .iter()
421                .map(f64::to_string)
422                .collect::<Vec<String>>()
423                .join(", ")
424        )
425    }
426}
427
428impl<const N: usize> Neg for Cartesian<N> {
429    type Output = Self;
430    #[inline]
431    fn neg(self) -> Self::Output {
432        Self {
433            coordinates: self.coordinates.map(|x| -x),
434        }
435    }
436}
437
438impl Cross for Cartesian<3> {
439    #[inline]
440    fn cross(&self, other: &Self) -> Self {
441        Cartesian::from((
442            self[1] * other[2] - self[2] * other[1],
443            self[2] * other[0] - self[0] * other[2],
444            self[0] * other[1] - self[1] * other[0],
445        ))
446    }
447}
448
449impl Wedge for Cartesian<3> {
450    type Bivector = Cartesian<3>;
451
452    #[inline]
453    /// Compute the wedge product of two vectors.
454    ///
455    /// ```math
456    /// \textbf{A}=\textbf{a}\wedge{\textbf{b}}
457    /// ```
458    ///
459    /// # Example
460    ///
461    /// ```
462    /// use hoomd_vector::{Cartesian, Wedge};
463    ///
464    /// let a = Cartesian::from([1.0, 0.0, 0.0]);
465    /// let b = Cartesian::from([0.0, 1.0, 0.0]);
466    /// assert_eq!(a.wedge(&b), [0.0, 0.0, 1.0].into());
467    /// ```
468    fn wedge(&self, other: &Self) -> Self::Bivector {
469        self.cross(other)
470    }
471}
472
473impl Cartesian<4> {
474    /// Compute the `N-1`-ary ([co-unary]) cross product of a four-dimensional vector.
475    ///
476    /// This function is a generalization of the cross product to general dimension,
477    /// identifying a unique vector perpendicular to `N-1` other vectors. In two
478    /// dimensions, this is [`Cartesian::<2>::perpendicular`], and in three dimensions this is
479    /// simply [`Cross`].
480    ///
481    /// In geometric algebra terms, this is the Hodge dual of the exterior (Wedge)
482    /// product. Geometrically, the dot product of the resulting vector with all inputs
483    /// is zero and the magnitude of the vector is the hypervolume of the
484    /// parallelepiped spanned by the inputs.
485    ///
486    /// [co-unary]: https://ncatlab.org/nlab/show/cross+product#counary
487    #[inline]
488    #[must_use]
489    pub fn counary_cross(vectors: &[Self; 3]) -> Self {
490        /// Calculate a 2x2 matrix determinant while compensating for errors.
491        /// <https://pharr.org/matt/blog/2019/11/03/difference-of-floats>
492        #[inline(always)]
493        fn diff_of_products(a: f64, b: f64, c: f64, d: f64) -> f64 {
494            let cd = c * d;
495            let err = (-c).mul_add(d, cd);
496            let dop = a.mul_add(b, -cd);
497            dop + err
498        }
499        let u = &vectors[0];
500        let v = &vectors[1];
501        let w = &vectors[2];
502
503        // 2x2 Minors via Kahan's Exact FMA
504        let m01 = diff_of_products(v[0], w[1], v[1], w[0]);
505        let m02 = diff_of_products(v[0], w[2], v[2], w[0]);
506        let m03 = diff_of_products(v[0], w[3], v[3], w[0]);
507        let m12 = diff_of_products(v[1], w[2], v[2], w[1]);
508        let m13 = diff_of_products(v[1], w[3], v[3], w[1]);
509        let m23 = diff_of_products(v[2], w[3], v[3], w[2]);
510
511        // Association is important here, as we can accumulate small sign errors if we
512        // do not correctly group terms!
513        [
514            (u[1] * m23 - u[2] * m13 + u[3] * m12),
515            -(u[0] * m23 - u[2] * m03 + u[3] * m02),
516            (u[0] * m13 - u[1] * m03 + u[3] * m01),
517            -(u[0] * m12 - u[1] * m02 + u[2] * m01),
518        ]
519        .into()
520    }
521}
522
523impl<const N: usize> Distribution<Cartesian<N>> for StandardUniform {
524    /// Sample a Cartesian vector from the uniform [-1, 1] hypercube.
525    ///
526    /// Each coordinate in the vector is in the closed range [-1, 1].
527    ///
528    /// # Example
529    ///
530    /// ```
531    /// use hoomd_vector::Cartesian;
532    /// use rand::{RngExt, SeedableRng, rngs::StdRng};
533    ///
534    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
535    /// let mut rng = StdRng::seed_from_u64(1);
536    /// let v: Cartesian<3> = rng.random();
537    /// # Ok(())
538    /// # }
539    /// ```
540    #[inline]
541    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> Cartesian<N> {
542        let uniform = Uniform::new_inclusive(-1.0, 1.0)
543            .expect("hard-coded range should form a valid distribution");
544        Cartesian {
545            coordinates: array::from_fn(|_| uniform.sample(rng)),
546        }
547    }
548}
549
550impl<const N: usize, T> Index<T> for Cartesian<N>
551where
552    T: Into<usize> + std::slice::SliceIndex<[f64], Output = f64>,
553{
554    type Output = f64;
555    /// Get the value of the vector at coordinate i.
556    ///
557    /// # Example
558    /// ```
559    /// use hoomd_vector::Cartesian;
560    ///
561    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
562    /// let v = Cartesian::<3>::try_from(3..6)?;
563    /// assert_eq!((v[0], v[1], v[2]), (3.0, 4.0, 5.0));
564    /// # Ok(())
565    /// # }
566    /// ```
567    #[inline]
568    fn index(&self, index: T) -> &Self::Output {
569        &self.coordinates[index]
570    }
571}
572
573impl Cartesian<2> {
574    /// Construct a 2-vector perpendicular to self.
575    ///
576    /// Given a vector $`(v_x, v_y)`$ `perpendicular` returns the vector
577    /// rotated by $`\pi/2`$:
578    /// ```math
579    /// (-v_y, v_x)
580    /// ```
581    ///
582    /// # Example
583    /// ```
584    /// use hoomd_vector::Cartesian;
585    ///
586    /// let v = Cartesian::from([1.0, -4.5]);
587    /// assert_eq!(v.perpendicular(), [4.5, 1.0].into());
588    /// ```
589    #[inline]
590    #[must_use]
591    pub fn perpendicular(self) -> Self {
592        Cartesian::from([-self[1], self[0]])
593    }
594}
595
596impl Wedge for Cartesian<2> {
597    type Bivector = f64;
598
599    /// Compute the wedge product of two vectors.
600    ///
601    /// ```math
602    /// \textbf{A}=\textbf{a}\wedge{\textbf{b}}
603    /// ```
604    ///
605    /// # Example
606    ///
607    /// ```
608    /// use hoomd_vector::{Cartesian, Wedge};
609    ///
610    /// let a = Cartesian::from([2.0, 1.0]);
611    /// let b = Cartesian::from([3.0, 1.0]);
612    ///
613    /// assert_eq!(a.wedge(&b), -1.0);
614    /// ```
615    #[inline]
616    fn wedge(&self, other: &Self) -> Self::Bivector {
617        self[0] * other[1] - self[1] * other[0]
618    }
619}
620
621impl<const N: usize> Outer for Cartesian<N> {
622    type Tensor = Matrix<N, N>;
623
624    #[inline]
625    fn outer(&self, other: &Self) -> Self::Tensor {
626        self.to_column_matrix().matmul(&other.to_row_matrix())
627    }
628}
629
630impl<const N: usize, T> IndexMut<T> for Cartesian<N>
631where
632    T: Into<usize> + std::slice::SliceIndex<[f64], Output = f64>,
633{
634    /// Get a mutable reference to the value of the vector at coordinate i.
635    ///
636    /// # Example
637    /// ```
638    /// use hoomd_vector::Cartesian;
639    ///
640    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
641    /// let mut v = Cartesian::<3>::try_from(3..6)?;
642    /// assert_eq!((v[0], v[1], v[2]), (3.0, 4.0, 5.0));
643    /// v[0] += 1.0;
644    /// assert_eq!(v[0], 4.0);
645    /// # Ok(())
646    /// # }
647    /// ```
648    #[inline]
649    fn index_mut(&mut self, index: T) -> &mut Self::Output {
650        &mut self.coordinates[index]
651    }
652}
653
654impl<const N: usize> Sum for Cartesian<N> {
655    #[inline]
656    fn sum<I>(iter: I) -> Self
657    where
658        I: Iterator<Item = Self>,
659    {
660        iter.fold(Cartesian::default(), |acc, x| acc + x)
661    }
662}
663
664/// Rotate vectors efficiently.
665///
666/// Construct a [`RotationMatrix`] to efficiently rotate many vectors by the same rotation.
667///
668/// See:
669/// * [`RotationMatrix::from<Angle>`]
670/// * [`RotationMatrix::from<Versor>`]
671///
672/// [`RotationMatrix`] _intentionally_ does not implement [`Rotation`](crate::Rotation).
673/// [`Angle`](crate::Angle) and [`Versor`](crate::Versor) are representations of
674/// rotations that are often the most effective and numerically stable to
675/// manipulate.
676#[serde_as]
677#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
678pub struct RotationMatrix<const N: usize> {
679    /// Rows of the rotation matrix.
680    #[serde_as(as = "[_; N]")]
681    pub(crate) rows: [Cartesian<N>; N],
682}
683
684impl<const N: usize> From<RotationMatrix<N>> for Matrix<N, N> {
685    #[inline]
686    fn from(value: RotationMatrix<N>) -> Self {
687        Self {
688            rows: value.rows().map(|arr| arr.coordinates),
689        }
690    }
691}
692
693impl<const N: usize> RotationMatrix<N> {
694    /// Get the rows of the rotation matrix.
695    ///
696    /// # Example
697    /// ```
698    /// use hoomd_vector::{Angle, InnerProduct, RotationMatrix, Vector};
699    /// use std::f64::consts::PI;
700    ///
701    /// let a = Angle::from(PI / 2.0);
702    ///
703    /// let matrix = RotationMatrix::from(a);
704    /// assert!(matrix.rows()[0].dot(&[0.0, -1.0].into()) > 0.99);
705    /// assert!(matrix.rows()[1].dot(&[1.0, 0.0].into()) > 0.99);
706    /// ```
707    #[inline]
708    #[must_use]
709    pub fn rows(&self) -> [Cartesian<N>; N] {
710        self.rows
711    }
712
713    /// Create a matrix that performs the inverse rotation.
714    ///
715    /// Matrix inversion is cheaper than [`Angle`](crate::Angle) ->
716    /// [`RotationMatrix`] and [`Versor`](crate::Versor) -> [`RotationMatrix`]
717    /// conversions. When you need both rotations, convert once and then invert.
718    ///
719    /// # Example
720    /// ```
721    /// use hoomd_vector::{Angle, InnerProduct, RotationMatrix, Vector};
722    /// use std::f64::consts::PI;
723    ///
724    /// let a = Angle::from(PI / 2.0);
725    ///
726    /// let matrix = RotationMatrix::from(a);
727    /// let inverted_matrix = matrix.inverted();
728    /// assert!(inverted_matrix.rows()[0].dot(&[0.0, 1.0].into()) > 0.99);
729    /// assert!(inverted_matrix.rows()[1].dot(&[-1.0, 0.0].into()) > 0.99);
730    /// ```
731    #[inline]
732    #[must_use]
733    pub fn inverted(self) -> Self {
734        Self {
735            rows: array::from_fn(|i| array::from_fn::<_, N, _>(|j| self.rows()[j][i]).into()),
736        }
737    }
738}
739
740impl<const N: usize> Default for RotationMatrix<N> {
741    /// Create an identity matrix.
742    ///
743    /// ```math
744    /// \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}
745    /// ```
746    /// ,
747    /// ```math
748    /// \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}
749    /// ```
750    /// , and so on.
751    ///
752    /// # Example
753    ///
754    /// ```
755    /// use hoomd_vector::RotationMatrix;
756    ///
757    /// let identity = RotationMatrix::<3>::default();
758    /// ```
759    #[inline]
760    fn default() -> RotationMatrix<N> {
761        RotationMatrix {
762            rows: array::from_fn(|i| array::from_fn(|j| if i == j { 1.0 } else { 0.0 }).into()),
763        }
764    }
765}
766
767impl<const N: usize> fmt::Display for RotationMatrix<N> {
768    #[inline]
769    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
770        write!(
771            f,
772            "[{}]",
773            self.rows
774                .iter()
775                .map(Cartesian::<N>::to_string)
776                .collect::<Vec<String>>()
777                .join("\n ")
778        )
779    }
780}
781
782impl<const N: usize> Rotate<Cartesian<N>> for RotationMatrix<N> {
783    type Matrix = RotationMatrix<N>;
784
785    #[inline]
786    /// Rotate a [`Cartesian<N>`] by a [`RotationMatrix`]
787    ///
788    /// # Examples
789    /// ```
790    /// use approxim::assert_relative_eq;
791    /// use hoomd_vector::{Angle, Cartesian, Rotate, RotationMatrix};
792    /// use std::f64::consts::PI;
793    ///
794    /// let v = Cartesian::from([-1.0, 0.0]);
795    /// let a = Angle::from(PI / 2.0);
796    ///
797    /// let matrix = RotationMatrix::from(a);
798    /// let rotated = matrix.rotate(&v);
799    /// assert_relative_eq!(rotated, [0.0, -1.0].into());
800    /// ```
801    ///
802    /// ```
803    /// use approxim::assert_relative_eq;
804    /// use hoomd_vector::{Cartesian, Rotate, RotationMatrix, Versor};
805    /// use std::f64::consts::PI;
806    ///
807    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
808    /// let a = Cartesian::from([-1.0, 0.0, 0.0]);
809    /// let v = Versor::from_axis_angle([0.0, 0.0, 1.0].try_into()?, PI / 2.0);
810    ///
811    /// let matrix = RotationMatrix::from(v);
812    /// let b = matrix.rotate(&a);
813    /// assert_relative_eq!(b, [0.0, -1.0, 0.0].into());
814    /// # Ok(())
815    /// # }
816    /// ```
817    fn rotate(&self, vector: &Cartesian<N>) -> Cartesian<N> {
818        let mut coordinates = [0.0; N];
819
820        for (result, row) in coordinates.iter_mut().zip(self.rows.iter()) {
821            *result = row.dot(vector);
822        }
823
824        Cartesian { coordinates }
825    }
826}
827
828impl<const N: usize, const K: usize> MatMul<Matrix<N, K>> for RotationMatrix<N> {
829    type Output = Matrix<N, K>;
830
831    #[inline]
832    fn matmul(&self, rhs: &Matrix<N, K>) -> Self::Output {
833        Matrix::from(*self).matmul(rhs)
834    }
835}
836
837impl<const N: usize> Cartesian<N> {
838    /// Convert a [`Cartesian<N>`] into a row matrix [`Matrix<1, N>`].
839    ///
840    /// # Example
841    /// ```
842    /// use hoomd_linear_algebra::matrix::Matrix;
843    /// use hoomd_vector::Cartesian;
844    ///
845    /// let a = Cartesian::from([1.0, -2.0, 3.0]);
846    ///
847    /// let b = a.to_row_matrix();
848    /// assert_eq!(b.rows, [[1.0, -2.0, 3.0]]);
849    /// ```
850    #[inline]
851    #[must_use]
852    pub fn to_row_matrix(self) -> Matrix<1, N> {
853        Matrix {
854            rows: [self.coordinates],
855        }
856    }
857    /// Convert a [`Cartesian<N>`] into a column matrix [`Matrix<N, 1>`].
858    ///
859    /// # Example
860    /// ```
861    /// use hoomd_linear_algebra::matrix::Matrix;
862    /// use hoomd_vector::Cartesian;
863    ///
864    /// let a = Cartesian::from([1.0, -2.0, 3.0]);
865    ///
866    /// let b = a.to_column_matrix();
867    /// assert_eq!(b.rows, [[1.0], [-2.0], [3.0]]);
868    /// ```
869    #[inline]
870    #[must_use]
871    pub fn to_column_matrix(self) -> Matrix<N, 1> {
872        Matrix {
873            rows: std::array::from_fn(|i| [self[i]]),
874        }
875    }
876}
877
878impl<const N: usize> From<Matrix<1, N>> for Cartesian<N> {
879    #[inline]
880    fn from(value: Matrix<1, N>) -> Self {
881        value.rows[0].into()
882    }
883}
884
885#[cfg(test)]
886mod tests {
887    use crate::{Angle, Rotation, Versor};
888
889    use super::*;
890    use approxim::assert_relative_eq;
891    use paste::paste;
892    use rand::{RngExt, SeedableRng, rngs::StdRng};
893    use rstest::rstest;
894
895    // Parameterize a test function over an array of vector lengths
896    macro_rules! parameterize_vector_length {
897        // macro with name as above that takes an identifier (fn) and an expression
898        // $(...),* matches 0 or more expressions (values) separated by commas
899        ($test_body:ident, [$($dim:expr),*]) => {
900
901            // Now, we repeat the test block 0 or more times, one for each $dim
902            $(
903                // paste package combines values in [< >] to form a new ident
904                paste! {
905                    #[test]
906                    fn [< $test_body "_" $dim>]() {
907                        const DIM: usize = $dim;
908                        $test_body::<DIM>();
909                    }
910                }
911            )*
912        };
913    }
914
915    /// Generate a pair of length N vectors.
916    /// The first vector ranges from [0, N-1] and the second ranges from [N, 2*N-1]
917    fn generate_vector_pair<const N: usize>() -> (Cartesian<N>, Cartesian<N>) {
918        (
919            Cartesian::try_from(0..N).unwrap(),
920            Cartesian::try_from(N..N * 2).unwrap(),
921        )
922    }
923
924    fn index<const N: usize>() {
925        let (_, b) = generate_vector_pair::<N>();
926        assert!(zip(0..N, b.coordinates.iter()).all(|(i, &x)| b[i] == x));
927    }
928    parameterize_vector_length!(index, [2, 3, 4, 8, 16, 32]);
929
930    fn index_mut<const N: usize>() {
931        let (a, mut b) = generate_vector_pair::<N>();
932        zip(0..N, b.coordinates.iter_mut()).for_each(|(i, x)| *x = a[i]);
933        assert_eq!(a, b);
934    }
935    parameterize_vector_length!(index_mut, [2, 3, 4, 8, 16, 32]);
936
937    fn add_explicit<const N: usize>() {
938        let (a, b) = generate_vector_pair::<N>();
939        let c = a.add(b);
940
941        let addition_answer: Vec<f64> = (0..(2 * N))
942            .step_by(2)
943            .map(|x| (x + N) as f64)
944            .collect::<Vec<_>>();
945
946        assert_eq!(c, Cartesian::try_from(addition_answer).unwrap());
947    }
948    parameterize_vector_length!(add_explicit, [2, 3, 4, 8, 16, 32]);
949
950    fn add_operator<const N: usize>() {
951        let (a, b) = generate_vector_pair::<N>();
952        let c = a + b;
953
954        let addition_answer: Vec<f64> = (0..(2 * N))
955            .step_by(2)
956            .map(|x| (x + N) as f64)
957            .collect::<Vec<_>>();
958
959        assert_eq!(c, Cartesian::try_from(addition_answer).unwrap());
960    }
961    parameterize_vector_length!(add_operator, [2, 3, 4, 8, 16, 32]);
962
963    fn add_assign<const N: usize>() {
964        let (a, b) = generate_vector_pair::<N>();
965        let mut c = a;
966        c += b;
967
968        let addition_answer: Vec<f64> = (0..(2 * N))
969            .step_by(2)
970            .map(|x| (x + N) as f64)
971            .collect::<Vec<_>>();
972
973        assert_eq!(c, Cartesian::try_from(addition_answer).unwrap());
974    }
975    parameterize_vector_length!(add_assign, [2, 3, 4, 8, 16, 32]);
976
977    fn sub_operator<const N: usize>() {
978        let (a, b) = generate_vector_pair::<N>();
979        let c = a - b;
980
981        let subtraction_answer = [-(N as f64); N];
982
983        assert_eq!(c, subtraction_answer.into());
984    }
985    parameterize_vector_length!(sub_operator, [2, 3, 4, 8, 16, 32]);
986
987    fn sub_assign<const N: usize>() {
988        let (a, b) = generate_vector_pair::<N>();
989        let mut c = a;
990        c -= b;
991
992        let subtraction_answer = [-(N as f64); N];
993
994        assert_eq!(c, subtraction_answer.into());
995    }
996
997    parameterize_vector_length!(sub_assign, [2, 3, 4, 8, 16, 32]);
998
999    fn mul_operator<const N: usize>() {
1000        let (a, _) = generate_vector_pair::<N>();
1001        let b = 12.0;
1002        let c = a * b;
1003
1004        let multiplication_answer: Vec<f64> = (0..N).map(|x| (x as f64) * b).collect::<Vec<_>>();
1005
1006        assert_eq!(c, Cartesian::try_from(multiplication_answer).unwrap());
1007    }
1008    parameterize_vector_length!(mul_operator, [2, 3, 4, 8, 16, 32]);
1009
1010    fn mul_assign<const N: usize>() {
1011        let (mut a, _) = generate_vector_pair::<N>();
1012        let b = 12.0;
1013        a *= b;
1014
1015        let multiplication_answer: Vec<f64> = (0..N).map(|x| (x as f64) * b).collect::<Vec<_>>();
1016
1017        assert_eq!(a, Cartesian::try_from(multiplication_answer).unwrap());
1018    }
1019
1020    parameterize_vector_length!(mul_assign, [2, 3, 4, 8, 16, 32]);
1021
1022    fn div_operator<const N: usize>() {
1023        let (a, _) = generate_vector_pair::<N>();
1024        let b = 12.0;
1025        let c = a / b;
1026
1027        let division_answer: Vec<f64> = (0..N).map(|x| (x as f64) / b).collect::<Vec<_>>();
1028
1029        assert_eq!(c, Cartesian::try_from(division_answer).unwrap());
1030    }
1031    parameterize_vector_length!(div_operator, [2, 3, 4, 8, 16, 32]);
1032
1033    fn div_assign<const N: usize>() {
1034        let (mut a, _) = generate_vector_pair::<N>();
1035        let b = 12.0;
1036        a /= b;
1037
1038        let division_answer: Vec<f64> = (0..N).map(|x| (x as f64) / b).collect::<Vec<_>>();
1039
1040        assert_eq!(a, Cartesian::try_from(division_answer).unwrap());
1041    }
1042
1043    parameterize_vector_length!(div_assign, [2, 3, 4, 8, 16, 32]);
1044
1045    fn compute_add_ref_ref<const N: usize>(a: &Cartesian<N>, b: &Cartesian<N>) -> Cartesian<N> {
1046        *a + *b
1047    }
1048
1049    fn compute_add_ref_type<const N: usize>(a: &Cartesian<N>, b: Cartesian<N>) -> Cartesian<N> {
1050        *a + b
1051    }
1052
1053    fn compute_add_type_ref<const N: usize>(a: Cartesian<N>, b: &Cartesian<N>) -> Cartesian<N> {
1054        a + *b
1055    }
1056
1057    fn add_with_refs<const N: usize>() {
1058        let (a, b) = generate_vector_pair::<N>();
1059
1060        let addition_answer = Cartesian::try_from(
1061            (0..(2 * N))
1062                .step_by(2)
1063                .map(|x| (x + N) as f64)
1064                .collect::<Vec<_>>(),
1065        )
1066        .unwrap();
1067
1068        let c = compute_add_ref_ref(&a, &b);
1069        assert_eq!(c, addition_answer);
1070
1071        let c = compute_add_ref_type(&a, b);
1072        assert_eq!(c, addition_answer);
1073
1074        let c = compute_add_type_ref(a, &b);
1075        assert_eq!(c, addition_answer);
1076    }
1077    parameterize_vector_length!(add_with_refs, [2, 3, 4, 8, 16, 32]);
1078
1079    fn dot<const N: usize>() {
1080        let (a, b) = generate_vector_pair::<N>();
1081        let c = a.dot(&b);
1082
1083        let n = N as f64;
1084
1085        // Analytical solution to sum of k * (n+k), k = 0 to n-1
1086        let dot_ans = (5.0 * n.powi(3) - 6.0 * n.powi(2) + n) / 6.0;
1087
1088        assert_eq!(c, dot_ans);
1089    }
1090    parameterize_vector_length!(dot, [2, 3, 4, 8, 16, 32]);
1091
1092    fn neg<const N: usize>() {
1093        let a: Cartesian<N> = Cartesian::try_from(0..N).unwrap();
1094        let b = -a;
1095        for (i, j) in zip(a.coordinates.iter(), b.coordinates.iter()) {
1096            assert_eq!(*i, -j);
1097        }
1098    }
1099    parameterize_vector_length!(neg, [2, 3, 4, 8, 16, 32]);
1100
1101    #[test]
1102    fn cross() {
1103        let (a, b) = generate_vector_pair::<3>();
1104        let c = a.cross(&b);
1105
1106        // Analytical solution
1107        let cross_ans = Cartesian::from((-3.0, 6.0, -3.0));
1108
1109        assert_eq!(c, cross_ans);
1110
1111        let a = Cartesian::from([1.0, 0.0, 0.0]);
1112        let b = Cartesian::from([0.0, 1.0, 0.0]);
1113        assert_eq!(a.cross(&b), [0.0, 0.0, 1.0].into());
1114
1115        let a = Cartesian::from([0.0, 3.0, 0.0]);
1116        let b = Cartesian::from([0.0, 0.0, 2.0]);
1117        assert_eq!(a.cross(&b), [6.0, 0.0, 0.0].into());
1118
1119        let a = Cartesian::from([2.0, 0.0, 0.0]);
1120        let b = Cartesian::from([0.0, 0.0, 4.0]);
1121        assert_eq!(a.cross(&b), [0.0, -8.0, 0.0].into());
1122    }
1123
1124    #[test]
1125    fn display() {
1126        let a = Cartesian::from([1.5, 2.125, -3.875]);
1127        let s = format!("{a}");
1128        assert_eq!(s, "[1.5, 2.125, -3.875]");
1129
1130        let a = Cartesian::from([10.0, 20.0, 30.0, 40.0]);
1131        let s = format!("{a}");
1132        assert_eq!(s, "[10, 20, 30, 40]");
1133    }
1134
1135    #[test]
1136    fn from_2_tuple() {
1137        let a = Cartesian::from((3.0, 0.125));
1138        assert_eq!(a.coordinates, [3.0, 0.125]);
1139    }
1140
1141    #[test]
1142    fn from_3_tuple() {
1143        let a = Cartesian::from((-0.5, 2.0, 18.125));
1144        assert_eq!(a.coordinates, [-0.5, 2.0, 18.125]);
1145    }
1146
1147    fn from_vec<const N: usize>() {
1148        let mut vec = Vec::with_capacity(N);
1149
1150        assert_eq!(
1151            Cartesian::<N>::try_from(vec.clone()),
1152            Err(Error::InvalidVectorLength)
1153        );
1154
1155        for i in 0..N {
1156            vec.push(i as f64 * 0.5);
1157        }
1158        let a = Cartesian::<N>::try_from(vec.clone()).unwrap();
1159
1160        assert_eq!(vec, Vec::from(a.coordinates));
1161
1162        vec.push(1.0);
1163        assert_eq!(
1164            Cartesian::<N>::try_from(vec.clone()),
1165            Err(Error::InvalidVectorLength)
1166        );
1167    }
1168    parameterize_vector_length!(from_vec, [2, 3, 4, 8, 16, 32]);
1169
1170    fn random_in_range<const N: usize>() {
1171        // Loosely verify we are drawing from the correct distribution
1172        let mut rng = StdRng::seed_from_u64(1);
1173        let a: Cartesian<N> = rng.random();
1174
1175        assert!(a.coordinates.iter().all(|&x| -1.0 < x && x < 1.0));
1176
1177        // This test will fail ~1e-3008 percent of the time - it's probably fine
1178        if N == 10_000 {
1179            assert!(a.coordinates.iter().any(|&x| x < 0.0));
1180        }
1181    }
1182
1183    parameterize_vector_length!(random_in_range, [2, 3, 4, 8, 16, 32, 10_000]);
1184
1185    #[test]
1186    fn to_unit() {
1187        let a = Cartesian::from((2.0, 0.0, 0.0));
1188        let (Unit(unit_a), _) = a
1189            .to_unit()
1190            .expect("hard-coded vector should have non-zero length");
1191        assert_eq!(unit_a, [1.0, 0.0, 0.0].into());
1192
1193        let (Unit(unit_a), _) = a.to_unit_unchecked();
1194        assert_eq!(unit_a, [1.0, 0.0, 0.0].into());
1195
1196        let Unit(unit_a) = Unit::<Cartesian<3>>::try_from([1.0, 0.0, 0.0])
1197            .expect("hard-coded vector should have non-zero length");
1198        assert_eq!(unit_a, [1.0, 0.0, 0.0].into());
1199
1200        let a = Cartesian::from((3.0, 0.0, 4.0));
1201        let (Unit(unit_a), _) = a
1202            .to_unit()
1203            .expect("hard-coded vector should have non-zero length");
1204        assert_eq!(unit_a, [3.0 / 5.0, 0.0, 4.0 / 5.0].into());
1205
1206        let (Unit(unit_a), _) = a.to_unit_unchecked();
1207        assert_eq!(unit_a, [3.0 / 5.0, 0.0, 4.0 / 5.0].into());
1208
1209        let Unit(unit_a) = Unit::<Cartesian<3>>::try_from([3.0, 0.0, 4.0])
1210            .expect("hard-coded vector should have non-zero length");
1211        assert_eq!(unit_a, [3.0 / 5.0, 0.0, 4.0 / 5.0].into());
1212
1213        let a = Cartesian::from((0.0, 0.0, 0.0));
1214        assert!(matches!(a.to_unit(), Err(Error::InvalidVectorMagnitude)));
1215        assert!(matches!(
1216            Unit::<Cartesian<3>>::try_from([0.0, 0.0, 0.0]),
1217            Err(Error::InvalidVectorMagnitude)
1218        ));
1219    }
1220
1221    #[test]
1222    fn sum() {
1223        let total: Cartesian<2> = [Cartesian::from((1.0, 2.0)), Cartesian::from((-2.0, -1.0))]
1224            .into_iter()
1225            .sum();
1226
1227        assert_eq!(total, [-1.0, 1.0].into());
1228    }
1229
1230    #[test]
1231    fn perpendicular() {
1232        let v = Cartesian::from([1.0, -4.5]);
1233        assert_eq!(v.perpendicular(), [4.5, 1.0].into());
1234    }
1235
1236    #[test]
1237    fn test_rotationmatrix_display() {
1238        let m = RotationMatrix::<3>::default();
1239        let s = format!("{m}");
1240        assert_eq!(
1241            s,
1242            "[[1, 0, 0]
1243 [0, 1, 0]
1244 [0, 0, 1]]"
1245        );
1246
1247        let m = RotationMatrix::<2>::default();
1248        let s = format!("{m}");
1249        assert_eq!(
1250            s,
1251            "[[1, 0]
1252 [0, 1]]"
1253        );
1254    }
1255
1256    #[rstest(
1257        angle => [
1258            Angle::default(),
1259            Angle::from(std::f64::consts::FRAC_PI_3),
1260            Angle::from(999.9 * std::f64::consts::PI / 1.234)
1261        ]
1262    )]
1263    fn test_rotationmatrix_invert_2d(angle: Angle) {
1264        let transposed = RotationMatrix::from(angle).inverted();
1265        let inverted = RotationMatrix::from(angle.inverted());
1266        for (a, b) in transposed.rows().iter().zip(inverted.rows().iter()) {
1267            assert_relative_eq!(a, b);
1268        }
1269    }
1270    #[rstest(
1271        versor => [
1272            Versor::default(),
1273            Versor::from_axis_angle(
1274                Unit::try_from([1.0, 0.0, 0.0]).unwrap(), std::f64::consts::PI / 1.234
1275            )
1276        ]
1277    )]
1278    fn test_rotationmatrix_invert_3d(versor: Versor) {
1279        let transposed = RotationMatrix::from(versor).inverted();
1280        let inverted = RotationMatrix::from(versor.inverted());
1281        for (a, b) in transposed.rows().iter().zip(inverted.rows().iter()) {
1282            assert_relative_eq!(a, b);
1283        }
1284    }
1285
1286    #[rstest]
1287    fn test_generalized_cross_product_combinations(
1288        #[values(1.0, -2.5, 4.0)] r1: f64,
1289        #[values(1.0, 3.0, -5.1)] r2: f64,
1290        #[values(2.0, -1.5, 0.0)] r3: f64,
1291        #[values(0.0, 33.0, -12.5, 0.04)] x: f64,
1292        #[values(0.0, 1.0, 3.0, -2.6)] y: f64,
1293        #[values(0.0, 2.0, -1.5, -0.1)] z: f64,
1294    ) {
1295        // Construct a lower-triangular matrix of vectors
1296        let vecs = [
1297            Cartesian::from([r1, 0.0, 0.0, 0.0]),
1298            Cartesian::from([x, r2, 0.0, 0.0]),
1299            Cartesian::from([y, z, r3, 0.0]),
1300        ];
1301
1302        let cross = Cartesian::<4>::counary_cross(&vecs);
1303
1304        // Validate the result is perpendicular to all inputs
1305        for v in &vecs {
1306            assert_relative_eq!(cross.dot(v), 0.0);
1307        }
1308
1309        // Validate the magnitude is exactly the spanned volume
1310        // For a triangular matrix, the volume is the product of the diagonals.
1311        let expected_volume_sq = (r1 * r2 * r3).powi(2);
1312        assert_relative_eq!(cross.norm_squared(), expected_volume_sq);
1313    }
1314}
hoomd_vector/cartesian.rs

hoomd_vector/
cartesian.rs
