1 {-# LANGUAGE ExistentialQuantification #-}
2 {-# LANGUAGE FlexibleContexts #-}
3 {-# LANGUAGE FlexibleInstances #-}
4 {-# LANGUAGE MultiParamTypeClasses #-}
5 {-# LANGUAGE NoMonomorphismRestriction #-}
6 {-# LANGUAGE ScopedTypeVariables #-}
7 {-# LANGUAGE TypeFamilies #-}
8 {-# LANGUAGE RebindableSyntax #-}
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.
18 import Data.List (intercalate)
20 import Data.Vector.Fixed (
29 import qualified Data.Vector.Fixed as V (
42 import Data.Vector.Fixed.Cont ( Arity, arity )
43 import Linear.Vector ( Vec, delete )
45 import Normed ( Normed(..) )
47 -- We want the "max" that works on Ord, not the one that only works on
48 -- Bool/Integer from the Lattice class!
49 import NumericPrelude hiding ( (*), abs, max)
50 import qualified NumericPrelude as NP ( (*) )
51 import qualified Algebra.Absolute as Absolute ( C )
52 import Algebra.Absolute ( abs )
53 import qualified Algebra.Additive as Additive ( C )
54 import qualified Algebra.Algebraic as Algebraic ( C )
55 import Algebra.Algebraic ( root )
56 import qualified Algebra.Field as Field ( C )
57 import qualified Algebra.Ring as Ring ( C )
58 import qualified Algebra.Module as Module ( C )
59 import qualified Algebra.RealRing as RealRing ( C )
60 import qualified Algebra.ToRational as ToRational ( C )
61 import qualified Algebra.Transcendental as Transcendental ( C )
62 import qualified Prelude as P ( map, max)
64 -- | Our main matrix type.
65 data Mat m n a = (Arity m, Arity n) => Mat (Vec m (Vec n a))
67 -- Type synonyms for n-by-n matrices.
68 type Mat0 a = Mat Z Z a
69 type Mat1 a = Mat N1 N1 a
70 type Mat2 a = Mat N2 N2 a
71 type Mat3 a = Mat N3 N3 a
72 type Mat4 a = Mat N4 N4 a
73 type Mat5 a = Mat N5 N5 a
74 type Mat6 a = Mat N6 N6 a
75 type Mat7 a = Mat N7 N7 a
77 -- * Type synonyms for 1-by-n row "vectors".
79 -- | Type synonym for row vectors expressed as 1-by-n matrices.
80 type Row n a = Mat N1 n a
82 type Row1 a = Row N1 a
83 type Row2 a = Row N2 a
84 type Row3 a = Row N3 a
85 type Row4 a = Row N4 a
86 type Row5 a = Row N5 a
88 -- * Type synonyms for n-by-1 column "vectors".
90 -- | Type synonym for column vectors expressed as n-by-1 matrices.
91 type Col n a = Mat n N1 a
94 type Col1 a = Col N1 a
95 type Col2 a = Col N2 a
96 type Col3 a = Col N3 a
97 type Col4 a = Col N4 a
98 type Col5 a = Col N5 a
99 type Col6 a = Col N6 a
100 type Col7 a = Col N7 a
101 type Col8 a = Col N8 a
102 type Col9 a = Col N9 a
103 type Col10 a = Col N10 a
104 type Col11 a = Col N11 a
105 type Col12 a = Col N12 a
106 type Col13 a = Col N13 a
107 type Col14 a = Col N14 a
108 type Col15 a = Col N15 a
109 type Col16 a = Col N16 a
110 type Col17 a = Col N17 a
111 type Col18 a = Col N18 a
112 type Col19 a = Col N19 a
113 type Col20 a = Col N20 a
114 type Col21 a = Col N21 a
115 type Col22 a = Col N22 a
116 type Col23 a = Col N23 a
117 type Col24 a = Col N24 a
118 type Col25 a = Col N25 a
119 type Col26 a = Col N26 a
120 type Col27 a = Col N27 a
121 type Col28 a = Col N28 a
122 type Col29 a = Col N29 a
123 type Col30 a = Col N30 a
124 type Col31 a = Col N31 a
125 type Col32 a = Col N32 a
128 instance (Eq a) => Eq (Mat m n a) where
129 -- | Compare a row at a time.
133 -- >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int
134 -- >>> let m2 = fromList [[1,2],[3,4]] :: Mat2 Int
135 -- >>> let m3 = fromList [[5,6],[7,8]] :: Mat2 Int
141 (Mat rows_one) == (Mat rows_two) =
142 V.and $ V.zipWith comp rows_one rows_two
144 -- Compare a row, one column at a time.
145 comp row1 row2 = V.and (V.zipWith (==) row1 row2)
148 instance (Show a) => Show (Mat m n a) where
149 -- | Display matrices and vectors as ordinary tuples. This is poor
150 -- practice, but these results are primarily displayed
151 -- interactively and convenience trumps correctness (said the guy
152 -- who insists his vector lengths be statically checked at
157 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
162 "(" ++ (intercalate "," (V.toList row_strings)) ++ ")"
164 row_strings = V.map show_vector rows
166 "(" ++ (intercalate "," element_strings) ++ ")"
169 element_strings = P.map show v1l
172 -- | Convert a matrix to a nested list.
173 toList :: Mat m n a -> [[a]]
174 toList (Mat rows) = map V.toList (V.toList rows)
177 -- | Create a matrix from a nested list.
178 fromList :: (Arity m, Arity n) => [[a]] -> Mat m n a
179 fromList vs = Mat (V.fromList $ map V.fromList vs)
182 -- | Unsafe indexing. Much faster than the safe indexing.
183 (!!!) :: (Arity m, Arity n) => Mat m n a -> (Int, Int) -> a
184 (!!!) (Mat rows) (i, j) = (rows ! i) ! j
191 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
219 (!!?) :: (Arity m, Arity n) => Mat m n a -> (Int, Int) -> Maybe a
221 ifoldl2 f Nothing matrix
223 f k l found cur = if (k,l) == idx then (Just cur) else found
226 -- | The number of rows in the matrix.
227 nrows :: forall m n a. (Arity m) => Mat m n a -> Int
228 nrows _ = arity (undefined :: m)
231 -- | The number of columns in the first row of the
232 -- matrix. Implementation stolen from Data.Vector.Fixed.length.
233 ncols :: forall m n a. (Arity n) => Mat m n a -> Int
234 ncols _ = arity (undefined :: n)
237 -- | Return the @i@th row of @m@ as a matrix. Unsafe.
238 row :: (Arity m, Arity n) => Mat m n a -> Int -> Row n a
242 lambda _ j = m !!! (i, j)
245 -- | Return the @j@th column of @m@ as a matrix. Unsafe.
246 column :: (Arity m, Arity n) => Mat m n a -> Int -> Col m a
250 lambda i _ = m !!! (i, j)
253 -- | Transpose @m@; switch it's columns and its rows. This is a dirty
254 -- implementation, but I don't see a better way.
256 -- TODO: Don't cheat with fromList.
260 -- >>> let m = fromList [[1,2], [3,4]] :: Mat2 Int
264 transpose :: (Arity m, Arity n) => Mat m n a -> Mat n m a
268 lambda i j = matrix !!! (j,i)
271 -- | Is @m@ symmetric?
275 -- >>> let m1 = fromList [[1,2], [2,1]] :: Mat2 Int
279 -- >>> let m2 = fromList [[1,2], [3,1]] :: Mat2 Int
283 symmetric :: (Eq a, Arity m) => Mat m m a -> Bool
288 -- | Construct a new matrix from a function @lambda@. The function
289 -- @lambda@ should take two parameters i,j corresponding to the
290 -- entries in the matrix. The i,j entry of the resulting matrix will
291 -- have the value returned by lambda i j.
295 -- >>> let lambda i j = i + j
296 -- >>> construct lambda :: Mat3 Int
297 -- ((0,1,2),(1,2,3),(2,3,4))
299 construct :: forall m n a. (Arity m, Arity n)
300 => (Int -> Int -> a) -> Mat m n a
301 construct lambda = Mat $ generate make_row
303 make_row :: Int -> Vec n a
304 make_row i = generate (lambda i)
307 -- | Create an identity matrix with the right dimensions.
311 -- >>> identity_matrix :: Mat3 Int
312 -- ((1,0,0),(0,1,0),(0,0,1))
313 -- >>> identity_matrix :: Mat3 Double
314 -- ((1.0,0.0,0.0),(0.0,1.0,0.0),(0.0,0.0,1.0))
316 identity_matrix :: (Arity m, Ring.C a) => Mat m m a
318 construct (\i j -> if i == j then (fromInteger 1) else (fromInteger 0))
321 -- | Given a positive-definite matrix @m@, computes the
322 -- upper-triangular matrix @r@ with (transpose r)*r == m and all
323 -- values on the diagonal of @r@ positive.
327 -- >>> let m1 = fromList [[20,-1], [-1,20]] :: Mat2 Double
328 -- >>> let r = cholesky m1
329 -- >>> frobenius_norm ((transpose r)*r - m1) < 1e-10
331 -- >>> is_upper_triangular r
334 -- >>> let k1 = [6, -3, 0, 0, 0, 0, 0] :: [Double]
335 -- >>> let k2 = [-3, 10.5, -7.5, 0, 0, 0, 0] :: [Double]
336 -- >>> let k3 = [0, -7.5, 12.5, 0, 0, 0, 0] :: [Double]
337 -- >>> let k4 = [0, 0, 0, 6, 0, 0, 0] :: [Double]
338 -- >>> let k5 = [0, 0, 0, 0, 6, 0, 0] :: [Double]
339 -- >>> let k6 = [0, 0, 0, 0, 0, 6, 0] :: [Double]
340 -- >>> let k7 = [0, 0, 0, 0, 0, 0, 15] :: [Double]
341 -- >>> let big_K = fromList [k1,k2,k3,k4,k5,k6,k7] :: Mat N7 N7 Double
343 -- >>> let e1 = [2.449489742783178,0,0,0,0,0,0] :: [Double]
344 -- >>> let e2 = [-1.224744871391589,3,0,0,0,0,0] :: [Double]
345 -- >>> let e3 = [0,-5/2,5/2,0,0,0,0] :: [Double]
346 -- >>> let e4 = [0,0,0,2.449489742783178,0,0,0] :: [Double]
347 -- >>> let e5 = [0,0,0,0,2.449489742783178,0,0] :: [Double]
348 -- >>> let e6 = [0,0,0,0,0,2.449489742783178,0] :: [Double]
349 -- >>> let e7 = [0,0,0,0,0,0,3.872983346207417] :: [Double]
350 -- >>> let expected = fromList [e1,e2,e3,e4,e5,e6,e7] :: Mat N7 N7 Double
352 -- >>> let r = cholesky big_K
353 -- >>> frobenius_norm (r - (transpose expected)) < 1e-12
356 cholesky :: forall m a. (Algebraic.C a, Arity m)
357 => (Mat m m a) -> (Mat m m a)
358 cholesky m = ifoldl2 f zero m
360 f :: Int -> Int -> (Mat m m a) -> a -> (Mat m m a)
361 f i j cur_R _ = set_idx cur_R (i,j) (r cur_R i j)
363 r :: (Mat m m a) -> Int -> Int -> a
365 | i == j = sqrt(m !!! (i,j) - sum [(cur_R !!! (k,i))^2 | k <- [0..i-1]])
366 | i < j = (((m !!! (i,j))
367 - sum [(cur_R !!! (k,i)) NP.* (cur_R !!! (k,j))
368 | k <- [0..i-1]]))/(cur_R !!! (i,i))
373 -- | Returns True if the given matrix is upper-triangular, and False
374 -- otherwise. The parameter @epsilon@ lets the caller choose a
379 -- >>> let m = fromList [[1,1],[1e-12,1]] :: Mat2 Double
380 -- >>> is_upper_triangular m
382 -- >>> is_upper_triangular' 1e-10 m
385 is_upper_triangular' :: forall m n a.
386 (Ord a, Ring.C a, Absolute.C a, Arity m, Arity n)
387 => a -- ^ The tolerance @epsilon@.
390 is_upper_triangular' epsilon matrix =
391 ifoldl2 f True matrix
393 f :: Int -> Int -> Bool -> a -> Bool
394 f _ _ False _ = False
397 -- use "less than or equal to" so zero is a valid epsilon
398 | otherwise = abs x <= epsilon
401 -- | Returns True if the given matrix is upper-triangular, and False
402 -- otherwise. We don't delegate to the general
403 -- 'is_upper_triangular'' here because it imposes additional
404 -- typeclass constraints throughout the library.
408 -- >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
409 -- >>> is_upper_triangular m
412 -- >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
413 -- >>> is_upper_triangular m
416 is_upper_triangular :: forall m n a.
417 (Eq a, Ring.C a, Arity m, Arity n)
419 is_upper_triangular matrix =
420 ifoldl2 f True matrix
422 f :: Int -> Int -> Bool -> a -> Bool
423 f _ _ False _ = False
430 -- | Returns True if the given matrix is lower-triangular, and False
435 -- >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
436 -- >>> is_lower_triangular m
439 -- >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
440 -- >>> is_lower_triangular m
443 is_lower_triangular :: (Eq a,
449 is_lower_triangular = is_upper_triangular . transpose
452 -- | Returns True if the given matrix is lower-triangular, and False
453 -- otherwise. The parameter @epsilon@ lets the caller choose a
458 -- >>> let m = fromList [[1,1e-12],[1,1]] :: Mat2 Double
459 -- >>> is_lower_triangular m
461 -- >>> is_lower_triangular' 1e-12 m
464 is_lower_triangular' :: (Ord a,
469 => a -- ^ The tolerance @epsilon@.
472 is_lower_triangular' epsilon = (is_upper_triangular' epsilon) . transpose
475 -- | Returns True if the given matrix is triangular, and False
480 -- >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
481 -- >>> is_triangular m
484 -- >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
485 -- >>> is_triangular m
488 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
489 -- >>> is_triangular m
492 is_triangular :: (Ord a,
499 is_triangular m = is_upper_triangular m || is_lower_triangular m
502 -- | Delete the @i@th row and @j@th column from the matrix. The name
503 -- \"preminor\" is made up, but is meant to signify that this is
504 -- usually used in the computationof a minor. A minor is simply the
505 -- determinant of a preminor in that case.
509 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
510 -- >>> preminor m 0 0 :: Mat2 Int
512 -- >>> preminor m 1 1 :: Mat2 Int
515 preminor :: (Arity m, Arity n)
520 preminor (Mat rows) i j = m
522 rows' = delete rows i
523 m = Mat $ V.map ((flip delete) j) rows'
526 -- | Compute the i,jth minor of a @matrix@.
530 -- >>> let m1 = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Double
534 minor :: (Arity m, Determined (Mat m m) a)
539 minor matrix i j = determinant (preminor matrix i j)
541 class (Eq a, Ring.C a) => Determined p a where
542 determinant :: (p a) -> a
544 instance (Eq a, Ring.C a) => Determined (Mat (S Z) (S Z)) a where
545 determinant = unscalar
551 Determined (Mat (S n) (S n)) a)
552 => Determined (Mat (S (S n)) (S (S n))) a where
553 -- | The recursive definition with a special-case for triangular matrices.
557 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
562 | is_triangular m = product [ m !!! (i,i) | i <- [0..(nrows m)-1] ]
563 | otherwise = determinant_recursive
567 determinant_recursive =
568 sum [ (-1)^(toInteger j) NP.* (m' 0 j) NP.* (minor m 0 j)
569 | j <- [0..(ncols m)-1] ]
573 -- | Matrix multiplication.
577 -- >>> let m1 = fromList [[1,2,3], [4,5,6]] :: Mat N2 N3 Int
578 -- >>> let m2 = fromList [[1,2],[3,4],[5,6]] :: Mat N3 N2 Int
583 (*) :: (Ring.C a, Arity m, Arity n, Arity p)
587 (*) m1 m2 = construct lambda
589 lambda i j = (transpose $ row m1 i) `dot` (column m2 j)
593 instance (Ring.C a, Arity m, Arity n) => Additive.C (Mat m n a) where
595 (Mat rows_one) + (Mat rows_two) =
596 Mat $ V.zipWith (V.zipWith (+)) rows_one rows_two
598 (Mat rows_one) - (Mat rows_two) =
599 Mat $ V.zipWith (V.zipWith (-)) rows_one rows_two
601 zero = Mat (V.replicate $ V.replicate (fromInteger 0))
604 instance (Ring.C a, Arity m, Arity n, m ~ n) => Ring.C (Mat (S m) (S n) a) where
605 -- The first * is ring multiplication, the second is matrix
607 one = identity_matrix
611 instance (Ring.C a, Arity m, Arity n) => Module.C a (Mat m n a) where
612 -- We can multiply a matrix by a scalar of the same type as its
614 x *> (Mat rows) = Mat $ V.map (V.map (NP.* x)) rows
617 instance (Absolute.C a,
621 => Normed (Col m a) where
622 -- | Generic p-norms for vectors in R^m that are represented as m-by-1
627 -- >>> let v1 = vec2d (3,4)
633 -- >>> let v1 = vec2d (-1,1) :: Col2 Double
634 -- >>> norm_p 1 v1 :: Double
637 -- >>> let v1 = vec0d :: Col0 Double
641 norm_p p (Mat rows) =
642 (root p') $ sum [fromRational' (toRational $ abs x)^p' | x <- xs]
645 xs = concat $ V.toList $ V.map V.toList rows
647 -- | The infinity norm.
651 -- >>> let v1 = vec3d (1,5,2)
655 norm_infty (Mat rows) =
656 fromRational' $ toRational
657 $ (V.foldl P.max 0) $ V.map (V.foldl P.max 0) rows
660 -- | Compute the Frobenius norm of a matrix. This essentially treats
661 -- the matrix as one long vector containing all of its entries (in
662 -- any order, it doesn't matter).
666 -- >>> let m = fromList [[1, 2, 3],[4,5,6],[7,8,9]] :: Mat3 Double
667 -- >>> frobenius_norm m == sqrt 285
670 -- >>> let m = fromList [[1, -1, 1],[-1,1,-1],[1,-1,1]] :: Mat3 Double
671 -- >>> frobenius_norm m == 3
674 frobenius_norm :: (Arity m, Arity n, Algebraic.C a, Ring.C a)
677 frobenius_norm matrix =
678 sqrt $ element_sum2 $ squares
680 squares = map2 (^2) matrix
683 -- Vector helpers. We want it to be easy to create low-dimension
684 -- column vectors, which are nx1 matrices.
686 -- | Convenient constructor for 2D vectors.
690 -- >>> import Roots.Simple
691 -- >>> let fst m = m !!! (0,0)
692 -- >>> let snd m = m !!! (1,0)
693 -- >>> let h = 0.5 :: Double
694 -- >>> let g1 m = 1.0 + h NP.* exp(-((fst m)^2))/(1.0 + (snd m)^2)
695 -- >>> let g2 m = 0.5 + h NP.* atan((fst m)^2 + (snd m)^2)
696 -- >>> let g u = vec2d ((g1 u), (g2 u))
697 -- >>> let u0 = vec2d (1.0, 1.0)
698 -- >>> let eps = 1/(10^9)
699 -- >>> fixed_point g eps u0
700 -- ((1.0728549599342185),(1.0820591495686167))
705 vec1d :: (a) -> Col1 a
706 vec1d (x) = Mat (mk1 (mk1 x))
708 vec2d :: (a,a) -> Col2 a
709 vec2d (x,y) = Mat (mk2 (mk1 x) (mk1 y))
711 vec3d :: (a,a,a) -> Col3 a
712 vec3d (x,y,z) = Mat (mk3 (mk1 x) (mk1 y) (mk1 z))
714 vec4d :: (a,a,a,a) -> Col4 a
715 vec4d (w,x,y,z) = Mat (mk4 (mk1 w) (mk1 x) (mk1 y) (mk1 z))
717 vec5d :: (a,a,a,a,a) -> Col5 a
718 vec5d (v,w,x,y,z) = Mat (mk5 (mk1 v) (mk1 w) (mk1 x) (mk1 y) (mk1 z))
721 -- Since we commandeered multiplication, we need to create 1x1
722 -- matrices in order to multiply things.
723 scalar :: a -> Mat1 a
724 scalar x = Mat (mk1 (mk1 x))
726 -- Get the scalar value out of a 1x1 matrix.
727 unscalar :: Mat1 a -> a
728 unscalar (Mat rows) = V.head $ V.head rows
731 dot :: (Ring.C a, Arity m)
735 v1 `dot` v2 = element_sum2 $ zipwith2 (NP.*) v1 v2
738 -- | The angle between @v1@ and @v2@ in Euclidean space.
742 -- >>> let v1 = vec2d (1.0, 0.0)
743 -- >>> let v2 = vec2d (0.0, 1.0)
744 -- >>> angle v1 v2 == pi/2.0
747 angle :: (Transcendental.C a,
758 theta = (recip norms) NP.* (v1 `dot` v2)
759 norms = (norm v1) NP.* (norm v2)
762 -- | Retrieve the diagonal elements of the given matrix as a \"column
763 -- vector,\" i.e. a m-by-1 matrix. We require the matrix to be
764 -- square to avoid ambiguity in the return type which would ideally
765 -- have dimension min(m,n) supposing an m-by-n matrix.
769 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
773 diagonal :: (Arity m) => Mat m m a -> Col m a
777 lambda i _ = matrix !!! (i,i)
780 -- | Given a square @matrix@, return a new matrix of the same size
781 -- containing only the on-diagonal entries of @matrix@. The
782 -- off-diagonal entries are set to zero.
786 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
787 -- >>> diagonal_part m
788 -- ((1,0,0),(0,5,0),(0,0,9))
790 diagonal_part :: (Arity m, Ring.C a)
793 diagonal_part matrix =
796 lambda i j = if i == j then matrix !!! (i,j) else 0
799 -- | Given a square @matrix@, return a new matrix of the same size
800 -- containing only the on-diagonal and below-diagonal entries of
801 -- @matrix@. The above-diagonal entries are set to zero.
805 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
807 -- ((1,0,0),(4,5,0),(7,8,9))
809 lt_part :: (Arity m, Ring.C a)
815 lambda i j = if i >= j then matrix !!! (i,j) else 0
818 -- | Given a square @matrix@, return a new matrix of the same size
819 -- containing only the below-diagonal entries of @matrix@. The on-
820 -- and above-diagonal entries are set to zero.
824 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
825 -- >>> lt_part_strict m
826 -- ((0,0,0),(4,0,0),(7,8,0))
828 lt_part_strict :: (Arity m, Ring.C a)
831 lt_part_strict matrix =
834 lambda i j = if i > j then matrix !!! (i,j) else 0
837 -- | Given a square @matrix@, return a new matrix of the same size
838 -- containing only the on-diagonal and above-diagonal entries of
839 -- @matrix@. The below-diagonal entries are set to zero.
843 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
845 -- ((1,2,3),(0,5,6),(0,0,9))
847 ut_part :: (Arity m, Ring.C a)
850 ut_part = transpose . lt_part . transpose
853 -- | Given a square @matrix@, return a new matrix of the same size
854 -- containing only the above-diagonal entries of @matrix@. The on-
855 -- and below-diagonal entries are set to zero.
859 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
860 -- >>> ut_part_strict m
861 -- ((0,2,3),(0,0,6),(0,0,0))
863 ut_part_strict :: (Arity m, Ring.C a)
866 ut_part_strict = transpose . lt_part_strict . transpose
869 -- | Compute the trace of a square matrix, the sum of the elements
870 -- which lie on its diagonal. We require the matrix to be
871 -- square to avoid ambiguity in the return type which would ideally
872 -- have dimension min(m,n) supposing an m-by-n matrix.
876 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
880 trace :: (Arity m, Ring.C a) => Mat m m a -> a
881 trace = element_sum2 . diagonal
885 -- | Zip together two matrices.
887 -- TODO: don't cheat with construct (map V.zips instead).
891 -- >>> let m1 = fromList [[1],[1],[1]] :: Col3 Int
892 -- >>> let m2 = fromList [[1],[2],[3]] :: Col3 Int
894 -- (((1,1)),((1,2)),((1,3)))
896 -- >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int
897 -- >>> let m2 = fromList [[1,1],[1,1]] :: Mat2 Int
899 -- (((1,1),(2,1)),((3,1),(4,1)))
901 zip2 :: (Arity m, Arity n) => Mat m n a -> Mat m n b -> Mat m n (a,b)
905 lambda i j = (m1 !!! (i,j), m2 !!! (i,j))
908 -- | Zip together three matrices.
910 -- TODO: don't cheat with construct (map V.zips instead).
914 -- >>> let m1 = fromList [[1],[1],[1]] :: Col3 Int
915 -- >>> let m2 = fromList [[1],[2],[3]] :: Col3 Int
916 -- >>> let m3 = fromList [[4],[5],[6]] :: Col3 Int
917 -- >>> zip2three m1 m2 m3
918 -- (((1,1,4)),((1,2,5)),((1,3,6)))
920 -- >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int
921 -- >>> let m2 = fromList [[1,1],[1,1]] :: Mat2 Int
922 -- >>> let m3 = fromList [[8,2],[6,3]] :: Mat2 Int
923 -- >>> zip2three m1 m2 m3
924 -- (((1,1,8),(2,1,2)),((3,1,6),(4,1,3)))
926 zip2three :: (Arity m, Arity n)
934 lambda i j = (m1 !!! (i,j), m2 !!! (i,j), m3 !!! (i,j))
937 -- | Zip together two matrices using the supplied function.
941 -- >>> let c1 = fromList [[1],[2],[3]] :: Col3 Integer
942 -- >>> let c2 = fromList [[4],[5],[6]] :: Col3 Integer
943 -- >>> zipwith2 (^) c1 c2
946 zipwith2 :: (Arity m, Arity n)
954 lambda i j = f (c1 !!! (i,j)) (c2 !!! (i,j))
957 -- | Map a function over a matrix of any dimensions.
961 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
965 map2 :: (a -> b) -> Mat m n a -> Mat m n b
972 -- | Fold over the entire matrix passing the coordinates @i@ and @j@
973 -- (of the row/column) to the accumulation function. The fold occurs
974 -- from top-left to bottom-right.
978 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
979 -- >>> ifoldl2 (\i j cur _ -> cur + i + j) 0 m
982 ifoldl2 :: forall a b m n.
983 (Int -> Int -> b -> a -> b)
987 ifoldl2 f initial (Mat rows) =
988 V.ifoldl row_function initial rows
990 -- | The order that we need this in (so that @g idx@ makes sense)
991 -- is a little funny. So that we don't need to pass weird
992 -- functions into ifoldl2, we swap the second and third
993 -- arguments of @f@ calling the result @g@.
994 g :: Int -> b -> Int -> a -> b
997 row_function :: b -> Int -> Vec n a -> b
998 row_function rowinit idx r = V.ifoldl (g idx) rowinit r
1001 -- | Left fold over the entries of a matrix (top-left to bottom-right).
1003 foldl2 :: forall a b m n.
1008 foldl2 f initial matrix =
1009 -- Use the index fold but ignore the index arguments.
1010 let g _ _ = f in ifoldl2 g initial matrix
1013 -- | Fold over the entire matrix passing the coordinates @i@ and @j@
1014 -- (of the row/column) to the accumulation function. The fold occurs
1015 -- from bottom-right to top-left.
1017 -- The order of the arguments in the supplied function are different
1018 -- from those in V.ifoldr; we keep them similar to ifoldl2.
1022 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
1023 -- >>> ifoldr2 (\i j cur _ -> cur + i + j) 0 m
1026 ifoldr2 :: forall a b m n.
1027 (Int -> Int -> b -> a -> b)
1031 ifoldr2 f initial (Mat rows) =
1032 V.ifoldr row_function initial rows
1034 -- | Swap the order of arguments in @f@ so that it agrees with the
1035 -- @f@ passed to ifoldl2.
1036 g :: Int -> Int -> a -> b -> b
1037 g w x y z = f w x z y
1039 row_function :: Int -> Vec n a -> b -> b
1040 row_function idx r rowinit = V.ifoldr (g idx) rowinit r
1043 -- | Map a function over a matrix of any dimensions, passing the
1044 -- coordinates @i@ and @j@ to the function @f@.
1048 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
1049 -- >>> imap2 (\i j _ -> i+j) m
1052 imap2 :: (Int -> Int -> a -> b) -> Mat m n a -> Mat m n b
1053 imap2 f (Mat rows) =
1059 -- | Reverse the order of elements in a matrix.
1063 -- >>> let m1 = fromList [[1,2,3]] :: Row3 Int
1067 -- >>> let m1 = vec3d (1,2,3 :: Int)
1071 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
1073 -- ((9,8,7),(6,5,4),(3,2,1))
1075 reverse2 :: (Arity m, Arity n) => Mat m n a -> Mat m n a
1076 reverse2 (Mat rows) = Mat $ V.reverse $ V.map V.reverse rows
1079 -- | Unsafely set the (i,j) element of the given matrix.
1083 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
1084 -- >>> set_idx m (1,1) 17
1085 -- ((1,2,3),(4,17,6),(7,8,9))
1087 set_idx :: forall m n a.
1093 set_idx matrix (i,j) newval =
1094 imap2 updater matrix
1096 updater :: Int -> Int -> a -> a
1097 updater k l existing =
1103 -- | Compute the i,jth cofactor of the given @matrix@. This simply
1104 -- premultiplues the i,jth minor by (-1)^(i+j).
1105 cofactor :: (Arity m, Determined (Mat m m) a)
1106 => Mat (S m) (S m) a
1110 cofactor matrix i j =
1111 (-1)^(toInteger i + toInteger j) NP.* (minor matrix i j)
1114 -- | Compute the inverse of a matrix using cofactor expansion
1115 -- (generalized Cramer's rule).
1119 -- >>> let m1 = fromList [[37,22],[17,54]] :: Mat2 Double
1120 -- >>> let e1 = [54/1624, -22/1624] :: [Double]
1121 -- >>> let e2 = [-17/1624, 37/1624] :: [Double]
1122 -- >>> let expected = fromList [e1, e2] :: Mat2 Double
1123 -- >>> let actual = inverse m1
1124 -- >>> frobenius_norm (actual - expected) < 1e-12
1127 inverse :: (Arity m,
1128 Determined (Mat (S m) (S m)) a,
1129 Determined (Mat m m) a,
1131 => Mat (S m) (S m) a
1132 -> Mat (S m) (S m) a
1134 (1 / (determinant matrix)) *> (transpose $ construct lambda)
1136 lambda i j = cofactor matrix i j
1140 -- | Retrieve the rows of a matrix as a column matrix. If the given
1141 -- matrix is m-by-n, the result would be an m-by-1 column whose
1142 -- entries are 1-by-n row matrices.
1146 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
1147 -- >>> (rows2 m) !!! (0,0)
1149 -- >>> (rows2 m) !!! (1,0)
1152 rows2 :: (Arity m, Arity n)
1156 Mat $ V.map (mk1. Mat . mk1) rows
1160 -- | Sum the elements of a matrix.
1164 -- >>> let m = fromList [[1,-1],[3,4]] :: Mat2 Int
1165 -- >>> element_sum2 m
1168 element_sum2 :: (Arity m, Arity n, Additive.C a) => Mat m n a -> a
1169 element_sum2 = foldl2 (+) zero