]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: add subalgebra() method.
[sage.d.git] / mjo / eja / eja_algebra.py
index ccc5006f3d79397d81ad98bbf45a7262e36dd6e9..8005412e65da8b6192072420dffe794c7975b5d0 100644 (file)
@@ -31,7 +31,8 @@ from sage.modules.free_module import FreeModule, VectorSpace
 from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
                             PolynomialRing,
                             QuadraticField)
-from mjo.eja.eja_element import FiniteDimensionalEJAElement
+from mjo.eja.eja_element import (CartesianProductEJAElement,
+                                 FiniteDimensionalEJAElement)
 from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
 from mjo.eja.eja_utils import _mat2vec
 
@@ -41,7 +42,15 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
     INPUT:
 
-      - basis -- a tuple of basis elements in their matrix form.
+      - basis -- a tuple of basis elements in "matrix form," which
+        must be the same form as the arguments to ``jordan_product``
+        and ``inner_product``. In reality, "matrix form" can be either
+        vectors, matrices, or a Cartesian product (ordered tuple)
+        of vectors or matrices. All of these would ideally be vector
+        spaces in sage with no special-casing needed; but in reality
+        we turn vectors into column-matrices and Cartesian products
+        `(a,b)` into column matrices `(a,b)^{T}` after converting
+        `a` and `b` themselves.
 
       - jordan_product -- function of two elements (in matrix form)
         that returns their jordan product in this algebra; this will
@@ -62,10 +71,10 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                  field=AA,
                  orthonormalize=True,
                  associative=False,
+                 cartesian_product=False,
                  check_field=True,
                  check_axioms=True,
-                 prefix='e',
-                 category=None):
+                 prefix='e'):
 
         if check_field:
             if not field.is_subring(RR):
@@ -76,7 +85,12 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         # If the basis given to us wasn't over the field that it's
         # supposed to be over, fix that. Or, you know, crash.
-        basis = tuple( b.change_ring(field) for b in basis )
+        if not cartesian_product:
+            # The field for a cartesian product algebra comes from one
+            # of its factors and is the same for all factors, so
+            # there's no need to "reapply" it on product algebras.
+            basis = tuple( b.change_ring(field) for b in basis )
+
 
         if check_axioms:
             # Check commutativity of the Jordan and inner-products.
@@ -94,12 +108,13 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                 raise ValueError("inner-product is not commutative")
 
 
-        if category is None:
-            category = MagmaticAlgebras(field).FiniteDimensional()
-            category = category.WithBasis().Unital()
-            if associative:
-                # Element subalgebras can take advantage of this.
-                category = category.Associative()
+        category = MagmaticAlgebras(field).FiniteDimensional()
+        category = category.WithBasis().Unital()
+        if associative:
+            # Element subalgebras can take advantage of this.
+            category = category.Associative()
+        if cartesian_product:
+            category = category.CartesianProducts()
 
         # Call the superclass constructor so that we can use its from_vector()
         # method to build our multiplication table.
@@ -117,10 +132,17 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # we see in things like x = 1*e1 + 2*e2.
         vector_basis = basis
 
+        def flatten(b):
+            # flatten a vector, matrix, or cartesian product of those
+            # things into a long list.
+            if cartesian_product:
+                return sum(( b_i.list() for b_i in b ), [])
+            else:
+                return b.list()
+
         degree = 0
         if n > 0:
-            # Works on both column and square matrices...
-            degree = len(basis[0].list())
+            degree = len(flatten(basis[0]))
 
         # Build an ambient space that fits our matrix basis when
         # written out as "long vectors."
@@ -134,7 +156,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             # Save a copy of the un-orthonormalized basis for later.
             # Convert it to ambient V (vector) coordinates while we're
             # at it, because we'd have to do it later anyway.
-            deortho_vector_basis = tuple( V(b.list()) for b in basis )
+            deortho_vector_basis = tuple( V(flatten(b)) for b in basis )
 
             from mjo.eja.eja_utils import gram_schmidt
             basis = tuple(gram_schmidt(basis, inner_product))
