]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: re-disable base ring to algebra coercion.
[sage.d.git] / mjo / eja / eja_algebra.py
index 71cdd6f133623d49469569f70ef9504f8b1efd10..689a3db016437d1e6eda5c6372e52a3513896671 100644 (file)
@@ -5,7 +5,7 @@ are used in optimization, and have some additional nice methods beyond
 what can be supported in a general Jordan Algebra.
 """
 
-from itertools import izip, repeat
+from itertools import repeat
 
 from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra
 from sage.categories.magmatic_algebras import MagmaticAlgebras
@@ -26,13 +26,30 @@ lazy_import('mjo.eja.eja_subalgebra',
 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 _coerce_map_from_base_ring(self):
+        """
+        Disable the map from the base ring into the algebra.
+
+        Performing a nonsense conversion like this automatically
+        is counterpedagogical. The fallback is to try the usual
+        element constructor, which should also fail.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import random_eja
+
+        TESTS::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: J(1)
+            Traceback (most recent call last):
+            ...
+            ValueError: not a naturally-represented algebra element
+
+        """
+        return None
 
     def __init__(self,
                  field,
@@ -94,8 +111,10 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         # 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 ]
+        self._multiplication_table = [
+            list(map(lambda x: self.from_vector(x), ls))
+            for ls in mult_table
+        ]
 
 
     def _element_constructor_(self, elt):
@@ -109,7 +128,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
-            ....:                                  RealCartesianProductEJA,
+            ....:                                  HadamardEJA,
             ....:                                  RealSymmetricEJA)
 
         EXAMPLES:
@@ -137,7 +156,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         vector representations) back and forth faithfully::
 
             sage: set_random_seed()
-            sage: J = RealCartesianProductEJA.random_instance()
+            sage: J = HadamardEJA.random_instance()
             sage: x = J.random_element()
             sage: J(x.to_vector().column()) == x
             True
@@ -147,15 +166,22 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             True
 
         """
+        msg = "not a naturally-represented algebra element"
         if elt == 0:
             # The superclass implementation of random_element()
             # needs to be able to coerce "0" into the algebra.
             return self.zero()
+        elif elt in self.base_ring():
+            # Ensure that no base ring -> algebra coercion is performed
+            # by this method. There's some stupidity in sage that would
+            # otherwise propagate to this method; for example, sage thinks
+            # that the integer 3 belongs to the space of 2-by-2 matrices.
+            raise ValueError(msg)
 
         natural_basis = self.natural_basis()
         basis_space = natural_basis[0].matrix_space()
         if elt not in basis_space:
-            raise ValueError("not a naturally-represented algebra element")
+            raise ValueError(msg)
 
         # Thanks for nothing! Matrix spaces aren't vector spaces in
         # Sage, so we have to figure out its natural-basis coordinates
@@ -407,7 +433,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         S = PolynomialRing(S, R.variable_names())
         t = S(t)
 
-        return sum( a[k]*(t**k) for k in xrange(len(a)) )
+        return sum( a[k]*(t**k) for k in range(len(a)) )
 
 
     def inner_product(self, x, y):
@@ -493,7 +519,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         """
         M = list(self._multiplication_table) # copy
-        for i in xrange(len(M)):
+        for i in range(len(M)):
             # M had better be "square"
             M[i] = [self.monomial(i)] + M[i]
         M = [["*"] + list(self.gens())] + M
@@ -580,12 +606,12 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import (RealCartesianProductEJA,
+            sage: from mjo.eja.eja_algebra import (HadamardEJA,
             ....:                                  random_eja)
 
         EXAMPLES::
 
-            sage: J = RealCartesianProductEJA(5)
+            sage: J = HadamardEJA(5)
             sage: J.one()
             e0 + e1 + e2 + e3 + e4
 
@@ -765,7 +791,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             True
 
         """
-        return  tuple( self.random_element() for idx in xrange(count) )
+        return  tuple( self.random_element() for idx in range(count) )
 
 
     def rank(self):
@@ -901,8 +927,7 @@ class KnownRankEJA(object):
         return cls(n, field, **kwargs)
 
 
-class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra,
-                              KnownRankEJA):
+class HadamardEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
     """
     Return the Euclidean Jordan Algebra corresponding to the set
     `R^n` under the Hadamard product.
@@ -913,13 +938,13 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra,
 
     SETUP::
 
-        sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
+        sage: from mjo.eja.eja_algebra import HadamardEJA
 
     EXAMPLES:
 
     This multiplication table can be verified by hand::
 
-        sage: J = RealCartesianProductEJA(3)
+        sage: J = HadamardEJA(3)
         sage: e0,e1,e2 = J.gens()
         sage: e0*e0
         e0
@@ -938,16 +963,16 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra,
 
     We can change the generator prefix::
 
-        sage: RealCartesianProductEJA(3, prefix='r').gens()
+        sage: HadamardEJA(3, prefix='r').gens()
         (r0, r1, r2)
 
     """
     def __init__(self, n, field=QQ, **kwargs):
         V = VectorSpace(field, n)
