# Some linear algebra functions

# Transpose of a matrix
SetHelp ("Transpose", "linear_algebra", "Transpose of a matrix")
function Transpose(M)=M.'
protect ("Transpose")

# Conjugate Transpose of a matrix
SetHelp ("ConjugateTranspose", "linear_algebra", "Conjugate transpose of a matrix (adjoint)")
function ConjugateTranspose(M)=M'
protect ("Conjugate Transpose")

#FIXME: We'll need more sophisticated tools for Tensors

# PivotColumns is now built-in for speed

# Non-pivot columns
SetHelp ("NonPivotColumns", "linear_algebra", "Return the columns that are not the pivot columns of a matrix")
function NonPivotColumns(M) =
  (
    non_pivot = [1:columns(M)];
    pivots = PivotColumns(M);
    if not IsNull (pivots) then
      SetMinus (non_pivot, pivots@(1,))
    else
      non_pivot
  )
protect ("NonPivotColumns")

# Delete a row from a matrix
# Shamelessly copied from DeleteRowColumn
SetHelp("DeleteRow", "matrix", "Delete a row of a matrix")
function DeleteRow(M,row) = (
        row_dim_M=rows(M);
        if row_dim_M==1 then return null;   # if delete the only row...
	     if row==1         then M@(2:row_dim_M,)
        else if row==row_dim_M then M@(1:row_dim_M-1,)
        else                       [M@(1:row-1,);M@(row+1:row_dim_M,)]
)
protect ("DeleteRow")

# Delete a column from a matrix
# Also shamelessly copied from DeleteRowColumn
SetHelp("DeleteColumn", "matrix", "Delete a column of a matrix")
function DeleteColumn(M,col) = (
        col_dim_M=columns(M);
        if col_dim_M==1 then return null;   # if delete the only column...
	     if col==1         then M@(,2:col_dim_M)
        else if col==col_dim_M then M@(,1:col_dim_M-1)
        else                       [M@(,1:col-1),M@(,col+1:col_dim_M)]
)
protect ("DeleteColumn")

