hoomd_manifold/minkowski.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 vector types in Minkowski space.
5
6use approxim::{approx_derive::RelativeEq, assert_relative_eq};
7use rand::{
8 Rng,
9 distr::{Distribution, StandardUniform, Uniform},
10};
11use serde::{Deserialize, Serialize};
12use serde_with::serde_as;
13use std::{
14 array,
15 f64::consts::PI,
16 fmt,
17 iter::zip,
18 ops::{Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign},
19};
20
21use crate::{Error, HyperbolicRotate};
22use hoomd_utility::valid::PositiveReal;
23use hoomd_vector::{Metric, Vector};
24
25/// A vector in N-dimensional Minkowski space.
26///
27/// [`Minkowski<N>`] implements (N-1,1)-dimensional Minkowski space with the
28/// metric signature $`(+ , \cdots , + , -)`$. [`Minkowski`] supports
29/// [`Vector`] operations such as vector addition and rescaling.
30///
31/// ## Constructing Minkowski vectors
32///
33/// Similar to [`Cartesian`], N-dimensional vectors can be
34/// constructed using an array of (real-valued) coordinates. Three- and
35/// four-dimensional vectors can also be constructed from tuples:
36/// ```
37/// use hoomd_manifold::Minkowski;
38///
39/// fn from_array() -> Minkowski<5> {
40/// Minkowski::from([1.0, 2.0, 3.0, 4.0, 5.0])
41/// }
42/// fn from_tuples() -> Minkowski<3> {
43/// Minkowski::from((6.0, 7.0, 8.0))
44/// }
45/// ```
46///
47/// ## Operating on Minkowski vectors
48///
49/// [`Minkowski`] implements everything from [`Vector`], which includes vector
50/// addition/subtraction and multiplication by a scalar.
51///
52/// ```
53/// use hoomd_manifold::Minkowski;
54///
55/// // Vector addition
56/// let mut a = Minkowski::from([1.0, 1.0, 1.0, 1.0]);
57/// let mut b = Minkowski::from([0.0, 0.0, 0.0, 2.0]);
58/// a += b;
59///
60/// // Multiplication by a scalar
61/// let mut c = a * 4.0;
62///
63/// // Division by a scalar
64/// c /= 2.0;
65///
66/// assert_eq!(c, [2.0, 2.0, 2.0, 6.0].into());
67/// ```
68///
69/// The distance metric on Minkowski space is given by the "spacetime interval"
70/// ```math
71/// d_M^2(\vec{u},\vec{v}) = (\vec{u}-\vec{v})^T \eta (\vec{u}-\vec{v})
72/// = (u_1-v_1)^2 +\cdots + (u_{N-1}-v_{N-1})^2 - (u_N - v_N)^2
73/// ```
74/// Note that because this metric is not positive-definite, [`Minkowski`] this
75/// "spacetime interval" is not a true inner-product, and therefore
76/// [`Minkowski`] does not implement the methods of [`InnerProduct`].
77///
78/// [`InnerProduct`]: hoomd_vector::InnerProduct
79/// [`Cartesian`]: hoomd_vector::Cartesian
80///
81/// ```
82/// use hoomd_manifold::Minkowski;
83/// use hoomd_vector::Metric;
84///
85/// let x = Minkowski::from([1.0, 0.0, 5.0]);
86/// let y = Minkowski::from([0.0, 0.0, 3.0]);
87/// assert_eq!(-3.0, x.distance_squared(&y));
88/// ```
89
90#[serde_as]
91#[derive(Clone, Copy, Debug, PartialEq, RelativeEq, Serialize, Deserialize)]
92#[approx(epsilon_type = f64)]
93pub struct Minkowski<const N: usize> {
94 /// The vector's coordinates.
95 ///
96 /// The final component is the one associated with a minus sign (-) in the metric.
97 #[serde_as(as = "[_; N]")]
98 #[approx(into_iter)]
99 pub coordinates: [f64; N],
100}
101
102impl<const N: usize> Default for Minkowski<N> {
103 /// Create a 0 vector in Minkowski space.
104 ///
105 /// # Example
106 /// ```
107 /// use hoomd_manifold::Minkowski;
108 ///
109 /// let v = Minkowski::<3>::default();
110 /// assert_eq!(v, [0.0; 3].into())
111 /// ```
112 #[inline]
113 fn default() -> Self {
114 Minkowski::from([0.0; N])
115 }
116}
117
118impl<const N: usize> From<[f64; N]> for Minkowski<N> {
119 /// Create a vector in Minkowski space with the given coordinates.
120 ///
121 /// The last component has a (-) signature, while the preceding coordinates
122 /// have (+) signatures in the metric.
123 ///
124 /// # Example
125 /// ```
126 /// use hoomd_manifold::Minkowski;
127 ///
128 /// let v = Minkowski::from([1.0, 2.0]);
129 /// ```
130 #[inline]
131 fn from(coordinates: [f64; N]) -> Self {
132 Self { coordinates }
133 }
134}
135
136impl From<(f64, f64)> for Minkowski<2> {
137 #[inline]
138 fn from(coordinates: (f64, f64)) -> Self {
139 Self {
140 coordinates: coordinates.into(),
141 }
142 }
143}
144
145impl From<(f64, f64, f64)> for Minkowski<3> {
146 #[inline]
147 fn from(coordinates: (f64, f64, f64)) -> Self {
148 Self {
149 coordinates: coordinates.into(),
150 }
151 }
152}
153
154impl From<(f64, f64, f64, f64)> for Minkowski<4> {
155 #[inline]
156 fn from(coordinates: (f64, f64, f64, f64)) -> Self {
157 Self {
158 coordinates: coordinates.into(),
159 }
160 }
161}
162
163impl<const N: usize> TryFrom<Vec<f64>> for Minkowski<N> {
164 type Error = Error;
165
166 /// Create a vector in Minkowski with coordinates given by a [`Vec<f64>`]
167 ///
168 /// # Example
169 /// ```
170 /// use hoomd_manifold::Minkowski;
171 ///
172 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
173 /// let v = Minkowski::<3>::try_from(vec![5.0, 4.0, 3.0])?;
174 /// assert_eq!(v, [5.0, 4.0, 3.0].into());
175 /// # Ok(())
176 /// # }
177 /// ```
178 /// <div class="warning">
179 ///
180 /// Use `Minkowski::From<[f64; N]>` in performance critical code.
181 ///
182 /// </div>
183 #[inline]
184 fn try_from(value: Vec<f64>) -> Result<Self, Self::Error> {
185 let coordinates = value.try_into().map_err(|_| Error::InvalidVectorLength)?;
186 Ok(Self { coordinates })
187 }
188}
189
190impl<const N: usize> TryFrom<std::ops::Range<usize>> for Minkowski<N> {
191 type Error = Error;
192
193 /// Create a vector in Minkowski space with coordinates given by a range.
194 ///
195 /// # Example
196 /// ```
197 /// use hoomd_manifold::Minkowski;
198 ///
199 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
200 /// let v = Minkowski::<3>::try_from(1..4)?;
201 /// assert_eq!(v, [1.0, 2.0, 3.0].into());
202 /// # Ok(())
203 /// # }
204 /// ```
205 #[inline]
206 fn try_from(value: std::ops::Range<usize>) -> Result<Self, Self::Error> {
207 if value.len() != N {
208 return Err(Error::InvalidVectorLength);
209 }
210
211 // The default value of 0 will never be used due to the above error check.
212 let mut iter = value;
213 let coordinates = array::from_fn(|_| iter.next().unwrap_or(0) as f64);
214 Ok(Self { coordinates })
215 }
216}
217
218impl<const N: usize> Metric for Minkowski<N> {
219 /// Computes the squared distance between two points in Minkowski space.
220 ///
221 /// Employs the "mostly plusses" metric signature (+ ... + -).
222 /// ```math
223 /// d^2_M(\vec{x},\vec{y}) = -(x_N-y_N)^2 + \sum_{i=1}^{N-1} (x_i - y_i)^2
224 /// ```
225 ///
226 /// # Example
227 /// ```
228 /// use hoomd_manifold::Minkowski;
229 /// use hoomd_vector::{Metric, Vector};
230 ///
231 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
232 /// let x = Minkowski::from([0.0, 2.0, 3.0]);
233 /// let y = Minkowski::from([1.0, 0.0, 0.0]);
234 /// assert_eq!(-4.0, x.distance_squared(&y));
235 /// # Ok(())
236 /// # }
237 /// ```
238 #[inline]
239 fn distance_squared(&self, other: &Self) -> f64 {
240 let last_component = -(self.coordinates[N - 1] - other.coordinates[N - 1]).powi(2);
241 zip(
242 self.coordinates[0..N - 1].iter(),
243 other.coordinates[0..N - 1].iter(),
244 )
245 .fold(last_component, |product, x| product + (x.0 - x.1).powi(2))
246 }
247 #[inline]
248 fn distance(&self, other: &Self) -> f64 {
249 ((self.distance_squared(other)).abs()).sqrt()
250 }
251 #[inline]
252 fn n_dimensions(&self) -> usize {
253 N
254 }
255}
256
257impl<const N: usize> Vector for Minkowski<N> {}
258
259impl<const N: usize> Add for Minkowski<N> {
260 type Output = Self;
261
262 #[inline]
263 fn add(self, rhs: Self) -> Self {
264 let mut coordinates = [0.0; N];
265
266 for (result, (a, b)) in coordinates
267 .iter_mut()
268 .zip(self.coordinates.iter().zip(rhs.coordinates.iter()))
269 {
270 *result = a + b;
271 }
272 Self { coordinates }
273 }
274}
275
276impl<const N: usize> AddAssign for Minkowski<N> {
277 #[inline]
278 fn add_assign(&mut self, rhs: Self) {
279 for (result, a) in self.coordinates.iter_mut().zip(rhs.coordinates.iter()) {
280 *result += a;
281 }
282 }
283}
284
285impl<const N: usize> Div<f64> for Minkowski<N> {
286 type Output = Self;
287
288 #[inline]
289 fn div(self, rhs: f64) -> Self {
290 let mut coordinates = [0.0; N];
291
292 for (result, a) in coordinates.iter_mut().zip(self.coordinates) {
293 *result = a / rhs;
294 }
295 Self { coordinates }
296 }
297}
298
299impl<const N: usize> DivAssign<f64> for Minkowski<N> {
300 #[inline]
301 fn div_assign(&mut self, rhs: f64) {
302 for result in &mut self.coordinates {
303 *result /= rhs;
304 }
305 }
306}
307
308impl<const N: usize> Mul<f64> for Minkowski<N> {
309 type Output = Self;
310
311 #[inline]
312 fn mul(self, rhs: f64) -> Self {
313 let mut coordinates = [0.0; N];
314 for (result, a) in coordinates.iter_mut().zip(self.coordinates) {
315 *result = a * rhs;
316 }
317 Self { coordinates }
318 }
319}
320
321impl<const N: usize> MulAssign<f64> for Minkowski<N> {
322 #[inline]
323 fn mul_assign(&mut self, rhs: f64) {
324 for result in &mut self.coordinates {
325 *result *= rhs;
326 }
327 }
328}
329
330impl<const N: usize> Sub for Minkowski<N> {
331 type Output = Self;
332
333 #[inline]
334 fn sub(self, rhs: Self) -> Self {
335 let mut coordinates = [0.0; N];
336 for (result, (a, b)) in coordinates
337 .iter_mut()
338 .zip(self.coordinates.iter().zip(rhs.coordinates.iter()))
339 {
340 *result = a - b;
341 }
342 Self { coordinates }
343 }
344}
345
346impl<const N: usize> SubAssign for Minkowski<N> {
347 #[inline]
348 fn sub_assign(&mut self, rhs: Self) {
349 for (result, a) in self.coordinates.iter_mut().zip(rhs.coordinates) {
350 *result -= a;
351 }
352 }
353}
354
355impl<const N: usize> fmt::Display for Minkowski<N> {
356 #[inline]
357 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
358 write!(
359 f,
360 "[{}]",
361 self.coordinates
362 .iter()
363 .map(f64::to_string)
364 .collect::<Vec<String>>()
365 .join(", ")
366 )
367 }
368}
369
370impl<const N: usize> Neg for Minkowski<N> {
371 type Output = Self;
372 #[inline]
373 fn neg(self) -> Self::Output {
374 let mut result = self;
375 result.coordinates.iter_mut().for_each(|x| *x = -*x);
376 result
377 }
378}
379
380impl<const N: usize, T> Index<T> for Minkowski<N>
381where
382 T: Into<usize> + std::slice::SliceIndex<[f64], Output = f64>,
383{
384 type Output = f64;
385 /// Get the value of the vector at coordinate i.
386 ///
387 /// # Example
388 /// ```
389 /// use hoomd_manifold::Minkowski;
390 ///
391 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
392 /// let v = Minkowski::<3>::try_from(4..7)?;
393 /// assert_eq!((v[0], v[1], v[2]), (4.0, 5.0, 6.0));
394 /// # Ok(())
395 /// # }
396 /// ```
397 #[inline]
398 fn index(&self, index: T) -> &Self::Output {
399 &self.coordinates[index]
400 }
401}
402
403impl<const N: usize, T> IndexMut<T> for Minkowski<N>
404where
405 T: Into<usize> + std::slice::SliceIndex<[f64], Output = f64>,
406{
407 /// Get a mutable reference to the value of the vector at coordinate `i`.
408 ///
409 /// # Example
410 /// ```
411 /// use hoomd_manifold::Minkowski;
412 ///
413 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
414 /// let mut v = Minkowski::<3>::try_from(4..7)?;
415 /// assert_eq!((v[0], v[1], v[2]), (4.0, 5.0, 6.0));
416 /// v[0] += 1.0;
417 /// assert_eq!(v[0], 5.0);
418 /// # Ok(())
419 /// # }
420 /// ```
421 #[inline]
422 fn index_mut(&mut self, index: T) -> &mut Self::Output {
423 &mut self.coordinates[index]
424 }
425}
426
427impl<const N: usize> Distribution<Minkowski<N>> for StandardUniform {
428 /// Sample a Minkowski vector from the uniform $` [-1, 1] `$ hypercube.
429 ///
430 /// # Example
431 /// ```
432 /// use hoomd_manifold::Minkowski;
433 /// use rand::{RngExt, SeedableRng, rngs::StdRng};
434 ///
435 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
436 /// let mut rng = StdRng::seed_from_u64(1);
437 /// let v: Minkowski<3> = rng.random();
438 /// # Ok(())
439 /// # }
440 /// ```
441 #[inline]
442 fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> Minkowski<N> {
443 let uniform = Uniform::new_inclusive(-1.0, 1.0)
444 .expect("hard-coded range should form a valid distribution");
445 Minkowski {
446 coordinates: array::from_fn(|_| uniform.sample(rng)),
447 }
448 }
449}
450
451/// Point on the top sheet of a Hyperboloid.
452///
453/// [`Hyperbolic`] implements an embedding of the top sheet of an
454/// (N-1)-dimensional two-sheeted hyperboloid in N-dimensional Minkowski space.
455/// This surface has constant negative curvature and therefore serves as a model
456/// of (N-1)-dimensional hyperbolic space.
457///
458/// Explicitly, for N-dimensional Minkowski space with metric $`\eta =
459/// \operatorname{diag}(+,\cdots,+,-)`$, the hyperboloid with skirt width $`R`$ is
460/// defined by the set of points with components satisfying
461/// ```math
462/// x_1^2 +\cdots x_{N-1}^2 - x_{N}^2 = -R^2
463/// ```
464/// Where the "top sheet" is defined by the $`x_N>0`$ solutions. In Minkowski
465/// space, the hyperboloid has a natural interpretation as the set of points
466/// with the same spacetime interval
467/// ```math
468/// \Delta s^2 = \vec{x}^T \eta \vec{x} = x_1^2 +\cdots x_{N-1}^2 - x_{N}^2
469/// ```
470///
471/// Note that the skirt width is fixed at $`R=1.0`$. In simulation, the global
472/// curvature may be tuned by changing the length scale of interactions.
473///
474/// [`Hyperbolic`] implements a [`Metric`] that computes the distance of the
475/// geodesic passing between two points on a hyperboloid with some given skirt
476/// width.
477///
478/// Two points on the hyperboloid with skirt width $` R = 1.0 `$:
479/// ```
480/// use hoomd_manifold::{Hyperbolic, Minkowski};
481/// use hoomd_vector::Metric;
482///
483/// let x = Hyperbolic::from_minkowski_coordinates([0.0, 0.0, 1.0].into());
484///
485/// let y = Hyperbolic::from_minkowski_coordinates(
486/// [0.0, 1.0, (2.0_f64).sqrt()].into(),
487/// );
488///
489/// assert_eq!(((2.0_f64).sqrt()).acosh(), x.distance(&y));
490/// ```
491///
492/// [`Metric`]: hoomd_vector::Metric
493#[derive(Clone, Copy, Debug, PartialEq, RelativeEq, Serialize, Deserialize)]
494pub struct Hyperbolic<const N: usize> {
495 /// A point on the surface of the upper sheet of a two-sheeted hyperboloid.
496 point: Minkowski<N>,
497}
498
499impl<const N: usize> Hyperbolic<N> {
500 /// Get the coordinates of the point on the hyperboloid.
501 #[must_use]
502 #[inline]
503 pub fn coordinates(&self) -> &[f64; N] {
504 &self.point.coordinates
505 }
506 /// Get the Minkowski point of the hyperboloid.
507 #[must_use]
508 #[inline]
509 pub fn point(&self) -> &Minkowski<N> {
510 &self.point
511 }
512 /// Create a Hyperbolic point from a Minkowski vector.
513 ///
514 /// # Example
515 /// ```
516 /// use hoomd_manifold::{Hyperbolic, Minkowski};
517 /// use hoomd_vector::Metric;
518 ///
519 /// let x = Hyperbolic::from_minkowski_coordinates([0.0, 0.0, 1.0].into());
520 /// ```
521 #[must_use]
522 #[inline]
523 pub fn from_minkowski_coordinates(point: Minkowski<N>) -> Hyperbolic<N> {
524 let skirt_squared = -point.distance_squared(&Minkowski::<N>::default());
525 assert_relative_eq!(skirt_squared, 1.0_f64, epsilon = 1e-8);
526 Hyperbolic { point }
527 }
528}
529
530impl Hyperbolic<3> {
531 /// Create a point on the surface of a three-dimensional hyperboloid from the polar representation.
532 #[must_use]
533 #[inline]
534 pub fn from_polar_coordinates(v: f64, theta: f64) -> Hyperbolic<3> {
535 let theta_mod = theta.rem_euclid(2.0 * PI);
536 let v_sinh = v.sinh();
537 let point = Minkowski::from([
538 v_sinh * (theta_mod.cos()),
539 v_sinh * (theta_mod.sin()),
540 v.cosh(),
541 ]);
542 Hyperbolic::from_minkowski_coordinates(point)
543 }
544}
545
546impl Hyperbolic<4> {
547 /// Create a point on the surface of a four-dimensional hyperboloid from the spherical representation.
548 #[must_use]
549 #[inline]
550 pub fn from_polar_coordinates(v: f64, theta: f64, phi: f64) -> Hyperbolic<4> {
551 let theta_mod = theta.rem_euclid(2.0 * PI);
552 let phi_mod = phi.rem_euclid(PI);
553 let v_sinh = v.sinh();
554 let point = Minkowski::from([
555 v_sinh * (theta_mod.cos()),
556 v_sinh * (theta_mod.sin()) * (phi_mod.cos()),
557 v_sinh * (theta_mod.sin()) * (phi_mod.sin()),
558 v.cosh(),
559 ]);
560 Hyperbolic::from_minkowski_coordinates(point)
561 }
562}
563
564impl<const N: usize> Hyperbolic<N> {
565 /// Compute the distance from a point to the cusp.
566 ///
567 /// Computes the length of the geodesic passing between the cusp
568 /// $`(0,\cdots,0,\rho)`$ and a given point on the hyperboloid.
569 ///
570 /// # Example
571 /// ```
572 /// use approxim::assert_relative_eq;
573 /// use hoomd_manifold::{Hyperbolic, Minkowski};
574 /// use hoomd_vector::Vector;
575 ///
576 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
577 /// let v: f64 = 4.2;
578 /// let x = Hyperbolic::from_minkowski_coordinates(
579 /// [v.sinh(), 0.0, v.cosh()].into(),
580 /// );
581 /// assert_relative_eq!(v, x.distance_from_cusp(), epsilon = 1e-12);
582 /// # Ok(())
583 /// # }
584 /// ```
585 #[inline]
586 #[must_use]
587 pub fn distance_from_cusp(&self) -> f64 {
588 (self.point.coordinates[N - 1]).acosh()
589 }
590
591 /// Projects points on the hyperboloid onto the Poincare disk/ball.
592 ///
593 /// # Example
594 /// ```
595 /// use approxim::assert_relative_eq;
596 /// use hoomd_manifold::{Hyperbolic, Minkowski};
597 /// use hoomd_vector::Vector;
598 ///
599 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
600 /// let v: f64 = 1.098612;
601 /// let x = Hyperbolic::from_minkowski_coordinates(
602 /// [v.sinh(), 0.0, v.cosh()].into(),
603 /// );
604 /// let projection = x.to_poincare();
605 /// assert_relative_eq!(
606 /// v.sinh() / (v.cosh() + 1.0),
607 /// projection[0],
608 /// epsilon = 1e-12
609 /// );
610 /// assert_relative_eq!(0.0, projection[1], epsilon = 1e-12);
611 /// # Ok(())
612 /// # }
613 /// ```
614 #[inline]
615 #[must_use]
616 pub fn to_poincare(&self) -> Vec<f64> {
617 (0..N - 1)
618 .collect::<Vec<usize>>()
619 .iter()
620 .map(|i| self.point.coordinates[*i] / (1.0 + self.point.coordinates[N - 1]))
621 .collect::<Vec<f64>>()
622 }
623}
624
625impl<const N: usize> Default for Hyperbolic<N> {
626 /// Construct a default point on a hyperboloid.
627 ///
628 /// The default `Hyperbolic<N>` point is on the cusp of a hyperboloid with
629 /// skirt width of 1.0 (i.e., the point $`(0.0, \cdots, 0.0, 1.0)`$).
630 #[inline]
631 fn default() -> Self {
632 let mut zero = Minkowski::<N>::default();
633 zero.coordinates[N - 1] = 1.0;
634 Hyperbolic { point: zero }
635 }
636}
637
638impl Metric for Hyperbolic<3> {
639 /// The distance between two [`Hyperbolic<3>`] points.
640 ///
641 /// Explicitly, the metric for two points $`\vec{u}`$ and $`\vec{v}`$ on a
642 /// hyperboloid
643 ///
644 /// ```math
645 /// d_{H_2}(\vec{u}, \vec{v}) = \operatorname{arccosh}\left[u_3v_3 - u_1v_1 - u_2v_2\right]
646 /// ```
647 /// This choice of metric furnishes a representation of 2-dimensional hyperbolic
648 /// space with Gaussian curvature $`K = -1`$.
649 #[inline(always)]
650 fn distance(&self, other: &Self) -> f64 {
651 let self_coords = self.point.coordinates;
652 let other_coords = other.point.coordinates;
653 (self_coords[2] * other_coords[2]
654 - self_coords[0] * other_coords[0]
655 - self_coords[1] * other_coords[1])
656 .acosh()
657 }
658
659 #[inline]
660 fn distance_squared(&self, other: &Self) -> f64 {
661 self.distance(other).powi(2)
662 }
663
664 #[inline]
665 fn n_dimensions(&self) -> usize {
666 2_usize
667 }
668}
669
670impl Metric for Hyperbolic<4> {
671 /// The distance between two [`Hyperbolic<4>`] points.
672 ///
673 /// Explicitly, the metric for two points $`\vec{u}`$ and $`\vec{v}`$ on a
674 /// hyperboloid is given by
675 ///
676 /// ```math
677 /// d_{H_3}(\vec{u}, \vec{v}) = \operatorname{arccosh}\left[u_4v_4 - u_1v_1 - u_2v_2 - u_3v_3\right]
678 /// ```
679 /// This choice of metric furnishes a representation of 3-dimensional hyperbolic
680 /// space with with Gaussian curvature $`K = -1`$.
681 #[inline]
682 fn distance(&self, other: &Self) -> f64 {
683 let self_coords = self.point.coordinates;
684 let other_coords = other.point.coordinates;
685 (self_coords[3] * other_coords[3]
686 - self_coords[0] * other_coords[0]
687 - self_coords[1] * other_coords[1]
688 - self_coords[2] * other_coords[2])
689 .acosh()
690 }
691
692 #[inline]
693 fn distance_squared(&self, other: &Self) -> f64 {
694 self.distance(other).powi(2)
695 }
696
697 #[inline]
698 fn n_dimensions(&self) -> usize {
699 3_usize
700 }
701}
702
703/// Hyperbolic rotations in Minkowski Space
704///
705/// Construct a [`HyperbolicRotationMatrix`] to apply SO(N-1, 1)
706/// transformations to N-dimensional Minkowski vectors. For Minkowski 4-vectors,
707/// [`Biquaternion`] should be used instead for numerical stability. See
708/// documentation in [`HyperbolicAngle`] for details on SO(2,1) transformations
709/// (i.e., two-dimensional hyperbolic space), and in [`Biquaternion`] for
710/// details on SO(3,1) transformations (i.e., three-dimensional hyperbolic
711/// space).
712///
713/// [`Biquaternion`]: crate::Biquaternion
714/// [`HyperbolicAngle`]: crate::HyperbolicAngle
715///
716/// In two dimensional hyperbolic space:
717/// ```
718/// use hoomd_manifold::{
719/// HyperbolicAngle, HyperbolicRotate, HyperbolicRotationMatrix, Minkowski,
720/// };
721/// use std::f64::consts::PI;
722///
723/// fn rotate_about_z(minkowski_vector: &Minkowski<3>) -> Minkowski<3> {
724/// let generators = HyperbolicAngle::from((PI, 0.0_f64, 0.0_f64));
725/// let rotation_matrix = HyperbolicRotationMatrix::from(generators);
726/// rotation_matrix.hyperbolic_rotate(&minkowski_vector)
727/// }
728///
729/// fn boost_in_x(minkowski_vector: &Minkowski<3>) -> Minkowski<3> {
730/// let generators = HyperbolicAngle::from((0.0_f64, 0.2_f64, 0.0_f64));
731/// let boost_matrix = HyperbolicRotationMatrix::from(generators);
732/// boost_matrix.hyperbolic_rotate(&minkowski_vector)
733/// }
734/// ```
735#[serde_as]
736#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
737pub struct HyperbolicRotationMatrix<const N: usize> {
738 /// Rows of the rotation matrix.
739 #[serde_as(as = "[_; N]")]
740 pub(crate) rows: [Minkowski<N>; N],
741}
742
743impl<const N: usize> HyperbolicRotate<Minkowski<N>> for HyperbolicRotationMatrix<N> {
744 type Matrix = HyperbolicRotationMatrix<N>;
745
746 #[inline]
747 /// Rotate a [`Minkowski<N>`] by a [`HyperbolicRotationMatrix`]
748 ///
749 /// # Examples
750 ///
751 /// Rotate point in 2D hyperbolic space about z-axis:
752 /// ```
753 /// use hoomd_manifold::{
754 /// HyperbolicAngle, HyperbolicRotate, HyperbolicRotationMatrix, Minkowski,
755 /// };
756 /// use std::f64::consts::PI;
757 ///
758 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
759 /// let v = Minkowski::from([1.0, 0.0, 1.0]);
760 /// let spatial_rotation = HyperbolicAngle::from((PI / 2.0, 0.0_f64, 0.0_f64));
761 /// let matrix = HyperbolicRotationMatrix::from(spatial_rotation);
762 /// let rotated = matrix.hyperbolic_rotate(&v);
763 /// let c = Minkowski::from([(PI / 2.0).cos(), (PI / 2.0).sin(), 1.0]);
764 /// assert_eq!(c, rotated);
765 /// # Ok(())
766 /// # }
767 /// ```
768 ///
769 /// Boost point in 2D hyperbolic space in x direction:
770 /// ```
771 /// use hoomd_manifold::{
772 /// HyperbolicAngle, HyperbolicRotate, HyperbolicRotationMatrix, Minkowski,
773 /// };
774 /// use num::complex::Complex;
775 /// use std::f64::consts::PI;
776 ///
777 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
778 /// let v = Minkowski::from([1.0, 0.0, 1.0]);
779 /// let small_boost = HyperbolicAngle::from((0.0_f64, 0.1_f64, 0.0_f64));
780 /// let matrix = HyperbolicRotationMatrix::from(small_boost);
781 /// let rotated = matrix.hyperbolic_rotate(&v);
782 /// let c = Minkowski::from([
783 /// (0.1_f64).sinh() + (0.1_f64).cosh(),
784 /// 0.0,
785 /// (0.1_f64).sinh() + (0.1_f64).cosh(),
786 /// ]);
787 /// assert_eq!(c, rotated);
788 /// # Ok(())
789 /// # }
790 /// ```
791 ///
792 /// Zero angles and rapidities does nothing:
793 /// ```
794 /// use hoomd_manifold::{
795 /// HyperbolicAngle, HyperbolicRotate, HyperbolicRotationMatrix, Minkowski,
796 /// };
797 /// use num::complex::Complex;
798 /// use std::f64::consts::PI;
799 ///
800 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
801 /// let v = Minkowski::from([1.0, 2.0, 1.0]);
802 /// let identity = HyperbolicAngle::from((0.0_f64, 0.0_f64, 0.0_f64));
803 /// let matrix = HyperbolicRotationMatrix::from(identity);
804 /// let rotated = matrix.hyperbolic_rotate(&v);
805 /// let c = Minkowski::from([1.0, 2.0, 1.0]);
806 /// assert_eq!(c, rotated);
807 /// # Ok(())
808 /// # }
809 /// ```
810 ///
811 /// Rotate point in 3D hyperbolic space about y axis using matrix representation:
812 /// ```
813 /// use approxim::assert_relative_eq;
814 /// use hoomd_manifold::{
815 /// Biquaternion, HyperbolicRotate, HyperbolicRotationMatrix, Minkowski,
816 /// UnitBiquaternion,
817 /// };
818 /// use num::complex::Complex;
819 /// use std::f64::consts::PI;
820 ///
821 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
822 /// let q = Biquaternion::from([
823 /// Complex::new(0.0, 0.0),
824 /// Complex::new((PI / 4.0).sin(), 0.0),
825 /// Complex::new(0.0, 0.0),
826 /// Complex::new((PI / 4.0).cos(), 0.0),
827 /// ]);
828 /// let v = q.to_unit()?;
829 /// let x = Minkowski::from([1.0, 0.0, 0.0, 1.0]);
830 /// let rotation = HyperbolicRotationMatrix::from(v);
831 /// let rotated = rotation.hyperbolic_rotate(&x);
832 /// assert_relative_eq!(rotated, [0.0, 0.0, -1.0, 1.0].into(), epsilon = 1e-12);
833 /// # Ok(())
834 /// # }
835 /// ```
836 ///
837 /// Boost point in 3D hyperbolic space in x direction using biquaternion algebra:
838 /// ```
839 /// use approxim::assert_relative_eq;
840 /// use hoomd_manifold::{
841 /// Biquaternion, HyperbolicRotate, Minkowski, UnitBiquaternion,
842 /// };
843 /// use num::complex::Complex;
844 /// use std::f64::consts::PI;
845 ///
846 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
847 /// let x = Minkowski::from([0.0, 0.0, 0.0, 1.0]);
848 /// let q = Biquaternion::from([
849 /// Complex::new(0.0, 0.25).sin(),
850 /// Complex::new(0.0, 0.0),
851 /// Complex::new(0.0, 0.0),
852 /// Complex::new(0.0, 0.25).cos(),
853 /// ]);
854 /// let v = q.to_unit()?;
855 /// let boosted = v.hyperbolic_rotate(&x);
856 /// assert_relative_eq!(
857 /// boosted,
858 /// [(0.5_f64).sinh(), 0.0, 0.0, (0.5_f64).cosh()].into(),
859 /// epsilon = 1e-12
860 /// );
861 /// # Ok(())
862 /// # }
863 /// ```
864 fn hyperbolic_rotate(&self, vector: &Minkowski<N>) -> Minkowski<N> {
865 let mut coordinates = [0.0; N];
866
867 for (result, row) in coordinates.iter_mut().zip(self.rows.iter()) {
868 *result = zip(row.coordinates.iter(), vector.coordinates.iter())
869 .fold(0.0, |product, x| product + x.0 * x.1);
870 }
871
872 Minkowski { coordinates }
873 }
874}
875
876/// Randomly distribute points locally on a hyperboloid.
877///
878/// [`HyperbolicDisk`] is a uniform distribution of points within distance `r`
879/// of a point on the 2-dimensional hyperboloid.
880///
881/// # Example
882///
883/// ```
884/// use hoomd_manifold::{
885/// Hyperbolic, HyperbolicAngle, HyperbolicDisk, HyperbolicRotate,
886/// HyperbolicRotationMatrix, Minkowski,
887/// };
888/// use hoomd_vector::Metric;
889/// use rand::{RngExt, SeedableRng, distr::Distribution, rngs::StdRng};
890///
891/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
892/// let mut rng = StdRng::seed_from_u64(12);
893///
894/// let v: HyperbolicAngle = rng.random();
895/// let matrix = HyperbolicRotationMatrix::from(v);
896/// let origin = Minkowski::from([0.0, 0.0, 1.0]);
897/// let random_point = Hyperbolic::from_minkowski_coordinates(
898/// matrix.hyperbolic_rotate(&origin),
899/// );
900///
901/// let r = 0.1;
902/// let mut rng_2 = StdRng::seed_from_u64(239);
903/// let disk = HyperbolicDisk {
904/// disk_radius: r.try_into()?,
905/// point: random_point,
906/// };
907/// let transformed_random_point: Hyperbolic<3> = disk.sample(&mut rng_2);
908///
909/// assert!(r > random_point.distance(&transformed_random_point));
910///
911/// # Ok(())
912/// # }
913/// ```
914#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
915pub struct HyperbolicDisk {
916 /// Max distance away from point.
917 pub disk_radius: PositiveReal,
918 /// The center of the disk.
919 pub point: Hyperbolic<3>,
920}
921
922impl Distribution<Hyperbolic<3>> for HyperbolicDisk {
923 /// Sample a random point in the hyperbolic disk.
924 ///
925 /// The implementation translates Minkowski 3-vector `point` along
926 /// the Hyperbolic by maximum distance of `disk_radius`. Note that because SO(2,1) is
927 /// non-Abelian, the point must be transformed to the cusp before a trial
928 /// move is applied (and then the point is transformed back). This ensures
929 /// that the max distance translated by the trial move does not exceed `disk_radius`.
930 #[inline]
931 fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> Hyperbolic<3> {
932 let max_boost = self.disk_radius.get();
933 let point = self.point;
934 let eta = (point.point.coordinates[2]).acosh();
935 let phi = point.point.coordinates[1].atan2(point.point.coordinates[0]);
936 let trial_boost = Uniform::new(0.0, 1.0).expect("r is positive and real");
937 let trial_rotation =
938 Uniform::new(-PI, PI).expect("hard-coded distribution should be valid");
939 let theta = trial_rotation.sample(rng);
940 let v1: f64 = trial_boost.sample(rng);
941 let v = v1.sqrt() * max_boost;
942 let (v_sinh, eta_sinh, eta_cosh, phi_sin, phi_cos) =
943 (v.sinh(), eta.sinh(), eta.cosh(), phi.sin(), phi.cos());
944 let trial_coords = [v_sinh * theta.cos(), v_sinh * theta.sin(), v.cosh()];
945 let transformed_point = Minkowski::from([
946 trial_coords[0] * (eta_cosh * (phi_cos.powi(2)) + phi_sin.powi(2))
947 + trial_coords[1] * phi_sin * phi_cos * (eta_cosh - 1.0)
948 + trial_coords[2] * eta_sinh * phi_cos,
949 trial_coords[0] * phi_sin * phi_cos * (eta_cosh - 1.0)
950 + trial_coords[1] * (eta_cosh * (phi_sin.powi(2)) + phi_cos.powi(2))
951 + trial_coords[2] * eta_sinh * phi_sin,
952 trial_coords[0] * eta_sinh * phi_cos
953 + trial_coords[1] * eta_sinh * phi_sin
954 + trial_coords[2] * eta_cosh,
955 ]);
956 Hyperbolic::from_minkowski_coordinates(transformed_point)
957 }
958}
959
960#[cfg(test)]
961mod tests {
962 use super::*;
963 use approxim::assert_relative_eq;
964 use paste::paste;
965 use rand::{RngExt, SeedableRng, rngs::StdRng};
966
967 // Parameterize a test function over an array of vector lengths
968 macro_rules! parameterize_vector_length {
969 // macro with name as above that takes an identifier (fn) and an expression
970 // $(...),* matches 0 or more expressions (values) separated by commas
971 ($test_body:ident, [$($dim:expr),*]) => {
972
973 // Now, we repeat the test block 0 or more times, one for each $dim
974 $(
975 // paste package combines values in [< >] to form a new ident
976 paste! {
977 #[test]
978 fn [< $test_body "_" $dim>]() {
979 const DIM: usize = $dim;
980 $test_body::<DIM>();
981 }
982 }
983 )*
984 };
985 }
986
987 /// Generate a pair of length N vectors.
988 /// The first vector ranges from [0, N-1] and the second ranges from [N, 2*N-1]
989 fn generate_vector_pair<const N: usize>() -> (Minkowski<N>, Minkowski<N>) {
990 (
991 Minkowski::try_from(0..N).unwrap(),
992 Minkowski::try_from(N..N * 2).unwrap(),
993 )
994 }
995
996 fn index<const N: usize>() {
997 let (_, b) = generate_vector_pair::<N>();
998 assert!(zip(0..N, b.coordinates.iter()).all(|(i, &x)| b[i] == x));
999 }
1000 parameterize_vector_length!(index, [2, 3, 4, 8, 16, 32]);
1001
1002 fn index_mut<const N: usize>() {
1003 let (a, mut b) = generate_vector_pair::<N>();
1004 zip(0..N, b.coordinates.iter_mut()).for_each(|(i, x)| *x = a[i]);
1005 assert_eq!(a, b);
1006 }
1007 parameterize_vector_length!(index_mut, [2, 3, 4, 8, 16, 32]);
1008
1009 fn add_explicit<const N: usize>() {
1010 let (a, b) = generate_vector_pair::<N>();
1011 let c = a.add(b);
1012
1013 let addition_answer: Vec<f64> = (0..(2 * N))
1014 .step_by(2)
1015 .map(|x| (x + N) as f64)
1016 .collect::<Vec<_>>();
1017
1018 assert_eq!(c, Minkowski::try_from(addition_answer).unwrap());
1019 }
1020 parameterize_vector_length!(add_explicit, [2, 3, 4, 8, 16, 32]);
1021
1022 fn add_operator<const N: usize>() {
1023 let (a, b) = generate_vector_pair::<N>();
1024 let c = a + b;
1025
1026 let addition_answer: Vec<f64> = (0..(2 * N))
1027 .step_by(2)
1028 .map(|x| (x + N) as f64)
1029 .collect::<Vec<_>>();
1030
1031 assert_eq!(c, Minkowski::try_from(addition_answer).unwrap());
1032 }
1033 parameterize_vector_length!(add_operator, [2, 3, 4, 8, 16, 32]);
1034
1035 fn add_assign<const N: usize>() {
1036 let (a, b) = generate_vector_pair::<N>();
1037 let mut c = a;
1038 c += b;
1039
1040 let addition_answer: Vec<f64> = (0..(2 * N))
1041 .step_by(2)
1042 .map(|x| (x + N) as f64)
1043 .collect::<Vec<_>>();
1044
1045 assert_eq!(c, Minkowski::try_from(addition_answer).unwrap());
1046 }
1047 parameterize_vector_length!(add_assign, [2, 3, 4, 8, 16, 32]);
1048
1049 fn sub_operator<const N: usize>() {
1050 let (a, b) = generate_vector_pair::<N>();
1051 let c = a - b;
1052
1053 let subtraction_answer = [-(N as f64); N];
1054
1055 assert_eq!(c, subtraction_answer.into());
1056 }
1057 parameterize_vector_length!(sub_operator, [2, 3, 4, 8, 16, 32]);
1058
1059 fn sub_assign<const N: usize>() {
1060 let (a, b) = generate_vector_pair::<N>();
1061 let mut c = a;
1062 c -= b;
1063
1064 let subtraction_answer = [-(N as f64); N];
1065
1066 assert_eq!(c, subtraction_answer.into());
1067 }
1068
1069 parameterize_vector_length!(sub_assign, [2, 3, 4, 8, 16, 32]);
1070
1071 fn mul_operator<const N: usize>() {
1072 let (a, _) = generate_vector_pair::<N>();
1073 let b = 12.0;
1074 let c = a * b;
1075
1076 let multiplication_answer: Vec<f64> = (0..N).map(|x| (x as f64) * b).collect::<Vec<_>>();
1077
1078 assert_eq!(c, Minkowski::try_from(multiplication_answer).unwrap());
1079 }
1080 parameterize_vector_length!(mul_operator, [2, 3, 4, 8, 16, 32]);
1081
1082 fn div_operator<const N: usize>() {
1083 let (a, _) = generate_vector_pair::<N>();
1084 let b = 12.0;
1085 let c = a / b;
1086
1087 let division_answer: Vec<f64> = (0..N).map(|x| (x as f64) / b).collect::<Vec<_>>();
1088
1089 assert_eq!(c, Minkowski::try_from(division_answer).unwrap());
1090 }
1091 parameterize_vector_length!(div_operator, [2, 3, 4, 8, 16, 32]);
1092
1093 fn div_assign<const N: usize>() {
1094 let (mut a, _) = generate_vector_pair::<N>();
1095 let b = 12.0;
1096 a /= b;
1097
1098 let division_answer: Vec<f64> = (0..N).map(|x| (x as f64) / b).collect::<Vec<_>>();
1099
1100 assert_eq!(a, Minkowski::try_from(division_answer).unwrap());
1101 }
1102
1103 parameterize_vector_length!(div_assign, [2, 3, 4, 8, 16, 32]);
1104
1105 fn compute_add_ref_ref<const N: usize>(a: &Minkowski<N>, b: &Minkowski<N>) -> Minkowski<N> {
1106 *a + *b
1107 }
1108
1109 fn compute_add_ref_type<const N: usize>(a: &Minkowski<N>, b: Minkowski<N>) -> Minkowski<N> {
1110 *a + b
1111 }
1112
1113 fn compute_add_type_ref<const N: usize>(a: Minkowski<N>, b: &Minkowski<N>) -> Minkowski<N> {
1114 a + *b
1115 }
1116
1117 fn add_with_refs<const N: usize>() {
1118 let (a, b) = generate_vector_pair::<N>();
1119
1120 let addition_answer = Minkowski::try_from(
1121 (0..(2 * N))
1122 .step_by(2)
1123 .map(|x| (x + N) as f64)
1124 .collect::<Vec<_>>(),
1125 )
1126 .unwrap();
1127
1128 let c = compute_add_ref_ref(&a, &b);
1129 assert_eq!(c, addition_answer);
1130
1131 let c = compute_add_ref_type(&a, b);
1132 assert_eq!(c, addition_answer);
1133
1134 let c = compute_add_type_ref(a, &b);
1135 assert_eq!(c, addition_answer);
1136 }
1137 parameterize_vector_length!(add_with_refs, [2, 3, 4, 8, 16, 32]);
1138
1139 #[test]
1140 fn display() {
1141 let a = Minkowski::from([1.67, -2.125, 42.01]);
1142 let s = format!("{a}");
1143 assert_eq!(s, "[1.67, -2.125, 42.01]");
1144
1145 let a = Minkowski::from([10.0, 20.0, 30.0, 40.0]);
1146 let s = format!("{a}");
1147 assert_eq!(s, "[10, 20, 30, 40]");
1148 }
1149
1150 #[test]
1151 fn from_2_tuple() {
1152 let a = Minkowski::from((13.0, 0.125));
1153 assert_eq!(a.coordinates, [13.0, 0.125]);
1154 }
1155
1156 #[test]
1157 fn from_3_tuple() {
1158 let a = Minkowski::from((-0.5, 2.0, 18.125));
1159 assert_eq!(a.coordinates, [-0.5, 2.0, 18.125]);
1160 }
1161
1162 fn from_vec<const N: usize>() {
1163 let mut vec = Vec::with_capacity(N);
1164
1165 assert_eq!(
1166 Minkowski::<N>::try_from(vec.clone()),
1167 Err(Error::InvalidVectorLength)
1168 );
1169
1170 for i in 0..N {
1171 vec.push(i as f64 * 0.5);
1172 }
1173 let a = Minkowski::<N>::try_from(vec.clone()).unwrap();
1174
1175 assert_eq!(vec, Vec::from(a.coordinates));
1176
1177 vec.push(1.0);
1178 assert_eq!(
1179 Minkowski::<N>::try_from(vec.clone()),
1180 Err(Error::InvalidVectorLength)
1181 );
1182 }
1183 parameterize_vector_length!(from_vec, [2, 3, 4, 8, 16, 32]);
1184
1185 fn random_in_range<const N: usize>() {
1186 // Loosely verify we are drawing from the correct distribution
1187 let mut rng = StdRng::seed_from_u64(1);
1188 let a: Minkowski<N> = rng.random();
1189
1190 assert!(a.coordinates.iter().all(|&x| -1.0 < x && x < 1.0));
1191
1192 // This test will fail ~1e-3008 percent of the time - it's probably fine
1193 if N == 10_000 {
1194 assert!(a.coordinates.iter().any(|&x| x < 0.0));
1195 }
1196 }
1197
1198 parameterize_vector_length!(random_in_range, [2, 3, 4, 8, 16, 32, 10_000]);
1199
1200 /// Generate a pair of points in 2-dimensional hyperbolic space
1201 fn generate_h2_pair() -> (Hyperbolic<3>, Hyperbolic<3>) {
1202 (
1203 Hyperbolic::<3>::from_polar_coordinates(3.2, 0.1),
1204 Hyperbolic::<3>::from_polar_coordinates(4.0, 3.1),
1205 )
1206 }
1207 /// Generate a pair of points in 3-dimensional hyperbolic space
1208 fn generate_h3_pair() -> (Hyperbolic<4>, Hyperbolic<4>) {
1209 (
1210 Hyperbolic::<4>::from_polar_coordinates(3.5, 0.4, 0.5),
1211 Hyperbolic::<4>::from_polar_coordinates(4.2, 2.7, 0.1),
1212 )
1213 }
1214
1215 #[test]
1216 fn hyperbolic_distance() {
1217 let (a, b) = generate_h2_pair();
1218 let ab_distance = a.distance(&b);
1219 let ab_numeric_answer = 7.194_993_724_795_472;
1220 assert_relative_eq!(ab_distance, ab_numeric_answer, epsilon = 1e-12);
1221
1222 let (c, d) = generate_h3_pair();
1223 let cd_distance = c.distance(&d);
1224 let cd_numeric_answer = 7.525_514_513_583_905;
1225 assert_relative_eq!(cd_distance, cd_numeric_answer, epsilon = 1e-12);
1226 }
1227
1228 #[test]
1229 fn poincare_projection() {
1230 let a = Hyperbolic::<3>::from_polar_coordinates(1.5, 1.5);
1231 let a_poincare = a.to_poincare();
1232 let a_numeric_poincare = [0.044_928_659_534_049_77, 0.633_557_895_753_136_3];
1233 assert_relative_eq![a_poincare[0], a_numeric_poincare[0], epsilon = 1e-12];
1234 assert_relative_eq![a_poincare[1], a_numeric_poincare[1], epsilon = 1e-12];
1235
1236 let b = Hyperbolic::<3>::from_polar_coordinates(0.5, 4.2);
1237 let b_poincare = b.to_poincare();
1238 let b_numeric_poincare = [-0.120_074_024_591_707_93, -0.213_465_172_363_015_63];
1239 assert_relative_eq![b_poincare[0], b_numeric_poincare[0], epsilon = 1e-12];
1240 assert_relative_eq![b_poincare[1], b_numeric_poincare[1], epsilon = 1e-12];
1241 }
1242
1243 #[test]
1244 fn specific_distances() {
1245 // Distance to the cusp
1246 let a = Hyperbolic::<3>::from_polar_coordinates(1.2, 3.2);
1247 let a_cusp_distance = a.distance_from_cusp();
1248 let a_cusp_distance_numeric = 1.2;
1249 assert_relative_eq!(a_cusp_distance, a_cusp_distance_numeric, epsilon = 1e-12);
1250
1251 let c = Hyperbolic::<4>::from_polar_coordinates(1.2, 3.2, 1.2);
1252 let c_cusp_distance = c.distance_from_cusp();
1253 let c_cusp_distance_numeric = 1.2;
1254 assert_relative_eq!(c_cusp_distance, c_cusp_distance_numeric, epsilon = 1e-12);
1255 }
1256
1257 #[test]
1258 fn random_hyperbolic() {
1259 // Generate ten random points on the Hyperbolic
1260 let mut rng = StdRng::seed_from_u64(42);
1261 let d = 0.1;
1262 let origin = Minkowski::from([0.0, 0.0, 1.0]);
1263 for _n in 0..10 {
1264 let disk = HyperbolicDisk {
1265 disk_radius: d.try_into().expect("hard-coded positive number"),
1266 point: Hyperbolic::<3>::from_minkowski_coordinates(origin),
1267 };
1268 let random_point: Hyperbolic<3> = disk.sample(&mut rng);
1269
1270 // check that points remain on Hyperbolic
1271 let rho = -random_point
1272 .point
1273 .distance_squared(&Minkowski::<3>::default());
1274 assert_relative_eq!(rho, 1.0, epsilon = 1e-12);
1275
1276 // check that points are within distance d of cusp
1277 let distance = random_point.distance_from_cusp();
1278 assert!(d > distance);
1279 }
1280 }
1281}