"""Projective quadric module.
Contains classes for quadrics and pencils of quadrics and utility functions.
"""
from itertools import combinations, combinations_with_replacement
from textwrap import indent
import numpy as np
import sympy
import ddg.geometry._subspaces as subspaces
import ddg.geometry.signatures
from ddg import nonexact
from ddg.abc import LinearTransformable
from ddg.geometry import signatures as sg
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
from ddg.math.linalg import nullspace
from ddg.nets._utils import compose, dehomogenize, homogenize
from ddg.nets.parametrizations.jacobi_elliptic_curve import jacobi_elliptic_curve
[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 : numpy.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.geometry.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 : ddg.geometry.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() == sg.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 : ddg.geometry.Subspace or numpy.ndarray
Subspace contained in self.subspace or Matrix representing such a
subspace (or ddg.geometry.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 v.T @ self.matrix @ 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 ddg.geometry.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.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, object_):
"""Alias for ``polarize(object_, self)``.
See also
--------
ddg.geometry.polarize
"""
return polarize(object_, self)
[docs] def polarize_oriented(self, object_):
"""Alias for ``polarize_oriented(object_, self)``.
See also
--------
polarize
"""
return polarize_oriented(object_, 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.geometry.signatures.Signature` and
:py:class:`ddg.geometry.signatures.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 : ddg.geometry.Quadric
"""
if not affine:
sgn, F = sg.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 = sg.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.Subspace of dimension ``corank - 1`` and the
same ambient dimension as the quadric.
"""
sgn, F = quadric_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.Subspace
Has dimension ``rank - 1`` and the same ambient dimension as the
quadric.
"""
sgn, F = quadric_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 dual(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 = quadric_normalization(self)
dual_kernel = self.singular_subspace.dual(ambient_space=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 : ddg.geometry.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.
Attributes
----------
Q1, Q2
ambient_dimension
subspace
Raises
------
ValueError
* If the geometries of the quadrics don't match
* If the quadrics are not in the same subspace
See also
--------
ddg.geometry.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 intersection(self):
"""
Intersection of all quadrics in the pencil
Returns
-------
ddg.geometry.QuadricIntersection
"""
Q1 = self.Q1
Q2 = self.Q2
return QuadricIntersection(Q1, Q2)
[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
-------
ddg.geometry.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 ddg.geometry.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 regular_matrix(self):
"""
Calculates a regular quadric Q2
This is achieved by choosing a lambda for ``lambda Q2 + Q1`` such that::
det(lambda Q2 +Q1) != 0
Together with ``Q1``, the new ``Q2`` spans the same pencil as the original
``Q1, Q2``.
Returns
-------
ddg.geometry.Quadric
Raises
------
NotImplementedError
if the pencil is degenerate and only has roots at infinity
"""
# extract the pencil spanning matrix Q2
Q2 = self.Q2
# check whether Q2 is degenerate (if not, we are done)
if np.linalg.det(Q2.matrix) == 0:
# calculate the roots of the pencil
roots = self.roots()
# we need at least one root to not be at infinity -> non-degenerate pencil
if roots == {(np.inf + 0j): 4}:
Q2 = self.quadric(-1)
if np.linalg.det(Q2.matrix) == 0:
raise NotImplementedError(
"We only treat non-degenerate pencils here."
)
else:
# set k to a random number in case there are no real roots
k = -1
# and then check if we can choose k as smaller than the first real root
for elem in roots.keys():
if elem.imag == 0:
k = elem.real - 1
break
# calculate new quadric Q2 in the pencil at point k
Q2 = self.quadric(k)
# make sure you didn't accidentally pick the Q1 quadric of the pencil
if Q2 == self.Q1:
Q2 = self.quadric(k - 1)
return Q2
[docs] def signature_sequence(self):
"""
Calculates the signature sequence for the pencil described by ``lmda*Q2+Q1``
Uses a representation via ``Q1`` and ``Q2`` of the pencil with ``Q2`` invertible
Returns
-------
ddg.geometry.signatures.SignatureSequence
"""
# choose -Q1 instead of Q1 to match the paper
Q1 = self.quadric([-1, 0])
pencil = Pencil(Q1, self.Q2)
# make sure we have a nonsingular pair for the pencil (Q2 is regular)
Q2 = pencil.regular_matrix()
# Use equivalent pencil spanned by a nonsingular pair
pencil = Pencil(Q1, Q2)
# calculate roots of the characteristic polynomial
roots = pencil.roots()
# collect necessary lists for signature sequence creation
indices = []
signatures = []
multiplicities = []
# for every root calculate index between roots and signature in root:
for elem in roots.keys():
if elem.imag == 0 and elem.real != np.inf:
# pick any element before the first root
if indices == []:
value = elem.real - 1
# after first step: find the median between two roots
else:
value = (elem.real + prev_elem) / 2
# define new previous element to find the median in the next step
prev_elem = elem.real
# calculate index for lambda=value and add to sequence
indices += [
signature(Quadric(value * pencil.Q2.matrix + pencil.Q1.matrix)).plus
]
# calculate signature in root (prev_elem) and add to sequence
signatures += [
signature(Quadric(prev_elem * pencil.Q2.matrix + pencil.Q1.matrix))
]
# add multiplicity of root to multiplicity sequence
multiplicities += [roots[elem]]
# if there are no real roots, create indices list for minimum sequence:
if indices == []:
indices += [2]
# add last index to sequence in the end
# here we know that the first and last index sum up to ambient dimension d
else:
d = len(Q1.matrix[0])
indices += [d - indices[0]]
# create signature sequence
sign_seq = ddg.geometry.signatures.SignatureSequence(
indices, signatures, multiplicities
)
return sign_seq
[docs] def eigenvalue_curve_function(self):
"""
Calculates the eigenvalue curve function
The eigenvalue curve function of a pencil represented by ``lambda Q2 + Q1``
is defined by ``f(u, lambda) = det(lmda*Q2+Q1-uI)``,
where ``I`` is the identity matrix and ``u, lambda`` are in ``RP¹``.
This is needed for special cases in the classification of pencils
by signature sequence.
Returns
-------
function calculating the eigenvalue curve
"""
# get a nonsingular pair of matrices from the pencil
Q1, Q2 = self.quadric([1, 0]), self.regular_matrix()
# calculates the eigenvalue curve as f(u, lambda) = det(lmda*Q2+Q1-uI)
# using the lagrange version of calculating the determinant
def polynomial(u, lmda):
return ddg.math.linalg.laplace_determinant(
lmda * Q2.matrix + Q1.matrix - np.diag([u, u, u, u])
)
return polynomial
def _special_cases_for_classification(
self, discriminant=True, atol=None, rtol=None
):
"""
Treats the special cases of the signature sequence::
seq = [2] and seq = [2 ((((1, 1)))) 2]
Parameters
----------
discriminant: bool
decides whether the discriminant is calculated
- True for case ``seq = [2]``
- False for case ``seq = [2 ((((1, 1)))) 2]``
Returns
-------
list
- first entry: str (case number)
- second entry: str (explanation of cases)
"""
# define a symbol for the eigenvalue curve function
u = sympy.Symbol("u")
# define some additional symbols for homogenisation
w = sympy.Symbol("w", positive=True)
v = sympy.Symbol("v", positive=True)
# calculate the eigenvalue curve
f = self.eigenvalue_curve_function()
# evaluate the eigenvalue polynomial for u=0
# to get the characteristic function of the pencil lambda*Q2 + Q1
characteristic_funct = f(0, u)
if discriminant:
# rounding done because of comparison to 0
discr = sympy.discriminant(characteristic_funct)
discr = float(discr)
# if the discriminant !=0, we have no multiple roots
if not ddg.nonexact.isclose(discr, 0, atol, rtol):
sequenceclass = [
"type4",
"Type 4: The Quadric Surface Intersection Curve"
+ " has two non-null-homotopic components.",
]
return sequenceclass
else:
pass
else:
pass
# calculate g as the squareroot of f(0,u)
# homogenise polynomial expression (x->v/w and then *w^4)
# to be able to call for the identity matrix and take the square root
gee = sympy.factor(
(characteristic_funct.subs({u: v / w}) * w**4) ** sympy.Rational(1, 2)
)
# calculates matrix Q2^-1*Q1 to calculate g(Q2⁻¹Q1)
M = sympy.Matrix(np.linalg.inv(self.Q2.matrix)) * sympy.Matrix(self.Q1.matrix)
value = sympy.powsimp(gee.subs({v: M, w: sympy.eye(4)}))
# turn matrix into unraveled numpy array to regard entries by themselves
value = np.array(value).ravel()
# case distinction: does the squareroot vanish for Q2^(-1)Q1?
if all([ddg.nonexact.isclose(float(entry), 0, atol, rtol) for entry in value]):
if discriminant:
sequenceclass = [
"type31",
"Type 31: The Quadric Surface Intersection Curve is planar with "
+ "one pair of skew real lines and one pair of two skew imaginary "
+ "lines. The two pairs of lines intersect at two pairs of "
+ "complex conjugate points.",
]
else:
sequenceclass = [
"type35",
"Type 35: The Quadric Surface Intersection Curve "
+ "is planar with one real double line and two other "
+ "skew real lines; each of the latter two lines intersects"
+ " the real double line.",
]
else:
if discriminant:
sequenceclass = [
"type10",
"Type 10: The Quadric Surface Intersection Curve has one real line "
+ "and one space cubic, intersecting at "
+ "two complex conjugate points.",
]
else:
sequenceclass = [
"type26",
"Type 26: The Quadric Surface Intersection Curve is planar with "
+ "one real conic and two real lines; these three components "
+ "intersect at a common real point.",
]
return sequenceclass
[docs] def classification(self):
"""
Classifies the quadric intersection curve of the pencil via the type
of signature sequence
Returns
-------
list
- first entry: str (case number)
- second entry: str (explanation of cases)
"""
# calculate signature sequence of pencil (in normalform)
seq = self.signature_sequence()
if seq.dimension != 4:
raise NotImplementedError(
"Classification only implemented for quadrics in RP³"
)
seq = seq.normalform()
seq_attributes = (seq.indices, seq.signatures, seq.multiplicities)
# determine the type from the dictionary
if seq_attributes in sg.signatures_to_types:
sequenceclass = sg.signatures_to_types[seq_attributes]
# cover special cases
if sequenceclass[0] == "specialcase1":
sequenceclass = self._special_cases_for_classification(
discriminant=True
)
elif sequenceclass[0] == "specialcase2":
sequenceclass = self._special_cases_for_classification(
discriminant=False
)
else:
raise ValueError(
"Not a valid pencil for classification. "
+ "The pencil needs to be nondegenerate in RP^3."
)
return sequenceclass
[docs] def index_sequence(self):
"""
Determines the index sequence of a pencil
Derived from the classification of the pencil via signature sequence
Returns
-------
ddg.geometry.signatures.IndexSequence
in normal form
"""
# classify the pencil and read off index sequence from there
case_number = self.classification()[0]
if case_number in sg.types_to_index_and_segre:
index_seq = sg.types_to_index_and_segre[case_number][0]
else:
raise Exception("Not a valid case")
return index_seq
[docs] def segre_symbol(self):
"""
Determines the segre symbol of a pencil
Derived from the classification of the pencil via signature sequence
Returns
-------
ddg.geometry.signatures.SegreSymbol
in normal form
"""
# classify the pencil and read off index sequence from there
case_number = self.classification()[0]
if case_number in sg.types_to_index_and_segre:
segre_symb = sg.types_to_index_and_segre[case_number][1]
else:
raise Exception("Not a valid case")
return segre_symb
[docs]class QuadricIntersection(Embeddable, LinearTransformable):
"""
Intersection of two quadrics.
Parameters
----------
Q1, Q2 : ddg.geometry.Quadric
Quadrics that define the intersection
Attributes
----------
Q1, Q2 : ddg.geometry.Quadric
dimension
dimension_complex
ambient_dimension
Raises
------
ValueError
* If the objects handed are not two Quadrics
* If the geometries of the quadrics don't match
* If the quadrics are not in the same subspace
"""
def __init__(self, Q1, Q2):
if (
not isinstance(Q1, Quadric)
or not isinstance(Q2, Quadric)
or Q1.ambient_dimension != Q2.ambient_dimension
):
raise ValueError(
"Please enter two quadrics with the same ambient dimension."
)
elif Q1.subspace != Q2.subspace:
raise ValueError("Quadrics are not in the same subspace.")
else:
self.Q1 = Q1
self.Q2 = Q2
@property
def ambient_dimension(self):
return self.Q1.ambient_dimension
@property
def dimension(self):
return self.Q1.dimension + self.Q2.dimension - self.ambient_dimension
@property
def dimension_complex(self):
return (
self.Q1.dimension_complex
+ self.Q2.dimension_complex
- self.ambient_dimension
)
[docs] def pencil(self):
"""Pencil whose intersection curve is represented and that the two matrices span
Returns
-------
ddg.geometry.Pencil
"""
return Pencil(self.Q1, self.Q2)
[docs] def embed(self, subspace=None):
return QuadricIntersection(
self.Q1.embed(subspace=subspace), self.Q2.embed(subspace=subspace)
)
[docs] def unembed(self, subspace=None):
return QuadricIntersection(
self.Q1.unembed(subspace=subspace), self.Q2.unembed(subspace=subspace)
)
[docs] def at_infinity(self):
return self.Q1.at_infinity() or self.Q2.at_infinity()
def __str__(self):
return (
"intersection of the two quadrics Q1 and Q2\n"
+ f"in {self.ambient_dimension}D projective space\n"
+ "Q1:\n"
+ str(self.Q1)
+ "\nQ2:\n"
+ str(self.Q2)
)
def __repr__(self):
return (
"QuadricIntersection(\n"
+ indent(repr(self.Q1), " ")
+ ",\n "
+ indent(repr(self.Q2), " ")
+ "\n)"
)
[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 ddg.geometry.Point
Point in homogeneous coordinates. Must lie in `quadric.subspace` if
in_subspace is True.
quadric : ddg.geometry.Quadric
in_subspace : bool (default=False)
Returns
-------
ddg.geometry.Quadric or ddg.geometry.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 : ddg.geometry.Point
radius : float
absolute : ddg.geometry.Quadric
Returns
-------
ddg.geometry.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 : ddg.geometry.Point
radius : float
absolute : ddg.geometry.Quadric
Returns
-------
ddg.geometry.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 cone_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 : ddg.geometry.Quadric
Returns
-------
ddg.geometry.Subspace
Raises
------
ValueError
If Quadric does not have the signature mentioned above.
"""
cone = cone.normalize(affine=True)
ad = cone.subspace.dimension
supported_signatures = [
sg.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 quadric_normalization(quadric, affine=False):
"""Signature and normalizing transformation.
Approximately satisfies::
F.T @ quadric.matrix @ F == sgn.matrix
See :py:class:`ddg.geometry.signatures.Signature` and
:py:class:`ddg.geometry.signatures.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 : ddg.geometry.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 = sg.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 = sg.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
def intersect_quadric_subspace(quadric, subspace):
"""Intersect quadric with subspace
Parameters
----------
quadric : ddg.geometry.Quadric
subspace : ddg.geometry.Subspace
Returns
-------
ddg.geometry.Quadric
ddg.geometry.
See also
--------
ddg.geometry.intersect
To intersect more than two objects of various types.
"""
meet = subspaces.intersect_subspaces(quadric.subspace, subspace)
gram = quadric.inner_product(meet.matrix, meet.matrix)
return Quadric(gram, subspace=meet)
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 : ddg.geometry.Quadric
subspace : ddg.geometry.Subspace
Returns
-------
ddg.geometry.Quadric or ddg.geometry.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.join
To join more than two objects of various types.
"""
if quadric.subspace <= subspace:
return subspace
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_with_diagonalisable_matrices(
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**2, 1, -1])
\_______ ________/ \_______ ________/
\/ \/
D1 D2
for some k in [-1, 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 : ddg.geometry.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**2, 1, -1]) iff
#
# a1 * diag1 + a2 * diag2 = [1, 1, 0, -1]
#
# for some a1, a2 and
#
# b1 * diag1 + b2 * diag2 - k**2 * [0, 1, 0, 0] = [0, 0, 1, -1]
#
# for some b1, b2 and some k in [-1, 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**2, 1, -1]) for some -1 <= k <= 1."
)
k = np.sqrt(-c2[-1])
if not -1 - atol <= k <= 1 + atol:
raise ValueError(f"k is {k}, but should be between -1 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, -1, atol=atol, rtol=rtol):
k = -1
curve = jacobi_elliptic_curve(k)
curve = homogenize(curve)
curve = compose(lambda x: quadric1.subspace.matrix @ x, curve)
if affine:
curve = dehomogenize(curve)
return curve
[docs]def polarize(object_, quadric):
"""Polarize a geometric object with respect to a quadric.
Parameters
----------
object_ : Geometric object
Currently supported: ddg.geometry.Subspace, Quadric.
quadric : ddg.geometry.Quadric
Returns
-------
polar : type depending on the input type
Polar object
Raises
------
TypeError
If the type of object is not supported.
See also
--------
polarize_oriented
Notes
-----
This function will not ensure any orientation of the resulting polar.
"""
if isinstance(object_, subspaces.Subspace):
return _polarize_subspace(object_, quadric)
elif isinstance(object_, Quadric):
return _polarize_quadric(object_, quadric)
else:
raise TypeError(
f"Polarization for {type(object_)} objects not supported (yet)."
)
def polarize_oriented(object_, quadric):
"""Polarize a geometric object with respect to a quadric.
Parameters
----------
object_ : Geometric object
Currently supported: Subspace
quadric : ddg.geometry.Quadric
Returns
-------
polar : type depending on the input type
Polar object
Raises
------
TypeError
If the type of object is not supported.
See also
--------
polarize
Notes
-----
The polar will be (positively) oriented with respect to the basis of the subspace
containing the quadric.
"""
if isinstance(object_, subspaces.Subspace):
return _polarize_subspace_oriented(object_, quadric)
else:
raise TypeError(
f"Polarization for {type(object_)} objects not " f"supported (yet)."
)
def _polarize_subspace_helper(polarize_function, subspace, quadric):
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 = polarize_function(polar_points)
polar_points = B @ polar_points
return subspaces.subspace_from_columns(polar_points)
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 : ddg.geometry.Subspace
quadric : ddg.geometry.Quadric
Returns
-------
ddg.geometry.Subspace
Raises
------
ValueError
If subspace is not contained in quadric.subspace
"""
return _polarize_subspace_helper(ip.col_complement, subspace, quadric)
def _polarize_subspace_oriented(subspace, quadric):
"""Oriented 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}.
Note that the orientation is based on the containing subspace of the quadric
and might be negative, if the subspace itself is negatively oriented.
Parameters
----------
subspace : ddg.geometry.Subspace
quadric : ddg.geometry.Quadric
Returns
-------
ddg.geometry.Subspace
Raises
------
ValueError
If subspace is not contained in quadric.subspace
"""
return _polarize_subspace_helper(ip.col_complement_oriented, subspace, quadric)
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 : ddg.geometry.Quadric
Quadric to polarize.
quadric2 : ddg.geometry.Quadric
Quadric with respect to which to polarize.
Returns
-------
ddg.geometry.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 = quadric_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 : ddg.geometry.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.geometry.signatures.Signature,
ddg.geometry.signatures.AffineSignature
"""
if subspace is None:
return quadric_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)
[docs]def intersect_quadrics(Q1, Q2):
"""
Creates an object representing the intersection curve of two quadrics
for visualisation purposes
Parameters
----------
Q1, Q2: ddg.geometry.Quadric
Two quadrics that span a pencil whose intersection curve is represented
Returns
-------
ddg.geometry.QuadricIntersection
"""
if not isinstance(Q1, Quadric) or not isinstance(Q2, Quadric):
raise ValueError("Please enter two quadrics.")
return Pencil(Q1, Q2).intersection()
[docs]def quadric_from_points(points):
"""Returns a quadric from homogeneous coordinates of a set of
points.
There is always a quadric going through (n^2 + 3n)/2 points in
n-dimensional projective space, e.g. nine points determine a
quadric in 3d projective space.
Parameters
----------
points : iterable of ddg.geometry.Point
Returns
-------
ddg.geometry.Quadric
Raises
------
ValueError
If the number of points does not match the dimension as
stated above.
ValueError
If the given points do not determine a unique quadric.
"""
point_coords = [p.points[0] for p in points]
d = len(point_coords[0])
k = len(point_coords)
m = int((d * (d + 1)) / 2 - 1)
if k != m:
raise ValueError(
f"Number of points {k} does not match the number required "
f"{m} to determine a quadric in {d - 1}-dimensional "
"projective space."
)
symmetric_matrices = []
for i, j in combinations_with_replacement(range(d), 2):
S = np.zeros(shape=(d, d))
S[i][j] = 1
S[j][i] = 1
symmetric_matrices.append(S)
A = np.zeros(shape=(k, k + 1))
for i, coord in enumerate(point_coords):
for j, S in enumerate(symmetric_matrices):
A[i][j] = coord.T @ S @ coord
kern = nullspace(A)
if kern.shape[1] != 1:
raise ValueError("The given points do not determine a unique quadric.")
quadric_matrix = np.zeros(shape=(d, d))
for a, b in zip(kern, symmetric_matrices):
quadric_matrix += a * b
return Quadric(quadric_matrix)
[docs]def quadric_from_three_skew_lines(lines):
"""Returns a quadric in from three skew lines in 3-dimensional
projective space.
Given three skew lines in 3-dimensional projective space, the
family of lines that meet all three lines sweeps out a unique
quadric.
Parameters
----------
lines : iterable of ddg.geometry.Subspace
Three skew lines in 3-dimensional projective space.
Returns
-------
ddg.geometry.Quadric
Raises
------
ValueError
If length of `lines` iterable is not exactly three.
Notes
-----
Skew lines are lines in 3-dimensional projective space that
are neither parallel nor intersecting.
"""
k = len(lines)
if k != 3:
raise ValueError(
f"Three skew lines determine a quadric in "
"3-dimensional projective space. Number of lines "
f"given is {k}."
)
line1, line2, line3 = lines
# Raise exception if given lines are not pairwise skew
for i, j in combinations(range(3), 2):
intersection = subspaces.intersect_subspaces(lines[i], lines[j])
if intersection.points != []:
if intersection.at_infinity():
raise ValueError(
"Given three lines are not pairwise skew. "
f"Line {i+1} and Line {j+1} are parallel."
)
else:
raise ValueError(
"Given three lines are not pairwise skew. "
f"Line {i+1} and Line {j+1} intersect."
)
# Generate sample points on the first line
base_point = line1.affine_points[0]
direction_vector = line1.affine_points[0] - line1.affine_points[1]
sample_points = [
base_point,
base_point + 2 * direction_vector,
base_point + 3 * direction_vector,
]
points = []
for point in sample_points:
# Get base point and direction vector on second line
q = line2.affine_points[0]
v = line2.affine_points[1] - line2.affine_points[0]
point_subspace = subspaces.Point(pm.homogenize(point))
# Generate plane from normal, intersect it with third line and take
# the span on intersection and given sample point.
plane = subspaces.hyperplane_from_normal(normal=np.cross(point - q, v), point=q)
intersection_point = subspaces.intersect_subspaces(plane, line3)
line = subspaces.join_subspaces(intersection_point, point_subspace)
# Append three points from a given sample point
points.append(subspaces.intersect_subspaces(line2, line))
points.append(subspaces.intersect_subspaces(line3, line))
points.append(point_subspace)
return quadric_from_points(points)