]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: fix the natural representation in trivial subalgebras.
[sage.d.git] / mjo / eja / eja_algebra.py
index 5e9c07c01be4bb5ac13644c5335a8b45304c6b7e..42b2474d6971d81c8a1fcda67b079892a3e8af25 100644 (file)
@@ -5,25 +5,33 @@ are used in optimization, and have some additional nice methods beyond
 what can be supported in a general Jordan Algebra.
 """
 
-#from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra import FiniteDimensionalAlgebra
 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
-from sage.categories.finite_dimensional_algebras_with_basis import FiniteDimensionalAlgebrasWithBasis
+from sage.categories.magmatic_algebras import MagmaticAlgebras
 from sage.combinat.free_module import CombinatorialFreeModule
 from sage.matrix.constructor import matrix
+from sage.matrix.matrix_space import MatrixSpace
 from sage.misc.cachefunc import cached_method
 from sage.misc.prandom import choice
-from sage.modules.free_module import VectorSpace
+from sage.misc.table import table
+from sage.modules.free_module import FreeModule, VectorSpace
 from sage.rings.integer_ring import ZZ
 from sage.rings.number_field.number_field import QuadraticField
 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
 from sage.rings.rational_field import QQ
 from sage.structure.element import is_Matrix
-from sage.structure.category_object import normalize_names
 
 from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement
 from mjo.eja.eja_utils import _mat2vec
 
 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
+    # This is an ugly hack needed to prevent the category framework
+    # from implementing a coercion from our base ring (e.g. the
+    # rationals) into the algebra. First of all -- such a coercion is
+    # nonsense to begin with. But more importantly, it tries to do so
+    # in the category of rings, and since our algebras aren't
+    # associative they generally won't be rings.
+    _no_generic_basering_coercion = True
+
     def __init__(self,
                  field,
                  mult_table,
@@ -50,9 +58,11 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         """
         self._rank = rank
         self._natural_basis = natural_basis
-        self._multiplication_table = mult_table
+
         if category is None:
-            category = FiniteDimensionalAlgebrasWithBasis(field).Unital()
+            category = MagmaticAlgebras(field).FiniteDimensional()
+            category = category.WithBasis().Unital()
+
         fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
         fda.__init__(field,
                      range(len(mult_table)),
@@ -60,6 +70,82 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
                      category=category)
         self.print_options(bracket='')
 
+        # The multiplication table we're given is necessarily in terms
+        # of vectors, because we don't have an algebra yet for
+        # anything to be an element of. However, it's faster in the
+        # long run to have the multiplication table be in terms of
+        # algebra elements. We do this after calling the superclass
+        # constructor so that from_vector() knows what to do.
+        self._multiplication_table = [ map(lambda x: self.from_vector(x), ls)
+                                       for ls in mult_table ]
+
+
+    def _element_constructor_(self, elt):
+        """
+        Construct an element of this algebra from its natural
+        representation.
+
+        This gets called only after the parent element _call_ method
+        fails to find a coercion for the argument.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            ....:                                  RealCartesianProductEJA,
+            ....:                                  RealSymmetricEJA)
+
+        EXAMPLES:
+
+        The identity in `S^n` is converted to the identity in the EJA::
+
+            sage: J = RealSymmetricEJA(3)
+            sage: I = matrix.identity(QQ,3)
+            sage: J(I) == J.one()
+            True
+
+        This skew-symmetric matrix can't be represented in the EJA::
+
+            sage: J = RealSymmetricEJA(3)
+            sage: A = matrix(QQ,3, lambda i,j: i-j)
+            sage: J(A)
+            Traceback (most recent call last):
+            ...
+            ArithmeticError: vector is not in free module
+
+        TESTS:
+
+        Ensure that we can convert any element of the two non-matrix
+        simple algebras (whose natural representations are their usual
+        vector representations) back and forth faithfully::
+
+            sage: set_random_seed()
+            sage: J = RealCartesianProductEJA(5)
+            sage: x = J.random_element()
+            sage: J(x.to_vector().column()) == x
+            True
+            sage: J = JordanSpinEJA(5)
+            sage: x = J.random_element()
+            sage: J(x.to_vector().column()) == x
+            True
+
+        """
+        if elt == 0:
+            # The superclass implementation of random_element()
+            # needs to be able to coerce "0" into the algebra.
+            return self.zero()
+
+        natural_basis = self.natural_basis()
+        if elt not in natural_basis[0].matrix_space():
+            raise ValueError("not a naturally-represented algebra element")
+
+        # Thanks for nothing! Matrix spaces aren't vector
+        # spaces in Sage, so we have to figure out its
+        # natural-basis coordinates ourselves.
+        V = VectorSpace(elt.base_ring(), elt.nrows()*elt.ncols())
+        W = V.span_of_basis( _mat2vec(s) for s in natural_basis )
+        coords =  W.coordinate_vector(_mat2vec(elt))
+        return self.from_vector(coords)
+
 
     def _repr_(self):
         """
@@ -74,20 +160,16 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         Ensure that it says what we think it says::
 
             sage: JordanSpinEJA(2, field=QQ)
-            Euclidean Jordan algebra of degree 2 over Rational Field
+            Euclidean Jordan algebra of dimension 2 over Rational Field
             sage: JordanSpinEJA(3, field=RDF)
-            Euclidean Jordan algebra of degree 3 over Real Double Field
+            Euclidean Jordan algebra of dimension 3 over Real Double Field
 
         """
