+
+
+-- | Given a square @matrix@, return a new matrix of the same size
+-- containing only the on-diagonal and below-diagonal entries of
+-- @matrix@. The above-diagonal entries are set to zero.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
+-- >>> lt_part m
+-- ((1,0,0),(4,5,0),(7,8,9))
+--
+lt_part :: (Arity m, Ring.C a)
+ => Mat m m a
+ -> Mat m m a
+lt_part matrix =
+ construct lambda
+ where
+ lambda i j = if i >= j then matrix !!! (i,j) else 0
+
+
+-- | Given a square @matrix@, return a new matrix of the same size
+-- containing only the below-diagonal entries of @matrix@. The on-
+-- and above-diagonal entries are set to zero.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
+-- >>> lt_part_strict m
+-- ((0,0,0),(4,0,0),(7,8,0))
+--
+lt_part_strict :: (Arity m, Ring.C a)
+ => Mat m m a
+ -> Mat m m a
+lt_part_strict matrix =
+ construct lambda
+ where
+ lambda i j = if i > j then matrix !!! (i,j) else 0
+
+
+-- | Given a square @matrix@, return a new matrix of the same size
+-- containing only the on-diagonal and above-diagonal entries of
+-- @matrix@. The below-diagonal entries are set to zero.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
+-- >>> ut_part m
+-- ((1,2,3),(0,5,6),(0,0,9))
+--
+ut_part :: (Arity m, Ring.C a)
+ => Mat m m a
+ -> Mat m m a
+ut_part = transpose . lt_part . transpose
+
+
+-- | Given a square @matrix@, return a new matrix of the same size
+-- containing only the above-diagonal entries of @matrix@. The on-
+-- and below-diagonal entries are set to zero.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
+-- >>> ut_part_strict m
+-- ((0,2,3),(0,0,6),(0,0,0))
+--
+ut_part_strict :: (Arity m, Ring.C a)
+ => Mat m m a
+ -> Mat m m a
+ut_part_strict = transpose . lt_part_strict . transpose
+
+
+-- | Compute the trace of a square matrix, the sum of the elements
+-- which lie on its diagonal. We require the matrix to be
+-- square to avoid ambiguity in the return type which would ideally
+-- have dimension min(m,n) supposing an m-by-n matrix.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
+-- >>> trace m
+-- 15
+--
+trace :: (Arity m, Ring.C a) => Mat m m a -> a
+trace = element_sum2 . diagonal
+
+
+
+-- | Zip together two matrices.
+--
+-- TODO: don't cheat with construct (map V.zips instead).
+--
+-- Examples:
+--
+-- >>> let m1 = fromList [[1],[1],[1]] :: Col3 Int
+-- >>> let m2 = fromList [[1],[2],[3]] :: Col3 Int
+-- >>> zip2 m1 m2
+-- (((1,1)),((1,2)),((1,3)))
+--
+-- >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int
+-- >>> let m2 = fromList [[1,1],[1,1]] :: Mat2 Int
+-- >>> zip2 m1 m2
+-- (((1,1),(2,1)),((3,1),(4,1)))
+--
+zip2 :: (Arity m, Arity n) => Mat m n a -> Mat m n b -> Mat m n (a,b)
+zip2 m1 m2 =
+ construct lambda
+ where
+ lambda i j = (m1 !!! (i,j), m2 !!! (i,j))
+
+
+-- | Zip together three matrices.
+--
+-- TODO: don't cheat with construct (map V.zips instead).
+--
+-- Examples:
+--
+-- >>> let m1 = fromList [[1],[1],[1]] :: Col3 Int
+-- >>> let m2 = fromList [[1],[2],[3]] :: Col3 Int
+-- >>> let m3 = fromList [[4],[5],[6]] :: Col3 Int
+-- >>> zip2three m1 m2 m3
+-- (((1,1,4)),((1,2,5)),((1,3,6)))
+--
+-- >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int
+-- >>> let m2 = fromList [[1,1],[1,1]] :: Mat2 Int
+-- >>> let m3 = fromList [[8,2],[6,3]] :: Mat2 Int
+-- >>> zip2three m1 m2 m3
+-- (((1,1,8),(2,1,2)),((3,1,6),(4,1,3)))
+--
+zip2three :: (Arity m, Arity n)
+ => Mat m n a
+ -> Mat m n a
+ -> Mat m n a
+ -> Mat m n (a,a,a)
+zip2three m1 m2 m3 =
+ construct lambda
+ where
+ lambda i j = (m1 !!! (i,j), m2 !!! (i,j), m3 !!! (i,j))
+
+
+-- | Zip together two matrices using the supplied function.
+--
+-- Examples:
+--
+-- >>> let c1 = fromList [[1],[2],[3]] :: Col3 Integer
+-- >>> let c2 = fromList [[4],[5],[6]] :: Col3 Integer
+-- >>> zipwith2 (^) c1 c2
+-- ((1),(32),(729))
+--
+zipwith2 :: (Arity m, Arity n)
+ => (a -> b -> c)
+ -> Mat m n a
+ -> Mat m n b
+ -> Mat m n c
+zipwith2 f c1 c2 =
+ construct lambda
+ where
+ lambda i j = f (c1 !!! (i,j)) (c2 !!! (i,j))
+
+
+-- | Map a function over a matrix of any dimensions.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
+-- >>> map2 (^2) m
+-- ((1,4),(9,16))
+--
+map2 :: (a -> b) -> Mat m n a -> Mat m n b
+map2 f (Mat rows) =
+ Mat $ V.map g rows
+ where
+ g = V.map f
+
+
+-- | Fold over the entire matrix passing the coordinates @i@ and @j@
+-- (of the row/column) to the accumulation function. The fold occurs
+-- from top-left to bottom-right.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
+-- >>> ifoldl2 (\i j cur _ -> cur + i + j) 0 m
+-- 18
+--
+ifoldl2 :: forall a b m n.
+ (Int -> Int -> b -> a -> b)
+ -> b
+ -> Mat m n a
+ -> b
+ifoldl2 f initial (Mat rows) =
+ V.ifoldl row_function initial rows
+ where
+ -- | The order that we need this in (so that @g idx@ makes sense)
+ -- is a little funny. So that we don't need to pass weird
+ -- functions into ifoldl2, we swap the second and third
+ -- arguments of @f@ calling the result @g@.
+ g :: Int -> b -> Int -> a -> b
+ g w x y = f w y x
+
+ row_function :: b -> Int -> Vec n a -> b
+ row_function rowinit idx r = V.ifoldl (g idx) rowinit r
+
+
+-- | Left fold over the entries of a matrix (top-left to bottom-right).
+--
+foldl2 :: forall a b m n.
+ (b -> a -> b)
+ -> b
+ -> Mat m n a
+ -> b
+foldl2 f initial matrix =
+ -- Use the index fold but ignore the index arguments.
+ let g _ _ = f in ifoldl2 g initial matrix
+
+
+-- | Fold over the entire matrix passing the coordinates @i@ and @j@
+-- (of the row/column) to the accumulation function. The fold occurs
+-- from bottom-right to top-left.
+--
+-- The order of the arguments in the supplied function are different
+-- from those in V.ifoldr; we keep them similar to ifoldl2.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
+-- >>> ifoldr2 (\i j cur _ -> cur + i + j) 0 m
+-- 18
+--
+ifoldr2 :: forall a b m n.
+ (Int -> Int -> b -> a -> b)
+ -> b
+ -> Mat m n a
+ -> b
+ifoldr2 f initial (Mat rows) =
+ V.ifoldr row_function initial rows
+ where
+ -- | Swap the order of arguments in @f@ so that it agrees with the
+ -- @f@ passed to ifoldl2.
+ g :: Int -> Int -> a -> b -> b
+ g w x y z = f w x z y
+
+ row_function :: Int -> Vec n a -> b -> b
+ row_function idx r rowinit = V.ifoldr (g idx) rowinit r
+
+
+-- | Map a function over a matrix of any dimensions, passing the
+-- coordinates @i@ and @j@ to the function @f@.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
+-- >>> imap2 (\i j _ -> i+j) m
+-- ((0,1),(1,2))
+--
+imap2 :: (Int -> Int -> a -> b) -> Mat m n a -> Mat m n b
+imap2 f (Mat rows) =
+ Mat $ V.imap g rows
+ where
+ g i = V.imap (f i)
+
+
+-- | Reverse the order of elements in a matrix.
+--
+-- Examples:
+--
+-- >>> let m1 = fromList [[1,2,3]] :: Row3 Int
+-- >>> reverse2 m1
+-- ((3,2,1))
+--
+-- >>> let m1 = vec3d (1,2,3 :: Int)
+-- >>> reverse2 m1
+-- ((3),(2),(1))
+--
+-- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
+-- >>> reverse2 m
+-- ((9,8,7),(6,5,4),(3,2,1))
+--
+reverse2 :: (Arity m, Arity n) => Mat m n a -> Mat m n a
+reverse2 (Mat rows) = Mat $ V.reverse $ V.map V.reverse rows
+
+
+-- | Unsafely set the (i,j) element of the given matrix.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
+-- >>> set_idx m (1,1) 17
+-- ((1,2,3),(4,17,6),(7,8,9))
+--
+set_idx :: forall m n a.
+ (Arity m, Arity n)
+ => Mat m n a
+ -> (Int, Int)
+ -> a
+ -> Mat m n a
+set_idx matrix (i,j) newval =
+ imap2 updater matrix
+ where
+ updater :: Int -> Int -> a -> a
+ updater k l existing =
+ if k == i && l == j
+ then newval
+ else existing
+
+
+-- | Compute the i,jth cofactor of the given @matrix@. This simply
+-- premultiplues the i,jth minor by (-1)^(i+j).
+cofactor :: (Arity m, Determined (Mat m m) a)
+ => Mat (S m) (S m) a
+ -> Int
+ -> Int
+ -> a
+cofactor matrix i j =
+ (-1)^(toInteger i + toInteger j) NP.* (minor matrix i j)
+
+
+-- | Compute the inverse of a matrix using cofactor expansion
+-- (generalized Cramer's rule).
+--
+-- Examples:
+--
+-- >>> let m1 = fromList [[37,22],[17,54]] :: Mat2 Double
+-- >>> let e1 = [54/1624, -22/1624] :: [Double]
+-- >>> let e2 = [-17/1624, 37/1624] :: [Double]
+-- >>> let expected = fromList [e1, e2] :: Mat2 Double
+-- >>> let actual = inverse m1
+-- >>> frobenius_norm (actual - expected) < 1e-12
+-- True
+--
+inverse :: (Arity m,
+ Determined (Mat (S m) (S m)) a,
+ Determined (Mat m m) a,
+ Field.C a)
+ => Mat (S m) (S m) a
+ -> Mat (S m) (S m) a
+inverse matrix =
+ (1 / (determinant matrix)) *> (transpose $ construct lambda)
+ where
+ lambda i j = cofactor matrix i j
+
+
+
+-- | Retrieve the rows of a matrix as a column matrix. If the given
+-- matrix is m-by-n, the result would be an m-by-1 column whose
+-- entries are 1-by-n row matrices.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
+-- >>> (rows2 m) !!! (0,0)
+-- ((1,2))
+-- >>> (rows2 m) !!! (1,0)
+-- ((3,4))
+--
+rows2 :: (Arity m, Arity n)
+ => Mat m n a
+ -> Col m (Row n a)
+rows2 (Mat rows) =
+ Mat $ V.map (mk1. Mat . mk1) rows
+
+
+
+-- | Sum the elements of a matrix.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,-1],[3,4]] :: Mat2 Int
+-- >>> element_sum2 m
+-- 7
+--
+element_sum2 :: (Arity m, Arity n, Additive.C a) => Mat m n a -> a
+element_sum2 = foldl2 (+) zero