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