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)
89 Matrix (V.fromList rows)
92 V.fromList [ sum [ (m1 ! (i,k)) * (m2 ! (k,j)) | k <- [0..(ncols m1)-1] ]
93 | j <- [0..(ncols m2)-1] ]
94 rows = [row i | i <- [0..(nrows m1)-1]]
96 -- | Is @m@ symmetric?
97 symmetric :: Eq a => (Matrix a) -> Bool
101 -- | Construct a new matrix from a function @lambda@. The function
102 -- @lambda@ should take two parameters i,j corresponding to the
103 -- entries in the matrix. The i,j entry of the resulting matrix will
104 -- have the value returned by lambda i j.
106 -- The @imax@ and @jmax@ parameters determine the size of the matrix.
108 construct :: Int -> Int -> (Int -> Int -> a) -> (Matrix a)
109 construct imax jmax lambda =
112 row i = V.fromList [ lambda i j | j <- [0..jmax] ]
113 rows = V.fromList [ row i | i <- [0..imax] ]
115 -- | Given a positive-definite matrix @m@, computes the
116 -- upper-triangular matrix @r@ with (transpose r)*r == m and all
117 -- values on the diagonal of @r@ positive.
118 cholesky :: forall a. RealFloat a => (Matrix a) -> (Matrix a)
120 construct (nrows m - 1) (ncols m - 1) r
123 r i j | i == j = sqrt(m ! (i,j) - sum [(r k i)**2 | k <- [0..i-1]])
125 (((m ! (i,j)) - sum [(r k i)*(r k j) | k <- [0..i-1]]))/(r i i)
128 -- | It's not correct to use Num here, but I really don't want to have
129 -- to define my own addition and subtraction.
130 instance Num a => Num (Matrix a) where
131 -- Standard componentwise addition.
132 (Matrix rows1) + (Matrix rows2) = Matrix (rows1 `rowsplus` rows2)
134 -- Standard componentwise subtraction.
135 (Matrix rows1) - (Matrix rows2) = Matrix (rows1 `rowsminus` rows2)
137 -- Matrix multiplication.
138 m1 * m2 = m1 `mtimes` m2
140 abs _ = error "absolute value of matrices is undefined"
142 signum _ = error "signum of matrices is undefined"
144 fromInteger _ = error "fromInteger of matrices is undefined"