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.Vector.Fixed as V
  11 import Data.Vector.Fixed.Boxed
  12 import Data.Vector.Fixed.Internal
  13
  14 import Normed
  15
  16 -- | The Vn newtype simply wraps (Vector v a) so that we avoid
  17 --   undecidable instances.
  18 newtype Vn a = Vn a
  19   deriving (Show)
  20
  21 -- | We would really like to say, "anything that is a vector of
  22 --   equatable things is itself equatable." The 'Vn' class
  23 --   allows us to express this without a GHC battle.
  24 --
  25 --   Examples:
  26 --
  27 --   >>> let v1 = make2d (1,2)
  28 --   >>> let v2 = make2d (1,2)
  29 --   >>> let v3 = make2d (3,4)
  30 --   >>> v1 == v2
  31 --   True
  32 --   >>> v1 == v3
  33 --   False
  34 --
  35 instance (Eq a, Vector v a, Vector v Bool) => Eq (Vn (v a)) where
  36   (Vn v1) == (Vn v2) = V.foldl (&&) True (V.zipWith (==) v1 v2)
  37
  38 -- | The use of 'Num' here is of course incorrect (otherwise, we
  39 --   wouldn't have to throw errors). But it's really nice to be able
  40 --   to use normal addition/subtraction.
  41 instance (Num a, Vector v a) => Num (Vn (v a)) where
  42   -- | Componentwise addition.
  43   --
  44   --   Examples:
  45   --
  46   --   >>> let v1 = make2d (1,2)
  47   --   >>> let v2 = make2d (3,4)
  48   --   >>> v1 + v2
  49   --   Vn fromList [4,6]
  50   --
  51   (Vn v1) + (Vn v2) = Vn $ V.zipWith (+) v1 v2
  52
  53   (Vn v1) - (Vn v2) = Vn $ V.zipWith (-) v1 v2
  54   fromInteger x = Vn $ V.replicate (fromInteger x)
  55   (*) = error "multiplication of vectors is undefined"
  56   abs = error "absolute value of vectors is undefined"
  57   signum = error "signum of vectors is undefined"
  58
  59 instance Functor Vn where
  60   fmap f (Vn v1) = Vn (f v1)
  61
  62 instance (RealFloat a, Ord a, Vector v a) => Normed (Vn (v a)) where
  63   -- We don't use V.maximum here because it relies on a type
  64   -- constraint that the vector be non-empty and I don't know how to
  65   -- pattern match it away.
  66   norm_infty (Vn v1) = fromRational $ toRational $ V.foldl max 0 v1
  67
  68   norm_p p (Vn v1) =
  69     fromRational $ toRational $ root $ V.sum $ V.map (exponentiate . abs) v1
  70     where
  71       exponentiate = (** (fromIntegral p))
  72       root = (** (recip (fromIntegral p)))
  73
  74 -- | Dot (standard inner) product.
  75 dot :: (Num a, Vector v a) => Vn (v a) -> Vn (v a) -> a
  76 dot (Vn v1) (Vn v2) = V.sum $ V.zipWith (*) v1 v2
  77
  78 -- | The angle between @v1@ and @v2@ in Euclidean space.
  79 angle :: (RealFloat a, Vector v a) => Vn (v a) -> Vn (v a) -> a
  80 angle v1 v2 =
  81   acos theta
  82   where
  83     theta = (v1 `dot` v2) / norms
  84     norms = (norm_p 2 v1) * (norm_p 2 v2)
  85
  86 -- | Convenience function for 2d vectors.
  87 make2d :: forall a. (a,a) -> Vn (Vec2 a)
  88 make2d (x,y) =
  89   Vn v1
  90   where
  91     v1 = vec $ con |> x |> y :: Vec2 a
  92
  93 -- | Convenience function for 3d vectors.
  94 make3d :: forall a. (a,a,a) -> Vn (Vec3 a)
  95 make3d (x,y,z) =
  96   Vn v1
  97   where
  98     v1 = vec $ con |> x |> y |> z :: Vec3 a