# Library for dealing with bilinear forms, and related topics like
# angles of vectors, GramSchmidt, etc.

# Evaluate (v,w) with respect to the bilinear form given by the matrix A.
#FIXME: check parameters
SetHelp ("BilinearForm", "linear_algebra", "Evaluate (v,w) with respect to the bilinear form given by the matrix A")
function BilinearForm(v,A,w) = (v.'*A*w)@(1)
protect ("BilinearForm")

# Evaluate (v,w) with respect to the sesquilinear form given by the matrix A.
#FIXME: check parameters
SetHelp ("SesquilinearForm", "linear_algebra", "Evaluate (v,w) with respect to the sesquilinear form given by the matrix A")
function SesquilinearForm(v,A,w) = (v'*A*w)@(1)
protect ("SesquilinearForm")

# The angle of two vectors, given an inner product
#FIXME: default inner product!
#FIXME: check parameters
SetHelp ("VectorAngle", "linear_algebra", "The angle of two vectors, given an inner product")
function VectorAngle(v,w,B...) =
  (
        # if you don't give anything, assume standard inner product
	if IsNull(B) or IsNull(B@(1)) then P=HermitianProduct
	# If it's a matrix, then use the sesquilinear form defined by that matrix
	else if IsMatrix(B@(1)) then (P = SesquilinearFormFunction(u,B@(1),w))
	# If B is just some function, assume that it's an inner product
        else if elements(B) > 1 or not IsFunction(B@(1)) then (error("VectorAngle: argument B not a function or Matrix");bailout)
	else P = B@(1);
 	
	acos(P(v,w)/(P(v,v)*P(w,w)))
)
protect ("VectorAngle")

#FIXME: check parameters
SetHelp ("BilinearFormFunction", "linear_algebra", "Return a function that evaluates two vectors with respect to the bilinear form given by A")
function BilinearFormFunction(A)=(`(v,w)=(v.'*A*w)@(1))
protect ("BilinearFormFunction")

SetHelp ("SesquilinearFormFunction", "linear_algebra", "Return a function that evaluates two vectors with respect to the sesquilinear form given by A")
function SesquilinearFormFunction(A)=(`(v,w)=(v'*A*w)@(1))
protect ("SesquilinearFormFunction")

# Projection onto a vector space
# Projection of vector v onto subspace W given a sesquilinear form B
SetHelp ("Projection", "linear_algebra", "Projection of vector v onto subspace W given a sesquilinear form B (if not given use hermitian product)")
function Projection(v,W,B...) =
  (
        # if you don't give anything, assume standard inner product
	if IsNull(B) or IsNull(B@(1)) then P=HermitianProduct
	# If it's a matrix, then use the sesquilinear form defined by that matrix
	else if IsMatrix(B@(1)) then (P = SesquilinearFormFunction(u,B@(1),w))
	# If B is just some function, assume that it's an inner product
        else if elements(B) > 1 or not IsFunction(B@(1)) then (error("Projection: argument B not a function or Matrix");bailout)
	else P = B@(1);

	sum c in ColumnsOf(W) do (
		(P(v,c)/P(c,c))*c
	)
  )
protect ("Projection")

# Gram-Schmidt Orthogonalization
# Described, for instance, in Hoffman & Kunze, ``Linear Algebra'' p.~280
SetHelp ("GramSchmidt", "linear_algebra", "Apply the Gram-Schmidt process (to the columns) with respect to inner product given by B (if not given use hermitian product)")
function GramSchmidt(v,B...) =
# Takes k column vectors v_1,...,v_k
# and returns a collection, orthonormal with respect to the inner product given by B(-,-): V x V -> R
# Note: if you plug in a bad function or a matrix that isn't symmetric and
# positive definite, then you'll get garbage as a result.
	(
        # if you don't give anything, assume standard inner product
	if IsNull(B) or IsNull(B@(1)) then P=HermitianProduct
	# If it's a matrix, then use the sesquilinear form defined by that matrix
	else if IsMatrix(B@(1)) then (function P(u,w)=SesquilinearForm(u,B@(1),w))
	# If B is just some function, assume that it's an inner product
        else if elements(B) > 1 or not IsFunction(B@(1)) then (error("GramSchmidt: argument B not a function or Matrix");bailout)
	else P = B@(1);

#FIXME: this is not the most stable algorithm
	w=v@(,1);			# First vector
	w=w/sqrt(P(w,w));
	for i = 2 to columns(v) do (	# now iterate
		ww = v@(,i)-Projection(v@(,i),w,P);
		w = [w, ww / sqrt (P(ww,ww))]
	);
	w
	)
protect ("GramSchmidt")