@@ -146,7 +168,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # Now create the vector space for the algebra, which will have
         # its own set of non-ambient coordinates (in terms of the
         # supplied basis).
-        vector_basis = tuple( V(b.list()) for b in basis )
+        vector_basis = tuple( V(flatten(b)) for b in basis )
         W = V.span_of_basis( vector_basis, check=check_axioms)
 
         if orthonormalize:
@@ -178,7 +200,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(elt.list()))
+                elt = W.coordinate_vector(V(flatten(elt)))
                 self._multiplication_table[i][j] = self.from_vector(elt)
 
                 if not orthonormalize:
@@ -283,22 +305,32 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: y = J.random_element()
             sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
             True
+
         """
         B = self._inner_product_matrix
         return (B*x.to_vector()).inner_product(y.to_vector())
 
 
-    def _is_commutative(self):
+    def is_associative(self):
         r"""
-        Whether or not this algebra's multiplication table is commutative.
+        Return whether or not this algebra's Jordan product is associative.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
+
+        EXAMPLES::
+
+            sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
+            sage: J.is_associative()
+            False
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by(orthonormalize=False)
+            sage: A.is_associative()
+            True
 
-        This method should of course always return ``True``, unless
-        this algebra was constructed with ``check_axioms=False`` and
-        passed an invalid multiplication table.
         """
-        return all( self.product_on_basis(i,j) == self.product_on_basis(i,j)
-                    for i in range(self.dimension())
-                    for j in range(self.dimension()) )
+        return "Associative" in self.category().axioms()
 
     def _is_jordanian(self):
         r"""
@@ -307,13 +339,13 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         We only check one arrangement of `x` and `y`, so for a
         ``True`` result to be truly true, you should also check
-        :meth:`_is_commutative`. This method should of course always
+        :meth:`is_commutative`. This method should of course always
         return ``True``, unless this algebra was constructed with
         ``check_axioms=False`` and passed an invalid multiplication table.
         """
-        return all( (self.monomial(i)**2)*(self.monomial(i)*self.monomial(j))
+        return all( (self.gens()[i]**2)*(self.gens()[i]*self.gens()[j])
                     ==
-                    (self.monomial(i))*((self.monomial(i)**2)*self.monomial(j))
+                    (self.gens()[i])*((self.gens()[i]**2)*self.gens()[j])
                     for i in range(self.dimension())
                     for j in range(self.dimension()) )
 
@@ -335,9 +367,9 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         for i in range(self.dimension()):
             for j in range(self.dimension()):
                 for k in range(self.dimension()):
-                    x = self.monomial(i)
-                    y = self.monomial(j)
-                    z = self.monomial(k)
+                    x = self.gens()[i]
+                    y = self.gens()[j]
+                    z = self.gens()[k]
                     diff = (x*y).inner_product(z) - x.inner_product(y*z)
 
                     if self.base_ring().is_exact():
@@ -660,8 +692,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         # And to each subsequent row, prepend an entry that belongs to
         # the left-side "header column."
-        M += [ [self.monomial(i)] + [ self.product_on_basis(i,j)
-                                      for j in range(n) ]
+        M += [ [self.gens()[i]] + [ self.product_on_basis(i,j)
+                                    for j in range(n) ]
                for i in range(n) ]
 
         return table(M, header_row=True, header_column=True, frame=True)
@@ -997,14 +1029,12 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         if not c.is_idempotent():
             raise ValueError("element is not idempotent: %s" % c)
 
-        from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
-
         # Default these to what they should be if they turn out to be
         # trivial, because eigenspaces_left() won't return eigenvalues
         # 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 = FiniteDimensionalEJASubalgebra(self, ())
+        trivial = self.subalgebra(())
         J0 = trivial                          # eigenvalue zero
         J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
         J1 = trivial                          # eigenvalue one
@@ -1014,9 +1044,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                 J5 = eigspace
             else:
                 gens = tuple( self.from_vector(b) for b in eigspace.basis() )
-                subalg = FiniteDimensionalEJASubalgebra(self,
-                                                        gens,
-                                                        check_axioms=False)
+                subalg = self.subalgebra(gens, check_axioms=False)
                 if eigval == 0:
                     J0 = subalg
                 elif eigval == 1:
@@ -1129,7 +1157,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         def L_x_i_j(i,j):
             # From a result in my book, these are the entries of the
             # basis representation of L_x.
-            return sum( vars[k]*self.monomial(k).operator().matrix()[i,j]
+            return sum( vars[k]*self.gens()[k].operator().matrix()[i,j]
                         for k in range(n) )
 
         L_x = matrix(F, n, n, L_x_i_j)
@@ -1235,6 +1263,17 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         return len(self._charpoly_coefficients())
 
 
+    def subalgebra(self, basis, **kwargs):
+        r"""
+        Create a subalgebra of this algebra from the given basis.
+
+        This is a simple wrapper around a subalgebra class constructor
+        that can be overridden in subclasses.
+        """
+        from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra
+        return FiniteDimensionalEJASubalgebra(self, basis, **kwargs)
+
+
     def vector_space(self):
         """
         Return the vector space that underlies this algebra.
