1 {-# LANGUAGE ScopedTypeVariables #-}
3 -- | A Matrix type using Data.Vector as the underlying type. In other
4 -- words, the size is not fixed, but at least we have safe indexing if
7 -- This should be replaced with a fixed-size implementation eventually!
12 import qualified Data.Vector as V
14 type Rows a = V.Vector (V.Vector a)
15 type Columns a = V.Vector (V.Vector a)
16 data Matrix a = Matrix (Rows a) deriving Eq
19 (!) :: (Matrix a) -> (Int, Int) -> a
20 (Matrix rows) ! (i, j) = (rows V.! i) V.! j
23 (!?) :: (Matrix a) -> (Int, Int) -> Maybe a
24 (Matrix rows) !? (i, j) = do
29 -- | Unsafe indexing without bounds checking
30 unsafeIndex :: (Matrix a) -> (Int, Int) -> a
31 (Matrix rows) `unsafeIndex` (i, j) =
32 (rows `V.unsafeIndex` i) `V.unsafeIndex` j
34 -- | Return the @i@th column of @m@. Unsafe!
35 column :: (Matrix a) -> Int -> (V.Vector a)
36 column (Matrix rows) i =
37 V.fromList [row V.! i | row <- V.toList rows]
39 -- | The number of rows in the matrix.
40 nrows :: (Matrix a) -> Int
41 nrows (Matrix rows) = V.length rows
43 -- | The number of columns in the first row of the matrix.
44 ncols :: (Matrix a) -> Int
46 | V.length rows == 0 = 0
47 | otherwise = V.length (rows V.! 0)
49 -- | Return the vector of @m@'s columns.
50 columns :: (Matrix a) -> (Columns a)
52 V.fromList [column m i | i <- [0..(ncols m)-1]]
54 -- | Transose @m@; switch it's columns and its rows.
55 transpose :: (Matrix a) -> (Matrix a)
59 instance Show a => Show (Matrix a) where
61 concat $ V.toList $ V.map show_row rows
62 where show_row r = "[" ++ (show r) ++ "]\n"
64 instance Functor Matrix where
65 f `fmap` (Matrix rows) = Matrix (V.map (fmap f) rows)
69 vplus :: Num a => (V.Vector a) -> (V.Vector a) -> (V.Vector a)
70 vplus xs ys = V.zipWith (+) xs ys
72 -- | Vector subtraction.
73 vminus :: Num a => (V.Vector a) -> (V.Vector a) -> (V.Vector a)
74 vminus xs ys = V.zipWith (-) xs ys
76 -- | Add two vectors of rows.
77 rowsplus :: Num a => (Rows a) -> (Rows a) -> (Rows a)
79 V.zipWith vplus rs1 rs2
81 -- | Subtract two vectors of rows.
82 rowsminus :: Num a => (Rows a) -> (Rows a) -> (Rows a)
84 V.zipWith vminus rs1 rs2
86 -- | Matrix multiplication.
87 mtimes :: Num a => (Matrix a) -> (Matrix a) -> (Matrix a)
88 mtimes m1@(Matrix rows1) m2@(Matrix rows2) =
89 Matrix (V.fromList rows)
91 row i = V.fromList [ sum [ (m1 ! (i,k)) * (m2 ! (k,j)) | k <- [0..(ncols m1)-1] ]
92 | j <- [0..(ncols m2)-1] ]
93 rows = [row i | i <- [0..(nrows m1)-1]]
95 -- | Is @m@ symmetric?
96 symmetric :: Eq a => (Matrix a) -> Bool
100 -- | Construct a new matrix from a function @lambda@. The function
101 -- @lambda@ should take two parameters i,j corresponding to the
102 -- entries in the matrix. The i,j entry of the resulting matrix will
103 -- have the value returned by lambda i j.
105 -- The @imax@ and @jmax@ parameters determine the size of the matrix.
107 construct :: Int -> Int -> (Int -> Int -> a) -> (Matrix a)
108 construct imax jmax lambda =
111 row i = V.fromList [ lambda i j | j <- [0..jmax] ]
112 rows = V.fromList [ row i | i <- [0..imax] ]
114 -- | Given a positive-definite matrix @m@, computes the
115 -- upper-triangular matrix @r@ with (transpose r)*r == m and all
116 -- values on the diagonal of @r@ positive.
117 cholesky :: forall a. RealFloat a => (Matrix a) -> (Matrix a)
119 construct (nrows m - 1) (ncols m - 1) r
123 | i == j = sqrt(m ! (i,j) - sum [(r k i)**2 | k <- [0..i-1]])
124 | i < j = (((m ! (i,j)) - sum [(r k i)*(r k j) | k <- [0..i-1]]))/(r i i)
126 -- | It's not correct to use Num here, but I really don't want to have
127 -- to define my own addition and subtraction.
128 instance Num a => Num (Matrix a) where
129 -- Standard componentwise addition.
130 (Matrix rows1) + (Matrix rows2) = Matrix (rows1 `rowsplus` rows2)
132 -- Standard componentwise subtraction.
133 (Matrix rows1) - (Matrix rows2) = Matrix (rows1 `rowsminus` rows2)
135 -- Matrix multiplication.
136 m1 * m2 = m1 `mtimes` m2
138 abs _ = error "absolute value of matrices is undefined"
140 signum _ = error "signum of matrices is undefined"
142 fromInteger x = error "fromInteger of matrices is undefined"