from textwrap import indent
import numpy as np
from scipy.linalg import orth
import ddg.math.inner_product as ip
import ddg.math.linalg as la
import ddg.math.projective as pmath
from ddg import nonexact
from ddg.abc import LinearTransformable
from ddg.geometry.abc import Embeddable
from ddg.math.euclidean import embed
from ddg.math.random import random_matrix
[docs]class Subspace(Embeddable, LinearTransformable):
"""Subspace of a projective space.
Represented by a number of vectors in homogeneous coordinates.
Parameters
----------
*points : array_like of list or numpy.ndarray of shape (n,)
Homogeneous coordinate vectors spanning the subspace. If the span is
1-dimensional, automatically casts to _subspaces.Point
Non-float dtypes will be converted to float.
atol, rtol : float (default=None)
Tolerances that will be used internally by the class. If these are set to
None, the global defaults will be used.
Attributes
----------
points
matrix
affine_points
affine_matrix
affine_point_and_directions
dimension
ambient_dimension
codimension
Methods
-------
at_infinity
See Also
--------
ddg.geometry.Point
Notes
-----
Subspaces can be empty. In this case, points are stored as an empty matrix
with shape (n,0) if the ambient dimension is n-1. Empty spaces must be
initialized with 0-vectors.
"""
def __init__(self, *points, atol=None, rtol=None):
self.atol = atol
self.rtol = rtol
# Exception blocks
if len(np.shape(points[0])) != 1:
raise ValueError(
"Vectors spanning the subspace must be 1D "
"arrays given as separate arguments."
)
if np.size(points[0]) == 0:
raise ValueError(
"Can't determine ambient dimension from empty coordinate vectors."
)
for p in points[1:]:
if np.shape(p) != np.shape(points[0]):
raise ValueError("Point vector length mismatch.")
point_matrix = np.array(points).T
if point_matrix.dtype == complex:
raise TypeError("Complex subspaces are not currently supported.")
# Convert integer matrices to float matrices to avoid confusing behavior
point_matrix = point_matrix.astype(float)
# Reduce linearly dependent vectors.
# Also converts np.zeros(n) to an empty array of shape (n,0).
temp = la.col_basis(point_matrix, atol=self.atol, rtol=self.rtol)
if np.shape(temp) != np.shape(point_matrix):
point_matrix = temp
self.matrix = point_matrix
# Cast to Point. This would bypass Point.__init__ if it existed.
if np.shape(self.matrix)[1] == 1:
self.__class__ = Point
@property
def points(self):
"""List of points spanning the subspace in homogeneous coordinates.
Returns
-------
list of numpy.ndarray of shape (ambient_dimension+1,)
"""
return la.matrix_to_arrays(self.matrix)
[docs] def embed(self, subspace=None):
if subspace is None:
subspace = coordinate_hyperplane(self.ambient_dimension + 1)
return subspace_from_columns(subspace.matrix @ self.matrix)
[docs] def unembed(self, subspace=None):
if subspace is None:
subspace = coordinate_hyperplane(self.ambient_dimension)
if not self <= subspace:
raise ValueError("subspace must contain self.")
return subspace_from_columns(la.coordinates(self.matrix, subspace.matrix))
@property
def affine_points(self):
"""List of points spanning the affine part of the subspace.
Computes points in affine coordinates that, when homogenized, form a
basis for the subspace.
Returns
-------
list of numpy.ndarray of shape (ambient_dimension,)
"""
return la.matrix_to_arrays(self.affine_matrix)
@property
def affine_matrix(self):
"""Matrix whose columns span the affine part of the subspace.
Computes points in affine coordinates that, when homogenized, form a
basis for the subspace. Returns them as the columns of a matrix.
Returns
-------
numpy.ndarray of shape (ambient_dimension, dimension+1)
Raises
------
ValueError
If the subspace is at infinity.
"""
if self.at_infinity():
raise ValueError("Subspace at infinity.")
points = self.matrix.copy()
basemask = nonexact.isclose(points[-1], 0, atol=self.atol, rtol=self.rtol)
pointmask = np.where(np.logical_not(basemask))[0]
directmask = np.where(basemask)[0]
points[:, pointmask] = points[:, pointmask] / points[-1, pointmask]
point = points.T[pointmask[0]]
point = point[:, np.newaxis]
norms = np.linalg.norm(points[:, directmask], axis=0)
points[:, directmask] = points[:, directmask] / norms
points[:, directmask] = points[:, directmask] + point
return np.delete(points, -1, axis=0)
@property
def affine_point_and_directions(self):
"""Return point in the subspace and directions.
Returns affine coordinate vectors p and d1, ..., dk such that
the original subspace is spanned by [p,1], [d1,0], ..., [dk,0].
The point is the one that is closest to the zero vector in affine
coordinates and the directions are orthonormal.
Returns
-------
numpy.ndarray, list of numpy.ndarray
See also
--------
ddg.geometry.subspace_from_affine_points_and_directions
"""
if self.at_infinity():
raise ValueError("Subspace at infinity.")
points = self.matrix.copy()
basemask = nonexact.isclose(points[-1], 0, atol=self.atol, rtol=self.rtol)
# dehomogenize the points that are not at infinity
pointmask = np.where(np.logical_not(basemask))[0]
points[:, pointmask] = points[:, pointmask] / points[-1, pointmask]
# isolate the point to be returned
pidx = pointmask[0]
pointmask = pointmask[1:]
point = points.T[pidx]
point = point[:, np.newaxis]
# Directions as differences. The other directions are points at
# infinity and we can just delete the 0.
points[:, pointmask] -= point
if pidx != 0:
points[:, (0, pidx)] = points[:, (pidx, 0)]
temp = np.delete(points, -1, axis=0)
directions = temp[:, 1:]
if np.size(directions) != 0:
directions = orth(directions)
directions = list(directions.T)
point = temp[:, 0]
# project so the point is the closest one to the origin
point -= sum(np.dot(point, d) * d for d in directions)
return point, directions
[docs] def at_infinity(self):
points = self.matrix
tmp = nonexact.isclose(points[-1], 0, atol=self.atol, rtol=self.rtol)
return tmp.all()
@property
def dimension(self):
return np.shape(self.matrix)[1] - 1
@property
def ambient_dimension(self):
return np.shape(self.matrix)[0] - 1
@property
def codimension(self):
return self.ambient_dimension - self.dimension
[docs] def dual(self, ambient_space=None):
"""Return dual subspace.
Parameters
----------
ambient_space : ddg.geometry.Subspace (default=None)
Other subspace containing this subspace that should be used as
ambient space when dualizing. None means whole space.
Returns
-------
ddg.geometry.Subspace
"""
return _dualize(ip.col_complement, self, ambient_space)
[docs] def dual_oriented(self, ambient_space=None):
"""Return oriented dual subspace.
Parameters
----------
ambient_space : ddg.geometry.Subspace (default=None)
Other subspace containing this subspace that should be used as
ambient space when dualizing. None means whole space.
Returns
-------
ddg.geometry.Subspace
"""
return _dualize(ip.col_complement_oriented, self, ambient_space)
[docs] def orthonormalize(self):
"""Returns subspace with a regularized basis. Useful for visualization.
See :py:func:`.orthonormalize_subspace` for full documentation.
See also
--------
ddg.geometry.orthonormalize_subspace
Equivalent function
"""
return orthonormalize_subspace(self)
[docs] def center(self, center, remove_index=None, insert_index=0):
"""Returns subspace with a point inserted into the basis. Useful for
visualization.
See :py:func:`.center_subspace` for full documentation.
See also
--------
ddg.geometry.center_subspace
Equivalent function
"""
return center_subspace(
self, center, remove_index=remove_index, insert_index=insert_index
)
[docs] def orthonormalize_and_center(self, center):
"""Returns subspace with a regularized basis and a point
inserted in the basis
Useful for visualization.
See :py:func:`.orthonormalize_and_center_subspace` for full documentation.
See also
--------
ddg.geometry.orthonormalize_and_center_subspace
Equivalent function
"""
return orthonormalize_and_center_subspace(self, center)
[docs] def dehomogenize(self):
"""Returns subspace with last entries of basis vectors normalized to 1.
Useful for visualization.
See :py:func:`.dehomogenize_subspace` for full documentation.
See also
--------
ddg.geometry.dehomogenize_subspace
Equivalent function
"""
return dehomogenize_subspace(self)
def __contains__(self, p):
"""
Checks whether a homogeneous coordinate vector or Point object is
contained in the subspace.
"""
if isinstance(p, Point):
p = p.point
if len(p) != self.ambient_dimension + 1:
raise ValueError(
"Dimension mismatch: Neither a homogeneous "
"coordinate vector nor Point object in "
f"{self.ambient_dimension}D space."
)
eps = [np.dot(p, q) for q in self.dual().points]
return all(
nonexact.isclose(e, 0.0, atol=self.atol, rtol=self.rtol) for e in eps
)
def __ge__(self, other):
"""
Checks if subspace contains `other`.
"""
if isinstance(other, Subspace):
if self.ambient_dimension != other.ambient_dimension:
raise ValueError(
"Ambient dimension mismatch: "
f"{self.ambient_dimension} and "
f"{other.ambient_dimension}."
)
points = other.matrix
dual_matrix = self.dual().matrix.T
eps = dual_matrix.dot(points)
return nonexact.isclose(
eps, np.zeros(eps.shape), atol=self.atol, rtol=self.rtol
).all()
else:
raise TypeError(repr(other) + " is not of ddg.geometry.Subspace class.")
def __le__(self, other):
"""
Checks if subspace is contained in `other`.
"""
return other.__ge__(self)
def __eq__(self, other):
if not isinstance(other, Subspace):
return False
return self <= other and self >= other
def __repr__(self):
if self.dimension != -1:
points_param = ", ".join([repr(p) for p in self.points])
else:
points_param = repr(np.zeros(self.ambient_dimension + 1))
return f"Subspace({points_param})"
def __str__(self):
if self.dimension == -1:
return f"Empty space in {self.ambient_dimension}D projective space"
elif self.codimension == 0:
s = "Whole space"
elif self.dimension == 1:
s = "Line"
elif self.dimension == 2:
s = "Plane"
elif self.codimension == 1:
s = "Hyperplane"
else:
s = f"{self.dimension}D Subspace"
return (
f"{s} in {self.ambient_dimension}D projective space\n"
"Homogeneous coordinates (columns):\n" + indent(str(self.matrix), " ")
)
[docs]class Point(Subspace):
"""Subclass for points (0-dim. projective subspaces).
Parameters
----------
point : list or numpy.ndarray of shape (n,)
atol, rtol : float (default=None)
Tolerances that will be used internally by the class. If these are set to
None, the global defaults will be used.
Attributes
----------
point
affine_point
See also
--------
ddg.geometry.Subspace
"""
def __init__(self, point, atol=None, rtol=None):
super().__init__(point, atol=atol, rtol=rtol)
@property
def point(self):
"""Homogeneous coordinate vector representing the point.
Returns
-------
numpy.ndarray (dimension+1,)
"""
return self.points[0]
@property
def affine_point(self):
"""Affine coordinate vector representing the point.
Returns
-------
numpy.ndarray of shape (dimension,)
Raises
------
ValueError
If the point is at infinity.
"""
return self.affine_points[0]
def __repr__(self):
return f"Point({repr(self.point)})"
def __str__(self):
return (
f"Point in {self.ambient_dimension}D projective space\n"
f"Homogeneous coordinates: {self.point}"
)
[docs]def subspace_from_rows(matrix, **kwargs):
"""Create subspace using the rows of a matrix as points.
Parameters
----------
matrix : numpy.ndarray of shape (k+1,n+1)
**kwargs : dict
Other arguments to be passed to the Subspace.
Returns
-------
ddg.geometry.Subspace
Notes
-----
An empty subspace can be created by giving a 0-matrix or a matrix with
shape (0,n+1).
"""
if np.size(matrix) == 0:
matrix = [np.zeros(np.shape(matrix)[1])]
return Subspace(*matrix, **kwargs)
[docs]def subspace_from_columns(matrix, **kwargs):
"""Create subspace using the columns of a matrix as points.
Parameters
----------
matrix : numpy.ndarray of shape (n+1,k+1)
**kwargs : dict
Other arguments to be passed to the Subspace.
Returns
-------
ddg.geometry.Subspace
Notes
-----
An empty subspace can be created by giving a 0-matrix or a matrix with
shape (n,0).
"""
return subspace_from_rows(matrix.T, **kwargs)
[docs]def hyperplane_from_normal(normal, level=None, point=None, **kwargs):
"""
An affine hyperplane defined by either
<x, normal> == level,
or the normal and a point (given in affine coordinates) on the plane.
"""
if (level is None) == (point is None):
raise ValueError(
"Hyperplane must be specified by a normal and either "
"the level or one point."
)
elif level is not None:
h = Point(np.concatenate([normal, [-level]]), **kwargs).dual_oriented()
else:
level = np.dot(normal, point)
h = hyperplane_from_normal(normal, level=level, **kwargs)
if "atol" in kwargs:
h.atol = kwargs["atol"]
if "rtol" in kwargs:
h.rtol = kwargs["rtol"]
return h
[docs]def normals(subspace):
"""Return all normal directions of a subspace.
If p1, p2 are two points in the subspace in affine coordinates and n is a
normal direction, then
<p1 - p2, n> = 0.
Parameters
----------
subspace : ddg.geometry.Subspace
Returns
-------
list of numpy.ndarray of shape (subspace.ambient_dimension,)
Normals in affine coordinates.
"""
if subspace.dimension < 1:
return whole_space(subspace.ambient_dimension)
dirs = subspace.affine_point_and_directions[1]
dirs = np.column_stack(dirs)
normals = ip.col_complement_oriented(dirs)
return list(normals.T)
[docs]def normal_with_level(hyperplane):
"""Unit normal vector and negative level of a hyperplane in one vector.
A hyperplane having normal `N` and level `l` means that in affine coordinates, the
hyperplane is the set of points `x` such that ::
<x, N> = l,
which is equivalent to ::
<x, N> - l = 0.
This function computes a representative vector `v = [N, -l]` of the dual space,
normalized such that `N` has unit length.
`-N` and `-l` define the same hyperplane. To fix this sign, an orientation would
have to be given on the hyperplane. This function does not distinguish between the
two choices, so the sign should be treated as random.
Parameters
----------
hyperplane : ddg.geometry.Subspace
Hyperplane in n-dimensional ambient space.
Returns
-------
Nl : numpy.ndarray of shape (n+1,)
The vector `[N, -l]`. The subspace `Point(Nl).dual()` is the original
hyperplane.
"""
if hyperplane.codimension != 1:
raise TypeError(
"This function only supports hyperplanes. For multiple normal directions, "
"use the function ddg.geometry._subspaces.normals."
)
if hyperplane.at_infinity():
raise ValueError("Hyperplane is at infinity.")
n = hyperplane.dual_oriented().point
return n / np.linalg.norm(n[:-1])
[docs]def normal(hyperplane):
"""Unit normal vector N of a hyperplane.
The hyperplane consists of all points `x` such that ::
<x, N> = l
where `l` is the level of the hyperplane.
Parameters
----------
hyperplane : ddg.geometry.Subspace
Hyperplane in n-dimensional ambient space.
Returns
-------
N : numpy.ndarray of shape (n,)
See also
--------
ddg.geometry.normal_with_level
ddg.geometry.level
"""
return normal_with_level(hyperplane)[:-1]
[docs]def level(hyperplane):
"""Level l of a hyperplane.
The hyperplane consists of all points `x` such that ::
<x, N> = l
where `N` is the unit normal of the hyperplane.
The level is the distance of the hyperplane from the origin in the normal direction,
i.e. the hyperplane contains l * N.
Parameters
----------
hyperplane : ddg.geometry.Subspace
Hyperplane in n-dimensional ambient space.
Returns
-------
l : float
See also
--------
ddg.geometry.normal_with_level
ddg.geometry.normal
"""
return -normal_with_level(hyperplane)[-1]
[docs]def subspace_from_affine_points(*points, **kwargs):
"""Create a subspace from affine points.
Parameters
----------
*points : iterable of numpy.ndarray or list
Affine coordinate vectors
**kwargs : dict
Other keyword arguments to be passed to the subspace.
Returns
-------
ddg.geometry.Subspace
Subspace spanned by the given points
"""
points = [pmath.homogenize(p) for p in points]
return Subspace(*points, **kwargs)
[docs]def subspace_from_affine_rows(matrix, **kwargs):
"""Create a subspace from the rows of a matrix, interpreted as affine
coordinates.
Parameters
----------
matrix : numpy.ndarray
Array whichs columns describe the affine points
**kwargs : dict
Other keyword arguments to be passed to the subspace.
Returns
-------
ddg.geometry.Subspace
Subspace spanned by the given points
"""
return subspace_from_affine_points(*matrix, **kwargs)
[docs]def subspace_from_affine_columns(matrix, **kwargs):
"""Create a subspace from the columns of a matrix, interpreted as affine
coordinates.
Parameters
----------
matrix : numpy.ndarray
Array whichs columns describe the affine points
**kwargs : dict
Other keyword arguments to be passed to the subspace.
Returns
-------
ddg.geometry.Subspace
Subspace spanned by the given points
"""
return subspace_from_affine_rows(matrix.T, **kwargs)
[docs]def subspace_from_affine_points_and_directions(points, directions, **kwargs):
"""Create a subspace from affine points and directions.
Parameters
----------
points : (list of numpy.array) or numpy.array
List of points or a single point
directions : (list of numpy.array) or numpy.array
List of directions or a single direction
**kwargs : dict
Other keyword arguments to be passed to the subspace.
Returns
-------
ddg.geometry.Subspace
Subspace spanned by the given points and directions
Raises
------
ValueError
If both points and directions are empty lists.
"""
if len(np.shape(points)) == 1 and len(points) != 0:
points = [points]
if len(np.shape(directions)) == 1 and len(directions) != 0:
directions = [directions]
num_pts = len(points)
pts_and_dirs = points + directions
if len(pts_and_dirs) == 0:
raise ValueError("Points and directions can not both be empty.")
for i, p in enumerate(pts_and_dirs[:num_pts]):
pts_and_dirs[i] = pmath.homogenize(p)
for i, d in enumerate(pts_and_dirs[num_pts:]):
pts_and_dirs[num_pts + i] = embed(d)
return Subspace(*pts_and_dirs, **kwargs)
[docs]def random_subspace(ambient_dimension, subspace_dimension=None, seed=None):
"""Create a random subspace.
Parameters
----------
ambient_dimension : int
subspace_dimension : int (default=None)
Dimension of the resulting subspace. If None, this is chosen randomly between 0
and *ambient_dimension* (inclusive).
seed : int, numpy.random.Generator or None (default=None)
Passes a fixed seed to default Generator class. If passed a Generator,
just uses that generator for random generation.
Returns
-------
ddg.geometry.Subspace
"""
gen = np.random.default_rng(seed=seed)
if subspace_dimension is None:
subspace_dimension = gen.integers(0, ambient_dimension + 1)
return subspace_from_columns(
random_matrix(ambient_dimension + 1, subspace_dimension + 1, seed=gen)
)
[docs]def random_contained_subspace(subspace, dimension=None, seed=None):
"""Create a random subspace contained in another subspace.
Parameters
----------
subspace : ddg.geometry.Subspace
The subspace that should contain the resulting subspace.
dimension : int (default=None)
Dimension of the resulting subspace. If None, this is chosen randomly between 0
and `subspace.dimension` (inclusive).
seed : int, numpy.random.Generator or None (default=None)
Passes a fixed seed to default Generator class. If passed a Generator,
just uses that generator for random generation.
Returns
-------
ddg.geometry.Subspace
Raises
------
ValueError
If dimension is greater than `subspace.dimension`.
"""
gen = np.random.default_rng(seed=seed)
if dimension is None:
dimension = gen.integers(0, subspace.dimension + 1)
elif dimension > subspace.dimension:
raise ValueError("dimension can not be greater than subspace.dimension.")
coords = random_matrix(subspace.dimension + 1, dimension + 1, seed=gen)
return subspace_from_columns(subspace.matrix @ coords)
[docs]def orthonormalize_subspace(subspace):
"""Regularize the basis of a subspace. Useful for visualization.
Returns an equal subspace with a basis of the form `[p, v1,...,vk]`, where
`p` is the point in the subspace closest to the origin with last entry 1
and `v1,...,vk` are orthonormal direction vectors at infinity. This is only
possible if the subspace is not at infinity.
When using the "convex" parametrization in to_smooth_net, the
parametrization of the subspace with this new basis will be ::
f(x1,...,xk) = p + x1 * v1 + ... + xk * vk
Parameters
----------
subspace : ddg.geometry.Subspace
Returns
-------
ddg.geometry.Subspace
The same subspace with the regularized basis.
Raises
------
ValueError
If subspace is at infinity.
"""
if subspace.at_infinity():
raise ValueError(
'Can not compute regularized "orthonormal" basis of subspace at infinity.'
)
return subspace_from_affine_points_and_directions(
*subspace.affine_point_and_directions, atol=subspace.atol, rtol=subspace.rtol
)
[docs]def center_subspace(subspace, center, remove_index=None, insert_index=0):
"""Insert a specific point into the basis of a subspace.
This is especially useful in preparation for visualization: When applying
this function and then using the "convex" parametrization in
`to_smooth_net`, the parametrization will be centered around the new
`center`, if it is a point not at infinity.
Parameters
----------
subspace : ddg.geometry.Subspace
center : ddg.geometry.Point or numpy.ndarray of shape (k,) or (k+1,)
Where `k` is the ambient dimension of the subspace. Point must lie in
subspace. Can be given as Point object or as array containing affine or
homogeneous coordinates.
remove_index : int (default=None)
Index of basis vector that should be replaced with `center`. If None,
we search for one.
insert_index : int (default=0)
Index where center should be inserted.
Returns
-------
ddg.geometry.Subspace
The same subspace with a basis that has `center` at position
`insert_index`.
Raises
------
ValueError
* If center does not lie in the subspace.
* If the vector at `remove_index` could not be replaced with `center`
while still having a basis.
"""
if isinstance(center, Point):
center = center.point
elif np.size(center) == subspace.ambient_dimension:
center = pmath.homogenize(center)
if center not in subspace:
raise ValueError("Center must lie in subspace.")
# Let v1,...,vk be a basis. If w = a1 * v1 + ... + ak * vk and ai != 0,
# then v1,...,w,...,vk is again a basis, where vi was exchanged with w.
coords = la.coordinates(
center, subspace.matrix, atol=subspace.atol, rtol=subspace.rtol
)
if remove_index is not None:
idx = remove_index
if nonexact.isclose(coords[idx], 0, atol=subspace.atol, rtol=subspace.rtol):
raise ValueError(
f"Could not replace vector {subspace.points[remove_index]} "
f"with center {center} and still have a basis."
)
else:
for idx, coord in enumerate(coords):
if not nonexact.isclose(coord, 0, atol=subspace.atol, rtol=subspace.rtol):
break
mask = np.full(subspace.dimension + 1, True)
mask[idx] = False
# Convert to non-negative index, because otherwise insert will act
# strangely.
if -(subspace.dimension + 1) <= insert_index <= -1:
insert_index += subspace.dimension + 1
return subspace_from_columns(
np.insert(subspace.matrix[:, mask], insert_index, center, axis=1),
atol=subspace.atol,
rtol=subspace.rtol,
)
[docs]def orthonormalize_and_center_subspace(subspace, center):
"""Regularize the basis of a subspace. Useful for visualization.
Returns an equal subspace with a basis of the form `[p, v1,...,vk]`, where
`p` is the center with last entry 1
and `v1,...,vk` are orthonormal direction vectors at infinity. This is only
possible if the subspace is not at infinity.
When using the "convex" parametrization in to_smooth_net, the
parametrization of the subspace with this new basis will be ::
f(x1,...,xk) = p + x1 * v1 + ... + xk * vk
Parameters
----------
subspace : ddg.geometry.Subspace
center : ddg.geometry.Point or numpy.ndarray of shape (k,) or (k+1,)
Where `k` is the ambient dimension of the subspace. Point must lie in
subspace. Can be given as Point object or as array containing affine or
homogeneous coordinates.
Returns
-------
ddg.geometry.Subspace
The same subspace with the regularized basis.
Raises
------
ValueError
* If center does not lie in the subspace.
* If subspace is at infinity.
"""
if isinstance(center, Point):
center = center.point
elif np.size(center) == subspace.ambient_dimension:
center = pmath.homogenize(center)
if center not in subspace:
raise ValueError("Center must lie in subspace.")
if subspace.at_infinity():
raise ValueError(
'Can not compute regularized "orthonormal" basis of subspace at infinity.'
)
return subspace_from_affine_points_and_directions(
pmath.dehomogenize(center),
*subspace.affine_point_and_directions[1:],
atol=subspace.atol,
rtol=subspace.rtol,
)
[docs]def dehomogenize_subspace(subspace):
"""Normalize last entries of basis vectors not at infinity to 1.
Parameters
----------
subspace : ddg.geometry.Subspace
Returns
-------
ddg.geometry.Subspace
The same subspace with normalized basis vectors.
"""
entries = subspace.matrix[-1]
diag = []
for e in entries:
if nonexact.isclose(e, 0, atol=subspace.atol, rtol=subspace.rtol):
diag.append(1)
else:
diag.append(1 / e)
return subspace_from_columns(
subspace.matrix @ np.diag(diag), atol=subspace.atol, rtol=subspace.rtol
)
def _dualize(f, subspace, ambient_space=None):
"""Helper function for subspace dualization.
Parameters
----------
f : Callable
Function taking a matrix whose columns describe a subspace as its
only argument and returns a similar matrix describing the dual.
subspace : ddg.geometry.Subspace
Subspace to dualize.
ambient_space : ddg.geometry.Subspace (default=None)
Other subspace containing this subspace that should be used as
ambient space when dualizing. None means whole space.
Returns
-------
ddg.geometry.Subspace
"""
if ambient_space is None:
return subspace_from_columns(f(subspace.matrix))
if not (subspace <= ambient_space):
raise ValueError("Subspace must be contained in ambient space.")
return intersect_subspaces(_dualize(f, subspace), ambient_space)
[docs]def whole_space(dimension, **kwargs):
"""Create whole space of a certain dimension as a trivial subspace.
Parameters
----------
Dimension : int
Dimension of the space.
**kwargs
Keyword arguments other than points to be passed to Subspace.
Returns
-------
ddg.geometry.Subspace
See also
--------
ddg.geometry.Subspace
"""
return subspace_from_columns(np.eye(dimension + 1), **kwargs)
[docs]def coordinate_hyperplane(n, zero_component=-1, origin_index=-1):
""" "Equatorial plane" in n-dimensional space.
Returns the coordinate hyperplane {x_i = 0}, where `i=zero_component`.
Parameters
----------
n : int
Ambient dimension
zero_component : int (default=-1)
Which component should be 0 in this hyperplane in affine coordinates.
origin_index : int (default=-1)
Where to insert the "origin" (the point not at infinity)
Returns
-------
ddg.geometry.Subspace
"""
zc = zero_component % n
oi = origin_index % n
H = np.eye(n)
H = np.delete(H, zc, axis=1)
H = np.insert(H, oi, np.zeros(n), axis=1)
H = np.vstack([H, la.e(oi, n)])
return subspace_from_columns(H)
def join_subspaces(subspace1, subspace2):
"""Construct the join of two _subspaces.
Parameters
----------
subspace1 : ddg.geometry.Subspace
subspace2 : ddg.geometry.Subspace
Returns
-------
ddg.geometry.Subspace
Raises
------
ValueError
If ambient_dimension of the subspaces does not match.
See also
--------
ddg.geometry.join
To join an arbitrary number of subspaces and other objects.
"""
if subspace1.ambient_dimension != subspace2.ambient_dimension:
raise ValueError(
"Ambient dimension mismatch: "
f"{subspace1.ambient_dimension} and "
f"{subspace2.ambient_dimension}."
)
atol, rtol = nonexact.combine_tols(subspace1, subspace2)
# The tolerances given here are used when checking linear dependence of basis
# vectors.
result = Subspace(*(subspace1.points + subspace2.points), atol=atol, rtol=rtol)
# We reset the tolerances because the general convention is that functions return
# objects with default tolerances.
result.atol = None
result.rtol = None
return result
def intersect_subspaces(subspace1, subspace2):
"""Construct the intersection of two linear _subspaces.
Parameters
----------
subspace1 : ddg.geometry.Subspace
subspace2 : ddg.geometry.Subspace
Returns
-------
ddg.geometry.Subspace
Raises
------
ValueError
If ambient_dimension of the subspaces does not match.
See also
--------
ddg.geometry.intersect
To intersect an arbitrary number of subspaces and other objects.
"""
if subspace1.ambient_dimension != subspace2.ambient_dimension:
raise ValueError(
"Ambient dimension mismatch: "
f"{subspace1.ambient_dimension} and "
f"{subspace2.ambient_dimension}."
)
atol, rtol = nonexact.combine_tols(subspace1, subspace2)
points1 = subspace1.dual().matrix
points2 = subspace2.dual().matrix
points = np.append(points1, points2, axis=1)
points = la.col_basis(points, atol=atol, rtol=rtol)
result = subspace_from_columns(points).dual()
return result
[docs]def least_square_subspace(points, k):
"""
Computes the k-dimensional subspace closest to the affine
coordinates of N points.
Parameters
----------
points : sequence of ddg.geometry.Point
Points containing homogeneous coordinates.
k : int
Dimension of the subspace.
Returns
-------
subspace : ddg.geometry.Subspace
k-dimensional projective subspace.
See Also
--------
:py:func:`ddg.geometry.least_square_subspace_from_affine_points`
"""
return least_square_subspace_from_affine_points([p.affine_point for p in points], k)
[docs]def least_square_subspace_from_affine_points(points, k):
"""
Computes the k-dimensional subspace closest to N affine
points. See below for more details.
Parameters
----------
points : sequence of numpy.ndarray, each of shape (n,)
Affine coordinates of N points.
k : int
Dimension of the subspace.
Returns
-------
subspace : ddg.geometry.Subspace
k-dimensional projective subspace.
Raises
------
ValueError
If k is strictly greater than the ambient projective
dimension n.
Notes
-----
The returned subspace is *best-fitting* k-dimensional subspace
to the N given points in the sense that it minimizes the sum of
the squares of the Euclidean distances of the points to it.
See Also
--------
:py:func:`least_square_subspace`
"""
points_span = subspace_from_affine_points(*points)
d = points_span.ambient_dimension
r = points_span.dimension
if k > d:
raise ValueError(
f"Can not return subspace with dimension {k} and ambient dimension {d}. The"
" dimension of subspace must be less than or equal the ambient dimension."
)
elif k == d:
subspace = whole_space(k)
elif k == r:
subspace = points_span
else:
# Compute barycenter of the given points
G = np.array(points).sum(axis=0) / np.array(points).shape[0]
# Run SVD
*_, vh = np.linalg.svd(points - G)
directions = la.matrix_to_arrays(vh[:k, :].T)
subspace = subspace_from_affine_points_and_directions(G, directions)
return subspace