From 55cf46834f181be01b17b4c5a02ecd772c4e3090 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Mon, 4 Feb 2013 20:23:33 -0500 Subject: [PATCH] Clean up imports in FixedMatrix.hs. Add matrix multiplication and Cholesky factorization to FixedMatrix.hs. Add some doctests to FixedMatrix.hs. --- src/FixedMatrix.hs | 162 +++++++++++++++++++++++++++++++++++---------- 1 file changed, 126 insertions(+), 36 deletions(-) diff --git a/src/FixedMatrix.hs b/src/FixedMatrix.hs index 7b34fee..e3a1e2b 100644 --- a/src/FixedMatrix.hs +++ b/src/FixedMatrix.hs @@ -7,8 +7,19 @@ module FixedMatrix where -import FixedVector as FV -import qualified Data.Vector.Fixed as V +import FixedVector +import Data.Vector.Fixed ( + Arity(..), + Dim, + Vector, + (!), + ) +import qualified Data.Vector.Fixed as V ( + fromList, + length, + map, + toList + ) import Data.Vector.Fixed.Internal (arity) type Mat v w a = Vn v (Vn w a) @@ -17,63 +28,72 @@ type Mat3 a = Mat Vec3D Vec3D a type Mat4 a = Mat Vec4D Vec4D a -- | Convert a matrix to a nested list. -toList :: (V.Vector v (Vn w a), V.Vector w a) => Mat v w a -> [[a]] -toList m = Prelude.map V.toList (V.toList m) +toList :: (Vector v (Vn w a), Vector w a) => Mat v w a -> [[a]] +toList m = map V.toList (V.toList m) -- | Create a matrix from a nested list. -fromList :: (V.Vector v (Vn w a), V.Vector w a) => [[a]] -> Mat v w a -fromList vs = V.fromList $ Prelude.map V.fromList vs +fromList :: (Vector v (Vn w a), Vector w a) => [[a]] -> Mat v w a +fromList vs = V.fromList $ map V.fromList vs -- | Unsafe indexing. -(!) :: (V.Vector v (Vn w a), V.Vector w a) => Mat v w a -> (Int, Int) -> a -(!) m (i, j) = (row m i) V.! j +(!!!) :: (Vector v (Vn w a), Vector w a) => Mat v w a -> (Int, Int) -> a +(!!!) m (i, j) = (row m i) ! j -- | Safe indexing. -(!?) :: (V.Vector v (Vn w a), V.Vector w a) => Mat v w a - -> (Int, Int) - -> Maybe a -(!?) m (i, j) +(!!?) :: (Vector v (Vn w a), Vector w a) => Mat v w a + -> (Int, Int) + -> Maybe a +(!!?) m (i, j) | i < 0 || j < 0 = Nothing | i > V.length m = Nothing | otherwise = if j > V.length (row m j) then Nothing - else Just $ (row m j) V.! j + else Just $ (row m j) ! j -- | The number of rows in the matrix. -nrows :: forall v w a. (V.Vector v (Vn w a), V.Vector w a) => Mat v w a -> Int +nrows :: forall v w a. (Vector v (Vn w a), Vector w a) => Mat v w a -> Int nrows = V.length -- | The number of columns in the first row of the -- matrix. Implementation stolen from Data.Vector.Fixed.length. -ncols :: forall v w a. (V.Vector v (Vn w a), V.Vector w a) => Mat v w a -> Int -ncols _ = arity (undefined :: V.Dim w) +ncols :: forall v w a. (Vector v (Vn w a), Vector w a) => Mat v w a -> Int +ncols _ = arity (undefined :: Dim w) -- | Return the @i@th row of @m@. Unsafe. -row :: (V.Vector v (Vn w a), V.Vector w a) => Mat v w a - -> Int - -> Vn w a -row m i = m V.! i +row :: (Vector v (Vn w a), Vector w a) => Mat v w a + -> Int + -> Vn w a +row m i = m ! i -- | Return the @j@th column of @m@. Unsafe. -column :: (V.Vector v a, V.Vector v (Vn w a), V.Vector w a) => Mat v w a - -> Int - -> Vn v a +column :: (Vector v a, Vector v (Vn w a), Vector w a) => Mat v w a + -> Int + -> Vn v a column m j = V.map (element j) m where - element = flip (V.!) + element = flip (!) --- | Transose @m@; switch it's columns and its rows. This is a dirty +-- | 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 :: (V.Vector v (Vn w a), - V.Vector w (Vn v a), - V.Vector v a, - V.Vector w a) +-- +-- TODO: Don't cheat with fromList. +-- +-- Examples: +-- +-- >>> let m = fromList [[1,2], [3,4]] :: Mat2 Int +-- >>> transpose m +-- ((1,3),(2,4)) +-- +transpose :: (Vector v (Vn w a), + Vector w (Vn v a), + Vector v a, + Vector w a) => Mat v w a -> Mat w v a transpose m = V.fromList column_list @@ -81,10 +101,21 @@ transpose m = V.fromList column_list column_list = [ column m i | i <- [0..(ncols m)-1] ] -- | Is @m@ symmetric? -symmetric :: (V.Vector v (Vn w a), - V.Vector w a, +-- +-- Examples: +-- +-- >>> let m1 = fromList [[1,2], [2,1]] :: Mat2 Int +-- >>> symmetric m1 +-- True +-- +-- >>> let m2 = fromList [[1,2], [3,1]] :: Mat2 Int +-- >>> symmetric m2 +-- False +-- +symmetric :: (Vector v (Vn w a), + Vector w a, v ~ w, - V.Vector w Bool, + Vector w Bool, Eq a) => Mat v w a -> Bool @@ -96,15 +127,74 @@ symmetric m = -- @lambda@ should take two parameters i,j corresponding to the -- 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)) +-- construct :: forall v w a. - (V.Vector v (Vn w a), - V.Vector w a) + (Vector v (Vn w a), + Vector w a) => (Int -> Int -> a) -> Mat v w a construct lambda = rows where -- The arity trick is used in Data.Vector.Fixed.length. - imax = (arity (undefined :: V.Dim v)) - 1 - jmax = (arity (undefined :: V.Dim w)) - 1 + 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] ] + +-- | 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. +-- +-- Examples: +-- +-- >>> let m1 = fromList [[20,-1], [-1,20]] :: Mat2 Double +-- >>> cholesky m1 +-- ((4.47213595499958,-0.22360679774997896),(0.0,4.466542286825459)) +-- >>> (transpose (cholesky m1)) `mult` (cholesky m1) +-- ((20.000000000000004,-1.0),(-1.0,20.0)) +-- +cholesky :: forall a v w. + (RealFloat a, + Vector v (Vn w a), + Vector w a) + => (Mat v w a) + -> (Mat v w a) +cholesky m = construct r + where + r :: Int -> Int -> a + r i j | i == j = sqrt(m !!! (i,j) - sum [(r k i)**2 | k <- [0..i-1]]) + | i < j = + (((m !!! (i,j)) - sum [(r k i)*(r k j) | k <- [0..i-1]]))/(r i i) + | otherwise = 0 + +-- | Matrix multiplication. Our 'Num' instance doesn't define one, and +-- we need additional restrictions on the result type anyway. +-- +-- Examples: +-- +-- >>> let m1 = fromList [[1,2,3], [4,5,6]] :: Mat Vec2D Vec3D Int +-- >>> let m2 = fromList [[1,2],[3,4],[5,6]] :: Mat Vec3D Vec2D Int +-- >>> m1 `mult` m2 +-- ((22,28),(49,64)) +-- +mult :: (Num a, + Vector v (Vn w a), + Vector w a, + Vector w (Vn z a), + Vector z a, + Vector v (Vn z a)) + => Mat v w a + -> Mat w z a + -> Mat v z a +mult m1 m2 = construct lambda + where + lambda i j = + sum [(m1 !!! (i,k)) * (m2 !!! (k,j)) | k <- [0..(ncols m1)-1] ] -- 2.44.2