mjo/eja/eja_utils.py

   1 from sage.functions.other import sqrt
   2 from sage.structure.element import is_Matrix
   3 from sage.matrix.constructor import matrix
   4 from sage.modules.free_module_element import vector
   5
   6 def _charpoly_sage_input(s):
   7     r"""
   8     Helper function that you can use on the string output from sage
   9     to convert a charpoly coefficient into the corresponding input
  10     to be cached.
  11
  12     SETUP::
  13
  14         sage: from mjo.eja.eja_algebra import JordanSpinEJA
  15         sage: from mjo.eja.eja_utils import _charpoly_sage_input
  16
  17     EXAMPLES::
  18
  19         sage: J = JordanSpinEJA(4,QQ)
  20         sage: a = J._charpoly_coefficients()
  21         sage: a[0]
  22         X1^2 - X2^2 - X3^2 - X4^2
  23         sage: _charpoly_sage_input(str(a[0]))
  24         'X[0]**2 - X[1]**2 - X[2]**2 - X[3]**2'
  25
  26     """
  27     import re
  28
  29     exponent_out = r"\^"
  30     exponent_in = r"**"
  31
  32     digit_out = r"X([0-9]+)"
  33
  34     def replace_digit(m):
  35         # m is a match object
  36         return "X[" + str(int(m.group(1)) - 1) + "]"
  37
  38     s = re.sub(exponent_out, exponent_in, s)
  39     return re.sub(digit_out, replace_digit, s)
  40
  41
  42 def _scale(x, alpha):
  43     r"""
  44     Scale the vector, matrix, or cartesian-product-of-those-things
  45     ``x`` by ``alpha``.
  46
  47     This works around the inability to scale certain elements of
  48     Cartesian product spaces, as reported in
  49
  50       https://trac.sagemath.org/ticket/31435
  51
  52     ..WARNING:
  53
  54         This will do the wrong thing if you feed it a tuple or list.
  55
  56     SETUP::
  57
  58         sage: from mjo.eja.eja_utils import _scale
  59
  60     EXAMPLES::
  61
  62         sage: v = vector(QQ, (1,2,3))
  63         sage: _scale(v,2)
  64         (2, 4, 6)
  65         sage: m = matrix(QQ, [[1,2],[3,4]])
  66         sage: M = cartesian_product([m.parent(), m.parent()])
  67         sage: _scale(M((m,m)), 2)
  68         ([2 4]
  69         [6 8], [2 4]
  70         [6 8])
  71
  72     """
  73     if hasattr(x, 'cartesian_factors'):
  74         P = x.parent()
  75         return P(tuple( _scale(x_i, alpha)
  76                         for x_i in x.cartesian_factors() ))
  77     else:
  78         return x*alpha
  79
  80
  81 def _all2list(x):
  82     r"""
  83     Flatten a vector, matrix, or cartesian product of those things
  84     into a long list.
  85
  86     If the entries of the matrix themselves belong to a real vector
  87     space (such as the complex numbers which can be thought of as
  88     pairs of real numbers), they will also be expanded in vector form
  89     and flattened into the list.
  90
  91     SETUP::
  92
  93         sage: from mjo.eja.eja_utils import _all2list
  94         sage: from mjo.hurwitz import (QuaternionMatrixAlgebra,
  95         ....:                          Octonions,
  96         ....:                          OctonionMatrixAlgebra)
  97
  98     EXAMPLES::
  99
 100         sage: _all2list([[1]])
 101         [1]
 102
 103     ::
 104
 105         sage: V1 = VectorSpace(QQ,2)
 106         sage: V2 = MatrixSpace(QQ,2)
 107         sage: x1 = V1([1,1])
 108         sage: x2 = V1([1,-1])
 109         sage: y1 = V2.one()
 110         sage: y2 = V2([0,1,1,0])
 111         sage: _all2list((x1,y1))
 112         [1, 1, 1, 0, 0, 1]
 113         sage: _all2list((x2,y2))
 114         [1, -1, 0, 1, 1, 0]
 115         sage: M = cartesian_product([V1,V2])
 116         sage: _all2list(M((x1,y1)))
 117         [1, 1, 1, 0, 0, 1]
 118         sage: _all2list(M((x2,y2)))
 119         [1, -1, 0, 1, 1, 0]
 120
 121     ::
 122
 123         sage: _all2list(Octonions().one())
 124         [1, 0, 0, 0, 0, 0, 0, 0]
 125         sage: _all2list(OctonionMatrixAlgebra(1).one())
 126         [1, 0, 0, 0, 0, 0, 0, 0]
 127
 128     ::
 129
 130         sage: _all2list(QuaternionAlgebra(QQ, -1, -1).one())
 131         [1, 0, 0, 0]
 132         sage: _all2list(QuaternionMatrixAlgebra(1).one())
 133         [1, 0, 0, 0]
 134
 135     ::
 136
 137         sage: V1 = VectorSpace(QQ,2)
 138         sage: V2 = OctonionMatrixAlgebra(1,field=QQ)
 139         sage: C = cartesian_product([V1,V2])
 140         sage: x1 = V1([3,4])
 141         sage: y1 = V2.one()
 142         sage: _all2list(C( (x1,y1) ))
 143         [3, 4, 1, 0, 0, 0, 0, 0, 0, 0]
 144
 145     """
 146     if hasattr(x, 'to_vector'):
 147         # This works on matrices of e.g. octonions directly, without
 148         # first needing to convert them to a list of octonions and
 149         # then recursing down into the list. It also avoids the wonky
 150         # list(x) when x is an element of a CFM. I don't know what it
 151         # returns but it aint the coordinates. We don't recurse
 152         # because vectors can only contain ring elements as entries.
 153         return x.to_vector().list()
 154
 155     if is_Matrix(x):
 156         # This sucks, but for performance reasons we don't want to
 157         # call _all2list recursively on the contents of a matrix
 158         # when we don't have to (they only contain ring elements
 159         # as entries)
 160         return x.list()
 161
 162     try:
 163         xl = list(x)
 164     except TypeError: # x is not iterable
 165         return [x]
 166
 167     if xl == [x]:
 168         # Avoid the retardation of list(QQ(1)) == [1].
 169         return [x]
 170
 171     return sum( map(_all2list, xl) , [])
 172
 173
 174
 175 def _mat2vec(m):
 176         return vector(m.base_ring(), m.list())
 177
 178 def _vec2mat(v):
 179         return matrix(v.base_ring(), sqrt(v.degree()), v.list())
 180
 181 def gram_schmidt(v, inner_product=None):
 182     """
 183     Perform Gram-Schmidt on the list ``v`` which are assumed to be
 184     vectors over the same base ring. Returns a list of orthonormalized
 185     vectors over the same base ring, which means that your base ring
 186     needs to contain the appropriate roots.
 187
 188     SETUP::
 189
 190         sage: from mjo.eja.eja_utils import gram_schmidt
 191
 192     EXAMPLES:
 193
 194     If you start with an orthonormal set, you get it back. We can use
 195     the rationals here because we don't need any square roots::
 196
 197         sage: v1 = vector(QQ, (1,0,0))
 198         sage: v2 = vector(QQ, (0,1,0))
 199         sage: v3 = vector(QQ, (0,0,1))
 200         sage: v = [v1,v2,v3]
 201         sage: gram_schmidt(v) == v
 202         True
 203
 204     The usual inner-product and norm are default::
 205
 206         sage: v1 = vector(AA,(1,2,3))
 207         sage: v2 = vector(AA,(1,-1,6))
 208         sage: v3 = vector(AA,(2,1,-1))
 209         sage: v = [v1,v2,v3]
 210         sage: u = gram_schmidt(v)
 211         sage: all( u_i.inner_product(u_i).sqrt() == 1 for u_i in u )
 212         True
 213         sage: bool(u[0].inner_product(u[1]) == 0)
 214         True
 215         sage: bool(u[0].inner_product(u[2]) == 0)
 216         True
 217         sage: bool(u[1].inner_product(u[2]) == 0)
 218         True
 219
 220
 221     But if you supply a custom inner product, the result is
 222     orthonormal with respect to that (and not the usual inner
 223     product)::
 224
 225         sage: v1 = vector(AA,(1,2,3))
 226         sage: v2 = vector(AA,(1,-1,6))
 227         sage: v3 = vector(AA,(2,1,-1))
 228         sage: v = [v1,v2,v3]
 229         sage: B = matrix(AA, [ [6, 4, 2],
 230         ....:                  [4, 5, 4],
 231         ....:                  [2, 4, 9] ])
 232         sage: ip = lambda x,y: (B*x).inner_product(y)
 233         sage: norm = lambda x: ip(x,x)
 234         sage: u = gram_schmidt(v,ip)
 235         sage: all( norm(u_i) == 1 for u_i in u )
 236         True
 237         sage: ip(u[0],u[1]).is_zero()
 238         True
 239         sage: ip(u[0],u[2]).is_zero()
 240         True
 241         sage: ip(u[1],u[2]).is_zero()
 242         True
 243
 244     This Gram-Schmidt routine can be used on matrices as well, so long
 245     as an appropriate inner-product is provided::
 246
 247         sage: E11 = matrix(AA, [ [1,0],
 248         ....:                    [0,0] ])
 249         sage: E12 = matrix(AA, [ [0,1],
 250         ....:                    [1,0] ])
 251         sage: E22 = matrix(AA, [ [0,0],
 252         ....:                    [0,1] ])
 253         sage: I = matrix.identity(AA,2)
 254         sage: trace_ip = lambda X,Y: (X*Y).trace()
 255         sage: gram_schmidt([E11,E12,I,E22], inner_product=trace_ip)
 256         [
 257         [1 0]  [                  0 0.7071067811865475?]  [0 0]
 258         [0 0], [0.7071067811865475?                   0], [0 1]
 259         ]
 260
 261     It even works on Cartesian product spaces whose factors are vector
 262     or matrix spaces::
 263
 264         sage: V1 = VectorSpace(AA,2)
 265         sage: V2 = MatrixSpace(AA,2)
 266         sage: M = cartesian_product([V1,V2])
 267         sage: x1 = V1([1,1])
 268         sage: x2 = V1([1,-1])
 269         sage: y1 = V2.one()
 270         sage: y2 = V2([0,1,1,0])
 271         sage: z1 = M((x1,y1))
 272         sage: z2 = M((x2,y2))
 273         sage: def ip(a,b):
 274         ....:     return a[0].inner_product(b[0]) + (a[1]*b[1]).trace()
 275         sage: U = gram_schmidt([z1,z2], inner_product=ip)
 276         sage: ip(U[0],U[1])
 277         0
 278         sage: ip(U[0],U[0])
 279         1
 280         sage: ip(U[1],U[1])
 281         1
 282
 283     TESTS:
 284
 285     Ensure that zero vectors don't get in the way::
 286
 287         sage: v1 = vector(AA,(1,2,3))
 288         sage: v2 = vector(AA,(1,-1,6))
 289         sage: v3 = vector(AA,(0,0,0))
 290         sage: v = [v1,v2,v3]
 291         sage: len(gram_schmidt(v)) == 2
 292         True
 293     """
 294     if inner_product is None:
 295         inner_product = lambda x,y: x.inner_product(y)
 296     def norm(x):
 297         ip = inner_product(x,x)
 298         # Don't expand the given field; the inner-product's codomain
 299         # is already correct. For example QQ(2).sqrt() returns sqrt(2)
 300         # in SR, and that will give you weird errors about symbolics
 301         # when what's really going wrong is that you're trying to
 302         # orthonormalize in QQ.
 303         return ip.parent()(ip.sqrt())
 304
 305     v = list(v) # make a copy, don't clobber the input
 306
 307     # Drop all zero vectors before we start.
 308     v = [ v_i for v_i in v if not v_i.is_zero() ]
 309
 310     if len(v) == 0:
 311         # cool
 312         return v
 313
 314     R = v[0].base_ring()
 315
 316     # Our "zero" needs to belong to the right space for sum() to work.
 317     zero = v[0].parent().zero()
 318
 319     sc = lambda x,a: a*x
 320     if hasattr(v[0], 'cartesian_factors'):
 321         # Only use the slow implementation if necessary.
 322         sc = _scale
 323
 324     def proj(x,y):
 325         return sc(x, (inner_product(x,y)/inner_product(x,x)))
 326
 327     # First orthogonalize...
 328     for i in range(1,len(v)):
 329         # Earlier vectors can be made into zero so we have to ignore them.
 330         v[i] -= sum( (proj(v[j],v[i])
 331                       for j in range(i)
 332                       if not v[j].is_zero() ),
 333                      zero )
 334
 335     # And now drop all zero vectors again if they were "orthogonalized out."
 336     v = [ v_i for v_i in v if not v_i.is_zero() ]
 337
 338     # Just normalize. If the algebra is missing the roots, we can't add
 339     # them here because then our subalgebra would have a bigger field
 340     # than the superalgebra.
 341     for i in range(len(v)):
 342         v[i] = sc(v[i], ~norm(v[i]))
 343
 344     return v