"""Projective quadric module.
Contains classes for quadrics and pencils of quadrics and utility functions.
"""
from textwrap import indent
import numpy as np
import sympy
import ddg.datastructures.nets.utils as nutils
import ddg.geometry.subspaces as subspaces
from ddg import nonexact
from ddg.abc import LinearTransformable
from ddg.datastructures.nets.net_generators.jacobi_elliptic_curve import (
jacobi_elliptic_curve,
)
from ddg.datastructures.nets.utils import compose
from ddg.geometry.abc import Embeddable
from ddg.math import inner_product as ip
from ddg.math import linalg as la
from ddg.math import projective as pm
from ddg.math import symmetric_matrices as sm
from ddg.math.euclidean import embed
[docs]class Quadric(Embeddable, LinearTransformable):
r"""Quadric in real projective space.
Let `q` be a quadratic form on **R**\ ^(n+1). The quadric corresponding to
`q` is the set ::
{ [x] in RP^n : q(x) = 0 }.
Parameters
----------
matrix : np.ndarray of shape (n, n)
Symmetric matrix representing the quadratic form in homogeneous
coordinates.
The matrix is interpreted as the Gram matrix with respect to the given
basis of `subspace`.
Non-float dtypes will be converted to float.
subspace : ddg.Subspace or iterable of numpy.ndarray (default=None)
Subspace containing the quadric. Can be given as a list of homogeneous
coordinate vectors.
If None is given, uses whole space with standard basis.
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
----------
matrix : numpy.ndarray of shape (n, n)
subspace : Subspace
ambient_dimension
dimension
dimension_complex
rank
corank
is_degenerate
singular_subspace
non_degenerate_subspace
atol, rtol
Raises
------
ValueError
* If `matrix` is not symmetric
* If the dimension of the subspace does not match the shape of
`matrix`.
Notes
-----
This class supports the following operators:
- ``==``, meaning equality of the quadratic forms.
- ``<=``, meaning that one quadratic form is a restriction of the other.
- ``in``, meaning set membership.
"""
def __init__(self, matrix, subspace=None, atol=None, rtol=None):
self.atol = atol
self.rtol = rtol
if matrix.dtype == complex:
raise TypeError("Complex quadrics are not currently supported.")
if not sm.is_symmetric(matrix, atol=self.atol, rtol=self.rtol):
raise ValueError("Not a symmetric matrix.")
# Convert integer matrices to float matrices to avoid confusing behavior
self.matrix = matrix.astype(float)
# None means whole space
if subspace is None:
subspace = subspaces.whole_space(np.shape(self.matrix)[0] - 1)
# Given as list of points
elif not isinstance(subspace, subspaces.Subspace):
subspace = subspaces.Subspace(*subspace)
if subspace.dimension != np.shape(self.matrix)[0] - 1:
raise ValueError(
"Dimension mismatch: Subspace is "
f"{subspace.dimension}D, but should be "
f"{np.shape(self.matrix)[0]-1}D."
)
self.subspace = subspace
@property
def ambient_dimension(self):
return self.subspace.ambient_dimension
@property
def dimension(self):
"""The (real) dimension of the quadric.
More precisely, the maximum of the dimensions of the manifolds
contained in the set of real points in the quadric.
Ignores affine coordinates! I.e. we don't distinguish between the
quadric lying at infinity or not.
Returns
-------
int
Notes
-----
This function was fully intuited. There is no guarantee that this gives
correct results.
"""
# Reasoning for why this might be correct:
# Let n be the dimension of the space the quadric is contained in. The
# dimension of a non-degenerate quadric is n-1, except if the quadratic
# form is positive definite, in which case it is empty with dimension
# -1. A degenerate quadric is the family of lines connecting all points
# in its singular subspace with all points in the quadric obtained by
# restricting the quadratic form to a subspace complementary to its
# kernel. Therefore, the dimension is the sum of the dimension of these
# two objects + 1.
# Fully degenerate, whole space. Treated separately because the
# restriction step below doesn't work in this case.
if self.rank == 0:
return self.subspace.dimension
# positive definite, empty
m = self.subspace.dimension
if self.signature() == sm.Signature(m + 1, 0):
return -1
if not self.is_degenerate:
return self.subspace.dimension - 1
restriction = intersect_quadric_subspace(self, self.non_degenerate_subspace)
rd = restriction.dimension
# This formula works even if the restriction is positive definite /
# empty. The -1 is canceled by the 1.
return rd + self.singular_subspace.dimension + 1
@property
def dimension_complex(self):
"""Complex dimension of the quadric.
This is the same as the real dimension, except that the positive
definite quadric is not empty and not a special case anymore. In fact,
there is only one non-degenerate quadric up to projective
transformation, the one with diagonal ``[1,...,1]``.
Returns
-------
int
"""
if self.rank == 0:
return self.subspace.dimension
return self.subspace.dimension - 1
[docs] def at_infinity(self):
return self.subspace.at_infinity()
def _coordinate_helper(self, v):
"""Get coordinates of vectors in terms of self.subspace's basis.
Parameters
----------
v : Subspace or numpy.ndarray
Subspace contained in self.subspace or Matrix representing such a
subspace (or Point).
Returns
-------
numpy.ndarray of shape (self.subspace.dimension + 1,)
"""
if isinstance(v, subspaces.Point):
v = v.point
elif isinstance(v, subspaces.Subspace):
v = v.matrix
return la.coordinates(v, self.subspace.matrix, atol=self.atol, rtol=self.rtol)
[docs] def embed(self, subspace=None):
if subspace is None:
subspace = subspaces.coordinate_hyperplane(self.ambient_dimension + 1)
return Quadric(
self.matrix,
self.subspace.embed(subspace),
atol=self.atol,
rtol=self.rtol,
)
[docs] def unembed(self, subspace=None):
if subspace is None:
subspace = subspaces.coordinate_hyperplane(self.ambient_dimension)
if not self.subspace <= subspace:
raise ValueError("subspace must contain self.")
return Quadric(self.matrix, subspace=self.subspace.unembed(subspace))
[docs] def inner_product(self, v, w):
"""Inner product defined by / defining the quadric.
Parameters
----------
v, w : numpy.ndarray of shape (ambient_dimension + 1,)
Two vectors that lie in the containing subspace.
Returns
-------
float
``a.T @ matrix @ b``, where `a` and `b` are the coordinates of `v`
and `w` in terms of the given basis of the containing subspace.
"""
v = self._coordinate_helper(v)
w = self._coordinate_helper(w)
return ip.inner_product_from_matrix(self.matrix)(v, w)
[docs] def cayley_klein_distance(self, v, w):
"""Cayley-Klein "distance".
This is defined as ::
b(v, w) ** 2 / (b(v, v) * b(w, w))
Where b is the inner product induced by the absolute. None of the
points can lie on the absolute.
Parameters
----------
v, w : numpy.ndarray or Point
Homogeneous coordinate vectors or Point instances
Returns
-------
float
Raises
------
ValueError
If either of the points is contained in the absolute quadric.
"""
if isinstance(v, subspaces.Point):
v = v.point
if isinstance(w, subspaces.Point):
w = w.point
if v in self or w in self:
raise ValueError(
"The Cayley-Klein distance can not be computed for points in "
"the absolute quadric."
)
b = self.inner_product
return b(v, w) ** 2 / (b(v, v) * b(w, w))
[docs] def conjugate(self, S1, S2):
"""Conjugacy of two subspaces.
Whether two subspaces are conjugate, i.e. if approximately
``inner_product(v, w) == 0`` for all v in S1 and w in S2.
Parameters
----------
S1, S2 : ddg.geometry.subspaces.Subspace or numpy.ndarray
Subspaces contained in self.subspace. Arrays will be interpreted as
matrices whose columns are points in homogeneous coordinates. In
particular, 1D arrays are simply points in homogeneous coordinates.
Returns
-------
bool
"""
S1 = self._coordinate_helper(S1)
S2 = self._coordinate_helper(S2)
eps = S1.T @ self.matrix @ S2
return nonexact.allclose(eps, 0.0, atol=self.atol, rtol=self.rtol)
[docs] def polarize(self, obj):
"""Alias for ``polarize(obj, self)``.
See also
--------
polarize
"""
return polarize(obj, self)
[docs] def signature(self, subspace=None, affine=False):
"""Alias for :py:func:`signature`."""
return signature(self, subspace=subspace, affine=affine)
[docs] def normalize(self, affine=False):
"""Return an equal quadric expressed in coordinates that normalize it.
Returns an equal quadric Q_normalized such that ::
Q_normalized.matrix == sgn.matrix
where sgn = self.signature(affine).
See :py:class:`ddg.math.symmetric_matrices.Signature` and
:py:class:`ddg.math.symmetric_matrices.AffineSignature` for more
information.
If ``affine=True``, ``Q_normalized.subspace.points`` will be in the
form [d1,...,dn, p], where n is the dimension of the subspace,
d1,...,dn are directions at infinity and p is a point in the subspace
not at infinity.
Parameters
----------
affine : bool (default=False)
Returns
-------
Q_normalized : Quadric
"""
if not affine:
sgn, F = sm.projective_normalization(
self.matrix, atol=self.atol, rtol=self.rtol
)
B_new = self.subspace.matrix @ F
else:
# Compute a regularized basis to ensure the transformation is
# affine
point, directions = self.subspace.affine_point_and_directions
B_regularized = np.column_stack(
[embed(d) for d in directions] + [pm.homogenize(point)]
)
# Express the quadric in the new coordinates
# B_regularized = self.subspace.matrix @ S
S = self._coordinate_helper(B_regularized)
sgn, F = sm.affine_normalization(
S.T @ self.matrix @ S, atol=self.atol, rtol=self.rtol
)
B_new = B_regularized @ F
return Quadric(sgn.matrix, subspace=B_new.T)
@property
def rank(self):
"""Rank of self.matrix.
Returns
-------
int
"""
return self.signature().rank
@property
def corank(self):
"""Corank of self.matrix.
Returns
-------
int
"""
return self.signature().zero
@property
def singular_subspace(self):
"""Singular subspace.
Subspace containing all singular points of the quadric (points in the
kernel of ``self.matrix``).
Returns
-------
ddg.geometry.subspaces.Subspace of dimension ``corank - 1`` and the
same ambient dimension as the quadric.
"""
sgn, F = normalization(self)
mask = np.full(np.shape(self.matrix)[0], True)
mask[: self.rank] = False
kernel = np.eye(np.shape(self.matrix)[0])[:, mask]
return subspaces.subspace_from_columns(
self.subspace.matrix @ F @ kernel,
)
@property
def non_degenerate_subspace(self):
"""Non-degenerate subspace.
Complementary subspace of the kernel of ``self.matrix``.
Returns
-------
ddg.geometry.subspaces.Subspace
Has dimension ``rank - 1`` and the same ambient dimension as the
quadric.
"""
sgn, F = normalization(self)
kernel_complement = np.eye(np.shape(self.matrix)[0])[:, : self.rank]
return subspaces.subspace_from_columns(
self.subspace.matrix @ F @ kernel_complement,
)
@property
def is_degenerate(self):
"""Whether the quadric is degenerate.
Returns
-------
bool
"""
return self.corank != 0
[docs] def dualize(self):
B = self.subspace.matrix
# B.T @ B is invertible because the columns of B are linearly
# independent.
dual_basis = B @ np.linalg.inv(B.T @ B)
if not self.is_degenerate:
return Quadric(np.linalg.inv(self.matrix), subspace=dual_basis.T)
else:
sgn, F = normalization(self)
dual_kernel = self.singular_subspace.dualize(self.subspace)
cone = Quadric(F @ sgn.matrix @ F.T, subspace=dual_basis.T)
return intersect_quadric_subspace(cone, dual_kernel)
def __contains__(self, v):
if v not in self.subspace:
return False
else:
return self.conjugate(v, v)
def __eq__(self, other):
if isinstance(other, np.ndarray):
# unsigned int, signed int, float or complex
if other.dtype.kind not in set("uifc"):
return False
if other.shape != self.matrix.shape:
return False
other = Quadric(other)
atol = self.atol
rtol = self.rtol
if isinstance(other, Quadric):
atol, rtol = nonexact.combine_tols(self, other)
else:
return NotImplemented
# Subspace.__eq__ raises an error instead of returning False when
# ambient dimensions don't match. This preserves the behavior of
# Quadric.__eq__. These should be made consistent.
if self.ambient_dimension != other.ambient_dimension:
return False
if self.subspace != other.subspace:
return False
M1 = self.matrix
M2 = other.matrix
# Coordinate transformation
S = la.coordinates(
self.subspace.matrix, other.subspace.matrix, atol=atol, rtol=rtol
)
M2 = S.T @ M2 @ S
# check if the matrices are the same up to a scalar multiple.
M1_nonzero = np.logical_not(nonexact.isclose(M1, 0, atol=atol, rtol=rtol))
M2_nonzero = np.logical_not(nonexact.isclose(M2, 0, atol=atol, rtol=rtol))
if (M1_nonzero != M2_nonzero).any():
return False
if not M1_nonzero.any(): # Need to do this to not get index error
return True
i, j = np.argwhere(M1_nonzero)[0]
return nonexact.allclose((M2[i, j] / M1[i, j]) * M1, M2, atol=atol, rtol=rtol)
def __le__(self, other):
if not isinstance(other, Quadric):
return NotImplemented
if not self.subspace <= other.subspace:
return False
restriction = intersect_quadric_subspace(other, self.subspace)
return self == restriction
def __repr__(self):
argstrings = [repr(self.matrix)]
if not np.array_equal(self.subspace.matrix, np.eye(self.ambient_dimension + 1)):
argstrings.append(f"subspace={repr(self.subspace.points)}")
argstrings = ",\n".join(argstrings)
return "Quadric(\n" + indent(argstrings, " ") + "\n)"
def __str__(self):
# Sections will be joined with newlines later.
sections = []
# First line says either "quadric in nD projective space" or
# "quadric contained in nD subspace in mD projective space".
# If n == 2, "quadric" is replaced by "conic".
first_line = []
if self.subspace.dimension == 2:
first_line.append("conic")
else:
first_line.append("quadric")
if self.subspace.codimension > 0:
first_line.append(f"contained in {self.subspace.dimension}D subspace")
first_line.append(f"in {self.ambient_dimension}D projective space")
sections.append(" ".join(first_line))
# "signature: (p, q)" or "signature: (p, q, r)".
sections.append(f"signature: {self.signature()}")
# matrix:
# [[...]
# ...
# [...]]
sections.append("matrix:\n" + indent(str(self.matrix), " "))
# basis of containing subspace (columns):
# [[...]
# ...
# [...]]
# Or:
# basis of ambient space (columns):
# [[...]
# ...
# [...]]
if not np.array_equal(self.subspace.matrix, np.eye(self.ambient_dimension + 1)):
if self.subspace.codimension > 0:
subspace_heading = "basis of containing subspace (columns):"
else:
subspace_heading = "basis of ambient space (columns):"
sections.append(
subspace_heading + "\n" + indent(str(self.subspace.matrix), " ")
)
return "\n".join(sections)
[docs]class Pencil(LinearTransformable):
"""Pencil of quadrics in a projective space.
Parameters
----------
Q1, Q2 : Quadric or numpy.ndarray
Quadrics spanning the pencil. If given as arrays, they will be
created.
**kwargs : dict
Keyword arguments to be passed to Quadric during creation of Q1 and Q2,
if they are given as arrays.
Raises
------
ValueError
* If the geometries of the quadrics don't match
* If the quadrics are not in the same subspace
See also
--------
ddg.geometry.quadrics.Quadric
"""
def __init__(self, Q1, Q2, **kwargs):
quadrics = [q for q in [Q1, Q2] if isinstance(q, Quadric)]
if len(quadrics) == 1:
kwargs.setdefault("subspace", quadrics[0].subspace)
if isinstance(Q1, np.ndarray):
Q1 = Quadric(Q1, **kwargs)
if isinstance(Q2, np.ndarray):
Q2 = Quadric(Q2, **kwargs)
if Q1.subspace != Q2.subspace:
raise ValueError("Quadrics are not in the same subspace.")
self.Q1 = Q1
self.Q2 = Q2
@property
def ambient_dimension(self):
return self.Q1.ambient_dimension
@property
def subspace(self):
return self.Q1.subspace
[docs] def matrix(self, u):
"""Return the matrix of a quadric in the pencil.
The pencil is parametrized as::
RP^1 -> pencil, [u1, u2] -> [u1*Q1 + u2*Q2].
We identify [1, a] <-> a and [0, 1] <-> inf.
Parameters
----------
u : array_like (float, float) or float
A finite float will be interpreted as (1, u) and inf will be
interpreted as (0, 1). -inf is interpreted as (0, -1).
Returns
-------
numpy.ndarray of shape (n, n)
"""
shape = np.shape(u)
if len(shape) == 0:
if np.isinf(u):
u1, u2 = 0.0, np.sign(u) * 1.0
else:
u1, u2 = 1.0, u
else:
u1, u2 = u
return u1 * self.Q1.matrix + u2 * self.Q2.matrix
[docs] def quadric(self, u):
"""Return a quadric in the pencil. Complex quadrics are not supported.
The pencil is parametrized as::
RP^1 -> pencil, [u1, u2] -> [u1*Q1 + u2*Q2].
We identify [1, a] <-> a and [0, 1] <-> inf.
Parameters
----------
u : array_like (float, float) or float including inf.
A finite float will be interpreted as (1, u) and inf will be
interpreted as (0, 1). -inf is interpreted as (0, -1).
Returns
-------
Quadric
"""
return Quadric(self.matrix(u), subspace=self.Q1.subspace)
def _roots(self, real=True):
"""
Calculate real or complex roots of the polynomial ``det(Q1 + lambda * Q2)``
without any infinity shenanigans.
Parameters
----------
real : bool
Returns
-------
dict {complex: int} or {float: int}
Dictionary mapping roots to their multiplicity.
"""
# For other possibilities to find roots of polynomials, see
# https://docs.sympy.org/dev/guides/solving/find-roots-polynomial.html
lambda_ = sympy.Symbol("lambda")
Q = sympy.Matrix(self.Q1.matrix) + lambda_ * sympy.Matrix(self.Q2.matrix)
Q = sympy.nsimplify(Q)
p = sympy.Poly(Q.det(), lambda_)
# float() / complex() converts sympy's (Complex)RootOf objects or expressions
# to float/complex.
if real:
roots = {
float(root): multiplicity
for root, multiplicity in p.real_roots(multiple=False)
}
else:
roots = {
complex(root): multiplicity
for root, multiplicity in p.all_roots(multiple=False)
}
return roots
[docs] def roots(self):
"""Find values for which quadrics in the pencil are degenerate.
This function calculates the complex roots of the polynomial
``P := det(Q1 + lambda * Q2)`` and their multiplicity. If the sum of the
multiplicities less than ``N := subspace.dimension + 1`` ,which means that `Q2`
is degenerate, ``complex('inf')`` is added and assigned the "remaining"
multiplicity ``N - deg(P)``.
Returns
-------
dict {complex: int}
Dictionary whose keys are the roots as complex numbers (including
infinity) and whose values are their multiplicity.
"""
roots = self._roots(real=False)
N = self.Q1.subspace.dimension + 1
sum_mult = sum(roots.values())
if N != sum_mult:
roots[complex("inf")] = N - sum_mult
return roots
[docs] def degenerate_quadrics(self):
"""Calculate real degenerate quadrics in the pencil.
Returns
-------
list of Quadric
"""
result = [self.quadric(root) for root in self._roots(real=True)]
if self.Q2.is_degenerate:
result.append(self.Q2)
return result
[docs]def touching_cone(p, quadric, in_subspace=False):
"""Calculate touching cone.
The touching cone of a quadric Q and a point P is the set of lines through
P that are tangent to Q.
Note that by default, this function looks at tangency in the ambient space
as opposed to within `quadric.subspace`. This means that if the quadric is
contained in a proper subspace, the function returns the join of `p` and
`quadric`.
Parameters
----------
p : numpy.ndarray of shape (quadric.ambient_dimension+1,) or Point
Point in homogeneous coordinates. Must lie in `quadric.subspace` if
in_subspace is True.
quadric : Quadric
in_subspace : bool (default=False)
Returns
-------
Quadric or Subspace
Raises
------
ValueError
If `in_subspace` is True and `p` is not contained in
`quadric.subspace`.
"""
if not isinstance(p, subspaces.Point):
p = subspaces.Point(p)
if not in_subspace and quadric.subspace.codimension > 0:
return join_quadric_subspace(quadric, p)
if in_subspace and p not in quadric.subspace:
raise ValueError(
f"Point {p} not contained in quadric.subspace with matrix:\n"
f"{quadric.subspace.matrix}"
)
return cayley_klein_sphere(p, 1, quadric)
[docs]def cayley_klein_sphere(center, radius, absolute):
"""Create a Cayley-Klein sphere from center, radius and absolute.
Let b be the inner product induced by the `absolute`. A Cayley-Klein sphere
is defined as the solution set of::
b(center, x) ** 2 - radius * b(center, center) * b(x, x) = 0
Note that if `center` lies on the `absolute`, you will just get the polar
hyperplane of `center`.
Parameters
----------
center : Point
radius : float
absolute : Quadric
Returns
-------
Quadric
Raises
------
ValueError
If `center` is not in `absolute.subspace`.
"""
# TODO maybe check if center is in absolute and set r=0 in that case to
# avoid numerical issues?
if center not in absolute.subspace:
raise ValueError(
"Center of Cayley-Klein sphere must lie in the same subspace as "
"the absolute quadric."
)
b = absolute.inner_product
c = center.point
return generalized_cayley_klein_sphere(center, radius * b(c, c), absolute)
[docs]def generalized_cayley_klein_sphere(center, radius, absolute):
"""Create Generalized Cayley-Klein sphere from center, "radius" and
absolute.
Let b be the inner product induced by the `absolute`. A Cayley-Klein sphere
is defined as the solution set of::
b(center, x) ** 2 - radius * b(x, x) = 0
This allows centers to lie on the absolute quadric, at the expense of
depending on the representative vector of `center`. The equation of a
regular Cayley-Klein sphere can still be obtained by using the new radius
``radius * b(center, center)``.
Parameters
----------
center : Point
radius : float
absolute : Quadric
Returns
-------
Quadric
Raises
------
ValueError
If `center` is not in `absolute.subspace`.
"""
if center not in absolute.subspace:
raise ValueError(
"Center of Cayley-Klein sphere must lie in the same subspace as "
"the absolute quadric."
)
A = absolute.matrix
c = absolute._coordinate_helper(center)
return Quadric(A.T @ np.outer(c, c) @ A - radius * A, subspace=absolute.subspace)
[docs]def axis(cone):
"""Find the axis of a non-parabolic cone.
More precisely: The axis of a non-parabolic quadric with signature
(n-1, 1, 1), where n is the dimension of the containing subspace.
Parameters
----------
cone : Quadric
Returns
-------
Subspace
Raises
------
ValueError
If Quadric does not have the signature mentioned above.
"""
cone = cone.normalize(affine=True)
ad = cone.subspace.dimension
supported_signatures = [
sm.AffineSignature(ad - 1, 1, 1, last_entry=le) for le in [1, -1, 0]
]
sgn = cone.signature(affine=True)
if sgn not in supported_signatures:
raise ValueError("Quadric is not a non-parabolic cone.")
return subspaces.Subspace(*cone.subspace.points[-2:])
[docs]def normalization(quadric, affine=False):
"""Signature and normalizing transformation.
Approximately satisfies::
F.T @ quadric.matrix @ F == sgn.matrix
See :py:class:`ddg.math.symmetric_matrices.Signature` and
:py:class:`ddg.math.symmetric_matrices.AffineSignature` for more
information.
Note that this is a transformation within the subspace as opposed to a
transformation of the ambient space and thus cannot be used with
``quadric.transform``.
Parameters
----------
quadric : Quadric
affine : bool (default=False)
Whether to normalize to projective or affine normal form.
Returns
-------
sgn : Signature or AffineSignature
F : numpy.ndarray of shape (ambient_dimension+1, dimension+2)
Transformation that normalizes the quadric
"""
if not affine:
sgn, F = sm.projective_normalization(
quadric.matrix, atol=quadric.atol, rtol=quadric.rtol
)
else:
# Compute a regularized basis to ensure the transformation is affine
point, directions = quadric.subspace.affine_point_and_directions
B_regularized = np.column_stack(
[embed(d) for d in directions] + [pm.homogenize(point)]
)
# Express the quadric in the new coordinates
# B_regularized = quadric.subspace.matrix @ S
S = la.coordinates(
B_regularized, quadric.subspace.matrix, atol=quadric.atol, rtol=quadric.rtol
)
sgn, F = sm.affine_normalization(
S.T @ quadric.matrix @ S, atol=quadric.atol, rtol=quadric.rtol
)
# Express the transformation in the original coordinates again
F = S @ F
return sgn, F
[docs]def intersect_quadric_subspace(quadric, subspace):
"""Intersect quadric with subspace
Parameters
----------
quadric : Quadric
subspace : Subspace
Returns
-------
Quadric
See also
--------
ddg.geometry.intersection.intersect
To intersect more than two objects of various types.
"""
meet = subspaces.intersect_subspaces(quadric.subspace, subspace)
gram = ip.gram_matrix(quadric.inner_product, meet.matrix)
return Quadric(gram, subspace=meet)
[docs]def join_quadric_subspace(quadric, subspace):
"""Join a quadric and a subspace.
The result will be the projective cone with top S and basis Q.
Parameters
----------
quadric : Quadric
subspace : Subspace
Returns
-------
Quadric or Subspace
Raises
------
ValueError
If `quadric.subspace` and `subspace` are not disjoint.
NotImplementedError
* If quadric.subspace and subspace are not disjoint and subspace is not
a point
* If the join is an object like a "solid cone", which we can't
represent.
See also
--------
ddg.geometry.intersection.join
To join more than two objects of various types.
"""
join = subspaces.join_subspaces(quadric.subspace, subspace)
jd = join.dimension
qsd = quadric.subspace.dimension
sd = subspace.dimension
# subspace and quadric.subpspace are disjoint
if jd == qsd + sd + 1:
# fmt: off
Q = np.block(
[[quadric.matrix, np.zeros((qsd + 1, jd - qsd))],
[np.zeros((jd - qsd, qsd + 1)), np.zeros((jd - qsd, jd - qsd))]]
)
# fmt: on
return Quadric(Q, subspace=(quadric.subspace.points + subspace.points))
if not isinstance(subspace, subspaces.Point):
# Not sure how hard this is to generalize
raise NotImplementedError(
"Join of a quadric with a subspace not disjoint from "
"quadric.subspace is currently only implemented for points."
)
if subspace not in quadric:
# Let a = subspace, Q = quadric.matrix, q the bilinear form. Then:
# join is NOT equal to quadric.subspace
# <=> There is a line through a that does not intersect quadric.
# The representative matrix of the touching cone is
# C = Q a a.T Q - (a.T Q a) Q, up to scaling. So Q is a subset of
# {x | x.T C x <= 0} or {x | x.T C x >= 0}, depending on the
# representative matrix. If there is a b with the opposite sign, the
# whole line [a + lambda b] lies in the complement of the respective
# set, so it does not intersect the quadric.
tc = touching_cone(subspace, quadric, in_subspace=True)
if tc.signature().is_indefinite:
raise NotImplementedError(
"The join is a 'solid' cone, which we can't represent. Use "
"touching_cone with in_subspace=True to obtain its boundary."
)
return quadric.subspace
else: # Point is on the quadric.
if subspace in quadric.singular_subspace:
# If the point (=: a) is in the kernel, we just get the quadric
# back:
# 1. join is a subset of quadric:
# A point in the join is on a line connecting a with a point on
# the quadric. Plug this into the quadratic form.
# 2. quadric is a subset of join:
# Every point p on the quadric lies on a line [mu a + nu p].
return quadric
else:
# This is technically not correct because the tangent plane is not
# contained in the join. But since it is the only thing not
# contained, we ignore it.
# Let q be the bilinear form of the quadric.
# join is the whole space
# <=> All lines [mu a + nu b] through a intersect quadric \ {a}
# (except lines in the polar hyperplane, which we ignore).
# <=> For all b there are mu, nu such that
# 2 mu nu q(a, b) + mu^2 q(b, b) = 0.
# We assume q(a, b) != 0 because a is not in ker(q) and b is
# not polar to a.
# Choose mu=1, nu=-q(b, b)/(2 * q(a, b)) to get the
# intersection point with quadric \ {a}.
return quadric.subspace
[docs]def intersect_quadrics(quadric1, quadric2, affine=False, atol=1e-9, rtol=0.0):
r"""Intersection curve of certain quadrics.
Currently only works for two quadrics that are both contained in a pencil
::
u1 * diag([1, 1, 0, -1]) + u2 * diag([0, k, 1, -1])
\_______ ________/ \_______ ________/
\/ \/
D1 D2
for some k in [0, 1]. The quadrics can be contained in any 3D subspace as
long as they are given w.r.t. the same basis.
The function checks if both quadrics are contained in a pencil like this.
If you *know* your two quadrics with matrices Q1, Q2 are contained in a
transformed version of such a pencil, i.e. ::
u1 * F.T @ D1 @ F + u2 * F.T @ D2 @ F
for some known invertible F, you can make it work as well: Let B be
quadric1.subspace.matrix. Then the quadrics ::
Quadric(inv(F.T) @ Q1 @ inv(F), subspace=(B @ inv(F)).T)
Quadric(inv(F.T) @ Q1 @ inv(F), subspace=(B @ inv(F)).T)
can be intersected with this function. Automatic detection of this is not
yet implemented.
Parameters
----------
quadric1, quadric2 : Quadric
affine : bool (default=False)
Whether the resulting curve should output affine or homogeneous
coordinates.
atol : float (default=1e-9)
rtol : float (default=0.0)
This function uses the global tolerance defaults if `atol` or `rtol` are set to
None. See :py:mod:`ddg.nonexact` for details.
Returns
-------
NetCollection
Raises
------
ValueError
* If the two quadrics are given in different coordinate systems.
* If quadric1 and quadric2 do not span the above-mentioned pencil.
Notes
-----
For visualization purposes, not all subspace bases work well, because the
curve might pass through infinity.
"""
if not nonexact.allclose(
quadric1.subspace.matrix, quadric2.subspace.matrix, atol=atol, rtol=rtol
):
raise ValueError("Both quadrics must be expressed in the same coordinates.")
Q1 = quadric1.matrix
Q2 = quadric2.matrix
diag1 = np.diagonal(Q1)
diag2 = np.diagonal(Q2)
if not nonexact.allclose(
Q1 - np.diag(diag1), 0, atol=atol, rtol=rtol
) or not nonexact.allclose(Q2 - np.diag(diag2), 0, atol=atol, rtol=rtol):
raise ValueError("Both quadrics must be given by diagonal matrices.")
# Q1 = diag(diag1) and Q2 = diag(diag2) span the same pencil as
# diag([1, 1, 0, -1]), diag([0, k, 1, -1]) iff
#
# a1 * diag1 + a2 * diag2 = [1, 1, 0, -1]
#
# for some a1, a2 and
#
# b1 * diag1 + b2 * diag2 - k * [0, 1, 0, 0] = [0, 0, 1, -1]
#
# for some b1, b2 and some k in [0, 1].
S1 = np.column_stack((diag1, diag2))
S2 = np.column_stack((diag1, diag2, [0, 1, 0, 0]))
try:
c1 = la.coordinates([1, 1, 0, -1], S1, atol=atol, rtol=rtol)
c2 = la.coordinates([0, 0, 1, -1], S2, atol=atol, rtol=rtol)
# one of the Points does not lie in the subspace
except ValueError:
raise ValueError(
"Pencil does not contain diag([1, 1, 0, -1]) or does not contain "
"diag([0, k, 1, -1]) for some 0 <= k <= 1."
)
k = -c2[-1]
if not 0 - atol <= k <= 1 + atol:
raise ValueError(f"k is {k}, but should be between 0 and 1.")
# Breaks outside of [0, 1] and doesn't work well close to 1.
# The tolerances are chosen so that the scipy.ellipj functions still give
# acceptable results.
if nonexact.isclose(k, 1, atol=atol, rtol=rtol):
k = 1
elif nonexact.isclose(k, 0, atol=atol, rtol=rtol):
k = 0
curve = jacobi_elliptic_curve(k)
curve = nutils.homogenize(curve)
curve = compose(lambda x: quadric1.subspace.matrix @ x, curve)
if affine:
curve = nutils.dehomogenize(curve)
return curve
[docs]def polarize(obj, quadric):
"""Polarize a geometric object with respect to a quadric.
Parameters
----------
obj : Geometric object
Currently supported: Subspace, Quadric.
quadric : Quadric
Returns
-------
type(obj)
Raises
------
TypeError
If the type of object is not supported.
"""
if isinstance(obj, subspaces.Subspace):
return _polarize_subspace(obj, quadric)
elif isinstance(obj, Quadric):
return _polarize_quadric(obj, quadric)
else:
raise TypeError(f"Polarization for {type(obj)} objects not supported (yet).")
def _polarize_subspace(subspace, quadric):
"""Polarize a subspace with respect to a quadric.
If b is the defining bilinear form of the quadric, the polar space of a
subspace S is the set {x | b(x,y) = 0 for all y in S}.
Parameters
----------
subspace : Subspace
quadric : Quadric
Returns
-------
Subspace
Raises
------
ValueError
If subspace is not contained in quadric.subspace
"""
if not subspace <= quadric.subspace:
raise ValueError("Subspace must be contained in quadric.subspace.")
B = quadric.subspace.matrix
Q = quadric.matrix
points = la.coordinates(subspace.matrix, B, atol=subspace.atol, rtol=subspace.rtol)
polar_points = Q @ points
polar_points = ip.get_col_complement(polar_points)
polar_points = B @ polar_points
return subspaces.subspace_from_columns(polar_points)
def _polarize_quadric(quadric1, quadric2):
"""Polarize a quadric with respect to another quadric.
The polar quadric is obtained by polarizing all tangent hyperplanes.
Parameters
----------
quadric1 : Quadric
Quadric to polarize.
quadric2 : Quadric
Quadric with respect to which to polarize.
Returns
-------
Quadric
Raises
------
ValueError
If quadric1 is degenerate.
"""
if quadric1.subspace != quadric2.subspace:
raise ValueError("Both quadrics must be contained in the same subspace.")
Q1 = quadric1.matrix
Q2 = quadric2.matrix
S = la.coordinates(
quadric1.subspace.matrix,
quadric2.subspace.matrix,
atol=quadric1.atol,
rtol=quadric1.rtol,
)
Q2 = S.T @ Q2 @ S
if not quadric1.is_degenerate:
Q_polar = Q2.T @ np.linalg.inv(Q1) @ Q2
return Quadric(Q_polar, subspace=quadric1.subspace.points)
else:
kernel = quadric1.singular_subspace
polar_kernel = quadric2.polarize(kernel)
sgn, F = normalization(quadric1)
cone = Quadric(
Q2.T @ F @ sgn.matrix @ F.T @ Q2, subspace=quadric1.subspace.points
)
return intersect_quadric_subspace(cone, polar_kernel)
[docs]def signature(quadric, subspace=None, affine=False):
"""Signature of a quadric, optionally restricted to a subspace.
Parameters
----------
subspace : Subspace or None (default=None)
Subspace contained in self.subspace
affine : bool (default=False)
Whether to return projective or affine signature
Returns
-------
Signature or AffineSignature
Raises
------
ValueError
If subspace is not contained in self.subspace.
See also
--------
ddg.math.symmetric_matrices.Signature,
ddg.math.symmetric_matrices.AffineSignature
"""
if subspace is None:
return normalization(quadric, affine=affine)[0]
elif not subspace <= quadric.subspace:
raise ValueError("Subspace must be contained in self.subspace.")
else:
restriction = intersect_quadric_subspace(quadric, subspace)
return signature(restriction, affine=affine)