X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FVector.hs;h=1aaf78daa42e31ca0eae5e106b52aa3827276882;hb=b5bc0c0bc71915093d0883e2fd923808ada4db6c;hp=97bb5d8aea2c94a8100ed981bff988d6a9df60ab;hpb=d3b5d9259dd9f0650d51447d7f57ed6d782dfdeb;p=numerical-analysis.git diff --git a/src/Vector.hs b/src/Vector.hs index 97bb5d8..1aaf78d 100644 --- a/src/Vector.hs +++ b/src/Vector.hs @@ -1,24 +1,305 @@ +{-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} --- | The 'Vector' class represents elements of a normed vector --- space. We define instances for all common numeric types. module Vector where -import Data.Number.BigFloat +import Data.List (intercalate) +import Data.Vector.Fixed ( + Dim, + Fun(..), + N2, + N3, + N4, + Vector(..), + construct, + inspect, + toList, + ) +import qualified Data.Vector.Fixed as V ( + eq, + foldl, + length, + map, + replicate, + sum, + zipWith + ) -class (Num a) => Vector a where - norm :: RealFrac b => a -> b +import Normed --- Define instances for common numeric types. -instance Vector Integer where - norm = fromInteger +-- | The Vn newtype simply wraps (Vector v a) so that we avoid +-- undecidable instances. +newtype Vn v a = Vn (v a) -instance Vector Rational where - norm = fromRational +-- | Declare the dimension of the wrapper to be the dimension of what +-- it contains. +type instance Dim (Vn v) = Dim v -instance Epsilon e => Vector (BigFloat e) where - norm = fromRational . toRational +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 -instance Vector Double where - norm = fromRational . toRational + -- | Defer to the inspect defined on the contained type. + inspect (Vn v1) = inspect v1 + +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 + + +-- | 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. +-- +-- Examples: +-- +-- >>> let v1 = make2d (1,2) +-- >>> let v2 = make2d (1,2) +-- >>> let v3 = make2d (3,4) +-- >>> v1 == v2 +-- True +-- >>> v1 == v3 +-- False +-- +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 + + -- | 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" + + +-- | 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 + -- | 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. + -- + -- Examples: + -- + -- >>> let v1 = make3d (1,5,2) + -- >>> norm_infty v1 + -- 5 + -- + norm_infty (Vn v1) = realToFrac $ V.foldl max 0 v1 + + -- | Generic p-norms. The usual norm in R^n is (norm_p 2). + -- + -- Examples: + -- + -- >>> let v1 = make2d (3,4) + -- >>> norm_p 1 v1 + -- 7.0 + -- >>> 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) + + +-- | 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: +-- +-- >>> import Roots.Simple +-- >>> let h = 0.5 :: Double +-- >>> let g1 (Vn (Vec2D x y)) = 1.0 + h*exp(-(x^2))/(1.0 + y^2) +-- >>> let g2 (Vn (Vec2D x y)) = 0.5 + h*atan(x^2 + y^2) +-- >>> let g u = make2d ((g1 u), (g2 u)) +-- >>> let u0 = make2d (1.0, 1.0) +-- >>> let eps = 1/(10^9) +-- >>> 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)