From: Michael Orlitzky Date: Tue, 5 Feb 2013 01:25:14 +0000 (-0500) Subject: Remove non-fixed Matrix module. X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=c00e5ae9829a358890922083267564ec13798061;p=numerical-analysis.git Remove non-fixed Matrix module. Drop the "Fixed" prefix for Vector/Matrix modules. Update ghci includes. --- diff --git a/.ghci b/.ghci index c9819c7..47792b8 100644 --- a/.ghci +++ b/.ghci @@ -6,10 +6,11 @@ :load src/Aliases.hs src/Integration/Simpson.hs src/Integration/Trapezoid.hs + src/Matrix.hs src/Misc.hs src/ODE/IVP.hs src/Roots/Simple.hs - src/TwoTuple.hs + src/Vector.hs :} -- Just for convenience. @@ -18,10 +19,11 @@ import Data.Number.BigFloat import Aliases import Integration.Simpson import Integration.Trapezoid +import Matrix import Misc import ODE.IVP import Roots.Simple -import TwoTuple +import Vector -- Use a calmer prompt. :set prompt "numerical-analysis> " diff --git a/src/FixedMatrix.hs b/src/FixedMatrix.hs deleted file mode 100644 index e3a1e2b..0000000 --- a/src/FixedMatrix.hs +++ /dev/null @@ -1,200 +0,0 @@ -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE TypeFamilies #-} - -module FixedMatrix -where - -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) -type Mat2 a = Mat Vec2D Vec2D a -type Mat3 a = Mat Vec3D Vec3D a -type Mat4 a = Mat Vec4D Vec4D a - --- | Convert a matrix to a nested list. -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 :: (Vector v (Vn w a), Vector w a) => [[a]] -> Mat v w a -fromList vs = V.fromList $ map V.fromList vs - - --- | Unsafe indexing. -(!!!) :: (Vector v (Vn w a), Vector w a) => Mat v w a -> (Int, Int) -> a -(!!!) m (i, j) = (row m i) ! j - --- | Safe indexing. -(!!?) :: (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) ! j - - --- | The number of rows in the matrix. -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. (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 :: (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 :: (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 (!) - - --- | 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. --- --- 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 - where - column_list = [ column m i | i <- [0..(ncols m)-1] ] - --- | Is @m@ symmetric? --- --- 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, - Vector w Bool, - Eq a) - => Mat v w a - -> Bool -symmetric m = - m == (transpose m) - - --- | Construct a new matrix from a function @lambda@. The function --- @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. - (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 :: 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] ] diff --git a/src/Matrix.hs b/src/Matrix.hs index 349b71a..ea44ae9 100644 --- a/src/Matrix.hs +++ b/src/Matrix.hs @@ -1,144 +1,200 @@ {-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE TypeFamilies #-} --- | A Matrix type using Data.Vector as the underlying type. In other --- words, the size is not fixed, but at least we have safe indexing if --- we want it. --- --- This should be replaced with a fixed-size implementation eventually! --- module Matrix where -import qualified Data.Vector as V - -type Rows a = V.Vector (V.Vector a) -type Columns a = V.Vector (V.Vector a) -data Matrix a = Matrix (Rows a) deriving Eq - --- | Unsafe indexing -(!) :: (Matrix a) -> (Int, Int) -> a -(Matrix rows) ! (i, j) = (rows V.! i) V.! j +import Vector +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) +type Mat2 a = Mat Vec2D Vec2D a +type Mat3 a = Mat Vec3D Vec3D a +type Mat4 a = Mat Vec4D Vec4D a + +-- | Convert a matrix to a nested list. +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 :: (Vector v (Vn w a), Vector w a) => [[a]] -> Mat v w a +fromList vs = V.fromList $ map V.fromList vs + + +-- | Unsafe indexing. +(!!!) :: (Vector v (Vn w a), Vector w a) => Mat v w a -> (Int, Int) -> a +(!!!) m (i, j) = (row m i) ! j + +-- | Safe indexing. +(!!?) :: (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) ! j --- | Safe indexing -(!?) :: (Matrix a) -> (Int, Int) -> Maybe a -(Matrix rows) !? (i, j) = do - row <- rows V.!? i - col <- row V.!? j - return col --- | Unsafe indexing without bounds checking -unsafeIndex :: (Matrix a) -> (Int, Int) -> a -(Matrix rows) `unsafeIndex` (i, j) = - (rows `V.unsafeIndex` i) `V.unsafeIndex` j +-- | The number of rows in the matrix. +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. (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 :: (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 :: (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 (!) --- | Return the @i@th column of @m@. Unsafe! -column :: (Matrix a) -> Int -> (V.Vector a) -column (Matrix rows) i = - V.fromList [row V.! i | row <- V.toList rows] --- | The number of rows in the matrix. -nrows :: (Matrix a) -> Int -nrows (Matrix rows) = V.length rows - --- | The number of columns in the first row of the matrix. -ncols :: (Matrix a) -> Int -ncols (Matrix rows) - | V.length rows == 0 = 0 - | otherwise = V.length (rows V.! 0) - --- | Return the vector of @m@'s columns. -columns :: (Matrix a) -> (Columns a) -columns m = - V.fromList [column m i | i <- [0..(ncols m)-1]] - --- | Transose @m@; switch it's columns and its rows. -transpose :: (Matrix a) -> (Matrix a) -transpose m = - Matrix (columns m) - -instance Show a => Show (Matrix a) where - show (Matrix rows) = - concat $ V.toList $ V.map show_row rows - where show_row r = "[" ++ (show r) ++ "]\n" - -instance Functor Matrix where - f `fmap` (Matrix rows) = Matrix (V.map (fmap f) rows) - - --- | Vector addition. -vplus :: Num a => (V.Vector a) -> (V.Vector a) -> (V.Vector a) -vplus xs ys = V.zipWith (+) xs ys - --- | Vector subtraction. -vminus :: Num a => (V.Vector a) -> (V.Vector a) -> (V.Vector a) -vminus xs ys = V.zipWith (-) xs ys - --- | Add two vectors of rows. -rowsplus :: Num a => (Rows a) -> (Rows a) -> (Rows a) -rowsplus rs1 rs2 = - V.zipWith vplus rs1 rs2 - --- | Subtract two vectors of rows. -rowsminus :: Num a => (Rows a) -> (Rows a) -> (Rows a) -rowsminus rs1 rs2 = - V.zipWith vminus rs1 rs2 - --- | Matrix multiplication. -mtimes :: Num a => (Matrix a) -> (Matrix a) -> (Matrix a) -mtimes m1 m2 = - Matrix (V.fromList rows) +-- | 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. +-- +-- 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 where - row i = - V.fromList [ sum [ (m1 ! (i,k)) * (m2 ! (k,j)) | k <- [0..(ncols m1)-1] ] - | j <- [0..(ncols m2)-1] ] - rows = [row i | i <- [0..(nrows m1)-1]] + column_list = [ column m i | i <- [0..(ncols m)-1] ] -- | Is @m@ symmetric? -symmetric :: Eq a => (Matrix a) -> Bool +-- +-- 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, + Vector w Bool, + Eq a) + => Mat v w a + -> Bool symmetric m = m == (transpose m) + -- | Construct a new matrix from a function @lambda@. The function -- @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. -- --- The @imax@ and @jmax@ parameters determine the size of the matrix. +-- TODO: Don't cheat with fromList. +-- +-- Examples: -- -construct :: Int -> Int -> (Int -> Int -> a) -> (Matrix a) -construct imax jmax lambda = - Matrix rows +-- >>> 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 (Vn w a), + Vector w a) + => (Int -> Int -> a) + -> Mat v w a +construct lambda = rows where - row i = V.fromList [ lambda i j | j <- [0..jmax] ] - rows = V.fromList [ row i | i <- [0..imax] ] + -- 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] ] -- | 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. -cholesky :: forall a. RealFloat a => (Matrix a) -> (Matrix a) -cholesky m = - construct (nrows m - 1) (ncols m - 1) r +-- +-- 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]]) + 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) + (((m !!! (i,j)) - sum [(r k i)*(r k j) | k <- [0..i-1]]))/(r i i) | otherwise = 0 --- | It's not correct to use Num here, but I really don't want to have --- to define my own addition and subtraction. -instance Num a => Num (Matrix a) where - -- Standard componentwise addition. - (Matrix rows1) + (Matrix rows2) = Matrix (rows1 `rowsplus` rows2) - - -- Standard componentwise subtraction. - (Matrix rows1) - (Matrix rows2) = Matrix (rows1 `rowsminus` rows2) - - -- Matrix multiplication. - m1 * m2 = m1 `mtimes` m2 - - abs _ = error "absolute value of matrices is undefined" - - signum _ = error "signum of matrices is undefined" - - fromInteger _ = error "fromInteger of matrices is undefined" +-- | 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] ] diff --git a/src/FixedVector.hs b/src/Vector.hs similarity index 99% rename from src/FixedVector.hs rename to src/Vector.hs index 1bb828d..b42932c 100644 --- a/src/FixedVector.hs +++ b/src/Vector.hs @@ -4,7 +4,7 @@ {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeFamilies #-} -module FixedVector +module Vector where import Data.List (intercalate)