hoomd_manifold/
hyperbolic_angle.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//! Implements a three-dimensional representation of SO(3,1) on Minkowski
5//! space.
6
7use num::complex::Complex;
8use rand::{
9    Rng,
10    distr::{Distribution, StandardUniform, Uniform},
11};
12use serde::{Deserialize, Serialize};
13use std::{f64::consts::PI, fmt};
14
15use crate::HyperbolicRotationMatrix;
16
17/// Combine rotations and boosts in hyperbolic space.
18///
19/// "Hyperbolic rotations" describe elements of the group SO(N,1), which
20/// preserve Hyperbolics embedded in [`Minkowski<N+1>`]. Transformations include
21/// pure spatial rotations as well as "boosts".
22///
23/// [`Minkowski<N+1>`]: crate::Minkowski
24///
25/// For two-dimensional hyperbolic surfaces, use [`HyperbolicAngle`] to
26/// implement elements of SO(2,1) which rotate points about the z-axis or boost
27///  points along the x- and y-axes. SO(2,1) transformations represent
28/// symmetries of two-dimensional hyperbolic space. The algebra so(2,1) is
29/// generated by the matrices
30/// ```math
31/// K_x = \begin{bmatrix} 0&0&1\\ 0&0&0\\1&0&0 \end{bmatrix} \qquad\qquad
32/// K_y = \begin{bmatrix} 0&0&0\\ 0&0&1\\0&1&0 \end{bmatrix}\qquad\qquad
33/// J_z = \begin{bmatrix} 0&-1&0\\ 1&0&0\\0&0&0 \end{bmatrix}
34/// ```
35/// which satisfy the commutation relations
36/// ```math
37/// [J_z,K_x] = K_y \qquad [K_x,K_y] = -J_z \qquad [J_z, K_y] = -K_x
38/// ```
39/// The group elements SO(2,1) can be obtained by exponentiating the generators
40/// near identity,
41/// i.e.,
42/// ```math
43/// G = \exp(aJ_z + b K_x + c K_y)
44/// ```
45/// In particular,
46/// ```math
47/// \exp(aJ_z) = \begin{bmatrix}\cos(a)&-\sin(a)&0\\ \sin(a)&\cos(a)&0 \\ 0&0&1\end{bmatrix}
48/// \\ \exp(bK_x) = \begin{bmatrix}\cosh(b)&0&\sinh(b)\\ 0&1&0 \\ \sinh(b)&0&\cosh(b)\end{bmatrix}
49/// \\ \exp(cK_y) = \begin{bmatrix}1&0&0\\ 0&\cosh(c)&\sinh(c) \\ 0&\sinh(c)&\cosh(c)\end{bmatrix}
50/// ```
51/// [`HyperbolicAngle`] is a struct of three real numbers $`(a,b,c)`$ which
52/// parameterize the generators in this exponential map. The first value $`a`$
53/// may be interpreted as a rotation angle, while the latter two numbers $`b, c`$
54/// may be interpreted as rapidities for boosts along the x- and y- directions,
55/// respectively. Use [`HyperbolicRotationMatrix`] to generate the matrix from
56/// the values defined by [`HyperbolicAngle`].
57///
58/// Rotation about z axis:
59/// ```
60/// use approxim::assert_relative_eq;
61/// use hoomd_manifold::{
62///     HyperbolicAngle, HyperbolicRotate, HyperbolicRotationMatrix, Minkowski,
63/// };
64/// use std::f64::consts::PI;
65///
66/// let v = Minkowski::from([1.0, 0.0, 1.0]);
67/// let rotation_about_z = HyperbolicAngle::from((PI / 2.0, 0.0_f64, 0.0_f64));
68/// let matrix = HyperbolicRotationMatrix::from(rotation_about_z);
69/// let rotated = matrix.hyperbolic_rotate(&v);
70/// assert_relative_eq!(rotated, [0.0, 1.0, 1.0].into(), epsilon = 1e-12);
71/// ```
72///
73/// Boost in the y direction:
74/// ```
75/// use approxim::assert_relative_eq;
76/// use hoomd_manifold::{
77///     HyperbolicAngle, HyperbolicRotate, HyperbolicRotationMatrix, Minkowski,
78/// };
79///
80/// let v = Minkowski::from([1.0, 0.0, 1.0]);
81/// let boost_in_y = HyperbolicAngle::from((0.0_f64, 0.0_f64, 0.5_f64));
82/// let matrix = HyperbolicRotationMatrix::from(boost_in_y);
83/// let boosted = matrix.hyperbolic_rotate(&v);
84/// assert_relative_eq!(
85///     boosted,
86///     [1.0, (0.5_f64).sinh(), (0.5_f64).cosh()].into(),
87///     epsilon = 1e-12
88/// );
89/// ```
90///
91/// The default value `HyperbolicAngle::Default` returns the tuple
92/// `(0.0,0.0,0.0)`, which corresponds to the identify transformation.
93#[derive(Clone, Copy, Debug, Default, PartialEq, Serialize, Deserialize)]
94pub struct HyperbolicAngle {
95    /// Get the rotation angle in radians.
96    pub angles: (f64, f64, f64),
97}
98
99impl HyperbolicAngle {
100    /// Modulo the second and third components of a `HyperbolicAngle` by $`2 \pi`$.
101    #[inline]
102    #[must_use]
103    pub fn to_reduced(self) -> Self {
104        Self {
105            angles: (
106                self.angles.0.rem_euclid(2.0 * PI),
107                self.angles.1,
108                self.angles.2,
109            ),
110        }
111    }
112}
113
114impl From<HyperbolicAngle> for HyperbolicRotationMatrix<3> {
115    /// Create a matrix belonging to SO(2,1) from the generators of the algebra so(2,1).
116    #[inline]
117    fn from(angle_list: HyperbolicAngle) -> HyperbolicRotationMatrix<3> {
118        let (a, b, c) = angle_list.angles;
119        if (a, b, c) == (0.0_f64, 0.0_f64, 0.0_f64) {
120            HyperbolicRotationMatrix {
121                rows: [
122                    [1.0, 0.0, 0.0].into(),
123                    [0.0, 1.0, 0.0].into(),
124                    [0.0, 0.0, 1.0].into(),
125                ],
126            }
127        } else {
128            let arg_sq = -a.powi(2) + b.powi(2) + c.powi(2);
129            let arg = Complex::new(arg_sq, 0.0);
130            let arg_sqrt = arg.sqrt();
131            let ch = arg_sqrt.cosh();
132            let sh = arg_sqrt.sinh();
133            let sh_c = Complex::new(
134                arg_sqrt.re * sh.re - arg_sqrt.im * sh.im,
135                arg_sqrt.re * sh.im + arg_sqrt.im * sh.re,
136            );
137
138            HyperbolicRotationMatrix {
139                rows: [
140                    [
141                        ((Complex::new(c.powi(2), 0.0) + ch.scale(b.powi(2) - a.powi(2)))
142                            .scale(1.0 / arg_sq))
143                        .re,
144                        ((Complex::new(b * c, 0.0) - ch.scale(b * c) + sh_c.scale(a))
145                            .scale(-1.0 / arg_sq))
146                        .re,
147                        ((Complex::new(-a * c, 0.0) + ch.scale(a * c) - sh_c.scale(b))
148                            .scale(-1.0 / arg_sq))
149                        .re,
150                    ]
151                    .into(),
152                    [
153                        ((Complex::new(b * c, 0.0) - ch.scale(b * c) - sh_c.scale(a))
154                            .scale(-1.0 / arg_sq))
155                        .re,
156                        ((Complex::new(b.powi(2), 0.0) + ch.scale(c.powi(2) - a.powi(2)))
157                            .scale(1.0 / arg_sq))
158                        .re,
159                        ((Complex::new(a * b, 0.0) - ch.scale(a * b) - sh_c.scale(c))
160                            .scale(-1.0 / arg_sq))
161                        .re,
162                    ]
163                    .into(),
164                    [
165                        ((Complex::new(a * c, 0.0) - ch.scale(a * c) - sh_c.scale(b))
166                            .scale(-1.0 / arg_sq))
167                        .re,
168                        ((Complex::new(-a * b, 0.0) + ch.scale(a * b) - sh_c.scale(c))
169                            .scale(-1.0 / arg_sq))
170                        .re,
171                        ((Complex::new(a.powi(2), 0.0) - ch.scale(b.powi(2) + c.powi(2)))
172                            .scale(-1.0 / arg_sq))
173                        .re,
174                    ]
175                    .into(),
176                ],
177            }
178        }
179    }
180}
181
182impl From<(f64, f64, f64)> for HyperbolicAngle {
183    /// A tuple of three real numbers (a,b,c) describing the coefficients for the
184    /// generators (`J_z`, `K_x`, `K_y`) of the algebra so(2,1).
185    #[inline]
186    fn from(value: (f64, f64, f64)) -> Self {
187        Self { angles: value }
188    }
189}
190
191impl fmt::Display for HyperbolicAngle {
192    /// Format a [`HyperbolicAngle`] as `[{v[0]}, {v[1]}, {v[2]}]`.
193    #[inline]
194    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
195        write!(
196            f,
197            "[{}, {}, {}]",
198            self.angles.0, self.angles.1, self.angles.2
199        )
200    }
201}
202
203impl Distribution<HyperbolicAngle> for StandardUniform {
204    /// Sample a random hyperbolic angle.
205    ///
206    /// Note that the distribution is not
207    /// truly uniform---the boost magnitude is sampled from a square-root
208    /// distribution (between zero and one) while the angle is sampled
209    /// from a uniform distribution between 0 and $`2\pi`$.
210    ///
211    /// Every point on the Hyperbolic can be generated by a boost in the
212    /// x-direction followed by a rotation. The [Baker-Campbell-Hausdorff Lemma]
213    /// gives an expression for the relationship between subsequent SO(2,1)
214    /// transformations and the so(2,1) algebra generators. Explicitly,
215    ///
216    /// [Baker-Campbell-Hausdorff Lemma]: https://en.wikipedia.org/wiki/Baker%E2%80%93Campbell%E2%80%93Hausdorff_formula
217    ///
218    /// ```math
219    /// \exp[\theta J_z]\exp[v K_x] = \exp[(\theta + \frac{1}{12}\theta v^2)J_z + (v-\frac{1}{12}\theta^2v)K_x + \frac{1}{2}\theta v K_y + \text{higher order}]
220    /// ```
221    /// To avoid clustering around the cusp of the Hyperbolic, the random
222    /// boost parameter $`v`$ is sampled from a uniform distribution between
223    /// 0 and 1 and then passed through a square root function.
224    ///
225    /// # Example
226    ///
227    /// ```
228    /// use approxim::assert_relative_eq;
229    /// use hoomd_manifold::{
230    ///     HyperbolicAngle, HyperbolicRotate, HyperbolicRotationMatrix, Minkowski,
231    /// };
232    /// use hoomd_vector::{Metric, Vector};
233    /// use rand::{RngExt, SeedableRng, rngs::StdRng};
234    ///
235    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
236    /// let mut rng = StdRng::seed_from_u64(1);
237    /// let v: HyperbolicAngle = rng.random();
238    ///
239    /// let matrix = HyperbolicRotationMatrix::from(v);
240    /// let origin = Minkowski::from([0.0, 0.0, 1.0]);
241    /// let point = matrix.hyperbolic_rotate(&origin);
242    /// assert_relative_eq!(
243    ///     -1.0,
244    ///     point.distance_squared(&Minkowski::<3>::default()),
245    ///     epsilon = 1e-12
246    /// );
247    /// # Ok(())
248    /// # }
249    /// ```
250    #[inline]
251    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> HyperbolicAngle {
252        let uniform_angle =
253            Uniform::new(0.0, 2.0 * PI).expect("hard-coded distribution should be valid");
254        let uniform_boost =
255            Uniform::new(0.0, 1.0).expect("hard-coded distribution should be valid");
256        let theta = uniform_angle.sample(rng);
257        let v: f64 = uniform_boost.sample(rng);
258        let v_sqrt = v.sqrt();
259        HyperbolicAngle::from((
260            theta + v * theta / 12.0,
261            v_sqrt - theta.powi(2) * v_sqrt / 12.0,
262            theta * v_sqrt / 2.0,
263        ))
264    }
265}
266
267#[cfg(test)]
268mod tests {
269    use super::*;
270    use crate::{HyperbolicRotate, HyperbolicRotationMatrix, Minkowski};
271    use approxim::assert_relative_eq;
272    use rand::{RngExt, SeedableRng, rngs::StdRng};
273    use rstest::rstest;
274    use std::f64::consts::PI;
275
276    // Test named cases of three input values (hyperbolic angle tuple, Minkowski input, and answer)
277    #[rstest]
278    #[case::z_pi((PI, 0.0_f64, 0.0_f64), (1.0,0.0,0.0), (-1.0,0.0,0.0))]
279    #[case::z_pi_halves((PI/2.0, 0.0_f64, 0.0_f64), (1.0,0.0,0.0), (0.0,1.0,0.0))]
280    #[case::x_one((0.0_f64, 1.0_f64, 0.0_f64), (1.0,0.0,0.0), ((1.0_f64).cosh(),0.0,(1.0_f64).sinh()))]
281    #[case::x_tenth((0.0_f64, 0.1_f64, 0.0_f64), (1.0,0.0,1.0), ((0.1_f64).cosh()+(0.1_f64).sinh(),0.0,(0.1_f64).cosh()+(0.1_f64).sinh()))]
282    #[case::x_negative_one((0.0_f64, -1.0_f64, 0.0_f64), (1.0,0.0,0.0), ((1.0_f64).cosh(),0.0,-(1.0_f64).sinh()))]
283    #[case::y_one((0.0_f64, 0.0_f64, 1.0_f64), (0.0,1.0,0.0), (0.0,(1.0_f64).cosh(),(1.0_f64).sinh()))]
284    #[case::y_tenth((0.0_f64, 0.0_f64, 0.1_f64), (1.0,0.0,1.0), (1.0,(0.1_f64).sinh(),(0.1_f64).cosh()))]
285    #[case::y_negative_one((0.0_f64, 0.0_f64, -1.0_f64), (0.0,1.0,0.0), (0.0,(1.0_f64).cosh(),-(1.0_f64).sinh()))]
286
287    fn rotate_hyperbolic(
288        #[case] hyperbolic_angle: (f64, f64, f64),
289        #[case] vec: (f64, f64, f64),
290        #[case] ans: (f64, f64, f64),
291    ) {
292        let hyperbolic_angle = HyperbolicAngle::from(hyperbolic_angle);
293        let vec = Minkowski::from(vec);
294        let ans = Minkowski::from(ans);
295
296        let matrix = HyperbolicRotationMatrix::from(hyperbolic_angle);
297        let transformed = matrix.hyperbolic_rotate(&vec);
298        assert_relative_eq!(transformed, ans, epsilon = 1e-12);
299    }
300
301    #[test]
302    fn display() {
303        let a = HyperbolicAngle::from((1.3, 1.1, 1.2));
304        let s = format!("[{}, {}, {}]", a.angles.0, a.angles.1, a.angles.2);
305        assert_eq!(s, "[1.3, 1.1, 1.2]");
306    }
307
308    #[test]
309    fn reduced() {
310        let a = HyperbolicAngle::from((2.0 * PI + 0.1, 1.2, 2.3)).to_reduced();
311        assert_relative_eq!(a.angles.0, 0.1, epsilon = 1e-12);
312        assert_relative_eq!(a.angles.1, 1.2, epsilon = 1e-12);
313        assert_relative_eq!(a.angles.2, 2.3, epsilon = 1e-12);
314
315        let b = HyperbolicAngle::from((2.0 * PI, 1.2, 2.3)).to_reduced();
316        assert_relative_eq!(b.angles.0, 0.0, epsilon = 1e-12);
317        assert_relative_eq!(b.angles.1, 1.2, epsilon = 1e-12);
318        assert_relative_eq!(b.angles.2, 2.3, epsilon = 1e-12);
319
320        let c = HyperbolicAngle::from((0.5, 1.2, 2.3)).to_reduced();
321        assert_relative_eq!(c.angles.0, 0.5, epsilon = 1e-12);
322        assert_relative_eq!(c.angles.1, 1.2, epsilon = 1e-12);
323        assert_relative_eq!(c.angles.2, 2.3, epsilon = 1e-12);
324    }
325
326    #[test]
327    fn random() {
328        let mut rng = StdRng::seed_from_u64(13);
329
330        for _ in 0..10000 {
331            let a_sample: HyperbolicAngle = rng.random();
332            let a = a_sample.to_reduced();
333            assert!(a.angles.0 >= 0.0 && a.angles.0 < 2.0 * PI);
334            assert!(a.angles.1 < 1.0);
335            assert!(a.angles.2 < PI);
336        }
337    }
338}
hoomd_manifold/hyperbolic_angle.rs

hoomd_manifold/
hyperbolic_angle.rs
