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 5bb5519664c5fa4028d284775493bb7da4f5598e..c6f4a83c71af1c5231113ec555fa79f2dcfcc038 100644 (file)
--- a/src/Linear/Matrix.hs
+++ b/src/Linear/Matrix.hs
@@ -2,38 +2,59 @@
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE MultiParamTypeClasses #-}

+{-# LANGUAGE NoMonomorphismRestriction #-}

 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE RebindableSyntax #-}
 
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE RebindableSyntax #-}
 

+-- | Boxed matrices; that is, boxed m-vectors of boxed n-vectors. We
+--   assume that the underlying representation is
+--   Data.Vector.Fixed.Boxed.Vec for simplicity. It was tried in
+--   generality and failed.
+--

 module Linear.Matrix
 where
 
 import Data.List (intercalate)
 
 import Data.Vector.Fixed (
 module Linear.Matrix
 where
 
 import Data.List (intercalate)
 
 import Data.Vector.Fixed (

-  Dim,
+  (!),

   N1,
   N1,

-  Vector
+  N2,
+  N3,
+  N4,
+  N5,
+  S,
+  Z,
+  generate,
+  mk1,
+  mk2,
+  mk3,
+  mk4,
+  mk5

   )
 import qualified Data.Vector.Fixed as V (
   and,
   fromList,
   )
 import qualified Data.Vector.Fixed as V (
   and,
   fromList,

+  head,

   length,
   map,
   length,
   map,

+  maximum,

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

-import Data.Vector.Fixed.Internal (Arity, arity, S)
+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


@@ -41,16 +62,14 @@ import qualified Algebra.ToRational as ToRational
 import qualified Algebra.Transcendental as Transcendental
 import qualified Prelude as P
 
 import qualified Algebra.Transcendental as Transcendental
 import qualified Prelude as P
 

-data Mat v w a = (Vector v (w a), Vector w a) => Mat (v (w a))
-type Mat1 a = Mat D1 D1 a
-type Mat2 a = Mat D2 D2 a
-type Mat3 a = Mat D3 D3 a
-type Mat4 a = Mat D4 D4 a
+data Mat m n a = (Arity m, Arity n) => Mat (Vec m (Vec n a))
+type Mat1 a = Mat N1 N1 a
+type Mat2 a = Mat N2 N2 a
+type Mat3 a = Mat N3 N3 a
+type Mat4 a = Mat N4 N4 a
+type Mat5 a = Mat N5 N5 a

 
 

--- We can't just declare that all instances of Vector are instances of
--- Eq unfortunately. We wind up with an overlapping instance for
--- w (w a).
-instance (Eq a, Vector v Bool, Vector w Bool) => Eq (Mat v w a) where
+instance (Eq a) => Eq (Mat m n a) where

   -- | Compare a row at a time.
   --
   --   Examples:
   -- | Compare a row at a time.
   --
   --   Examples:


@@ -70,7 +89,7 @@ instance (Eq a, Vector v Bool, Vector w Bool) => Eq (Mat v w a) where
       comp row1 row2 = V.and (V.zipWith (==) row1 row2)
 
 
       comp row1 row2 = V.and (V.zipWith (==) row1 row2)
 
 

-instance (Show a, Vector v String, Vector w String) => Show (Mat v w a) where
+instance (Show a) => Show (Mat m n a) where

   -- | Display matrices and vectors as ordinary tuples. This is poor
   --   practice, but these results are primarily displayed
   --   interactively and convenience trumps correctness (said the guy
   -- | Display matrices and vectors as ordinary tuples. This is poor
   --   practice, but these results are primarily displayed
   --   interactively and convenience trumps correctness (said the guy


@@ -94,22 +113,21 @@ instance (Show a, Vector v String, Vector w String) => Show (Mat v w a) where
           element_strings = P.map show v1l
 
 
           element_strings = P.map show v1l
 
 

-

 -- | Convert a matrix to a nested list.
 -- | Convert a matrix to a nested list.

-toList :: Mat v w a -> [[a]]
+toList :: Mat m n a -> [[a]]

 toList (Mat rows) = map V.toList (V.toList rows)
 
 -- | Create a matrix from a nested list.
 toList (Mat rows) = map V.toList (V.toList rows)
 
 -- | Create a matrix from a nested list.

-fromList :: (Vector v (w a), Vector w a, Vector v a) => [[a]] -> Mat v w a
+fromList :: (Arity m, Arity n) => [[a]] -> Mat m n a

 fromList vs = Mat (V.fromList $ map V.fromList vs)
 
 
 -- | Unsafe indexing.
 fromList vs = Mat (V.fromList $ map V.fromList vs)
 
 
 -- | Unsafe indexing.

-(!!!) :: (Vector w a) => Mat v w a -> (Int, Int) -> a
+(!!!) :: (Arity m, Arity n) => Mat m n a -> (Int, Int) -> a

 (!!!) m (i, j) = (row m i) ! j
 
 -- | Safe indexing.
 (!!!) m (i, j) = (row m i) ! j
 
 -- | Safe indexing.

-(!!?) :: Mat v w a -> (Int, Int) -> Maybe a
+(!!?) :: Mat m n a -> (Int, Int) -> Maybe a

 (!!?) m@(Mat rows) (i, j)
   | i < 0 || j < 0 = Nothing
   | i > V.length rows = Nothing
 (!!?) m@(Mat rows) (i, j)
   | i < 0 || j < 0 = Nothing
   | i > V.length rows = Nothing


@@ -119,27 +137,30 @@ fromList vs = Mat (V.fromList $ map V.fromList vs)
 
 
 -- | The number of rows in the matrix.
 
 
 -- | The number of rows in the matrix.

-nrows :: Mat v w a -> Int
-nrows (Mat rows) = V.length rows
+nrows :: forall m n a. (Arity m) => Mat m n a -> Int
+nrows _ = arity (undefined :: m)

 
 -- | The number of columns in the first row of the
 --   matrix. Implementation stolen from Data.Vector.Fixed.length.
 
 -- | The number of columns in the first row of the
 --   matrix. Implementation stolen from Data.Vector.Fixed.length.

-ncols :: forall v w a. (Vector w a) => Mat v w a -> Int
-ncols _ = (arity (undefined :: Dim w))
+ncols :: forall m n a. (Arity n) => Mat m n a -> Int
+ncols _ = arity (undefined :: n)
+

 
 -- | Return the @i@th row of @m@. Unsafe.
 
 -- | Return the @i@th row of @m@. Unsafe.

-row :: Mat v w a -> Int -> w a
+row :: Mat m n a -> Int -> (Vec n a)

 row (Mat rows) i = rows ! i
 
 
 -- | Return the @j@th column of @m@. Unsafe.
 row (Mat rows) i = rows ! i
 
 
 -- | Return the @j@th column of @m@. Unsafe.

-column :: (Vector v a) => Mat v w a -> Int -> v a
+column :: Mat m n a -> Int -> (Vec m a)

 column (Mat rows) j =
   V.map (element j) rows
   where
     element = flip (!)
 
 
 column (Mat rows) j =
   V.map (element j) rows
   where
     element = flip (!)
 
 

+
+

 -- | Transpose @m@; switch it's columns and its rows. This is a dirty
 --   implementation.. it would be a little cleaner to use imap, but it
 --   doesn't seem to work.
 -- | Transpose @m@; switch it's columns and its rows. This is a dirty
 --   implementation.. it would be a little cleaner to use imap, but it
 --   doesn't seem to work.


@@ -152,11 +173,7 @@ column (Mat rows) j =
 --   >>> transpose m
 --   ((1,3),(2,4))
 --
 --   >>> transpose m
 --   ((1,3),(2,4))
 --

-transpose :: (Vector w (v a),
-              Vector v a,
-              Vector w a)
-             => Mat v w a
-             -> Mat w v a
+transpose :: (Arity m, Arity n) => Mat m n a -> Mat n m a

 transpose m = Mat $ V.fromList column_list
   where
     column_list = [ column m i | i <- [0..(ncols m)-1] ]
 transpose m = Mat $ V.fromList column_list
   where
     column_list = [ column m i | i <- [0..(ncols m)-1] ]


@@ -174,13 +191,7 @@ transpose m = Mat $ V.fromList column_list
 --   >>> symmetric m2
 --   False
 --
 --   >>> symmetric m2
 --   False
 --

-symmetric :: (Vector v (w a),
-              Vector w a,
-              v ~ w,
-              Vector w Bool,
-              Eq a)
-             => Mat v w a
-             -> Bool
+symmetric :: (Eq a, Arity m) => Mat m m a -> Bool

 symmetric m =
   m == (transpose m)
 
 symmetric m =
   m == (transpose m)
 


@@ -190,26 +201,32 @@ symmetric m =
 --   entries in the matrix. The i,j entry of the resulting matrix will
 --   have the value returned by lambda i j.
 --
 --   entries in the matrix. The i,j entry of the resulting matrix will
 --   have the value returned by lambda i j.
 --

---   TODO: Don't cheat with fromList.
---

 --   Examples:
 --
 --   >>> let lambda i j = i + j
 --   >>> construct lambda :: Mat3 Int
 --   ((0,1,2),(1,2,3),(2,3,4))
 --
 --   Examples:
 --
 --   >>> let lambda i j = i + j
 --   >>> construct lambda :: Mat3 Int
 --   ((0,1,2),(1,2,3),(2,3,4))
 --

-construct :: forall v w a.
-             (Vector v (w a),
-              Vector w a)
-             => (Int -> Int -> a)
-             -> Mat v w a
-construct lambda = Mat rows
+construct :: forall m n a. (Arity m, Arity n)
+          => (Int -> Int -> a) -> Mat m n a
+construct lambda = Mat $ generate make_row

   where
   where

-    -- The arity trick is used in Data.Vector.Fixed.length.
-    imax = (arity (undefined :: Dim v)) - 1
-    jmax = (arity (undefined :: Dim w)) - 1
-    row' i = V.fromList [ lambda i j | j <- [0..jmax] ]
-    rows = V.fromList [ row' i | i <- [0..imax] ]
+    make_row :: Int -> Vec n a
+    make_row i = generate (lambda i)
+
+
+-- | Create an identity matrix with the right dimensions.
+--
+--   Examples:
+--
+--   >>> identity_matrix :: Mat3 Int
+--   ((1,0,0),(0,1,0),(0,0,1))
+--   >>> identity_matrix :: Mat3 Double
+--   ((1.0,0.0,0.0),(0.0,1.0,0.0),(0.0,0.0,1.0))
+--
+identity_matrix :: (Arity m, Ring.C a) => Mat m m a
+identity_matrix =
+  construct (\i j -> if i == j then (fromInteger 1) else (fromInteger 0))

 
 -- | Given a positive-definite matrix @m@, computes the
 --   upper-triangular matrix @r@ with (transpose r)*r == m and all
 
 -- | Given a positive-definite matrix @m@, computes the
 --   upper-triangular matrix @r@ with (transpose r)*r == m and all


@@ -223,13 +240,8 @@ construct lambda = Mat rows
 --   >>> (transpose (cholesky m1)) * (cholesky m1)
 --   ((20.000000000000004,-1.0),(-1.0,20.0))
 --
 --   >>> (transpose (cholesky m1)) * (cholesky m1)
 --   ((20.000000000000004,-1.0),(-1.0,20.0))
 --

-cholesky :: forall a v w.
-            (Algebraic.C a,
-             Vector v (w a),
-             Vector w a,
-             Vector v a)
-            => (Mat v w a)
-            -> (Mat v w a)
+cholesky :: forall m n a. (Algebraic.C a, Arity m, Arity n)
+         => (Mat m n a) -> (Mat m n a)

 cholesky m = construct r
   where
     r :: Int -> Int -> a
 cholesky m = construct r
   where
     r :: Int -> Int -> a


@@ -240,20 +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, Vector w a) => Mat v w 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] ]


@@ -261,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:
 --


@@ -277,16 +320,39 @@ 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,

-                        Vector w a,
-                        Vector w (v a),
-                        Vector v a)
-                    => Mat v w a
+                        Absolute.C a,
+                        Arity m,
+                        Arity n)
+                    => Mat m n a

                     -> Bool
 is_lower_triangular = is_upper_triangular . transpose
 
 
                     -> Bool
 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.
 --


@@ -304,12 +370,12 @@ 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,

-                  Vector w a,
-                  Vector w (v a),
-                  Vector v a)
-               => Mat v w a
+                  Absolute.C a,
+                  Arity m,
+                  Arity n)
+              => Mat m n a

               -> Bool
 is_triangular m = is_upper_triangular m || is_lower_triangular m
 
               -> Bool
 is_triangular m = is_upper_triangular m || is_lower_triangular m
 


