X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FLinear%2FMatrix.hs;h=a5fbb8403b0832becf04b0d48a68be4f0725345c;hb=9ef91ad4ec3a5c0966f0850d40310722b6c38b68;hp=34920b4f025e7e2027588ab07c14d6772b679ab2;hpb=4e464a486bef07db44de9c3d3fae0c8094401b09;p=numerical-analysis.git

diff --git a/src/Linear/Matrix.hs b/src/Linear/Matrix.hs
index 34920b4..a5fbb84 100644
--- a/src/Linear/Matrix.hs
+++ b/src/Linear/Matrix.hs
@@ -37,6 +37,7 @@ import qualified Data.Vector.Fixed as V (
   and,
   fromList,
   head,
+  ifoldl,
   length,
   map,
   maximum,
@@ -91,6 +92,11 @@ type Col3 a = Col N3 a
 type Col4 a = Col N4 a
 type Col5 a = Col N5 a
 
+-- We need a big column for Gaussian quadrature.
+type N10 = S (S (S (S (S N5))))
+type Col10 a = Col N10 a
+
+
 instance (Eq a) => Eq (Mat m n a) where
   -- | Compare a row at a time.
   --
@@ -146,60 +152,47 @@ fromList vs = Mat (V.fromList $ map V.fromList vs)
 
 -- | Unsafe indexing.
 (!!!) :: (Arity m, Arity n) => Mat m n a -> (Int, Int) -> a
-(!!!) m (i, j) = (row m i) ! j
+(!!!) (Mat rows) (i, j) = (rows ! i) ! j
+
 
 -- | Safe indexing.
-(!!?) :: Mat m n a -> (Int, Int) -> Maybe a
-(!!?) m@(Mat rows) (i, j)
+(!!?) :: (Arity m, Arity n) => Mat m n a -> (Int, Int) -> Maybe a
+(!!?) matrix (i, j)
   | i < 0 || j < 0 = Nothing
-  | i > V.length rows = Nothing
-  | otherwise = if j > V.length (row m j)
-                then Nothing
-                else Just $ (row m j) ! j
+  | i > (nrows matrix) - 1 = Nothing
+  | j > (ncols matrix) - 1 = Nothing
+  | otherwise = Just $ matrix !!! (i,j)
 
 
 -- | The number of rows in the matrix.
 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.
 ncols :: forall m n a. (Arity n) => Mat m n a -> Int
 ncols _ = arity (undefined :: n)
 
 
--- | Return the @i@th row of @m@. Unsafe.
-row :: Mat m n a -> Int -> (Vec n a)
-row (Mat rows) i = rows ! i
-
-
 -- | Return the @i@th row of @m@ as a matrix. Unsafe.
-row' :: (Arity m, Arity n) => Mat m n a -> Int -> Row n a
-row' m i =
+row :: (Arity m, Arity n) => Mat m n a -> Int -> Row n a
+row m i =
   construct lambda
   where
     lambda _ j = m !!! (i, j)
 
 
--- | Return the @j@th column of @m@. Unsafe.
-column :: Mat m n a -> Int -> (Vec m a)
-column (Mat rows) j =
-  V.map (element j) rows
-  where
-    element = flip (!)
-
-
 -- | Return the @j@th column of @m@ as a matrix. Unsafe.
-column' :: (Arity m, Arity n) => Mat m n a -> Int -> Col m a
-column' m j =
+column :: (Arity m, Arity n) => Mat m n a -> Int -> Col m a
+column m j =
   construct lambda
   where
     lambda i _ = m !!! (i, j)
 
 
 -- | 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.
+--   implementation, but I don't see a better way.
 --
 --   TODO: Don't cheat with fromList.
 --
@@ -210,9 +203,10 @@ column' m j =
 --   ((1,3),(2,4))
 --
 transpose :: (Arity m, Arity n) => Mat m n a -> Mat n m a
-transpose m = Mat $ V.fromList column_list
+transpose matrix =
+  construct lambda
   where
-    column_list = [ column m i | i <- [0..(ncols m)-1] ]
+    lambda i j = matrix !!! (j,i)
 
 
 -- | Is @m@ symmetric?
@@ -264,6 +258,7 @@ 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
 --   values on the diagonal of @r@ positive.
@@ -767,3 +762,84 @@ trace matrix =
   let (Mat rows) = diagonal matrix
   in
     element_sum $ V.map V.head rows
+
+
+-- | Zip together two column matrices.
+--
+--   Examples:
+--
+--   >>> let m1 = fromList [[1],[1],[1]] :: Col3 Int
+--   >>> let m2 = fromList [[1],[2],[3]] :: Col3 Int
+--   >>> colzip m1 m2
+--   (((1,1)),((1,2)),((1,3)))
+--
+colzip :: Arity m => Col m a -> Col m a -> Col m (a,a)
+colzip c1 c2 =
+  construct lambda
+  where
+    lambda i j = (c1 !!! (i,j), c2 !!! (i,j))
+
+
+-- | Zip together two column matrices using the supplied function.
+--
+--   Examples:
+--
+--   >>> let c1 = fromList [[1],[2],[3]] :: Col3 Integer
+--   >>> let c2 = fromList [[4],[5],[6]] :: Col3 Integer
+--   >>> colzipwith (^) c1 c2
+--   ((1),(32),(729))
+--
+colzipwith :: Arity m
+           => (a -> a -> b)
+           -> Col m a
+           -> Col m a
+           -> Col m b
+colzipwith f c1 c2 =
+  construct lambda
+  where
+    lambda i j = f (c1 !!! (i,j)) (c2 !!! (i,j))
+
+
+-- | Map a function over a matrix of any dimensions.
+--
+--   Examples:
+--
+--   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
+--   >>> map2 (^2) m
+--   ((1,4),(9,16))
+--
+map2 :: (a -> b) -> Mat m n a -> Mat m n b
+map2 f (Mat rows) =
+  Mat $ V.map g rows
+  where
+    g = V.map f
+
+
+-- | Fold over the entire matrix passing the coordinates @i@ and @j@
+--   (of the row/column) to the accumulation function.
+--
+--   Examples:
+--
+--   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
+--   >>>  ifoldl2 (\i j cur _ -> cur + i + j) 0 m
+--   18
+--
+ifoldl2 :: forall a b m n.
+           (Int -> Int -> b -> a -> b)
+        -> b
+        -> Mat m n a
+        -> b
+ifoldl2 f initial (Mat rows) =
+  V.ifoldl row_function initial rows
+  where
+    -- | The order that we need this in (so that @g idx@ makes sense)
+    --   is a little funny. So that we don't need to pass weird
+    --   functions into ifoldl2, we swap the second and third
+    --   arguments of @f@ calling the result @g@.
+    g :: Int -> b -> Int -> a -> b
+    g w x y = f w y x
+
+    row_function :: b -> Int -> Vec n a -> b
+    row_function rowinit idx r = V.ifoldl (g idx) rowinit r
+
+