X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_algebra.py;h=79187594472d82e24a0b700305b7ddfc7b192536;hb=b259821a73cb75a6d5a81d13256802be023b0fa9;hp=8b37a83602ef9b00b5a14c55785c2163b44df32b;hpb=0cf1bd4fb459733e559f2040089e1905fc6af9ca;p=sage.d.git diff --git a/mjo/eja/eja_algebra.py b/mjo/eja/eja_algebra.py index 8b37a83..7918759 100644 --- a/mjo/eja/eja_algebra.py +++ b/mjo/eja/eja_algebra.py @@ -170,6 +170,17 @@ from mjo.eja.eja_element import FiniteDimensionalEJAElement from mjo.eja.eja_operator import FiniteDimensionalEJAOperator from mjo.eja.eja_utils import _all2list, _mat2vec +def EuclideanJordanAlgebras(field): + r""" + The category of Euclidean Jordan algebras over ``field``, which + must be a subfield of the real numbers. For now this is just a + convenient wrapper around all of the other category axioms that + apply to all EJAs. + """ + category = MagmaticAlgebras(field).FiniteDimensional() + category = category.WithBasis().Unital().Commutative() + return category + class FiniteDimensionalEJA(CombinatorialFreeModule): r""" A finite-dimensional Euclidean Jordan algebra. @@ -228,6 +239,26 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): """ Element = FiniteDimensionalEJAElement + @staticmethod + def _check_input_field(field): + if not field.is_subring(RR): + # Note: this does return true for the real algebraic + # field, the rationals, and any quadratic field where + # we've specified a real embedding. + raise ValueError("scalar field is not real") + + @staticmethod + def _check_input_axioms(basis, jordan_product, inner_product): + if not all( jordan_product(bi,bj) == jordan_product(bj,bi) + for bi in basis + for bj in basis ): + raise ValueError("Jordan product is not commutative") + + if not all( inner_product(bi,bj) == inner_product(bj,bi) + for bi in basis + for bj in basis ): + raise ValueError("inner-product is not commutative") + def __init__(self, basis, jordan_product, @@ -244,30 +275,14 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): n = len(basis) if check_field: - if not field.is_subring(RR): - # Note: this does return true for the real algebraic - # field, the rationals, and any quadratic field where - # we've specified a real embedding. - raise ValueError("scalar field is not real") + self._check_input_field(field) if check_axioms: # Check commutativity of the Jordan and inner-products. # This has to be done before we build the multiplication # and inner-product tables/matrices, because we take # advantage of symmetry in the process. - if not all( jordan_product(bi,bj) == jordan_product(bj,bi) - for bi in basis - for bj in basis ): - raise ValueError("Jordan product is not commutative") - - if not all( inner_product(bi,bj) == inner_product(bj,bi) - for bi in basis - for bj in basis ): - raise ValueError("inner-product is not commutative") - - - category = MagmaticAlgebras(field).FiniteDimensional() - category = category.WithBasis().Unital().Commutative() + self._check_input_axioms(basis, jordan_product, inner_product) if n <= 1: # All zero- and one-dimensional algebras are just the real @@ -286,6 +301,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): for bj in basis for bk in basis) + category = EuclideanJordanAlgebras(field) + if associative: # Element subalgebras can take advantage of this. category = category.Associative() @@ -309,7 +326,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # as well as a subspace W of V spanned by those (vectorized) # basis elements. The W-coordinates are the coefficients that # we see in things like x = 1*b1 + 2*b2. - vector_basis = basis degree = 0 if n > 0: @@ -319,9 +335,11 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # written out as "long vectors." V = VectorSpace(field, degree) - # The matrix that will hole the orthonormal -> unorthonormal - # coordinate transformation. - self._deortho_matrix = None + # The matrix that will hold the orthonormal -> unorthonormal + # coordinate transformation. Default to an identity matrix of + # the appropriate size to avoid special cases for None + # everywhere. + self._deortho_matrix = matrix.identity(field,n) if orthonormalize: # Save a copy of the un-orthonormalized basis for later. @@ -346,17 +364,22 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # its own set of non-ambient coordinates (in terms of the # supplied basis). vector_basis = tuple( V(_all2list(b)) for b in basis ) - W = V.span_of_basis( vector_basis, check=check_axioms) + + # Save the span of our matrix basis (when written out as long + # vectors) because otherwise we'll have to reconstruct it + # every time we want to coerce a matrix into the algebra. + self._matrix_span = V.span_of_basis( vector_basis, check=check_axioms) if orthonormalize: - # Now "W" is the vector space of our algebra coordinates. The - # variables "X1", "X2",... refer to the entries of vectors in - # W. Thus to convert back and forth between the orthonormal - # coordinates and the given ones, we need to stick the original - # basis in W. + # Now "self._matrix_span" is the vector space of our + # algebra coordinates. The variables "X1", "X2",... refer + # to the entries of vectors in self._matrix_span. Thus to + # convert back and forth between the orthonormal + # coordinates and the given ones, we need to stick the + # original basis in self._matrix_span. U = V.span_of_basis( deortho_vector_basis, check=check_axioms) - self._deortho_matrix = matrix( U.coordinate_vector(q) - for q in vector_basis ) + self._deortho_matrix = matrix.column( U.coordinate_vector(q) + for q in vector_basis ) # Now we actually compute the multiplication and inner-product @@ -377,7 +400,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # The jordan product returns a matrixy answer, so we # have to convert it to the algebra coordinates. elt = jordan_product(q_i, q_j) - elt = W.coordinate_vector(V(_all2list(elt))) + elt = self._matrix_span.coordinate_vector(V(_all2list(elt))) self._multiplication_table[i][j] = self.from_vector(elt) if not orthonormalize: @@ -684,8 +707,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): def _element_constructor_(self, elt): """ - Construct an element of this algebra from its vector or matrix - representation. + Construct an element of this algebra or a subalgebra from its + EJA element, vector, or matrix representation. This gets called only after the parent element _call_ method fails to find a coercion for the argument. @@ -724,6 +747,16 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) ) b1 + b5 + Subalgebra elements are embedded into the superalgebra:: + + sage: J = JordanSpinEJA(3) + sage: J.one() + b0 + sage: x = sum(J.gens()) + sage: A = x.subalgebra_generated_by() + sage: J(A.one()) + b0 + TESTS: Ensure that we can convert any element back and forth @@ -748,6 +781,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): Traceback (most recent call last): ... ValueError: not an element of this algebra + """ msg = "not an element of this algebra" if elt in self.base_ring(): @@ -757,13 +791,16 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # that the integer 3 belongs to the space of 2-by-2 matrices. raise ValueError(msg) - try: - # Try to convert a vector into a column-matrix... + if hasattr(elt, 'superalgebra_element'): + # Handle subalgebra elements + if elt.parent().superalgebra() == self: + return elt.superalgebra_element() + + if hasattr(elt, 'sparse_vector'): + # Convert a vector into a column-matrix. We check for + # "sparse_vector" and not "column" because matrices also + # have a "column" method. elt = elt.column() - except (AttributeError, TypeError): - # and ignore failure, because we weren't really expecting - # a vector as an argument anyway. - pass if elt not in self.matrix_space(): raise ValueError(msg) @@ -780,15 +817,10 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # is that we're already converting everything to long vectors, # and that strategy works for tuples as well. # - # We pass check=False because the matrix basis is "guaranteed" - # to be linearly independent... right? Ha ha. - elt = _all2list(elt) - V = VectorSpace(self.base_ring(), len(elt)) - W = V.span_of_basis( (V(_all2list(s)) for s in self.matrix_basis()), - check=False) + elt = self._matrix_span.ambient_vector_space()(_all2list(elt)) try: - coords = W.coordinate_vector(V(elt)) + coords = self._matrix_span.coordinate_vector(elt) except ArithmeticError: # vector is not in free module raise ValueError(msg) @@ -1394,7 +1426,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule): # corresponding to trivial spaces (e.g. it returns only the # eigenspace corresponding to lambda=1 if you take the # decomposition relative to the identity element). - trivial = self.subalgebra(()) + trivial = self.subalgebra((), check_axioms=False) J0 = trivial # eigenvalue zero J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half J1 = trivial # eigenvalue one @@ -1753,13 +1785,6 @@ class RationalBasisEJA(FiniteDimensionalEJA): a = ( a_i.change_ring(self.base_ring()) for a_i in self._rational_algebra._charpoly_coefficients() ) - if self._deortho_matrix is None: - # This can happen if our base ring was, say, AA and we - # chose not to (or didn't need to) orthonormalize. It's - # still faster to do the computations over QQ even if - # the numbers in the boxes stay the same. - return tuple(a) - # Otherwise, convert the coordinate variables back to the # deorthonormalized ones. R = self.coordinate_polynomial_ring() @@ -1833,11 +1858,36 @@ class ConcreteEJA(FiniteDimensionalEJA): def random_instance(cls, max_dimension=None, *args, **kwargs): """ Return a random instance of this type of algebra whose dimension - is less than or equal to ``max_dimension``. If the dimension bound - is omitted, then the ``_max_random_instance_dimension()`` is used - to get a suitable bound. + is less than or equal to the lesser of ``max_dimension`` and + the value returned by ``_max_random_instance_dimension()``. If + the dimension bound is omitted, then only the + ``_max_random_instance_dimension()`` is used as a bound. This method should be implemented in each subclass. + + SETUP:: + + sage: from mjo.eja.eja_algebra import ConcreteEJA + + TESTS: + + Both the class bound and the ``max_dimension`` argument are upper + bounds on the dimension of the algebra returned:: + + sage: from sage.misc.prandom import choice + sage: eja_class = choice(ConcreteEJA.__subclasses__()) + sage: class_max_d = eja_class._max_random_instance_dimension() + sage: J = eja_class.random_instance(max_dimension=20, + ....: field=QQ, + ....: orthonormalize=False) + sage: J.dimension() <= class_max_d + True + sage: J = eja_class.random_instance(max_dimension=2, + ....: field=QQ, + ....: orthonormalize=False) + sage: J.dimension() <= 2 + True + """ from sage.misc.prandom import choice eja_class = choice(cls.__subclasses__()) @@ -1845,7 +1895,7 @@ class ConcreteEJA(FiniteDimensionalEJA): # These all bubble up to the RationalBasisEJA superclass # constructor, so any (kw)args valid there are also valid # here. - return eja_class.random_instance(*args, **kwargs) + return eja_class.random_instance(max_dimension, *args, **kwargs) class MatrixEJA(FiniteDimensionalEJA): @@ -2077,14 +2127,15 @@ class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/2 - 1/2)) @classmethod - def random_instance(cls, max_dimension=None, **kwargs): + def random_instance(cls, max_dimension=None, *args, **kwargs): """ Return a random instance of this type of algebra. """ - if max_dimension is None: - max_dimension = cls._max_random_instance_dimension() - max_size = cls._max_random_instance_size(max_dimension) + 1 - n = ZZ.random_element(max_size) + class_max_d = cls._max_random_instance_dimension() + if (max_dimension is None or max_dimension > class_max_d): + max_dimension = class_max_d + max_size = cls._max_random_instance_size(max_dimension) + n = ZZ.random_element(max_size + 1) return cls(n, **kwargs) def __init__(self, n, field=AA, **kwargs): @@ -2201,13 +2252,15 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): return ZZ(int(ZZ(max_dimension).sqrt())) @classmethod - def random_instance(cls, max_dimension=None, **kwargs): + def random_instance(cls, max_dimension=None, *args, **kwargs): """ Return a random instance of this type of algebra. """ - if max_dimension is None: - max_dimension = cls._max_random_instance_dimension() - n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1) + class_max_d = cls._max_random_instance_dimension() + if (max_dimension is None or max_dimension > class_max_d): + max_dimension = class_max_d + max_size = cls._max_random_instance_size(max_dimension) + n = ZZ.random_element(max_size + 1) return cls(n, **kwargs) @@ -2297,13 +2350,15 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/4 + 1/4)) @classmethod - def random_instance(cls, max_dimension=None, **kwargs): + def random_instance(cls, max_dimension=None, *args, **kwargs): """ Return a random instance of this type of algebra. """ - if max_dimension is None: - max_dimension = cls._max_random_instance_dimension() - n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1) + class_max_d = cls._max_random_instance_dimension() + if (max_dimension is None or max_dimension > class_max_d): + max_dimension = class_max_d + max_size = cls._max_random_instance_size(max_dimension) + n = ZZ.random_element(max_size + 1) return cls(n, **kwargs) class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): @@ -2410,13 +2465,15 @@ class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA): return 0 @classmethod - def random_instance(cls, max_dimension=None, **kwargs): + def random_instance(cls, max_dimension=None, *args, **kwargs): """ Return a random instance of this type of algebra. """ - if max_dimension is None: - max_dimension = cls._max_random_instance_dimension() - n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1) + class_max_d = cls._max_random_instance_dimension() + if (max_dimension is None or max_dimension > class_max_d): + max_dimension = class_max_d + max_size = cls._max_random_instance_size(max_dimension) + n = ZZ.random_element(max_size + 1) return cls(n, **kwargs) def __init__(self, n, field=AA, **kwargs): @@ -2552,13 +2609,15 @@ class HadamardEJA(RationalBasisEJA, ConcreteEJA): return max_dimension @classmethod - def random_instance(cls, max_dimension=None, **kwargs): + def random_instance(cls, max_dimension=None, *args, **kwargs): """ Return a random instance of this type of algebra. """ - if max_dimension is None: - max_dimension = cls._max_random_instance_dimension() - n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1) + class_max_d = cls._max_random_instance_dimension() + if (max_dimension is None or max_dimension > class_max_d): + max_dimension = class_max_d + max_size = cls._max_random_instance_size(max_dimension) + n = ZZ.random_element(max_size + 1) return cls(n, **kwargs) @@ -2713,13 +2772,16 @@ class BilinearFormEJA(RationalBasisEJA, ConcreteEJA): return max_dimension @classmethod - def random_instance(cls, max_dimension=None, **kwargs): + def random_instance(cls, max_dimension=None, *args, **kwargs): """ Return a random instance of this algebra. """ - if max_dimension is None: - max_dimension = cls._max_random_instance_dimension() - n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1) + class_max_d = cls._max_random_instance_dimension() + if (max_dimension is None or max_dimension > class_max_d): + max_dimension = class_max_d + max_size = cls._max_random_instance_size(max_dimension) + n = ZZ.random_element(max_size + 1) + if n.is_zero(): B = matrix.identity(ZZ, n) return cls(B, **kwargs) @@ -2730,6 +2792,7 @@ class BilinearFormEJA(RationalBasisEJA, ConcreteEJA): alpha = ZZ.zero() while alpha.is_zero(): alpha = ZZ.random_element().abs() + B22 = M.transpose()*M + alpha*I from sage.matrix.special import block_matrix @@ -2803,15 +2866,17 @@ class JordanSpinEJA(BilinearFormEJA): super().__init__(B, *args, **kwargs) @classmethod - def random_instance(cls, max_dimension=None, **kwargs): + def random_instance(cls, max_dimension=None, *args, **kwargs): """ Return a random instance of this type of algebra. Needed here to override the implementation for ``BilinearFormEJA``. """ - if max_dimension is None: - max_dimension = cls._max_random_instance_dimension() - n = ZZ.random_element(cls._max_random_instance_size(max_dimension) + 1) + class_max_d = cls._max_random_instance_dimension() + if (max_dimension is None or max_dimension > class_max_d): + max_dimension = class_max_d + max_size = cls._max_random_instance_size(max_dimension) + n = ZZ.random_element(max_size + 1) return cls(n, **kwargs) @@ -2868,9 +2933,12 @@ class TrivialEJA(RationalBasisEJA, ConcreteEJA): self.one.set_cache( self.zero() ) @classmethod - def random_instance(cls, **kwargs): + def random_instance(cls, max_dimension=None, *args, **kwargs): # We don't take a "size" argument so the superclass method is - # inappropriate for us. + # inappropriate for us. The ``max_dimension`` argument is + # included so that if this method is called generically with a + # ``max_dimension=`` argument, we don't try to pass + # it on to the algebra constructor. return cls(**kwargs) @@ -3086,6 +3154,13 @@ class CartesianProductEJA(FiniteDimensionalEJA): check_field=False, check_axioms=False) + # Since we don't (re)orthonormalize the basis, the FDEJA + # constructor is going to set self._deortho_matrix to the + # identity matrix. Here we set it to the correct value using + # the deortho matrices from our factors. + self._deortho_matrix = matrix.block_diagonal( [J._deortho_matrix + for J in factors] ) + self.rank.set_cache(sum(J.rank() for J in factors)) ones = tuple(J.one().to_matrix() for J in factors) self.one.set_cache(self(ones)) @@ -3416,6 +3491,21 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA, RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA def random_eja(max_dimension=None, *args, **kwargs): + r""" + + SETUP:: + + sage: from mjo.eja.eja_algebra import random_eja + + TESTS:: + + sage: set_random_seed() + sage: n = ZZ.random_element(1,5) + sage: J = random_eja(max_dimension=n, field=QQ, orthonormalize=False) + sage: J.dimension() <= n + True + + """ # Use the ConcreteEJA default as the total upper bound (regardless # of any whether or not any individual factors set a lower limit). if max_dimension is None: