    -- Copyright (c) 2008-2015, David Amos. All rights reserved.

    -- |A module providing elementary operations involving scalars, vectors, and matrices
    -- over a ring or field. Vectors are represented as [a], matrices as [[a]].
    -- (No distinction is made between row and column vectors.)
    -- It is the caller's responsibility to ensure that the lists have the correct number of elements.
    --
    -- The mnemonic for many of the arithmetic operations is that the number of angle brackets
    -- on each side indicates the dimension of the argument on that side. For example,
    -- v \<*\>\> m is multiplication of a vector on the left by a matrix on the right.
    module Math.Algebra.LinearAlgebra where

    import Prelude hiding ( (*>), (<*>) )

    import qualified Data.List as L
    import Math.Core.Field -- not actually used in this module

    infixr 8 *>, *>>
    infixr 7 <<*>
    infixl 7 <.>, <*>, <<*>>, <*>>
    infixl 6 <+>, <->, <<+>>, <<->>

    -- The mnemonic for these operations is that the number of angle brackets on each side indicates the dimension of the argument on that side

    -- vector operations

    -- |u \<+\> v returns the sum u+v of vectors
    (<+>) :: (Num a) => [a] -> [a] -> [a]
    u <+> v = zipWith (+) u v

    -- |u \<-\> v returns the difference u-v of vectors
    (<->) :: (Num a) => [a] -> [a] -> [a]
    u <-> v = zipWith (-) u v

    -- |k *\> v returns the product k*v of the scalar k and the vector v
    (*>) :: (Num a) => a -> [a] -> [a]
    k *> v = map (k*) v

    -- |u \<.\> v returns the dot product of vectors (also called inner or scalar product)
    (<.>) :: (Num a) => [a] -> [a] -> a
    u <.> v = sum (zipWith (*) u v)

    -- |u \<*\> v returns the tensor product of vectors (also called outer or matrix product)
    (<*>) :: (Num a) => [a] -> [a] -> [[a]]
    u <*> v = [ [a*b | b <- v] | a <- u]

    -- matrix operations

    -- |a \<\<+\>\> b returns the sum a+b of matrices
    (<<+>>) :: (Num a) => [[a]] -> [[a]] -> [[a]]
    a <<+>> b = (zipWith . zipWith) (+) a b

    -- |a \<\<-\>\> b returns the difference a-b of matrices
    (<<->>) :: (Num a) => [[a]] -> [[a]] -> [[a]]
    a <<->> b = (zipWith . zipWith) (-) a b

    -- |a \<\<*\>\> b returns the product a*b of matrices
    (<<*>>) :: (Num a) => [[a]] -> [[a]] -> [[a]]
    a <<*>> b = [ [u <.> v | v <- L.transpose b] | u <- a]

    -- |k *\>\> m returns the product k*m of the scalar k and the matrix m
    (*>>) :: (Num a) => a -> [[a]] -> [[a]]
    k *>> m = (map . map) (k*) m

    -- |m \<\<*\> v is multiplication of a vector by a matrix on the left
    (<<*>) :: (Num a) => [[a]] -> [a] -> [a]
    m <<*> v = map (<.> v) m

    -- |v \<*\>\> m is multiplication of a vector by a matrix on the right
    (<*>>) :: (Num a) => [a] -> [[a]] -> [a]
    v <*>> m = map (v <.>) (L.transpose m)

    fMatrix n f = [[f i j | j <- [1..n]] | i <- [1..n]]

    -- version with indices from zero
    fMatrix' n f = [[f i j | j <- [0..n-1]] | i <- [0..n-1]]

    -- idMx n = fMatrix n (\i j -> if i == j then 1 else 0)

    idMx n = idMxs !! n where
        idMxs = map snd $ iterate next (0,[])
        next (j,m) = (j+1, (1 : replicate j 0) : map (0:) m)

    -- |iMx n is the n*n identity matrix
    iMx :: (Num t) => Int -> [[t]]
    iMx n = idMx n

    -- |jMx n is the n*n matrix of all 1s
    jMx :: (Num t) => Int -> [[t]]
    jMx n = replicate n (replicate n 1)

    -- |zMx n is the n*n matrix of all 0s
    zMx :: (Num t) => Int -> [[t]]
    zMx n = replicate n (replicate n 0)

    {-
    -- VECTORS

    data Vector d k = V [k] deriving (Eq,Ord,Show)

    instance (IntegerAsType d, Num k) => Num (Vector d k) where
        V a + V b = V $ a <+> b
        V a - V b = V $ a <-> b
        negate (V a) = V $ map negate a
        fromInteger 0 = V $ replicate d' 0 where d' = fromInteger $ value (undefined :: d)

    V v <>> M m = V $ v <*>> m

    M m <<> V v = V $ m <<*> v

    k |> V v = V $ k *> v
    -}

    -- MATRICES

    {-
    -- Square matrices of dimension d over field k
    data Matrix d k = M [[k]] deriving (Eq,Ord,Show)

    instance (IntegerAsType d, Num k) => Num (Matrix d k) where
        M a + M b = M $ a <<+>> b
        M a - M b = M $ a <<->> b
        negate (M a) = M $ (map . map) negate a
        M a * M b = M $ a <<*>> b
        fromInteger 0 = M $ zMx d' where d' = fromInteger $ value (undefined :: d)
        fromInteger 1 = M $ idMx d' where d' = fromInteger $ value (undefined :: d)

    instance (IntegerAsType d, Fractional a) => Fractional (Matrix d a) where
    	recip (M a) = case inverse a of
    		Nothing -> error "Matrix.recip: matrix is singular"
    		Just a' -> M a'
    -}

    -- |The inverse of a matrix (over a field), if it exists
    inverse :: (Eq a, Fractional a) => [[a]] -> Maybe [[a]]
    inverse m =
        let d = length m -- the dimension
            i = idMx d
            m' = zipWith (++) m i
            i1 = inverse1 m'
            i2 = inverse2 i1
        in if length i1 == d
           then Just i2
           else Nothing

    -- given (M|I), use row operations to get to (U|A), where U is upper triangular with 1s on diagonal
    inverse1 [] = []
    inverse1 ((x:xs):rs) =
        if x /= 0
        then let r' = (1/x) *> xs
             in (1:r') : inverse1 [ys <-> y *> r' | (y:ys) <- rs]
        else case filter (\r' -> head r' /= 0) rs of
             [] -> [] -- early termination, which will be detected in calling function
             r:_ -> inverse1 (((x:xs) <+> r) : rs)
    -- This is basically row echelon form

    -- given (U|A), use row operations to get to M^-1
    inverse2 [] = []
    inverse2 ((1:r):rs) = inverse2' r rs : inverse2 rs where
        inverse2' xs [] = xs
        inverse2' (x:xs) ((1:r):rs) = inverse2' (xs <-> x *> r) rs

    xs ! i = xs !! (i-1) -- ie, a 1-based list lookup instead of 0-based

    rowEchelonForm [] = []
    rowEchelonForm ((x:xs):rs) =
        if x /= 0
        then let r' = (1/x) *> xs
             in (1:r') : map (0:) (rowEchelonForm [ys <-> y *> r' | (y:ys) <- rs])
        else case filter (\r' -> head r' /= 0) rs of
             [] -> map (0:) (rowEchelonForm $ xs : map tail rs)
             r:_ -> rowEchelonForm (((x:xs) <+> r) : rs)
    rowEchelonForm zs@([]:_) = zs

    reducedRowEchelonForm :: (Eq a, Fractional a) => [[a]] -> [[a]]
    reducedRowEchelonForm m = reverse $ reduce $ reverse $ rowEchelonForm m where
        reduce (r:rs) = let r':rs' = reduceStep (r:rs) in r' : reduce rs' -- is this scanl or similar?
        reduce [] = []
        reduceStep ((1:xs):rs) = (1:xs) : [ 0: (ys <-> y *> xs) | y:ys <- rs]
        reduceStep rs@((0:_):_) = zipWith (:) (map head rs) (reduceStep $ map tail rs)
        reduceStep rs = rs

    -- Given a matrix m and (column) vector b, either find (column vector) x such that m x == b,
    -- or indicate that there is none
    solveLinearSystem m b =
        let augmented = zipWith (\r x -> r ++ [x]) m b -- augmented matrix
            trisystem = inverse1 augmented -- upper triangular form
            solution = reverse $ solveTriSystem $ reverse $ map reverse trisystem
        in if length solution == length b then Just solution else Nothing
        where solveTriSystem ([v,c]:rs) =
                  let x = v/c -- the first row tells us that cx == v
                      rs' = map (\(v':c':r) -> (v'-c'*x):r) rs
                  in x : solveTriSystem rs'
              solveTriSystem [] = []
              solveTriSystem _ = [] -- abnormal termination - m wasn't invertible

    isZero v = all (==0) v

    -- inSpanRE m v returns whether the vector v is in the span of the matrix m, where m is required to be in row echelon form
    inSpanRE ((1:xs):bs) (y:ys) = inSpanRE (map tail bs) (if y == 0 then ys else ys <-> y *> xs)
    inSpanRE ((0:xs):bs) (y:ys) = if y == 0 then inSpanRE (xs : map tail bs) ys else False
    inSpanRE _ ys = isZero ys

    rank m = length $ filter (not . isZero) $ rowEchelonForm m

    -- kernel of a matrix
    -- returns basis for vectors v s.t m <<*> v == 0
    kernel m = kernelRRE $ reducedRowEchelonForm m

    kernelRRE m =
        let nc = length $ head m -- the number of columns
            is = findLeadingCols 1 (L.transpose m) -- these are the indices of the columns which have a leading 1
            js = [1..nc] L.\\ is
            freeCols = let m' = take (length is) m -- discard zero rows
                       in zip is $ L.transpose [map (negate . (!j)) m' | j <- js]
            boundCols = zip js (idMx $ length js)
        in L.transpose $ map snd $ L.sort $ freeCols ++ boundCols
        where
        findLeadingCols i (c@(1:_):cs) = i : findLeadingCols (i+1) (map tail cs)
        findLeadingCols i (c@(0:_):cs) = findLeadingCols (i+1) cs
        findLeadingCols _ _ = []

    -- m ^- n = recip m ^ n

    -- t (M m) = M (L.transpose m)

    -- |The determinant of a matrix (over a field)
    det :: (Eq a, Fractional a) => [[a]] -> a
    det [[x]] = x
    det ((x:xs):rs) =
        if x /= 0
        then let r' = (1/x) *> xs
             in x * det [ys <-> y *> r' | (y:ys) <- rs]
        else case filter (\r' -> head r' /= 0) rs of
             [] -> 0
             r:_ -> det (((x:xs) <+> r) : rs)

    {-
    class IntegerAsType a where
        value :: a -> Integer

    data Z
    instance IntegerAsType Z where
        value _ = 0

    data S a
    instance IntegerAsType a => IntegerAsType (S a) where
        value _ = value (undefined :: a) + 1
    -}