mjo/eja/eja_utils.py

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