# Column Reduced Echelon Form
SetHelp("cref", "linear_algebra", "Compute the Column Reduced Echelon Form")
function cref(M)=(rref(M.').')
protect ("cref")
SetHelpAlias ("cref", "CREF")
CREF = cref
protect ("CREF")
SetHelpAlias ("cref", "ColumnReducedEchelonForm")
ColumnReducedEchelonForm = cref
protect ("ColumnReducedEchelonForm")

SetHelp ("StripZeroRows", "matrix", "Removes any all-zero rows of M")
function StripZeroRows (M) =
(
  if IsNull (M) then return null
  else if not IsMatrix (M) then
    (error("StripZeroColumns: argument not a matrix");bailout);
  StripZeroColumns (M.').'
)
protect ("StripZeroRows")

# Cross product of two vectors (in R^3)
##<refentry id="fun-linear-algebra-cross-product">
##<indexterm zone="fun-linear-algebra-cross-product">
##<primary>CrossProduct</primary>
##</indexterm>
##<indexterm>
##<primary>cross product</primary>
##<see>CrossProduct</see>
##</indexterm>
##
##  <refmeta>
##    <refentrytitle>CrossProduct</refentrytitle>
##  </refmeta>
##
##  <refnamediv>
##    <refname>CrossProduct</refname>
##    <refpurpose>compute the cross product of two vectors</refpurpose>
##  </refnamediv>
##
##      <funcprototype>
##        <funcdef>Vector <function>CrossProduct</function></funcdef>
##        <paramdef>Vector <parameter>v</parameter></paramdef>
##        <paramdef>Vector <parameter>w</parameter></paramdef>
##      </funcprototype>
##    </funcsynopsis>
##  </refsynopsisdiv>
##
##  <refsect1>
##    <title>Description</title>
##    <para><function>CrossProduct</function> computes the cross product
##    of the vectors <parameter>v</parameter> and <parameter>w</parameter>,
##    which need to be vectors in R<sup>3</sup>. FIXME: generalize to n-1
##    vectors in R<sup>n</sup>.</para>
##    <para>FIXME: define it, talk about it</para>
##  </refsect1>
##  <refsect1>
##    <title>Examples</title>
##    <para>Some examples...FIXME</para>
##  </refsect1>
##</refentry>
SetHelp ("CrossProduct", "linear_algebra", "CrossProduct of two vectors in R^3")
function CrossProduct(v,w) =
 (
  M=[v;w;1,1,1];
  [Minor(M,3,1),Minor(M,3,2),Minor(M,3,3)]
 )
protect ("CrossProduct")

# Direct sum of a vector of matrices
SetHelp ("DirectSumMatrixVector", "linear_algebra", "Direct sum of a vector of matrices")
function DirectSumMatrixVector(v) = (
  if not IsMatrix (v) or rows (v) > 1 then
    (error("DirectSumMatrixVector: argument not a horizontal vector");bailout);
  if columns (v) == 1 then
    [v@(1)]
  else if columns (v) == 2 then
    [v@(1),0;0,v@(2)]
  else
    [v@(1),0;0,DirectSumMatrixVector(v@(1,2:columns(v)))] 
)
protect ("DirectSumMatrixVector");

# Direct sum of matrices
SetHelp ("DirectSum", "linear_algebra", "Direct sum of matrices")
function DirectSum(M,N...) = (
  if IsNull (N) then
    [M] 
  else if columns(N) == 1 then
    [M,0;0,N@(1)] 
  else
    [M,0;0,DirectSumMatrixVector(N)] 
)
protect ("DirectSum");

# (Kronecker) Tensor product of two matrices
#FIXME!
# function KroneckerProduct(M,N) = TensorProduct(M,N)
# function TensorProduct(M,N) = Null

# Checks to see if a matrix/value is invertible
# Note that given an integer matrix, returns true iff it is invertible OVER THE INTEGERS,
# so for instance [2,0;0,1] is not invertible
# FIXME: We need good coercion rules!
SetHelp ("IsInvertible", "linear_algebra", "Is a matrix (or number) invertible (Integer matrix is invertible iff it's invertible over the integers)")
function IsInvertible(n) =
(
  if IsNull (n) then return false
 else if (IsInteger(n)) then return (|n|==1)     # Only invertible integers are 1, -1
 else if (IsValue(n)) then return (n!=0)    # Elements of fields (Q,R,C) are invertible iff != 0
 else if (IsMatrix(n)) then (if IsMatrixSquare(n) then return (IsInvertible(det(n))) else return false);
                                             # Matrices are invertible iff they are square and their
                                             #  determinant is a unit
 error ("IsInvertible: argument not a value or matrix");bailout
)
protect ("IsInvertible")

SetHelp ("IsInvertibleField", "linear_algebra", "Is a matrix (or number) invertible over a field")
function IsInvertibleField(n) =
(
  if IsNull (n) then return false
 else if (IsValue(n)) then return (n!=0)    # Elements of fields (Q,R,C) are invertible iff != 0
 else if (IsMatrix(n)) then (if IsMatrixSquare(n) then return (IsInvertible(det(n))) else return true);
                                             # Matrices are invertible iff they are square and their
                                             #  determinant is a unit
 error ("IsInvertibleField: argument not a value or matrix");bailout
)
protect ("IsInvertibleField")

# SmithNormal Form for fields (will end up with 1's on the diagonal)
SetHelp("SmithNormalFormField", "linear_algebra", "Smith Normal Form for fields (will end up with 1's on the diagonal)")
function SmithNormalFormField(A) =
(
  ref (ref (A).').'
)
protect ("SmithNormalFormField")

SetHelp("SmithNormalFormInteger", "linear_algebra", "Smith Normal Form for square integer matrices (not its characteristic)")
function SmithNormalFormInteger(M) =
(
	if not IsMatrix (M) or not IsMatrixInteger(M) or not IsMatrixSquare(M) then
		(error ("SmithNormalFormInteger: M must be a square integer matrix");bailout);
	d = DeterminantalDivisorsInteger (M);
	S = zeros (rows (M), columns (M));
	S@(1,1) = d@(1);
	for i=2 to rows(M) do
		S@(i,i) = d@(i)/d@(i-1);
	S
)
protect ("SmithNormalFormInteger")

# AuxilliaryUnitMatrix
SetHelp("AuxilliaryUnitMatrix", "linear_algebra", "Get the auxilliary unit matrix of size n")
function AuxilliaryUnitMatrix(n) =
(
  if not IsInteger (n) or n < 2 then
    (error ("AuxilliaryUnitMatrix: n must be an integer greater than 1");bailout)
  else
    [0,I(n-1);0,0]
)
protect ("AuxilliaryUnitMatrix")

# AuxilliaryUnitMatrix
SetHelp("JordanBlock", "linear_algebra", "Get the jordan block corresponding to lambda and n")
function JordanBlock(n,lambda) =
(
  if not IsInteger (n) or n < 1 then
    (error ("JordanBlock: n must be an integer greater than 0");bailout)
  else if not IsValue (n) or n < 1 then
    (error ("JordanBlock: lambda must be a value");bailout)
  else if n == 1 then
    [lambda]
  else
    lambda * I(n) + AuxilliaryUnitMatrix(n)
)
protect ("JordanBlock")
SetHelpAlias("JordanBlock", "J")
J = JordanBlock
protect ("J")

# Dot product
SetHelp("DotProduct", "matrix", "Get the dot product of two vectors (no conjugates)")
function DotProduct(u,v) =
(
  if not IsVector (u) or not IsVector (v) then
    (error ("DotProduct: u,v must be vectors");bailout)
  else if columns (u) == 1 and columns (v) == 1 then
    ( u .' * v )@(1)
  else if rows (u) == 1 and columns (v) == 1 then
    ( u * v )@(1)
  else if rows (u) == 1 and rows (v) == 1 then
    ( u * (v.') )@(1)
  else # if columns (u) == 1 and rows (v) == 1 then
    ( ( u .' ) * (v.') )@(1)
)
protect ("DotProduct")

# Outer product
SetHelp("OuterProduct", "matrix", "Get the outer product of two vectors")
function OuterProduct(u,v) =
(
  if not IsVector (u) or not IsVector (v) then
    (error ("OuterProduct: u,v must be vectors");bailout)
  else if columns (u) == 1 and columns (v) == 1 then
    u * (v')
  else if rows (u) == 1 and columns (v) == 1 then
    (u.') * (v')
  else if rows (u) == 1 and rows (v) == 1 then
    (u.') * conj(v)
  else # if columns (u) == 1 and rows (v) == 1 then
    u * conj(v)
)
protect ("OuterProduct")

SetHelp("InfNorm", "linear_algebra", "Get the Inf Norm of a vector");
function InfNorm(v) = (
  if not IsVector (v) then
    (error ("InfNorm: v must be a vector");bailout)
  else
    max (|v|)
)
protect("InfNorm")

SetHelp("Norm", "linear_algebra", "Get the p Norm (or 2 Norm if no p is supplied) of a vector");
function Norm(v,p...) = (
  if not IsVector (v) then
    (error ("Norm: v must be a vector");bailout)
  else if IsNull(p) then
    sqrt (HermitianProduct (v,v))
  else if elements(p) > 1 or not IsInteger(p@(1)) then
    (error ("Norm: p must be an integer or not specified");bailout)
  else
    (sum k in v do (|k|^p@(1)))^(1/p@(1))
)
protect("Norm")
SetHelpAlias ("Norm", "norm")
norm = Norm
protect("norm")

SetHelp ("CharacteristicPolynomial", "linear_algebra", "Get the characteristic polynomial as a vector")
function CharacteristicPolynomial(M) = (
	if not IsMatrix (M) or not IsMatrixSquare(M) then
		(error ("CharacteristicPolynomial: M must be a square matrix");bailout);

	n = columns(M);
	v = [ det (M) ];
	for i=2 to n-1 do
		v@(i) = (-1) ^ (i-1) * sum gamma in Combinations(n-i+1,n) do det (Submatrix (M, gamma, gamma));
	v@(n) = (-1) ^ (n-1) * trace (M);
	v@(n+1) = (-1) ^ n;
	v
)
protect("CharacteristicPolynomial")
SetHelpAlias ("CharacteristicPolynomial", "CharPoly")
CharPoly = CharacteristicPolynomial
protect("CharPoly")

SetHelp ("CharacteristicPolynomialFunction", "linear_algebra", "Get the characteristic polynomial as a function")
function CharacteristicPolynomialFunction(M) = PolyToFunction (CharacteristicPolynomial (M))
protect("CharacteristicPolynomialFunction")

SetHelp ("DeterminantalDivisorsInteger", "linear_algebra", "Get the determinantal divisors of an integer matrix (not its characteristic)")
function DeterminantalDivisorsInteger(M) = (
	if not IsMatrix (M) or not IsMatrixInteger(M) or not IsMatrixSquare(M) then
		(error ("DeterminantalDivisorsInteger: M must be a square integer matrix");bailout);
	v = [gcd(M)];
	n = rows (M);
	for k=2 to n-1 do (
		g = 0;
		for gamma in Combinations(k,n) do
			for omega in Combinations(k,n) do
				g = gcd (g, det (Submatrix (M, gamma, omega)));
		v@(k) = g
	);
	v@(n) = |det (M)|;
	v
)
protect("DeterminantalDivisorsInteger")

SetHelp("InvariantFactorsInteger", "linear_algebra", "Get the invariant factors of a square integer matrix (not its characteristic)")
function InvariantFactorsInteger(M) =
(
	if not IsMatrix (M) or not IsMatrixInteger(M) or not IsMatrixSquare(M) then
		(error ("InvariantFactorsInteger: M must be a square integer matrix");bailout);
	d = DeterminantalDivisorsInteger (M);
	H = [d@(1)];
	for i=2 to rows(M) do
		H@(i) = d@(i)/d@(i-1);
	H
)
protect ("InvariantFactorsInteger")

SetHelp("IsNormal", "linear_algebra", "Is a matrix normal")
function IsNormal(M) = (
	if not IsMatrix(M) or not IsMatrixSquare(M) then
		(error("IsNormal: argument not a square matrix");bailout);
	M' * M == M * M'
)
protect ("IsNormal")

SetHelp("IsUnitary", "linear_algebra", "Is a matrix unitary")
function IsUnitary(M) = (
	if not IsMatrix(M) or not IsMatrixSquare(M) then
		(error("IsUnitary: argument not a square matrix");bailout);
	M' * M == I(columns(M))
)
protect ("IsUnitary")

SetHelp("IsHermitian", "linear_algebra", "Is a matrix hermitian")
function IsHermitian(M) = (
	if not IsMatrix(M) or not IsMatrixSquare(M) then
		(error("IsHermitian: argument not a square matrix");bailout);
	M' == M
)
protect ("IsHermitian")

SetHelp("IsSkewHermitian", "linear_algebra", "Is a matrix skew-hermitian")
function IsSkewHermitian(M) = (
	if not IsMatrix(M) or not IsMatrixSquare(M) then
		(error("IsSkewHermitian: argument not a square matrix");bailout);
	M' == -M
)
protect ("IsSkewHermitian")

SetHelp("IsPositiveDefinite", "linear_algebra", "Is a matrix positive definite")
function IsPositiveDefinite(M) = (
	if not IsMatrix(M) or not IsMatrixSquare(M) then
		(error("IsPositiveDefinite: argument not a square matrix");bailout);
	if M' == M then (
		n = rows(M);
		for k=1 to n do (
			for c in Combinations(k,n) do (
				if det(Submatrix (M,c,c)) <= 0 then
					return false
			)
		);
		true
	) else
		false
)
protect ("IsPositiveDefinite")

SetHelp("IsPositiveSemidefinite", "linear_algebra", "Is a matrix positive semidefinite")
function IsPositiveSemidefinite(M) = (
	if not IsMatrix(M) or not IsMatrixSquare(M) then
		(error("IsPositiveSemidefinite: argument not a square matrix");bailout);
	if M' == M then (
		n = rows(M);
		for k=1 to n do (
			for c in Combinations(k,n) do (
				if det(Submatrix (M,c,c)) < 0 then
					return false
			)
		);
		true
	) else
		false
)
protect ("IsPositiveSemidefinite")

SetHelp("RayleighQuotient", "linear_algebra", "Return the Rayleigh quotient of a matrix and a vector")
function RayleighQuotient(A,x) = (
	# FIXME: argument checking
	if columns (x) > 1 then x = x.';
	( ( x' * A * x ) / ( x' * x ) ) @(1)
)
protect ("RayleighQuotient")

SetHelp("Eigenvalues", "linear_algebra", "Get the eigenvalues of a matrix (Currently only for up to 4x4 or triangular matrices)")
function Eigenvalues(M) = (
	if not IsMatrix (M) or not IsMatrixSquare (M) then
		(error("Eigenvalues: argument not a square matrix");bailout);
	if IsUpperTriangular (M) or IsLowerTriangular (M) then (
		DiagonalOf (M)
	) else if columns (M) == 2 then (
		b = - trace (M) ; c = det (M);
		[ ( -b - sqrt (b^2 - 4*c) ) / 2 ; ( -b + sqrt (b^2 - 4*c) ) / 2 ]
	) else if columns (M) == 3 then (
		CubicFormula (CharacteristicPolynomial(M))
	) else if columns (M) == 4 then (
		QuarticFormula (CharacteristicPolynomial(M))
	) else
		(error("Eigenvalues: Currently only implemented for up to 4x4 or triangular matrices");bailout)
)
protect ("Eigenvalues")
eig = Eigenvalues;
SetHelpAlias ("Eigenvalues", "eig")
protect ("eig")

# Written by David W. Hutchison, dahutchi@indiana.edu for Genius 0.7.1
# Copyright (C) 2004 David W. Hutchison (GPL)
# LU decomposition of a matrix, aka Crout and/or Cholesky reduction.
# (ISBN 0-201-11577-8 pp.99-103)
# This function will compute the lower and upper materices of a square matrix
# 'a.' The upper triangular matrix features a diagonal of values 1 (one).  This
# is not Doolittle's Method which features the 1's diagonal on the lower
# matrix. 
SetHelp("LUDecomposition", "linear_algebra", "Get the LU decomposition of A and store the result in the L and U which should be references.  If not possible returns false.")
function LUDecomposition(A,L,U) = ( 
	if not IsMatrix(A) or not IsMatrixSquare(A) then
        	(error("LUDecomposition: A must be a square matrix");bailout);
	if not IsFunctionRef(L) or not IsFunctionRef(U) then
        	(error("LUDecomposition: L and U must be references");bailout);
	*L = *U = null;
	if A@(1,1) == 0 then (
		return false
	);
	n = rows(A);

	lower = zeros(n,n);
	lower@(,1) = A@(,1);
	upper = I(n);
	upper@(1,) = A@(1,) / lower@(1,1);

	for j=2 to n do (
		for i=j to n do (
			lower@(i,j) = A@(i,j) -
				(sum k=1 to j-1 do
					(lower@(i,k)*upper@(k,j)))
		);
		ljj = lower@(j,j);
		if ljj == 0 then (
			return false
		);
		for i = j+1 to n do (
			upper@(j,i) = (A@(j,i) -
				(sum k=1 to j-1 do
					(lower@(j,k)*upper@(k,i)))) /
				ljj
		)
	);
	*L = lower;
	*U = upper;
	true
);
protect("LUDecomposition");

# FIXME: Should use more stable algorithm
SetHelp("QRDecomposition", "linear_algebra", "Get the QR decomposition of A, returns R and Q can be a reference")
function QRDecomposition(A,Q) = ( 
	if not IsMatrix(A) or not IsMatrixSquare(A) then
        	(error("QRDecomposition: A must be a square matrix");bailout);
	if not IsFunctionRef(Q) and not IsNull(Q) then
        	(error("QRDecomposition: Q must be a reference or null");bailout);

	QQ = GramSchmidt (A);

	if IsFunctionRef(Q) then *Q = QQ;
	
	UpperTriangular (QQ' * A)
)
protect("QRDecomposition")

# See for example: http://www.cs.utk.edu/~dongarra/etemplates/node97.html
SetHelp("RayleighQuotientIteration", "linear_algebra", "Compute an eigenvalue using the Rayleigh Quotient Iteration method until we are epsilon from eigenvalue or for maxiter iterations")
function RayleighQuotientIteration(A,x,epsilon,maxiter,vecref) = (
	if not IsMatrix(A) or not IsMatrixSquare(A) then
        	(error("RayleighQuotientIteration: A must be a square matrix");bailout);
	if not IsVector(x) or elements(x) != columns(A) then
        	(error("RayleighQuotientIteration: x must be a vector of same size as A");bailout);
	if not IsValue(epsilon) or not IsPositiveInteger(maxiter) then
        	(error("RayleighQuotientIteration: epsilon and maxiter must be numbers");bailout);
	if not IsFunctionRef(vecref) and not IsNull(vecref) then
        	(error("RayleighQuotientIteration: vecref must be a reference or null");bailout);

	if columns (x) > 1 then x = x.';

	x = x / sqrt (HermitianProduct (x,x));

	p = ( ( x' * A * x ) / ( x' * x ) ) @(1);

	for k = 1 to maxiter do (
		y = SolveLinearSystem ((A-p),x);
		if IsNull(y) then (
		        if IsFunctionRef(vecref) then *vecref = x;
			return p
		);

		nsq = HermitianProduct (y,y);

		p = p + (y'*x)@(1) / nsq;

		x = y/sqrt(nsq);

		if nsq >= 1/epsilon then (
		        if IsFunctionRef(vecref) then *vecref = x;
			return p
		)
	);
	null
)
protect ("RayleighQuotientIteration")