Make safe indexing slower but not dependent on the number of rows/columns.

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

diff --git a/src/Linear/Matrix.hs b/src/Linear/Matrix.hs
index 054bb6e3afac2777ef16e200ff64b9c60a9fd3fd..611545fc9dc44e97bdf3c30f0333c7a7b7064ce6 100644 (file)
--- a/src/Linear/Matrix.hs
+++ b/src/Linear/Matrix.hs
@@ -38,7 +38,6 @@ import qualified Data.Vector.Fixed as V (
   fromList,
   head,
   ifoldl,
-  length,
   map,
   maximum,
   replicate,
@@ -145,67 +144,85 @@ instance (Show a) => Show (Mat m n a) where
 toList :: Mat m n a -> [[a]]
 toList (Mat rows) = map V.toList (V.toList rows)
 
+
 -- | Create a matrix from a nested list.
 fromList :: (Arity m, Arity n) => [[a]] -> Mat m n a
 fromList vs = Mat (V.fromList $ map V.fromList vs)
 
 
--- | Unsafe indexing.
+-- | Unsafe indexing. Much faster than the safe 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)
-  | 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
+--
+--   Examples:
+--
+--   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
+--   >>> m !!? (-1,-1)
+--   Nothing
+--   >>>  m !!? (-1,0)
+--   Nothing
+--   >>>  m !!? (-1,1)
+--   Nothing
+--   >>>  m !!? (0,-1)
+--   Nothing
+--   >>>  m !!? (0,0)
+--   Just 1
+--   >>>  m !!? (0,1)
+--   Just 2
+--   >>>  m !!? (1,-1)
+--   Nothing
+--   >>>  m !!? (1,0)
+--   Just 3
+--   >>>  m !!? (1,1)
+--   Just 4
+--   >>>  m !!? (2,-1)
+--   Nothing
+--   >>>  m !!? (2,0)
+--   Nothing
+--   >>>  m !!? (2,1)
+--   Nothing
+--   >>>  m !!? (2,2)
+--   Nothing
+--
+(!!?) :: (Arity m, Arity n) => Mat m n a -> (Int, Int) -> Maybe a
+(!!?) matrix idx =
+  ifoldl2 f Nothing matrix
+  where
+    f k l found cur = if (k,l) == idx then (Just cur) else found
 
 
 -- | 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.
 --
@@ -216,9 +233,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?
@@ -270,6 +288,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.
@@ -852,5 +871,3 @@ ifoldl2 f initial (Mat rows) =
 
     row_function :: b -> Int -> Vec n a -> b
     row_function rowinit idx r = V.ifoldl (g idx) rowinit r
-
-