@@ -2362,7 +2401,11 @@ class HadamardEJA(ConcreteEJA):
         if "check_axioms" not in kwargs: kwargs["check_axioms"] = False
 
         column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() )
-        super().__init__(column_basis, jordan_product, inner_product, **kwargs)
+        super().__init__(column_basis,
+                         jordan_product,
+                         inner_product,
+                         associative=True,
+                         **kwargs)
         self.rank.set_cache(n)
 
         if n == 0:
@@ -2730,6 +2773,52 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
         Real Field (+) Euclidean Jordan algebra of dimension 6 over
         Algebraic Real Field
 
+    Rank is additive on a Cartesian product::
+
+        sage: J1 = HadamardEJA(1)
+        sage: J2 = RealSymmetricEJA(2)
+        sage: J = cartesian_product([J1,J2])
+        sage: J1.rank.clear_cache()
+        sage: J2.rank.clear_cache()
+        sage: J.rank.clear_cache()
+        sage: J.rank()
+        3
+        sage: J.rank() == J1.rank() + J2.rank()
+        True
+
+    The same rank computation works over the rationals, with whatever
+    basis you like::
+
+        sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False)
+        sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False)
+        sage: J = cartesian_product([J1,J2])
+        sage: J1.rank.clear_cache()
+        sage: J2.rank.clear_cache()
+        sage: J.rank.clear_cache()
+        sage: J.rank()
+        3
+        sage: J.rank() == J1.rank() + J2.rank()
+        True
+
+    The product algebra will be associative if and only if all of its
+    components are associative::
+
+        sage: J1 = HadamardEJA(2)
+        sage: J1.is_associative()
+        True
+        sage: J2 = HadamardEJA(3)
+        sage: J2.is_associative()
+        True
+        sage: J3 = RealSymmetricEJA(3)
+        sage: J3.is_associative()
+        False
+        sage: CP1 = cartesian_product([J1,J2])
+        sage: CP1.is_associative()
+        True
+        sage: CP2 = cartesian_product([J1,J3])
+        sage: CP2.is_associative()
+        False
+
     TESTS:
 
     All factors must share the same base field::
@@ -2741,45 +2830,55 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
         ...
         ValueError: all factors must share the same base field
 
-    The "cached" Jordan and inner products are the componentwise
-    ones::
+    The cached unit element is the same one that would be computed::
 
