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 (
35 import qualified Data.Vector.Fixed as V (
48 import Data.Vector.Fixed.Cont ( Arity, arity )
49 import Linear.Vector ( Vec, delete, element_sum )
50 import Normed ( Normed(..) )
52 import NumericPrelude hiding ( (*), abs )
53 import qualified NumericPrelude as NP ( (*) )
54 import qualified Algebra.Absolute as Absolute ( C )
55 import Algebra.Absolute ( abs )
56 import qualified Algebra.Additive as Additive ( C )
57 import qualified Algebra.Algebraic as Algebraic ( C )
58 import Algebra.Algebraic ( root )
59 import qualified Algebra.Ring as Ring ( C )
60 import qualified Algebra.Module as Module ( C )
61 import qualified Algebra.RealRing as RealRing ( C )
62 import qualified Algebra.ToRational as ToRational ( C )
63 import qualified Algebra.Transcendental as Transcendental ( C )
64 import qualified Prelude as P ( map )
66 -- | Our main matrix type.
67 data Mat m n a = (Arity m, Arity n) => Mat (Vec m (Vec n a))
69 -- Type synonyms for n-by-n matrices.
70 type Mat1 a = Mat N1 N1 a
71 type Mat2 a = Mat N2 N2 a
72 type Mat3 a = Mat N3 N3 a
73 type Mat4 a = Mat N4 N4 a
74 type Mat5 a = Mat N5 N5 a
76 -- | Type synonym for row vectors expressed as 1-by-n matrices.
77 type Row n a = Mat N1 n a
79 -- Type synonyms for 1-by-n row "vectors".
80 type Row1 a = Row N1 a
81 type Row2 a = Row N2 a
82 type Row3 a = Row N3 a
83 type Row4 a = Row N4 a
84 type Row5 a = Row N5 a
86 -- | Type synonym for column vectors expressed as n-by-1 matrices.
87 type Col n a = Mat n N1 a
89 -- Type synonyms for n-by-1 column "vectors".
90 type Col1 a = Col N1 a
91 type Col2 a = Col N2 a
92 type Col3 a = Col N3 a
93 type Col4 a = Col N4 a
94 type Col5 a = Col N5 a
96 -- We need a big column for Gaussian quadrature.
97 type N10 = S (S (S (S (S N5))))
98 type Col10 a = Col N10 a
101 instance (Eq a) => Eq (Mat m n a) where
102 -- | Compare a row at a time.
106 -- >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int
107 -- >>> let m2 = fromList [[1,2],[3,4]] :: Mat2 Int
108 -- >>> let m3 = fromList [[5,6],[7,8]] :: Mat2 Int
114 (Mat rows1) == (Mat rows2) =
115 V.and $ V.zipWith comp rows1 rows2
117 -- Compare a row, one column at a time.
118 comp row1 row2 = V.and (V.zipWith (==) row1 row2)
121 instance (Show a) => Show (Mat m n a) where
122 -- | Display matrices and vectors as ordinary tuples. This is poor
123 -- practice, but these results are primarily displayed
124 -- interactively and convenience trumps correctness (said the guy
125 -- who insists his vector lengths be statically checked at
130 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
135 "(" ++ (intercalate "," (V.toList row_strings)) ++ ")"
137 row_strings = V.map show_vector rows
139 "(" ++ (intercalate "," element_strings) ++ ")"
142 element_strings = P.map show v1l
145 -- | Convert a matrix to a nested list.
146 toList :: Mat m n a -> [[a]]
147 toList (Mat rows) = map V.toList (V.toList rows)
150 -- | Create a matrix from a nested list.
151 fromList :: (Arity m, Arity n) => [[a]] -> Mat m n a
152 fromList vs = Mat (V.fromList $ map V.fromList vs)
155 -- | Unsafe indexing. Much faster than the safe indexing.
156 (!!!) :: (Arity m, Arity n) => Mat m n a -> (Int, Int) -> a
157 (!!!) (Mat rows) (i, j) = (rows ! i) ! j
164 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
192 (!!?) :: (Arity m, Arity n) => Mat m n a -> (Int, Int) -> Maybe a
194 ifoldl2 f Nothing matrix
196 f k l found cur = if (k,l) == idx then (Just cur) else found
199 -- | The number of rows in the matrix.
200 nrows :: forall m n a. (Arity m) => Mat m n a -> Int
201 nrows _ = arity (undefined :: m)
204 -- | The number of columns in the first row of the
205 -- matrix. Implementation stolen from Data.Vector.Fixed.length.
206 ncols :: forall m n a. (Arity n) => Mat m n a -> Int
207 ncols _ = arity (undefined :: n)
210 -- | Return the @i@th row of @m@ as a matrix. Unsafe.
211 row :: (Arity m, Arity n) => Mat m n a -> Int -> Row n a
215 lambda _ j = m !!! (i, j)
218 -- | Return the @j@th column of @m@ as a matrix. Unsafe.
219 column :: (Arity m, Arity n) => Mat m n a -> Int -> Col m a
223 lambda i _ = m !!! (i, j)
226 -- | Transpose @m@; switch it's columns and its rows. This is a dirty
227 -- implementation, but I don't see a better way.
229 -- TODO: Don't cheat with fromList.
233 -- >>> let m = fromList [[1,2], [3,4]] :: Mat2 Int
237 transpose :: (Arity m, Arity n) => Mat m n a -> Mat n m a
241 lambda i j = matrix !!! (j,i)
244 -- | Is @m@ symmetric?
248 -- >>> let m1 = fromList [[1,2], [2,1]] :: Mat2 Int
252 -- >>> let m2 = fromList [[1,2], [3,1]] :: Mat2 Int
256 symmetric :: (Eq a, Arity m) => Mat m m a -> Bool
261 -- | Construct a new matrix from a function @lambda@. The function
262 -- @lambda@ should take two parameters i,j corresponding to the
263 -- entries in the matrix. The i,j entry of the resulting matrix will
264 -- have the value returned by lambda i j.
268 -- >>> let lambda i j = i + j
269 -- >>> construct lambda :: Mat3 Int
270 -- ((0,1,2),(1,2,3),(2,3,4))
272 construct :: forall m n a. (Arity m, Arity n)
273 => (Int -> Int -> a) -> Mat m n a
274 construct lambda = Mat $ generate make_row
276 make_row :: Int -> Vec n a
277 make_row i = generate (lambda i)
280 -- | Create an identity matrix with the right dimensions.
284 -- >>> identity_matrix :: Mat3 Int
285 -- ((1,0,0),(0,1,0),(0,0,1))
286 -- >>> identity_matrix :: Mat3 Double
287 -- ((1.0,0.0,0.0),(0.0,1.0,0.0),(0.0,0.0,1.0))
289 identity_matrix :: (Arity m, Ring.C a) => Mat m m a
291 construct (\i j -> if i == j then (fromInteger 1) else (fromInteger 0))
294 -- | Given a positive-definite matrix @m@, computes the
295 -- upper-triangular matrix @r@ with (transpose r)*r == m and all
296 -- values on the diagonal of @r@ positive.
300 -- >>> let m1 = fromList [[20,-1], [-1,20]] :: Mat2 Double
301 -- >>> let r = cholesky m1
302 -- >>> frobenius_norm ((transpose r)*r - m1) < 1e-10
304 -- >>> is_upper_triangular r
307 -- >>> import Naturals ( N7 )
308 -- >>> let k1 = [6, -3, 0, 0, 0, 0, 0] :: [Double]
309 -- >>> let k2 = [-3, 10.5, -7.5, 0, 0, 0, 0] :: [Double]
310 -- >>> let k3 = [0, -7.5, 12.5, 0, 0, 0, 0] :: [Double]
311 -- >>> let k4 = [0, 0, 0, 6, 0, 0, 0] :: [Double]
312 -- >>> let k5 = [0, 0, 0, 0, 6, 0, 0] :: [Double]
313 -- >>> let k6 = [0, 0, 0, 0, 0, 6, 0] :: [Double]
314 -- >>> let k7 = [0, 0, 0, 0, 0, 0, 15] :: [Double]
315 -- >>> let big_K = fromList [k1,k2,k3,k4,k5,k6,k7] :: Mat N7 N7 Double
317 -- >>> let e1 = [2.449489742783178,0,0,0,0,0,0] :: [Double]
318 -- >>> let e2 = [-1.224744871391589,3,0,0,0,0,0] :: [Double]
319 -- >>> let e3 = [0,-5/2,5/2,0,0,0,0] :: [Double]
320 -- >>> let e4 = [0,0,0,2.449489742783178,0,0,0] :: [Double]
321 -- >>> let e5 = [0,0,0,0,2.449489742783178,0,0] :: [Double]
322 -- >>> let e6 = [0,0,0,0,0,2.449489742783178,0] :: [Double]
323 -- >>> let e7 = [0,0,0,0,0,0,3.872983346207417] :: [Double]
324 -- >>> let expected = fromList [e1,e2,e3,e4,e5,e6,e7] :: Mat N7 N7 Double
326 -- >>> let r = cholesky big_K
327 -- >>> frobenius_norm (r - (transpose expected)) < 1e-12
330 cholesky :: forall m n a. (Algebraic.C a, Arity m, Arity n)
331 => (Mat m n a) -> (Mat m n a)
332 cholesky m = construct r
335 r i j | i == j = sqrt(m !!! (i,j) - sum [(r k i)^2 | k <- [0..i-1]])
337 (((m !!! (i,j)) - sum [(r k i) NP.* (r k j) | k <- [0..i-1]]))/(r i i)
341 -- | Returns True if the given matrix is upper-triangular, and False
342 -- otherwise. The parameter @epsilon@ lets the caller choose a
347 -- >>> let m = fromList [[1,1],[1e-12,1]] :: Mat2 Double
348 -- >>> is_upper_triangular m
350 -- >>> is_upper_triangular' 1e-10 m
353 is_upper_triangular' :: forall m n a.
354 (Ord a, Ring.C a, Absolute.C a, Arity m, Arity n)
355 => a -- ^ The tolerance @epsilon@.
358 is_upper_triangular' epsilon matrix =
359 ifoldl2 f True matrix
361 f :: Int -> Int -> Bool -> a -> Bool
362 f _ _ False _ = False
365 -- use "less than or equal to" so zero is a valid epsilon
366 | otherwise = abs x <= epsilon
369 -- | Returns True if the given matrix is upper-triangular, and False
370 -- otherwise. We don't delegate to the general
371 -- 'is_upper_triangular'' here because it imposes additional
372 -- typeclass constraints throughout the library.
376 -- >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
377 -- >>> is_upper_triangular m
380 -- >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
381 -- >>> is_upper_triangular m
384 is_upper_triangular :: forall m n a.
385 (Eq a, Ring.C a, Arity m, Arity n)
387 is_upper_triangular matrix =
388 ifoldl2 f True matrix
390 f :: Int -> Int -> Bool -> a -> Bool
391 f _ _ False _ = False
398 -- | Returns True if the given matrix is lower-triangular, and False
403 -- >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
404 -- >>> is_lower_triangular m
407 -- >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
408 -- >>> is_lower_triangular m
411 is_lower_triangular :: (Eq a,
417 is_lower_triangular = is_upper_triangular . transpose
420 -- | Returns True if the given matrix is lower-triangular, and False
421 -- otherwise. The parameter @epsilon@ lets the caller choose a
426 -- >>> let m = fromList [[1,1e-12],[1,1]] :: Mat2 Double
427 -- >>> is_lower_triangular m
429 -- >>> is_lower_triangular' 1e-12 m
432 is_lower_triangular' :: (Ord a,
437 => a -- ^ The tolerance @epsilon@.
440 is_lower_triangular' epsilon = (is_upper_triangular' epsilon) . transpose
443 -- | Returns True if the given matrix is triangular, and False
448 -- >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
449 -- >>> is_triangular m
452 -- >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
453 -- >>> is_triangular m
456 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
457 -- >>> is_triangular m
460 is_triangular :: (Ord a,
467 is_triangular m = is_upper_triangular m || is_lower_triangular m
470 -- | Return the (i,j)th minor of m.
474 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
475 -- >>> minor m 0 0 :: Mat2 Int
477 -- >>> minor m 1 1 :: Mat2 Int
488 minor (Mat rows) i j = m
490 rows' = delete rows i
491 m = Mat $ V.map ((flip delete) j) rows'
494 class (Eq a, Ring.C a) => Determined p a where
495 determinant :: (p a) -> a
497 instance (Eq a, Ring.C a) => Determined (Mat (S Z) (S Z)) a where
498 determinant (Mat rows) = (V.head . V.head) rows
504 Determined (Mat (S n) (S n)) a)
505 => Determined (Mat (S (S n)) (S (S n))) a where
506 -- | The recursive definition with a special-case for triangular matrices.
510 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
515 | is_triangular m = product [ m !!! (i,i) | i <- [0..(nrows m)-1] ]
516 | otherwise = determinant_recursive
520 det_minor i j = determinant (minor m i j)
522 determinant_recursive =
523 sum [ (-1)^(toInteger j) NP.* (m' 0 j) NP.* (det_minor 0 j)
524 | j <- [0..(ncols m)-1] ]
528 -- | Matrix multiplication.
532 -- >>> let m1 = fromList [[1,2,3], [4,5,6]] :: Mat N2 N3 Int
533 -- >>> let m2 = fromList [[1,2],[3,4],[5,6]] :: Mat N3 N2 Int
538 (*) :: (Ring.C a, Arity m, Arity n, Arity p)
542 (*) m1 m2 = construct lambda
545 sum [(m1 !!! (i,k)) NP.* (m2 !!! (k,j)) | k <- [0..(ncols m1)-1] ]
549 instance (Ring.C a, Arity m, Arity n) => Additive.C (Mat m n a) where
551 (Mat rows1) + (Mat rows2) =
552 Mat $ V.zipWith (V.zipWith (+)) rows1 rows2
554 (Mat rows1) - (Mat rows2) =
555 Mat $ V.zipWith (V.zipWith (-)) rows1 rows2
557 zero = Mat (V.replicate $ V.replicate (fromInteger 0))
560 instance (Ring.C a, Arity m, Arity n, m ~ n) => Ring.C (Mat m n a) where
561 -- The first * is ring multiplication, the second is matrix
566 instance (Ring.C a, Arity m, Arity n) => Module.C a (Mat m n a) where
567 -- We can multiply a matrix by a scalar of the same type as its
569 x *> (Mat rows) = Mat $ V.map (V.map (NP.* x)) rows
572 instance (Algebraic.C a,
575 => Normed (Mat (S m) N1 a) where
576 -- | Generic p-norms for vectors in R^n that are represented as nx1
581 -- >>> let v1 = vec2d (3,4)
587 norm_p p (Mat rows) =
588 (root p') $ sum [fromRational' (toRational x)^p' | x <- xs]
591 xs = concat $ V.toList $ V.map V.toList rows
593 -- | The infinity norm.
597 -- >>> let v1 = vec3d (1,5,2)
601 norm_infty (Mat rows) =
602 fromRational' $ toRational $ V.maximum $ V.map V.maximum rows
605 -- | Compute the Frobenius norm of a matrix. This essentially treats
606 -- the matrix as one long vector containing all of its entries (in
607 -- any order, it doesn't matter).
611 -- >>> let m = fromList [[1, 2, 3],[4,5,6],[7,8,9]] :: Mat3 Double
612 -- >>> frobenius_norm m == sqrt 285
615 -- >>> let m = fromList [[1, -1, 1],[-1,1,-1],[1,-1,1]] :: Mat3 Double
616 -- >>> frobenius_norm m == 3
619 frobenius_norm :: (Algebraic.C a, Ring.C a) => Mat m n a -> a
620 frobenius_norm (Mat rows) =
621 sqrt $ element_sum $ V.map row_sum rows
623 -- | Square and add up the entries of a row.
624 row_sum = element_sum . V.map (^2)
627 -- Vector helpers. We want it to be easy to create low-dimension
628 -- column vectors, which are nx1 matrices.
630 -- | Convenient constructor for 2D vectors.
634 -- >>> import Roots.Simple
635 -- >>> let fst m = m !!! (0,0)
636 -- >>> let snd m = m !!! (1,0)
637 -- >>> let h = 0.5 :: Double
638 -- >>> let g1 m = 1.0 + h NP.* exp(-((fst m)^2))/(1.0 + (snd m)^2)
639 -- >>> let g2 m = 0.5 + h NP.* atan((fst m)^2 + (snd m)^2)
640 -- >>> let g u = vec2d ((g1 u), (g2 u))
641 -- >>> let u0 = vec2d (1.0, 1.0)
642 -- >>> let eps = 1/(10^9)
643 -- >>> fixed_point g eps u0
644 -- ((1.0728549599342185),(1.0820591495686167))
646 vec1d :: (a) -> Col1 a
647 vec1d (x) = Mat (mk1 (mk1 x))
649 vec2d :: (a,a) -> Col2 a
650 vec2d (x,y) = Mat (mk2 (mk1 x) (mk1 y))
652 vec3d :: (a,a,a) -> Col3 a
653 vec3d (x,y,z) = Mat (mk3 (mk1 x) (mk1 y) (mk1 z))
655 vec4d :: (a,a,a,a) -> Col4 a
656 vec4d (w,x,y,z) = Mat (mk4 (mk1 w) (mk1 x) (mk1 y) (mk1 z))
658 vec5d :: (a,a,a,a,a) -> Col5 a
659 vec5d (v,w,x,y,z) = Mat (mk5 (mk1 v) (mk1 w) (mk1 x) (mk1 y) (mk1 z))
661 -- Since we commandeered multiplication, we need to create 1x1
662 -- matrices in order to multiply things.
663 scalar :: a -> Mat1 a
664 scalar x = Mat (mk1 (mk1 x))
666 -- Get the scalar value out of a 1x1 matrix.
667 unscalar :: Mat1 a -> a
668 unscalar (Mat rows) = V.head $ V.head rows
671 dot :: (Ring.C a, Arity m)
675 v1 `dot` v2 = unscalar $ ((transpose v1) * v2)
678 -- | The angle between @v1@ and @v2@ in Euclidean space.
682 -- >>> let v1 = vec2d (1.0, 0.0)
683 -- >>> let v2 = vec2d (0.0, 1.0)
684 -- >>> angle v1 v2 == pi/2.0
687 angle :: (Transcendental.C a,
698 theta = (recip norms) NP.* (v1 `dot` v2)
699 norms = (norm v1) NP.* (norm v2)
702 -- | Retrieve the diagonal elements of the given matrix as a \"column
703 -- vector,\" i.e. a m-by-1 matrix. We require the matrix to be
704 -- square to avoid ambiguity in the return type which would ideally
705 -- have dimension min(m,n) supposing an m-by-n matrix.
709 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
713 diagonal :: (Arity m) => Mat m m a -> Col m a
717 lambda i _ = matrix !!! (i,i)
720 -- | Given a square @matrix@, return a new matrix of the same size
721 -- containing only the on-diagonal entries of @matrix@. The
722 -- off-diagonal entries are set to zero.
726 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
727 -- >>> diagonal_part m
728 -- ((1,0,0),(0,5,0),(0,0,9))
730 diagonal_part :: (Arity m, Ring.C a)
733 diagonal_part matrix =
736 lambda i j = if i == j then matrix !!! (i,j) else 0
739 -- | Given a square @matrix@, return a new matrix of the same size
740 -- containing only the on-diagonal and below-diagonal entries of
741 -- @matrix@. The above-diagonal entries are set to zero.
745 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
747 -- ((1,0,0),(4,5,0),(7,8,9))
749 lt_part :: (Arity m, Ring.C a)
755 lambda i j = if i >= j then matrix !!! (i,j) else 0
758 -- | Given a square @matrix@, return a new matrix of the same size
759 -- containing only the below-diagonal entries of @matrix@. The on-
760 -- and above-diagonal entries are set to zero.
764 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
765 -- >>> lt_part_strict m
766 -- ((0,0,0),(4,0,0),(7,8,0))
768 lt_part_strict :: (Arity m, Ring.C a)
771 lt_part_strict matrix =
774 lambda i j = if i > j then matrix !!! (i,j) else 0
777 -- | Given a square @matrix@, return a new matrix of the same size
778 -- containing only the on-diagonal and above-diagonal entries of
779 -- @matrix@. The below-diagonal entries are set to zero.
783 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
785 -- ((1,2,3),(0,5,6),(0,0,9))
787 ut_part :: (Arity m, Ring.C a)
790 ut_part = transpose . lt_part . transpose
793 -- | Given a square @matrix@, return a new matrix of the same size
794 -- containing only the above-diagonal entries of @matrix@. The on-
795 -- and below-diagonal entries are set to zero.
799 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
800 -- >>> ut_part_strict m
801 -- ((0,2,3),(0,0,6),(0,0,0))
803 ut_part_strict :: (Arity m, Ring.C a)
806 ut_part_strict = transpose . lt_part_strict . transpose
809 -- | Compute the trace of a square matrix, the sum of the elements
810 -- which lie on its diagonal. We require the matrix to be
811 -- square to avoid ambiguity in the return type which would ideally
812 -- have dimension min(m,n) supposing an m-by-n matrix.
816 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
820 trace :: (Arity m, Ring.C a) => Mat m m a -> a
822 let (Mat rows) = diagonal matrix
824 element_sum $ V.map V.head rows
827 -- | Zip together two matrices.
829 -- TODO: don't cheat with construct (map V.zips instead).
833 -- >>> let m1 = fromList [[1],[1],[1]] :: Col3 Int
834 -- >>> let m2 = fromList [[1],[2],[3]] :: Col3 Int
836 -- (((1,1)),((1,2)),((1,3)))
838 -- >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int
839 -- >>> let m2 = fromList [[1,1],[1,1]] :: Mat2 Int
841 -- (((1,1),(2,1)),((3,1),(4,1)))
843 zip2 :: (Arity m, Arity n) => Mat m n a -> Mat m n a -> Mat m n (a,a)
847 lambda i j = (m1 !!! (i,j), m2 !!! (i,j))
850 -- | Zip together three matrices.
852 -- TODO: don't cheat with construct (map V.zips instead).
856 -- >>> let m1 = fromList [[1],[1],[1]] :: Col3 Int
857 -- >>> let m2 = fromList [[1],[2],[3]] :: Col3 Int
858 -- >>> let m3 = fromList [[4],[5],[6]] :: Col3 Int
859 -- >>> zip2three m1 m2 m3
860 -- (((1,1,4)),((1,2,5)),((1,3,6)))
862 -- >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int
863 -- >>> let m2 = fromList [[1,1],[1,1]] :: Mat2 Int
864 -- >>> let m3 = fromList [[8,2],[6,3]] :: Mat2 Int
865 -- >>> zip2three m1 m2 m3
866 -- (((1,1,8),(2,1,2)),((3,1,6),(4,1,3)))
868 zip2three :: (Arity m, Arity n)
876 lambda i j = (m1 !!! (i,j), m2 !!! (i,j), m3 !!! (i,j))
879 -- | Zip together two matrices using the supplied function.
883 -- >>> let c1 = fromList [[1],[2],[3]] :: Col3 Integer
884 -- >>> let c2 = fromList [[4],[5],[6]] :: Col3 Integer
885 -- >>> zipwith2 (^) c1 c2
896 lambda i j = f (c1 !!! (i,j)) (c2 !!! (i,j))
899 -- | Map a function over a matrix of any dimensions.
903 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
907 map2 :: (a -> b) -> Mat m n a -> Mat m n b
914 -- | Fold over the entire matrix passing the coordinates @i@ and @j@
915 -- (of the row/column) to the accumulation function. The fold occurs
916 -- from top-left to bottom-right.
920 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
921 -- >>> ifoldl2 (\i j cur _ -> cur + i + j) 0 m
924 ifoldl2 :: forall a b m n.
925 (Int -> Int -> b -> a -> b)
929 ifoldl2 f initial (Mat rows) =
930 V.ifoldl row_function initial rows
932 -- | The order that we need this in (so that @g idx@ makes sense)
933 -- is a little funny. So that we don't need to pass weird
934 -- functions into ifoldl2, we swap the second and third
935 -- arguments of @f@ calling the result @g@.
936 g :: Int -> b -> Int -> a -> b
939 row_function :: b -> Int -> Vec n a -> b
940 row_function rowinit idx r = V.ifoldl (g idx) rowinit r
943 -- | Fold over the entire matrix passing the coordinates @i@ and @j@
944 -- (of the row/column) to the accumulation function. The fold occurs
945 -- from bottom-right to top-left.
947 -- The order of the arguments in the supplied function are different
948 -- from those in V.ifoldr; we keep them similar to ifoldl2.
952 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
953 -- >>> ifoldr2 (\i j cur _ -> cur + i + j) 0 m
956 ifoldr2 :: forall a b m n.
957 (Int -> Int -> b -> a -> b)
961 ifoldr2 f initial (Mat rows) =
962 V.ifoldr row_function initial rows
964 -- | Swap the order of arguments in @f@ so that it agrees with the
965 -- @f@ passed to ifoldl2.
966 g :: Int -> Int -> a -> b -> b
967 g w x y z = f w x z y
969 row_function :: Int -> Vec n a -> b -> b
970 row_function idx r rowinit = V.ifoldr (g idx) rowinit r
973 -- | Map a function over a matrix of any dimensions, passing the
974 -- coordinates @i@ and @j@ to the function @f@.
978 -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
979 -- >>> imap2 (\i j _ -> i+j) m
982 imap2 :: (Int -> Int -> a -> b) -> Mat m n a -> Mat m n b
989 -- | Reverse the order of elements in a matrix.
993 -- >>> let m1 = fromList [[1,2,3]] :: Row3 Int
997 -- >>> let m1 = vec3d (1,2,3 :: Int)
1001 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
1003 -- ((9,8,7),(6,5,4),(3,2,1))
1005 reverse2 :: (Arity m, Arity n) => Mat m n a -> Mat m n a
1006 reverse2 (Mat rows) = Mat $ V.reverse $ V.map V.reverse rows
1009 -- | Unsafely set the (i,j) element of the given matrix.
1013 -- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
1014 -- >>> set_idx m (1,1) 17
1015 -- ((1,2,3),(4,17,6),(7,8,9))
1017 set_idx :: forall m n a.
1023 set_idx matrix (i,j) newval =
1024 imap2 updater matrix
1026 updater :: Int -> Int -> a -> a
1027 updater k l existing =