Unclassified
LS,
here an haskell file that implements overloaded matrix operations. Such
overloading is not found
in Jan Skibinski's Orthogonals.lhs program. I would like to make this
code public domain
and add this file to the "hslibs", does anyone know how to go about it?
Greetings Frank Mulder
-- overloaded matrix operations
--
-- by F.W. Mulder
-- fwmulder(a)signaal.nl
-- Hollandse Signaalapparaten B.V.
--
module Matrix (
t,Mat(M,V,VT,S),Mtx((./),(.*),(.~),det,fromScalar,fromVec,
fromMat,vecLength,matDims,matShow,solve,lu,ldl,zeros,
ones,diag,trace,distL,isCovMat),
SM,sm_find, sm_mk_spm, sm_mk_mat, sm_size,sm_rows,sm_cols)
where
import List
import Maybe
-------------------- indexless linear algebra dense matrices
------------------
t = -1000
infixr 9 .~
data Mat = M [[Double]] | V [Double] | VT [Double] | S Double
deriving (Eq,Show,Read)
class Mtx a where
(./) :: a -> Double -> a
(.*) :: a -> Double -> a
(.~) :: a -> Integer -> a
det :: a -> Double
fromScalar :: a -> Double
fromVec :: a -> [Double]
fromMat :: a -> [[Double]]
vecLength :: a -> Int
matDims :: a -> (Int,Int)
matShow :: a -> String
solve :: (a -> (a,a)) -> a -> a -> a
lu :: a -> (a,a)
ldl :: a -> (a,a)
zeros :: Int -> a
ones :: Int -> a
diag :: a -> a
trace :: a -> Double
distL :: a -> a -> Double
isCovMat :: a -> Bool
instance Mtx Mat where
(./) (M a) b = M [[e / b | e <- r] | r <- a]
(./) (V a) b = V [e / b | e <- a]
(.*) (M a) b = M [[e * b | e <- r] | r <- a]
(.*) (V a) b = V [e * b | e <- a]
(.~) (M a) n | n == -1 = M (inv' a)
| n == t = M (transpose a)
| otherwise = error "Matrix:inv"
(.~) (V a) n | n == t = VT a
| otherwise = error "Matrix:inv"
(.~) (VT a) n | n == t = V a
| otherwise = error "Matrix:inv"
det (M a) = det' a
det _ = error "Matrix:det"
fromScalar (S a) = a
fromVec (V a) = a
fromVec (VT a) = a
fromMat (M a) = a
vecLength (V x) = length x
matDims (M a) = (length a, length (head a))
matDims _ = error "Matrix:matDims"
matShow (M b) = "\n\n" ++ concat [show r ++ "\n" | r <- b] ++ ";"
solve fac m = solve' fac m
lu m = lu' m
zeros n = V (take n (repeat 0.0))
ones n = V (take n (repeat 1.0))
diag (M a) = V (diag' a)
trace (M a) = sum (diag' a)
distL a b = distL' a b
isCovMat a = isCovMat' a
instance Num a => Num [a] where
(+) a b = zipWith (+) a b
(-) a b = zipWith (-) a b
instance Num Mat where
(+) (M a) (M b) = M (a + b)
(+) (V a) (V b) = V (a + b)
(+) (VT a) (VT b) = VT (a + b)
(-) (M a) (M b) = M (a - b)
(-) (V a) (V b) = V (a - b)
(-) (VT a) (VT b) = VT (a - b)
(*) (M a) (M b) = M (a `xmm` b)
(*) (M a) (V b) = V (a `xmv` b)
(*) (VT a) (V b) = S (inprod a b)
(*) (V a) (VT b) = M (outprod a b)
(*) (VT a) (M b) = VT (a `xvm` b)
diag' [r] = [head r]
diag' ([]:_) = error "Matrix:diag'"
diag' (r:rs) = head r : diag' [tail r | r <- rs]
dims' :: [[Double]] -> (Int,Int)
dims' [r] = (length r, 0)
dims' (r:rs) = (length r, length rs + 1)
isCovMat' b = n == m && all (\e -> e > 0.0) (fromVec (diag b))
where (n,m) = matDims b
------------------- standard determinand, inverse
-----------------------------
det' [[x]] = x
det' a@(v:_) = inprod v [cofactor a 0 j | j <- [0..n-1]]
where n = length v
cofactor a i j = fromIntegral ((-1)^(i+j)) * det' (minor a i j)
minor a i j = [omit j v | v <- omit i a]
omit 0 (x:xs) = xs
omit _ [] = []
omit n (x:xs) = (x:omit (n-1) xs)
inv' [[x]] = [[1.0/x]]
inv' a = [[(cofactor a j i) / d | j <- [0..n]] | i <- [0..n]]
where n = length a - 1
d = det' a
----------------------- vector products
---------------------------------------
inprod [] _ = 0.0
inprod (x:xs) (y:ys) = x * y + inprod xs ys
outprod a b = [[e * f | f <- b] | e <- a]
----------------------- matrix products
---------------------------------------
xmm a b = [[inprod r c | c <- transpose b] | r <- a]
xmv a b = [inprod r b | r <- a]
xvm a b = [inprod a c | c <- b]
----------------------- LU decomposition
--------------------------------------
lu' (M a) = factLU a [head a] []
factLU [[a]] u l = (M (transpose (l ++ [[1.0]])), M u)
factLU a u l = factLU a' (u ++ [head a']) (l ++ [[1.0] ++ cm])
where (c,cs) = unzip [(head r, tail r) | r <- a]
cm = [e / (head c) | e <- (tail c)]
rm = head cs
a' = tail cs --- \\ outprod cm rm
----------------------- choleski decomposition
--------------------------------
ldl' (M (r:rs)) = (M rt, M (transpose rt))
where rt = ldl2 1 ([],r,rs)
ldl2 j (u,r,[]) = [factCH r u j 1 []]
ldl2 j (u,r,l) = (ro : ldl2 (j+1) (u++[ro],head l, tail l))
where ro = factCH r u j 1 []
factCH [] _ _ _ _ = []
factCH (a:as) u j i lj | i > j = []
| i == j = (el1 : factCH as u j (i+1) (lj++[el1]))
| otherwise = (el2 : factCH as u j (i+1) (lj++[el2]))
where
el1 = sqrt (a - inprod lj lj)
el2 = (a - inprod li lj) / lii
(li,(lii:_)) = splitAt (i-1) (u!!(i-1))
-- forward and backward substitution for lower and upper triang.
matrices -----
forwL l b = subst_ l b []
backU u b = reverse (subst_ (map reverse (reverse u)) (reverse b) [])
subst_ [] _ s = s
subst_ (l:ls) (b:bs) s = subst_ ls bs (s ++ [xi])
where xi = (b - inprod (init l) s) / (last l)
solve' fac a (V b) = let (M l,M u) = fac a in V (backU u (forwL l b))
distL' (V x) (M r) = inprod x (backU (transpose r) (forwL r x))
-- Matrix.hs
-- overloaded matrix operations
--
-- by F.W. Mulder
-- fwmulder(a)signaal.nl
-- Hollandse Signaalapparaten B.V.
--
module Matrix ( t,Mat(M,V,VT,S),Mtx((./),(.*),(.~),det,fromScalar,fromVec,
fromMat,vecLength,matDims,matShow,solve,lu,ldl,zeros,
ones,diag,trace,distL,isCovMat),
SM,sm_find, sm_mk_spm, sm_mk_mat, sm_size,sm_rows,sm_cols)
where
import List
import Maybe
-------------------- indexless linear algebra dense matrices ------------------
t = -1000
infixr 9 .~
data Mat = M [[Double]] | V [Double] | VT [Double] | S Double
deriving (Eq,Show,Read)
class Mtx a where
(./) :: a -> Double -> a
(.*) :: a -> Double -> a
(.~) :: a -> Integer -> a
det :: a -> Double
fromScalar :: a -> Double
fromVec :: a -> [Double]
fromMat :: a -> [[Double]]
vecLength :: a -> Int
matDims :: a -> (Int,Int)
matShow :: a -> String
solve :: (a -> (a,a)) -> a -> a -> a
lu :: a -> (a,a)
ldl :: a -> (a,a)
zeros :: Int -> a
ones :: Int -> a
diag :: a -> a
trace :: a -> Double
distL :: a -> a -> Double
isCovMat :: a -> Bool
instance Mtx Mat where
(./) (M a) b = M [[e / b | e <- r] | r <- a]
(./) (V a) b = V [e / b | e <- a]
(.*) (M a) b = M [[e * b | e <- r] | r <- a]
(.*) (V a) b = V [e * b | e <- a]
(.~) (M a) n | n == -1 = M (inv' a)
| n == t = M (transpose a)
| otherwise = error "Matrix:inv"
(.~) (V a) n | n == t = VT a
| otherwise = error "Matrix:inv"
(.~) (VT a) n | n == t = V a
| otherwise = error "Matrix:inv"
det (M a) = det' a
det _ = error "Matrix:det"
fromScalar (S a) = a
fromVec (V a) = a
fromVec (VT a) = a
fromMat (M a) = a
vecLength (V x) = length x
matDims (M a) = (length a, length (head a))
matDims _ = error "Matrix:matDims"
matShow (M b) = "\n\n" ++ concat [show r ++ "\n" | r <- b] ++ ";"
solve fac m = solve' fac m
lu m = lu' m
zeros n = V (take n (repeat 0.0))
ones n = V (take n (repeat 1.0))
diag (M a) = V (diag' a)
trace (M a) = sum (diag' a)
distL a b = distL' a b
isCovMat a = isCovMat' a
instance Num a => Num [a] where
(+) a b = zipWith (+) a b
(-) a b = zipWith (-) a b
instance Num Mat where
(+) (M a) (M b) = M (a + b)
(+) (V a) (V b) = V (a + b)
(+) (VT a) (VT b) = VT (a + b)
(-) (M a) (M b) = M (a - b)
(-) (V a) (V b) = V (a - b)
(-) (VT a) (VT b) = VT (a - b)
(*) (M a) (M b) = M (a `xmm` b)
(*) (M a) (V b) = V (a `xmv` b)
(*) (VT a) (V b) = S (inprod a b)
(*) (V a) (VT b) = M (outprod a b)
(*) (VT a) (M b) = VT (a `xvm` b)
diag' [r] = [head r]
diag' ([]:_) = error "Matrix:diag'"
diag' (r:rs) = head r : diag' [tail r | r <- rs]
dims' :: [[Double]] -> (Int,Int)
dims' [r] = (length r, 0)
dims' (r:rs) = (length r, length rs + 1)
isCovMat' b = n == m && all (\e -> e > 0.0) (fromVec (diag b))
where (n,m) = matDims b
------------------- standard determinand, inverse -----------------------------
det' [[x]] = x
det' a@(v:_) = inprod v [cofactor a 0 j | j <- [0..n-1]]
where n = length v
cofactor a i j = fromIntegral ((-1)^(i+j)) * det' (minor a i j)
minor a i j = [omit j v | v <- omit i a]
omit 0 (x:xs) = xs
omit _ [] = []
omit n (x:xs) = (x:omit (n-1) xs)
inv' [[x]] = [[1.0/x]]
inv' a = [[(cofactor a j i) / d | j <- [0..n]] | i <- [0..n]]
where n = length a - 1
d = det' a
----------------------- vector products ---------------------------------------
inprod [] _ = 0.0
inprod (x:xs) (y:ys) = x * y + inprod xs ys
outprod a b = [[e * f | f <- b] | e <- a]
----------------------- matrix products ---------------------------------------
xmm a b = [[inprod r c | c <- transpose b] | r <- a]
xmv a b = [inprod r b | r <- a]
xvm a b = [inprod a c | c <- b]
----------------------- LU decomposition --------------------------------------
lu' (M a) = factLU a [head a] []
factLU [[a]] u l = (M (transpose (l ++ [[1.0]])), M u)
factLU a u l = factLU a' (u ++ [head a']) (l ++ [[1.0] ++ cm])
where (c,cs) = unzip [(head r, tail r) | r <- a]
cm = [e / (head c) | e <- (tail c)]
rm = head cs
a' = tail cs --- \\ outprod cm rm
----------------------- choleski decomposition --------------------------------
ldl' (M (r:rs)) = (M rt, M (transpose rt))
where rt = ldl2 1 ([],r,rs)
ldl2 j (u,r,[]) = [factCH r u j 1 []]
ldl2 j (u,r,l) = (ro : ldl2 (j+1) (u++[ro],head l, tail l))
where ro = factCH r u j 1 []
factCH [] _ _ _ _ = []
factCH (a:as) u j i lj | i > j = []
| i == j = (el1 : factCH as u j (i+1) (lj++[el1]))
| otherwise = (el2 : factCH as u j (i+1) (lj++[el2]))
where
el1 = sqrt (a - inprod lj lj)
el2 = (a - inprod li lj) / lii
(li,(lii:_)) = splitAt (i-1) (u!!(i-1))
-- forward and backward substitution for lower and upper triang. matrices -----
forwL l b = subst_ l b []
backU u b = reverse (subst_ (map reverse (reverse u)) (reverse b) [])
subst_ [] _ s = s
subst_ (l:ls) (b:bs) s = subst_ ls bs (s ++ [xi])
where xi = (b - inprod (init l) s) / (last l)
solve' fac a (V b) = let (M l,M u) = fac a in V (backU u (forwL l b))
distL' (V x) (M r) = inprod x (backU (transpose r) (forwL r x))
---------------------- sparse matrices ----------------------------------------
type SM = ((Int,Int),Double)
sm_size [[]] = (0,0)
sm_size (r:rs) = (length rs + 1, length r)
sm_find ix sm = if (cc /= Nothing) then (Just (ix,fromJust cc)) else Nothing
where cc = lookup ix sm
sm_mk_spm b = ([((i,j),e) | (i,r) <- zip [1..n] b, (j,e) <- zip [1..m] r],n,m)
where (n,m) = sm_size b
sm_mk_mat (sm,n,m) large =
[[mkc (sm_find (i,j) sm) | j <- [1..m]] | i <- [1..n]]
where mkc Nothing = large
mkc (Just (_,c)) = c
sm_rows [] = []
sm_rows [e] = [[e]]
sm_rows (e@((i1,_),_):es) = (e : rw) : sm_rows rws
where (rw,rws) = partition (\((i,_),_) -> i == i1) es
sm_cols [] = []
sm_cols [e] = [[e]]
sm_cols (e@((_,j1),_):es) = (e : cl) : sm_cols cls
where (cl,cls) = partition (\((_,j),_) -> j == j1) es
