packages feed

vector-space 0.1.3 → 0.2.0

raw patch · 9 files changed

+1346/−231 lines, 9 filesdep +OpenGLdep +old-timenew-component:exe:Perf

Dependencies added: OpenGL, old-time

Files

+ src/Data/AdditiveGroup.hs view
@@ -0,0 +1,71 @@+----------------------------------------------------------------------+-- |+-- Module      :   Data.AdditiveGroup+-- Copyright   :  (c) Conal Elliott and Andy J Gill 2008+-- License     :  BSD3+-- +-- Maintainer  :  conal@conal.net, andygill@ku.edu+-- Stability   :  experimental+-- +-- Groups: zero, addition, and negation (additive inverse)+----------------------------------------------------------------------++module Data.AdditiveGroup+  ( +    AdditiveGroup(..), (^-^)+  ) where++import Control.Applicative+import Data.Complex hiding (magnitude)++infixl 6 ^+^, ^-^++-- | Additive group @v@.+class AdditiveGroup v where+  -- | The zero element: identity for '(^+^)'+  zeroV :: v+  -- | Add vectors+  (^+^) :: v -> v -> v+  -- | Additive inverse+  negateV :: v -> v++-- | Group subtraction+(^-^) :: AdditiveGroup v => v -> v -> v+v ^-^ v' = v ^+^ negateV v'++instance AdditiveGroup Double where+  zeroV   = 0.0+  (^+^)   = (+)+  negateV = negate++instance AdditiveGroup Float where+  zeroV   = 0.0+  (^+^)   = (+)+  negateV = negate++instance (RealFloat v, AdditiveGroup v) => AdditiveGroup (Complex v) where+  zeroV   = zeroV :+ zeroV+  (^+^)   = (+)+  negateV = negate++-- Hm.  The 'RealFloat' constraint is unfortunate here.  It's due to a+-- questionable decision to place 'RealFloat' into the definition of the+-- 'Complex' /type/, rather than in functions and instances as needed.++instance (AdditiveGroup u,AdditiveGroup v) => AdditiveGroup (u,v) where+  zeroV             = (zeroV,zeroV)+  (u,v) ^+^ (u',v') = (u^+^u',v^+^v')+  negateV (u,v)     = (negateV u, negateV v)++instance (AdditiveGroup u,AdditiveGroup v,AdditiveGroup w)+    => AdditiveGroup (u,v,w) where+  zeroV                  = (zeroV,zeroV,zeroV)+  (u,v,w) ^+^ (u',v',w') = (u^+^u',v^+^v',w^+^w')+  negateV (u,v,w)        = (negateV u, negateV v, negateV w)+++-- Standard instance for an applicative functor applied to a vector space.+instance AdditiveGroup v => AdditiveGroup (a->v) where+  zeroV   = pure   zeroV+  (^+^)   = liftA2 (^+^)+  negateV = fmap   negateV
+ src/Data/Cross.hs view
@@ -0,0 +1,110 @@+{-# LANGUAGE FlexibleInstances, TypeOperators, UndecidableInstances+           , TypeSynonymInstances+  #-}+{-# OPTIONS_GHC -Wall #-}+----------------------------------------------------------------------+-- |+-- Module      :  Data.Cross+-- Copyright   :  (c) Conal Elliott 2008+-- License     :  BSD3+-- +-- Maintainer  :  conal@conal.net+-- Stability   :  experimental+-- +-- Cross products and normals+----------------------------------------------------------------------++module Data.Cross+  (+    HasNormal(..), normal+  , One, Two, Three+  , HasCross2(..), HasCross3(..)+  ) where++import Data.VectorSpace+import Data.LinearMap+import Data.Derivative++-- | Thing with a normal vector (not necessarily normalized).+class HasNormal v where normalVec :: v -> v++-- | Normalized normal vector.  See also 'cross'.+normal :: (HasNormal v, InnerSpace v s, Floating s) => v -> v+normal = normalized . normalVec+++-- | Singleton+type One   s = s++-- | Homogeneous pair+type Two   s = (s,s)++-- | Homogeneous triple+type Three s = (s,s,s)++-- | Cross product of various forms of 2D vectors+class HasCross2 v where cross2 :: v -> v++instance Num s => HasCross2 (s,s) where+  cross2 (x,y) = (-y,x)  -- or @(y,-x)@?++-- TODO: Eliminate the 'Num' constraint by using negateV.++-- "Variable occurs more often in a constraint than in the instance+-- head".  Hence UndecidableInstances.++instance (LMapDom a s, VectorSpace v s, HasCross2 v) => HasCross2 (a:>v) where+  -- 2d cross-product is linear+  cross2 = fmapD cross2++instance (Num s, LMapDom s s) => HasNormal (One s :> Two s) where+  normalVec v = cross2 (derivativeAt v 1)++-- Does this problem come from the choice of 'VectorSpace' instance?++instance (Num s, LMapDom s s, VectorSpace s s)+    => HasNormal (Two (One s :> s)) where+  normalVec = unpairD . normalVec . pairD+++-- | Cross product of various forms of 3D vectors+class HasCross3 v where cross3 :: v -> v -> v++instance Num s => HasCross3 (s,s,s) where+  (ax,ay,az) `cross3` (bx,by,bz) = ( ay * bz - az * by+                                   , az * bx - ax * bz+                                   , ax * by - ay * bx )++-- TODO: Eliminate the 'Num' constraint by using 'VectorSpace' operations.++instance (LMapDom a s, VectorSpace v s, HasCross3 v) => HasCross3 (a:>v) where+  -- 3D cross-product is bilinear (curried linear)+  cross3 = distrib cross3++instance (Num s, LMapDom s s) => HasNormal (Two s :> Three s) where+  normalVec v = d (1,0) `cross3` d (0,1)+   where+     d = derivativeAt v++instance (Num s, VectorSpace s s, LMapDom s s) => HasNormal (Three (Two s :> s)) where+  normalVec = untripleD . normalVec . tripleD++---- Could go elsewhere++pairD :: (LMapDom a s, VectorSpace b s, VectorSpace c s) =>+         (a:>b,a:>c) -> a:>(b,c)+pairD (u,v) = liftD2 (,) u v++tripleD :: (LMapDom a s, VectorSpace b s, VectorSpace c s, VectorSpace d s) =>+           (a:>b,a:>c,a:>d) -> a:>(b,c,d)+tripleD (u,v,w) = liftD3 (,,) u v w++unpairD :: (LMapDom a s, VectorSpace a s, VectorSpace b s, VectorSpace c s) =>+           (a :> (b,c)) -> (a:>b, a:>c)+unpairD d = (fst <$>> d, snd <$>> d)++untripleD :: ( LMapDom a s , VectorSpace a s, VectorSpace b s+             , VectorSpace c s, VectorSpace d s) =>+             (a :> (b,c,d)) -> (a:>b, a:>c, a:>d)+untripleD d =+  ((\ (a,_,_) -> a) <$>> d, (\ (_,b,_) -> b) <$>> d, (\ (_,_,c) -> c) <$>> d)
src/Data/Derivative.hs view
@@ -1,5 +1,12 @@-{-# LANGUAGE TypeOperators, MultiParamTypeClasses, UndecidableInstances #-}-{-# OPTIONS_GHC -Wall #-}+{-# LANGUAGE TypeOperators, MultiParamTypeClasses, UndecidableInstances+           , TypeSynonymInstances, FlexibleInstances, FunctionalDependencies+           , FlexibleContexts+           , GeneralizedNewtypeDeriving, StandaloneDeriving+  #-}++-- TODO: remove FlexibleContexts++-- {-# OPTIONS_GHC -Wall #-} ---------------------------------------------------------------------- -- | -- Module      :  Data.Derivative@@ -9,171 +16,178 @@ -- Maintainer  :  conal@conal.net -- Stability   :  experimental -- --- Infinite derivative towers via linear maps.  See blog posts--- <http://conal.net/blog/tag/derivatives/>+-- This module is a wrapper around Data.Maclaurin or Data.Horner, to+-- change the 'VectorSpace' instance for '(:>)'. ----------------------------------------------------------------------  module Data.Derivative   (-    (:>)(..), (:~>), dZero, dConst+    (:>), powVal, derivative, derivativeAt+  , (:~>), dZero, pureD+  , fmapD, (<$>>){-, (<*>>)-}, liftD2, liftD3   , idD, fstD, sndD-  , linearD, distribD+  , linearD, distrib   , (@.), (>-<)   ) where -import Control.Applicative+import Data.Function (on)  import Data.VectorSpace import Data.NumInstances ()+import Data.LinearMap +import qualified Data.Maclaurin as D+-- import qualified Data.Horner as D -infixr 9 `D`, @.+-- TODO: sync Data.Horner interface with Data.Maclaurin.  Better: refactor+-- Maclaurin into a module that depends on the representation (e.g.,+-- Horner vs Maclaurin) and the rest that doesn't.++infixr 9 @.+infixl 4 {-<*>>,-} <$>> infix  0 >-<   -- | Tower of derivatives.--- --- Warning, the 'Applicative' instance is missing its 'pure' (due to a--- 'VectorSpace' type constraint).  Use 'dConst' instead.-data a :> b = D { dVal :: b, dDeriv :: a :-* (a :> b) }+newtype (a :> b) = T { unT :: a D.:> b } --- data a :> b = D b (a :-* (a :> b))+deriving instance (VectorSpace b b, LMapDom a b, Num b       ) => Num           (a :> b)+deriving instance (VectorSpace b b, LMapDom a b, Fractional b) => Fractional    (a :> b)+deriving instance (VectorSpace b b, LMapDom a b, Floating b  ) => Floating      (a :> b)+deriving instance (VectorSpace b s, LMapDom a s              ) => AdditiveGroup (a :> b) ++inT :: ((a D.:> b) -> (c D.:> d))+    -> ((a   :> b) -> (c   :> d))+inT  = (T .).(. unT)++inT2 :: ((a D.:> b) -> (c D.:> d) -> (e D.:> f))+     -> ((a   :> b) -> (c   :> d) -> (e   :> f))+inT2  = (inT .).(. unT)++inT3 :: ((a D.:> b) -> (c D.:> d) -> (e D.:> f) -> (g D.:> h))+     -> ((a   :> b) -> (c   :> d) -> (e   :> f) -> (g   :> h))+inT3  = (inT2 .).(. unT)++-- | Extract the value from a derivative tower+powVal :: (a :> b) -> b+powVal = D.powVal . unT++-- | Extract the derivative from a derivative tower+derivative :: (VectorSpace b s, LMapDom a s) =>+              (a :> b) -> (a :-* (a :> b))+derivative = fmapL T . D.derivative . unT++-- | Sampled derivative.  For avoiding an awkward typing problem related+-- to the two required 'VectorSpace' instances.+derivativeAt :: (LMapDom a s, VectorSpace b s) =>+                (a :> b) -> a -> (a :> b)+derivativeAt d = T . D.derivativeAt (unT d)++-- The definition of 'derivativeAt' takes care to share partial+-- applications of 'D.derivativeAt', which is useful in power series+-- representations for which 'derivative' is not free (Horner).+ -- | Infinitely differentiable functions type a :~> b = a -> (a:>b) --- So we could define--- ---   data a :> b = D b (a :~> b)--- --- with the restriction that the a :~> b is linear+-- -- Handy for missing methods.+-- noOv :: String -> a+-- noOv op = error (op ++ ": not defined on a :> b") -instance Functor ((:>) a) where-  fmap f (D b b') = D (f b) ((fmap.fmap) f b') --- I think fmap will be meaningful only with *linear* functions.+-- | Derivative tower full of 'zeroV'.+dZero :: (LMapDom a s, VectorSpace b s) => a:>b+dZero = T D.dZero --- Handy for missing methods.-noOv :: String -> a-noOv op = error (op ++ ": not defined on a :> b") -instance Applicative ((:>) a) where-    -- pure = dConst    -- not!  see below.-    pure = noOv "pure.  use dConst instead."-    D f f' <*> D b b' = D (f b) (liftA2 (<*>) f' b')+-- | Constant derivative tower.+pureD :: (LMapDom a s, VectorSpace b s) => b -> a:>b+pureD = fmap T D.pureD --- Why can't we define 'pure' as 'dConst'?  Because of the extra type--- constraint that @VectorSpace b@ (not @a@).  Oh well.  Be careful not to--- use 'pure', okay?  Alternatively, I could define the '(<*>)' (naming it--- something else) and then say @foo <$> p <*^> q <*^> ...@.+-- | Map a /linear/ function over a derivative tower.+fmapD, (<$>>) :: (LMapDom a s, VectorSpace b s) =>+                 (b -> c) -> (a :> b) -> (a :> c)+fmapD = fmap inT D.fmapD --- | Constant derivative tower.-dConst :: VectorSpace b s => b -> a:>b-dConst b = b `D` const dZero+(<$>>) = fmapD --- | Derivative tower full of 'zeroV'.-dZero :: VectorSpace b s => a:>b-dZero = dConst zeroV+-- -- | Like '(<*>)' for derivative towers.+-- (<*>>) :: (LMapDom a s, VectorSpace b s, VectorSpace c s) =>+--           (a :> (b -> c)) -> (a :> b) -> (a :> c)+-- (<*>>) = inT2 (D.<*>>) +-- | Apply a /linear/ binary function over derivative towers.+liftD2 :: (VectorSpace b s, LMapDom a s, VectorSpace c s, VectorSpace d s) =>+          (b -> c -> d) -> (a :> b) -> (a :> c) -> (a :> d)+liftD2 = fmap inT2 D.liftD2+++-- | Apply a /linear/ ternary function over derivative towers.+liftD3 :: ( LMapDom a s+          , VectorSpace b s, VectorSpace c s+          , VectorSpace d s, VectorSpace e s ) =>+          (b -> c -> d -> e)+       -> (a :> b) -> (a :> c) -> (a :> d) -> (a :> e)+liftD3 = fmap inT3 D.liftD3+ -- | Differentiable identity function.  Sometimes called "the -- derivation variable" or similar, but it's not really a variable.-idD :: VectorSpace u s => u :~> u-idD = linearD id---- or---   dId v = D v dConst+idD :: (LMapDom u s, VectorSpace u s) => u :~> u+idD = fmap T D.idD  -- | Every linear function has a constant derivative equal to the function -- itself (as a linear map).-linearD :: VectorSpace v s => (u :-* v) -> (u :~> v)-linearD f u = D (f u) (dConst . f)-+linearD :: (LMapDom u s, VectorSpace v s) => (u -> v) -> (u :~> v)+linearD = (fmap.fmap) T D.linearD  -- Other examples of linear functions  -- | Differentiable version of 'fst'-fstD :: VectorSpace a s => (a,b) :~> a-fstD = linearD fst+fstD :: (VectorSpace a s, LMapDom b s, LMapDom a s) => (a,b) :~> a+fstD = fmap T D.fstD  -- | Differentiable version of 'snd'-sndD :: VectorSpace b s => (a,b) :~> b-sndD = linearD snd+sndD :: (VectorSpace b s, LMapDom b s, LMapDom a s) => (a,b) :~> b+sndD = fmap T D.sndD  -- | Derivative tower for applying a binary function that distributes over--- addition, such as multiplication.-distribD :: (VectorSpace u s) =>-             (b -> c -> u) -> ((a :> b) -> (a :> c) -> (a :> u))-          -> (a :> b) -> (a :> c) -> (a :> u)-distribD op opD u@(D u0 u') v@(D v0 v') =-  D (u0 `op` v0) ((u `opD`) . v' ^+^ (`opD` v) . u')+-- addition, such as multiplication.  A bit weaker assumption than+-- bilinearity.+distrib :: (LMapDom a s, VectorSpace b s, VectorSpace c s, VectorSpace u s) =>+           (b -> c -> u) -> (a :> b) -> (a :> c) -> (a :> u)+distrib = fmap inT2 D.distrib --- Equivalently,--- ---   distribD op opD u@(D u0 u') v@(D v0 v') =---     D (u0 `op` v0) (\ da -> (u `opD` v' da) ^+^ (u' da `opD` v))+-- I'm not sure about the next three, which discard information +instance Show b => Show (a :> b) where show    = show     .   unT+instance Eq   b => Eq   (a :> b) where (==)    = (==)    `on` unT+instance Ord  b => Ord  (a :> b) where compare = compare `on` unT --- I'm not sure about the next three, which discard information-instance Show b => Show (a :> b) where show    = noOv "show"-instance Eq   b => Eq   (a :> b) where (==)    = noOv "(==)"-instance Ord  b => Ord  (a :> b) where compare = noOv "compare"+-- These next two instances are the reason for having this module.  The+-- scalar field differs from the one used in the underlying+-- representation.  I can't have both instances.  First, there is a+-- functional dependency on 'VectorSpace', which says that a vector space+-- type determines its scalar field type.  I experimented with dropping+-- the fundep, and then managing the resulting ambiguity.  However,+-- the two 'VectorSpace' instances overlap considerably. -instance VectorSpace u s => VectorSpace (a :> u) (a :> s) where-  zeroV   = dConst   zeroV    -- or dZero-  (*^)    = distribD (*^) (*^)-  negateV = fmap     negateV-  (^+^)   = liftA2   (^+^)+instance (LMapDom a s, VectorSpace u s, VectorSpace s s)+    => VectorSpace (a :> u) (a :> s) where+  (*^)    = inT2 (D.**^) -- not '(D.*^)' -instance (InnerSpace u s, InnerSpace s s') =>+instance (InnerSpace u s, InnerSpace s s', VectorSpace s s, LMapDom a s) =>      InnerSpace (a :> u) (a :> s) where-  (<.>) = distribD (<.>) (<.>)+  (<.>)    = inT2 (D.<*.>) -- not '(D.<.>)' --- | Chain rule.-(@.) :: (b :~> c) -> (a :~> b) -> (a :~> c)-(h @. g) a0 = D c0 (c' @. b')-  where-    D b0 b' = g a0-    D c0 c' = h b0 +-- | Chain rule.+(@.) :: (LMapDom b s, LMapDom a s, VectorSpace c s) =>+        (b :~> c) -> (a :~> b) -> (a :~> c)+h @. g = T . ((unT . h) D.@. (unT . g))  -- | Specialized chain rule.-(>-<) :: VectorSpace u s => (u -> u) -> ((a :> u) -> (a :> s))+(>-<) :: (LMapDom a s, VectorSpace s s, VectorSpace u s) =>+         (u -> u) -> ((a :> u) -> (a :> s))       -> (a :> u) -> (a :> u)--f >-< f' = \ u@(D u0 u') -> D (f u0) ((f' u *^) . u')---- Equivalently:--- ---   f >-< f' = \ u@(D u0 u') -> D (f u0) (\ da -> f' u *^ u' da)--instance (Num b, VectorSpace b b) => Num (a:>b) where-  fromInteger = dConst . fromInteger-  (+) = liftA2   (+)-  (-) = liftA2   (-)-  (*) = distribD (*) (*)-  negate = negate >-< -1-  abs    = abs    >-< signum-  signum = signum >-< 0  -- derivative wrong at zero--instance (Fractional b, VectorSpace b b) => Fractional (a:>b) where-  fromRational = dConst . fromRational-  recip        = recip >-< recip sqr--sqr :: Num a => a -> a-sqr x = x*x--instance (Floating b, VectorSpace b b) => Floating (a:>b) where-  pi    = dConst pi-  exp   = exp   >-< exp-  log   = log   >-< recip-  sqrt  = sqrt  >-< recip (2 * sqrt)-  sin   = sin   >-< cos-  cos   = cos   >-< - sin-  sinh  = sinh  >-< cosh-  cosh  = cosh  >-< sinh-  asin  = asin  >-< recip (sqrt (1-sqr))-  acos  = acos  >-< recip (- sqrt (1-sqr))-  atan  = atan  >-< recip (1+sqr)-  asinh = asinh >-< recip (sqrt (1+sqr))-  acosh = acosh >-< recip (- sqrt (sqr-1))-  atanh = atanh >-< recip (1-sqr)+f >-< f' = inT (f D.>-< (unT . f' . T))
+ src/Data/Horner.hs view
@@ -0,0 +1,220 @@+{-# LANGUAGE TypeOperators, MultiParamTypeClasses, UndecidableInstances+           , TypeSynonymInstances, FlexibleInstances, FunctionalDependencies+  #-}+{-# OPTIONS_GHC -Wall #-}+----------------------------------------------------------------------+-- |+-- Module      :  Data.Horner+-- Copyright   :  (c) Conal Elliott 2008+-- License     :  BSD3+-- +-- Maintainer  :  conal@conal.net+-- Stability   :  experimental+-- +-- Infinite derivative towers via linear maps, using the Horner+-- representation.  See blog posts <http://conal.net/blog/tag/derivatives/>.+----------------------------------------------------------------------++module Data.Horner+  (+    (:>), powVal, derivative, integral+  , (:~>), dZero, dConst+  , idD, fstD, sndD+  , linearD, distrib+  , (@.), (>-<)+  -- , HasDeriv(..)+  )+  where++import Control.Applicative++import Data.VectorSpace+import Data.LinearMap+import Data.NumInstances ()++infixr 9 `H`, @.+infix  0 >-<++-- | Power series+-- +-- Warning, the 'Applicative' instance is missing its 'pure' (due to a+-- 'VectorSpace' type constraint).  Use 'dConst' instead.+data a :> b = H b (a :-* (a :> b))++-- | The plain-old (0th order) value+powVal :: (a :> b) -> b+powVal (H b _) = b++-- Apply successive functions to successive values+apPow :: [b -> c] -> (a :> b) -> (a :> c)+apPow [] _ = error "apPow: finite function list"+apPow (f : fs) (b0 `H` bt) = H (f b0) (apPow fs . bt)++-- Count.  Avoids the 'Enum' requirement of [1..]+from :: Num s => s -> [s]+from n = n : from (n+1) ++-- | Derivative of a power series+derivative :: (VectorSpace b s, Num s) =>+         (a :> b) -> (a :-* (a :> b))+derivative (H _ bt) = apPow ((*^) <$> from 1) . bt++-- | Integral of a power series+integral :: (VectorSpace b s, Fractional s) =>+            b -> (a :-* (a :> b)) -> (a :> b)+integral b0 bt = H b0 (apPow (((*^).recip) <$> from 1) . bt)++-- | Infinitely differentiable functions+type a :~> b = a -> (a:>b)++-- So we could define+-- +--   data a :> b = H b (a :~> b)+-- +-- with the restriction that the a :~> b is linear++instance Functor ((:>) a) where+  fmap f (H b b') = H (f b) ((fmap.fmap) f b')++-- I think fmap will be meaningful only with *linear* functions.++-- Handy for missing methods.+noOv :: String -> a+noOv op = error (op ++ ": not defined on a :> b")++instance Applicative ((:>) a) where+  -- pure = dConst    -- not!  see below.+  pure = noOv "pure"  -- use dConst instead+  H f f' <*> H b b' = H (f b) (liftA2 (<*>) f' b')++-- Why can't we define 'pure' as 'dConst'?  Because of the extra type+-- constraint that @VectorSpace b@ (not @a@).  Oh well.  Be careful not to+-- use 'pure', okay?  Alternatively, I could define the '(<*>)' (naming it+-- something else) and then say @foo <$> p <*^> q <*^> ...@.++-- | Constant derivative tower.+dConst :: VectorSpace b s => b -> a:>b+dConst b = b `H` const dZero++-- | Derivative tower full of 'zeroV'.+dZero :: VectorSpace b s => a:>b+dZero = dConst zeroV++-- | Differentiable identity function.  Sometimes called "the+-- derivation variable" or similar, but it's not really a variable.+idD :: VectorSpace u s => u :~> u+idD = linearD id++-- or+--   dId v = H v dConst++-- | Every linear function has a constant derivative equal to the function+-- itself (as a linear map).+linearD :: VectorSpace v s => (u :-* v) -> (u :~> v)+linearD f u = H (f u) (dConst . f)+++-- Other examples of linear functions++-- | Differentiable version of 'fst'+fstD :: VectorSpace a s => (a,b) :~> a+fstD = linearD fst++-- | Differentiable version of 'snd'+sndD :: VectorSpace b s => (a,b) :~> b+sndD = linearD snd++-- | Derivative tower for applying a binary function that distributes over+-- addition, such as multiplication.  A bit weaker assumption than+-- bilinearity.+distrib :: (VectorSpace u s) =>+           (b -> c -> u) -> (a :> b) -> (a :> c) -> (a :> u)+distrib op = opD+ where+   opD (H u0 ut) v@(H v0 vt) =+     H (u0 `op` v0) (fmap (u0 `op`) . vt ^+^ (`opD` v) . ut)+++-- Equivalently,+-- +--   distrib op = opD+--    where+--      opD u@(H u0 u') v@(H v0 v') =+--        H (u0 `op` v0) (\ da -> ((u0 `op`) <$> v' da) ^+^ (u' da `opD` v))++++-- I'm not sure about the next three, which discard information++instance Show b => Show (a :> b) where show    = noOv "show"+instance Eq   b => Eq   (a :> b) where (==)    = noOv "(==)"+instance Ord  b => Ord  (a :> b) where compare = noOv "compare"++instance (LMapDom a s, VectorSpace u s) => AdditiveGroup (a :> u) where+  zeroV   = pureD  zeroV    -- or dZero+  negateV = fmapD  negateV+  (^+^)   = liftD2 (^+^)++instance (LMapDom a s, VectorSpace u s) => VectorSpace (a :> u) s where+  (*^) s = fmapD  ((*^) s)++(**^) :: (VectorSpace c s, VectorSpace s s, LMapDom a s) =>+         (a :> s) -> (a :> c) -> (a :> c)+(**^) = distrib (*^)++-- | Chain rule.+(@.) :: (VectorSpace b s, VectorSpace c s, Num s) =>+        (b :~> c) -> (a :~> b) -> (a :~> c)+(h @. g) a0 = H c0 (derivative c @. derivative b)+  where+    b@(H b0 _) = g a0+    c@(H c0 _) = h b0+++-- | Specialized chain rule.+(>-<) :: (VectorSpace u s, Fractional s) => (u -> u) -> ((a :> u) -> (a :> s))+      -> (a :> u) -> (a :> u)++-- f >-< f' = \ u@(D u0 u') -> D (f u0) ((f' u *^) . u')++f >-< f' = \ u@(H u0 _) -> integral (f u0) ((f' u *^) . derivative u)++-- TODO: consider eliminating @Num s@.  I just need a multiplicative unit.++-- Equivalently:+-- +--   f >-< f' = \ u@(H u0 u') -> H (f u0) (\ da -> f' u *^ u' da)++instance (Fractional b, VectorSpace b b) => Num (a:>b) where+  fromInteger = dConst . fromInteger+  (+) = liftA2  (+)+  (-) = liftA2  (-)+  (*) = distrib (*)+  +  negate = negate >-< -1+  abs    = abs    >-< signum+  signum = signum >-< 0  -- derivative wrong at zero++instance (Fractional b, VectorSpace b b) => Fractional (a:>b) where+  fromRational = dConst . fromRational+  recip        = recip >-< recip sqr++sqr :: Num a => a -> a+sqr x = x*x++instance (Floating b, VectorSpace b b) => Floating (a:>b) where+  pi    = dConst pi+  exp   = exp   >-< exp+  log   = log   >-< recip+  sqrt  = sqrt  >-< recip (2 * sqrt)+  sin   = sin   >-< cos+  cos   = cos   >-< - sin+  sinh  = sinh  >-< cosh+  cosh  = cosh  >-< sinh+  asin  = asin  >-< recip (sqrt (1-sqr))+  acos  = acos  >-< recip (- sqrt (1-sqr))+  atan  = atan  >-< recip (1+sqr)+  asinh = asinh >-< recip (sqrt (1+sqr))+  acosh = acosh >-< recip (- sqrt (sqr-1))+  atanh = atanh >-< recip (1-sqr)+
+ src/Data/LinearMap.hs view
@@ -0,0 +1,315 @@+{-# LANGUAGE TypeOperators, TypeFamilies, UndecidableInstances+           , MultiParamTypeClasses, FlexibleInstances+           , FunctionalDependencies+  #-}+{-# OPTIONS_GHC -Wall -fno-warn-orphans -fglasgow-exts #-}+----------------------------------------------------------------------+-- |+-- Module      :  Data.LinearMap+-- Copyright   :  (c) Conal Elliott 2008+-- License     :  BSD3+-- +-- Maintainer  :  conal@conal.net+-- Stability   :  experimental+-- +-- Linear maps+----------------------------------------------------------------------++module Data.LinearMap+  (+    LMapDom(..), inL, inL2, inL3+  , linearK, (.*)+  , pureL --, (<*>*)+  , fmapL, (<$>*), liftL2, liftL3, idL+  ) where++-- -fglasgow-exts above enables the RULES pragma++import Control.Applicative+import Data.Function++import Data.VectorSpace++-- Temporary+import Graphics.Rendering.OpenGL.GL.CoordTrans+  (Vector2(..),Vector3(..))+++infixr 9 :-*+infixr 9 .*+infixl 9 `lapply`+infixl 4 <$>* -- , <*>*+++-- | Domain of a linear map.+class VectorSpace a s => LMapDom a s | a -> s where+  -- | Linear map type+  data (:-*) a :: * -> *+  -- | Linear map as function+  lapply :: VectorSpace b s => (a :-* b) -> (a -> b)+  -- | Function (assumed linear) as linear map.+  linear :: (a -> b) -> (a :-* b)++-- Neither 'VectorSpace' nor even 'AdditiveGroup' is really required as a+-- 'LMapDom' superclass.  Instead, we could have additional constraints in+-- some 'LMapDom' instances and related functions.+++{-# RULES+"linear.lapply"   forall m. linear (lapply m) = m+"lapply.linear"   forall f. lapply (linear f) = f+ #-}+++-- | Transform a linear map by transforming a linear function.+inL :: (LMapDom c s, VectorSpace b s', LMapDom a s') =>+        ((a -> b) -> (c -> d)) -> ((a :-* b) -> (c :-* d))+{-# INLINE inL #-}+inL h = linear . h . lapply++-- | Transform a linear maps by transforming linear functions.+inL2 :: ( LMapDom c s, VectorSpace b s', LMapDom a s'+        , LMapDom e s, VectorSpace d s ) =>+        ((a -> b) -> (c -> d) -> (e -> f))+     -> ((a :-* b) -> (c :-* d) -> (e :-* f))+{-# INLINE inL2 #-}+inL2 h = inL . h . lapply++-- inL2 h m n+--   = (inL . h . lapply) m n+--   = inL (h (lapply m)) n+--   = (linear . (h (lapply m)) . lapply) n+--   = linear (h (lapply m) (lapply n)+--   = linear (h (lapply m) (lapply n))+++-- | Transform a linear maps by transforming linear functions.+inL3 :: ( LMapDom a s, VectorSpace b s+        , VectorSpace f s , LMapDom p s, LMapDom c s'+        , VectorSpace d s', LMapDom e s ) =>+        ((a -> b) -> (c -> d) -> (e -> f) -> (p -> q))+     -> ((a :-* b) -> (c :-* d) -> (e :-* f) -> (p :-* q))+{-# INLINE inL3 #-}+inL3 h = inL2 . h . lapply+++-- TODO: go through this whole module and relax constraints on scalar+-- fields.  See if it helps with the scalar dilemma in Mac2+++-- | Constant value as a linear map+pureL :: LMapDom a s => b -> (a :-* b)+pureL b = linear (const b)++-- -- | Like '(<*>)' for linear maps.+-- (<*>*) :: (LMapDom a s, VectorSpace b s, VectorSpace c s) =>+--         (a :-* (b -> c)) -> (a :-* b) -> (a :-* c)+-- (<*>*) = inL2 (<*>)++-- | Map a /linear/ function over a linear map.+fmapL, (<$>*) :: (LMapDom a s, VectorSpace b s) =>+                 (b -> c) -> (a :-* b) -> (a :-* c)+{-# INLINE fmapL #-}+fmapL = inL . fmap++-- fmapL f+--   = inL (f .)+--   = linear . (f .) . lapply+--   = \ m -> linear (fmap f (lapply m))++(<$>*) = fmapL+++{-# RULES+"fmapL.fmapL"  forall f g m. fmapL f (fmapL g m) = fmapL (f.g) m+ #-}++-- | Apply a /linear/ binary function over linear maps.+liftL2 :: ( LMapDom a s, VectorSpace b s+           , VectorSpace c s, VectorSpace d s) =>+           (b -> c -> d) -> (a :-* b) -> (a :-* c) -> (a :-* d)+liftL2 = inL2 . liftA2++-- liftL2 f a b+--   = inL2 (liftA2 f) a b+--   = linear (liftA2 f (lapply a) (lapply b))++-- I expected the following definition to be equivalent, thanks to rewriting:+-- +--   liftL2 f b c = fmapL f b <*>* c+--     = linear (f . lapply b) <*>* c+--     = linear (lapply (linear (f . lapply b)) <*> lapply c)+--     = linear ((f . lapply b) <*> lapply c)+--     = linear (liftA2 f (lapply b) (lapply c))+--     = (inL2.liftA2) f b c++-- The rewrite isn't happening, however.  And the '(<*>*)' definition+-- yields an incredibly slow implementation.+-- +-- Here's that derivation again, in slo-mo:++--   liftL2 f b c+--     = fmapL f b <*>* c                                  -- inline liftL2+--     = inL (f .) b <*>* c                                -- inline fmapL+--     = (linear . (f .) lapply) b <*>* c                  -- inline inL+--     = linear (f . lapply b) <*>* c                      -- inline (.)+--     = inL2 (<*>) (linear (f . lapply b)) c              -- inline (<*>*)+--     = (inL . (<*>) . lapply) (linear (f . lapply b)) c  -- inline inL2+--     = inL ((<*>) (lapply (linear (f . lapply b)))) c    -- inline (.)+--     = inL ((<*>) (f . lapply b)) c                      -- RULE lapply.linear+--     = (linear . ((<*>) (f . lapply b)) .lapply) c       -- inline inL+--     = linear ((<*>) (f . lapply b) (lapply c))          -- inline (.)+--     = linear ((f . lapply b) <*> (lapply c))            -- infix <*>+--     = linear ((f <$> lapply b) <*> (lapply c))          -- <$> on (a ->)+--     = linear (liftA2 f (lapply b) (lapply c))           -- uninline liftA2++-- When I compile this module, I don't get any firings of lapply.linear+++-- | Apply a /linear/ ternary function over linear maps.+liftL3 :: ( LMapDom a s, VectorSpace b s, VectorSpace c s+           , VectorSpace d s, VectorSpace e s) =>+           (b -> c -> d -> e)+        -> (a :-* b) -> (a :-* c) -> (a :-* d) -> (a :-* e)+liftL3 = inL3 . liftA3+++-- TODO: Get clear about the linearity requirements of 'apL', 'liftL2',+-- and 'liftL3'.++-- | Identity linear map+idL :: LMapDom a s => a :-* a+idL = linear id++-- | Compose linear maps+(.*) :: (VectorSpace c s, LMapDom b s, LMapDom a s) =>+        (b :-* c) -> (a :-* b) -> (a :-* c)+(.*) = inL2 (.)+++instance (VectorSpace v s, LMapDom u s) => AdditiveGroup (u :-* v) where+  zeroV   = pureL  zeroV+  (^+^)   = liftL2 (^+^)+  negateV = fmapL  negateV++instance (VectorSpace v s, LMapDom u s) => VectorSpace (u :-* v) s where+  (*^) s  = fmapL  ((*^) s)++-- Or possibly the following non-standard definition:++-- instance (VectorSpace v s, VectorSpace s s, LMapDom u s)+--       => VectorSpace (u :-* v) (u :-* s) where+--   zeroV   = pureL  zeroV+--   (*^)    = liftL2 (*^)+--   (^+^)   = liftL2 (^+^)+--   negateV = fmapL  negateV++-- Alternatively, add some methods to 'LMapDom' and use them as follows.+-- May be more efficient.++-- -- Linear maps form a vector space+-- instance (VectorSpace o s, LinearMap a o s) => VectorSpace (a :--> o) s where+--   zeroV   = zeroL+--   (^+^)   = addL+--   (*^)    = scaleL+--   negateV = negateL++++--- Instances of LMapDom++instance LMapDom Float Float where+  data Float :-* o   = FloatL o+  lapply (FloatL o)  = (*^ o)+  linear f           = FloatL (f 1)++instance LMapDom Double Double where+  data Double :-* o  = DoubleL o+  lapply (DoubleL o) = (*^ o)+  linear f           = DoubleL (f 1)+++-- | Convenience function for 'linear' definitions.  Both functions are+-- assumed linear.+linearK :: (LMapDom a s) => (a -> b) -> (b -> c) -> a :-* c+linearK k f = linear (f . k)++-- instance LMapDom (Double,Double) Double where+--   data (Double,Double) :-* o = PairD o o+--   PairD ao bo `lapply` (a,b) = a *^ ao ^+^ b *^ bo+--   linear f = PairD (f (0,1)) (f (1,0))++instance (LMapDom a s, LMapDom b s) => LMapDom (a,b) s where+  data (a,b) :-* o = PairL (a :-* o) (b :-* o)+  PairL ao bo `lapply` (a,b) = ao `lapply` a ^+^ bo `lapply` b+  linear f = PairL (linear (\ a -> f (a,zeroV)))+                   (linear (\ b -> f (zeroV,b)))++--   linear = liftA2 PairL+--              (linearK (\ a -> (a,zeroV)))+--              (linearK (\ b -> (zeroV,b)))++instance (LMapDom a s, LMapDom b s, LMapDom c s) => LMapDom (a,b,c) s where+  data (a,b,c) :-* o = TripleL (a :-* o) (b :-* o) (c :-* o)+  TripleL ao bo co `lapply` (a,b,c) =+    ao `lapply` a ^+^ bo `lapply` b ^+^ co `lapply` c+  linear = liftA3 TripleL+             (linearK (\ a -> (a,zeroV,zeroV)))+             (linearK (\ b -> (zeroV,b,zeroV)))+             (linearK (\ c -> (zeroV,zeroV,c)))+++++-- TODO: unfst, unsnd, pair, unpair+++---- OpenGL stuff.++-- I'd rather this code be in a different package.  It's here as a+-- temporary bug workaround.  In ghc-6.8.2, I get the following error+-- message if the 'LMapDom' instance (below) is compiled in a separate+-- module.+-- +--     Type indexes must match class instance head
+--     Found o but expected Vector2 u
+--     In the associated type instance for `:-*'
+--     In the instance declaration for `LMapDom (Vector2 u) s'
++++-- TODO: is UndecidableInstances still necessary?++instance AdditiveGroup u => AdditiveGroup (Vector2 u) where+  zeroV                         = Vector2 zeroV zeroV+  Vector2 u v ^+^ Vector2 u' v' = Vector2 (u^+^u') (v^+^v')+  negateV (Vector2 u v)         = Vector2 (negateV u) (negateV v)++instance (VectorSpace u s) => VectorSpace (Vector2 u) s where+  s *^ Vector2 u v            = Vector2 (s*^u) (s*^v)++instance (InnerSpace u s, AdditiveGroup s)+    => InnerSpace (Vector2 u) s where+  Vector2 u v <.> Vector2 u' v' = u<.>u' ^+^ v<.>v'++instance LMapDom u s => LMapDom (Vector2 u) s where+  data Vector2 u :-* o = VecL (u :-* o) (u :-* o)+  VecL ao bo `lapply` Vector2 a b = ao `lapply` a ^+^ bo `lapply` b+  linear = liftA2 VecL+             (linearK (\ a -> Vector2 a zeroV))+             (linearK (\ b -> Vector2 zeroV b))+++instance AdditiveGroup u => AdditiveGroup (Vector3 u) where+  zeroV                   = Vector3 zeroV zeroV zeroV+  Vector3 u v w ^+^ Vector3 u' v' w'+                          = Vector3 (u^+^u') (v^+^v') (w^+^w')+  negateV (Vector3 u v w) = Vector3 (negateV u) (negateV v) (negateV w)++instance VectorSpace u s => VectorSpace (Vector3 u) s where+  s *^ Vector3 u v w    = Vector3 (s*^u) (s*^v) (s*^w)++instance (InnerSpace u s, AdditiveGroup s)+    => InnerSpace (Vector3 u) s where+  Vector3 u v w <.> Vector3 u' v' w' = u<.>u' ^+^ v<.>v' ^+^ w<.>w'+
+ src/Data/Maclaurin.hs view
@@ -0,0 +1,248 @@+{-# LANGUAGE TypeOperators, MultiParamTypeClasses, UndecidableInstances+           , TypeSynonymInstances, FlexibleInstances, FunctionalDependencies+           , FlexibleContexts+  #-}++-- TODO: remove FlexibleContexts++{-# OPTIONS_GHC -Wall #-}+----------------------------------------------------------------------+-- |+-- Module      :  Data.Maclaurin+-- Copyright   :  (c) Conal Elliott 2008+-- License     :  BSD3+-- +-- Maintainer  :  conal@conal.net+-- Stability   :  experimental+-- +-- Infinite derivative towers via linear maps, using the Maclaurin+-- representation.  See blog posts <http://conal.net/blog/tag/derivatives/>.+----------------------------------------------------------------------++module Data.Maclaurin+  (+    (:>), powVal, derivative, derivativeAt+  , (:~>), dZero, pureD+  , fmapD, (<$>>){-, (<*>>)-}, liftD2, liftD3+  , idD, fstD, sndD+  , linearD, distrib+  , (@.), (>-<)+  ,(**^), (<*.>)+  -- , HasDeriv(..)+  -- experimental+  -- , liftD3+  ) where++-- import Control.Applicative++import Data.VectorSpace+import Data.NumInstances ()+import Data.LinearMap+++infixr 9 `D`, @.+infixl 4 {-<*>>,-} <$>>+infix  0 >-<+++-- | Tower of derivatives.+data a :> b = D { powVal :: b, derivative :: a :-* (a :> b) }++-- | Infinitely differentiable functions+type a :~> b = a -> (a:>b)++-- | Sampled derivative.  For avoiding an awkward typing problem related+-- to the two required 'VectorSpace' instances.+derivativeAt :: (VectorSpace b s, LMapDom a s) =>+                (a :> b) -> a -> (a :> b)+derivativeAt d = lapply (derivative d)++-- The crucial point here is for '($*)' to be interpreted with respect to+-- the 'VectorSpace' instance in this module, not Mac.++-- The argument order for 'derivativeAt' allows partial evaluation, which+-- is useful in power series representations for which 'derivative' is not+-- free (Horner).++-- Handy for missing methods.+noOv :: String -> a+noOv op = error (op ++ ": not defined on a :> b")++-- | Derivative tower full of 'zeroV'.+dZero :: (LMapDom a s, AdditiveGroup b) => a:>b+dZero = pureD zeroV++-- | Constant derivative tower.+pureD :: (LMapDom a s, AdditiveGroup b) => b -> a:>b+pureD b = b `D` pureL dZero++-- | Map a /linear/ function over a derivative tower.+fmapD, (<$>>) :: (LMapDom a s, VectorSpace b s) =>+                 (b -> c) -> (a :> b) -> (a :> c)+fmapD f (D b0 b') = D (f b0) ((fmapL.fmapD) f b')++(<$>>) = fmapD++-- -- | Like '(<*>)' for derivative towers.+-- (<*>>) :: (LMapDom a s, VectorSpace b s, VectorSpace c s) =>+--           (a :> (b -> c)) -> (a :> b) -> (a :> c)+-- D f0 f' <*>> D x0 x' = D (f0 x0) (liftL2 (<*>>) f' x')++-- | Apply a /linear/ binary function over derivative towers.+liftD2 :: (VectorSpace b s, LMapDom a s, VectorSpace c s, VectorSpace d s) =>+          (b -> c -> d) -> (a :> b) -> (a :> c) -> (a :> d)+liftD2 f (D b0 b') (D c0 c') = D (f b0 c0) (liftL2 (liftD2 f) b' c')++-- | Apply a /linear/ ternary function over derivative towers.+liftD3 :: ( LMapDom a s+          , VectorSpace b s, VectorSpace c s+          , VectorSpace d s, VectorSpace e s ) =>+          (b -> c -> d -> e)+       -> (a :> b) -> (a :> c) -> (a :> d) -> (a :> e)+liftD3 f (D b0 b') (D c0 c') (D d0 d') = D (f b0 c0 d0) (liftL3 (liftD3 f) b' c' d')++-- | Differentiable identity function.  Sometimes called "the+-- derivation variable" or similar, but it's not really a variable.+idD :: (LMapDom u s, VectorSpace u s) => u :~> u+idD = linearD id++-- or+--   dId v = D v pureD++-- | Every linear function has a constant derivative equal to the function+-- itself (as a linear map).+linearD :: (LMapDom u s, VectorSpace v s) => (u -> v) -> (u :~> v)+linearD f u = D (f u) (linear (pureD . f))++-- Other examples of linear functions++-- | Differentiable version of 'fst'+fstD :: (VectorSpace a s, LMapDom b s, LMapDom a s) => (a,b) :~> a+fstD = linearD fst++-- | Differentiable version of 'snd'+sndD :: (VectorSpace b s, LMapDom b s, LMapDom a s) => (a,b) :~> b+sndD = linearD snd++-- | Derivative tower for applying a binary function that distributes over+-- addition, such as multiplication.  A bit weaker assumption than+-- bilinearity.+distrib :: (LMapDom a s, VectorSpace b s, VectorSpace c s, VectorSpace u s) =>+           (b -> c -> u) -> (a :> b) -> (a :> c) -> (a :> u)+distrib op = opD+ where+   opD u@(D u0 u') v@(D v0 v') =+     D (u0 `op` v0) (linear (\ da -> u `opD` (v' `lapply` da) ^+^+                                     (u' `lapply` da) `opD` v))++-- Equivalently:+-- +--    opD u@(D u0 u') v@(D v0 v') =+--      D (u0 `op` v0) (linear ((u `opD`) . lapply v' ^+^ (`opD` v) . lapply u'))+-- +-- or+-- +--    opD u@(D u0 u') v@(D v0 v') =+--      D (u0 `op` v0) ( linear ((u `opD`) . lapply v') ^+^+--                       linear ((`opD` v) . lapply u') )+-- or even+-- +--    opD u@(D u0 u') v@(D v0 v') =+--      D (u0 `op` v0) ( inL ((u `opD`) .) v' ^+^ inL ((`opD` v) .) u' )++++++-- TODO: look for a simpler definition of distrib.  this definition almost+-- fits liftLM2.+++-- I'm not sure about the next three, which discard information++instance Show b => Show (a :> b) where show    = noOv "show"+instance Eq   b => Eq   (a :> b) where (==)    = noOv "(==)"+instance Ord  b => Ord  (a :> b) where compare = noOv "compare"++instance (LMapDom a s, VectorSpace u s) => AdditiveGroup (a :> u) where+  zeroV   = pureD  zeroV    -- or dZero+  negateV = fmapD  negateV+  (^+^)   = liftD2 (^+^)++instance (LMapDom a s, VectorSpace u s) => VectorSpace (a :> u) s where+  (*^) s = fmapD  ((*^) s)++(**^) :: (VectorSpace c s, VectorSpace s s, LMapDom a s) =>+         (a :> s) -> (a :> c) -> (a :> c)+(**^) = distrib (*^)++-- ouch!  InnerSpace one won't work at all, for the same reason as for functions.+                  +-- instance (InnerSpace u s) => InnerSpace (a :> u) s where+--   (<.>) = distrib (<.>)++(<*.>) :: (LMapDom a s, InnerSpace b s, VectorSpace s s) =>+          (a :> b) -> (a :> b) -> (a :> s)+(<*.>) s = distrib (<.>) s+++-- The instances below are the one I think we'll want externally.+-- However, the ones above allow the definition of @a:>b@ to work out.+-- The module "Data.Mac" rewraps to provide the alternate instances.++-- instance (LMapDom a s, VectorSpace u s, VectorSpace s s)+--     => VectorSpace (a :> u) (a :> s) where+--   (*^) = (**^)++-- instance (InnerSpace u s, InnerSpace s s', VectorSpace s s, LMapDom a s) =>+--      InnerSpace (a :> u) (a :> s) where+--   (<.>) = (<*.>)+++-- | Chain rule.  See also '(>-<)'.+(@.) :: (LMapDom b s, LMapDom a s, VectorSpace c s) =>+        (b :~> c) -> (a :~> b) -> (a :~> c)+(h @. g) a0 = D c0 (inL2 (@.) c' b')+  where+    D b0 b' = g a0+    D c0 c' = h b0++-- | Specialized chain rule.  See also '(\@.)'+(>-<) :: (LMapDom a s, VectorSpace s s, VectorSpace u s) =>+         (u -> u) -> ((a :> u) -> (a :> s))+      -> (a :> u) -> (a :> u)+f >-< f' = \ u@(D u0 u') -> D (f u0) ((f' u **^) <$>* u')++-- TODO: express '(>-<)' in terms of '(@.)'.  If I can't, then understand why not.++instance (LMapDom a b, Num b, VectorSpace b b) => Num (a:>b) where+  fromInteger = pureD . fromInteger+  (+) = liftD2  (+)+  (-) = liftD2  (-)+  (*) = distrib (*)+  negate = negate >-< -1+  abs    = abs    >-< signum+  signum = signum >-< 0  -- derivative wrong at zero++instance (LMapDom a b, Fractional b, VectorSpace b b) => Fractional (a:>b) where+  fromRational = pureD . fromRational+  recip        = recip >-< recip sqr++sqr :: Num a => a -> a+sqr x = x*x++instance (LMapDom a b, Floating b, VectorSpace b b) => Floating (a:>b) where+  pi    = pureD pi+  exp   = exp   >-< exp+  log   = log   >-< recip+  sqrt  = sqrt  >-< recip (2 * sqrt)+  sin   = sin   >-< cos+  cos   = cos   >-< - sin+  sinh  = sinh  >-< cosh+  cosh  = cosh  >-< sinh+  asin  = asin  >-< recip (sqrt (1-sqr))+  acos  = acos  >-< recip (- sqrt (1-sqr))+  atan  = atan  >-< recip (1+sqr)+  asinh = asinh >-< recip (sqrt (1+sqr))+  acosh = acosh >-< recip (- sqrt (sqr-1))+  atanh = atanh >-< recip (1-sqr)
src/Data/VectorSpace.hs view
@@ -1,5 +1,5 @@ {-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies -           , FlexibleInstances, FlexibleContexts, UndecidableInstances+           , FlexibleInstances, UndecidableInstances  #-} ---------------------------------------------------------------------- -- |@@ -15,54 +15,30 @@  module Data.VectorSpace   ( -    VectorSpace(..), (^-^), (^/), (^*)-  , InnerSpace(..) --, Scalar+    module Data.AdditiveGroup+  , VectorSpace(..), (^/), (^*)+  , InnerSpace(..)   , lerp, magnitudeSq, magnitude, normalized-  , (:-*)   ) where -import Control.Applicative import Data.Complex hiding (magnitude) -infixr 9 :-*+import Data.AdditiveGroup+ infixr 7 *^, ^/, <.> infixl 7 ^*-infixl 6 ^+^, ^-^ --- | Vector space @v@ over a scalar field @s@-class VectorSpace v s | v -> s where-  -- | The zero vector-  zeroV :: v+-- | Vector space @v@ over a scalar field @s@.  Extends 'AdditiveGroup'+-- with scalar multiplication.+class AdditiveGroup v => VectorSpace v s | v -> s where   -- | Scale a vector   (*^)  :: s -> v -> v-  -- | Add vectors-  (^+^) :: v -> v -> v-  -- | Additive inverse-  negateV :: v -> v --- | Adds inner (dot) products-class VectorSpace v s => InnerSpace v s | v -> s where+-- | Adds inner (dot) products.+class VectorSpace v s => InnerSpace v s where   -- | Inner/dot product   (<.>) :: v -> v -> s --- | Convenience.  Maybe add methods later.--- class VectorSpace s s => Scalar s---- TODO: consider replacing v with a type constructor argument:--- --- class VectorSpace v where---   zeroV :: v s---   (*^)  :: s -> v s -> v s---   (^+^) :: v s -> v s -> v s---   (<.>)   :: v s -> v s -> s--- --- Perhaps with constraints on s.  We couldn't then define instances for--- doubles & floats.---- | Vector subtraction-(^-^) :: VectorSpace v s => v -> v -> v-v ^-^ v' = v ^+^ negateV v'- -- | Vector divided by scalar (^/) :: (Fractional s, VectorSpace v s) => v -> s -> v v ^/ s = (1/s) *^ v@@ -77,7 +53,7 @@  -- | Square of the length of a vector.  Sometimes useful for efficiency. -- See also 'magnitude'.-magnitudeSq :: InnerSpace v s =>  v -> s+magnitudeSq :: InnerSpace v s => v -> s magnitudeSq v = v <.> v  -- | Length of a vector.   See also 'magnitudeSq'.@@ -87,34 +63,16 @@ -- | Vector in same direction as given one but with length of one.  If -- given the zero vector, then return it. normalized :: (InnerSpace v s, Floating s) =>  v -> v-normalized v | mag /= 0  = v ^/ mag-             | otherwise = v-  where-    mag = magnitude v--instance VectorSpace Double Double where-  zeroV   = 0.0-  (*^)    = (*)-  (^+^)   = (+)-  negateV = negate--instance InnerSpace Double Double where- (<.>) = (*)+normalized v = v ^/ magnitude v -instance VectorSpace Float Float where-  zeroV   = 0.0-  (*^)    = (*)-  (^+^)   = (+)-  negateV = negate+instance VectorSpace Double Double where (*^)  = (*)+instance InnerSpace  Double Double where (<.>) = (*) -instance InnerSpace Float Float where-  (<.>) = (*)+instance VectorSpace Float  Float  where (*^)  = (*)+instance InnerSpace  Float  Float  where (<.>) = (*)  instance (RealFloat v, VectorSpace v s) => VectorSpace (Complex v) s where-  zeroV       = zeroV :+ zeroV   s*^(u :+ v) = s*^u :+ s*^v-  (^+^)       = (+)-  negateV     = negate  instance (RealFloat v, InnerSpace v s, VectorSpace s s')      => InnerSpace (Complex v) s where@@ -129,10 +87,7 @@ --   Coverage Condition fails for one of the functional dependencies ...)  instance (VectorSpace u s,VectorSpace v s) => VectorSpace (u,v) s where-  zeroV             = (zeroV,zeroV)-  s *^ (u,v)        = (s*^u,s*^v)-  (u,v) ^+^ (u',v') = (u^+^u',v^+^v')-  negateV (u,v)     = (negateV u, negateV v)+  s *^ (u,v) = (s*^u,s*^v)  instance (InnerSpace u s,InnerSpace v s, VectorSpace s s')     => InnerSpace (u,v) s where@@ -143,48 +98,15 @@  instance (VectorSpace u s,VectorSpace v s,VectorSpace w s)     => VectorSpace (u,v,w) s where-  zeroV                  = (zeroV,zeroV,zeroV)-  s *^ (u,v,w)           = (s*^u,s*^v,s*^w)-  (u,v,w) ^+^ (u',v',w') = (u^+^u',v^+^v',w^+^w')-  negateV (u,v,w)        = (negateV u, negateV v, negateV w)+  s *^ (u,v,w) = (s*^u,s*^v,s*^w)  instance (InnerSpace u s,InnerSpace v s,InnerSpace w s, VectorSpace s s')     => InnerSpace (u,v,w) s where   (u,v,w) <.> (u',v',w') = u<.>u' ^+^ v<.>v' ^+^ w<.>w'  --- Standard instance for an applicative functor applied to a vector space.-instance VectorSpace v s => VectorSpace (a->v) s where-  zeroV   = pure   zeroV-  (*^) s  = fmap   (s *^)-  (^+^)   = liftA2 (^+^)-  negateV = fmap   negateV---- I don't know how to make the InnerSpace class work out, because the--- inner product would have to combine two vector *functions* into a--- scalar value.--- ---   instance InnerSpace v s => InnerSpace (a->v) s where---     (<.>) = ???---- Alternatively, we could use (a->s) as the scalar field:--- ---   -- Standard instance for an applicative functor applied to a vector space.---   instance VectorSpace v s => VectorSpace (a->v) (a->s) where---     zeroV   = pure   zeroV---     (*^)    = liftA2 (*^)---     (^+^)   = liftA2 (^+^)---     negateV = fmap negateV--- ---   instance InnerSpace v s => InnerSpace (a->v) (a->s) where---     (<.>) = liftA2 (<.>)--- --- This definition, however, doesn't fit the standard notion of linear--- maps as vector spaces.----- | Linear transformations/maps.  For now, represented as simple--- functions.  The 'VectorSpace' instance for functions gives the usual--- meaning for a vector space of linear transformations.+-- -- Standard instance for an applicative functor applied to a vector space.+-- instance VectorSpace v s => VectorSpace (a->v) s where+--   (*^) s = fmap (s *^) -type a :-* b = a -> b+-- No 'InnerSpace' instance for @(a->v)@.
+ tests/src/Perf.hs view
@@ -0,0 +1,197 @@+{-# LANGUAGE TypeOperators, MultiParamTypeClasses, UndecidableInstances, FlexibleInstances+  #-}+++-- This module tests *performance* of the vector-space operations, such that it is possible to catch performance regressions.+++module Main where++import Control.Applicative+import System.Time+import Data.List++import Data.NumInstances ()+import Data.VectorSpace+import Data.Cross+import Data.Derivative+import Data.LinearMap++type Surf s        = (s,s) -> (s,s,s)+type HeightField s = (s,s) -> s+type Curve2 s      = s -> (s,s)++type Warp1 s        = s -> s+type Warp2 s        = (s,s) -> (s,s)+type Warp3 s        = (s,s,s) -> (s,s,s)++type R = Double++cosU, sinU :: Floating s => s -> s+cosU = cos . mul2pi+sinU = sin . mul2pi++mul2pi :: Floating s => s -> s+mul2pi = (* (2*pi))++torus :: (Floating s, VectorSpace s s) => s -> s -> Surf s+torus sr cr = revolve (\ s -> (sr,0) ^+^ cr *^ circle s)++-- Try use rules to optimize?+-- # RULES "sphere" sphere1 = spec_sphere1+sphere1 :: Floating s => Surf s+sphere1 = revolve semiCircle++spec_sphere1 :: Surf ((Double,Double) :> Double)+spec_sphere1 = sphere1++semiCircle :: Floating s => Curve2 s+semiCircle = circle . (/ 2)++circle :: Floating s => Curve2 s+circle = liftA2 (,) cosU sinU++revolveG :: Floating s => (s -> Curve2 s) -> Surf s+revolveG curveF = \ (u,v) -> onXY (rotate (-2*pi*v)) (addY (curveF v) u)++revolve :: Floating s => Curve2 s -> Surf s+revolve curve = revolveG (const curve)++rotate :: Floating s => s -> Warp2 s+rotate theta = \ (x,y) -> (x * c - y * s, y * c + x * s)+ where c = cos theta+       s = sin theta++addX, addY, addZ :: Num s => (a -> Two s) -> (a -> Three s)+addX = fmap (\ (y,z) -> (0,y,z))+addY = fmap (\ (x,z) -> (x,0,z))+addZ = fmap (\ (x,y) -> (x,y,0))++addYZ,addXZ,addXY :: Num s => (a -> One s) -> (a -> Three s)+addYZ = fmap (\ x -> (x,0,0))+addXZ = fmap (\ y -> (0,y,0))+addXY = fmap (\ z -> (0,0,z))++onX,onY,onZ :: Warp1 s -> Warp3 s+onX f (x,y,z) = (f x, y, z)+onY f (x,y,z) = (x, f y, z)+onZ f (x,y,z) = (x, y, f z)++onXY,onYZ,onXZ :: Warp2 s -> Warp3 s+onXY f (x,y,z) = (x',y',z ) where (x',y') = f (x,y)+onXZ f (x,y,z) = (x',y ,z') where (x',z') = f (x,z)+onYZ f (x,y,z) = (x ,y',z') where (y',z') = f (y,z)+++onX',onY',onZ' :: Warp1 s -> (a -> Three s) -> (a -> Three s)+onX' = fmap fmap onX+onY' = fmap fmap onY+onZ' = fmap fmap onZ++onXY',onXZ',onYZ' :: Warp2 s -> (a -> Three s) -> (a -> Three s)+onXY' = fmap fmap onXY+onXZ' = fmap fmap onXZ+onYZ' = fmap fmap onYZ++displace :: (InnerSpace v s, Floating s, HasNormal v, Applicative f) =>+            f v -> f s -> f v+displace = liftA2 displaceV++displaceV :: (InnerSpace v s, Floating s, HasNormal v) =>+             v -> s -> v+displaceV v s = v ^+^ s *^ normal v++------------------------------------------------------------------------------++surfs3 :: [(Surf ((Double,Double) :> Double),String)]+surfs3 = [ (displace surf hmap,m1 ++ " `displace` " ++ m2) +	 | (surf,m1) <- surfs2+	 , (hmap,m2) <- hmaps+	 ]++surfs2 :: [(Surf ((Double,Double) :> Double),String)]+surfs2 = [ (displace surf hmap,m1 ++ " `displace` " ++ m2) +	 | (surf,m1) <- surfs+	 , (hmap,m2) <- hmaps+	 ]++surfs :: [(Surf ((Double,Double) :> Double),String)]+surfs =+  [ (torus 1 (1/2) ,"torus")+  , (sphere1,"sphere")+  ]++hmaps :: [(HeightField ((Double,Double) :> Double),String)]+hmaps = +  [ (\ (_,_) -> 0,"flat")+  , (\ (u,v) -> cosU u * sinU v,"eggcrate")+  ]++main :: IO ()+main = do +	let loop msg fun t count (points:pss) = do+		sequence_ [ p1 `seq` p2 `seq` p3 `seq` n1 `seq` n2 `seq` n3 `seq` return ()+        	          | (x,y) <- points+		          , let ((p1,p2,p3),(n1,n2,n3)) = vsurf fun (x,y) ]+		diff <- currRelTime t+--		print diff+		if diff > 2+		  then do let count' = count + length points+			  putStrLn $ "Sample count rate for " ++ msg ++ " is " ++ show (fromIntegral count' / diff) ++ " (total count = " ++ show count' ++ ")"+			  return ()+		  else loop msg fun t (count + length points) pss+	    loop _ _ _ _ _ = return ()++	sequence_ [ do t <- getClockTime+		       loop msg fun t 0 (progressive_filter samples_2d)+		  | (fun,msg) <- concat [ surfs, surfs, surfs, surfs2, surfs3 ]+	 	  ]++currRelTime :: ClockTime -> IO Double+currRelTime (TOD sec0 pico0) = fmap delta getClockTime+ where+   delta (TOD sec pico) =+     fromIntegral (sec-sec0) + 1.0e-12 * fromIntegral (pico-pico0)++------------------------------------------------------------------------------++vsurf :: Surf ((R,R) :> R) -> (R,R) -> ((R,R,R),(R,R,R))+vsurf surf = toVN3 . vector3D . surf . unvector2D . idD++type SurfPt s = (s,s) :> (s,s,s)++toVN3 :: (LMapDom s s,Floating s, InnerSpace s s) => SurfPt s -> ((s,s,s),(s,s,s))+toVN3 v = ( powVal v+	  , powVal (normal v)+	  )+vector3D :: (LMapDom a s,VectorSpace s s) => (a :> s,a :> s,a :> s) -> (a :> (s,s,s))+vector3D (u,v,w) = liftD3 (,,) u v w+unvector2D :: (Data.LinearMap.LMapDom a s,VectorSpace s s) => (a :> (s,s)) -> (a :> s,a :> s) +unvector2D d = ( (\ (x,_) -> x) <$>> d+	       , (\ (_,y) -> y) <$>> d+	       )++------------------------------------------------------------------------------++between :: [Double] -> [Double]+between xs = [ (n + m) / 2 | (n,m) <- zip xs (tail xs) ]++samples_1d :: [[Double]]+samples_1d = fn [0,1]+     where+	fn :: [Double] -> [[Double]]+	fn points = points : fn (sort (points ++ between points))++samples_2d :: [[(Double,Double)]]+samples_2d =  [ [ (a,b) +		| a <- sam+		, b <- sam+		]+  	      | sam <- samples_1d+	      ]++-- only allows new points through.+progressive_filter :: (Ord a) => [[a]] -> [[a]]+progressive_filter xs = head sorted_xs : [ y \\ x | (x,y) <- zip sorted_xs (tail sorted_xs) ]+  where+	sorted_xs = map sort xs
vector-space.cabal view
@@ -1,6 +1,7 @@ Name:                vector-space-Version:             0.1.3-Synopsis: 	     Vector & affine spaces, plus derivatives+Version:             0.2.0+Cabal-Version:       >= 1.2+Synopsis:            Vector & affine spaces, plus derivatives Category:            math Description:   vector-space provides classes and generic operations for vector@@ -18,18 +19,35 @@ Author:              Conal Elliott  Maintainer:          conal@conal.net Homepage:            http://haskell.org/haskellwiki/vector-space-Package-Url:	     http://code.haskell.org/vector-space+Package-Url:         http://code.haskell.org/vector-space Copyright:           (c) 2007-2008 by Conal Elliott License:             BSD3 Stability:           experimental-build-type:	     Simple-Hs-Source-Dirs:      src-Extensions:          -Build-Depends:       base-Exposed-Modules:     -		     Data.VectorSpace-		     Data.Derivative-		     Data.AffineSpace-		     Data.NumInstances-		     -ghc-options:         -Wall -O2+build-type:          Simple+Build-Depends:       base, OpenGL+++Library+  Hs-Source-Dirs:      src+  Extensions:          +  Exposed-Modules:     +                     Data.AdditiveGroup+                     Data.VectorSpace+                     Data.LinearMap+                     Data.Maclaurin+                     Data.Derivative+                     Data.Cross+                     Data.AffineSpace+                     Data.NumInstances+  ghc-options:         -Wall -O2+  ghc-prof-options:    -prof -auto-all ++Executable Perf+  main-is:           Perf.hs+  build-depends:     base, OpenGL, old-time+  Hs-Source-Dirs:    src, tests/src+  ghc-options:       -Wall -O2+  ghc-prof-options:    -prof -auto-all   ++--                   Data.Horner+-- The OpenGL dependency is a temporary workaround for a ghc-6.8.2 type family bug.