1 {-# LANGUAGE FlexibleContexts #-}
2 {-# LANGUAGE FlexibleInstances #-}
3 {-# LANGUAGE MultiParamTypeClasses #-}
4 {-# LANGUAGE ScopedTypeVariables #-}
5 {-# LANGUAGE TypeFamilies #-}
10 import Data.List (intercalate)
11 import Data.Vector.Fixed (
22 import qualified Data.Vector.Fixed as V (
34 -- | The Vn newtype simply wraps (Vector v a) so that we avoid
35 -- undecidable instances.
36 newtype Vn v a = Vn (v a)
38 -- | Declare the dimension of the wrapper to be the dimension of what
40 type instance Dim (Vn v) = Dim v
42 instance (Vector v a) => Vector (Vn v) a where
43 -- | Fortunately, 'Fun' is an instance of 'Functor'. The
44 -- 'construct' defined on our contained type will return a
45 -- 'Fun', and we simply slap our constructor on top with fmap.
46 construct = fmap Vn construct
48 -- | Defer to the inspect defined on the contained type.
49 inspect (Vn v1) = inspect v1
51 instance (Show a, Vector v a) => Show (Vn v a) where
52 -- | Display vectors as ordinary tuples. This is poor practice, but
53 -- these results are primarily displayed interactively and
54 -- convenience trumps correctness (said the guy who insists his
55 -- vector lengths be statically checked at compile-time).
59 -- >>> let v1 = make2d (1,2)
64 "(" ++ (intercalate "," element_strings) ++ ")"
67 element_strings = Prelude.map show v1l
70 -- | We would really like to say, "anything that is a vector of
71 -- equatable things is itself equatable." The 'Vn' class
72 -- allows us to express this without a GHC battle.
76 -- >>> let v1 = make2d (1,2)
77 -- >>> let v2 = make2d (1,2)
78 -- >>> let v3 = make2d (3,4)
84 instance (Eq a, Vector v a) => Eq (Vn v a) where
85 (Vn v1) == (Vn v2) = v1 `V.eq` v2
88 -- | The use of 'Num' here is of course incorrect (otherwise, we
89 -- wouldn't have to throw errors). But it's really nice to be able
90 -- to use normal addition/subtraction.
91 instance (Num a, Vector v a) => Num (Vn v a) where
92 -- | Componentwise addition.
96 -- >>> let v1 = make2d (1,2)
97 -- >>> let v2 = make2d (3,4)
101 (Vn v1) + (Vn v2) = Vn $ V.zipWith (+) v1 v2
103 -- | Componentwise subtraction.
107 -- >>> let v1 = make2d (1,2)
108 -- >>> let v2 = make2d (3,4)
112 (Vn v1) - (Vn v2) = Vn $ V.zipWith (-) v1 v2
114 -- | Create an n-vector whose components are all equal to the given
115 -- integer. The result type must be specified since otherwise the
116 -- length n would be unknown.
120 -- >>> let v1 = fromInteger 17 :: Vn Vec3 Int
123 fromInteger x = Vn $ V.replicate (fromInteger x)
124 (*) = error "multiplication of vectors is undefined"
125 abs = error "absolute value of vectors is undefined"
126 signum = error "signum of vectors is undefined"
129 -- | This is probably useless, since the vectors we usually contain
130 -- aren't functor instances.
131 instance (Functor v) => Functor (Vn v) where
132 fmap f (Vn v1) = Vn (f `fmap` v1)
135 instance (RealFloat a, Ord a, Vector v a) => Normed (Vn v a) where
136 -- | The infinity norm. We don't use V.maximum here because it
137 -- relies on a type constraint that the vector be non-empty and I
138 -- don't know how to pattern match it away.
142 -- >>> let v1 = make3d (1,5,2)
146 norm_infty (Vn v1) = realToFrac $ V.foldl max 0 v1
148 -- | Generic p-norms. The usual norm in R^n is (norm_p 2).
152 -- >>> let v1 = make2d (3,4)
159 realToFrac $ root $ V.sum $ V.map (exponentiate . abs) v1
161 exponentiate = (** (fromIntegral p))
162 root = (** (recip (fromIntegral p)))
164 -- | Dot (standard inner) product.
168 -- >>> let v1 = make3d (1,2,3)
169 -- >>> let v2 = make3d (4,5,6)
173 dot :: (Num a, Vector v a) => Vn v a -> Vn v a -> a
174 dot (Vn v1) (Vn v2) = V.sum $ V.zipWith (*) v1 v2
177 -- | The angle between @v1@ and @v2@ in Euclidean space.
181 -- >>> let v1 = make2d (1.0, 0.0)
182 -- >>> let v2 = make2d (0.0, 1.0)
183 -- >>> angle v1 v2 == pi/2.0
186 angle :: (RealFloat a, Vector v a) => Vn v a -> Vn v a -> a
190 theta = (v1 `dot` v2) / norms
191 norms = (norm v1) * (norm v2)
194 -- | Unsafe indexing.
198 -- >>> let v1 = make2d (1,2)
202 (!) :: (Vector v a) => v a -> Int -> a
203 (!) v1 idx = (toList v1) !! idx
209 -- >>> let v1 = make3d (1,2,3)
215 (!?) :: (Vector v a) => v a -> Int -> Maybe a
217 | idx < 0 || idx >= V.length v1 = Nothing
218 | otherwise = Just $ v1 ! idx
223 -- * Low-dimension vector wrappers.
225 -- These wrappers are instances of 'Vector', so they inherit all of
226 -- the userful instances defined above. But, they use fixed
227 -- constructors, so you can pattern match out the individual
230 -- | Convenient constructor for 2D vectors.
234 -- >>> import Roots.Simple
235 -- >>> let h = 0.5 :: Double
236 -- >>> let g1 (Vn (Vec2D x y)) = 1.0 + h*exp(-(x^2))/(1.0 + y^2)
237 -- >>> let g2 (Vn (Vec2D x y)) = 0.5 + h*atan(x^2 + y^2)
238 -- >>> let g u = make2d ((g1 u), (g2 u))
239 -- >>> let u0 = make2d (1.0, 1.0)
240 -- >>> let eps = 1/(10^9)
241 -- >>> fixed_point g eps u0
242 -- (1.0728549599342185,1.0820591495686167)
244 data Vec2D a = Vec2D a a
245 type instance Dim Vec2D = N2
246 instance Vector Vec2D a where
247 inspect (Vec2D x y) (Fun f) = f x y
248 construct = Fun Vec2D
250 data Vec3D a = Vec3D a a a
251 type instance Dim Vec3D = N3
252 instance Vector Vec3D a where
253 inspect (Vec3D x y z) (Fun f) = f x y z
254 construct = Fun Vec3D
256 data Vec4D a = Vec4D a a a a
257 type instance Dim Vec4D = N4
258 instance Vector Vec4D a where
259 inspect (Vec4D w x y z) (Fun f) = f w x y z
260 construct = Fun Vec4D
263 -- | Convenience function for creating 2d vectors.
267 -- >>> let v1 = make2d (1,2)
270 -- >>> let Vn (Vec2D x y) = v1
274 make2d :: forall a. (a,a) -> Vn Vec2D a
275 make2d (x,y) = Vn (Vec2D x y)
278 -- | Convenience function for creating 3d vectors.
282 -- >>> let v1 = make3d (1,2,3)
285 -- >>> let Vn (Vec3D x y z) = v1
289 make3d :: forall a. (a,a,a) -> Vn Vec3D a
290 make3d (x,y,z) = Vn (Vec3D x y z)
293 -- | Convenience function for creating 4d vectors.
297 -- >>> let v1 = make4d (1,2,3,4)
300 -- >>> let Vn (Vec4D w x y z) = v1
304 make4d :: forall a. (a,a,a,a) -> Vn Vec4D a
305 make4d (w,x,y,z) = Vn (Vec4D w x y z)