From 303c5e7bba583f08e59bc6c848be8e75c1155a3b Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Wed, 20 Feb 2013 21:47:50 -0500 Subject: [PATCH] A huge pile of crap upon Matrix/Vector. --- src/Linear/Matrix.hs | 260 +++++++++++++++++++++++++++++++++------- src/Linear/Vector.hs | 278 +++++++++---------------------------------- 2 files changed, 270 insertions(+), 268 deletions(-) diff --git a/src/Linear/Matrix.hs b/src/Linear/Matrix.hs index 63c0348..39576dc 100644 --- a/src/Linear/Matrix.hs +++ b/src/Linear/Matrix.hs @@ -1,78 +1,146 @@ -{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE ExistentialQuantification #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE RebindableSyntax #-} module Linear.Matrix where +import Data.List (intercalate) + import Data.Vector.Fixed ( Dim, Vector ) import qualified Data.Vector.Fixed as V ( + Fun(..), + N1, + and, + eq, + foldl, fromList, length, map, - toList + maximum, + replicate, + toList, + zipWith ) -import Data.Vector.Fixed.Internal (arity) - +import Data.Vector.Fixed.Internal (Arity, arity, S, Dim) import Linear.Vector +import Normed + +import NumericPrelude hiding (abs) +import qualified Algebra.Algebraic as Algebraic +import qualified Algebra.Absolute as Absolute +import qualified Algebra.Additive as Additive +import qualified Algebra.Ring as Ring +import Algebra.Absolute (abs) +import qualified Algebra.Field as Field +import qualified Algebra.RealField as RealField +import qualified Algebra.RealRing as RealRing +import qualified Algebra.ToRational as ToRational +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 + +-- 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 + -- | Compare a row at a time. + -- + -- Examples: + -- + -- >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int + -- >>> let m2 = fromList [[1,2],[3,4]] :: Mat2 Int + -- >>> let m3 = fromList [[5,6],[7,8]] :: Mat2 Int + -- >>> m1 == m2 + -- True + -- >>> m1 == m3 + -- False + -- + (Mat rows1) == (Mat rows2) = + V.and $ V.zipWith comp rows1 rows2 + where + -- Compare a row, one column at a time. + comp row1 row2 = V.and (V.zipWith (==) row1 row2) + + +instance (Show a, Vector v String, Vector w String) => Show (Mat v w 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 + -- who insists his vector lengths be statically checked at + -- compile-time). + -- + -- Examples: + -- + -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int + -- >>> show m + -- ((1,2),(3,4)) + -- + show (Mat rows) = + "(" ++ (intercalate "," (V.toList row_strings)) ++ ")" + where + row_strings = V.map show_vector rows + show_vector v1 = + "(" ++ (intercalate "," element_strings) ++ ")" + where + v1l = V.toList v1 + element_strings = P.map show v1l + -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) +toList :: Mat v w a -> [[a]] +toList (Mat rows) = map V.toList (V.toList rows) -- | 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 +fromList :: (Vector v (w a), Vector w a, Vector v a) => [[a]] -> Mat v w a +fromList vs = Mat (V.fromList $ map V.fromList vs) -- | Unsafe indexing. -(!!!) :: (Vector v (Vn w a), Vector w a) => Mat v w a -> (Int, Int) -> 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) +(!!?) :: Mat v w a -> (Int, Int) -> Maybe a +(!!?) m@(Mat rows) (i, j) | i < 0 || j < 0 = Nothing - | i > V.length m = Nothing + | i > V.length rows = 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 +nrows :: Mat v w a -> Int +nrows (Mat rows) = V.length rows -- | 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) +ncols :: forall v 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 +row :: Mat v w a -> Int -> w a +row (Mat rows) i = rows ! 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 +column :: (Vector v a) => Mat v w a -> Int -> v a +column (Mat rows) j = + V.map (element j) rows where element = flip (!) @@ -89,16 +157,16 @@ column m j = -- >>> transpose m -- ((1,3),(2,4)) -- -transpose :: (Vector v (Vn w a), - Vector w (Vn v a), +transpose :: (Vector w (v a), Vector v a, Vector w a) => Mat v w a -> Mat w v a -transpose m = V.fromList column_list +transpose m = Mat $ V.fromList column_list where column_list = [ column m i | i <- [0..(ncols m)-1] ] + -- | Is @m@ symmetric? -- -- Examples: @@ -111,7 +179,7 @@ transpose m = V.fromList column_list -- >>> symmetric m2 -- False -- -symmetric :: (Vector v (Vn w a), +symmetric :: (Vector v (w a), Vector w a, v ~ w, Vector w Bool, @@ -136,11 +204,11 @@ symmetric m = -- ((0,1,2),(1,2,3),(2,3,4)) -- construct :: forall v w a. - (Vector v (Vn w a), + (Vector v (w a), Vector w a) => (Int -> Int -> a) -> Mat v w a -construct lambda = rows +construct lambda = Mat rows where -- The arity trick is used in Data.Vector.Fixed.length. imax = (arity (undefined :: Dim v)) - 1 @@ -161,19 +229,21 @@ construct lambda = rows -- ((20.000000000000004,-1.0),(-1.0,20.0)) -- cholesky :: forall a v w. - (RealFloat a, - Vector v (Vn w a), - Vector w a) + (Algebraic.C a, + Vector v (w a), + Vector w a, + Vector v 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) | otherwise = 0 + -- | Matrix multiplication. Our 'Num' instance doesn't define one, and -- we need additional restrictions on the result type anyway. -- @@ -184,12 +254,11 @@ cholesky m = construct r -- >>> m1 `mult` m2 -- ((22,28),(49,64)) -- -mult :: (Num a, - Vector v (Vn w a), +mult :: (Ring.C a, + Vector v a, Vector w a, - Vector w (Vn z a), Vector z a, - Vector v (Vn z a)) + Vector v (z a)) => Mat v w a -> Mat w z a -> Mat v z a @@ -197,3 +266,104 @@ mult m1 m2 = construct lambda where lambda i j = sum [(m1 !!! (i,k)) * (m2 !!! (k,j)) | k <- [0..(ncols m1)-1] ] + + + +instance (Ring.C a, + Vector v (w a), + Vector w a) + => Additive.C (Mat v w a) where + + (Mat rows1) + (Mat rows2) = + Mat $ V.zipWith (V.zipWith (+)) rows1 rows2 + + (Mat rows1) - (Mat rows2) = + Mat $ V.zipWith (V.zipWith (-)) rows1 rows2 + + 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 + one = Mat (V.replicate $ V.replicate (fromInteger 1)) + m1 * m2 = m1 `mult` m2 + + +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 + -- Treat the matrix as a big vector. + norm_p p (Mat rows) = + sqrt $ sum [(fromRational' $ toRational x)^2 | x <- xs] + where + xs = concat $ V.toList $ V.map V.toList rows + + 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 + + + + + +-- Vector helpers. We want it to be easy to create low-dimension +-- column vectors. +type Vec a b = Mat a D1 b + +vec2d :: (a,a) -> Mat D2 D1 a +vec2d (x,y) = Mat (D2 (D1 x) (D1 y)) + +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 D4 D1 a +vec4d (w,x,y,z) = Mat (D4 (D1 w) (D1 x) (D1 y) (D1 z)) + +dot :: (RealRing.C a, + Dim w ~ V.N1, + Vector v a, + Vector w a, + Vector w (v a), + Vector w (w a)) + => Mat v w a + -> Mat v w a + -> a +v1 `dot` v2 = ((transpose v1) `mult` v2) !!! (0, 0) + + +-- | The angle between @v1@ and @v2@ in Euclidean space. +-- +-- Examples: +-- +-- >>> let v1 = make2d (1.0, 0.0) +-- >>> let v2 = make2d (0.0, 1.0) +-- >>> angle v1 v2 == pi/2.0 +-- True +-- +angle :: (Transcendental.C a, + RealRing.C a, + Dim w ~ V.N1, + Vector w (w a), + Vector v [a], + Vector v a, + Vector w a, + Vector v (w a), + Vector w (v a), + ToRational.C a) + => Mat v w a + -> Mat v w a + -> a +angle v1 v2 = + acos theta + where + theta = (recip norms) * (v1 `dot` v2) + norms = (norm v1) * (norm v2) diff --git a/src/Linear/Vector.hs b/src/Linear/Vector.hs index 7cc5e00..9774dcd 100644 --- a/src/Linear/Vector.hs +++ b/src/Linear/Vector.hs @@ -11,6 +11,7 @@ import Data.List (intercalate) import Data.Vector.Fixed ( Dim, Fun(..), + N1, N2, N3, N4, @@ -20,119 +21,74 @@ import Data.Vector.Fixed ( toList, ) import qualified Data.Vector.Fixed as V ( - eq, - foldl, length, - map, - replicate, - sum, - zipWith ) import Normed --- | The Vn newtype simply wraps (Vector v a) so that we avoid --- undecidable instances. -newtype Vn v a = Vn (v a) --- | Declare the dimension of the wrapper to be the dimension of what --- it contains. -type instance Dim (Vn v) = Dim v +-- * Low-dimension vector wrappers. +-- +-- These wrappers are instances of 'Vector', so they inherit all of +-- the userful instances defined above. But, they use fixed +-- constructors, so you can pattern match out the individual +-- components. -instance (Vector v a) => Vector (Vn v) a where - -- | Fortunately, 'Fun' is an instance of 'Functor'. The - -- 'construct' defined on our contained type will return a - -- 'Fun', and we simply slap our constructor on top with fmap. - construct = fmap Vn construct +data D1 a = D1 a +type instance Dim D1 = N1 +instance Vector D1 a where + inspect (D1 x) (Fun f) = f x + construct = Fun D1 - -- | Defer to the inspect defined on the contained type. - inspect (Vn v1) = inspect v1 +data D2 a = D2 a a +type instance Dim D2 = N2 +instance Vector D2 a where + inspect (D2 x y) (Fun f) = f x y + construct = Fun D2 -instance (Show a, Vector v a) => Show (Vn v a) where - -- | Display vectors as ordinary tuples. This is poor practice, but - -- these results are primarily displayed interactively and - -- convenience trumps correctness (said the guy who insists his - -- vector lengths be statically checked at compile-time). - -- - -- Examples: - -- - -- >>> let v1 = make2d (1,2) - -- >>> show v1 - -- (1,2) - -- - show (Vn v1) = - "(" ++ (intercalate "," element_strings) ++ ")" - where - v1l = toList v1 - element_strings = Prelude.map show v1l +data D3 a = D3 a a a +type instance Dim D3 = N3 +instance Vector D3 a where + inspect (D3 x y z) (Fun f) = f x y z + construct = Fun D3 + +data D4 a = D4 a a a a +type instance Dim D4 = N4 +instance Vector D4 a where + inspect (D4 w x y z) (Fun f) = f w x y z + construct = Fun D4 --- | We would really like to say, "anything that is a vector of --- equatable things is itself equatable." The 'Vn' class --- allows us to express this without a GHC battle. +-- | Unsafe indexing. -- -- Examples: -- --- >>> let v1 = make2d (1,2) --- >>> let v2 = make2d (1,2) --- >>> let v3 = make2d (3,4) --- >>> v1 == v2 --- True --- >>> v1 == v3 --- False +-- >>> let v1 = Vec2D 1 2 +-- >>> v1 ! 1 +-- 2 -- -instance (Eq a, Vector v a) => Eq (Vn v a) where - (Vn v1) == (Vn v2) = v1 `V.eq` v2 - - --- | The use of 'Num' here is of course incorrect (otherwise, we --- wouldn't have to throw errors). But it's really nice to be able --- to use normal addition/subtraction. -instance (Num a, Vector v a) => Num (Vn v a) where - -- | Componentwise addition. - -- - -- Examples: - -- - -- >>> let v1 = make2d (1,2) - -- >>> let v2 = make2d (3,4) - -- >>> v1 + v2 - -- (4,6) - -- - (Vn v1) + (Vn v2) = Vn $ V.zipWith (+) v1 v2 - - -- | Componentwise subtraction. - -- - -- Examples: - -- - -- >>> let v1 = make2d (1,2) - -- >>> let v2 = make2d (3,4) - -- >>> v1 - v2 - -- (-2,-2) - -- - (Vn v1) - (Vn v2) = Vn $ V.zipWith (-) v1 v2 +(!) :: (Vector v a) => v a -> Int -> a +(!) v1 idx = (toList v1) !! idx - -- | Create an n-vector whose components are all equal to the given - -- integer. The result type must be specified since otherwise the - -- length n would be unknown. - -- - -- Examples: - -- - -- >>> let v1 = fromInteger 17 :: Vn Vec3 Int - -- (17,17,17) - -- - fromInteger x = Vn $ V.replicate (fromInteger x) - (*) = error "multiplication of vectors is undefined" - abs = error "absolute value of vectors is undefined" - signum = error "signum of vectors is undefined" +-- | Safe indexing. +-- +-- Examples: +-- +-- >>> let v1 = Vec3D 1 2 3 +-- >>> v1 !? 2 +-- Just 3 +-- >>> v1 !? 3 +-- Nothing +-- +(!?) :: (Vector v a) => v a -> Int -> Maybe a +(!?) v1 idx + | idx < 0 || idx >= V.length v1 = Nothing + | otherwise = Just $ v1 ! idx --- | This is probably useless, since the vectors we usually contain --- aren't functor instances. -instance (Functor v) => Functor (Vn v) where - fmap f (Vn v1) = Vn (f `fmap` v1) -instance (RealFloat a, Ord a, Vector v a) => Normed (Vn v a) where +--instance (RealFloat a, Ord a, Vector v a) => Normed (Vn v a) where -- | 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. @@ -143,7 +99,7 @@ instance (RealFloat a, Ord a, Vector v a) => Normed (Vn v a) where -- >>> norm_infty v1 -- 5 -- - norm_infty (Vn v1) = realToFrac $ V.foldl max 0 v1 +-- norm_infty (Vn v1) = realToFrac $ V.foldl max 0 v1 -- | Generic p-norms. The usual norm in R^n is (norm_p 2). -- @@ -155,78 +111,16 @@ instance (RealFloat a, Ord a, Vector v a) => Normed (Vn v a) where -- >>> norm_p 2 v1 -- 5.0 -- - norm_p p (Vn v1) = - realToFrac $ root $ V.sum $ V.map (exponentiate . abs) v1 - where - exponentiate = (** (fromIntegral p)) - root = (** (recip (fromIntegral p))) - --- | Dot (standard inner) product. --- --- Examples: --- --- >>> let v1 = make3d (1,2,3) --- >>> let v2 = make3d (4,5,6) --- >>> dot v1 v2 --- 32 --- -dot :: (Num a, Vector v a) => Vn v a -> Vn v a -> a -dot (Vn v1) (Vn v2) = V.sum $ V.zipWith (*) v1 v2 - - --- | The angle between @v1@ and @v2@ in Euclidean space. --- --- Examples: --- --- >>> let v1 = make2d (1.0, 0.0) --- >>> let v2 = make2d (0.0, 1.0) --- >>> angle v1 v2 == pi/2.0 --- True --- -angle :: (RealFloat a, Vector v a) => Vn v a -> Vn v a -> a -angle v1 v2 = - acos theta - where - theta = (v1 `dot` v2) / norms - norms = (norm v1) * (norm v2) +-- norm_p p (Vn v1) = +-- realToFrac $ root $ V.sum $ V.map (exponentiate . abs) v1 +-- where +-- exponentiate = (** (fromIntegral p)) +-- root = (** (recip (fromIntegral p))) --- | Unsafe indexing. --- --- Examples: --- --- >>> let v1 = make2d (1,2) --- >>> v1 ! 1 --- 2 --- -(!) :: (Vector v a) => v a -> Int -> a -(!) v1 idx = (toList v1) !! idx --- | Safe indexing. --- --- Examples: --- --- >>> let v1 = make3d (1,2,3) --- >>> v1 !? 2 --- Just 3 --- >>> v1 !? 3 --- Nothing --- -(!?) :: (Vector v a) => v a -> Int -> Maybe a -(!?) v1 idx - | idx < 0 || idx >= V.length v1 = Nothing - | otherwise = Just $ v1 ! idx - - --- * Low-dimension vector wrappers. --- --- These wrappers are instances of 'Vector', so they inherit all of --- the userful instances defined above. But, they use fixed --- constructors, so you can pattern match out the individual --- components. - -- | Convenient constructor for 2D vectors. -- -- Examples: @@ -241,65 +135,3 @@ angle v1 v2 = -- >>> fixed_point g eps u0 -- (1.0728549599342185,1.0820591495686167) -- -data Vec2D a = Vec2D a a -type instance Dim Vec2D = N2 -instance Vector Vec2D a where - inspect (Vec2D x y) (Fun f) = f x y - construct = Fun Vec2D - -data Vec3D a = Vec3D a a a -type instance Dim Vec3D = N3 -instance Vector Vec3D a where - inspect (Vec3D x y z) (Fun f) = f x y z - construct = Fun Vec3D - -data Vec4D a = Vec4D a a a a -type instance Dim Vec4D = N4 -instance Vector Vec4D a where - inspect (Vec4D w x y z) (Fun f) = f w x y z - construct = Fun Vec4D - - --- | Convenience function for creating 2d vectors. --- --- Examples: --- --- >>> let v1 = make2d (1,2) --- >>> v1 --- (1,2) --- >>> let Vn (Vec2D x y) = v1 --- >>> (x,y) --- (1,2) --- -make2d :: forall a. (a,a) -> Vn Vec2D a -make2d (x,y) = Vn (Vec2D x y) - - --- | Convenience function for creating 3d vectors. --- --- Examples: --- --- >>> let v1 = make3d (1,2,3) --- >>> v1 --- (1,2,3) --- >>> let Vn (Vec3D x y z) = v1 --- >>> (x,y,z) --- (1,2,3) --- -make3d :: forall a. (a,a,a) -> Vn Vec3D a -make3d (x,y,z) = Vn (Vec3D x y z) - - --- | Convenience function for creating 4d vectors. --- --- Examples: --- --- >>> let v1 = make4d (1,2,3,4) --- >>> v1 --- (1,2,3,4) --- >>> let Vn (Vec4D w x y z) = v1 --- >>> (w,x,y,z) --- (1,2,3,4) --- -make4d :: forall a. (a,a,a,a) -> Vn Vec4D a -make4d (w,x,y,z) = Vn (Vec4D w x y z) -- 2.44.2