-        # TODO: change this to say "dimension" and fix all the tests.
-        fmt = "Euclidean Jordan algebra of degree {} over {}"
+        fmt = "Euclidean Jordan algebra of dimension {} over {}"
         return fmt.format(self.dimension(), self.base_ring())
 
     def product_on_basis(self, i, j):
-        ei = self.basis()[i]
-        ej = self.basis()[j]
-        Lei = self._multiplication_table[i]
-        return self.from_vector(Lei*ej.to_vector())
+        return self._multiplication_table[i][j]
 
     def _a_regular_element(self):
         """
@@ -126,8 +208,13 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         determinant).
         """
         z = self._a_regular_element()
-        V = self.vector_space()
-        V1 = V.span_of_basis( (z**k).to_vector() for k in range(self.rank()) )
+        # Don't use the parent vector space directly here in case this
+        # happens to be a subalgebra. In that case, we would be e.g.
+        # two-dimensional but span_of_basis() would expect three
+        # coordinates.
+        V = VectorSpace(self.base_ring(), self.vector_space().dimension())
+        basis = [ (z**k).to_vector() for k in range(self.rank()) ]
+        V1 = V.span_of_basis( basis )
         b =  (V1.basis() + V1.complement().basis())
         return V.span_of_basis(b)
 
@@ -178,15 +265,40 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         r = self.rank()
         n = self.dimension()
 
-        # Construct a new algebra over a multivariate polynomial ring...
+        # Turn my vector space into a module so that "vectors" can
+        # have multivatiate polynomial entries.
         names = tuple('X' + str(i) for i in range(1,n+1))
         R = PolynomialRing(self.base_ring(), names)
-        J = FiniteDimensionalEuclideanJordanAlgebra(
-              R,
-              tuple(self._multiplication_table),
-              r)
 
-        idmat = matrix.identity(J.base_ring(), n)
+        # Using change_ring() on the parent's vector space doesn't work
+        # here because, in a subalgebra, that vector space has a basis
+        # and change_ring() tries to bring the basis along with it. And
+        # that doesn't work unless the new ring is a PID, which it usually
+        # won't be.
+        V = FreeModule(R,n)
+
+        # Now let x = (X1,X2,...,Xn) be the vector whose entries are
+        # indeterminates...
+        x = V(names)
+
+        # And figure out the "left multiplication by x" matrix in
+        # that setting.
+        lmbx_cols = []
+        monomial_matrices = [ self.monomial(i).operator().matrix()
+                              for i in range(n) ] # don't recompute these!
+        for k in range(n):
+            ek = self.monomial(k).to_vector()
+            lmbx_cols.append(
+              sum( x[i]*(monomial_matrices[i]*ek)
+                   for i in range(n) ) )
+        Lx = matrix.column(R, lmbx_cols)
+
+        # Now we can compute powers of x "symbolically"
+        x_powers = [self.one().to_vector(), x]
+        for d in range(2, r+1):
+            x_powers.append( Lx*(x_powers[-1]) )
+
+        idmat = matrix.identity(R, n)
 
         W = self._charpoly_basis_space()
         W = W.change_ring(R.fraction_field())
