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

[numerical-analysis.git] / src / Linear / Matrix.hs

diff --git a/src/Linear/Matrix.hs b/src/Linear/Matrix.hs
index ef0e9b6ff3ccba65c310b81126ea915d543e0411..5bb5519664c5fa4028d284775493bb7da4f5598e 100644 (file)
--- a/src/Linear/Matrix.hs
+++ b/src/Linear/Matrix.hs
@@ -13,23 +13,19 @@ import Data.List (intercalate)
 
 import Data.Vector.Fixed (
   Dim,
 
 import Data.Vector.Fixed (
   Dim,

+  N1,

   Vector
   )
 import qualified Data.Vector.Fixed as V (
   Vector
   )
 import qualified Data.Vector.Fixed as V (

-  Fun(..),
-  N1,

   and,
   and,

-  eq,
-  foldl,

   fromList,
   length,
   map,
   fromList,
   length,
   map,

-  maximum,

   replicate,
   toList,
   zipWith
   )
   replicate,
   toList,
   zipWith
   )

-import Data.Vector.Fixed.Internal (Arity, arity, S, Dim)
+import Data.Vector.Fixed.Internal (Arity, arity, S)

 import Linear.Vector
 import Normed
 
 import Linear.Vector
 import Normed
 


@@ -37,13 +33,9 @@ import NumericPrelude hiding ((*), abs)
 import qualified NumericPrelude as NP ((*))
 import qualified Algebra.Algebraic as Algebraic
 import Algebra.Algebraic (root)
 import qualified NumericPrelude as NP ((*))
 import qualified Algebra.Algebraic as Algebraic
 import Algebra.Algebraic (root)

-import qualified Algebra.Absolute as Absolute

 import qualified Algebra.Additive as Additive
 import qualified Algebra.Ring as Ring
 import qualified Algebra.Additive as Additive
 import qualified Algebra.Ring as Ring

-import Algebra.Absolute (abs)
-import qualified Algebra.Field as Field

 import qualified Algebra.Module as Module
 import qualified Algebra.Module as Module

-import qualified Algebra.RealField as RealField

 import qualified Algebra.RealRing as RealRing
 import qualified Algebra.ToRational as ToRational
 import qualified Algebra.Transcendental as Transcendental
 import qualified Algebra.RealRing as RealRing
 import qualified Algebra.ToRational as ToRational
 import qualified Algebra.Transcendental as Transcendental


@@ -247,6 +239,139 @@ cholesky m = construct r
           | otherwise = 0
 
 
           | 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.
 --
 -- | Matrix multiplication. Our 'Num' instance doesn't define one, and
 --   we need additional restrictions on the result type anyway.
 --


@@ -294,7 +419,7 @@ instance (Ring.C a,
          => Ring.C (Mat v w a) where
   -- The first * is ring multiplication, the second is matrix
   -- multiplication.
          => Ring.C (Mat v w a) where
   -- The first * is ring multiplication, the second is matrix
   -- multiplication.

-  m1 * m2 = m1 * m1
+  m1 * m2 = m1 * m2

 
 
 instance (Ring.C a,
 
 
 instance (Ring.C a,


@@ -359,14 +484,17 @@ instance (Algebraic.C a,
 --
 --   >>> import Roots.Simple
 --   >>> let h = 0.5 :: Double
 --
 --   >>> import Roots.Simple
 --   >>> let h = 0.5 :: Double

---   >>> let g1 (Mat (D2 (D1 x) (D1 y))) = 1.0 + h*exp(-(x^2))/(1.0 + y^2)
---   >>> let g2 (Mat (D2 (D1 x) (D1 y))) = 0.5 + h*atan(x^2 + y^2)
+--   >>> let g1 (Mat (D2 (D1 x) (D1 y))) = 1.0 + h NP.* exp(-(x^2))/(1.0 + y^2)
+--   >>> let g2 (Mat (D2 (D1 x) (D1 y))) = 0.5 + h NP.* atan(x^2 + y^2)

 --   >>> let g u = vec2d ((g1 u), (g2 u))
 --   >>> let u0 = vec2d (1.0, 1.0)
 --   >>> let eps = 1/(10^9)
 --   >>> fixed_point g eps u0
 --   >>> let g u = vec2d ((g1 u), (g2 u))
 --   >>> let u0 = vec2d (1.0, 1.0)
 --   >>> let eps = 1/(10^9)
 --   >>> fixed_point g eps u0

---   (1.0728549599342185,1.0820591495686167)
+--   ((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))
 
 vec2d :: (a,a) -> Mat D2 D1 a
 vec2d (x,y) = Mat (D2 (D1 x) (D1 y))
 


@@ -382,7 +510,7 @@ scalar :: a -> Mat D1 D1 a
 scalar x = Mat (D1 (D1 x))
 
 dot :: (RealRing.C 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,
         Dim v ~ S n,
         Vector v a,
         Vector w a,


@@ -405,7 +533,7 @@ v1 `dot` v2 = ((transpose v1) * v2) !!! (0, 0)
 --
 angle :: (Transcendental.C a,
           RealRing.C a,
 --
 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],
           Dim v ~ S n,
           Vector w (w a),
           Vector v [a],