Add triangular functions, determinant, minor (all half-baked) to Linear.Matrix.

diff --git a/src/Linear/Matrix.hs b/src/Linear/Matrix.hs
index 4f2e845da6b97380b703441c1fc6963c568cb943..5bb5519664c5fa4028d284775493bb7da4f5598e 100644 (file)
--- a/src/Linear/Matrix.hs
+++ b/src/Linear/Matrix.hs
@@ -13,10 +13,10 @@ import Data.List (intercalate)
 
 import Data.Vector.Fixed (
   Dim,
+  N1,
   Vector
   )
 import qualified Data.Vector.Fixed as V (
-  N1,
   and,
   fromList,
   length,
@@ -239,6 +239,139 @@ cholesky m = construct r
           | otherwise = 0
 
 
+-- | Returns True if the given matrix is upper-triangular, and False
+--   otherwise.
+--
+--   Examples:
+--
+--   >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
+--   >>> is_upper_triangular m
+--   False
+--
+--   >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
+--   >>> is_upper_triangular m
+--   True
+--
+is_upper_triangular :: (Eq a, Ring.C a, Vector w a) => Mat v w a -> Bool
+is_upper_triangular m =
+  and $ concat results
+  where
+    results = [[ test i j | i <- [0..(nrows m)-1]] | j <- [0..(ncols m)-1] ]
+
+    test :: Int -> Int -> Bool
+    test i j
+      | i <= j = True
+      | otherwise = m !!! (i,j) == 0
+
+
+-- | Returns True if the given matrix is lower-triangular, and False
+--   otherwise.
+--
+--   Examples:
+--
+--   >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
+--   >>> is_lower_triangular m
+--   True
+--
+--   >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
+--   >>> is_lower_triangular m
+--   False
+--
+is_lower_triangular :: (Eq a,
+                        Ring.C a,
+                        Vector w a,
+                        Vector w (v a),
+                        Vector v a)
+                    => Mat v w a
+                    -> Bool
+is_lower_triangular = is_upper_triangular . transpose
+
+
+-- | Returns True if the given matrix is triangular, and False
+--   otherwise.
+--
+--   Examples:
+--
+--   >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
+--   >>> is_triangular m
+--   True
+--
+--   >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
+--   >>> is_triangular m
+--   True
+--
+--   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
+--   >>> is_triangular m
+--   False
+--
+is_triangular :: (Eq a,
+                  Ring.C a,
+                  Vector w a,
+                  Vector w (v a),
+                  Vector v a)
+               => Mat v w a
+              -> Bool
+is_triangular m = is_upper_triangular m || is_lower_triangular m
+
+
+-- | Return the (i,j)th minor of m.
+--
+--   Examples:
+--
+--   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
+--   >>> minor m 0 0 :: Mat2 Int
+--   ((5,6),(8,9))
+--   >>> minor m 1 1 :: Mat2 Int
+--   ((1,3),(7,9))
+--
+minor :: (Dim v ~ S (Dim u),
+          Dim w ~ S (Dim z),
+          Vector z a,
+          Vector u (w a),
+          Vector u (z a))
+      => Mat v w a
+      -> Int
+      -> Int
+      -> Mat u z a
+minor (Mat rows) i j = m
+  where
+    rows' = delete rows i
+    m = Mat $ V.map ((flip delete) j) rows'
+
+
+determinant :: (Eq a,
+                Ring.C a,
+                Vector w a,
+                Vector w (v a),
+                Vector v a,
+                Dim v ~ S r,
+                Dim w ~ S t)
+             => Mat v w a
+             -> a
+determinant m
+  | is_triangular m = product [ m !!! (i,i) | i <- [0..(nrows m)-1] ]
+  | otherwise = undefined --determinant_recursive m
+
+{-
+determinant_recursive :: forall v w a r c.
+                         (Eq a,
+                          Ring.C a,
+                          Vector w a)
+                       => Mat (v r) (w c) a
+                       -> a
+determinant_recursive m
+  | (ncols m) == 0 || (nrows m) == 0 = error "don't do that"
+  | (ncols m) == 1 && (nrows m) == 1 = m !!! (0,0) -- Base case
+  | otherwise =
+      sum [ (-1)^(1+(toInteger j)) NP.* (m' 1 j) NP.* (det_minor 1 j)
+            | j <- [0..(ncols m)-1] ]
+      where
+        m' i j = m !!! (i,j)
+
+        det_minor :: Int -> Int -> a
+        det_minor i j = determinant (minor m i j)
+-}
+
 -- | Matrix multiplication. Our 'Num' instance doesn't define one, and
 --   we need additional restrictions on the result type anyway.
 --
@@ -359,6 +492,9 @@ instance (Algebraic.C a,
 --   >>> fixed_point g eps u0
 --   ((1.0728549599342185),(1.0820591495686167))
 --
+vec1d :: (a) -> Mat D1 D1 a
+vec1d (x) = Mat (D1 (D1 x))
+
 vec2d :: (a,a) -> Mat D2 D1 a
 vec2d (x,y) = Mat (D2 (D1 x) (D1 y))
 
@@ -374,7 +510,7 @@ scalar :: a -> Mat D1 D1 a
 scalar x = Mat (D1 (D1 x))
 
 dot :: (RealRing.C a,
-        Dim w ~ V.N1,
+        Dim w ~ N1,
         Dim v ~ S n,
         Vector v a,
         Vector w a,
@@ -397,7 +533,7 @@ v1 `dot` v2 = ((transpose v1) * v2) !!! (0, 0)
 --
 angle :: (Transcendental.C a,
           RealRing.C a,
-          Dim w ~ V.N1,
+          Dim w ~ N1,
           Dim v ~ S n,
           Vector w (w a),
           Vector v [a],