@@ -324,73 +390,68 @@ is_triangular m = is_upper_triangular m || is_lower_triangular m
 --   >>> minor m 1 1 :: Mat2 Int
 --   ((1,3),(7,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
+minor :: (m ~ S r,
+          n ~ S t,
+          Arity r,
+          Arity t)
+      => Mat m n a

       -> Int
       -> Int
       -> Int
       -> Int

-      -> Mat u z a
+      -> Mat r t a

 minor (Mat rows) i j = m
   where
     rows' = delete rows i
     m = Mat $ V.map ((flip delete) j) rows'
 
 
 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.
+class (Eq a, Ring.C a) => Determined p a where
+  determinant :: (p a) -> a
+
+instance (Eq a, Ring.C a) => Determined (Mat (S Z) (S Z)) a where
+  determinant (Mat rows) = (V.head . V.head) rows
+
+instance (Ord 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
+  -- | The recursive definition with a special-case for triangular matrices.
+  --
+  --   Examples:
+  --
+  --   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
+  --   >>> determinant m
+  --   -1
+  --
+  determinant m
+    | is_triangular m = product [ m !!! (i,i) | i <- [0..(nrows m)-1] ]
+    | otherwise = determinant_recursive
+        where
+          m' i j = m !!! (i,j)
+
+          det_minor i j = determinant (minor m i j)
+
+          determinant_recursive =
+            sum [ (-1)^(toInteger j) NP.* (m' 0 j) NP.* (det_minor 0 j)
+              | j <- [0..(ncols m)-1] ]
+
+
+
+-- | Matrix multiplication.

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

---   >>> let m1 = fromList [[1,2,3], [4,5,6]]  :: Mat D2 D3 Int
---   >>> let m2 = fromList [[1,2],[3,4],[5,6]] :: Mat D3 D2 Int
+--   >>> let m1 = fromList [[1,2,3], [4,5,6]]  :: Mat N2 N3 Int
+--   >>> let m2 = fromList [[1,2],[3,4],[5,6]] :: Mat N3 N2 Int

 --   >>> m1 * m2
 --   ((22,28),(49,64))
 --
 infixl 7 *
 --   >>> m1 * m2
 --   ((22,28),(49,64))
 --
 infixl 7 *

-(*) :: (Ring.C a,
-         Vector v a,
-         Vector w a,
-         Vector z a,
-         Vector v (z a))
-        => Mat v w a
-        -> Mat w z a
-        -> Mat v z a
+(*) :: (Ring.C a, Arity m, Arity n, Arity p)
+        => Mat m n a
+        -> Mat n p a
+        -> Mat m p a

 (*) m1 m2 = construct lambda
   where
     lambda i j =
 (*) m1 m2 = construct lambda
   where
     lambda i j =


@@ -398,10 +459,7 @@ infixl 7 *
 
 
 
 
 
 

-instance (Ring.C a,
-          Vector v (w a),
-          Vector w a)
-         => Additive.C (Mat v w a) where
+instance (Ring.C a, Arity m, Arity n) => Additive.C (Mat m n a) where

 
   (Mat rows1) + (Mat rows2) =
     Mat $ V.zipWith (V.zipWith (+)) rows1 rows2
 
   (Mat rows1) + (Mat rows2) =
     Mat $ V.zipWith (V.zipWith (+)) rows1 rows2


@@ -412,20 +470,13 @@ instance (Ring.C a,
   zero = Mat (V.replicate $ V.replicate (fromInteger 0))
 
 
   zero = Mat (V.replicate $ V.replicate (fromInteger 0))
 
 

-instance (Ring.C a,
-          Vector v (w a),
-          Vector w a,
-          v ~ w)
-         => Ring.C (Mat v w a) where
+instance (Ring.C a, Arity m, Arity n, m ~ n) => Ring.C (Mat m n a) where

   -- The first * is ring multiplication, the second is matrix
   -- multiplication.
   m1 * m2 = m1 * m2
 
 
   -- The first * is ring multiplication, the second is matrix
   -- multiplication.
   m1 * m2 = m1 * m2
 
 

-instance (Ring.C a,
-          Vector v (w a),
-          Vector w a)
-         => Module.C a (Mat v w a) where
+instance (Ring.C a, Arity m, Arity n) => Module.C a (Mat m n a) where

   -- We can multiply a matrix by a scalar of the same type as its
   -- elements.
   x *> (Mat rows) = Mat $ V.map (V.map (NP.* x)) rows
   -- We can multiply a matrix by a scalar of the same type as its
   -- elements.
   x *> (Mat rows) = Mat $ V.map (V.map (NP.* x)) rows


@@ -433,13 +484,10 @@ instance (Ring.C a,
 
 instance (Algebraic.C a,
           ToRational.C a,
 
 instance (Algebraic.C a,
           ToRational.C a,

-          Vector v (w a),
-          Vector w a,
-          Vector v a,
-          Vector v [a])
-         => Normed (Mat v w a) where
-  -- | Generic p-norms. The usual norm in R^n is (norm_p 2). We treat
-  --   all matrices as big vectors.
+          Arity m)
+         => Normed (Mat (S m) N1 a) where
+  -- | Generic p-norms for vectors in R^n that are represented as nx1
+  --   matrices.

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


@@ -450,14 +498,12 @@ instance (Algebraic.C a,
   --   5.0
   --
   norm_p p (Mat rows) =
   --   5.0
   --
   norm_p p (Mat rows) =

-    (root p') $ sum [(fromRational' $ toRational x)^p' | x <- xs]
+    (root p') $ sum [fromRational' (toRational x)^p' | x <- xs]

     where
       p' = toInteger p
       xs = concat $ V.toList $ V.map V.toList rows
 
     where
       p' = toInteger p
       xs = concat $ V.toList $ V.map V.toList rows
 

-  -- | The infinity norm. We don't use V.maximum here because it
-  --   relies on a type constraint that the vector be non-empty and I
-  --   don't know how to pattern match it away.
+  -- | The infinity norm.

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


@@ -465,14 +511,30 @@ instance (Algebraic.C a,
   --   >>> norm_infty v1
   --   5
   --
   --   >>> norm_infty v1
   --   5
   --

-  norm_infty m@(Mat rows)
-    | nrows m == 0 || ncols m == 0 = 0
-    | otherwise =
-        fromRational' $ toRational $
-        P.maximum $ V.toList $ V.map (P.maximum . V.toList) rows
-
+  norm_infty (Mat rows) =
+    fromRational' $ toRational $ V.maximum $ V.map V.maximum rows

 
 
 
 

+-- | Compute the Frobenius norm of a matrix. This essentially treats
+--   the matrix as one long vector containing all of its entries (in
+--   any order, it doesn't matter).
+--
+--   Examples:
+--
+--   >>> let m = fromList [[1, 2, 3],[4,5,6],[7,8,9]] :: Mat3 Double
+--   >>> frobenius_norm m == sqrt 285
+--   True
+--
+--   >>> let m = fromList [[1, -1, 1],[-1,1,-1],[1,-1,1]] :: Mat3 Double
+--   >>> frobenius_norm m == 3
+--   True
+--
+frobenius_norm :: (Algebraic.C a, Ring.C a) => Mat m n a -> a
+frobenius_norm (Mat rows) =
+  sqrt $ element_sum $ V.map row_sum rows
+  where
+    -- | Square and add up the entries of a row.
+    row_sum = element_sum . V.map (^2)

 
 
 -- Vector helpers. We want it to be easy to create low-dimension
 
 
 -- Vector helpers. We want it to be easy to create low-dimension


@@ -483,41 +545,40 @@ instance (Algebraic.C a,
 --   Examples:
 --
 --   >>> import Roots.Simple
 --   Examples:
 --
 --   >>> import Roots.Simple

+--   >>> let fst m = m !!! (0,0)
+--   >>> let snd m = m !!! (1,0)

 --   >>> let h = 0.5 :: Double
 --   >>> let h = 0.5 :: Double

---   >>> 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 g1 m = 1.0 + h NP.* exp(-((fst m)^2))/(1.0 + (snd m)^2)
+--   >>> let g2 m = 0.5 + h NP.* atan((fst m)^2 + (snd m)^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
 --   ((1.0728549599342185),(1.0820591495686167))
 --
 --   >>> 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))
 --

-vec1d :: (a) -> Mat D1 D1 a
-vec1d (x) = Mat (D1 (D1 x))
+vec1d :: (a) -> Mat N1 N1 a
+vec1d (x) = Mat (mk1 (mk1 x))
+
+vec2d :: (a,a) -> Mat N2 N1 a
+vec2d (x,y) = Mat (mk2 (mk1 x) (mk1 y))

 
 

-vec2d :: (a,a) -> Mat D2 D1 a
-vec2d (x,y) = Mat (D2 (D1 x) (D1 y))
+vec3d :: (a,a,a) -> Mat N3 N1 a
+vec3d (x,y,z) = Mat (mk3 (mk1 x) (mk1 y) (mk1 z))

 
 

-vec3d :: (a,a,a) -> Mat D3 D1 a
-vec3d (x,y,z) = Mat (D3 (D1 x) (D1 y) (D1 z))
+vec4d :: (a,a,a,a) -> Mat N4 N1 a
+vec4d (w,x,y,z) = Mat (mk4 (mk1 w) (mk1 x) (mk1 y) (mk1 z))

 
 

-vec4d :: (a,a,a,a) -> Mat D4 D1 a
-vec4d (w,x,y,z) = Mat (D4 (D1 w) (D1 x) (D1 y) (D1 z))
+vec5d :: (a,a,a,a,a) -> Mat N5 N1 a
+vec5d (v,w,x,y,z) = Mat (mk5 (mk1 v) (mk1 w) (mk1 x) (mk1 y) (mk1 z))

 
 -- Since we commandeered multiplication, we need to create 1x1
 -- matrices in order to multiply things.
 
 -- Since we commandeered multiplication, we need to create 1x1
 -- matrices in order to multiply things.

-scalar :: a -> Mat D1 D1 a
-scalar x = Mat (D1 (D1 x))
-
-dot :: (RealRing.C a,
-        Dim w ~ N1,
-        Dim v ~ S n,
-        Vector v a,
-        Vector w a,
-        Vector w (v a),
-        Vector w (w a))
-       => Mat v w a
-       -> Mat v w a
+scalar :: a -> Mat N1 N1 a
+scalar x = Mat (mk1 (mk1 x))
+
+dot :: (RealRing.C a, n ~ N1, m ~ S t, Arity t)
+       => Mat m n a
+       -> Mat m n a

        -> a
 v1 `dot` v2 = ((transpose v1) * v2) !!! (0, 0)
 
        -> a
 v1 `dot` v2 = ((transpose v1) * v2) !!! (0, 0)
 


@@ -533,20 +594,105 @@ v1 `dot` v2 = ((transpose v1) * v2) !!! (0, 0)
 --
 angle :: (Transcendental.C a,
           RealRing.C a,
 --
 angle :: (Transcendental.C a,
           RealRing.C a,

-          Dim w ~ N1,
-          Dim v ~ S n,
-          Vector w (w a),
-          Vector v [a],
-          Vector v a,
-          Vector w a,
-          Vector v (w a),
-          Vector w (v a),
+          n ~ N1,
+          m ~ S t,
+          Arity t,

           ToRational.C a)
           ToRational.C a)

-          => Mat v w a
-          -> Mat v w a
+          => Mat m n a
+          -> Mat m n a

           -> a
 angle v1 v2 =
   acos theta
   where
    theta = (recip norms) NP.* (v1 `dot` v2)
    norms = (norm v1) NP.* (norm v2)
           -> a
 angle v1 v2 =
   acos theta
   where
    theta = (recip norms) NP.* (v1 `dot` v2)
    norms = (norm v1) NP.* (norm v2)

+
+
+
+-- | Given a square @matrix@, return a new matrix of the same size
+--   containing only the on-diagonal entries of @matrix@. The
+--   off-diagonal entries are set to zero.
+--
+--   Examples:
+--
+--   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
+--   >>> diagonal_part m
+--   ((1,0,0),(0,5,0),(0,0,9))
+--
+diagonal_part :: (Arity m, Ring.C a)
+         => Mat m m a
+         -> Mat m m a
+diagonal_part matrix =
+  construct lambda
+  where
+    lambda i j = if i == j then matrix !!! (i,j) else 0
+
+
+-- | Given a square @matrix@, return a new matrix of the same size
+--   containing only the on-diagonal and below-diagonal entries of
+--   @matrix@. The above-diagonal entries are set to zero.
+--
+--   Examples:
+--
+--   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
+--   >>> lt_part m
+--   ((1,0,0),(4,5,0),(7,8,9))
+--
+lt_part :: (Arity m, Ring.C a)
+        => Mat m m a
+        -> Mat m m a
+lt_part matrix =
+  construct lambda
+  where
+    lambda i j = if i >= j then matrix !!! (i,j) else 0
+
+
+-- | Given a square @matrix@, return a new matrix of the same size
+--   containing only the below-diagonal entries of @matrix@. The on-
+--   and above-diagonal entries are set to zero.
+--
+--   Examples:
+--
+--   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
+--   >>> lt_part_strict m
+--   ((0,0,0),(4,0,0),(7,8,0))
+--
+lt_part_strict :: (Arity m, Ring.C a)
+        => Mat m m a
+        -> Mat m m a
+lt_part_strict matrix =
+  construct lambda
+  where
+    lambda i j = if i > j then matrix !!! (i,j) else 0
+
+
+-- | Given a square @matrix@, return a new matrix of the same size
+--   containing only the on-diagonal and above-diagonal entries of
+--   @matrix@. The below-diagonal entries are set to zero.
+--
+--   Examples:
+--
+--   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
+--   >>> ut_part m
+--   ((1,2,3),(0,5,6),(0,0,9))
+--
+ut_part :: (Arity m, Ring.C a)
+        => Mat m m a
+        -> Mat m m a
+ut_part = transpose . lt_part . transpose
+
+
+-- | Given a square @matrix@, return a new matrix of the same size
+--   containing only the above-diagonal entries of @matrix@. The on-
+--   and below-diagonal entries are set to zero.
+--
+--   Examples:
+--
+--   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
+--   >>> ut_part_strict m
+--   ((0,2,3),(0,0,6),(0,0,0))
+--
+ut_part_strict :: (Arity m, Ring.C a)
+        => Mat m m a
+        -> Mat m m a
+ut_part_strict = transpose . lt_part_strict . transpose