From 3c1c9170143a6412b050602cd79a383e3da6c821 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Fri, 27 Nov 2020 09:16:11 -0500
Subject: [PATCH] eja: allow non-standard inner product in gram_schmidt.

---
 mjo/eja/eja_utils.py | 42 ++++++++++++++++++++++++++++++++++++------
 1 file changed, 36 insertions(+), 6 deletions(-)

diff --git a/mjo/eja/eja_utils.py b/mjo/eja/eja_utils.py
index 49e3078..4d70e06 100644
--- a/mjo/eja/eja_utils.py
+++ b/mjo/eja/eja_utils.py
@@ -11,7 +11,7 @@ def _mat2vec(m):
 def _vec2mat(v):
         return matrix(v.base_ring(), sqrt(v.degree()), v.list())
 
-def gram_schmidt(v):
+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
@@ -22,7 +22,9 @@ def gram_schmidt(v):
 
         sage: from mjo.eja.eja_utils import gram_schmidt
 
-    EXAMPLES::
+    EXAMPLES:
+
+    The usual inner-product and norm are default::
 
         sage: v1 = vector(QQ,(1,2,3))
         sage: v2 = vector(QQ,(1,-1,6))
@@ -38,6 +40,30 @@ def gram_schmidt(v):
         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
+
     TESTS:
 
     Ensure that zero vectors don't get in the way::
@@ -50,8 +76,12 @@ def gram_schmidt(v):
         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 (y.inner_product(x)/x.inner_product(x))*x
+        return (inner_product(x,y)/inner_product(x,x))*x
 
     v = list(v) # make a copy, don't clobber the input
 
@@ -65,7 +95,7 @@ def gram_schmidt(v):
     R = v[0].base_ring()
 
     # First orthogonalize...
-    for i in xrange(1,len(v)):
+    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() )
 
@@ -75,7 +105,7 @@ def gram_schmidt(v):
     # 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 xrange(len(v)):
-        v[i] = v[i] / v[i].norm()
+    for i in range(len(v)):
+        v[i] = v[i] / norm(v[i])
 
     return v
-- 
2.44.2