@@ -206,18 +318,10 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         # We want the middle equivalent thing in our matrix, but use
         # the first equivalent thing instead so that we can pass in
         # standard coordinates.
-        x = J(W(R.gens()))
-
-        # Handle the zeroth power separately, because computing
-        # the unit element in J is mathematically suspect.
-        x0 = W.coordinate_vector(self.one().to_vector())
-        l1  = [ x0.column() ]
-        l1 += [ W.coordinate_vector((x**k).to_vector()).column()
-                for k in range(1,r) ]
-        l2 = [idmat.column(k-1).column() for k in range(r+1, n+1)]
-        A_of_x = matrix.block(R, 1, n, (l1 + l2))
-        xr = W.coordinate_vector((x**r).to_vector())
-        return (A_of_x, x, xr, A_of_x.det())
+        x_powers = [ W.coordinate_vector(xp) for xp in x_powers ]
+        l2 = [idmat.column(k-1) for k in range(r+1, n+1)]
+        A_of_x = matrix.column(R, n, (x_powers[:r] + l2))
+        return (A_of_x, x, x_powers[r], A_of_x.det())
 
 
     @cached_method
@@ -312,6 +416,63 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return x.trace_inner_product(y)
 
 
+    def is_trivial(self):
+        """
+        Return whether or not this algebra is trivial.
+
+        A trivial algebra contains only the zero element.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
+
+        EXAMPLES::
+
+            sage: J = ComplexHermitianEJA(3)
+            sage: J.is_trivial()
+            False
+            sage: A = J.zero().subalgebra_generated_by()
+            sage: A.is_trivial()
+            True
+
+        """
+        return self.dimension() == 0
+
+
+    def multiplication_table(self):
+        """
+        Return a visual representation of this algebra's multiplication
+        table (on basis elements).
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import JordanSpinEJA
+
+        EXAMPLES::
+
+            sage: J = JordanSpinEJA(4)
+            sage: J.multiplication_table()
+            +----++----+----+----+----+
+            | *  || e0 | e1 | e2 | e3 |
+            +====++====+====+====+====+
+            | e0 || e0 | e1 | e2 | e3 |
+            +----++----+----+----+----+
+            | e1 || e1 | e0 | 0  | 0  |
+            +----++----+----+----+----+
+            | e2 || e2 | 0  | e0 | 0  |
+            +----++----+----+----+----+
+            | e3 || e3 | 0  | 0  | e0 |
+            +----++----+----+----+----+
+
+        """
+        M = list(self._multiplication_table) # copy
+        for i in range(len(M)):
+            # M had better be "square"
+            M[i] = [self.monomial(i)] + M[i]
+        M = [["*"] + list(self.gens())] + M
+        return table(M, header_row=True, header_column=True, frame=True)
+
+
     def natural_basis(self):
         """
         Return a more-natural representation of this algebra's basis.
@@ -354,11 +515,23 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         """
         if self._natural_basis is None:
-            return tuple( b.to_vector().column() for b in self.basis() )
+            M = self.natural_basis_space()
+            return tuple( M(b.to_vector()) for b in self.basis() )
         else:
             return self._natural_basis
 
 
+    def natural_basis_space(self):
+        """
+        Return the matrix space in which this algebra's natural basis
+        elements live.
+        """
+        if self._natural_basis is None or len(self._natural_basis) == 0:
+            return MatrixSpace(self.base_ring(), self.dimension(), 1)
+        else:
+            return self._natural_basis[0].matrix_space()
+
+
     @cached_method
     def one(self):
         """
@@ -375,7 +548,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: J.one()
             e0 + e1 + e2 + e3 + e4
 
-        TESTS::
+        TESTS:
 
         The identity element acts like the identity::
 
@@ -420,6 +593,15 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return self.linear_combination(zip(self.gens(), coeffs))
 
 
+    def random_element(self):
+        # Temporary workaround for https://trac.sagemath.org/ticket/28327
+        if self.is_trivial():
+            return self.zero()
+        else:
+            s = super(FiniteDimensionalEuclideanJordanAlgebra, self)
+            return s.random_element()
+
+
     def rank(self):
         """
         Return the rank of this EJA.
