1 {-# LANGUAGE ExistentialQuantification #-}
2 {-# LANGUAGE FlexibleContexts #-}
3 {-# LANGUAGE FlexibleInstances #-}
4 {-# LANGUAGE MultiParamTypeClasses #-}
5 {-# LANGUAGE ScopedTypeVariables #-}
6 {-# LANGUAGE TypeFamilies #-}
7 {-# LANGUAGE RebindableSyntax #-}
9 -- | Boxed matrices; that is, boxed m-vectors of boxed n-vectors. We
10 -- assume that the underlying representation is
11 -- Data.Vector.Fixed.Boxed.Vec for simplicity. It was tried in
12 -- generality and failed.
17 import Data.List (intercalate)
19 import Data.Vector.Fixed (
33 import qualified Data.Vector.Fixed as V (
44 import Data.Vector.Fixed.Boxed (Vec)
45 import Data.Vector.Fixed.Internal (Arity, arity)
49 import NumericPrelude hiding ((*), abs)
50 import qualified NumericPrelude as NP ((*))
51 import qualified Algebra.Algebraic as Algebraic
52 import Algebra.Algebraic (root)
53 import qualified Algebra.Additive as Additive
54 import qualified Algebra.Ring as Ring
55 import qualified Algebra.Module as Module
56 import qualified Algebra.RealRing as RealRing
57 import qualified Algebra.ToRational as ToRational
58 import qualified Algebra.Transcendental as Transcendental
59 import qualified Prelude as P
61 data Mat m n a = (Arity m, Arity n) => Mat (Vec m (Vec n a))
62 type Mat1 a = Mat N1 N1 a
63 type Mat2 a = Mat N2 N2 a
64 type Mat3 a = Mat N3 N3 a
65 type Mat4 a = Mat N4 N4 a
66 type Mat5 a = Mat N5 N5 a
68 instance (Eq a) => Eq (Mat m n a) where
69 -- | Compare a row at a time.
73 -- >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int
74 -- >>> let m2 = fromList [[1,2],[3,4]] :: Mat2 Int
75 -- >>> let m3 = fromList [[5,6],[7,8]] :: Mat2 Int
81 (Mat rows1) == (Mat rows2) =
82 V.and $ V.zipWith comp rows1 rows2
84 -- Compare a row, one column at a time.
85 comp row1 row2 = V.and (V.zipWith (==) row1 row2)
88 instance (Show a) => Show (Mat m n a) where
89 -- | Display matrices and vectors as ordinary tuples. This is poor
90 -- practice, but these results are primarily displayed
91 -- interactively and convenience trumps correctness (said the guy
92 -- who insists his vector lengths be statically checked at
97 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
102 "(" ++ (intercalate "," (V.toList row_strings)) ++ ")"
104 row_strings = V.map show_vector rows
106 "(" ++ (intercalate "," element_strings) ++ ")"
109 element_strings = P.map show v1l
112 -- | Convert a matrix to a nested list.
113 toList :: Mat m n a -> [[a]]
114 toList (Mat rows) = map V.toList (V.toList rows)
116 -- | Create a matrix from a nested list.
117 fromList :: (Arity m, Arity n) => [[a]] -> Mat m n a
118 fromList vs = Mat (V.fromList $ map V.fromList vs)
121 -- | Unsafe indexing.
122 (!!!) :: (Arity m, Arity n) => Mat m n a -> (Int, Int) -> a
123 (!!!) m (i, j) = (row m i) ! j
126 (!!?) :: Mat m n a -> (Int, Int) -> Maybe a
127 (!!?) m@(Mat rows) (i, j)
128 | i < 0 || j < 0 = Nothing
129 | i > V.length rows = Nothing
130 | otherwise = if j > V.length (row m j)
132 else Just $ (row m j) ! j
135 -- | The number of rows in the matrix.
136 nrows :: forall m n a. (Arity m) => Mat m n a -> Int
137 nrows _ = arity (undefined :: m)
139 -- | The number of columns in the first row of the
140 -- matrix. Implementation stolen from Data.Vector.Fixed.length.
141 ncols :: forall m n a. (Arity n) => Mat m n a -> Int
142 ncols _ = arity (undefined :: n)
145 -- | Return the @i@th row of @m@. Unsafe.
146 row :: Mat m n a -> Int -> (Vec n a)
147 row (Mat rows) i = rows ! i
150 -- | Return the @j@th column of @m@. Unsafe.
151 column :: Mat m n a -> Int -> (Vec m a)
152 column (Mat rows) j =
153 V.map (element j) rows
160 -- | Transpose @m@; switch it's columns and its rows. This is a dirty
161 -- implementation.. it would be a little cleaner to use imap, but it
162 -- doesn't seem to work.
164 -- TODO: Don't cheat with fromList.
168 -- >>> let m = fromList [[1,2], [3,4]] :: Mat2 Int
172 transpose :: (Arity m, Arity n) => Mat m n a -> Mat n m a
173 transpose m = Mat $ V.fromList column_list
175 column_list = [ column m i | i <- [0..(ncols m)-1] ]
178 -- | Is @m@ symmetric?
182 -- >>> let m1 = fromList [[1,2], [2,1]] :: Mat2 Int
186 -- >>> let m2 = fromList [[1,2], [3,1]] :: Mat2 Int
190 symmetric :: (Eq a, Arity m) => Mat m m a -> Bool
195 -- | Construct a new matrix from a function @lambda@. The function
196 -- @lambda@ should take two parameters i,j corresponding to the
197 -- entries in the matrix. The i,j entry of the resulting matrix will
198 -- have the value returned by lambda i j.
200 -- TODO: Don't cheat with fromList.
204 -- >>> let lambda i j = i + j
205 -- >>> construct lambda :: Mat3 Int
206 -- ((0,1,2),(1,2,3),(2,3,4))
208 construct :: forall m n a. (Arity m, Arity n)
209 => (Int -> Int -> a) -> Mat m n a
210 construct lambda = Mat rows
212 -- The arity trick is used in Data.Vector.Fixed.length.
213 imax = (arity (undefined :: m)) - 1
214 jmax = (arity (undefined :: n)) - 1
215 row' i = V.fromList [ lambda i j | j <- [0..jmax] ]
216 rows = V.fromList [ row' i | i <- [0..imax] ]
219 -- | Given a positive-definite matrix @m@, computes the
220 -- upper-triangular matrix @r@ with (transpose r)*r == m and all
221 -- values on the diagonal of @r@ positive.
225 -- >>> let m1 = fromList [[20,-1], [-1,20]] :: Mat2 Double
227 -- ((4.47213595499958,-0.22360679774997896),(0.0,4.466542286825459))
228 -- >>> (transpose (cholesky m1)) * (cholesky m1)
229 -- ((20.000000000000004,-1.0),(-1.0,20.0))
231 cholesky :: forall m n a. (Algebraic.C a, Arity m, Arity n)
232 => (Mat m n a) -> (Mat m n a)
233 cholesky m = construct r
236 r i j | i == j = sqrt(m !!! (i,j) - sum [(r k i)^2 | k <- [0..i-1]])
238 (((m !!! (i,j)) - sum [(r k i) NP.* (r k j) | k <- [0..i-1]]))/(r i i)
242 -- | Returns True if the given matrix is upper-triangular, and False
247 -- >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
248 -- >>> is_upper_triangular m
251 -- >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
252 -- >>> is_upper_triangular m
255 is_upper_triangular :: (Eq a, Ring.C a, Arity m, Arity n)
257 is_upper_triangular m =
260 results = [[ test i j | i <- [0..(nrows m)-1]] | j <- [0..(ncols m)-1] ]
262 test :: Int -> Int -> Bool
265 | otherwise = m !!! (i,j) == 0
268 -- | Returns True if the given matrix is lower-triangular, and False
273 -- >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
274 -- >>> is_lower_triangular m
277 -- >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
278 -- >>> is_lower_triangular m
281 is_lower_triangular :: (Eq a,
287 is_lower_triangular = is_upper_triangular . transpose
290 -- | Returns True if the given matrix is triangular, and False
295 -- >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
296 -- >>> is_triangular m
299 -- >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
300 -- >>> is_triangular m
303 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
304 -- >>> is_triangular m
307 is_triangular :: (Eq a,
313 is_triangular m = is_upper_triangular m || is_lower_triangular m
316 -- | Return the (i,j)th minor of m.
320 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
321 -- >>> minor m 0 0 :: Mat2 Int
323 -- >>> minor m 1 1 :: Mat2 Int
334 minor (Mat rows) i j = m
336 rows' = delete rows i
337 m = Mat $ V.map ((flip delete) j) rows'
340 class (Eq a, Ring.C a) => Determined p a where
341 determinant :: (p a) -> a
343 instance (Eq a, Ring.C a) => Determined (Mat (S Z) (S Z)) a where
344 determinant (Mat rows) = (V.head . V.head) rows
349 Determined (Mat (S n) (S n)) a)
350 => Determined (Mat (S (S n)) (S (S n))) a where
351 -- | The recursive definition with a special-case for triangular matrices.
355 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
360 | is_triangular m = product [ m !!! (i,i) | i <- [0..(nrows m)-1] ]
361 | otherwise = determinant_recursive
365 det_minor i j = determinant (minor m i j)
367 determinant_recursive =
368 sum [ (-1)^(toInteger j) NP.* (m' 0 j) NP.* (det_minor 0 j)
369 | j <- [0..(ncols m)-1] ]
373 -- | Matrix multiplication.
377 -- >>> let m1 = fromList [[1,2,3], [4,5,6]] :: Mat N2 N3 Int
378 -- >>> let m2 = fromList [[1,2],[3,4],[5,6]] :: Mat N3 N2 Int
383 (*) :: (Ring.C a, Arity m, Arity n, Arity p)
387 (*) m1 m2 = construct lambda
390 sum [(m1 !!! (i,k)) NP.* (m2 !!! (k,j)) | k <- [0..(ncols m1)-1] ]
394 instance (Ring.C a, Arity m, Arity n) => Additive.C (Mat m n a) where
396 (Mat rows1) + (Mat rows2) =
397 Mat $ V.zipWith (V.zipWith (+)) rows1 rows2
399 (Mat rows1) - (Mat rows2) =
400 Mat $ V.zipWith (V.zipWith (-)) rows1 rows2
402 zero = Mat (V.replicate $ V.replicate (fromInteger 0))
405 instance (Ring.C a, Arity m, Arity n, m ~ n) => Ring.C (Mat m n a) where
406 -- The first * is ring multiplication, the second is matrix
411 instance (Ring.C a, Arity m, Arity n) => Module.C a (Mat m n a) where
412 -- We can multiply a matrix by a scalar of the same type as its
414 x *> (Mat rows) = Mat $ V.map (V.map (NP.* x)) rows
417 instance (Algebraic.C a,
421 => Normed (Mat (S m) (S n) a) where
422 -- | Generic p-norms. The usual norm in R^n is (norm_p 2). We treat
423 -- all matrices as big vectors.
427 -- >>> let v1 = vec2d (3,4)
433 norm_p p (Mat rows) =
434 (root p') $ sum [(fromRational' $ toRational x)^p' | x <- xs]
437 xs = concat $ V.toList $ V.map V.toList rows
439 -- | The infinity norm.
443 -- >>> let v1 = vec3d (1,5,2)
447 norm_infty (Mat rows) =
448 fromRational' $ toRational $ V.maximum $ V.map V.maximum rows
454 -- Vector helpers. We want it to be easy to create low-dimension
455 -- column vectors, which are nx1 matrices.
457 -- | Convenient constructor for 2D vectors.
461 -- >>> import Roots.Simple
462 -- >>> let fst m = m !!! (0,0)
463 -- >>> let snd m = m !!! (1,0)
464 -- >>> let h = 0.5 :: Double
465 -- >>> let g1 m = 1.0 + h NP.* exp(-((fst m)^2))/(1.0 + (snd m)^2)
466 -- >>> let g2 m = 0.5 + h NP.* atan((fst m)^2 + (snd m)^2)
467 -- >>> let g u = vec2d ((g1 u), (g2 u))
468 -- >>> let u0 = vec2d (1.0, 1.0)
469 -- >>> let eps = 1/(10^9)
470 -- >>> fixed_point g eps u0
471 -- ((1.0728549599342185),(1.0820591495686167))
473 vec1d :: (a) -> Mat N1 N1 a
474 vec1d (x) = Mat (mk1 (mk1 x))
476 vec2d :: (a,a) -> Mat N2 N1 a
477 vec2d (x,y) = Mat (mk2 (mk1 x) (mk1 y))
479 vec3d :: (a,a,a) -> Mat N3 N1 a
480 vec3d (x,y,z) = Mat (mk3 (mk1 x) (mk1 y) (mk1 z))
482 vec4d :: (a,a,a,a) -> Mat N4 N1 a
483 vec4d (w,x,y,z) = Mat (mk4 (mk1 w) (mk1 x) (mk1 y) (mk1 z))
485 vec5d :: (a,a,a,a,a) -> Mat N5 N1 a
486 vec5d (v,w,x,y,z) = Mat (mk5 (mk1 v) (mk1 w) (mk1 x) (mk1 y) (mk1 z))
488 -- Since we commandeered multiplication, we need to create 1x1
489 -- matrices in order to multiply things.
490 scalar :: a -> Mat N1 N1 a
491 scalar x = Mat (mk1 (mk1 x))
493 dot :: (RealRing.C a, n ~ N1, m ~ S t, Arity t)
497 v1 `dot` v2 = ((transpose v1) * v2) !!! (0, 0)
500 -- | The angle between @v1@ and @v2@ in Euclidean space.
504 -- >>> let v1 = vec2d (1.0, 0.0)
505 -- >>> let v2 = vec2d (0.0, 1.0)
506 -- >>> angle v1 v2 == pi/2.0
509 angle :: (Transcendental.C a,
521 theta = (recip norms) NP.* (v1 `dot` v2)
522 norms = (norm v1) NP.* (norm v2)