import math
import numpy as np
import numpy.linalg as la
import ddg.math.linalg as linalg
[docs]def get_col_complement(A, inner_product=None, ensure_pos_det=True):
"""
Returns the orthogonal complement of the columns of a given
array with respect to the given inner product.
Parameters
----------
A : numpy.array
Array for which to return the column complement
inner_product : Callable (optional, default=None)
Inner product to be used. When the argument is None, the standard
inner product will be used.
ensure_pos_det : bool (optional, default=True)
When set to True, the function will make sure that the resulting basis
will be positively oriented.
Returns
-------
numpy.array
Column complement of the given array
"""
if inner_product is None:
inner_product = inner_product_from_matrix(np.eye(np.shape(A)[0]))
B = linalg.col_basis(A)
# Return whole space. Avoids index error
if (B == 0).all(): # 0 matrix or empty matrix
return np.eye(np.shape(B)[0])
AL = np.dot(get_matrix(inner_product, len(B[:, 0])), B)
Q, R = la.qr(AL, mode="complete")
ACN = orthonormalize(Q[:, len(B[0]) :], inner_product)
mat = np.concatenate((A, ACN), axis=1)
if (
ACN.shape[1] != 0
and mat.shape[0] == mat.shape[1]
and ensure_pos_det
and np.linalg.slogdet(mat)[0] < 0
):
ACN.T[-1] = -ACN.T[-1]
return ACN # A complement normalized
[docs]def get_row_complement(A, inner_product=None, ensure_pos_det=True):
"""
Returns the orthogonal complement of the rows of a given
array with respect to the given inner product.
Parameters
----------
A : numpy.array
Array for which to return the column complement
inner_product : Callable (optional, default=None)
Inner product to be used. When the argument is None, the standard
inner product will be used.
ensure_pos_det : bool (optional, default=True)
When set to True, the function will make sure that the resulting basis
will be positively oriented.
Returns
-------
numpy.array
Column complement of the given array
"""
return get_col_complement(A.T, inner_product, ensure_pos_det).T
[docs]def orthonormalize(U, inner_product, tol=1e-6):
"""
Turn U into an orthonormal form with respect to the inner_product
:param U:
:param inner_product:
:return:
"""
UN = U
for i in range(len(U[0])):
u = U[:, i]
# print(Projective.inner_product(u, u))
uu = inner_product(u, u)
if math.fabs(uu) < tol:
for j in range(i + 1, len(U[0])):
v = U[:, j]
vv = inner_product(v, v)
if math.fabs(vv) >= tol:
U[:, i] = v
U[:, j] = u
break
elif math.fabs(inner_product(u, v)) > 1e-6:
U[:, i] = u + v
U[:, j] = u - v
UN[:, i] = normalize(u, inner_product)
# print(str(i) + ":" + str(UN[:, i]))
for j in range(i + 1, len(U[0])):
UN[:, j] = (
UN[:, j]
- inner_product(UN[:, j], UN[:, i])
/ inner_product(UN[:, i], UN[:, i])
* UN[:, i]
)
return UN
[docs]def normalize(v, inner_product):
vv = inner_product(v, v)
if vv != 0:
return v / math.sqrt(math.fabs(vv))
else:
return v
[docs]def gram_matrix(inner_product, matrix_of_vectors):
"""
The vectors are supposed to be the columns of the matrix
"""
n = len(matrix_of_vectors[0])
g = np.empty([n, n])
for i, v in enumerate(matrix_of_vectors.T):
for j, w in enumerate(matrix_of_vectors.T):
g[i, j] = inner_product(v, w)
return g
[docs]def signature(inner_product, matrix):
g = gram_matrix(inner_product, matrix)
s, v = la.eig(g)
signature = np.empty(len(s))
for i, l in enumerate(s):
if l == 0:
signature[i] = 0
else:
signature[i] = l / math.fabs(l)
return np.sort(signature)
[docs]def reflect(normal, source, inner_product):
n2 = inner_product(normal, normal)
return source - 2 * inner_product(source, normal) / n2 * normal
[docs]def get_matrix(inner_product, n):
eye = np.eye(n)
return gram_matrix(inner_product, eye)
[docs]def inner_product_from_matrix(matrix):
def inner_product(v, w):
return np.dot(np.dot(matrix, v), w)
return inner_product
[docs]def project_onto_complement(normal, source, inner_product):
n2 = inner_product(normal, normal)
return source - inner_product(source, normal) / n2 * normal
[docs]def polarity_matrix(inner_product, n):
return get_matrix(inner_product, n).T
[docs]def points_on_quadric(subspace, inner_product):
return points_on_quadric_from_onb(orthonormalize(subspace.T, inner_product).T)
[docs]def points_on_quadric_from_onb(onb, inner_product):
"""
Compute possible points on the quadric from an orthonormal basis by suitable
linear combinations of space- and time-like vectors. Light-like vectors will
be appended to the list.
Parameters
----------
onb : numpy.array
rows of the matrix form an orthonormal basis
Returns
-------
numpy.array
sum and difference of the space- and time-like vectors, and the
light-like vectors
"""
spacelike = np.array([c for c in onb if inner_product(c, c) > 0])
timelike = np.array([c for c in onb if inner_product(c, c) < 0])
lightlike = np.array([c for c in onb if inner_product(c, c) == 0])
s = np.empty((len(spacelike) * len(timelike) * 2 + len(lightlike), len(onb[0])))
i = 0
for space in spacelike:
for time in timelike:
s[i, :] = space + time
i = i + 1
s[i, :] = space - time
i = i + 1
for light in lightlike:
s[i, :] = light
i = i + 1
return s