import numpy as np
from textwrap import indent
from ddg.abc import NonExact
import ddg.math.linalg as la
from ddg.geometry.geometries import Geometry, EuclideanGeometry, get_geometry
from ddg.math.projective import homogenize
from ddg.geometry.subspaces import (
Subspace,
Point,
whole_space,
subspace_from_affine_points_and_directions,
subspace_from_columns
)
import functools
[docs]@functools.total_ordering
class Sphere(NonExact):
"""Basic k-sphere in Euclidean space or any Cayley-Klein geometry.
Parameters
----------
center : ddg.Point or numpy.ndarray of shape (n+1,)
Center of Sphere given in homogeneous coordinates.
radius : float
Radius of Sphere.
subspace : ddg.Subspace or iterable of numpy.ndarray of shape (n+1,), optional
A (k+1)-dimensional Subspace of n-dimensional projective space
used to define k-Sphere.
If given as an iterable of numpy.ndarray, its elements will be
interpreted as the spanning points in homogeneous coordinates.
geometry : ddg.Geometry (default='euclidean')
Geometry to which Sphere belongs, can be set to any
CayleyKleinGeometry.
atol : float (default=None)
Absolute tolerance.
rtol : float (default=None)
Relative tolerance.
Attributes
----------
center
radius
subspace
ambient_dimension : int
Dimension of ambient projective space.
dimension : int
Dimension of Sphere.
Notes
-----
This class supports comparison with <, >, <=, >= and ==. These implement
the subset relation "⊆". This class also supports the "in" keyword, which
implements the element relation "∈".
"""
def __init__(self, center, radius, subspace=None, geometry='euclidean',
atol=None, rtol=None):
if not isinstance(center, Point):
center = Point(center)
self.center = center
self.radius = radius
if subspace is not None and subspace:
if not isinstance(subspace, Subspace):
points = np.array(subspace)
subspace = Subspace(*points)
if subspace.ambient_dimension != self.center.ambient_dimension:
raise ValueError(f"Dimension mismatch. Can not define Sphere in"
f" {self.center.ambient_dimension}-dimensional"
f" space with Subspace in {subspace.ambient_dimension}-"
f"dimensional space.")
else:
subspace = whole_space(self.center.ambient_dimension)
self.subspace = subspace
if self.center not in self.subspace:
raise ValueError("Center of sphere is not in containing subspace.")
self.set_geometry(geometry)
super().__init__(atol=atol, rtol=rtol)
@property
def ambient_dimension(self):
return self.center.ambient_dimension
@property
def dimension(self):
return self.subspace.dimension - 1
[docs] def cayley_klein_radius(self):
"""Return radius in terms of the Cayley-Klein distance.
Returns
-------
float
"""
return self.geometry.metric_to_cayley_klein_distance(self.radius)
[docs] def is_hypersphere(self):
"""Checks if the sphere is a hypersphere.
Returns : bool
True, if dimension of Sphere is one less than
dimension of ambient space."""
return self.dimension == self.ambient_dimension - 1
[docs] def is_circle(self):
"""Checks if the sphere is a circle.
Returns : bool
True, if Sphere is 1-dimensional Circle."""
return self.dimension == 1
[docs] def set_geometry(self, geometry):
"""Resets geometry attribute after initialization.
Parameters
----------
geometry : ddg.Geometry, str or tuple (str, int)
Geometry of the Sphere, set to 'euclidean'
at initialization.
"""
if isinstance(geometry, Geometry):
new_geo = geometry
elif isinstance(geometry, str):
new_geo = get_geometry(geometry, dimension=self.ambient_dimension)
else:
raise TypeError('Unknown identifier for geometries. Use either a '
'Geometry object or a string.')
if new_geo.dimension != self.ambient_dimension:
raise ValueError(f'Dimension mismatch: Dimension of the geometry '
f'({new_geo.dimension}) does not match ambient '
f'dimension of the subspace '
f'({self.ambient_dimension})')
self.geometry = new_geo
def __eq__(self, other):
"""
Two Spheres are equal iff they lie in same subspace in the same ambient
space, have same dimension, center point and radius.
Equality of floats is tested using asymmetric(!) np.isclose(a,b)
where b is the reference value.
"""
if not isinstance(other, Sphere):
return NotImplemented
if self.ambient_dimension != other.ambient_dimension:
return False
if self.subspace != other.subspace:
return False
atol = max(self.atol, other.atol)
rtol = max(self.rtol, other.rtol)
return (
np.isclose(self.radius, other.radius, atol=atol, rtol=rtol)
and self.center == other.center
)
def __le__(self, other):
"""Implements the subset relation "⊆".
This is currently only implemented for Euclidean spheres and works by
lifting them to Möbius geometry and checking containment of the
corresponding subspaces.
"""
if not isinstance(other, Sphere):
return NotImplemented
if not isinstance(self.geometry, EuclideanGeometry):
return NotImplemented
# I did not manage to resolve circular imports the normal way (changing
# "from a import b; b()" statements to "import a; a.b()" statements).
# register(a.b) for singledispatch functions still complained about
# circular imports.
import ddg.geometry.conversion as geoconv
# We want to embed into the standard plane, but perhaps given in a
# different basis: e_n+2 should be at position
# EuclideanGeometry.affine_component.
n = self.ambient_dimension
# First n basis vectors in ambient dimension n+1
H = np.eye(n+2, n)
i = self.geometry.affine_component
# Convert to non-negative index for np.insert
if -n <= i <= -1:
i += n + 1
H = np.insert(H, i, la.e(-1, n+2), axis=1)
H = subspace_from_columns(H)
S1 = geoconv._euclidean_to_moebius(self, subspace=H)
S2 = geoconv._euclidean_to_moebius(other, subspace=H)
# S1, S2 are quadrics. It suffices to compare subspaces (which involves
# fewer calculations), since both S1 and S2 are intersections of
# subspaces with the Möbius quadric.
return S1.subspace <= S2.subspace
def __contains__(self, p):
"""implements the element relation "∈".
A point is contained in a sphere if and only if it is contained in
`self.subspace` and the distance from `center` is (almost) equal to
`self.radius`. The distance is measured using `self.geometry.d`, if the
geometry provides it.
"""
if isinstance(p, Point):
p = p.point
if np.size(p) != self.center.ambient_dimension + 1:
return False
if p not in self.subspace:
return False
dist = self.geometry.d(self.center, p)
return np.isclose(self.radius, dist, atol=self.atol, rtol=self.rtol)
def __repr__(self):
if self.is_circle():
if self.is_hypersphere():
return (
f"Circle({repr(self.center)}, {self.radius},"
f"geometry={self.geometry})"
)
else:
return (
f"Circle({repr(self.center)}, {self.radius},"
f"subspace={repr(self.subspace)},"
f"geometry={self.geometry})"
)
elif self.is_hypersphere():
return (
f"Hypersphere({repr(self.center)}, {self.radius},"
f"geometry={self.geometry})"
)
else:
return (
f"Sphere({repr(self.center)}, {self.radius},"
f"subspace={repr(self.subspace)},"
f"geometry={self.geometry})"
)
def __str__(self):
geometry_str = self.geometry.name.split(" ", 1)[0]
if self.is_circle():
if self.is_hypersphere():
return (
f"Circle in {self.ambient_dimension}D "
f"{geometry_str.capitalize()} Geometry\n"
f"Center: {self.center.point},\n"
f"Radius: {self.radius}."
)
else:
return (
f"Circle in {self.ambient_dimension}D "
f"{geometry_str.capitalize()} Geometry\n"
f"Center: {self.center.point},\n"
f"Radius: {self.radius},\n"
f"Subspace coordinates:\n"
+ indent(str(self.subspace.matrix), ' ')
)
elif self.is_hypersphere():
return (
f"Hypersphere in {self.ambient_dimension}D "
f"{geometry_str.capitalize()} Geometry\n"
f"Center: {self.center.point},\n"
f"Radius: {self.radius}."
)
else:
return (
f"{self.dimension}-Sphere in {self.ambient_dimension}D "
f"{geometry_str.capitalize()} Geometry\n"
f"Center: {self.center.point},\n"
f"Radius: {self.radius},\n"
f"Subspace coordinates:\n"
+ indent(str(self.subspace.matrix), ' ')
)
[docs]def sphere_from_affine_point_and_normals(
point, radius, normals=[], affine_component=-1
):
"""Create an Euclidean sphere from an affine point and normal directions.
Parameters
----------
point : numpy.ndarray of shape (n,)
Affine point in n-dimensional space.
radius : float
Radius of Sphere.
normals : list of numpy.ndarray of shape (n,).
Used to define (n-k) directions as a basis for orthogonal
complement of k-dimensional span of list of normals of length k.
If given as 2-dimensional numpy.ndarray of shape (k,n), the
rows(!) of matrix will be interpreted as the normal vectors.
affine_component : int, default=-1.
Returns
-------
ddg.Sphere
An Euclidean (n-k-1)-sphere contained in a (n-k)-dimensional
projective subspace spanned by affine point and directions.
"""
if isinstance(normals, list):
normals = np.array(normals)
center = homogenize(point)
if normals.size != 0:
directions = list(la.nullspace(normals).T)
subspace = subspace_from_affine_points_and_directions(
point, directions, affine_component
)
return Sphere(center, radius, subspace=subspace, geometry="euclidean")
else:
return Sphere(center, radius, geometry="euclidean")