-- Version of 10/28/06
module Matrix where

import List

type Vector = [Int]
type Matrix = [Vector]

--basic constructions for vectors

zeroVector :: Int -> Vector
zeroVector n = replicate n 0

--basic operations for vectors

dotProduct :: Vector -> Vector -> Int
dotProduct v w = sum ( zipWith (*) v w )

vectorSum :: Vector -> Vector -> Vector
vectorSum = zipWith (+)

vectorScalarProduct :: Int -> Vector -> Vector
vectorScalarProduct n vec = [ n * x | x <- vec ]

--basic constructions for matrices

elemMatrix :: Int -> Int -> Int -> Int -> Matrix
-- elemMatrix n i j v   is the n-by-n elementary matrix with 
-- entry  v  in the (i,j) place
elemMatrix n i j v
  = [ [ entry row column | column <- [1..n] ] | row <- [1..n] ]
  where
  entry x y
    | x == y           = 1
    | x == i && y == j = v
    | otherwise        = 0

idMatrix :: Int -> Matrix
idMatrix n = elemMatrix n 1 1 1

zeroMatrix :: Int -> Int -> Matrix
zeroMatrix i j = replicate i (zeroVector j)

--basic operations for matrices

matrixSum :: Matrix -> Matrix -> Matrix
matrixSum = zipWith vectorSum

matrixScalarProduct :: Int -> Matrix -> Matrix
matrixScalarProduct n m = [ vectorScalarProduct n row | row <- m ]

matrixProduct :: Matrix -> Matrix -> Matrix
matrixProduct m n = [ map (dotProduct r) (transpose n) | r <- m ]

{- The determinant and inverse functions given here are only for examples
of Haskell syntax.  Efficient versions using row operations are implemented
in RowOperations.hs .

--determinant using cofactors

remove :: Matrix -> Int -> Int -> Matrix
remove m i j  
  | m == [] || i < 1 || i > numRows m || j < 1 || j > numColumns m
    = error "(i,j) out of range"
  | otherwise = transpose ( cut (transpose ( cut m i ) ) j )

determinant :: Matrix -> Int
determinant [] = error "determinant: 0-by-0 matrix"
determinant [[n]] = n
determinant m = sum [ (-1)^(j+1) * (head m)!!(j-1) * determinant (remove m 1 j) | 
  j <- [1..(numColumns m) ] ]

--inverse

cofactor :: Matrix -> Int -> Int -> Int
cofactor m i j = (-1)^(i+j) * determinant (remove m i j)

cofactorMatrix :: Matrix -> Matrix
cofactorMatrix m = [ [ (cofactor m i j) | j <- [1..n] ] | i <- [1..n] ]
  where
  n = length m

inverse :: Matrix -> Matrix
inverse m = transpose [ [ quot x ( determinant m) | 
  x <- row ] | row <- (cofactorMatrix m) ]

-}

----------------------------------------------------------
-- matrix utilities

numRows :: Matrix -> Int
numRows = length

numColumns :: Matrix -> Int
numColumns = length . head

----------------------------------------------------------
-- other utilities

cut :: [a] -> Int -> [a]
cut [] n = []
cut xs n
  | n < 1 || n > (length xs) = xs
  | otherwise = (take (n-1) xs) ++ drop n xs

showMatrix :: Matrix -> IO()
showMatrix m = putStr (rowsWithEndlines m)
  where
  rowsWithEndlines m = concat ( map (\x -> (show x) ++ "\n") m )

----------------------------------------------------------
-- test data

mat1, mat2, mat3, mat4, mat5, mat6, mat7, mat8, mat9, mat10, 
  mat11, mat12, mat13 :: Matrix
mat1 = [ [1,2,3], [4,5,6], [7,8,9] ]
mat2 = [ [1,2,3], [4,5,6], [7,8,10] ]
mat3 = [ [1,2,3, 0], [-6, 4,5,6], [7,8,10, 10], [-3, -2, -1, 4] ]
mat4 = [ [1,2], [1,1]]
mat5 = [ [1, 1, 0], [0, 1, 0], [0, 0, 1] ]
mat6 = [[1,2,0, 3,0],[-6,4,0, 5,6],[7,8,72, 10,10],[-3,-2,-1,4, -103]]
mat7 = [[2,0, 3,0],[-4,0, 5,6],[8,72, 10,10],[-2,-1,4, -103],[8,9,0,-14]]
mat8 =[ [0,0,0], [0,0,0], [0,0,0]]
mat9 = [ [2, 0], [0, 1] ]
mat10 = [ [1, 0], [ 0, 2] ]
mat11 = [ [2, 0, 0], [0, 3, 0], [0, 0, 1] ]
mat12 = [[0,0,0,0,0],[1,1,-1,-1,0],[0,0,1,1,-1],[0,1,-1,-1,0],[1,0,0,1,-1]]
mat13 = [[-1,0,0,-1,0,0,0,-1,0,0,0,-1],[1,0,1,1,1,1,1,1,1,1,1,1],[-1,0,-1,0,-1,0,0,0,-1,0,0,0],[0,-1,0,-1,0,-1,0,0,0,-1,0,0],[0,0,0,1,-1,0,0,0,0,0,0,0],[-1,-1,-1,-1,0,-1,0,0,0,0,0,0],[1,0,0,0,0,0,-1,0,0,0,0,0],[0,1,0,0,0,0,0,-1,0,0,0,0],[0,0,0,0,0,0,0,1,-1,0,0,0],[0,0,0,0,-1,-1,-1,-1,0,-1,0,0],[0,0,0,0,1,0,0,0,0,0,-1,0],[0,0,0,0,0,1,0,0,0,0,0,-1]]
----------------------------------------------------------
{-
-- test of the inverse routine

eProduct :: Int -> [(Int,Int,Int)] -> Matrix
-- eProduct n [(Int,Int,Int)] is the n-by-n matrix which is a product
-- of the elementary n-by-n matrices given by the list triples
eProduct n [] = idMatrix n
eProduct n ((i,j,v):rest) = matrixProduct ( elemMatrix n i j v) (eProduct n rest)

minSize :: [(Int,Int,Int)] -> Int
-- smallest size of matrix for which all elementary matrices are defined
minSize list = maximum (concat [ [i,j] | (i,j,value) <- list ] )

checkInverse :: [(Int,Int,Int)] -> String
checkInverse list =
  "\n   M       = " ++ (show m) ++ 
  "\nInverse(M) = " ++ (show (inverse m)) ++ 
  if matrixProduct m (inverse m) == idMatrix n then "\nOK.\n" else "\nError.\n"
  where
  m = eProduct n list
  n = minSize list

list1 :: [(Int,Int,Int)]
list1 =  [(1,2,1), (1,3,-1), (1,2,1), (3,2,-2), (3,1,-3)]

list2 :: [(Int,Int,Int)]
list2 =  [(1,2,4), (4,2,-1), (4,1,-2), (4,1,1), (1,3,-3), (2,3,2), 
  (1,2,2), (1,4,-3), (1,3,-1), (3,2,-1), (3,1,-1)]

test :: IO()
test = putStr ( 
  ( checkInverse list1 ) ++ 
  ( checkInverse list2 )
  )

----------------------------------------------------------

--Here is the output of the test routine:

Main> test

[[1,-2,1],[0,1,0],[3,-4,4]]
[[1,-2,1],[0,1,0],[3,-4,4]]
Equal.


[[4,-15,8,3],[0,1,-2,0],[4,-14,7,3],[1,-3,0,1]]
[[4,-15,8,3],[0,1,-2,0],[4,-14,7,3],[1,-3,0,1]]
Equal.

-}