-        sage: set_random_seed()
-        sage: J1 = random_eja()
-        sage: J2 = random_eja()
-        sage: J = cartesian_product([J1,J2])
-        sage: x,y = J.random_elements(2)
-        sage: x*y == J.cartesian_jordan_product(x,y)
-        True
-        sage: x.inner_product(y) == J.cartesian_inner_product(x,y)
+        sage: set_random_seed()              # long time
+        sage: J1 = random_eja()              # long time
+        sage: J2 = random_eja()              # long time
+        sage: J = cartesian_product([J1,J2]) # long time
+        sage: actual = J.one()               # long time
+        sage: J.one.clear_cache()            # long time
+        sage: expected = J.one()             # long time
+        sage: actual == expected             # long time
         True
 
     """
-    def __init__(self, modules, **kwargs):
+    def __init__(self, algebras, **kwargs):
         CombinatorialFreeModule_CartesianProduct.__init__(self,
-                                                          modules,
+                                                          algebras,
                                                           **kwargs)
-        field = modules[0].base_ring()
-        if not all( J.base_ring() == field for J in modules ):
+        field = algebras[0].base_ring()
+        if not all( J.base_ring() == field for J in algebras ):
             raise ValueError("all factors must share the same base field")
 
-        basis = tuple( b.to_vector().column() for b in self.basis() )
+        associative = all( m.is_associative() for m in algebras )
+
+        # The definition of matrix_space() and self.basis() relies
+        # only on the stuff in the CFM_CartesianProduct class, which
+        # we've already initialized.
+        Js = self.cartesian_factors()
+        m = len(Js)
+        MS = self.matrix_space()
+        basis = tuple(
+            MS(tuple( self.cartesian_projection(i)(b).to_matrix()
+                      for i in range(m) ))
+            for b in self.basis()
+        )
 
-        # Define jordan/inner products that operate on thbasis.
-        def jordan_product(x_mat,y_mat):
-            x = self.from_vector(_mat2vec(x_mat))
-            y = self.from_vector(_mat2vec(y_mat))
-            return self.cartesian_jordan_product(x,y).to_vector().column()
+        # Define jordan/inner products that operate on that matrix_basis.
+        def jordan_product(x,y):
+            return MS(tuple(
+                (Js[i](x[i])*Js[i](y[i])).to_matrix() for i in range(m)
+            ))
 
-        def inner_product(x_mat, y_mat):
-            x = self.from_vector(_mat2vec(x_mat))
-            y = self.from_vector(_mat2vec(y_mat))
-            return self.cartesian_inner_product(x,y)
+        def inner_product(x, y):
+            return sum(
+                Js[i](x[i]).inner_product(Js[i](y[i])) for i in range(m)
+            )
 
-        # Use whatever category the superclass came up with. Usually
-        # some join of the EJA and Cartesian product
-        # categories. There's no need to check the field since it
-        # already came from an EJA.
+        # There's no need to check the field since it already came
+        # from an EJA. Likewise the axioms are guaranteed to be
+        # satisfied, unless the guy writing this class sucks.
         #
         # If you want the basis to be orthonormalized, orthonormalize
         # the factors.
@@ -2789,11 +2888,39 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
                                       inner_product,
                                       field=field,
                                       orthonormalize=False,
+                                      associative=associative,
+                                      cartesian_product=True,
                                       check_field=False,
-                                      check_axioms=False,
-                                      category=self.category())
+                                      check_axioms=False)
+
+        ones = tuple(J.one() for J in algebras)
+        self.one.set_cache(self._cartesian_product_of_elements(ones))
+        self.rank.set_cache(sum(J.rank() for J in algebras))
+
+    def matrix_space(self):
+        r"""
+        Return the space that our matrix basis lives in as a Cartesian
+        product.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
+            ....:                                  RealSymmetricEJA)
+
+        EXAMPLES::
 
-        self.rank.set_cache(sum(J.rank() for J in modules))
+            sage: J1 = HadamardEJA(1)
+            sage: J2 = RealSymmetricEJA(2)
+            sage: J = cartesian_product([J1,J2])
+            sage: J.matrix_space()
+            The Cartesian product of (Full MatrixSpace of 1 by 1 dense
+            matrices over Algebraic Real Field, Full MatrixSpace of 2
+            by 2 dense matrices over Algebraic Real Field)
+
+        """
+        from sage.categories.cartesian_product import cartesian_product
+        return cartesian_product( [J.matrix_space()
+                                   for J in self.cartesian_factors()] )
 
     @cached_method
     def cartesian_projection(self, i):
@@ -2979,84 +3106,43 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
         return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix())
 
 
-    def cartesian_jordan_product(self, x, y):
+    def _element_constructor_(self, elt):
         r"""