-        mult_table = [ [ V.gen(i)*(i == j) for j in xrange(n) ]
-                       for i in xrange(n) ]
+        mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ]
+                       for i in range(n) ]
 
-        fdeja = super(RealCartesianProductEJA, self)
+        fdeja = super(HadamardEJA, self)
         return fdeja.__init__(field, mult_table, rank=n, **kwargs)
 
     def inner_product(self, x, y):
@@ -956,7 +981,7 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra,
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
+            sage: from mjo.eja.eja_algebra import HadamardEJA
 
         TESTS:
 
@@ -964,7 +989,7 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra,
         over `R^n`::
 
             sage: set_random_seed()
-            sage: J = RealCartesianProductEJA.random_instance()
+            sage: J = HadamardEJA.random_instance()
             sage: x,y = J.random_elements(2)
             sage: X = x.natural_representation()
             sage: Y = y.natural_representation()
@@ -1032,7 +1057,7 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
                 basis = tuple( s.change_ring(field) for s in basis )
             self._basis_normalizers = tuple(
                 ~(self.natural_inner_product(s,s).sqrt()) for s in basis )
-            basis = tuple(s*c for (s,c) in izip(basis,self._basis_normalizers))
+            basis = tuple(s*c for (s,c) in zip(basis,self._basis_normalizers))
 
         Qs = self.multiplication_table_from_matrix_basis(basis)
 
@@ -1055,8 +1080,8 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
             # with had entries in a nice field.
             return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coeff(i)
         else:
-            basis = ( (b/n) for (b,n) in izip(self.natural_basis(),
-                                              self._basis_normalizers) )
+            basis = ( (b/n) for (b,n) in zip(self.natural_basis(),
+                                             self._basis_normalizers) )
 
             # Do this over the rationals and convert back at the end.
             J = MatrixEuclideanJordanAlgebra(QQ,
@@ -1066,7 +1091,7 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
             (_,x,_,_) = J._charpoly_matrix_system()
             p = J._charpoly_coeff(i)
             # p might be missing some vars, have to substitute "optionally"
-            pairs = izip(x.base_ring().gens(), self._basis_normalizers)
+            pairs = zip(x.base_ring().gens(), self._basis_normalizers)
             substitutions = { v: v*c for (v,c) in pairs }
             result = p.subs(substitutions)
 
@@ -1099,9 +1124,9 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         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 xrange(n)] for i in xrange(n)]
-        for i in xrange(n):
-            for j in xrange(n):
+        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))
 
@@ -1272,8 +1297,8 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra, KnownRankEJA):
         # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
         # coordinates.
         S = []
-        for i in xrange(n):
-            for j in xrange(i+1):
+        for i in range(n):
+            for j in range(i+1):
                 Eij = matrix(field, n, lambda k,l: k==i and l==j)
                 if i == j:
                     Sij = Eij
@@ -1403,8 +1428,8 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         # Go top-left to bottom-right (reading order), converting every
         # 2-by-2 block we see to a single complex element.
         elements = []
-        for k in xrange(n/2):
-            for j in xrange(n/2):
+        for k in range(n/2):
+            for j in range(n/2):
                 submat = M[2*k:2*k+2,2*j:2*j+2]
                 if submat[0,0] != submat[1,1]:
                     raise ValueError('bad on-diagonal submatrix')
@@ -1555,8 +1580,8 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra, KnownRankEJA):
         #   * The diagonal will (as a result) be real.
         #
         S = []
-        for i in xrange(n):
-            for j in xrange(i+1):
+        for i in range(n):
+            for j in range(i+1):
                 Eij = matrix(F, n, lambda k,l: k==i and l==j)
                 if i == j:
                     Sij = cls.real_embed(Eij)
@@ -1693,8 +1718,8 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
         # quaternion block.
         elements = []
-        for l in xrange(n/4):
-            for m in xrange(n/4):
+        for l in range(n/4):
+            for m in range(n/4):
                 submat = ComplexMatrixEuclideanJordanAlgebra.real_unembed(
                     M[4*l:4*l+4,4*m:4*m+4] )
                 if submat[0,0] != submat[1,1].conjugate():
@@ -1846,8 +1871,8 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
         #   * The diagonal will (as a result) be real.
         #
         S = []
