"""Hyperbolic geometry module.
Contains model classes and functions for conversion between the models.
"""
from functools import singledispatch
import numpy as np
from ddg.geometry import euclidean_models, moebius_models
from ddg.geometry.conversion import quadric_to_subspaces
from ddg.geometry.geometry_model_templates import (
MetricCayleyKleinGeometry,
MetricGeometry,
)
from ddg.geometry.intersection import join, meet
from ddg.geometry.projections import (
central_project,
inverse_stereographic_project,
lift_sphere_to_quadric,
stereographic_project,
)
from ddg.geometry.quadrics import Quadric, touching_cone
from ddg.geometry.spheres import SphereLike
from ddg.geometry.subspaces import (
Point,
Subspace,
coordinate_hyperplane,
subspace_from_columns,
)
__all__ = [
"ProjectiveModel",
"PoincareDiskModel",
"HalfSpaceModel",
"HemisphereModel",
"projective_to_poincare",
"poincare_to_projective",
"projective_to_hemisphere",
"hemisphere_to_projective",
"hemisphere_to_poincare",
"poincare_to_hemisphere",
"hemisphere_to_half_space",
"half_space_to_hemisphere",
]
[docs]class ProjectiveModel(MetricCayleyKleinGeometry):
"""Projective model of hyperbolic geometry, aka Klein model.
.. rubric:: Model space
The model space can be thought of as the
inside of the unit ball. More precisely, it is the set ::
{ [x] : b(x, x) < 0 }
where `b` is the standard Lorentz scalar product represented by the matrix
`diag([1,...,1, -1])`.
.. rubric:: Representation of objects
- Hyperbolic subspaces are represented by subspace objects
- Spheres, distance curves and horospheres are sphere objects.
Note that a hyperbolic subspace is a subspace intersected with the model
space, but we use the whole subspace to represent it.
Parameters
----------
dimension : int
Attributes
----------
dimension : int
Notes
-----
This class supports comparison with `==`. Two geometries are equal
if and only if they have the same type and `dimension`.
This class also supports the `in` operator to check whether a point is in
the model space.
"""
@property
def absolute(self) -> Quadric:
"""The absolute quadric with matrix ``diag([1,...,1, -1])``.
Returns
-------
Quadric
"""
matrix = np.eye(self.dimension + 1)
matrix[-1, -1] = -1
return Quadric(matrix)
[docs] @staticmethod
def metric_to_cayley_klein_distance(d: float) -> float:
"""Return Metric distance converted to Cayley-Klein "distance"::
K = cosh(d) ** 2
Parameters
----------
d : float
Returns
-------
K : float
"""
return np.cosh(d) ** 2
[docs] @staticmethod
def cayley_klein_distance_to_metric(K: float) -> float:
"""Return Metric distance converted to Cayley-Klein "distance"::
d = arccosh(sqrt(K))
Parameters
----------
K : float
Returns
-------
d : float
Raises
------
ValueError
If K < 1
"""
if K < 1:
raise ValueError(
"Cayley-Klein distance must be in [1, inf[ to correspond to a "
"metric distance in hyperbolic geometry."
)
return np.arccosh(np.sqrt(K))
[docs] def geodesic(self, p1, p2):
"""Return geodesic connecting two points.
Parameters
----------
p1, p2 : Point or array_like of shape (self.dimension + 1,)
Returns
-------
Subspace
"""
if isinstance(p1, Point):
p1 = p1.point
if isinstance(p2, Point):
p2 = p2.point
return Subspace(p1, p2).dehomogenize()
[docs] def perpendicular_subspace(self, subspace, contained_subspace):
"""Get a subspace that is perpendicular to another subspace.
More precisely: The function returns the largest subspace S such that the
orthogonal projection of S onto `subspace` is equal to the orthogonal projection
of `contained_subspace` onto `subspace`.
Common special case: If `subspace` is a hyperplane and `contained_subspace` is
a point, this is the line that is orthogonal to the hyperplane and contains a
specified point.
This is computed as the join of `contained_subspace` and the polar subspace of
`subspace`.
Parameters
----------
subspace, contained_subspace : Subspace
Returns
-------
S : Subspace
"""
return join(self.absolute.polarize(subspace), contained_subspace)
[docs] def common_perpendicular(self, subspace1, subspace2) -> Subspace:
"""Subspace that is orthogonal to two given non-intersecting subspaces.
Common special case: The line that is perpendicular to two parallel hyperplanes.
This is computed as the join of the polarized subspaces. If the two subspaces
intersect, this join will lie outside of the unit sphere.
Parameters
----------
subspace1, subspace2 : Subspace
Returns
-------
S : Subspace
"""
return join(
self.absolute.polarize(subspace1), self.absolute.polarize(subspace2)
)
[docs] def d_point_hyperplane(self, point, subspace) -> float:
"""Perpendicular distance of a point to a hyperplane.
Parameters
----------
point : Point
subspace : Subspace
Returns
-------
float
"""
ptp = self.perpendicular_subspace(subspace, point)
point_on_line = meet(subspace, ptp)
return self.d(point, point_on_line)
[docs] def d_hyperplanes(self, subspace1, subspace2) -> float:
"""Perpendicular distance between two parallel hyperplanes.
Parameters
----------
subspace1, subspace2 : Subspace
Returns
-------
float
"""
cp = self.common_perpendicular(subspace1, subspace2)
p1 = meet(cp, subspace1)
p2 = meet(cp, subspace2)
return self.d(p1, p2)
[docs] def angle(self, subspace1, subspace2) -> float:
"""Angle between two hyperplanes.
Parameters
----------
subspace1, subspace2 : Subspace
Returns
-------
float
"""
n1 = self.absolute.polarize(subspace1).point
n2 = self.absolute.polarize(subspace2).point
return np.arccos(np.sqrt(abs(self.cayley_klein_distance(n1, n2))))
def __contains__(self, point):
return self.inner_product(point, point) < 0
[docs]class HemisphereModel(moebius_models._ProjectiveSubgeometry):
"""Hemisphere model of hyperbolic space.
Also known as the Möbius model because this is a subgeometry of Möbius
geometry. It is also a transformed version of the hyperboloid model.
.. rubric:: Model space
The model space is the upper open hemisphere of the unit sphere. More
precisely::
{ [x] in Quadric(diag([1,...,1, -1])) : dehomogenize([x])[-1] > 0 }.
.. rubric:: Representation of objects
All geometric objects (except for points) are represented by Möbius
spheres, i.e. Quadric objects contained in the Möbius quadric. Note that of
course the actual objects are quadrics intersected with the upper
hemisphere, but we use the whole quadric to represent them.
Parameters
----------
dimension : int
Attributes
----------
dimension : int
Notes
-----
This class supports comparison with `==`. Two geometries are equal
if and only if they have the same type, and `dimension`.
This class also supports the `in` operator to check whether a point is in
the model space.
"""
@property
def fixed_point(self):
"""Point outside the Möbius quadric used for computation of the metric.
This is the point with representative [0,...,0, 1, 0].
Returns
-------
Point
"""
p = np.zeros(self.ambient_dimension + 1)
p[-2] = 1
return Point(p)
[docs] def d(self, v, w):
"""Hyperbolic metric in Möbius model
Given by the formula ::
<x, y> / (<x, p> <y, p>) = -2 * sinh^2(d([x],[y]) / 2)
Where `p` is `self.fixed_point.point`.
Parameters
----------
v, w : Point or array_like of shape (n,)
Returns
-------
float
"""
if isinstance(v, Point):
v = v.point
if isinstance(w, Point):
w = w.point
p = self.fixed_point.point
b = self.inner_product
return 2 * np.arcsinh(np.sqrt(-0.5 * b(v, w) / (b(v, p) * b(w, p))))
def __contains__(self, point):
return super().__contains__(point) and point.affine_point[-1] > 0
[docs]class HalfSpaceModel(MetricGeometry):
"""Poincaré Half-space model.
.. rubric:: Model space
The model space is the half-space {xn > 0} in R^n, where `n` is
the `dimension`.
.. rubric:: Representation of objects
Objects are represented by their Euclidean counterparts:
- Subspaces, including other objects like horospheres that look like
subspaces, are represented by subspaces.
- Objects that look like spheres are represented by Euclidean spheres. This
means that the center of the spheres will not match the actual hyperbolic
center.
Note that we represent intersections of the upper half-space with an object
`obj` by the whole `obj`. For example, geodesics are semicircles, but we
represent them by the whole circle.
Parameters
----------
dimension : int
Attributes
----------
dimension : int
Notes
-----
This class supports comparison with `==`. Two geometries are equal
if and only if they have the same type, and `dimension`.
This class also supports the `in` operator to check whether a point is in
the model space.
"""
[docs] def d(self, v, w):
"""Metric in Poincaré half-space model.
Defined by::
arccosh(1 + |v-w|^2 / (2 * v[n] * w[n])),
where `n` is the `dimension`.
Parameters
----------
v, w : Point or array_like of shape(n,)
Returns
-------
float
Raises
------
ValueError
If the points are in different half-spaces.
Notes
-----
Source: Wikipedia
"""
if not isinstance(v, Point):
v = Point(v)
if not isinstance(w, Point):
w = Point(w)
v = v.affine_point
w = w.affine_point
# If this is not true, we can't apply Arcosh.
if np.sign(v[-1]) != np.sign(w[-1]):
raise ValueError("Points are in different half-spaces.")
return np.arccosh(1 + np.dot(v - w, v - w) / (2 * v[-1] * w[-1]))
@property
def absolute(self):
"""The model embedded in (n+1)-space.
Returns the hyperplane {x0 = 0}. The model is the upper half of this hyperplane.
This plane is projected to the Möbius quadric to convert between the
models.
Returns
-------
Subspace
"""
return coordinate_hyperplane(self.dimension + 1, zero_component=0)
def __contains__(self, point):
return point.ambient_dimension == self.dimension and point.affine_point[-1] > 0
[docs]class PoincareDiskModel(MetricGeometry):
"""Poincaré disk model.
.. rubric:: Model space
The model space is the open unit ball in R^n.
.. rubric:: Representation of objects
Geometric objects are represented by their Euclidean equivalents, mainly
spheres. Note:
- If an object is the intersection of the unit ball with a larger object,
we use the whole object to represent it. For example, geodesics are
circular arcs, but we represent them by the whole circle.
- Because we use Euclidean spheres to represent hyperbolic spheres which
are the same viewed as sets, the center of a sphere will not actually
match the center with respect to the hyperbolic metric. The hyperbolic
center can be obtained by representing the sphere in a different model
and converting the center to this model.
Parameters
----------
dimension : int
Attributes
----------
dimension : int
Notes
-----
This class supports comparison with `==`. Two geometries are equal
if and only if they have the same type, and `dimension`.
This class also supports the `in` operator to check whether a point is in
the model space.
"""
[docs] def d(self, v, w):
"""Metric in the Poincaré disk model.
Parameters
----------
v, w : Point or array_like of shape (n,)
Returns
-------
float
"""
if not isinstance(v, Point):
v = Point(v)
if not isinstance(w, Point):
w = Point(w)
o = np.zeros(self.dimension + 1)
o[-1] = 1
o = Point(o)
euc = euclidean_models.ProjectiveModel(self.dimension)
numer = 2 * euc.d(v, w) ** 2
denom = (1 - euc.d(o, v) ** 2) * (1 - euc.d(o, w) ** 2)
return np.arccosh(1 + numer / denom)
def __contains__(self, point):
unit_ball = ProjectiveModel(self.dimension).absolute
return unit_ball.inner_product(point, point) < 0
[docs]@singledispatch
def projective_to_hemisphere(object_):
"""Projective to hemisphere model conversion.hyperconv
This is a dispatcher, see type-specific functions for details.
"""
raise TypeError(
f"Conversion of {type(object_)} objects from projective "
f"to Möbius model is not supported."
)
@projective_to_hemisphere.register
def projective_to_hemisphere_subspace(subspace: Subspace, embedded=False):
"""Convert Subspace from Projective model to Möbius/hemisphere model.
Works by optionally embedding `subspace` into the equatorial plane (the
polar plane of [e_n+1]), then projecting to the Möbius quadric via [e_n+1].
The Möbius model is only one hemisphere, so the image is a semicircle, but
this will return the whole circle.
Parameters
----------
subspace : Subspace
embedded : bool (default=False)
If False, object will be embedded into the equatorial plane before
projecting.
Returns
-------
Quadric
"""
n = subspace.ambient_dimension
if not embedded:
subspace = subspace.embed()
geo = HemisphereModel(n)
else:
geo = HemisphereModel(n - 1)
return central_project(subspace, geo.absolute, geo.fixed_point)
@projective_to_hemisphere.register
def projective_to_hemisphere_point(point: Point, embedded=False):
"""Convert Point from projective model to Möbius/hemisphere model.
Parameters
----------
point : Point
embedded : bool (default=False)
If False, object will be embedded into the equatorial plane before
projecting.
Returns
-------
Tuple (Point, Point)
These are the central projection images in both hemispheres. The first
object in the tuple will always be the point that lies in the
half-space {x_n+1 > 0}.
"""
as_quadric = projective_to_hemisphere_subspace(point, embedded=embedded)
p1, p2 = quadric_to_subspaces(as_quadric)
if p1.affine_point[-1] > 0:
return p1, p2
return p2, p1
@projective_to_hemisphere.register
def projective_to_hemisphere_cayley_klein_sphere(sphere: SphereLike, embedded=False):
"""Convert Cayley-Klein sphere from Projective to Möbius/hemisphere model.
Works by embedding `quadric` into the equatorial plane (the polar plane of
[e_n+1]) and projecting to the Möbius quadric via [e_n+1].
Parameters
----------
sphere : Cayley-Klein sphere
embedded : bool (default=False)
If False, object will be embedded into the equatorial plane before
projecting.
Returns
-------
Tuple (Quadric, Quadric)
These are the central projection images in both hemispheres. The first
object in the tuple will always be the quadric that lies (more) in the
half-space {x_n+1 > 0}.
"""
n = sphere.ambient_dimension
if not embedded:
moeb = HemisphereModel(n)
sphere = sphere.embed()
else:
moeb = HemisphereModel(n - 1)
lift1, lift2 = lift_sphere_to_quadric(sphere, moeb.absolute, moeb.fixed_point)
# Select the image in the upper hemisphere. The quadric is (more) in the
# upper half-space if and only if its corresponding pole is in the upper
# half-space.
if moeb.absolute.polarize(lift1.subspace).affine_point[-1] > 0:
return lift1, lift2
return lift2, lift1
[docs]def hemisphere_to_projective(object_, embedded=False):
"""Convert from Möbius/hemisphere model to projective model.
Works by centrally projecting from the Möbius quadric to the equatorial
plane (the polar plane of [e_n+1]) via [e_n+1] and computing coordinates.
Parameters
----------
object_ : Quadric or Point
embedded : bool (default=False)
Whether to return the object in the equatorial plane or in n-dim. ambient space.
returns
-------
Quadric, Point or Subspace
"""
n = object_.ambient_dimension
H = coordinate_hyperplane(n)
projection_point = HemisphereModel(n - 1).fixed_point
# Case where we can't central project because the join of projection_point
# with the quadric would be a "solid cone". This is the case where the
# quadric is a circle orthogonal to the equator, which means it represents
# a geodesic.
if isinstance(object_, Quadric) and projection_point in object_.subspace:
tc = touching_cone(projection_point, object_, in_subspace=True)
p1, p2 = quadric_to_subspaces(meet(tc, H))
img = Subspace(p1.point, p2.point).dehomogenize()
else:
img = central_project(object_, H, projection_point)
if not embedded:
img = img.unembed()
return img
[docs]def projective_to_poincare(object_, embedded=False):
"""Convert from projective model to Poincaré disk model.
Just the composition ::
Projective -> Möbius -> Poincaré
"""
in_proj = projective_to_hemisphere(object_, embedded=embedded)
if isinstance(in_proj, tuple):
# Get the image in the lower hemisphere
in_proj = in_proj[1]
return hemisphere_to_poincare(in_proj, embedded=embedded)
[docs]def poincare_to_projective(object_, embedded=False):
"""Convert from Poincaré disk model to projective model.
Just the composition ::
Poincaré -> Möbius -> Projective
"""
return hemisphere_to_projective(
poincare_to_hemisphere(object_, embedded=embedded), embedded=embedded
)
[docs]def hemisphere_to_poincare(object_, embedded=False):
"""Hemisphere to Poincaré disk model conversion.
Alias for :py:func:`.moebius_models.projective_to_euclidean`
"""
return moebius_models.projective_to_euclidean(object_, embedded=embedded)
[docs]def poincare_to_hemisphere(object_, embedded=False):
"""Poincaré disk to hemisphere model conversion.
Alias for :py:func:`.moebius_models.euclidean_to_projective`
"""
return moebius_models.euclidean_to_projective(object_, embedded=embedded)
[docs]def hemisphere_to_half_space(object_, embedded=False):
"""Convert from Möbius/hemisphere model to Poincaré Half-plane model.
This function stereographically projects from the point [1, 0,...,0, 1] on
the Möbius quadric to the hyperplane {x1 = 0}, then computes coordinates.
Warnings
--------
No choice is made whether the upper or lower half-space is used.
Make sure your object resides in the correct hemisphere.
Parameters
----------
object_ : Point or Quadric
embedded : bool (default=False)
Whether to return the object in the equatorial plane or in n-dim. ambient space.
Returns
-------
Point or SphereLike
"""
n = object_.ambient_dimension
projection_point = np.zeros(n + 1)
projection_point[[0, -1]] = 1
projection_point = Point(projection_point)
plane = subspace_from_columns(np.eye(n + 1)[:, 1:])
img = stereographic_project(object_, plane, projection_point)
if isinstance(img, Quadric):
euc = euclidean_models.ProjectiveModel(img.ambient_dimension)
img = euc.quadric_to_sphere(img)
if not embedded:
img = img.unembed(plane)
return img
[docs]def half_space_to_hemisphere(object_, embedded=False):
"""Convert from Poincaré-Half plane model to Möbius/hemisphere model.
This function optionally embeds into the plane {x1 = 0}, then
stereographically projects from that plane to the Möbius quadric via the
point [1, 0,...,0, 1].
Warnings
--------
No choice is made whether the upper or lower hemisphere is used.
Make sure your object resides in the correct half-space.
Parameters
----------
object_ : Subspace or SphereLike
embedded : bool (default=False)
If False, embeds into the the plane {x1 = 0} before projecting.
Returns
-------
Subspace or Quadric
"""
if not embedded:
n = object_.ambient_dimension
else:
n = object_.ambient_dimension - 1
projection_point = np.zeros(n + 2)
projection_point[[0, -1]] = 1
projection_point = Point(projection_point)
Q = HemisphereModel(n).absolute
plane = subspace_from_columns(np.eye(n + 2)[:, 1:])
if not embedded:
object_ = object_.embed(plane)
return inverse_stereographic_project(object_, Q, projection_point)