From: Michael Orlitzky <michael@orlitzky.com>
Date: Sun, 3 Feb 2013 17:16:43 +0000 (-0500)
Subject: Begin implementation of normed, fixed-length vectors.
X-Git-Url: http://gitweb.michael.orlitzky.com/?p=numerical-analysis.git;a=commitdiff_plain;h=74e79199dcfec0639133ae9990dc33a2c5a095f0

Begin implementation of normed, fixed-length vectors.
---

diff --git a/src/FixedVector.hs b/src/FixedVector.hs
new file mode 100644
index 0000000..c01ae60
--- /dev/null
+++ b/src/FixedVector.hs
@@ -0,0 +1,98 @@
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeFamilies #-}
+
+module FixedVector
+where
+
+import Data.Vector.Fixed as V
+import Data.Vector.Fixed.Boxed
+import Data.Vector.Fixed.Internal
+
+import Normed
+
+-- | The Vn newtype simply wraps (Vector v a) so that we avoid
+--   undecidable instances.
+newtype Vn a = Vn a
+  deriving (Show)
+
+-- | 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, Vector v Bool) => Eq (Vn (v a)) where
+  (Vn v1) == (Vn v2) = V.foldl (&&) True (V.zipWith (==) v1 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
+  --   Vn fromList [4,6]
+  --
+  (Vn v1) + (Vn v2) = Vn $ V.zipWith (+) v1 v2
+
+  (Vn v1) - (Vn v2) = Vn $ V.zipWith (-) v1 v2
+  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"
+
+instance Functor Vn where
+  fmap f (Vn v1) = Vn (f v1)
+
+instance (RealFloat a, Ord a, Vector v a) => Normed (Vn (v a)) where
+  -- 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.
+  norm_infty (Vn v1) = fromRational $ toRational $ V.foldl max 0 v1
+
+  norm_p p (Vn v1) =
+    fromRational $ toRational $ root $ V.sum $ V.map (exponentiate . abs) v1
+    where
+      exponentiate = (** (fromIntegral p))
+      root = (** (recip (fromIntegral p)))
+
+-- | Dot (standard inner) product.
+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.
+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_p 2 v1) * (norm_p 2 v2)
+
+-- | Convenience function for 2d vectors.
+make2d :: forall a. (a,a) -> Vn (Vec2 a)
+make2d (x,y) =
+  Vn v1
+  where
+    v1 = vec $ con |> x |> y :: Vec2 a
+
+-- | Convenience function for 3d vectors.
+make3d :: forall a. (a,a,a) -> Vn (Vec3 a)
+make3d (x,y,z) =
+  Vn v1
+  where
+    v1 = vec $ con |> x |> y |> z :: Vec3 a
diff --git a/src/Normed.hs b/src/Normed.hs
new file mode 100644
index 0000000..b60c2b1
--- /dev/null
+++ b/src/Normed.hs
@@ -0,0 +1,29 @@
+{-# LANGUAGE FlexibleInstances #-}
+
+-- | The 'Normed' class represents elements of a normed vector
+--   space. We define instances for all common numeric types.
+module Normed
+where
+
+import Data.Number.BigFloat
+
+class Normed a where
+  norm_p :: (Integral c, RealFrac b) => c -> a -> b
+  norm_infty :: RealFrac b => a -> b
+
+-- Define instances for common numeric types.
+instance Normed Integer where
+  norm_p _ = fromInteger
+  norm_infty = fromInteger
+
+instance Normed Rational where
+  norm_p _ = fromRational
+  norm_infty = fromRational
+
+instance Epsilon e => Normed (BigFloat e) where
+  norm_p _ = fromRational . toRational
+  norm_infty = fromRational . toRational
+
+instance Normed Double where
+  norm_p _ = fromRational . toRational
+  norm_infty = fromRational . toRational