-        for i in xrange(n):
-            for j in xrange(i+1):
+        for i in range(n):
+            for j in range(i+1):
                 Eij = matrix(Q, n, lambda k,l: k==i and l==j)
                 if i == j:
                     Sij = cls.real_embed(Eij)
@@ -1874,11 +1899,131 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
         super(QuaternionHermitianEJA,self).__init__(field, basis, n, **kwargs)
 
 
+class BilinearFormEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
+    r"""
+    The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
+    with the half-trace inner product and jordan product ``x*y =
+    (x0*y0 + <B*x_bar,y_bar>, x0*y_bar + y0*x_bar)`` where ``B`` is a
+    symmetric positive-definite "bilinear form" matrix. It has
+    dimension `n` over the reals, and reduces to the ``JordanSpinEJA``
+    when ``B`` is the identity matrix of order ``n-1``.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import (BilinearFormEJA,
+        ....:                                  JordanSpinEJA)
+
+    EXAMPLES:
+
+    When no bilinear form is specified, the identity matrix is used,
+    and the resulting algebra is the Jordan spin algebra::
+
+        sage: J0 = BilinearFormEJA(3)
+        sage: J1 = JordanSpinEJA(3)
+        sage: J0.multiplication_table() == J0.multiplication_table()
+        True
+
+    TESTS:
+
+    We can create a zero-dimensional algebra::
+
+        sage: J = BilinearFormEJA(0)
+        sage: J.basis()
+        Finite family {}
+
+    We can check the multiplication condition given in the Jordan, von
+    Neumann, and Wigner paper (and also discussed on my "On the
+    symmetry..." paper). Note that this relies heavily on the standard
+    choice of basis, as does anything utilizing the bilinear form matrix::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(5)
+        sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular')
+        sage: B = M.transpose()*M
+        sage: J = BilinearFormEJA(n, B=B)
+        sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
+        sage: V = J.vector_space()
+        sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
+        ....:         for ei in eis ]
+        sage: actual = [ sis[i]*sis[j]
+        ....:            for i in range(n-1)
+        ....:            for j in range(n-1) ]
+        sage: expected = [ J.one() if i == j else J.zero()
+        ....:              for i in range(n-1)
+        ....:              for j in range(n-1) ]
+        sage: actual == expected
+        True
+    """
+    def __init__(self, n, field=QQ, B=None, **kwargs):
+        if B is None:
+            self._B = matrix.identity(field, max(0,n-1))
+        else:
+            self._B = B
+
+        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:]
+                z0 = x0*y0 + (self._B*xbar).inner_product(ybar)
+                zbar = y0*xbar + x0*ybar
+                z = V([z0] + zbar.list())
+                mult_table[i][j] = z
+
+        # The rank of this algebra is two, unless we're in a
+        # one-dimensional ambient space (because the rank is bounded
+        # by the ambient dimension).
+        fdeja = super(BilinearFormEJA, self)
+        return fdeja.__init__(field, mult_table, rank=min(n,2), **kwargs)
+
+    def inner_product(self, x, y):
+        r"""
+        Half of the trace inner product.
+
+        This is defined so that the special case of the Jordan spin
+        algebra gets the usual inner product.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import BilinearFormEJA
+
+        TESTS:
+
+        Ensure that this is one-half of the trace inner-product::
+
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(5)
+            sage: M = matrix.random(QQ, n-1, algorithm='unimodular')
+            sage: B = M.transpose()*M
+            sage: J = BilinearFormEJA(n, B=B)
+            sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis()
+            sage: V = J.vector_space()
+            sage: sis = [ J.from_vector(V([0] + (M.inverse()*ei).list()))
+            ....:         for ei in eis ]
+            sage: actual = [ sis[i]*sis[j]
+            ....:            for i in range(n-1)
+            ....:            for j in range(n-1) ]
+            sage: expected = [ J.one() if i == j else J.zero()
+            ....:              for i in range(n-1)
+            ....:              for j in range(n-1) ]
+
+        """
+        xvec = x.to_vector()
+        xbar = xvec[1:]
+        yvec = y.to_vector()
+        ybar = yvec[1:]
+        return x[0]*y[0] + (self._B*xbar).inner_product(ybar)
+
 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
     """
     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
     with the usual inner product and jordan product ``x*y =
-    (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
+    (<x,y>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
     the reals.
 
     SETUP::
@@ -1914,9 +2059,9 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra, KnownRankEJA):
     """
     def __init__(self, n, field=QQ, **kwargs):
         V = VectorSpace(field, n)
-        mult_table = [[V.zero() for j in xrange(n)] for i in xrange(n)]
-        for i in xrange(n):
-            for j in xrange(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]