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