X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Feja%2Feja_utils.py;h=29edf5b8a339b073e1426a12a98a6143e7af5069;hb=667e0df9c07589c03616ad8cf42eebe5c86de50b;hp=b6b0a0327c19842ce654023b00bf6e0b9af539a1;hpb=af22a5320eb5a915ff50bc6f90f87195d6cd615d;p=sage.d.git

diff --git a/mjo/eja/eja_utils.py b/mjo/eja/eja_utils.py
index b6b0a03..29edf5b 100644
--- a/mjo/eja/eja_utils.py
+++ b/mjo/eja/eja_utils.py
@@ -1,4 +1,141 @@
+from sage.functions.other import sqrt
+from sage.matrix.constructor import matrix
 from sage.modules.free_module_element import vector
 
+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())
+
+def _vec2mat(v):
+        return matrix(v.base_ring(), sqrt(v.degree()), v.list())
+
+def gram_schmidt(v, inner_product=None):
+    """
+    Perform Gram-Schmidt on the list ``v`` which are assumed to be
+    vectors over the same base ring. Returns a list of orthonormalized
+    vectors over the smallest extention ring containing the necessary
+    roots.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_utils import gram_schmidt
+
+    EXAMPLES:
+
+    The usual inner-product and norm are default::
+
+        sage: v1 = vector(QQ,(1,2,3))
+        sage: v2 = vector(QQ,(1,-1,6))
+        sage: v3 = vector(QQ,(2,1,-1))
+        sage: v = [v1,v2,v3]
+        sage: u = gram_schmidt(v)
+        sage: all( u_i.inner_product(u_i).sqrt() == 1 for u_i in u )
+        True
+        sage: bool(u[0].inner_product(u[1]) == 0)
+        True
+        sage: bool(u[0].inner_product(u[2]) == 0)
+        True
+        sage: bool(u[1].inner_product(u[2]) == 0)
+        True
+
+
+    But if you supply a custom inner product, the result is
+    orthonormal with respect to that (and not the usual inner
+    product)::
+
+        sage: v1 = vector(QQ,(1,2,3))
+        sage: v2 = vector(QQ,(1,-1,6))
+        sage: v3 = vector(QQ,(2,1,-1))
+        sage: v = [v1,v2,v3]
+        sage: B = matrix(QQ, [ [6, 4, 2],
+        ....:                  [4, 5, 4],
+        ....:                  [2, 4, 9] ])
+        sage: ip = lambda x,y: (B*x).inner_product(y)
+        sage: norm = lambda x: ip(x,x)
+        sage: u = gram_schmidt(v,ip)
+        sage: all( norm(u_i) == 1 for u_i in u )
+        True
+        sage: ip(u[0],u[1]).is_zero()
+        True
+        sage: ip(u[0],u[2]).is_zero()
+        True
+        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::
+
+        sage: v1 = vector(QQ,(1,2,3))
+        sage: v2 = vector(QQ,(1,-1,6))
+        sage: v3 = vector(QQ,(0,0,0))
+        sage: v = [v1,v2,v3]
+        sage: len(gram_schmidt(v)) == 2
+        True
+
+    """
+    if inner_product is 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.
+    v = [ v_i for v_i in v if not v_i.is_zero() ]
+
+    if len(v) == 0:
+        # cool
+        return v
+
+    R = v[0].base_ring()
+
+    # 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() )
+
+    # 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() ]
+
+    # Just normalize. If the algebra is missing the roots, we can't add
+    # 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])
+
+    return v