"""Quadrics"""
from textwrap import indent
from functools import total_ordering
import numpy as np
import sympy
from ddg.math import symmetric_matrices as sm
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.euclidean import embed
import ddg.geometry.subspaces as subspaces
import ddg.datastructures.nets.utils as nutils
from ddg.datastructures.nets.net_generators.jacobi_elliptic_curve import (
jacobi_elliptic_curve
)
from ddg.abc import LinearTransformable, NonExact
[docs]@total_ordering
class Quadric(LinearTransformable, NonExact):
"""Quadric in a projective space.
Parameters
----------
Q : np.ndarray
Symmetric matrix representing the quadratic form in homogeneous
coordinates.
The matrix is interpreted as the Gram matrix with respect to the given
basis of the containing subspace.
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.
Raises
------
ValueError
* If Q is not symmetric
* If the dimension of the subspace does not match the shape of Q.
"""
def __init__(self, Q, subspace=None, atol=None, rtol=None):
NonExact.__init__(self, atol=atol, rtol=rtol)
if not sm.is_symmetric(Q, atol=self.atol, rtol=self.rtol):
raise ValueError("Not a symmetric matrix.")
self.matrix = Q
# 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(f"Dimension mismatch: Subspace is "
f"{subspace.dimension}D, but should be "
f"{np.shape(self.matrix)[0]-1}D.")
self.subspace = subspace
LinearTransformable.__init__(self, self.ambient_dimension+1)
@property
def ambient_dimension(self):
return self.subspace.ambient_dimension
[docs] def dimension(self, real=True):
"""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
elif (self.signature() == sm.Signature(np.shape(self.matrix)[0], 0)
and real):
return -1
elif not self.is_degenerate:
return self.subspace.dimension - 1
restriction = intersect_quadric_subspace(
self,
self.non_degenerate_subspace
)
rd = restriction.dimension(real=real)
# 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
[docs] def at_infinity(self, affine_component=-1):
"""Whether or not quadric is at infinity.
Parameters
----------
affine_component : int (default=-1)
Returns
-------
bool
"""
return self.subspace.at_infinity(affine_component)
[docs] def coordinates(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.Subspace):
v = v.matrix
return la.coordinates(v, self.subspace.matrix, atol=self.atol,
rtol=self.rtol)
[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.coordinates(v)
w = self.coordinates(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.coordinates(S1)
S2 = self.coordinates(S2)
eps = S1.T @ self.matrix @ S2
return np.allclose(eps, 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, affine_component=-1):
"""Alias for :py:func:`signature`.
"""
return signature(self, subspace=subspace, affine=affine,
affine_component=affine_component)
[docs] def normalize(self, affine=False, affine_component=-1):
"""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, affine_component).
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 (according to ``affine_component``). This means
that the returned AffineSignature instance will always have
affine_component -1.
Parameters
----------
affine : bool (default=False)
affine_component : int (default=-1)
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:
i = affine_component
# Compute a regularized basis to ensure the transformation is
# affine
point, directions = self.subspace.affine_point_and_directions(i)
B_regularized = np.column_stack(
[embed(d, affine_component=i) for d in directions]
+ [pm.homogenize(point, affine_component=i)]
)
# Express the quadric in the new coordinates
# B_regularized = self.subspace.matrix @ S
S = self.coordinates(B_regularized)
# affine component always -1
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):
mat1 = np.array(self.matrix)
if isinstance(other, Quadric):
if self.subspace != other.subspace:
return False
mat2 = other.matrix
atol = max(self.atol, other.atol)
rtol = max(self.rtol, other.rtol)
# Coordinate transformation
S = la.coordinates(self.subspace.matrix, other.subspace.matrix,
atol=atol, rtol=rtol)
mat2 = S.T @ mat2 @ S
else:
mat2 = np.array(other)
atol = self.atol
rtol = self.rtol
if mat1.shape != mat2.shape:
return False
# check if the matrices are the same up to a scalar multiple
# find a nonzero entry of both
idx = np.nonzero(np.bitwise_and(np.where(mat1!=0, True, False),
np.where(mat2!=0, True, False)))
i,j = np.transpose(idx)[0]
mat1 = mat1 * (mat2[i,j]/mat1[i,j])
return np.isclose(mat1, mat2, atol=atol, rtol=rtol).all()
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):
degen_prefix = "degenerate " if self.is_degenerate else ""
if self.matrix.shape == (3, 3):
s = "conic"
else:
s = "quadric"
sections = []
sections.append(
f"{self.dimension()}D {degen_prefix}{s} in "
f"{self.ambient_dimension}D projective space"
)
sections.append(f"signature: {self.signature()}")
sections.append(f"matrix:\n" + indent(str(self.matrix), ' '))
if not np.array_equal(self.subspace.matrix,
np.eye(self.ambient_dimension + 1)):
sections.append(f"basis of containing subspace:\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.
If one is given as a quadric and the other as an array, the newly
created one will have matching affine component to the
given one, unless overridden by `kwargs`.
**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
super().__init__(Q1.ambient_dimension+1)
@property
def ambient_dimension(self):
return self.Q1.ambient_dimension
[docs] def quadric(self, u):
"""Return a quadric in the pencil.
The pencil is parametrized as [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).
Returns
-------
Quadric
"""
shape = np.shape(u)
if len(shape) == 0:
if np.isinf(u):
u1, u2 = 0., 1.
else:
u1, u2 = 1., u
else:
u1, u2 = u
Q = np.multiply(u1, self.Q1.matrix) + np.multiply(u2, self.Q2.matrix)
return Quadric(Q, subspace=self.Q1.subspace)
[docs] def roots(self):
"""Find values for which quadrics in the pencil are degenerate.
This function calculates the roots of the polynomial
`det(Q1 + lambda*Q2)` and their multiplicity.
Returns
-------
dict {complex or inf: int}
Dictionary whose keys are the roots as complex numbers (or
infinity) and whose values are their multiplicity.
"""
lamda = sympy.Symbol('lamda')
Q = sympy.Matrix(self.Q1.matrix) + lamda * sympy.Matrix(self.Q2.matrix)
Q = sympy.nsimplify(Q)
p = sympy.Poly(Q.det(), lamda, domain='QQ')
roots = sympy.roots(p)
# alternatively use ground_roots:
# roots = sympy.polys.polytools.ground_roots(Q.det())
roots = {complex(root):mult for root,mult in roots.items()}
N = np.shape(self.Q1.matrix)[0]
sum_mult = sum(roots.values())
if N != sum_mult:
roots[np.inf] = N - sum_mult
return roots
[docs] def degenerate_quadrics(self, real=True):
"""Calculate degenerate quadrics in the pencil.
Parameters
----------
real : bool (default=True)
If True, only use real roots.
Returns
-------
list of Quadric
"""
roots = self.roots().keys()
if real:
roots = sorted(
[np.real(root) for root in roots if np.isreal(root)]
)
return [
self.quadric([1, root]) if not np.isinf(root) else self.quadric([0,1])
for root in roots
]
[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` lies on `absolute`.
"""
# TODO maybe check if center is in absolute and set r=0 in that case to
# avoid numerical issues?
b = absolute.inner_product
c = center.point
return cayley_klein_horosphere(center, radius * b(c, c), absolute)
[docs]def cayley_klein_horosphere(center, radius, absolute):
"""Create Cayley-Klein Horosphere 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
"""
A = absolute.matrix
c = absolute.coordinates(center)
return Quadric(
A.T @ np.outer(c, c) @ A - radius * A,
subspace=absolute.subspace
)
[docs]def axis(cone, affine_component=-1):
"""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
affine_component : int (default=-1)
Returns
-------
Subspace
Raises
------
ValueError
If Quadric does not have the signature mentioned above.
"""
cone = cone.normalize(affine=True, affine_component=affine_component)
ad = cone.subspace.dimension
supported_signatures = [
sm.AffineSignature(ad - 1, 1, 1, affine_component_entry=ac)
for ac in [1, -1, 0]
]
sgn = cone.signature(affine=True, affine_component=affine_component)
if sgn not in supported_signatures:
raise ValueError(
f"Quadric is not a non-parabolic cone."
)
return subspaces.Subspace(*cone.subspace.points[-2:])
[docs]def normalization(quadric, affine=False, affine_component=-1):
"""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.
affine_component : int (default=-1)
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:
i = affine_component
# Compute a regularized basis to ensure the transformation is affine
point, directions = quadric.subspace.affine_point_and_directions(i)
B_regularized = np.column_stack(
[embed(d, affine_component=i) for d in directions]
+ [pm.homogenize(point, affine_component=i)]
)
# 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)
# affine component always -1
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
"""
meet = subspaces.meet(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.
"""
join = subspaces.join(quadric.subspace, subspace)
jd = join.dimension
qsd = quadric.subspace.dimension
sd = subspace.dimension
# subspace and quadric.subpspace are disjoint
if jd == qsd + sd + 1:
Q = np.block([
[quadric.matrix, np.zeros((qsd + 1, jd - qsd)) ],
[np.zeros((jd - qsd, qsd + 1)), np.zeros((jd - qsd, jd - qsd))]
])
return Quadric(Q, subspace=join)
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]@NonExact.nonexact_function
def intersect_quadrics(quadric1, quadric2, affine=False, affine_component=-1,
atol=1e-9, rtol=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.
affine_component : int (default=-1)
Used to dehomogenize the output curve.
atol : float (default=1e-9)
rtol : float(default=0.0)
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.
This function uses the global tolerance defaults if `atol` or `rtol` are
set to None. See ddg.abc.NonExact for details.
"""
if not np.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 np.allclose(Q1 - np.diag(diag1), 0, atol=atol, rtol=rtol)
or not np.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 np.isclose(k, 1, atol=atol, rtol=rtol):
k = 1
elif np.isclose(k, 0, atol=atol, rtol=rtol):
k = 0
curve = jacobi_elliptic_curve(k)
nutils.homogenize(curve)
curve.transform(quadric1.subspace.matrix)
if affine:
nutils.dehomogenize(curve, affine_component=affine_component)
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 "
f"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, affine_component=-1):
"""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
affine_component : int (default=-1)
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,
affine_component=affine_component)[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,
affine_component=affine_component)