eja: make gram_schmidt work in Cartesian product algebras.

[sage.d.git] / mjo / eja / eja_utils.py

diff --git a/mjo/eja/eja_utils.py b/mjo/eja/eja_utils.py
index 4d70e062c44e98ba90584fb0f97d3bbfef91e491..e8ed4db2ea23def442af71402cac98c5f5dfccf6 100644 (file)
--- a/mjo/eja/eja_utils.py
+++ b/mjo/eja/eja_utils.py
@@ -1,9 +1,22 @@
 from sage.functions.other import sqrt
 from sage.matrix.constructor import matrix
 from sage.modules.free_module_element import vector
-from sage.rings.number_field.number_field import NumberField
-from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
-from sage.rings.real_lazy import RLF
+
+def _all2list(x):
+    r"""
+    Flatten a vector, matrix, or cartesian product of those things
+    into a long list.
+    """
+    if hasattr(x, 'list'):
+        # Easy case...
+        return x.list()
+    if hasattr(x, 'cartesian_factors'):
+        # If it's a formal cartesian product space element, then
+        # we also know what to do...
+        return sum(( x_i.list() for x_i in x ), [])
+    else:
+        # But what if it's a tuple or something else?
+        return sum( map(_all2list,x), [] )
 
 def _mat2vec(m):
         return vector(m.base_ring(), m.list())
@@ -64,6 +77,23 @@ def gram_schmidt(v, inner_product=None):
         sage: ip(u[1],u[2]).is_zero()
         True
 
+    This Gram-Schmidt routine can be used on matrices as well, so long
+    as an appropriate inner-product is provided::
+
+        sage: E11 = matrix(QQ, [ [1,0],
+        ....:                    [0,0] ])
+        sage: E12 = matrix(QQ, [ [0,1],
+        ....:                    [1,0] ])
+        sage: E22 = matrix(QQ, [ [0,0],
+        ....:                    [0,1] ])
+        sage: I = matrix.identity(QQ,2)
+        sage: trace_ip = lambda X,Y: (X*Y).trace()
+        sage: gram_schmidt([E11,E12,I,E22], inner_product=trace_ip)
+        [
+        [1 0]  [          0 1/2*sqrt(2)]  [0 0]
+        [0 0], [1/2*sqrt(2)           0], [0 1]
+        ]
+
     TESTS:
 
     Ensure that zero vectors don't get in the way::
@@ -80,9 +110,6 @@ def gram_schmidt(v, inner_product=None):
         inner_product = lambda x,y: x.inner_product(y)
     norm = lambda x: inner_product(x,x).sqrt()
 
-    def proj(x,y):
-        return (inner_product(x,y)/inner_product(x,x))*x
-
     v = list(v) # make a copy, don't clobber the input
 
     # Drop all zero vectors before we start.
@@ -94,10 +121,26 @@ def gram_schmidt(v, inner_product=None):
 
     R = v[0].base_ring()
 
+    # Define a scaling operation that can be used on tuples.
+    # Oh and our "zero" needs to belong to the right space.
+    scale = lambda x,alpha: x*alpha
+    zero = v[0].parent().zero()
+    if hasattr(v[0], 'cartesian_factors'):
+        P = v[0].parent()
+        scale = lambda x,alpha: P(tuple( x_i*alpha
+                                         for x_i in x.cartesian_factors() ))
+
+
+    def proj(x,y):
+        return scale(x, (inner_product(x,y)/inner_product(x,x)))
+
     # First orthogonalize...
     for i in range(1,len(v)):
         # Earlier vectors can be made into zero so we have to ignore them.
-        v[i] -= sum( proj(v[j],v[i]) for j in range(i) if not v[j].is_zero() )
+        v[i] -= sum( (proj(v[j],v[i])
+                      for j in range(i)
+                      if not v[j].is_zero() ),
+                     zero )
 
     # And now drop all zero vectors again if they were "orthogonalized out."
     v = [ v_i for v_i in v if not v_i.is_zero() ]
@@ -106,6 +149,6 @@ def gram_schmidt(v, inner_product=None):
     # them here because then our subalgebra would have a bigger field
     # than the superalgebra.
     for i in range(len(v)):
-        v[i] = v[i] / norm(v[i])
+        v[i] = scale(v[i], ~norm(v[i]))
 
     return v