src/FixedVector.hs

   1 {-# LANGUAGE FlexibleContexts #-}
   2 {-# LANGUAGE FlexibleInstances #-}
   3 {-# LANGUAGE MultiParamTypeClasses #-}
   4 {-# LANGUAGE ScopedTypeVariables #-}
   5 {-# LANGUAGE TypeFamilies #-}
   6
   7 module FixedVector
   8 where
   9
  10 import Data.List (intercalate)
  11 import Data.Vector.Fixed (
  12   Dim,
  13   Fun(..),
  14   N1,
  15   N2,
  16   N3,
  17   N4,
  18   Vector(..),
  19   (!),
  20   construct,
  21   inspect,
  22   toList,
  23   )
  24 import qualified Data.Vector.Fixed as V (
  25   foldl,
  26   length,
  27   map,
  28   replicate,
  29   sum,
  30   zipWith
  31   )
  32
  33 import Normed
  34
  35 -- | The Vn newtype simply wraps (Vector v a) so that we avoid
  36 --   undecidable instances.
  37 newtype Vn v a = Vn (v a)
  38
  39 -- | Declare the dimension of the wrapper to be the dimension of what
  40 --   it contains.
  41 type instance Dim (Vn v) = Dim v
  42
  43 instance (Vector v a) => Vector (Vn v) a where
  44   -- | Fortunately, 'Fun' is an instance of 'Functor'. The
  45   --   'construct' defined on our contained type will return a
  46   --   'Fun', and we simply slap our constructor on top with fmap.
  47   construct = fmap Vn construct
  48
  49   -- | Defer to the inspect defined on the contained type.
  50   inspect (Vn v1) = inspect v1
  51
  52 instance (Show a, Vector v a) => Show (Vn v a) where
  53   -- | Display vectors as ordinary tuples. This is poor practice, but
  54   --   these results are primarily displayed interactively and
  55   --   convenience trumps correctness (said the guy who insists his
  56   --   vector lengths be statically checked at compile-time).
  57   --
  58   --   Examples:
  59   --
  60   --   >>> let v1 = make2d (1,2)
  61   --   >>> show v1
  62   --   (1,2)
  63   --
  64   show (Vn v1) =
  65     "(" ++ (intercalate "," element_strings) ++ ")"
  66     where
  67       v1l = toList v1
  68       element_strings = Prelude.map show v1l
  69
  70
  71 -- | We would really like to say, "anything that is a vector of
  72 --   equatable things is itself equatable." The 'Vn' class
  73 --   allows us to express this without a GHC battle.
  74 --
  75 --   Examples:
  76 --
  77 --   >>> let v1 = make2d (1,2)
  78 --   >>> let v2 = make2d (1,2)
  79 --   >>> let v3 = make2d (3,4)
  80 --   >>> v1 == v2
  81 --   True
  82 --   >>> v1 == v3
  83 --   False
  84 --
  85 instance (Eq a, Vector v a, Vector v Bool) => Eq (Vn v a) where
  86   (Vn v1) == (Vn v2) = V.foldl (&&) True (V.zipWith (==) v1 v2)
  87
  88
  89 -- | The use of 'Num' here is of course incorrect (otherwise, we
  90 --   wouldn't have to throw errors). But it's really nice to be able
  91 --   to use normal addition/subtraction.
  92 instance (Num a, Vector v a) => Num (Vn v a) where
  93   -- | Componentwise addition.
  94   --
  95   --   Examples:
  96   --
  97   --   >>> let v1 = make2d (1,2)
  98   --   >>> let v2 = make2d (3,4)
  99   --   >>> v1 + v2
 100   --   (4,6)
 101   --
 102   (Vn v1) + (Vn v2) = Vn $ V.zipWith (+) v1 v2
 103
 104   -- | Componentwise subtraction.
 105   --
 106   --   Examples:
 107   --
 108   --   >>> let v1 = make2d (1,2)
 109   --   >>> let v2 = make2d (3,4)
 110   --   >>> v1 - v2
 111   --   (-2,-2)
 112   --
 113   (Vn v1) - (Vn v2) = Vn $ V.zipWith (-) v1 v2
 114
 115   -- | Create an n-vector whose components are all equal to the given
 116   --   integer. The result type must be specified since otherwise the
 117   --   length n would be unknown.
 118   --
 119   --   Examples:
 120   --
 121   --   >>> let v1 = fromInteger 17 :: Vn Vec3 Int
 122   --   (17,17,17)
 123   --
 124   fromInteger x = Vn $ V.replicate (fromInteger x)
 125   (*) = error "multiplication of vectors is undefined"
 126   abs = error "absolute value of vectors is undefined"
 127   signum = error "signum of vectors is undefined"
 128
 129
 130 -- | This is probably useless, since the vectors we usually contain
 131 --   aren't functor instances.
 132 instance (Functor v) => Functor (Vn v) where
 133   fmap f (Vn v1) = Vn (f `fmap` v1)
 134
 135
 136 instance (RealFloat a, Ord a, Vector v a) => Normed (Vn v a) where
 137   -- | The infinity norm. We don't use V.maximum here because it
 138   --   relies on a type constraint that the vector be non-empty and I
 139   --   don't know how to pattern match it away.
 140   --
 141   --   Examples:
 142   --
 143   --   >>> let v1 = make3d (1,5,2)
 144   --   >>> norm_infty v1
 145   --   5
 146   --
 147   norm_infty (Vn v1) = fromRational $ toRational $ V.foldl max 0 v1
 148
 149   -- | Generic p-norms. The usual norm in R^n is (norm_p 2).
 150   --
 151   --   Examples:
 152   --
 153   --   >>> let v1 = make2d (3,4)
 154   --   >>> norm_p 1 v1
 155   --   7.0
 156   --   >>> norm_p 2 v1
 157   --   5.0
 158   --
 159   norm_p p (Vn v1) =
 160     fromRational $ toRational $ root $ V.sum $ V.map (exponentiate . abs) v1
 161     where
 162       exponentiate = (** (fromIntegral p))
 163       root = (** (recip (fromIntegral p)))
 164
 165 -- | Dot (standard inner) product.
 166 --
 167 --   Examples:
 168 --
 169 --   >>> let v1 = make3d (1,2,3)
 170 --   >>> let v2 = make3d (4,5,6)
 171 --   >>> dot v1 v2
 172 --   32
 173 --
 174 dot :: (Num a, Vector v a) => Vn v a -> Vn v a -> a
 175 dot (Vn v1) (Vn v2) = V.sum $ V.zipWith (*) v1 v2
 176
 177
 178 -- | The angle between @v1@ and @v2@ in Euclidean space.
 179 --
 180 --   Examples:
 181 --
 182 --   >>> let v1 = make2d (1.0, 0.0)
 183 --   >>> let v2 = make2d (0.0, 1.0)
 184 --   >>> angle v1 v2 == pi/2.0
 185 --   True
 186 --
 187 angle :: (RealFloat a, Vector v a) => Vn v a -> Vn v a -> a
 188 angle v1 v2 =
 189   acos theta
 190   where
 191     theta = (v1 `dot` v2) / norms
 192     norms = (norm v1) * (norm v2)
 193
 194
 195 -- | Safe indexing.
 196 --
 197 --   Examples:
 198 --
 199 --   >>> let v1 = make3d (1,2,3)
 200 --   >>> v1 !? 2
 201 --   Just 3
 202 --   >>> v1 !? 3
 203 --   Nothing
 204 --
 205 (!?) :: (Vector v a) => v a -> Int -> Maybe a
 206 (!?) v1 idx
 207   | idx < 0 || idx >= V.length v1 = Nothing
 208   | otherwise                     = Just $ v1 ! idx
 209
 210
 211
 212
 213 -- * Low-dimension vector wrappers.
 214 --
 215 -- These wrappers are instances of 'Vector', so they inherit all of
 216 -- the userful instances defined above. But, they use fixed
 217 -- constructors, so you can pattern match out the individual
 218 -- components.
 219
 220 data Vec2D a = Vec2D a a
 221 type instance Dim Vec2D = N2
 222 instance Vector Vec2D a where
 223   inspect (Vec2D x y) (Fun f) = f x y
 224   construct = Fun Vec2D
 225
 226 data Vec3D a = Vec3D a a a
 227 type instance Dim Vec3D = N3
 228 instance Vector Vec3D a where
 229   inspect (Vec3D x y z) (Fun f) = f x y z
 230   construct = Fun Vec3D
 231
 232 data Vec4D a = Vec4D a a a a
 233 type instance Dim Vec4D = N4
 234 instance Vector Vec4D a where
 235   inspect (Vec4D w x y z) (Fun f) = f w x y z
 236   construct = Fun Vec4D
 237
 238
 239 -- | Convenience function for creating 2d vectors.
 240 --
 241 --   Examples:
 242 --
 243 --   >>> let v1 = make2d (1,2)
 244 --   >>> v1
 245 --   (1,2)
 246 --   >>> let Vn (Vec2D x y) = v1
 247 --   >>> (x,y)
 248 --   (1,2)
 249 --
 250 make2d :: forall a. (a,a) -> Vn Vec2D a
 251 make2d (x,y) = Vn (Vec2D x y)
 252
 253
 254 -- | Convenience function for creating 3d vectors.
 255 --
 256 --   Examples:
 257 --
 258 --   >>> let v1 = make3d (1,2,3)
 259 --   >>> v1
 260 --   (1,2,3)
 261 --   >>> let Vn (Vec3D x y z) = v1
 262 --   >>> (x,y,z)
 263 --   (1,2,3)
 264 --
 265 make3d :: forall a. (a,a,a) -> Vn Vec3D a
 266 make3d (x,y,z) = Vn (Vec3D x y z)
 267
 268
 269 -- | Convenience function for creating 4d vectors.
 270 --
 271 --   Examples:
 272 --
 273 --   >>> let v1 = make4d (1,2,3,4)
 274 --   >>> v1
 275 --   (1,2,3,4)
 276 --   >>> let Vn (Vec4D w x y z) = v1
 277 --   >>> (w,x,y,z)
 278 --   (1,2,3,4)
 279 --
 280 make4d :: forall a. (a,a,a,a) -> Vn Vec4D a
 281 make4d (w,x,y,z) = Vn (Vec4D w x y z)