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