@@ -533,18 +715,21 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra):
         sage: e2*e2
         e2
 
+    TESTS:
+
+    We can change the generator prefix::
+
+        sage: RealCartesianProductEJA(3, prefix='r').gens()
+        (r0, r1, r2)
+
     """
-    def __init__(self, n, field=QQ):
-        # The superclass constructor takes a list of matrices, the ith
-        # representing right multiplication by the ith basis element
-        # in the vector space. So if e_1 = (1,0,0), then right
-        # (Hadamard) multiplication of x by e_1 picks out the first
-        # component of x; and likewise for the ith basis element e_i.
-        Qs = [ matrix(field, n, n, lambda k,j: 1*(k == j == i))
-               for i in xrange(n) ]
+    def __init__(self, n, field=QQ, **kwargs):
+        V = VectorSpace(field, n)
+        mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ]
+                       for i in range(n) ]
 
         fdeja = super(RealCartesianProductEJA, self)
-        return fdeja.__init__(field, Qs, rank=n)
+        return fdeja.__init__(field, mult_table, rank=n, **kwargs)
 
     def inner_product(self, x, y):
         return _usual_ip(x,y)
@@ -583,7 +768,7 @@ def random_eja():
     TESTS::
 
         sage: random_eja()
-        Euclidean Jordan algebra of degree...
+        Euclidean Jordan algebra of dimension...
 
     """
 
@@ -715,10 +900,7 @@ def _multiplication_table_from_matrix_basis(basis):
     multiplication on the right is matrix multiplication. Given a basis
     for the underlying matrix space, this function returns a
     multiplication table (obtained by looping through the basis
-    elements) for an algebra of those matrices. A reordered copy
-    of the basis is also returned to work around the fact that
-    the ``span()`` in this function will change the order of the basis
-    from what we think it is, to... something else.
+    elements) for an algebra of those matrices.
     """
     # In S^2, for example, we nominally have four coordinates even
     # though the space is of dimension three only. The vector space V
@@ -730,20 +912,14 @@ def _multiplication_table_from_matrix_basis(basis):
 
     V = VectorSpace(field, dimension**2)
     W = V.span_of_basis( _mat2vec(s) for s in basis )
+    n = len(basis)
+    mult_table = [[W.zero() for j in range(n)] for i in range(n)]
+    for i in range(n):
+        for j in range(n):
+            mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
+            mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
 
-    Qs = []
-    for s in basis:
-        # Brute force the multiplication-by-s matrix by looping
-        # through all elements of the basis and doing the computation
-        # to find out what the corresponding row should be.
-        Q_cols = []
-        for t in basis:
-            this_col = _mat2vec((s*t + t*s)/2)
-            Q_cols.append(W.coordinates(this_col))
-        Q = matrix.column(field, W.dimension(), Q_cols)
-        Qs.append(Q)
-
-    return Qs
+    return mult_table
 
 
 def _embed_complex_matrix(M):
@@ -1008,12 +1184,12 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
 
     TESTS:
 
-    The degree of this algebra is `(n^2 + n) / 2`::
+    The dimension of this algebra is `(n^2 + n) / 2`::
 
         sage: set_random_seed()
         sage: n = ZZ.random_element(1,5)
         sage: J = RealSymmetricEJA(n)
-        sage: J.degree() == (n^2 + n)/2
+        sage: J.dimension() == (n^2 + n)/2
         True
 
     The Jordan multiplication is what we think it is::
@@ -1032,8 +1208,13 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
         sage: J(expected) == x*y
         True
 
+    We can change the generator prefix::
+
+        sage: RealSymmetricEJA(3, prefix='q').gens()
+        (q0, q1, q2, q3, q4, q5)
+
     """
-    def __init__(self, n, field=QQ):
+    def __init__(self, n, field=QQ, **kwargs):
         S = _real_symmetric_basis(n, field=field)
         Qs = _multiplication_table_from_matrix_basis(S)
 
@@ -1041,7 +1222,8 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
         return fdeja.__init__(field,
                               Qs,
                               rank=n,
-                              natural_basis=S)
+                              natural_basis=S,
+                              **kwargs)
 
     def inner_product(self, x, y):
         return _matrix_ip(x,y)
@@ -1060,12 +1242,12 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
 
     TESTS:
 