-        The componentwise Jordan product.
+        Construct an element of this algebra from an ordered tuple.
 
-        We project ``x`` and ``y`` onto our factors, and add up the
-        Jordan products from the subalgebras. This may still be useful
-        after (if) the default Jordan product in the Cartesian product
-        algebra is overridden.
+        We just apply the element constructor from each of our factors
+        to the corresponding component of the tuple, and package up
+        the result.
 
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (HadamardEJA,
-            ....:                                  JordanSpinEJA)
+            ....:                                  RealSymmetricEJA)
 
-        EXAMPLE::
+        EXAMPLES::
 
             sage: J1 = HadamardEJA(3)
-            sage: J2 = JordanSpinEJA(3)
-            sage: J = cartesian_product([J1,J2])
-            sage: x1 = J1.from_vector(vector(QQ,(1,2,1)))
-            sage: y1 = J1.from_vector(vector(QQ,(1,0,2)))
-            sage: x2 = J2.from_vector(vector(QQ,(1,2,3)))
-            sage: y2 = J2.from_vector(vector(QQ,(1,1,1)))
-            sage: z1 = J.from_vector(vector(QQ,(1,2,1,1,2,3)))
-            sage: z2 = J.from_vector(vector(QQ,(1,0,2,1,1,1)))
-            sage: (x1*y1).to_vector()
-            (1, 0, 2)
-            sage: (x2*y2).to_vector()
-            (6, 3, 4)
-            sage: J.cartesian_jordan_product(z1,z2).to_vector()
-            (1, 0, 2, 6, 3, 4)
-
-        """
-        m = len(self.cartesian_factors())
-        projections = ( self.cartesian_projection(i) for i in range(m) )
-        products = ( P(x)*P(y) for P in projections )
-        return self._cartesian_product_of_elements(tuple(products))
-
-    def cartesian_inner_product(self, x, y):
-        r"""
-        The standard componentwise Cartesian inner-product.
-
-        We project ``x`` and ``y`` onto our factors, and add up the
-        inner-products from the subalgebras. This may still be useful
-        after (if) the default inner product in the Cartesian product
-        algebra is overridden.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import (HadamardEJA,
-            ....:                                  QuaternionHermitianEJA)
-
-        EXAMPLE::
-
-            sage: J1 = HadamardEJA(3,field=QQ)
-            sage: J2 = QuaternionHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: J2 = RealSymmetricEJA(2)
             sage: J = cartesian_product([J1,J2])
-            sage: x1 = J1.one()
-            sage: x2 = x1
-            sage: y1 = J2.one()
-            sage: y2 = y1
-            sage: x1.inner_product(x2)
-            3
-            sage: y1.inner_product(y2)
-            2
-            sage: z1 = J._cartesian_product_of_elements((x1,y1))
-            sage: z2 = J._cartesian_product_of_elements((x2,y2))
-            sage: J.cartesian_inner_product(z1,z2)
-            5
-
+            sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
+            e(0, 1) + e(1, 2)
         """
         m = len(self.cartesian_factors())
-        projections = ( self.cartesian_projection(i) for i in range(m) )
-        return sum( P(x).inner_product(P(y)) for P in projections )
-
+        try:
+            z = tuple( self.cartesian_factors()[i](elt[i]) for i in range(m) )
+            return self._cartesian_product_of_elements(z)
+        except:
+            raise ValueError("not an element of this algebra")
 
-    Element = FiniteDimensionalEJAElement
+    Element = CartesianProductEJAElement
 
 
 FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA
+
 random_eja = ConcreteEJA.random_instance
+#def random_eja(*args, **kwargs):
+#    from sage.categories.cartesian_product import cartesian_product
+#    J1 = HadamardEJA(1, **kwargs)
+#    J2 = RealSymmetricEJA(2, **kwargs)
+#    J =  cartesian_product([J1,J2])
+#    return J