]> gitweb.michael.orlitzky.com - numerical-analysis.git/blob - Linear/Matrix.hs
numerical-analysis.cabal,doc: relicense to AGPL-3.0+
[numerical-analysis.git] / Linear / Matrix.hs
1 {-# LANGUAGE ExistentialQuantification #-}
2 {-# LANGUAGE FlexibleContexts #-}
3 {-# LANGUAGE FlexibleInstances #-}
4 {-# LANGUAGE MultiParamTypeClasses #-}
5 {-# LANGUAGE NoMonomorphismRestriction #-}
6 {-# LANGUAGE ScopedTypeVariables #-}
7 {-# LANGUAGE TypeFamilies #-}
8 {-# LANGUAGE RebindableSyntax #-}
9
10 -- | Boxed matrices; that is, boxed m-vectors of boxed n-vectors. We
11 -- assume that the underlying representation is
12 -- Data.Vector.Fixed.Boxed.Vec for simplicity. It was tried in
13 -- generality and failed.
14 --
15 module Linear.Matrix
16 where
17
18 import Data.List (intercalate)
19
20 import Data.Vector.Fixed (
21 (!),
22 generate,
23 mk1,
24 mk2,
25 mk3,
26 mk4,
27 mk5 )
28 import qualified Data.Vector.Fixed as V (
29 and,
30 fromList,
31 head,
32 ifoldl,
33 ifoldr,
34 imap,
35 map,
36 maximum,
37 replicate,
38 reverse,
39 toList,
40 zipWith )
41 import Data.Vector.Fixed.Cont ( Arity, arity )
42 import Linear.Vector ( Vec, delete, element_sum )
43 import Naturals ( N1, N2, N3, N4, N5, N6, N7, N8, N9, N10, S, Z )
44 import Normed ( Normed(..) )
45
46 import NumericPrelude hiding ( (*), abs )
47 import qualified NumericPrelude as NP ( (*) )
48 import qualified Algebra.Absolute as Absolute ( C )
49 import Algebra.Absolute ( abs )
50 import qualified Algebra.Additive as Additive ( C )
51 import qualified Algebra.Algebraic as Algebraic ( C )
52 import Algebra.Algebraic ( root )
53 import qualified Algebra.Field as Field ( C )
54 import qualified Algebra.Ring as Ring ( C )
55 import qualified Algebra.Module as Module ( C )
56 import qualified Algebra.RealRing as RealRing ( C )
57 import qualified Algebra.ToRational as ToRational ( C )
58 import qualified Algebra.Transcendental as Transcendental ( C )
59 import qualified Prelude as P ( map )
60
61 -- | Our main matrix type.
62 data Mat m n a = (Arity m, Arity n) => Mat (Vec m (Vec n a))
63
64 -- Type synonyms for n-by-n matrices.
65 type Mat1 a = Mat N1 N1 a
66 type Mat2 a = Mat N2 N2 a
67 type Mat3 a = Mat N3 N3 a
68 type Mat4 a = Mat N4 N4 a
69 type Mat5 a = Mat N5 N5 a
70
71 -- * Type synonyms for 1-by-n row "vectors".
72
73 -- | Type synonym for row vectors expressed as 1-by-n matrices.
74 type Row n a = Mat N1 n a
75
76 type Row1 a = Row N1 a
77 type Row2 a = Row N2 a
78 type Row3 a = Row N3 a
79 type Row4 a = Row N4 a
80 type Row5 a = Row N5 a
81
82 -- * Type synonyms for n-by-1 column "vectors".
83
84 -- | Type synonym for column vectors expressed as n-by-1 matrices.
85 type Col n a = Mat n N1 a
86
87 type Col1 a = Col N1 a
88 type Col2 a = Col N2 a
89 type Col3 a = Col N3 a
90 type Col4 a = Col N4 a
91 type Col5 a = Col N5 a
92 type Col6 a = Col N6 a
93 type Col7 a = Col N7 a
94 type Col8 a = Col N8 a
95 type Col9 a = Col N9 a
96 type Col10 a = Col N10 a -- We need a big column for Gaussian quadrature.
97
98
99 instance (Eq a) => Eq (Mat m n a) where
100 -- | Compare a row at a time.
101 --
102 -- Examples:
103 --
104 -- >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int
105 -- >>> let m2 = fromList [[1,2],[3,4]] :: Mat2 Int
106 -- >>> let m3 = fromList [[5,6],[7,8]] :: Mat2 Int
107 -- >>> m1 == m2
108 -- True
109 -- >>> m1 == m3
110 -- False
111 --
112 (Mat rows1) == (Mat rows2) =
113 V.and $ V.zipWith comp rows1 rows2
114 where
115 -- Compare a row, one column at a time.
116 comp row1 row2 = V.and (V.zipWith (==) row1 row2)
117
118
119 instance (Show a) => Show (Mat m n a) where
120 -- | Display matrices and vectors as ordinary tuples. This is poor
121 -- practice, but these results are primarily displayed
122 -- interactively and convenience trumps correctness (said the guy
123 -- who insists his vector lengths be statically checked at
124 -- compile-time).
125 --
126 -- Examples:
127 --
128 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
129 -- >>> show m
130 -- ((1,2),(3,4))
131 --
132 show (Mat rows) =
133 "(" ++ (intercalate "," (V.toList row_strings)) ++ ")"
134 where
135 row_strings = V.map show_vector rows
136 show_vector v1 =
137 "(" ++ (intercalate "," element_strings) ++ ")"
138 where
139 v1l = V.toList v1
140 element_strings = P.map show v1l
141
142
143 -- | Convert a matrix to a nested list.
144 toList :: Mat m n a -> [[a]]
145 toList (Mat rows) = map V.toList (V.toList rows)
146
147
148 -- | Create a matrix from a nested list.
149 fromList :: (Arity m, Arity n) => [[a]] -> Mat m n a
150 fromList vs = Mat (V.fromList $ map V.fromList vs)
151
152
153 -- | Unsafe indexing. Much faster than the safe indexing.
154 (!!!) :: (Arity m, Arity n) => Mat m n a -> (Int, Int) -> a
155 (!!!) (Mat rows) (i, j) = (rows ! i) ! j
156
157
158 -- | Safe indexing.
159 --
160 -- Examples:
161 --
162 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
163 -- >>> m !!? (-1,-1)
164 -- Nothing
165 -- >>> m !!? (-1,0)
166 -- Nothing
167 -- >>> m !!? (-1,1)
168 -- Nothing
169 -- >>> m !!? (0,-1)
170 -- Nothing
171 -- >>> m !!? (0,0)
172 -- Just 1
173 -- >>> m !!? (0,1)
174 -- Just 2
175 -- >>> m !!? (1,-1)
176 -- Nothing
177 -- >>> m !!? (1,0)
178 -- Just 3
179 -- >>> m !!? (1,1)
180 -- Just 4
181 -- >>> m !!? (2,-1)
182 -- Nothing
183 -- >>> m !!? (2,0)
184 -- Nothing
185 -- >>> m !!? (2,1)
186 -- Nothing
187 -- >>> m !!? (2,2)
188 -- Nothing
189 --
190 (!!?) :: (Arity m, Arity n) => Mat m n a -> (Int, Int) -> Maybe a
191 (!!?) matrix idx =
192 ifoldl2 f Nothing matrix
193 where
194 f k l found cur = if (k,l) == idx then (Just cur) else found
195
196
197 -- | The number of rows in the matrix.
198 nrows :: forall m n a. (Arity m) => Mat m n a -> Int
199 nrows _ = arity (undefined :: m)
200
201
202 -- | The number of columns in the first row of the
203 -- matrix. Implementation stolen from Data.Vector.Fixed.length.
204 ncols :: forall m n a. (Arity n) => Mat m n a -> Int
205 ncols _ = arity (undefined :: n)
206
207
208 -- | Return the @i@th row of @m@ as a matrix. Unsafe.
209 row :: (Arity m, Arity n) => Mat m n a -> Int -> Row n a
210 row m i =
211 construct lambda
212 where
213 lambda _ j = m !!! (i, j)
214
215
216 -- | Return the @j@th column of @m@ as a matrix. Unsafe.
217 column :: (Arity m, Arity n) => Mat m n a -> Int -> Col m a
218 column m j =
219 construct lambda
220 where
221 lambda i _ = m !!! (i, j)
222
223
224 -- | Transpose @m@; switch it's columns and its rows. This is a dirty
225 -- implementation, but I don't see a better way.
226 --
227 -- TODO: Don't cheat with fromList.
228 --
229 -- Examples:
230 --
231 -- >>> let m = fromList [[1,2], [3,4]] :: Mat2 Int
232 -- >>> transpose m
233 -- ((1,3),(2,4))
234 --
235 transpose :: (Arity m, Arity n) => Mat m n a -> Mat n m a
236 transpose matrix =
237 construct lambda
238 where
239 lambda i j = matrix !!! (j,i)
240
241
242 -- | Is @m@ symmetric?
243 --
244 -- Examples:
245 --
246 -- >>> let m1 = fromList [[1,2], [2,1]] :: Mat2 Int
247 -- >>> symmetric m1
248 -- True
249 --
250 -- >>> let m2 = fromList [[1,2], [3,1]] :: Mat2 Int
251 -- >>> symmetric m2
252 -- False
253 --
254 symmetric :: (Eq a, Arity m) => Mat m m a -> Bool
255 symmetric m =
256 m == (transpose m)
257
258
259 -- | Construct a new matrix from a function @lambda@. The function
260 -- @lambda@ should take two parameters i,j corresponding to the
261 -- entries in the matrix. The i,j entry of the resulting matrix will
262 -- have the value returned by lambda i j.
263 --
264 -- Examples:
265 --
266 -- >>> let lambda i j = i + j
267 -- >>> construct lambda :: Mat3 Int
268 -- ((0,1,2),(1,2,3),(2,3,4))
269 --
270 construct :: forall m n a. (Arity m, Arity n)
271 => (Int -> Int -> a) -> Mat m n a
272 construct lambda = Mat $ generate make_row
273 where
274 make_row :: Int -> Vec n a
275 make_row i = generate (lambda i)
276
277
278 -- | Create an identity matrix with the right dimensions.
279 --
280 -- Examples:
281 --
282 -- >>> identity_matrix :: Mat3 Int
283 -- ((1,0,0),(0,1,0),(0,0,1))
284 -- >>> identity_matrix :: Mat3 Double
285 -- ((1.0,0.0,0.0),(0.0,1.0,0.0),(0.0,0.0,1.0))
286 --
287 identity_matrix :: (Arity m, Ring.C a) => Mat m m a
288 identity_matrix =
289 construct (\i j -> if i == j then (fromInteger 1) else (fromInteger 0))
290
291
292 -- | Given a positive-definite matrix @m@, computes the
293 -- upper-triangular matrix @r@ with (transpose r)*r == m and all
294 -- values on the diagonal of @r@ positive.
295 --
296 -- Examples:
297 --
298 -- >>> let m1 = fromList [[20,-1], [-1,20]] :: Mat2 Double
299 -- >>> let r = cholesky m1
300 -- >>> frobenius_norm ((transpose r)*r - m1) < 1e-10
301 -- True
302 -- >>> is_upper_triangular r
303 -- True
304 --
305 -- >>> let k1 = [6, -3, 0, 0, 0, 0, 0] :: [Double]
306 -- >>> let k2 = [-3, 10.5, -7.5, 0, 0, 0, 0] :: [Double]
307 -- >>> let k3 = [0, -7.5, 12.5, 0, 0, 0, 0] :: [Double]
308 -- >>> let k4 = [0, 0, 0, 6, 0, 0, 0] :: [Double]
309 -- >>> let k5 = [0, 0, 0, 0, 6, 0, 0] :: [Double]
310 -- >>> let k6 = [0, 0, 0, 0, 0, 6, 0] :: [Double]
311 -- >>> let k7 = [0, 0, 0, 0, 0, 0, 15] :: [Double]
312 -- >>> let big_K = fromList [k1,k2,k3,k4,k5,k6,k7] :: Mat N7 N7 Double
313 --
314 -- >>> let e1 = [2.449489742783178,0,0,0,0,0,0] :: [Double]
315 -- >>> let e2 = [-1.224744871391589,3,0,0,0,0,0] :: [Double]
316 -- >>> let e3 = [0,-5/2,5/2,0,0,0,0] :: [Double]
317 -- >>> let e4 = [0,0,0,2.449489742783178,0,0,0] :: [Double]
318 -- >>> let e5 = [0,0,0,0,2.449489742783178,0,0] :: [Double]
319 -- >>> let e6 = [0,0,0,0,0,2.449489742783178,0] :: [Double]
320 -- >>> let e7 = [0,0,0,0,0,0,3.872983346207417] :: [Double]
321 -- >>> let expected = fromList [e1,e2,e3,e4,e5,e6,e7] :: Mat N7 N7 Double
322 --
323 -- >>> let r = cholesky big_K
324 -- >>> frobenius_norm (r - (transpose expected)) < 1e-12
325 -- True
326 --
327 cholesky :: forall m n a. (Algebraic.C a, Arity m, Arity n)
328 => (Mat m n a) -> (Mat m n a)
329 cholesky m = construct r
330 where
331 r :: Int -> Int -> a
332 r i j | i == j = sqrt(m !!! (i,j) - sum [(r k i)^2 | k <- [0..i-1]])
333 | i < j =
334 (((m !!! (i,j)) - sum [(r k i) NP.* (r k j) | k <- [0..i-1]]))/(r i i)
335 | otherwise = 0
336
337
338 -- | Returns True if the given matrix is upper-triangular, and False
339 -- otherwise. The parameter @epsilon@ lets the caller choose a
340 -- tolerance.
341 --
342 -- Examples:
343 --
344 -- >>> let m = fromList [[1,1],[1e-12,1]] :: Mat2 Double
345 -- >>> is_upper_triangular m
346 -- False
347 -- >>> is_upper_triangular' 1e-10 m
348 -- True
349 --
350 is_upper_triangular' :: forall m n a.
351 (Ord a, Ring.C a, Absolute.C a, Arity m, Arity n)
352 => a -- ^ The tolerance @epsilon@.
353 -> Mat m n a
354 -> Bool
355 is_upper_triangular' epsilon matrix =
356 ifoldl2 f True matrix
357 where
358 f :: Int -> Int -> Bool -> a -> Bool
359 f _ _ False _ = False
360 f i j True x
361 | i <= j = True
362 -- use "less than or equal to" so zero is a valid epsilon
363 | otherwise = abs x <= epsilon
364
365
366 -- | Returns True if the given matrix is upper-triangular, and False
367 -- otherwise. We don't delegate to the general
368 -- 'is_upper_triangular'' here because it imposes additional
369 -- typeclass constraints throughout the library.
370 --
371 -- Examples:
372 --
373 -- >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
374 -- >>> is_upper_triangular m
375 -- False
376 --
377 -- >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
378 -- >>> is_upper_triangular m
379 -- True
380 --
381 is_upper_triangular :: forall m n a.
382 (Eq a, Ring.C a, Arity m, Arity n)
383 => Mat m n a -> Bool
384 is_upper_triangular matrix =
385 ifoldl2 f True matrix
386 where
387 f :: Int -> Int -> Bool -> a -> Bool
388 f _ _ False _ = False
389 f i j True x
390 | i <= j = True
391 | otherwise = x == 0
392
393
394
395 -- | Returns True if the given matrix is lower-triangular, and False
396 -- otherwise.
397 --
398 -- Examples:
399 --
400 -- >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
401 -- >>> is_lower_triangular m
402 -- True
403 --
404 -- >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
405 -- >>> is_lower_triangular m
406 -- False
407 --
408 is_lower_triangular :: (Eq a,
409 Ring.C a,
410 Arity m,
411 Arity n)
412 => Mat m n a
413 -> Bool
414 is_lower_triangular = is_upper_triangular . transpose
415
416
417 -- | Returns True if the given matrix is lower-triangular, and False
418 -- otherwise. The parameter @epsilon@ lets the caller choose a
419 -- tolerance.
420 --
421 -- Examples:
422 --
423 -- >>> let m = fromList [[1,1e-12],[1,1]] :: Mat2 Double
424 -- >>> is_lower_triangular m
425 -- False
426 -- >>> is_lower_triangular' 1e-12 m
427 -- True
428 --
429 is_lower_triangular' :: (Ord a,
430 Ring.C a,
431 Absolute.C a,
432 Arity m,
433 Arity n)
434 => a -- ^ The tolerance @epsilon@.
435 -> Mat m n a
436 -> Bool
437 is_lower_triangular' epsilon = (is_upper_triangular' epsilon) . transpose
438
439
440 -- | Returns True if the given matrix is triangular, and False
441 -- otherwise.
442 --
443 -- Examples:
444 --
445 -- >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
446 -- >>> is_triangular m
447 -- True
448 --
449 -- >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
450 -- >>> is_triangular m
451 -- True
452 --
453 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
454 -- >>> is_triangular m
455 -- False
456 --
457 is_triangular :: (Ord a,
458 Ring.C a,
459 Absolute.C a,
460 Arity m,
461 Arity n)
462 => Mat m n a
463 -> Bool
464 is_triangular m = is_upper_triangular m || is_lower_triangular m
465
466
467 -- | Delete the @i@th row and @j@th column from the matrix. The name
468 -- \"preminor\" is made up, but is meant to signify that this is
469 -- usually used in the computationof a minor. A minor is simply the
470 -- determinant of a preminor in that case.
471 --
472 -- Examples:
473 --
474 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
475 -- >>> preminor m 0 0 :: Mat2 Int
476 -- ((5,6),(8,9))
477 -- >>> preminor m 1 1 :: Mat2 Int
478 -- ((1,3),(7,9))
479 --
480 preminor :: (Arity m, Arity n)
481 => Mat (S m) (S n) a
482 -> Int
483 -> Int
484 -> Mat m n a
485 preminor (Mat rows) i j = m
486 where
487 rows' = delete rows i
488 m = Mat $ V.map ((flip delete) j) rows'
489
490
491 -- | Compute the i,jth minor of a @matrix@.
492 --
493 -- Examples:
494 --
495 -- >>> let m1 = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Double
496 -- >>> minor m1 1 1
497 -- -12.0
498 --
499 minor :: (Arity m, Determined (Mat m m) a)
500 => Mat (S m) (S m) a
501 -> Int
502 -> Int
503 -> a
504 minor matrix i j = determinant (preminor matrix i j)
505
506 class (Eq a, Ring.C a) => Determined p a where
507 determinant :: (p a) -> a
508
509 instance (Eq a, Ring.C a) => Determined (Mat (S Z) (S Z)) a where
510 determinant = unscalar
511
512 instance (Ord a,
513 Ring.C a,
514 Absolute.C a,
515 Arity n,
516 Determined (Mat (S n) (S n)) a)
517 => Determined (Mat (S (S n)) (S (S n))) a where
518 -- | The recursive definition with a special-case for triangular matrices.
519 --
520 -- Examples:
521 --
522 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
523 -- >>> determinant m
524 -- -1
525 --
526 determinant m
527 | is_triangular m = product [ m !!! (i,i) | i <- [0..(nrows m)-1] ]
528 | otherwise = determinant_recursive
529 where
530 m' i j = m !!! (i,j)
531
532 determinant_recursive =
533 sum [ (-1)^(toInteger j) NP.* (m' 0 j) NP.* (minor m 0 j)
534 | j <- [0..(ncols m)-1] ]
535
536
537
538 -- | Matrix multiplication.
539 --
540 -- Examples:
541 --
542 -- >>> let m1 = fromList [[1,2,3], [4,5,6]] :: Mat N2 N3 Int
543 -- >>> let m2 = fromList [[1,2],[3,4],[5,6]] :: Mat N3 N2 Int
544 -- >>> m1 * m2
545 -- ((22,28),(49,64))
546 --
547 infixl 7 *
548 (*) :: (Ring.C a, Arity m, Arity n, Arity p)
549 => Mat m n a
550 -> Mat n p a
551 -> Mat m p a
552 (*) m1 m2 = construct lambda
553 where
554 lambda i j =
555 sum [(m1 !!! (i,k)) NP.* (m2 !!! (k,j)) | k <- [0..(ncols m1)-1] ]
556
557
558
559 instance (Ring.C a, Arity m, Arity n) => Additive.C (Mat m n a) where
560
561 (Mat rows1) + (Mat rows2) =
562 Mat $ V.zipWith (V.zipWith (+)) rows1 rows2
563
564 (Mat rows1) - (Mat rows2) =
565 Mat $ V.zipWith (V.zipWith (-)) rows1 rows2
566
567 zero = Mat (V.replicate $ V.replicate (fromInteger 0))
568
569
570 instance (Ring.C a, Arity m, Arity n, m ~ n) => Ring.C (Mat m n a) where
571 -- The first * is ring multiplication, the second is matrix
572 -- multiplication.
573 m1 * m2 = m1 * m2
574
575
576 instance (Ring.C a, Arity m, Arity n) => Module.C a (Mat m n a) where
577 -- We can multiply a matrix by a scalar of the same type as its
578 -- elements.
579 x *> (Mat rows) = Mat $ V.map (V.map (NP.* x)) rows
580
581
582 instance (Absolute.C a,
583 Algebraic.C a,
584 ToRational.C a,
585 Arity m)
586 => Normed (Col (S m) a) where
587 -- | Generic p-norms for vectors in R^n that are represented as n-by-1
588 -- matrices.
589 --
590 -- Examples:
591 --
592 -- >>> let v1 = vec2d (3,4)
593 -- >>> norm_p 1 v1
594 -- 7.0
595 -- >>> norm_p 2 v1
596 -- 5.0
597 --
598 -- >>> let v1 = vec2d (-1,1) :: Col2 Double
599 -- >>> norm_p 1 v1 :: Double
600 -- 2.0
601 --
602 norm_p p (Mat rows) =
603 (root p') $ sum [fromRational' (toRational $ abs x)^p' | x <- xs]
604 where
605 p' = toInteger p
606 xs = concat $ V.toList $ V.map V.toList rows
607
608 -- | The infinity norm.
609 --
610 -- Examples:
611 --
612 -- >>> let v1 = vec3d (1,5,2)
613 -- >>> norm_infty v1
614 -- 5
615 --
616 norm_infty (Mat rows) =
617 fromRational' $ toRational $ V.maximum $ V.map V.maximum rows
618
619
620 -- | Compute the Frobenius norm of a matrix. This essentially treats
621 -- the matrix as one long vector containing all of its entries (in
622 -- any order, it doesn't matter).
623 --
624 -- Examples:
625 --
626 -- >>> let m = fromList [[1, 2, 3],[4,5,6],[7,8,9]] :: Mat3 Double
627 -- >>> frobenius_norm m == sqrt 285
628 -- True
629 --
630 -- >>> let m = fromList [[1, -1, 1],[-1,1,-1],[1,-1,1]] :: Mat3 Double
631 -- >>> frobenius_norm m == 3
632 -- True
633 --
634 frobenius_norm :: (Algebraic.C a, Ring.C a) => Mat m n a -> a
635 frobenius_norm (Mat rows) =
636 sqrt $ element_sum $ V.map row_sum rows
637 where
638 -- | Square and add up the entries of a row.
639 row_sum = element_sum . V.map (^2)
640
641
642 -- Vector helpers. We want it to be easy to create low-dimension
643 -- column vectors, which are nx1 matrices.
644
645 -- | Convenient constructor for 2D vectors.
646 --
647 -- Examples:
648 --
649 -- >>> import Roots.Simple
650 -- >>> let fst m = m !!! (0,0)
651 -- >>> let snd m = m !!! (1,0)
652 -- >>> let h = 0.5 :: Double
653 -- >>> let g1 m = 1.0 + h NP.* exp(-((fst m)^2))/(1.0 + (snd m)^2)
654 -- >>> let g2 m = 0.5 + h NP.* atan((fst m)^2 + (snd m)^2)
655 -- >>> let g u = vec2d ((g1 u), (g2 u))
656 -- >>> let u0 = vec2d (1.0, 1.0)
657 -- >>> let eps = 1/(10^9)
658 -- >>> fixed_point g eps u0
659 -- ((1.0728549599342185),(1.0820591495686167))
660 --
661 vec1d :: (a) -> Col1 a
662 vec1d (x) = Mat (mk1 (mk1 x))
663
664 vec2d :: (a,a) -> Col2 a
665 vec2d (x,y) = Mat (mk2 (mk1 x) (mk1 y))
666
667 vec3d :: (a,a,a) -> Col3 a
668 vec3d (x,y,z) = Mat (mk3 (mk1 x) (mk1 y) (mk1 z))
669
670 vec4d :: (a,a,a,a) -> Col4 a
671 vec4d (w,x,y,z) = Mat (mk4 (mk1 w) (mk1 x) (mk1 y) (mk1 z))
672
673 vec5d :: (a,a,a,a,a) -> Col5 a
674 vec5d (v,w,x,y,z) = Mat (mk5 (mk1 v) (mk1 w) (mk1 x) (mk1 y) (mk1 z))
675
676
677 -- Since we commandeered multiplication, we need to create 1x1
678 -- matrices in order to multiply things.
679 scalar :: a -> Mat1 a
680 scalar x = Mat (mk1 (mk1 x))
681
682 -- Get the scalar value out of a 1x1 matrix.
683 unscalar :: Mat1 a -> a
684 unscalar (Mat rows) = V.head $ V.head rows
685
686
687 dot :: (Ring.C a, Arity m)
688 => Col (S m) a
689 -> Col (S m) a
690 -> a
691 v1 `dot` v2 = unscalar $ ((transpose v1) * v2)
692
693
694 -- | The angle between @v1@ and @v2@ in Euclidean space.
695 --
696 -- Examples:
697 --
698 -- >>> let v1 = vec2d (1.0, 0.0)
699 -- >>> let v2 = vec2d (0.0, 1.0)
700 -- >>> angle v1 v2 == pi/2.0
701 -- True
702 --
703 angle :: (Transcendental.C a,
704 RealRing.C a,
705 m ~ S t,
706 Arity t,
707 ToRational.C a)
708 => Col m a
709 -> Col m a
710 -> a
711 angle v1 v2 =
712 acos theta
713 where
714 theta = (recip norms) NP.* (v1 `dot` v2)
715 norms = (norm v1) NP.* (norm v2)
716
717
718 -- | Retrieve the diagonal elements of the given matrix as a \"column
719 -- vector,\" i.e. a m-by-1 matrix. We require the matrix to be
720 -- square to avoid ambiguity in the return type which would ideally
721 -- have dimension min(m,n) supposing an m-by-n matrix.
722 --
723 -- Examples:
724 --
725 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
726 -- >>> diagonal m
727 -- ((1),(5),(9))
728 --
729 diagonal :: (Arity m) => Mat m m a -> Col m a
730 diagonal matrix =
731 construct lambda
732 where
733 lambda i _ = matrix !!! (i,i)
734
735
736 -- | Given a square @matrix@, return a new matrix of the same size
737 -- containing only the on-diagonal entries of @matrix@. The
738 -- off-diagonal entries are set to zero.
739 --
740 -- Examples:
741 --
742 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
743 -- >>> diagonal_part m
744 -- ((1,0,0),(0,5,0),(0,0,9))
745 --
746 diagonal_part :: (Arity m, Ring.C a)
747 => Mat m m a
748 -> Mat m m a
749 diagonal_part matrix =
750 construct lambda
751 where
752 lambda i j = if i == j then matrix !!! (i,j) else 0
753
754
755 -- | Given a square @matrix@, return a new matrix of the same size
756 -- containing only the on-diagonal and below-diagonal entries of
757 -- @matrix@. The above-diagonal entries are set to zero.
758 --
759 -- Examples:
760 --
761 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
762 -- >>> lt_part m
763 -- ((1,0,0),(4,5,0),(7,8,9))
764 --
765 lt_part :: (Arity m, Ring.C a)
766 => Mat m m a
767 -> Mat m m a
768 lt_part matrix =
769 construct lambda
770 where
771 lambda i j = if i >= j then matrix !!! (i,j) else 0
772
773
774 -- | Given a square @matrix@, return a new matrix of the same size
775 -- containing only the below-diagonal entries of @matrix@. The on-
776 -- and above-diagonal entries are set to zero.
777 --
778 -- Examples:
779 --
780 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
781 -- >>> lt_part_strict m
782 -- ((0,0,0),(4,0,0),(7,8,0))
783 --
784 lt_part_strict :: (Arity m, Ring.C a)
785 => Mat m m a
786 -> Mat m m a
787 lt_part_strict matrix =
788 construct lambda
789 where
790 lambda i j = if i > j then matrix !!! (i,j) else 0
791
792
793 -- | Given a square @matrix@, return a new matrix of the same size
794 -- containing only the on-diagonal and above-diagonal entries of
795 -- @matrix@. The below-diagonal entries are set to zero.
796 --
797 -- Examples:
798 --
799 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
800 -- >>> ut_part m
801 -- ((1,2,3),(0,5,6),(0,0,9))
802 --
803 ut_part :: (Arity m, Ring.C a)
804 => Mat m m a
805 -> Mat m m a
806 ut_part = transpose . lt_part . transpose
807
808
809 -- | Given a square @matrix@, return a new matrix of the same size
810 -- containing only the above-diagonal entries of @matrix@. The on-
811 -- and below-diagonal entries are set to zero.
812 --
813 -- Examples:
814 --
815 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
816 -- >>> ut_part_strict m
817 -- ((0,2,3),(0,0,6),(0,0,0))
818 --
819 ut_part_strict :: (Arity m, Ring.C a)
820 => Mat m m a
821 -> Mat m m a
822 ut_part_strict = transpose . lt_part_strict . transpose
823
824
825 -- | Compute the trace of a square matrix, the sum of the elements
826 -- which lie on its diagonal. We require the matrix to be
827 -- square to avoid ambiguity in the return type which would ideally
828 -- have dimension min(m,n) supposing an m-by-n matrix.
829 --
830 -- Examples:
831 --
832 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
833 -- >>> trace m
834 -- 15
835 --
836 trace :: (Arity m, Ring.C a) => Mat m m a -> a
837 trace matrix =
838 let (Mat rows) = diagonal matrix
839 in
840 element_sum $ V.map V.head rows
841
842
843 -- | Zip together two matrices.
844 --
845 -- TODO: don't cheat with construct (map V.zips instead).
846 --
847 -- Examples:
848 --
849 -- >>> let m1 = fromList [[1],[1],[1]] :: Col3 Int
850 -- >>> let m2 = fromList [[1],[2],[3]] :: Col3 Int
851 -- >>> zip2 m1 m2
852 -- (((1,1)),((1,2)),((1,3)))
853 --
854 -- >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int
855 -- >>> let m2 = fromList [[1,1],[1,1]] :: Mat2 Int
856 -- >>> zip2 m1 m2
857 -- (((1,1),(2,1)),((3,1),(4,1)))
858 --
859 zip2 :: (Arity m, Arity n) => Mat m n a -> Mat m n a -> Mat m n (a,a)
860 zip2 m1 m2 =
861 construct lambda
862 where
863 lambda i j = (m1 !!! (i,j), m2 !!! (i,j))
864
865
866 -- | Zip together three matrices.
867 --
868 -- TODO: don't cheat with construct (map V.zips instead).
869 --
870 -- Examples:
871 --
872 -- >>> let m1 = fromList [[1],[1],[1]] :: Col3 Int
873 -- >>> let m2 = fromList [[1],[2],[3]] :: Col3 Int
874 -- >>> let m3 = fromList [[4],[5],[6]] :: Col3 Int
875 -- >>> zip2three m1 m2 m3
876 -- (((1,1,4)),((1,2,5)),((1,3,6)))
877 --
878 -- >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int
879 -- >>> let m2 = fromList [[1,1],[1,1]] :: Mat2 Int
880 -- >>> let m3 = fromList [[8,2],[6,3]] :: Mat2 Int
881 -- >>> zip2three m1 m2 m3
882 -- (((1,1,8),(2,1,2)),((3,1,6),(4,1,3)))
883 --
884 zip2three :: (Arity m, Arity n)
885 => Mat m n a
886 -> Mat m n a
887 -> Mat m n a
888 -> Mat m n (a,a,a)
889 zip2three m1 m2 m3 =
890 construct lambda
891 where
892 lambda i j = (m1 !!! (i,j), m2 !!! (i,j), m3 !!! (i,j))
893
894
895 -- | Zip together two matrices using the supplied function.
896 --
897 -- Examples:
898 --
899 -- >>> let c1 = fromList [[1],[2],[3]] :: Col3 Integer
900 -- >>> let c2 = fromList [[4],[5],[6]] :: Col3 Integer
901 -- >>> zipwith2 (^) c1 c2
902 -- ((1),(32),(729))
903 --
904 zipwith2 :: Arity m
905 => (a -> a -> b)
906 -> Col m a
907 -> Col m a
908 -> Col m b
909 zipwith2 f c1 c2 =
910 construct lambda
911 where
912 lambda i j = f (c1 !!! (i,j)) (c2 !!! (i,j))
913
914
915 -- | Map a function over a matrix of any dimensions.
916 --
917 -- Examples:
918 --
919 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
920 -- >>> map2 (^2) m
921 -- ((1,4),(9,16))
922 --
923 map2 :: (a -> b) -> Mat m n a -> Mat m n b
924 map2 f (Mat rows) =
925 Mat $ V.map g rows
926 where
927 g = V.map f
928
929
930 -- | Fold over the entire matrix passing the coordinates @i@ and @j@
931 -- (of the row/column) to the accumulation function. The fold occurs
932 -- from top-left to bottom-right.
933 --
934 -- Examples:
935 --
936 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
937 -- >>> ifoldl2 (\i j cur _ -> cur + i + j) 0 m
938 -- 18
939 --
940 ifoldl2 :: forall a b m n.
941 (Int -> Int -> b -> a -> b)
942 -> b
943 -> Mat m n a
944 -> b
945 ifoldl2 f initial (Mat rows) =
946 V.ifoldl row_function initial rows
947 where
948 -- | The order that we need this in (so that @g idx@ makes sense)
949 -- is a little funny. So that we don't need to pass weird
950 -- functions into ifoldl2, we swap the second and third
951 -- arguments of @f@ calling the result @g@.
952 g :: Int -> b -> Int -> a -> b
953 g w x y = f w y x
954
955 row_function :: b -> Int -> Vec n a -> b
956 row_function rowinit idx r = V.ifoldl (g idx) rowinit r
957
958
959 -- | Fold over the entire matrix passing the coordinates @i@ and @j@
960 -- (of the row/column) to the accumulation function. The fold occurs
961 -- from bottom-right to top-left.
962 --
963 -- The order of the arguments in the supplied function are different
964 -- from those in V.ifoldr; we keep them similar to ifoldl2.
965 --
966 -- Examples:
967 --
968 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
969 -- >>> ifoldr2 (\i j cur _ -> cur + i + j) 0 m
970 -- 18
971 --
972 ifoldr2 :: forall a b m n.
973 (Int -> Int -> b -> a -> b)
974 -> b
975 -> Mat m n a
976 -> b
977 ifoldr2 f initial (Mat rows) =
978 V.ifoldr row_function initial rows
979 where
980 -- | Swap the order of arguments in @f@ so that it agrees with the
981 -- @f@ passed to ifoldl2.
982 g :: Int -> Int -> a -> b -> b
983 g w x y z = f w x z y
984
985 row_function :: Int -> Vec n a -> b -> b
986 row_function idx r rowinit = V.ifoldr (g idx) rowinit r
987
988
989 -- | Map a function over a matrix of any dimensions, passing the
990 -- coordinates @i@ and @j@ to the function @f@.
991 --
992 -- Examples:
993 --
994 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
995 -- >>> imap2 (\i j _ -> i+j) m
996 -- ((0,1),(1,2))
997 --
998 imap2 :: (Int -> Int -> a -> b) -> Mat m n a -> Mat m n b
999 imap2 f (Mat rows) =
1000 Mat $ V.imap g rows
1001 where
1002 g i = V.imap (f i)
1003
1004
1005 -- | Reverse the order of elements in a matrix.
1006 --
1007 -- Examples:
1008 --
1009 -- >>> let m1 = fromList [[1,2,3]] :: Row3 Int
1010 -- >>> reverse2 m1
1011 -- ((3,2,1))
1012 --
1013 -- >>> let m1 = vec3d (1,2,3 :: Int)
1014 -- >>> reverse2 m1
1015 -- ((3),(2),(1))
1016 --
1017 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
1018 -- >>> reverse2 m
1019 -- ((9,8,7),(6,5,4),(3,2,1))
1020 --
1021 reverse2 :: (Arity m, Arity n) => Mat m n a -> Mat m n a
1022 reverse2 (Mat rows) = Mat $ V.reverse $ V.map V.reverse rows
1023
1024
1025 -- | Unsafely set the (i,j) element of the given matrix.
1026 --
1027 -- Examples:
1028 --
1029 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
1030 -- >>> set_idx m (1,1) 17
1031 -- ((1,2,3),(4,17,6),(7,8,9))
1032 --
1033 set_idx :: forall m n a.
1034 (Arity m, Arity n)
1035 => Mat m n a
1036 -> (Int, Int)
1037 -> a
1038 -> Mat m n a
1039 set_idx matrix (i,j) newval =
1040 imap2 updater matrix
1041 where
1042 updater :: Int -> Int -> a -> a
1043 updater k l existing =
1044 if k == i && l == j
1045 then newval
1046 else existing
1047
1048
1049 -- | Compute the i,jth cofactor of the given @matrix@. This simply
1050 -- premultiplues the i,jth minor by (-1)^(i+j).
1051 cofactor :: (Arity m, Determined (Mat m m) a)
1052 => Mat (S m) (S m) a
1053 -> Int
1054 -> Int
1055 -> a
1056 cofactor matrix i j =
1057 (-1)^(toInteger i + toInteger j) NP.* (minor matrix i j)
1058
1059
1060 -- | Compute the inverse of a matrix using cofactor expansion
1061 -- (generalized Cramer's rule).
1062 --
1063 -- Examples:
1064 --
1065 -- >>> let m1 = fromList [[37,22],[17,54]] :: Mat2 Double
1066 -- >>> let e1 = [54/1624, -22/1624] :: [Double]
1067 -- >>> let e2 = [-17/1624, 37/1624] :: [Double]
1068 -- >>> let expected = fromList [e1, e2] :: Mat2 Double
1069 -- >>> let actual = inverse m1
1070 -- >>> frobenius_norm (actual - expected) < 1e-12
1071 -- True
1072 --
1073 inverse :: (Arity m,
1074 Determined (Mat (S m) (S m)) a,
1075 Determined (Mat m m) a,
1076 Field.C a)
1077 => Mat (S m) (S m) a
1078 -> Mat (S m) (S m) a
1079 inverse matrix =
1080 (1 / (determinant matrix)) *> (transpose $ construct lambda)
1081 where
1082 lambda i j = cofactor matrix i j
1083