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