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 two-sheeted 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. For
473/// example, rescaling distances by $`r \rightarrow r/\ell`$ has the same effect
474/// as running simulations with skirt length $`R = 1/\ell`$ or Gauss curvature
475/// $`K = -\ell^2`$.
476///
477/// [`Hyperbolic`] implements a [`Metric`] that computes the distance of the
478/// geodesic passing between two points on a hyperboloid with some given skirt
479/// width.
480///
481/// Two points on the hyperboloid with skirt width $` R = 1.0 `$:
482/// ```
483/// use hoomd_manifold::{Hyperbolic, Minkowski};
484/// use hoomd_vector::Metric;
485///
486/// let x = Hyperbolic::from_minkowski_coordinates([0.0, 0.0, 1.0].into());
487///
488/// let y = Hyperbolic::from_minkowski_coordinates(
489/// [0.0, 1.0, (2.0_f64).sqrt()].into(),
490/// );
491///
492/// assert_eq!(((2.0_f64).sqrt()).acosh(), x.distance(&y));
493/// ```
494///
495/// ## Remark on Numeric Stability
496///
497/// The hyperboloid model of hyperbolic space performs best when points are
498/// located nearby the cusp (i.e., Minkowski coordinates $`(0,0,1)`$).
499/// Because of the way trial moves are validated, Monte Carlo simulations
500/// become numerically unstable when any bodies are far away from the cusp.
501/// Ideally, simulations should keep bodies within a distance of $`\approx 8.0`$
502/// from the cusp. All of the boundary conditions implemented for `Hyperbolic`
503/// satisfy this condition and perform well in MC simulations.
504///
505/// Nevertheless, it is often necessary to run larger simulations. `Hyperbolic`
506/// methods will still run simulations with bodies located far away from the
507/// cusp, but with a reduced numerical tolerance. Note also that using too
508/// small of a maximum translation step size is inherently unstable.
509/// Simulations will still run if the specified maximum step size is too small,
510/// but `Hyperbolic` methods would no longer guarantee that translation moves
511/// are not translated more than the specified maximum.
512///
513/// A rough method selecting maximum step sizes: If $`r_{max}`$ is the
514/// distance between the cusp and the most distant body, avoid setting the
515/// maximum step size to be smaller than
516/// $`10^{\left(\frac{1}{2}\log_{10}\left[\cosh^2(r_{max})/2\right]-6\right)}`$. For example,
517/// if the farthest body has Minkowski coordinates
518/// $`(10^5, 10^5, \sqrt{2\cdot 10^{10} + 1})`$, then $`r_{max} = 12.5526`$ and
519/// therefore step sizes smaller than $`0.1`$ are unstable.
520///
521/// [`Metric`]: hoomd_vector::Metric
522#[derive(Clone, Copy, Debug, PartialEq, RelativeEq, Serialize, Deserialize)]
523pub struct Hyperbolic<const N: usize> {
524 /// A point on the surface of the upper sheet of a two-sheeted hyperboloid.
525 point: Minkowski<N>,
526}
527
528impl<const N: usize> Hyperbolic<N> {
529 /// Get the coordinates of the point on the hyperboloid.
530 #[must_use]
531 #[inline]
532 pub fn coordinates(&self) -> &[f64; N] {
533 &self.point.coordinates
534 }
535 /// Get the Minkowski point of the hyperboloid.
536 #[must_use]
537 #[inline]
538 pub fn point(&self) -> &Minkowski<N> {
539 &self.point
540 }
541 /// Create a Hyperbolic point from a Minkowski vector.
542 ///
543 /// # Example
544 /// ```
545 /// use hoomd_manifold::{Hyperbolic, Minkowski};
546 /// use hoomd_vector::Metric;
547 ///
548 /// let x = Hyperbolic::from_minkowski_coordinates([0.0, 0.0, 1.0].into());
549 /// ```
550 #[must_use]
551 #[inline]
552 #[allow(
553 clippy::cast_possible_truncation,
554 reason = "truncation to relax precision"
555 )]
556 pub fn from_minkowski_coordinates(point: Minkowski<N>) -> Hyperbolic<N> {
557 let pt = point.coordinates;
558 let min_coord_index = pt[0..N - 1]
559 .iter()
560 .enumerate()
561 .min_by(|(_, a), (_, b)| a.total_cmp(b))
562 .map_or(0, |(i, _)| i);
563 let min_coord = pt[min_coord_index];
564 let lhs = min_coord.mul_add(min_coord, 1.0);
565 let rhs = pt[0..N - 1]
566 .iter()
567 .enumerate()
568 .filter(|(j, _)| *j != min_coord_index)
569 .fold(pt[N - 1] * pt[N - 1], |sum, (_, x)| {
570 f64::mul_add(-x, *x, sum)
571 });
572 let tolerance = 9_i32 - 2_i32 * (point.coordinates[N - 1].log10().trunc() as i32);
573 assert_relative_eq!(lhs, rhs, epsilon = 10.0_f64.powi(-tolerance));
574 Hyperbolic { point }
575 }
576}
577
578impl Hyperbolic<3> {
579 /// Create a point on the surface of a three-dimensional hyperboloid from the polar representation.
580 #[must_use]
581 #[inline]
582 pub fn from_polar_coordinates(v: f64, theta: f64) -> Hyperbolic<3> {
583 let theta_mod = theta.rem_euclid(2.0 * PI);
584 let v_sinh = v.sinh();
585 let point = Minkowski::from([
586 v_sinh * (theta_mod.cos()),
587 v_sinh * (theta_mod.sin()),
588 v.cosh(),
589 ]);
590 Hyperbolic { point }
591 }
592}
593
594impl Hyperbolic<4> {
595 /// Create a point on the surface of a four-dimensional hyperboloid from the spherical representation.
596 #[must_use]
597 #[inline]
598 pub fn from_polar_coordinates(v: f64, theta: f64, phi: f64) -> Hyperbolic<4> {
599 let theta_mod = theta.rem_euclid(2.0 * PI);
600 let phi_mod = phi.rem_euclid(PI);
601 let v_sinh = v.sinh();
602 let point = Minkowski::from([
603 v_sinh * (theta_mod.cos()),
604 v_sinh * (theta_mod.sin()) * (phi_mod.cos()),
605 v_sinh * (theta_mod.sin()) * (phi_mod.sin()),
606 v.cosh(),
607 ]);
608 Hyperbolic { point }
609 }
610}
611
612impl<const N: usize> Hyperbolic<N> {
613 /// Compute the distance from a point to the cusp.
614 ///
615 /// Computes the length of the geodesic passing between the cusp
616 /// $`(0,\cdots,0,\rho)`$ and a given point on the hyperboloid.
617 ///
618 /// # Example
619 /// ```
620 /// use approxim::assert_relative_eq;
621 /// use hoomd_manifold::{Hyperbolic, Minkowski};
622 /// use hoomd_vector::Vector;
623 ///
624 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
625 /// let v: f64 = 4.2;
626 /// let x = Hyperbolic::from_minkowski_coordinates(
627 /// [v.sinh(), 0.0, v.cosh()].into(),
628 /// );
629 /// assert_relative_eq!(v, x.distance_from_cusp(), epsilon = 1e-12);
630 /// # Ok(())
631 /// # }
632 /// ```
633 #[inline]
634 #[must_use]
635 pub fn distance_from_cusp(&self) -> f64 {
636 (self.point.coordinates[N - 1]).acosh()
637 }
638
639 /// Projects points on the hyperboloid onto the Poincaré disk/ball.
640 ///
641 /// # Example
642 /// ```
643 /// use approxim::assert_relative_eq;
644 /// use hoomd_manifold::{Hyperbolic, Minkowski};
645 /// use hoomd_vector::Vector;
646 ///
647 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
648 /// let v: f64 = 1.098612;
649 /// let x = Hyperbolic::from_minkowski_coordinates(
650 /// [v.sinh(), 0.0, v.cosh()].into(),
651 /// );
652 /// let projection = x.to_poincare();
653 /// assert_relative_eq!(
654 /// v.sinh() / (v.cosh() + 1.0),
655 /// projection[0],
656 /// epsilon = 1e-12
657 /// );
658 /// assert_relative_eq!(0.0, projection[1], epsilon = 1e-12);
659 /// # Ok(())
660 /// # }
661 /// ```
662 #[inline]
663 #[must_use]
664 pub fn to_poincare(&self) -> Vec<f64> {
665 (0..N - 1)
666 .collect::<Vec<usize>>()
667 .iter()
668 .map(|i| self.point.coordinates[*i] / (1.0 + self.point.coordinates[N - 1]))
669 .collect::<Vec<f64>>()
670 }
671}
672
673impl<const N: usize> Default for Hyperbolic<N> {
674 /// Construct a default point on a hyperboloid.
675 ///
676 /// The default `Hyperbolic<N>` point is on the cusp of a hyperboloid with
677 /// skirt width of 1.0 (i.e., the point $`(0.0, \cdots, 0.0, 1.0)`$).
678 #[inline]
679 fn default() -> Self {
680 let mut zero = Minkowski::<N>::default();
681 zero.coordinates[N - 1] = 1.0;
682 Hyperbolic { point: zero }
683 }
684}
685
686impl Metric for Hyperbolic<3> {
687 /// The distance between two [`Hyperbolic<3>`] points.
688 ///
689 /// Explicitly, the metric for two points $`\vec{u}`$ and $`\vec{v}`$ on a
690 /// hyperboloid
691 ///
692 /// ```math
693 /// d_{H_2}(\vec{u}, \vec{v}) = \operatorname{arccosh}\left[u_3v_3 - u_1v_1 - u_2v_2\right]
694 /// ```
695 /// This choice of metric furnishes a representation of 2-dimensional hyperbolic
696 /// space with Gaussian curvature $`K = -1`$.
697 #[inline(always)]
698 fn distance(&self, other: &Self) -> f64 {
699 let self_coords = self.point.coordinates;
700 let other_coords = other.point.coordinates;
701 let arg = self_coords[2] * other_coords[2]
702 - self_coords[0] * other_coords[0]
703 - self_coords[1] * other_coords[1];
704 let arg_clipped = if arg <= 1.0 { 1.0 } else { arg };
705 arg_clipped.acosh()
706 }
707
708 #[inline]
709 fn distance_squared(&self, other: &Self) -> f64 {
710 self.distance(other).powi(2)
711 }
712
713 #[inline]
714 fn n_dimensions(&self) -> usize {
715 2_usize
716 }
717}
718
719impl Metric for Hyperbolic<4> {
720 /// The distance between two [`Hyperbolic<4>`] points.
721 ///
722 /// Explicitly, the metric for two points $`\vec{u}`$ and $`\vec{v}`$ on a
723 /// hyperboloid is given by
724 ///
725 /// ```math
726 /// d_{H_3}(\vec{u}, \vec{v}) = \operatorname{arccosh}\left[u_4v_4 - u_1v_1 - u_2v_2 - u_3v_3\right]
727 /// ```
728 /// This choice of metric furnishes a representation of 3-dimensional hyperbolic
729 /// space with with Gaussian curvature $`K = -1`$.
730 #[inline]
731 fn distance(&self, other: &Self) -> f64 {
732 let self_coords = self.point.coordinates;
733 let other_coords = other.point.coordinates;
734 let arg = self_coords[3] * other_coords[3]
735 - self_coords[0] * other_coords[0]
736 - self_coords[1] * other_coords[1]
737 - self_coords[2] * other_coords[2];
738 let arg_clipped = if arg <= 1.0 { 1.0 } else { arg };
739 arg_clipped.acosh()
740 }
741
742 #[inline]
743 fn distance_squared(&self, other: &Self) -> f64 {
744 self.distance(other).powi(2)
745 }
746
747 #[inline]
748 fn n_dimensions(&self) -> usize {
749 3_usize
750 }
751}
752
753/// Hyperbolic rotations in Minkowski Space
754///
755/// Construct a [`HyperbolicRotationMatrix`] to apply SO(N-1, 1)
756/// transformations to N-dimensional Minkowski vectors. For Minkowski 4-vectors,
757/// [`Biquaternion`] should be used instead for numerical stability. See
758/// documentation in [`HyperbolicAngle`] for details on SO(2,1) transformations
759/// (i.e., two-dimensional hyperbolic space), and in [`Biquaternion`] for
760/// details on SO(3,1) transformations (i.e., three-dimensional hyperbolic
761/// space).
762///
763/// [`Biquaternion`]: crate::Biquaternion
764/// [`HyperbolicAngle`]: crate::HyperbolicAngle
765///
766/// In two dimensional hyperbolic space:
767/// ```
768/// use hoomd_manifold::{
769/// HyperbolicAngle, HyperbolicRotate, HyperbolicRotationMatrix, Minkowski,
770/// };
771/// use std::f64::consts::PI;
772///
773/// fn rotate_about_z(minkowski_vector: &Minkowski<3>) -> Minkowski<3> {
774/// let generators = HyperbolicAngle::from((PI, 0.0_f64, 0.0_f64));
775/// let rotation_matrix = HyperbolicRotationMatrix::from(generators);
776/// rotation_matrix.hyperbolic_rotate(&minkowski_vector)
777/// }
778///
779/// fn boost_in_x(minkowski_vector: &Minkowski<3>) -> Minkowski<3> {
780/// let generators = HyperbolicAngle::from((0.0_f64, 0.2_f64, 0.0_f64));
781/// let boost_matrix = HyperbolicRotationMatrix::from(generators);
782/// boost_matrix.hyperbolic_rotate(&minkowski_vector)
783/// }
784/// ```
785#[serde_as]
786#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
787pub struct HyperbolicRotationMatrix<const N: usize> {
788 /// Rows of the rotation matrix.
789 #[serde_as(as = "[_; N]")]
790 pub(crate) rows: [Minkowski<N>; N],
791}
792
793impl<const N: usize> HyperbolicRotate<Minkowski<N>> for HyperbolicRotationMatrix<N> {
794 type Matrix = HyperbolicRotationMatrix<N>;
795
796 #[inline]
797 /// Rotate a [`Minkowski<N>`] by a [`HyperbolicRotationMatrix`]
798 ///
799 /// # Examples
800 ///
801 /// Rotate point in 2D hyperbolic space about z-axis:
802 /// ```
803 /// use hoomd_manifold::{
804 /// HyperbolicAngle, HyperbolicRotate, HyperbolicRotationMatrix, Minkowski,
805 /// };
806 /// use std::f64::consts::PI;
807 ///
808 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
809 /// let v = Minkowski::from([1.0, 0.0, 1.0]);
810 /// let spatial_rotation = HyperbolicAngle::from((PI / 2.0, 0.0_f64, 0.0_f64));
811 /// let matrix = HyperbolicRotationMatrix::from(spatial_rotation);
812 /// let rotated = matrix.hyperbolic_rotate(&v);
813 /// let c = Minkowski::from([(PI / 2.0).cos(), (PI / 2.0).sin(), 1.0]);
814 /// assert_eq!(c, rotated);
815 /// # Ok(())
816 /// # }
817 /// ```
818 ///
819 /// Boost point in 2D hyperbolic space in x direction:
820 /// ```
821 /// use hoomd_manifold::{
822 /// HyperbolicAngle, HyperbolicRotate, HyperbolicRotationMatrix, Minkowski,
823 /// };
824 /// use num::complex::Complex;
825 /// use std::f64::consts::PI;
826 ///
827 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
828 /// let v = Minkowski::from([1.0, 0.0, 1.0]);
829 /// let small_boost = HyperbolicAngle::from((0.0_f64, 0.1_f64, 0.0_f64));
830 /// let matrix = HyperbolicRotationMatrix::from(small_boost);
831 /// let rotated = matrix.hyperbolic_rotate(&v);
832 /// let c = Minkowski::from([
833 /// (0.1_f64).sinh() + (0.1_f64).cosh(),
834 /// 0.0,
835 /// (0.1_f64).sinh() + (0.1_f64).cosh(),
836 /// ]);
837 /// assert_eq!(c, rotated);
838 /// # Ok(())
839 /// # }
840 /// ```
841 ///
842 /// Zero angles and rapidities does nothing:
843 /// ```
844 /// use hoomd_manifold::{
845 /// HyperbolicAngle, HyperbolicRotate, HyperbolicRotationMatrix, Minkowski,
846 /// };
847 /// use num::complex::Complex;
848 /// use std::f64::consts::PI;
849 ///
850 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
851 /// let v = Minkowski::from([1.0, 2.0, 1.0]);
852 /// let identity = HyperbolicAngle::from((0.0_f64, 0.0_f64, 0.0_f64));
853 /// let matrix = HyperbolicRotationMatrix::from(identity);
854 /// let rotated = matrix.hyperbolic_rotate(&v);
855 /// let c = Minkowski::from([1.0, 2.0, 1.0]);
856 /// assert_eq!(c, rotated);
857 /// # Ok(())
858 /// # }
859 /// ```
860 ///
861 /// Rotate point in 3D hyperbolic space about y axis using matrix representation:
862 /// ```
863 /// use approxim::assert_relative_eq;
864 /// use hoomd_manifold::{
865 /// Biquaternion, HyperbolicRotate, HyperbolicRotationMatrix, Minkowski,
866 /// UnitBiquaternion,
867 /// };
868 /// use num::complex::Complex;
869 /// use std::f64::consts::PI;
870 ///
871 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
872 /// let q = Biquaternion::from([
873 /// Complex::new(0.0, 0.0),
874 /// Complex::new((PI / 4.0).sin(), 0.0),
875 /// Complex::new(0.0, 0.0),
876 /// Complex::new((PI / 4.0).cos(), 0.0),
877 /// ]);
878 /// let v = q.to_unit()?;
879 /// let x = Minkowski::from([1.0, 0.0, 0.0, 1.0]);
880 /// let rotation = HyperbolicRotationMatrix::from(v);
881 /// let rotated = rotation.hyperbolic_rotate(&x);
882 /// assert_relative_eq!(rotated, [0.0, 0.0, -1.0, 1.0].into(), epsilon = 1e-12);
883 /// # Ok(())
884 /// # }
885 /// ```
886 ///
887 /// Boost point in 3D hyperbolic space in x direction using biquaternion algebra:
888 /// ```
889 /// use approxim::assert_relative_eq;
890 /// use hoomd_manifold::{
891 /// Biquaternion, HyperbolicRotate, Minkowski, UnitBiquaternion,
892 /// };
893 /// use num::complex::Complex;
894 /// use std::f64::consts::PI;
895 ///
896 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
897 /// let x = Minkowski::from([0.0, 0.0, 0.0, 1.0]);
898 /// let q = Biquaternion::from([
899 /// Complex::new(0.0, 0.25).sin(),
900 /// Complex::new(0.0, 0.0),
901 /// Complex::new(0.0, 0.0),
902 /// Complex::new(0.0, 0.25).cos(),
903 /// ]);
904 /// let v = q.to_unit()?;
905 /// let boosted = v.hyperbolic_rotate(&x);
906 /// assert_relative_eq!(
907 /// boosted,
908 /// [(0.5_f64).sinh(), 0.0, 0.0, (0.5_f64).cosh()].into(),
909 /// epsilon = 1e-12
910 /// );
911 /// # Ok(())
912 /// # }
913 /// ```
914 fn hyperbolic_rotate(&self, vector: &Minkowski<N>) -> Minkowski<N> {
915 let mut coordinates = [0.0; N];
916
917 for (result, row) in coordinates.iter_mut().zip(self.rows.iter()) {
918 *result = zip(row.coordinates.iter(), vector.coordinates.iter())
919 .fold(0.0, |product, x| product + x.0 * x.1);
920 }
921
922 Minkowski { coordinates }
923 }
924}
925
926/// Randomly distribute points locally on a hyperboloid.
927///
928/// [`HyperbolicDisk`] is a uniform distribution of points within distance `r`
929/// of a point on the 2-dimensional hyperboloid.
930///
931/// # Example
932///
933/// ```
934/// use hoomd_manifold::{
935/// Hyperbolic, HyperbolicAngle, HyperbolicDisk, HyperbolicRotate,
936/// HyperbolicRotationMatrix, Minkowski,
937/// };
938/// use hoomd_vector::Metric;
939/// use rand::{RngExt, SeedableRng, distr::Distribution, rngs::StdRng};
940///
941/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
942/// let mut rng = StdRng::seed_from_u64(12);
943///
944/// let v: HyperbolicAngle = rng.random();
945/// let matrix = HyperbolicRotationMatrix::from(v);
946/// let origin = Minkowski::from([0.0, 0.0, 1.0]);
947/// let random_point = Hyperbolic::from_minkowski_coordinates(
948/// matrix.hyperbolic_rotate(&origin),
949/// );
950///
951/// let r = 0.1;
952/// let mut rng_2 = StdRng::seed_from_u64(239);
953/// let disk = HyperbolicDisk {
954/// disk_radius: r.try_into()?,
955/// point: random_point,
956/// };
957/// let transformed_random_point: Hyperbolic<3> = disk.sample(&mut rng_2);
958///
959/// assert!(r > random_point.distance(&transformed_random_point));
960///
961/// # Ok(())
962/// # }
963/// ```
964#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
965pub struct HyperbolicDisk {
966 /// Max distance away from point.
967 pub disk_radius: PositiveReal,
968 /// The center of the disk.
969 pub point: Hyperbolic<3>,
970}
971
972impl Distribution<Hyperbolic<3>> for HyperbolicDisk {
973 /// Sample a random point in the hyperbolic disk.
974 ///
975 /// The implementation translates Minkowski 3-vector `point` along
976 /// the Hyperbolic by maximum distance of `disk_radius`. Note that because SO(2,1) is
977 /// non-Abelian, the point must be transformed to the cusp before a trial
978 /// move is applied (and then the point is transformed back). This ensures
979 /// that the max distance translated by the trial move does not exceed `disk_radius`.
980 #[inline]
981 fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> Hyperbolic<3> {
982 let max_boost = self.disk_radius.get();
983 let point = self.point;
984 let eta = (point.point.coordinates[2]).acosh();
985 let phi = point.point.coordinates[1].atan2(point.point.coordinates[0]);
986 let trial_boost = Uniform::new(0.0, 1.0).expect("r is positive and real");
987 let trial_rotation =
988 Uniform::new(-PI, PI).expect("hard-coded distribution should be valid");
989 let theta = trial_rotation.sample(rng);
990 let v1: f64 = trial_boost.sample(rng);
991 let v = v1.sqrt() * max_boost;
992 let (v_sinh, eta_sinh, eta_cosh, phi_sin, phi_cos) =
993 (v.sinh(), eta.sinh(), eta.cosh(), phi.sin(), phi.cos());
994 let trial_coords = [v_sinh * theta.cos(), v_sinh * theta.sin(), v.cosh()];
995 let transformed_point = Minkowski::from([
996 trial_coords[0] * (eta_cosh * (phi_cos.powi(2)) + phi_sin.powi(2))
997 + trial_coords[1] * phi_sin * phi_cos * (eta_cosh - 1.0)
998 + trial_coords[2] * eta_sinh * phi_cos,
999 trial_coords[0] * phi_sin * phi_cos * (eta_cosh - 1.0)
1000 + trial_coords[1] * (eta_cosh * (phi_sin.powi(2)) + phi_cos.powi(2))
1001 + trial_coords[2] * eta_sinh * phi_sin,
1002 trial_coords[0] * eta_sinh * phi_cos
1003 + trial_coords[1] * eta_sinh * phi_sin
1004 + trial_coords[2] * eta_cosh,
1005 ]);
1006 Hyperbolic::from_minkowski_coordinates(transformed_point)
1007 }
1008}
1009
1010#[cfg(test)]
1011mod tests {
1012 use super::*;
1013 use approxim::assert_relative_eq;
1014 use paste::paste;
1015 use rand::{RngExt, SeedableRng, rngs::StdRng};
1016
1017 // Parameterize a test function over an array of vector lengths
1018 macro_rules! parameterize_vector_length {
1019 // macro with name as above that takes an identifier (fn) and an expression
1020 // $(...),* matches 0 or more expressions (values) separated by commas
1021 ($test_body:ident, [$($dim:expr),*]) => {
1022
1023 // Now, we repeat the test block 0 or more times, one for each $dim
1024 $(
1025 // paste package combines values in [< >] to form a new ident
1026 paste! {
1027 #[test]
1028 fn [< $test_body "_" $dim>]() {
1029 const DIM: usize = $dim;
1030 $test_body::<DIM>();
1031 }
1032 }
1033 )*
1034 };
1035 }
1036
1037 /// Generate a pair of length N vectors.
1038 /// The first vector ranges from [0, N-1] and the second ranges from [N, 2*N-1]
1039 fn generate_vector_pair<const N: usize>() -> (Minkowski<N>, Minkowski<N>) {
1040 (
1041 Minkowski::try_from(0..N).unwrap(),
1042 Minkowski::try_from(N..N * 2).unwrap(),
1043 )
1044 }
1045
1046 fn index<const N: usize>() {
1047 let (_, b) = generate_vector_pair::<N>();
1048 assert!(zip(0..N, b.coordinates.iter()).all(|(i, &x)| b[i] == x));
1049 }
1050 parameterize_vector_length!(index, [2, 3, 4, 8, 16, 32]);
1051
1052 fn index_mut<const N: usize>() {
1053 let (a, mut b) = generate_vector_pair::<N>();
1054 zip(0..N, b.coordinates.iter_mut()).for_each(|(i, x)| *x = a[i]);
1055 assert_eq!(a, b);
1056 }
1057 parameterize_vector_length!(index_mut, [2, 3, 4, 8, 16, 32]);
1058
1059 fn add_explicit<const N: usize>() {
1060 let (a, b) = generate_vector_pair::<N>();
1061 let c = a.add(b);
1062
1063 let addition_answer: Vec<f64> = (0..(2 * N))
1064 .step_by(2)
1065 .map(|x| (x + N) as f64)
1066 .collect::<Vec<_>>();
1067
1068 assert_eq!(c, Minkowski::try_from(addition_answer).unwrap());
1069 }
1070 parameterize_vector_length!(add_explicit, [2, 3, 4, 8, 16, 32]);
1071
1072 fn add_operator<const N: usize>() {
1073 let (a, b) = generate_vector_pair::<N>();
1074 let c = a + b;
1075
1076 let addition_answer: Vec<f64> = (0..(2 * N))
1077 .step_by(2)
1078 .map(|x| (x + N) as f64)
1079 .collect::<Vec<_>>();
1080
1081 assert_eq!(c, Minkowski::try_from(addition_answer).unwrap());
1082 }
1083 parameterize_vector_length!(add_operator, [2, 3, 4, 8, 16, 32]);
1084
1085 fn add_assign<const N: usize>() {
1086 let (a, b) = generate_vector_pair::<N>();
1087 let mut c = a;
1088 c += b;
1089
1090 let addition_answer: Vec<f64> = (0..(2 * N))
1091 .step_by(2)
1092 .map(|x| (x + N) as f64)
1093 .collect::<Vec<_>>();
1094
1095 assert_eq!(c, Minkowski::try_from(addition_answer).unwrap());
1096 }
1097 parameterize_vector_length!(add_assign, [2, 3, 4, 8, 16, 32]);
1098
1099 fn sub_operator<const N: usize>() {
1100 let (a, b) = generate_vector_pair::<N>();
1101 let c = a - b;
1102
1103 let subtraction_answer = [-(N as f64); N];
1104
1105 assert_eq!(c, subtraction_answer.into());
1106 }
1107 parameterize_vector_length!(sub_operator, [2, 3, 4, 8, 16, 32]);
1108
1109 fn sub_assign<const N: usize>() {
1110 let (a, b) = generate_vector_pair::<N>();
1111 let mut c = a;
1112 c -= b;
1113
1114 let subtraction_answer = [-(N as f64); N];
1115
1116 assert_eq!(c, subtraction_answer.into());
1117 }
1118
1119 parameterize_vector_length!(sub_assign, [2, 3, 4, 8, 16, 32]);
1120
1121 fn mul_operator<const N: usize>() {
1122 let (a, _) = generate_vector_pair::<N>();
1123 let b = 12.0;
1124 let c = a * b;
1125
1126 let multiplication_answer: Vec<f64> = (0..N).map(|x| (x as f64) * b).collect::<Vec<_>>();
1127
1128 assert_eq!(c, Minkowski::try_from(multiplication_answer).unwrap());
1129 }
1130 parameterize_vector_length!(mul_operator, [2, 3, 4, 8, 16, 32]);
1131
1132 fn div_operator<const N: usize>() {
1133 let (a, _) = generate_vector_pair::<N>();
1134 let b = 12.0;
1135 let c = a / b;
1136
1137 let division_answer: Vec<f64> = (0..N).map(|x| (x as f64) / b).collect::<Vec<_>>();
1138
1139 assert_eq!(c, Minkowski::try_from(division_answer).unwrap());
1140 }
1141 parameterize_vector_length!(div_operator, [2, 3, 4, 8, 16, 32]);
1142
1143 fn div_assign<const N: usize>() {
1144 let (mut a, _) = generate_vector_pair::<N>();
1145 let b = 12.0;
1146 a /= b;
1147
1148 let division_answer: Vec<f64> = (0..N).map(|x| (x as f64) / b).collect::<Vec<_>>();
1149
1150 assert_eq!(a, Minkowski::try_from(division_answer).unwrap());
1151 }
1152
1153 parameterize_vector_length!(div_assign, [2, 3, 4, 8, 16, 32]);
1154
1155 fn compute_add_ref_ref<const N: usize>(a: &Minkowski<N>, b: &Minkowski<N>) -> Minkowski<N> {
1156 *a + *b
1157 }
1158
1159 fn compute_add_ref_type<const N: usize>(a: &Minkowski<N>, b: Minkowski<N>) -> Minkowski<N> {
1160 *a + b
1161 }
1162
1163 fn compute_add_type_ref<const N: usize>(a: Minkowski<N>, b: &Minkowski<N>) -> Minkowski<N> {
1164 a + *b
1165 }
1166
1167 fn add_with_refs<const N: usize>() {
1168 let (a, b) = generate_vector_pair::<N>();
1169
1170 let addition_answer = Minkowski::try_from(
1171 (0..(2 * N))
1172 .step_by(2)
1173 .map(|x| (x + N) as f64)
1174 .collect::<Vec<_>>(),
1175 )
1176 .unwrap();
1177
1178 let c = compute_add_ref_ref(&a, &b);
1179 assert_eq!(c, addition_answer);
1180
1181 let c = compute_add_ref_type(&a, b);
1182 assert_eq!(c, addition_answer);
1183
1184 let c = compute_add_type_ref(a, &b);
1185 assert_eq!(c, addition_answer);
1186 }
1187 parameterize_vector_length!(add_with_refs, [2, 3, 4, 8, 16, 32]);
1188
1189 #[test]
1190 fn display() {
1191 let a = Minkowski::from([1.67, -2.125, 42.01]);
1192 let s = format!("{a}");
1193 assert_eq!(s, "[1.67, -2.125, 42.01]");
1194
1195 let a = Minkowski::from([10.0, 20.0, 30.0, 40.0]);
1196 let s = format!("{a}");
1197 assert_eq!(s, "[10, 20, 30, 40]");
1198 }
1199
1200 #[test]
1201 fn from_2_tuple() {
1202 let a = Minkowski::from((13.0, 0.125));
1203 assert_eq!(a.coordinates, [13.0, 0.125]);
1204 }
1205
1206 #[test]
1207 fn from_3_tuple() {
1208 let a = Minkowski::from((-0.5, 2.0, 18.125));
1209 assert_eq!(a.coordinates, [-0.5, 2.0, 18.125]);
1210 }
1211
1212 fn from_vec<const N: usize>() {
1213 let mut vec = Vec::with_capacity(N);
1214
1215 assert_eq!(
1216 Minkowski::<N>::try_from(vec.clone()),
1217 Err(Error::InvalidVectorLength)
1218 );
1219
1220 for i in 0..N {
1221 vec.push(i as f64 * 0.5);
1222 }
1223 let a = Minkowski::<N>::try_from(vec.clone()).unwrap();
1224
1225 assert_eq!(vec, Vec::from(a.coordinates));
1226
1227 vec.push(1.0);
1228 assert_eq!(
1229 Minkowski::<N>::try_from(vec.clone()),
1230 Err(Error::InvalidVectorLength)
1231 );
1232 }
1233 parameterize_vector_length!(from_vec, [2, 3, 4, 8, 16, 32]);
1234
1235 fn random_in_range<const N: usize>() {
1236 // Loosely verify we are drawing from the correct distribution
1237 let mut rng = StdRng::seed_from_u64(1);
1238 let a: Minkowski<N> = rng.random();
1239
1240 assert!(a.coordinates.iter().all(|&x| -1.0 < x && x < 1.0));
1241
1242 // This test will fail ~1e-3008 percent of the time - it's probably fine
1243 if N == 10_000 {
1244 assert!(a.coordinates.iter().any(|&x| x < 0.0));
1245 }
1246 }
1247
1248 parameterize_vector_length!(random_in_range, [2, 3, 4, 8, 16, 32, 10_000]);
1249
1250 /// Generate a pair of points in 2-dimensional hyperbolic space
1251 fn generate_h2_pair() -> (Hyperbolic<3>, Hyperbolic<3>) {
1252 (
1253 Hyperbolic::<3>::from_polar_coordinates(3.2, 0.1),
1254 Hyperbolic::<3>::from_polar_coordinates(4.0, 3.1),
1255 )
1256 }
1257 /// Generate a pair of points in 3-dimensional hyperbolic space
1258 fn generate_h3_pair() -> (Hyperbolic<4>, Hyperbolic<4>) {
1259 (
1260 Hyperbolic::<4>::from_polar_coordinates(3.5, 0.4, 0.5),
1261 Hyperbolic::<4>::from_polar_coordinates(4.2, 2.7, 0.1),
1262 )
1263 }
1264
1265 #[test]
1266 fn hyperbolic_distance() {
1267 let (a, b) = generate_h2_pair();
1268 let ab_distance = a.distance(&b);
1269 let ab_numeric_answer = 7.194_993_724_795_472;
1270 assert_relative_eq!(ab_distance, ab_numeric_answer, epsilon = 1e-12);
1271
1272 let (c, d) = generate_h3_pair();
1273 let cd_distance = c.distance(&d);
1274 let cd_numeric_answer = 7.525_514_513_583_905;
1275 assert_relative_eq!(cd_distance, cd_numeric_answer, epsilon = 1e-12);
1276 }
1277
1278 #[test]
1279 fn poincare_projection() {
1280 let a = Hyperbolic::<3>::from_polar_coordinates(1.5, 1.5);
1281 let a_poincare = a.to_poincare();
1282 let a_numeric_poincare = [0.044_928_659_534_049_77, 0.633_557_895_753_136_3];
1283 assert_relative_eq![a_poincare[0], a_numeric_poincare[0], epsilon = 1e-12];
1284 assert_relative_eq![a_poincare[1], a_numeric_poincare[1], epsilon = 1e-12];
1285
1286 let b = Hyperbolic::<3>::from_polar_coordinates(0.5, 4.2);
1287 let b_poincare = b.to_poincare();
1288 let b_numeric_poincare = [-0.120_074_024_591_707_93, -0.213_465_172_363_015_63];
1289 assert_relative_eq![b_poincare[0], b_numeric_poincare[0], epsilon = 1e-12];
1290 assert_relative_eq![b_poincare[1], b_numeric_poincare[1], epsilon = 1e-12];
1291 }
1292
1293 #[test]
1294 fn specific_distances() {
1295 // Distance to the cusp
1296 let a = Hyperbolic::<3>::from_polar_coordinates(1.2, 3.2);
1297 let a_cusp_distance = a.distance_from_cusp();
1298 let a_cusp_distance_numeric = 1.2;
1299 assert_relative_eq!(a_cusp_distance, a_cusp_distance_numeric, epsilon = 1e-12);
1300
1301 let c = Hyperbolic::<4>::from_polar_coordinates(1.2, 3.2, 1.2);
1302 let c_cusp_distance = c.distance_from_cusp();
1303 let c_cusp_distance_numeric = 1.2;
1304 assert_relative_eq!(c_cusp_distance, c_cusp_distance_numeric, epsilon = 1e-12);
1305 }
1306
1307 #[test]
1308 fn random_hyperbolic() {
1309 // Generate ten random points on the Hyperbolic
1310 let mut rng = StdRng::seed_from_u64(42);
1311 let d = 0.1;
1312 let origin = Minkowski::from([0.0, 0.0, 1.0]);
1313 for _n in 0..10 {
1314 let disk = HyperbolicDisk {
1315 disk_radius: d.try_into().expect("hard-coded positive number"),
1316 point: Hyperbolic::<3>::from_minkowski_coordinates(origin),
1317 };
1318 let random_point: Hyperbolic<3> = disk.sample(&mut rng);
1319
1320 // check that points remain on Hyperbolic
1321 let rho = -random_point
1322 .point
1323 .distance_squared(&Minkowski::<3>::default());
1324 assert_relative_eq!(rho, 1.0, epsilon = 1e-12);
1325
1326 // check that points are within distance d of cusp
1327 let distance = random_point.distance_from_cusp();
1328 assert!(d > distance);
1329 }
1330 }
1331}