"""Elliptic geometry module.
This is the metric Cayley-Klein geometry with empty / positive definite absolute.
"""
from functools import singledispatch
import numpy as np
import ddg
from ddg.geometry import moebius_models
from ddg.geometry._intersection import intersect
from ddg.geometry._projections import central_project, lift_sphere_to_quadric
from ddg.geometry._quadrics import Quadric
from ddg.geometry._subspaces import Point, Subspace
from ddg.geometry.geometry_model_templates import MetricCayleyKleinGeometry
from ddg.geometry.spheres import SphereLike
__all__ = [
"ProjectiveModel",
"SphericalModel",
"spherical_to_projective",
"projective_to_spherical",
]
def _digit(num, k):
"""Find the k-th binary digit of num."""
return (num & (1 << k)) >> k
[docs]class ProjectiveModel(MetricCayleyKleinGeometry):
"""Elliptic geometry.
This is a Cayley-Klein geometry with empty absolute, so the model space is
all of RPn.
Parameters
----------
dimension : int
Attributes
----------
dimension : int
"""
@property
def absolute(self):
"""The absolute quadric with matrix ``diag([1,...,1])``.
Returns
-------
Quadric
"""
return Quadric(np.diag(np.ones(self.dimension + 1)))
[docs] @staticmethod
def metric_to_cayley_klein_distance(d):
"""Return Metric distance converted to Cayley-Klein "distance".
`K` is given by ::
K = cos(d) ** 2
Parameters
----------
d : float
Returns
-------
K : float
"""
return np.cos(d) ** 2
[docs] @staticmethod
def cayley_klein_distance_to_metric(K):
"""Return Cayley-Klein distance converted to metric distance.
`d` is given by the equation ::
d = arccos(sqrt(K))
Parameters
----------
K : float
Returns
-------
d : float
Raises
------
ValueError
If K is not in [0, 1]
"""
if ddg.nonexact.isclose(K, 0):
return np.pi / 2
if ddg.nonexact.isclose(K, 1):
return 0.0
if not 0 <= K <= 1:
raise ValueError(
"Cayley-Klein distance must be in [0, 1] to correspond to a "
"metric distance in elliptic geometry."
)
return np.arccos(np.sqrt(K))
[docs] def moebius(self):
"""Corresponding projective model of Moebius geometry (Spherical model)."""
return SphericalModel(self.dimension)
[docs] def to_moebius(self, object_, embedded=False):
"""Alias for :py:func:`.projective_to_spherical`"""
return projective_to_spherical(object_, embedded=embedded)
[docs] def from_moebius(self, object_, embedded=False):
"""Alias for :py:func:`.spherical_to_projective`"""
return spherical_to_projective(object_, embedded=embedded)
[docs] def sphere_through_points(self, *points):
"""Find the k-spheres through given points.
Returns all the spheres of minimal dimension that contain the given points.
There might be up to `2 ** (len(points) - 1)` of these which are returned
as a list ordered by increasing radius.
The spheres might degenerate to subspaces which are
considered to have infinite radius.
Generically, n points define an (n-2)-sphere.
Parameters
----------
*points : iterable of type `ddg.geometry.Point`
The points for which to find the sphere.
Returns
-------
list of ddg.geometry.spheres.MetricCaleyKleinSphere or ddg.geometry.Subspace
The length of this list varies from 1 to `2 ** (len(points) - 1)`.
Raises
------
ValueError
- If the points are not contained in a sphere.
- If a given point is not an element of the according geometry.
See also
--------
ddg.geometry.euclidean_models.ProjectiveModel.sphere_through_points,
ddg.geometry.hyperbolic_models.ProjectiveModel.sphere_through_points
"""
sph = SphericalModel(self.dimension)
# Lift the points to the spherical model
lifted_points = []
for point in points:
if point in self:
lifted_points.append(projective_to_spherical(point))
else:
raise ValueError(f"A given point is not an element of {self}.")
spheres_and_radii = []
for i in range(2 ** (len(lifted_points) - 1)):
subspace = ddg.geometry.join(
*[lifted_points[j][_digit(i, j)] for j in range(len(lifted_points))]
)
if subspace.dimension > self.dimension:
continue
quadric = ddg.geometry.intersect(subspace, sph.absolute)
sphere = spherical_to_projective(quadric)
if sphere not in [sph for sph, _ in spheres_and_radii]:
if isinstance(sphere, ddg.geometry.Subspace):
spheres_and_radii.append((sphere, np.inf))
else:
spheres_and_radii.append((sphere, sphere.metric_radius()))
if spheres_and_radii == []:
raise ValueError("The given points are not contained in a sphere.")
spheres = [sph for sph, _ in sorted(spheres_and_radii, key=lambda x: x[1])]
return spheres
[docs] def orthogonal_sphere(self, *spheres):
"""Find the set of minimal orthogonal spheres for the given `*spheres`.
If the output of this function is not an empty list, then any sphere
of positive radius in the output list possesses the following property:
Any higher-dimensional sphere containing the former sphere is orthogonal
to all spheres given in the input.
Thus, it is an orthogonal sphere of minimal dimension, if it is unique.
This function returns all such spheres of minimal dimension. The list is
sorted by increasing radius. Subspaces can also be inserted and returned
as they can be thought of as spheres with maximal radius. Returns the
empty list if the set is empty (including the case when there
is no orthogonal spheres).
Works by lifting `*spheres` to SphericalModel. Since each lift
is a pair of antipodal spheres, there are 2**len(spheres) possible lifts
(the complementary lifts are antipodal and then their ort spheres
are projected back to the same ones, so the maximal possible
output consist of 2**(len(spheres) - 1) spheres),
each of which could provide a minimal orthogonal sphere or not.
A minimal orthogonal sphere for a lift is defined as
the intersection of all spheres orthogonal to all lifted ones.
This definition matches
`ddg.geometry.euclidean_models.ProjectiveModel.orthogonal_sphere`.
For example, m non-intersecting hyperspheres
in general position in n-dimensional elliptic space
have 2**(m - 1) orthogonal (m-2)-spheres,
which are returned by this function.
Parameters
----------
*spheres: iterable of type ddg.geometry.spheres.MetricCayleyKleinSphere or
ddg.geometry.Subspace
Cayley Klein spheres or subspaces
Returns
-------
list of ddg.geometry.spheres.MetricCayleyKleinSphere or ddg.geometry.Subspace
List of spheres orthogonal to the given spheres
sorted by increasing radius.
The length of the list varies from 0 to 2**(len(spheres) - 1).
See also
--------
ddg.geometry.euclidean_models.ProjectiveModel.orthogonal_sphere,
ddg.geometry.hyperbolic_models.ProjectiveModel.orthogonal_sphere
"""
sph = SphericalModel(self.dimension)
# Lift the spheres to spherical model.
# We have two types of lifts: tuples (lift of a sphere)
# and unique spheres (lift of a subspace)
sph_tuples = []
sph_spheres = []
for s in spheres:
l_s = projective_to_spherical(s)
if isinstance(l_s, tuple):
sph_tuples.append(l_s)
else:
sph_spheres.append(l_s)
ort_spheres_and_radii = []
# For each possible lift find the orthogonal sphere.
# First find the intersection of non-tuple subspaces
if len(sph_spheres) == 0:
core = sph.absolute.subspace
else:
core = intersect(*[s.subspace for s in sph_spheres])
for i in range(2 ** (len(sph_tuples) - 1)):
current_polar = core
for j in range(len(sph_tuples)):
# We use j-th binary digit of i to choose a lift.
# Note that the (len(sph_tuples) - 1)-th digit is always 0
current_polar = intersect(
current_polar, sph_tuples[j][_digit(i, j)].subspace
)
if current_polar.dimension < 0:
continue
current_ort = intersect(sph.absolute, sph.absolute.polarize(current_polar))
if current_ort.dimension < 0:
continue
# Project it back
pr_current_ort = spherical_to_projective(current_ort)
if pr_current_ort not in [sph for sph, _ in ort_spheres_and_radii]:
if isinstance(pr_current_ort, Subspace):
# Note that the "radius" of a subspace is considered
# to be infinite.
ort_spheres_and_radii.append((pr_current_ort, np.inf))
else:
ort_spheres_and_radii.append(
(pr_current_ort, pr_current_ort.metric_radius())
)
ort_spheres = [
sph for sph, _ in sorted(ort_spheres_and_radii, key=lambda x: x[1])
]
return ort_spheres
def __contains__(self, point):
return point.ambient_dimension == self.ambient_dimension
def __str__(self):
return f"Projective model of {self.dimension}D elliptic geometry"
[docs]class SphericalModel(moebius_models._ProjectiveSubgeometry):
"""Spherical model of elliptic geometry.
Also known as the Möbius model because this is a subgeometry of Möbius
geometry.
.. rubric:: Model space
The model space is the unit sphere where antipodal points are identified. More
precisely::
S^n/{±1} = { {[x], [y]} ⊆ S^n : dehomogenize([x]) = − dehomogenize([y]) }
where `S^n = Quadric(diag([1,...,1, -1]))` is the Möbius quadric.
.. 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.
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 inside the Möbius quadric used for computation of the metric.
This is the point with representative [0,...,0, 1].
Returns
-------
ddg.geometry.Point
"""
return self.elliptic_point
[docs] def d(self, v, w):
"""Elliptic metric in spherical model
Given by the formula ::
| <x, y> + <x, p> <y, p> |
| ---------------------- | = cos(d([x],[y])
| <x, p> <y, p> |
Where `p` is `self.fixed_point.point`.
Parameters
----------
v, w : ddg.geometry.Point or array_like of shape (n,)
Returns
-------
float
"""
if isinstance(v, ddg.geometry.Point):
v = v.point
if isinstance(w, ddg.geometry.Point):
w = w.point
p = self.fixed_point.point
b = self.absolute.inner_product
c = np.abs((b(v, w) + b(v, p) * b(w, p)) / (b(v, p) * b(w, p)))
if ddg.nonexact.isclose(c, 1):
# Due to numerical instabilities, the value c may be slightly
# larger than 1. Then it cannot be inserted into arccos.
return 0.0
return np.arccos(c)
def __str__(self):
return (
f"Spherical model of {self.dimension}D elliptic geometry\n"
+ f"in the projective model of {self.dimension}D Moebius geometry"
)
[docs]def spherical_to_projective(object_, embedded=False):
"""Convert from spherical model to projective model.
Works by centrally projecting from the Möbius quadric to the hyperplane at
infinity and computing coordinates.
If a quadric is given which is contained in the Möbius quadric, the result
will be converted to a corresponding Caley-Klein sphere.
Parameters
----------
object_ : ddg.geometry.Quadric or ddg.geometry.Point
embedded : bool (default=False)
Whether to return the object in the hyperplane at infinity
or in n-dim. ambient space.
returns
-------
ddg.geometry.Quadric,
ddg.geometry.Point,
ddg.geometry.Subspace or
ddg.geometry.spheres.MetricCayleyKleinSphere
"""
n = object_.ambient_dimension - 1
sph = SphericalModel(n)
projection_point = sph.fixed_point
H = ddg.geometry.Subspace(*np.eye(n + 1, n + 2))
if isinstance(object_, Quadric) and projection_point in object_.subspace:
img = intersect(object_.subspace, H)
elif isinstance(object_, ddg.geometry.Quadric) and object_ <= sph.absolute:
ell = ProjectiveModel(n)
polar_subspace = sph.absolute.polarize(object_.subspace)
projected_subspace = central_project(object_.subspace, H, projection_point)
projected_polar_subspace = central_project(polar_subspace, H, projection_point)
center = ddg.geometry.intersect(projected_subspace, projected_polar_subspace)
center_coords = center.unembed(H).point
point_on_quadric = ddg.geometry.Point(
ddg.math.inner_product.light_like_vectors_from_onb(
ddg.math.inner_product.orthonormalize(
np.asarray(object_.subspace.points).T, object_.inner_product
).T,
object_.inner_product,
)[0, :]
)
point_on_sphere_coords = (
central_project(point_on_quadric, H, projection_point).unembed(H).point
)
cayley_klein_radius = ell.absolute.inner_product(
center_coords, point_on_sphere_coords
) ** 2 / (
ell.absolute.inner_product(center_coords, center_coords)
* ell.absolute.inner_product(point_on_sphere_coords, point_on_sphere_coords)
)
img_unembedded = ell.sphere(
center_coords,
ell.cayley_klein_distance_to_metric(cayley_klein_radius),
projected_subspace.unembed(H),
)
img = img_unembedded.embed(H)
else:
img = central_project(object_, H, projection_point)
if not embedded:
img = img.unembed(H)
return img
[docs]@singledispatch
def projective_to_spherical(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 spherical model is not supported."
)
@projective_to_spherical.register
def projective_to_spherical_subspace(subspace: Subspace, embedded=False):
"""Convert Subspace from Projective model to spherical model.
Works by optionally embedding `subspace` into the hyperplane at infinity (the
polar plane of [e_n+2]), then projecting to the Möbius quadric via [e_n+2].
Parameters
----------
subspace : ddg.geometry.Subspace
embedded : bool (default=False)
If False, object will be embedded into the hyperplane at infinity before
projecting.
Returns
-------
ddg.geometry.Quadric
"""
n = subspace.ambient_dimension
if not embedded:
sph = SphericalModel(n)
H = ddg.geometry.Subspace(*np.eye(n + 1, n + 2))
subspace = subspace.embed(H)
else:
sph = SphericalModel(n - 1)
return central_project(subspace, sph.absolute, sph.fixed_point)
@projective_to_spherical.register
def projective_to_spherical_point(point: Point, embedded=False):
"""Convert Point from projective model to spherical model.
Works by optionally embedding `subspace` into the hyperplane at infinity (the
polar plane of [e_n+2]), then projecting to the Möbius quadric via [e_n+2].
There will always be two antipodal points resulting from this projection.
Both will be returned in no particular order.
Parameters
----------
point : ddg.geometry.Point
embedded : bool (default=False)
If False, object will be embedded into the hyperplane at infinity before
projecting.
Returns
-------
Tuple (ddg.geometry.Point, ddg.geometry.Point)
These are the two antipodal central projection images. They are returned
in no particular order.
"""
as_quadric = projective_to_spherical_subspace(point, embedded=embedded)
p1, p2 = ddg.geometry.quadric_to_subspaces(as_quadric)
return p1, p2
@projective_to_spherical.register
def projective_to_spherical_cayley_klein_sphere(sphere: SphereLike, embedded=False):
"""Convert Cayley-Klein sphere from projective model to spherical model.
Works by optionally embedding `quadric` into the hyperplane at infinity (the
polar plane of [e_n+2]), then projecting to the Möbius quadric via [e_n+2].
There will always be two antipodal quadrics resulting from this projection.
Both will be returned in no particular order.
Parameters
----------
sphere : Cayley-Klein sphere
embedded : bool (default=False)
If False, object will be embedded into the hyperplane at infinity before
projecting.
Returns
-------
Tuple (ddg.geometry.Quadric, ddg.geometry.Quadric)
These are the two antipoldal central projection images. They are returned
in no particular order.
"""
n = sphere.ambient_dimension
if not embedded:
sph = SphericalModel(n)
H = ddg.geometry.Subspace(*np.eye(n + 1, n + 2))
sphere = sphere.embed(H)
else:
sph = SphericalModel(n - 1)
return lift_sphere_to_quadric(sphere, sph.absolute, sph.fixed_point)