+-- | 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