]>
gitweb.michael.orlitzky.com - sage.d.git/blob - 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
5 def _charpoly_sage_input(s
):
7 Helper function that you can use on the string output from sage
8 to convert a charpoly coefficient into the corresponding input
13 sage: from mjo.eja.eja_algebra import JordanSpinEJA
14 sage: from mjo.eja.eja_utils import _charpoly_sage_input
18 sage: J = JordanSpinEJA(4,QQ)
19 sage: a = J._charpoly_coefficients()
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'
31 digit_out
= r
"X([0-9]+)"
35 return "X[" + str(int(m
.group(1)) - 1) + "]"
37 s
= re
.sub(exponent_out
, exponent_in
, s
)
38 return re
.sub(digit_out
, replace_digit
, s
)
43 Scale the vector, matrix, or cartesian-product-of-those-things
46 This works around the inability to scale certain elements of
47 Cartesian product spaces, as reported in
49 https://trac.sagemath.org/ticket/31435
53 This will do the wrong thing if you feed it a tuple or list.
57 sage: from mjo.eja.eja_utils import _scale
61 sage: v = vector(QQ, (1,2,3))
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)
72 if hasattr(x
, 'cartesian_factors'):
74 return P(tuple( _scale(x_i
, alpha
)
75 for x_i
in x
.cartesian_factors() ))
82 Flatten a vector, matrix, or cartesian product of those things
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.
92 sage: from mjo.eja.eja_utils import _all2list
93 sage: from mjo.hurwitz import (QuaternionMatrixAlgebra,
95 ....: OctonionMatrixAlgebra)
99 sage: _all2list([[1]])
104 sage: V1 = VectorSpace(QQ,2)
105 sage: V2 = MatrixSpace(QQ,2)
107 sage: x2 = V1([1,-1])
109 sage: y2 = V2([0,1,1,0])
110 sage: _all2list((x1,y1))
112 sage: _all2list((x2,y2))
114 sage: M = cartesian_product([V1,V2])
115 sage: _all2list(M((x1,y1)))
117 sage: _all2list(M((x2,y2)))
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]
129 sage: _all2list(QuaternionAlgebra(QQ, -1, -1).one())
131 sage: _all2list(QuaternionMatrixAlgebra(1).one())
136 sage: V1 = VectorSpace(QQ,2)
137 sage: V2 = OctonionMatrixAlgebra(1,field=QQ)
138 sage: C = cartesian_product([V1,V2])
141 sage: _all2list(C( (x1,y1) ))
142 [3, 4, 1, 0, 0, 0, 0, 0, 0, 0]
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())
156 except TypeError: # x is not iterable
160 # Avoid the retardation of list(QQ(1)) == [1].
163 return sum(list( map(_all2list
, xl
) ), [])
168 return vector(m
.base_ring(), m
.list())
171 return matrix(v
.base_ring(), sqrt(v
.degree()), v
.list())
173 def gram_schmidt(v
, inner_product
=None):
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.
182 sage: from mjo.eja.eja_utils import gram_schmidt
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::
189 sage: v1 = vector(QQ, (1,0,0))
190 sage: v2 = vector(QQ, (0,1,0))
191 sage: v3 = vector(QQ, (0,0,1))
193 sage: gram_schmidt(v) == v
196 The usual inner-product and norm are default::
198 sage: v1 = vector(AA,(1,2,3))
199 sage: v2 = vector(AA,(1,-1,6))
200 sage: v3 = vector(AA,(2,1,-1))
202 sage: u = gram_schmidt(v)
203 sage: all( u_i.inner_product(u_i).sqrt() == 1 for u_i in u )
205 sage: bool(u[0].inner_product(u[1]) == 0)
207 sage: bool(u[0].inner_product(u[2]) == 0)
209 sage: bool(u[1].inner_product(u[2]) == 0)
213 But if you supply a custom inner product, the result is
214 orthonormal with respect to that (and not the usual inner
217 sage: v1 = vector(AA,(1,2,3))
218 sage: v2 = vector(AA,(1,-1,6))
219 sage: v3 = vector(AA,(2,1,-1))
221 sage: B = matrix(AA, [ [6, 4, 2],
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 )
229 sage: ip(u[0],u[1]).is_zero()
231 sage: ip(u[0],u[2]).is_zero()
233 sage: ip(u[1],u[2]).is_zero()
236 This Gram-Schmidt routine can be used on matrices as well, so long
237 as an appropriate inner-product is provided::
239 sage: E11 = matrix(AA, [ [1,0],
241 sage: E12 = matrix(AA, [ [0,1],
243 sage: E22 = matrix(AA, [ [0,0],
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)
249 [1 0] [ 0 0.7071067811865475?] [0 0]
250 [0 0], [0.7071067811865475? 0], [0 1]
253 It even works on Cartesian product spaces whose factors are vector
256 sage: V1 = VectorSpace(AA,2)
257 sage: V2 = MatrixSpace(AA,2)
258 sage: M = cartesian_product([V1,V2])
260 sage: x2 = V1([1,-1])
262 sage: y2 = V2([0,1,1,0])
263 sage: z1 = M((x1,y1))
264 sage: z2 = M((x2,y2))
266 ....: return a[0].inner_product(b[0]) + (a[1]*b[1]).trace()
267 sage: U = gram_schmidt([z1,z2], inner_product=ip)
277 Ensure that zero vectors don't get in the way::
279 sage: v1 = vector(AA,(1,2,3))
280 sage: v2 = vector(AA,(1,-1,6))
281 sage: v3 = vector(AA,(0,0,0))
283 sage: len(gram_schmidt(v)) == 2
286 if inner_product
is None:
287 inner_product
= lambda x
,y
: x
.inner_product(y
)
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())
297 v
= list(v
) # make a copy, don't clobber the input
299 # Drop all zero vectors before we start.
300 v
= [ v_i
for v_i
in v
if not v_i
.is_zero() ]
308 # Our "zero" needs to belong to the right space for sum() to work.
309 zero
= v
[0].parent().zero()
312 if hasattr(v
[0], 'cartesian_factors'):
313 # Only use the slow implementation if necessary.
317 return sc(x
, (inner_product(x
,y
)/inner_product(x
,x
)))
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
])
324 if not v
[j
].is_zero() ),
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() ]
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
]))