Add tolerant versions of is_{upper,lower}_triangular.

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

diff --git a/src/Linear/Matrix.hs b/src/Linear/Matrix.hs
index 69a74b50eeacf7ea944c790d3858dd0cba12c266..c6f4a83c71af1c5231113ec555fa79f2dcfcc038 100644 (file)
--- a/src/Linear/Matrix.hs
+++ b/src/Linear/Matrix.hs
@@ -48,11 +48,13 @@ import Data.Vector.Fixed.Cont (Arity, arity)
 import Linear.Vector
 import Normed
 
 import Linear.Vector
 import Normed
 

-import NumericPrelude hiding ((*), abs)
-import qualified NumericPrelude as NP ((*))
+import NumericPrelude hiding ( (*), abs )
+import qualified NumericPrelude as NP ( (*) )
+import qualified Algebra.Absolute as Absolute ( C )
+import Algebra.Absolute ( abs )
+import qualified Algebra.Additive as Additive

 import qualified Algebra.Algebraic as Algebraic
 import Algebra.Algebraic (root)
 import qualified Algebra.Algebraic as Algebraic
 import Algebra.Algebraic (root)

-import qualified Algebra.Additive as Additive

 import qualified Algebra.Ring as Ring
 import qualified Algebra.Module as Module
 import qualified Algebra.RealRing as RealRing
 import qualified Algebra.Ring as Ring
 import qualified Algebra.Module as Module
 import qualified Algebra.RealRing as RealRing


@@ -250,21 +252,26 @@ cholesky m = construct r
 
 
 -- | Returns True if the given matrix is upper-triangular, and False
 
 
 -- | Returns True if the given matrix is upper-triangular, and False

---   otherwise.
+--   otherwise. The parameter @epsilon@ lets the caller choose a
+--   tolerance.

 --
 --   Examples:
 --
 --
 --   Examples:
 --

---   >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
+--   >>> let m = fromList [[1,1],[1e-12,1]] :: Mat2 Double

 --   >>> is_upper_triangular m
 --   False
 --   >>> is_upper_triangular m
 --   False

---
---   >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
---   >>> is_upper_triangular m
+--   >>> is_upper_triangular' 1e-10 m

 --   True
 --
 --   True
 --

-is_upper_triangular :: (Eq a, Ring.C a, Arity m, Arity n)
-                    => Mat m n a -> Bool
-is_upper_triangular m =
+--   TODO:
+--
+--     1. Don't cheat with lists.
+--
+is_upper_triangular' :: (Ord a, Ring.C a, Absolute.C a, Arity m, Arity n)
+                    => a -- ^ The tolerance @epsilon@.
+                    -> Mat m n a
+                    -> Bool
+is_upper_triangular' epsilon m =

   and $ concat results
   where
     results = [[ test i j | i <- [0..(nrows m)-1]] | j <- [0..(ncols m)-1] ]
   and $ concat results
   where
     results = [[ test i j | i <- [0..(nrows m)-1]] | j <- [0..(ncols m)-1] ]


@@ -272,11 +279,36 @@ is_upper_triangular m =
     test :: Int -> Int -> Bool
     test i j
       | i <= j = True
     test :: Int -> Int -> Bool
     test i j
       | i <= j = True

-      | otherwise = m !!! (i,j) == 0
+      -- use "less than or equal to" so zero is a valid epsilon
+      | otherwise = abs (m !!! (i,j)) <= epsilon
+
+
+-- | Returns True if the given matrix is upper-triangular, and False
+--   otherwise. A specialized version of 'is_upper_triangular\'' with
+--   @epsilon = 0@.
+--
+--   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
+--
+--   TODO:
+--
+--     1. The Ord constraint is too strong here, Eq would suffice.
+--
+is_upper_triangular :: (Ord a, Ring.C a, Absolute.C a, Arity m, Arity n)
+                    => Mat m n a -> Bool
+is_upper_triangular = is_upper_triangular' 0

 
 
 -- | Returns True if the given matrix is lower-triangular, and False
 
 
 -- | Returns True if the given matrix is lower-triangular, and False

---   otherwise.
+--   otherwise. This is a specialized version of 'is_lower_triangular\''
+--   with @epsilon = 0@.

 --
 --   Examples:
 --
 --
 --   Examples:
 --


@@ -288,8 +320,9 @@ is_upper_triangular m =
 --   >>> is_lower_triangular m
 --   False
 --
 --   >>> is_lower_triangular m
 --   False
 --

-is_lower_triangular :: (Eq a,
+is_lower_triangular :: (Ord a,

                         Ring.C a,
                         Ring.C a,

+                        Absolute.C a,

                         Arity m,
                         Arity n)
                     => Mat m n a
                         Arity m,
                         Arity n)
                     => Mat m n a


@@ -297,6 +330,29 @@ is_lower_triangular :: (Eq a,
 is_lower_triangular = is_upper_triangular . transpose
 
 
 is_lower_triangular = is_upper_triangular . transpose
 
 

+-- | Returns True if the given matrix is lower-triangular, and False
+--   otherwise. The parameter @epsilon@ lets the caller choose a
+--   tolerance.
+--
+--   Examples:
+--
+--   >>> let m = fromList [[1,1e-12],[1,1]] :: Mat2 Double
+--   >>> is_lower_triangular m
+--   False
+--   >>> is_lower_triangular' 1e-12 m
+--   True
+--
+is_lower_triangular' :: (Ord a,
+                         Ring.C a,
+                         Absolute.C a,
+                         Arity m,
+                         Arity n)
+                    => a -- ^ The tolerance @epsilon@.
+                    -> Mat m n a
+                    -> Bool
+is_lower_triangular' epsilon = (is_upper_triangular' epsilon) . transpose
+
+

 -- | Returns True if the given matrix is triangular, and False
 --   otherwise.
 --
 -- | Returns True if the given matrix is triangular, and False
 --   otherwise.
 --


@@ -314,8 +370,9 @@ is_lower_triangular = is_upper_triangular . transpose
 --   >>> is_triangular m
 --   False
 --
 --   >>> is_triangular m
 --   False
 --

-is_triangular :: (Eq a,
+is_triangular :: (Ord a,

                   Ring.C a,
                   Ring.C a,

+                  Absolute.C a,

                   Arity m,
                   Arity n)
               => Mat m n a
                   Arity m,
                   Arity n)
               => Mat m n a


@@ -353,8 +410,9 @@ class (Eq a, Ring.C a) => Determined p a where
 instance (Eq a, Ring.C a) => Determined (Mat (S Z) (S Z)) a where
   determinant (Mat rows) = (V.head . V.head) rows
 
 instance (Eq a, Ring.C a) => Determined (Mat (S Z) (S Z)) a where
   determinant (Mat rows) = (V.head . V.head) rows
 

-instance (Eq a,
+instance (Ord a,

           Ring.C a,
           Ring.C a,

+          Absolute.C a,

           Arity n,
           Determined (Mat (S n) (S n)) a)
          => Determined (Mat (S (S n)) (S (S n))) a where
           Arity n,
           Determined (Mat (S n) (S n)) a)
          => Determined (Mat (S (S n)) (S (S n))) a where