Matrix.hs (with Vectors instead of Lists, and Either monad)

-- Matrix.hs (with Vectors instead of Lists, and Either monad)
module Matrix where

import qualified Data.Vector as BV
import qualified Data.Vector.Unboxed as UV
import qualified System.Random as R
import qualified Control.Monad as CM

isMatrixValid matrix =
 if BV.length lengths == 0 then
  False
 else
  (BV.all (== BV.head lengths) lengths)
 where
  lengths = BV.map UV.length matrix

newMatrix matrix =
 if isMatrixValid matrix then Right matrix else Left "Matrix is not valid."

fromListOfLists ll = newMatrix $ BV.fromList (map UV.fromList ll)

boxedToUnboxedVector boxed = UV.fromList $ BV.toList boxed

getMatrixSize matrix =
 let
  rows = BV.length matrix
 in
  if isMatrixValid matrix then
   Right $ if rows == 0 then (0, 0) else (rows, UV.length $ BV.head matrix)
  else
   Left "Matrix cannot be measured because it is invalid."

at matrix (row,column) = (matrix BV.! row) UV.! column

generateMatrixI (rows,columns) f =
 BV.generate rows (generateRow columns)
 where
  generateRow columns row =
   UV.generate columns (\column -> f (row,column))

transpose matrix =
 getMatrixSize matrix >>= \(rows,columns) ->
 Right $ generateMatrixI (columns,rows) (\(column,row) -> at matrix (row,column))

newRowVector vector =
 let rowVector = BV.singleton vector
 in newMatrix rowVector

newColumnVector vector =
 Right $ generateMatrixI (UV.length vector,1) (\(r,c) -> vector UV.! r)

mapMatrix f matrix =
 newMatrix matrix >>= \matrix0 ->
 Right $ BV.map (UV.map f) matrix0

zipMatrixesWith f matrix0 matrix1 =
 getMatrixSize matrix0 >>= \size0 ->
 getMatrixSize matrix1 >>= \size1 ->
 if size0 /= size1 then
  Left "Matrixes cannot be zipped because their sizes don't match."
 else
  Right $ BV.zipWith (UV.zipWith f) matrix0 matrix1

mapMatrixWithIndexes f matrix =
 getMatrixSize matrix >>= \(rows,columns) ->
 Right $ generateMatrixI (rows,columns) (\p -> f p (at matrix p))

zipMatrixesWithIndexes f matrix0 matrix1 =
 getMatrixSize matrix0 >>= \size0 ->
 getMatrixSize matrix1 >>= \size1 ->
 if size0 /= size1 then
  Left "Matrixes cannot be zipped with indexes because their sizes don't match."
 else
  let (rows,columns) = size0
  in Right $ generateMatrixI size0 (\p -> f p (at matrix0 p) (at matrix1 p))

sumMatrixes matrix0 matrix1 = zipMatrixesWith (+) matrix0 matrix1

subtractMatrixes matrix0 matrix1 = zipMatrixesWith (-) matrix0 matrix1

scaleMatrix value matrix = mapMatrix (* value) matrix

multiplyMatrixes matrix0 matrix1 =
 getMatrixSize matrix0 >>= \size0 ->
 getMatrixSize matrix1 >>= \size1 ->
 let size = (fst size0,snd size1)
 in
  if snd size0 /= fst size1 then
   Left "Matrixes cannot be multiplied because their dimensions don't match."
  else
   transpose matrix1 >>= \tMatrix1 ->
   let
    pointProduct v0 v1 = UV.sum $ UV.zipWith (*) v0 v1
    generateValue (row,column) =
     pointProduct (matrix0 BV.! row) (tMatrix1 BV.! column)
   in Right $ generateMatrixI size (\p -> generateValue p)

newMatrixWithValue size value = generateMatrixI size (\p -> value)

pivot matrix (i,j) =
 getMatrixSize matrix >>= \(rows,columns) ->
 Right $ generateMatrixI (rows-1,columns-1) generateValue
 where
  generateValue (row,column) =
   at matrix (r,c)
   where
    r = if row < i then row else row + 1
    c = if column < j then column else column + 1

cofactorSign (i,j) = if (mod (i + j) 2) == 0 then 1.0 else -1.0

determinant matrix =
 getMatrixSize matrix >>= \(rows,columns) ->
 if rows /= columns then
  Left "Only square matrixes have determinant."
 else
  if rows == 1 then
   Right $ at matrix (0,0)
  else
   if rows == 2 then
    Right $ determinant2x2 matrix
   else
    determinantNxN rows matrix
 where
  determinant2x2 matrix =
   (at matrix (0,0)) * (at matrix (1,1)) - (at matrix (0,1)) * (at matrix (1,0))

determinantNxN n matrix =
 CM.foldM accum 0 (map detTerm [0..n-1])
 where
  accum acc term =
   case term of
    Left error -> Left error
    Right value -> Right (acc + value)
  detTerm j =
   pivot matrix p >>= \detPivot ->
   determinant detPivot >>= \det ->
   Right ((cofactorSign p) * (at matrix p) * det)
   where p = (0,j)

cofactorsMatrix matrix =
 getMatrixSize matrix >>= \(rows,columns) ->
 if rows /= columns then
  Left "Only square matrixes have cofactors matrix."
 else
  Right $ generateMatrixI (rows,rows) generateValue
 where
  generateValue p =
   case (pivot matrix p >>= \piv -> determinant piv) of
    Left error -> 0
    Right detValue -> (cofactorSign p) * detValue

inverse matrix =
 cofactorsMatrix matrix >>= \cofactors ->
 determinant matrix >>= \det ->
 if det == 0.0 then
  Left "This matrix cannot be inverted because its determinant is zero."
 else
  transpose cofactors >>= \transposedCofactors ->
  scaleMatrix (1.0/det) transposedCofactors

identity n one =
 generateMatrixI (n,n) (\(i,j) -> if i == j then one else zero)
 where zero = one - one

areNumbersAlmostEqual x0 x1 precision =
 (abs (x0 - x1)) <= precision

areMatrixesAlmostEqual matrix0 matrix1 precision =
 getMatrixSize matrix0 >>= \size0 ->
 getMatrixSize matrix1 >>= \size1 ->
 if size0 /= size1 then
  Left "The matrixes cannot be compared because their dimensions don't match."
 else
  let
   booleans = generateMatrixI size0 generateValue
   generateValue p = areNumbersAlmostEqual (at matrix0 p) (at matrix1 p) precision
  in Right $ BV.all (\row -> UV.all id row) booleans