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