-    The degree of this algebra is `n^2`::
+    The dimension of this algebra is `n^2`::
 
         sage: set_random_seed()
         sage: n = ZZ.random_element(1,5)
         sage: J = ComplexHermitianEJA(n)
-        sage: J.degree() == n^2
+        sage: J.dimension() == n^2
         True
 
     The Jordan multiplication is what we think it is::
@@ -1084,8 +1266,13 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
         sage: J(expected) == x*y
         True
 
+    We can change the generator prefix::
+
+        sage: ComplexHermitianEJA(2, prefix='z').gens()
+        (z0, z1, z2, z3)
+
     """
-    def __init__(self, n, field=QQ):
+    def __init__(self, n, field=QQ, **kwargs):
         S = _complex_hermitian_basis(n)
         Qs = _multiplication_table_from_matrix_basis(S)
 
@@ -1093,7 +1280,8 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
         return fdeja.__init__(field,
                               Qs,
                               rank=n,
-                              natural_basis=S)
+                              natural_basis=S,
+                              **kwargs)
 
 
     def inner_product(self, x, y):
@@ -1120,12 +1308,12 @@ class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
 
     TESTS:
 
-    The degree of this algebra is `n^2`::
+    The dimension of this algebra is `n^2`::
 
         sage: set_random_seed()
         sage: n = ZZ.random_element(1,5)
         sage: J = QuaternionHermitianEJA(n)
-        sage: J.degree() == 2*(n^2) - n
+        sage: J.dimension() == 2*(n^2) - n
         True
 
     The Jordan multiplication is what we think it is::
@@ -1144,8 +1332,13 @@ class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
         sage: J(expected) == x*y
         True
 
+    We can change the generator prefix::
+
+        sage: QuaternionHermitianEJA(2, prefix='a').gens()
+        (a0, a1, a2, a3, a4, a5)
+
     """
-    def __init__(self, n, field=QQ):
+    def __init__(self, n, field=QQ, **kwargs):
         S = _quaternion_hermitian_basis(n)
         Qs = _multiplication_table_from_matrix_basis(S)
 
@@ -1153,7 +1346,8 @@ class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
         return fdeja.__init__(field,
                               Qs,
                               rank=n,
-                              natural_basis=S)
+                              natural_basis=S,
+                              **kwargs)
 
     def inner_product(self, x, y):
         # Since a+bi+cj+dk on the diagonal is represented as
@@ -1200,26 +1394,34 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
         sage: e2*e3
         0
 
+    We can change the generator prefix::
+
+        sage: JordanSpinEJA(2, prefix='B').gens()
+        (B0, B1)
+
     """
-    def __init__(self, n, field=QQ):
-        Qs = []
-        id_matrix = matrix.identity(field, n)
-        for i in xrange(n):
-            ei = id_matrix.column(i)
-            Qi = matrix.zero(field, n)
-            Qi.set_row(0, ei)
-            Qi.set_column(0, ei)
-            Qi += matrix.diagonal(n, [ei[0]]*n)
-            # The addition of the diagonal matrix adds an extra ei[0] in the
-            # upper-left corner of the matrix.
-            Qi[0,0] = Qi[0,0] * ~field(2)
-            Qs.append(Qi)
+    def __init__(self, n, field=QQ, **kwargs):
+        V = VectorSpace(field, n)
+        mult_table = [[V.zero() for j in range(n)] for i in range(n)]
+        for i in range(n):
+            for j in range(n):
+                x = V.gen(i)
+                y = V.gen(j)
+                x0 = x[0]
+                xbar = x[1:]
+                y0 = y[0]
+                ybar = y[1:]
+                # z = x*y
+                z0 = x.inner_product(y)
+                zbar = y0*xbar + x0*ybar
+                z = V([z0] + zbar.list())
+                mult_table[i][j] = z
 
         # The rank of the spin algebra is two, unless we're in a
         # one-dimensional ambient space (because the rank is bounded by
         # the ambient dimension).
         fdeja = super(JordanSpinEJA, self)
-        return fdeja.__init__(field, Qs, rank=min(n,2))
+        return fdeja.__init__(field, mult_table, rank=min(n,2), **kwargs)
 
     def inner_product(self, x, y):
         return _usual_ip(x,y)