vector-space 0.5.9 → 0.19
raw patch · 14 files changed
Files
- COPYING +25/−0
- Makefile +0/−9
- src/Data/AdditiveGroup.hs +131/−22
- src/Data/AffineSpace.hs +129/−12
- src/Data/Basis.hs +56/−23
- src/Data/Cross.hs +6/−8
- src/Data/Horner.hs +0/−220
- src/Data/LinearMap.hs +176/−60
- src/Data/Maclaurin.hs +28/−41
- src/Data/NumInstances.hs +0/−174
- src/Data/VectorSpace.hs +103/−29
- src/Data/VectorSpace/Generic.hs +20/−0
- tests/src/Perf.hs +0/−203
- vector-space.cabal +21/−18
+ COPYING view
@@ -0,0 +1,25 @@+Copyright (c) 2009 Conal Elliott+All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:+1. Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.+2. Redistributions in binary form must reproduce the above copyright+ notice, this list of conditions and the following disclaimer in the+ documentation and/or other materials provided with the distribution.+3. The names of the authors may not be used to endorse or promote products+ derived from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR+IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES+OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.+IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT,+INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT+NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF+THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.+
− Makefile
@@ -1,9 +0,0 @@-# For special configuration, especially for docs. Otherwise see README.--server = code.haskell.org-server-dir = /srv/code-server-url-dir =--# extra-configure-args += --enable-library-profiling --enable-executable-profiling--include ../my-cabal-make.inc
src/Data/AdditiveGroup.hs view
@@ -1,47 +1,74 @@-{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE TypeOperators, CPP #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE DefaultSignatures #-}+{-# LANGUAGE ScopedTypeVariables #-} ---------------------------------------------------------------------- -- | -- 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(..), (^-^), sumV+ (+ AdditiveGroup(..), sumV , Sum(..), inSum, inSum2 ) where +import Prelude hiding (foldr)+ import Control.Applicative+#if !(MIN_VERSION_base(4,8,0)) import Data.Monoid (Monoid(..))+import Data.Foldable (Foldable)+#endif+import Data.Foldable (foldr) import Data.Complex hiding (magnitude)+import Data.Ratio+#if !(MIN_VERSION_base(4,11,0))+import Data.Semigroup (Semigroup(..))+#endif+import Foreign.C.Types (CSChar, CInt, CShort, CLong, CLLong, CIntMax, CFloat, CDouble) import Data.MemoTrie +import Data.VectorSpace.Generic+import qualified GHC.Generics as Gnrx+import GHC.Generics (Generic, (:*:)(..))+ infixl 6 ^+^, ^-^ -- | Additive group @v@. class AdditiveGroup v where -- | The zero element: identity for '(^+^)' zeroV :: v+ default zeroV :: (Generic v, AdditiveGroup (VRep v)) => v+ zeroV = Gnrx.to (zeroV :: VRep v)+ {-# INLINE zeroV #-} -- | Add vectors (^+^) :: v -> v -> v+ default (^+^) :: (Generic v, AdditiveGroup (VRep v)) => v -> v -> v+ v ^+^ v' = Gnrx.to (Gnrx.from v ^+^ Gnrx.from v' :: VRep v)+ {-# INLINE (^+^) #-} -- | Additive inverse negateV :: v -> v---- | Group subtraction-(^-^) :: AdditiveGroup v => v -> v -> v-v ^-^ v' = v ^+^ negateV v'+ default negateV :: (Generic v, AdditiveGroup (VRep v)) => v -> v+ negateV v = Gnrx.to (negateV $ Gnrx.from v :: VRep v)+ {-# INLINE negateV #-}+ -- | Group subtraction+ (^-^) :: v -> v -> v+ v ^-^ v' = v ^+^ negateV v' -- | Sum over several vectors-sumV :: AdditiveGroup v => [v] -> v+sumV :: (Foldable f, AdditiveGroup v) => f v -> v sumV = foldr (^+^) zeroV-+{-# INLINE sumV #-} instance AdditiveGroup () where zeroV = ()@@ -49,15 +76,28 @@ negateV = id -- For 'Num' types:--- +-- -- instance AdditiveGroup n where {zeroV=0; (^+^) = (+); negateV = negate} -instance AdditiveGroup Int where {zeroV=0; (^+^) = (+); negateV = negate}-instance AdditiveGroup Integer where {zeroV=0; (^+^) = (+); negateV = negate}-instance AdditiveGroup Float where {zeroV=0; (^+^) = (+); negateV = negate}-instance AdditiveGroup Double where {zeroV=0; (^+^) = (+); negateV = negate}+#define ScalarTypeCon(con,t) \+ instance con => AdditiveGroup (t) where {zeroV=0; (^+^) = (+); negateV = negate} +#define ScalarType(t) ScalarTypeCon((),t) +ScalarType(Int)+ScalarType(Integer)+ScalarType(Float)+ScalarType(Double)+ScalarType(CSChar)+ScalarType(CInt)+ScalarType(CShort)+ScalarType(CLong)+ScalarType(CLLong)+ScalarType(CIntMax)+ScalarType(CFloat)+ScalarType(CDouble)+ScalarTypeCon(Integral a,Ratio a)+ instance (RealFloat v, AdditiveGroup v) => AdditiveGroup (Complex v) where zeroV = zeroV :+ zeroV (^+^) = (+)@@ -100,7 +140,37 @@ Just a' ^+^ Just b' = Just (a' ^+^ b') negateV = fmap negateV +{- +Alexey Khudyakov wrote:++ I looked through vector-space package and found lawless instance. Namely Maybe's AdditiveGroup instance++ It's group so following relation is expected to hold. Otherwise it's not a group.+ > x ^+^ negateV x == zeroV++ Here is counterexample:++ > let x = Just 2 in x ^+^ negateV x == zeroV+ False++ I think it's not possible to sensibly define group instance for+ Maybe a at all.+++I see that the problem here is in distinguishing 'Just zeroV' from+Nothing. I could fix the Just + Just line to use Nothing instead of Just+zeroV when a' ^+^ b' == zeroV, although doing so would require Eq a and+hence lose some generality. Even so, the abstraction leak would probably+show up elsewhere.++Hm.++-}++++ -- Memo tries instance (HasTrie u, AdditiveGroup v) => AdditiveGroup (u :->: v) where zeroV = pure zeroV@@ -115,6 +185,7 @@ instance Functor Sum where fmap f (Sum a) = Sum (f a)+ {-# INLINE fmap #-} -- instance Applicative Sum where -- pure a = Sum a@@ -122,32 +193,40 @@ instance Applicative Sum where pure = Sum+ {-# INLINE pure #-} (<*>) = inSum2 ($)+ {-# INLINE (<*>) #-} +instance AdditiveGroup a => Semigroup (Sum a) where+ (<>) = liftA2 (^+^)+ {-# INLINE (<>) #-}+ instance AdditiveGroup a => Monoid (Sum a) where mempty = Sum zeroV- mappend = liftA2 (^+^)-+#if !(MIN_VERSION_base(4,11,0))+ mappend = (<>)+#endif -- | Application a unary function inside a 'Sum' inSum :: (a -> b) -> (Sum a -> Sum b) inSum = getSum ~> Sum+{-# INLINE inSum #-} -- | Application a binary function inside a 'Sum' inSum2 :: (a -> b -> c) -> (Sum a -> Sum b -> Sum c) inSum2 = getSum ~> inSum-+{-# INLINE inSum2 #-} instance AdditiveGroup a => AdditiveGroup (Sum a) where- zeroV = mempty- (^+^) = mappend+ zeroV = Sum zeroV+ (^+^) = (<>) negateV = inSum negateV - ---- to go elsewhere (~>) :: (a' -> a) -> (b -> b') -> ((a -> b) -> (a' -> b')) (i ~> o) f = o . f . i+{-# INLINE (~>) #-} -- result :: (b -> b') -> ((a -> b) -> (a -> b')) -- result = (.)@@ -156,3 +235,33 @@ -- argument = flip (.) -- g ~> f = result g . argument f++++instance AdditiveGroup a => AdditiveGroup (Gnrx.Rec0 a s) where+ zeroV = Gnrx.K1 zeroV+ {-# INLINE zeroV #-}+ negateV (Gnrx.K1 v) = Gnrx.K1 $ negateV v+ {-# INLINE negateV #-}+ Gnrx.K1 v ^+^ Gnrx.K1 w = Gnrx.K1 $ v ^+^ w+ {-# INLINE (^+^) #-}+ Gnrx.K1 v ^-^ Gnrx.K1 w = Gnrx.K1 $ v ^-^ w+ {-# INLINE (^-^) #-}+instance AdditiveGroup (f p) => AdditiveGroup (Gnrx.M1 i c f p) where+ zeroV = Gnrx.M1 zeroV+ {-# INLINE zeroV #-}+ negateV (Gnrx.M1 v) = Gnrx.M1 $ negateV v+ {-# INLINE negateV #-}+ Gnrx.M1 v ^+^ Gnrx.M1 w = Gnrx.M1 $ v ^+^ w+ {-# INLINE (^+^) #-}+ Gnrx.M1 v ^-^ Gnrx.M1 w = Gnrx.M1 $ v ^-^ w+ {-# INLINE (^-^) #-}+instance (AdditiveGroup (f p), AdditiveGroup (g p)) => AdditiveGroup ((f :*: g) p) where+ zeroV = zeroV :*: zeroV+ {-# INLINE zeroV #-}+ negateV (x:*:y) = negateV x :*: negateV y+ {-# INLINE negateV #-}+ (x:*:y) ^+^ (ξ:*:υ) = (x^+^ξ) :*: (y^+^υ)+ {-# INLINE (^+^) #-}+ (x:*:y) ^-^ (ξ:*:υ) = (x^-^ξ) :*: (y^-^υ)+ {-# INLINE (^-^) #-}
src/Data/AffineSpace.hs view
@@ -1,4 +1,10 @@-{-# LANGUAGE MultiParamTypeClasses, FlexibleContexts, TypeFamilies #-}+{-# LANGUAGE MultiParamTypeClasses, FlexibleContexts, TypeFamilies, CPP #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE DefaultSignatures #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE DeriveGeneric #-} ---------------------------------------------------------------------- -- | -- Module : Data.AffineSpace@@ -13,25 +19,51 @@ module Data.AffineSpace (- AffineSpace(..), (.-^), distanceSq, distance, alerp+ AffineSpace(..), (.-^), distanceSq, distance, alerp, affineCombo ) where-+#if !MIN_VERSION_base(4,10,0) import Control.Applicative (liftA2)+#endif+import Data.Ratio+import Foreign.C.Types (CSChar, CInt, CShort, CLong, CLLong, CIntMax, CFloat, CDouble)+import Control.Arrow(first) import Data.VectorSpace+import Data.Basis -infix 4 .+^, .-^, .-.+import Data.VectorSpace.Generic+import qualified GHC.Generics as Gnrx+import GHC.Generics (Generic, (:*:)(..)) +-- Through 0.8.4, I used the following fixities.+-- +-- infix 4 .+^, .-^, .-.+-- +-- Changed in 0.8.5 to match precedence of + and -, and to associate usefully.+-- Thanks to Ben Gamari for suggesting left-associativity.++infixl 6 .+^, .-^+infix 6 .-.++ -- TODO: Convert AffineSpace from fundep to associated type, and eliminate -- FunctionalDependencies above. class AdditiveGroup (Diff p) => AffineSpace p where -- | Associated vector space type Diff p+ type Diff p = GenericDiff p -- | Subtract points (.-.) :: p -> p -> Diff p+ default (.-.) :: ( Generic p, Diff p ~ GenericDiff p, AffineSpace (VRep p) )+ => p -> p -> Diff p+ p .-. q = GenericDiff+ $ (Gnrx.from p .-. (Gnrx.from q :: VRep p)) -- | Point plus vector (.+^) :: p -> Diff p -> p+ default (.+^) :: ( Generic p, Diff p ~ GenericDiff p, AffineSpace (VRep p) )+ => p -> Diff p -> p+ p .+^ GenericDiff q = Gnrx.to (Gnrx.from p .+^ q :: VRep p) -- | Point minus vector (.-^) :: AffineSpace p => p -> Diff p -> p@@ -55,17 +87,46 @@ p -> p -> Scalar (Diff p) -> p alerp p p' s = p .+^ (s *^ (p' .-. p)) -instance AffineSpace Double where- type Diff Double = Double- (.-.) = (-)- (.+^) = (+)+-- | Compute an affine combination (weighted average) of points.+-- The first element is used as origin and is weighted+-- such that all coefficients sum to 1. For example,+--+-- > affineCombo a [(0.3,b), (0.2,c)]+--+-- is equal to+--+-- > a .+^ (0.3 *^ (b .-. a) ^+^ 0.2 *^ (c .-. a))+--+-- and if @a@, @b@, and @c@ were in a vector space would also be equal to+--+-- > 0.5 *^ a ^+^ 0.3 *^ b ^+^ 0.2 *^ c+--+-- See also 'linearCombo' (on vector spaces).+affineCombo :: (AffineSpace p, v ~ Diff p, VectorSpace v) => p -> [(p,Scalar v)] -> p+affineCombo z l = z .+^ linearCombo (map (first (.-. z)) l) -instance AffineSpace Float where- type Diff Float = Float- (.-.) = (-)- (.+^) = (+)+#define ScalarTypeCon(con,t) \+ instance con => AffineSpace (t) where \+ { type Diff (t) = t \+ ; (.-.) = (-) \+ ; (.+^) = (+) } +#define ScalarType(t) ScalarTypeCon((),t) +ScalarType(Int)+ScalarType(Integer)+ScalarType(Double)+ScalarType(Float)+ScalarType(CSChar)+ScalarType(CInt)+ScalarType(CShort)+ScalarType(CLong)+ScalarType(CLLong)+ScalarType(CIntMax)+ScalarType(CDouble)+ScalarType(CFloat)+ScalarTypeCon(Integral a,Ratio a)+ instance (AffineSpace p, AffineSpace q) => AffineSpace (p,q) where type Diff (p,q) = (Diff p, Diff q) (p,q) .-. (p',q') = (p .-. p', q .-. q')@@ -81,3 +142,59 @@ type Diff (a -> p) = a -> Diff p (.-.) = liftA2 (.-.) (.+^) = liftA2 (.+^)++++newtype GenericDiff p = GenericDiff (Diff (VRep p))+ deriving (Generic)++instance AdditiveGroup (Diff (VRep p)) => AdditiveGroup (GenericDiff p)+instance VectorSpace (Diff (VRep p)) => VectorSpace (GenericDiff p)+instance (AdditiveGroup (Scalar (Diff (VRep p))), InnerSpace (Diff (VRep p))) => InnerSpace (GenericDiff p)+instance HasBasis (Diff (VRep p)) => HasBasis (GenericDiff p)++data AffineDiffProductSpace f g p = AffineDiffProductSpace+ !(Diff (f p)) !(Diff (g p)) deriving (Generic)+instance (AffineSpace (f p), AffineSpace (g p))+ => AdditiveGroup (AffineDiffProductSpace f g p)+instance ( AffineSpace (f p), AffineSpace (g p)+ , VectorSpace (Diff (f p)), VectorSpace (Diff (g p))+ , Scalar (Diff (f p)) ~ Scalar (Diff (g p)) )+ => VectorSpace (AffineDiffProductSpace f g p)+instance ( AdditiveGroup (Scalar (Diff (g p)))+ , AffineSpace (f p), AffineSpace (g p)+ , InnerSpace (Diff (f p)), InnerSpace (Diff (g p))+ , Scalar (Diff (f p)) ~ Scalar (Diff (g p))+ , Num (Scalar (Diff (f p))) )+ => InnerSpace (AffineDiffProductSpace f g p)+instance (AffineSpace (f p), AffineSpace (g p))+ => AffineSpace (AffineDiffProductSpace f g p) where+ type Diff (AffineDiffProductSpace f g p) = AffineDiffProductSpace f g p+ (.+^) = (^+^)+ (.-.) = (^-^)+instance ( AffineSpace (f p), AffineSpace (g p)+ , HasBasis (Diff (f p)), HasBasis (Diff (g p))+ , Scalar (Diff (f p)) ~ Scalar (Diff (g p)) )+ => HasBasis (AffineDiffProductSpace f g p) where+ type Basis (AffineDiffProductSpace f g p) = Either (Basis (Diff (f p)))+ (Basis (Diff (g p)))+ basisValue (Left bf) = AffineDiffProductSpace (basisValue bf) zeroV+ basisValue (Right bg) = AffineDiffProductSpace zeroV (basisValue bg)+ decompose (AffineDiffProductSpace vf vg)+ = map (first Left) (decompose vf) ++ map (first Right) (decompose vg)+ decompose' (AffineDiffProductSpace vf _) (Left bf) = decompose' vf bf+ decompose' (AffineDiffProductSpace _ vg) (Right bg) = decompose' vg bg+++instance AffineSpace a => AffineSpace (Gnrx.Rec0 a s) where+ type Diff (Gnrx.Rec0 a s) = Diff a+ Gnrx.K1 v .+^ w = Gnrx.K1 $ v .+^ w+ Gnrx.K1 v .-. Gnrx.K1 w = v .-. w+instance AffineSpace (f p) => AffineSpace (Gnrx.M1 i c f p) where+ type Diff (Gnrx.M1 i c f p) = Diff (f p)+ Gnrx.M1 v .+^ w = Gnrx.M1 $ v .+^ w+ Gnrx.M1 v .-. Gnrx.M1 w = v .-. w+instance (AffineSpace (f p), AffineSpace (g p)) => AffineSpace ((f :*: g) p) where+ type Diff ((f:*:g) p) = AffineDiffProductSpace f g p+ (x:*:y) .+^ AffineDiffProductSpace ξ υ = (x.+^ξ) :*: (y.+^υ)+ (x:*:y) .-. (ξ:*:υ) = AffineDiffProductSpace (x.-.ξ) (y.-.υ)
src/Data/Basis.hs view
@@ -1,10 +1,8 @@--- WARNING: this module depends on type families working fairly well, and--- requires ghc version at least 6.9. I didn't find a way to specify that--- dependency in the .cabal.--- {-# LANGUAGE TypeOperators, TypeFamilies, UndecidableInstances- , FlexibleInstances, MultiParamTypeClasses- #-}+ , FlexibleInstances, MultiParamTypeClasses, CPP #-}+{-# LANGUAGE DefaultSignatures #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE ScopedTypeVariables #-} {-# OPTIONS_GHC -Wall -fno-warn-orphans #-} ---------------------------------------------------------------------- -- |@@ -23,37 +21,51 @@ -- import Control.Applicative ((<$>)) import Control.Arrow (first)-import Data.Either+import Data.Ratio+import Foreign.C.Types (CFloat, CDouble)+import Data.Kind+-- import Data.Either import Data.VectorSpace +import Data.VectorSpace.Generic+import qualified GHC.Generics as Gnrx+import GHC.Generics (Generic, (:*:)(..))+ -- using associated data type instead of associated type synonym to work -- around ghc bug <http://hackage.haskell.org/trac/ghc/ticket/3038> class VectorSpace v => HasBasis v where -- | Representation of the canonical basis for @v@- type Basis v :: *+ type Basis v :: Type+ type Basis v = Basis (VRep v) -- | Interpret basis rep as a vector basisValue :: Basis v -> v+ default basisValue :: (Generic v, HasBasis (VRep v), Basis (VRep v) ~ Basis v)+ => Basis v -> v+ basisValue b = Gnrx.to (basisValue b :: VRep v) -- | Extract coordinates decompose :: v -> [(Basis v, Scalar v)]+ default decompose :: ( Generic v, HasBasis (VRep v)+ , Scalar (VRep v) ~ Scalar v, Basis (VRep v) ~ Basis v )+ => v -> [(Basis v, Scalar v)]+ decompose v = decompose (Gnrx.from v :: VRep v) -- | Experimental version. More elegant definitions, and friendly to -- infinite-dimensional vector spaces. decompose' :: v -> (Basis v -> Scalar v)+ default decompose' :: ( Generic v, HasBasis (VRep v)+ , Scalar (VRep v) ~ Scalar v, Basis (VRep v) ~ Basis v )+ => v -> Basis v -> Scalar v+ decompose' v = decompose' (Gnrx.from v :: VRep v) -- Defining property: recompose . decompose == id --- | Linear combination-linearCombo :: VectorSpace v => [(v,Scalar v)] -> v-linearCombo ps = sumV [s *^ v | (v,s) <- ps]- -- Turn a basis decomposition back into a vector. recompose :: HasBasis v => [(Basis v, Scalar v)] -> v recompose = linearCombo . fmap (first basisValue) -- recompose ps = linearCombo (first basisValue <$> ps) - -- I don't know how to define -- -- recompose' :: HasBasis v => (Basis v -> Scalar v) -> v@@ -61,18 +73,21 @@ -- However, I don't seem to use recompose anywhere. -- I don't even use basisValue or decompose. -instance HasBasis Float where- type Basis Float = ()- basisValue () = 1- decompose s = [((),s)]- decompose' s = const s+#define ScalarTypeCon(con,t) \+ instance con => HasBasis (t) where \+ { type Basis (t) = () \+ ; basisValue () = 1 \+ ; decompose s = [((),s)] \+ ; decompose' s = const s } -instance HasBasis Double where- type Basis Double = ()- basisValue () = 1- decompose s = [((),s)]- decompose' s = const s+#define ScalarType(t) ScalarTypeCon((),t) +ScalarType(Float)+ScalarType(CFloat)+ScalarType(Double)+ScalarType(CDouble)+ScalarTypeCon(Integral a, Ratio a)+ instance ( HasBasis u, s ~ Scalar u , HasBasis v, s ~ Scalar v ) => HasBasis (u,v) where@@ -136,3 +151,21 @@ t4 = basisValue (Right (Left ())) :: (Float,Double,Float) -}++instance HasBasis a => HasBasis (Gnrx.Rec0 a s) where+ type Basis (Gnrx.Rec0 a s) = Basis a+ basisValue = Gnrx.K1 . basisValue+ decompose = decompose . Gnrx.unK1+ decompose' = decompose' . Gnrx.unK1+instance HasBasis (f p) => HasBasis (Gnrx.M1 i c f p) where+ type Basis (Gnrx.M1 i c f p) = Basis (f p)+ basisValue = Gnrx.M1 . basisValue+ decompose = decompose . Gnrx.unM1+ decompose' = decompose' . Gnrx.unM1+instance (HasBasis (f p), HasBasis (g p), Scalar (f p) ~ Scalar (g p))+ => HasBasis ((f :*: g) p) where+ type Basis ((f:*:g) p) = Either (Basis (f p)) (Basis (g p))+ basisValue (Left bf) = basisValue bf :*: zeroV+ basisValue (Right bg) = zeroV :*: basisValue bg+ decompose (u:*:v) = decomp2 Left u ++ decomp2 Right v+ decompose' (u:*:v) = decompose' u `either` decompose' v
src/Data/Cross.hs view
@@ -1,6 +1,6 @@ {-# LANGUAGE FlexibleInstances, FlexibleContexts, TypeOperators- , TypeFamilies, TypeSynonymInstances- #-}+ , TypeFamilies, TypeSynonymInstances + , UndecidableInstances #-} {-# OPTIONS_GHC -Wall #-} ---------------------------------------------------------------------- -- |@@ -49,8 +49,7 @@ instance AdditiveGroup u => HasCross2 (u,u) where cross2 (x,y) = (negateV y,x) -- or @(y,-x)@? -instance ( HasBasis a, HasTrie (Basis a)- , VectorSpace v, HasCross2 v) => HasCross2 (a:>v) where+instance (HasTrie (Basis a), HasCross2 v) => HasCross2 (a:>v) where -- 2d cross-product is linear cross2 = fmapD cross2 @@ -74,8 +73,7 @@ -- l `atB` b = maybe zeroV (`untrie` b) l -instance ( Num s, VectorSpace s- , HasBasis s, HasTrie (Basis s), Basis s ~ ())+instance (VectorSpace s, HasBasis s, HasTrie (Basis s), Basis s ~ ()) => HasNormal (Two (One s :> s)) where normalVec = unpairD . normalVec . pairD @@ -102,7 +100,7 @@ where d = derivAtBasis v -instance ( Num s, VectorSpace s, HasBasis s, HasTrie (Basis s)- , HasNormal (Two s :> Three s))+instance ( VectorSpace s, HasBasis s, HasTrie (Basis s)+ , HasNormal (Two s :> Three s) ) => HasNormal (Three (Two s :> s)) where normalVec = untripleD . normalVec . tripleD
− src/Data/Horner.hs
@@ -1,220 +0,0 @@-{-# LANGUAGE TypeOperators, MultiParamTypeClasses- , TypeSynonymInstances, FlexibleInstances- #-}-{-# 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
@@ -1,34 +1,40 @@-{-# LANGUAGE TypeOperators, FlexibleContexts, TypeFamilies #-}+{-# LANGUAGE TupleSections #-}+{-# LANGUAGE CPP, TypeOperators, FlexibleContexts, TypeFamilies+ , GeneralizedNewtypeDeriving, StandaloneDeriving, UndecidableInstances #-} {-# OPTIONS_GHC -Wall -fno-warn-orphans #-}--- {-# OPTIONS_GHC -funbox-strict-fields #-}--- {-# OPTIONS_GHC -ddump-simpl-stats -ddump-simpl #-} ---------------------------------------------------------------------- -- | -- Module : Data.LinearMap--- Copyright : (c) Conal Elliott 2008+-- Copyright : (c) Conal Elliott 2008-2016 -- License : BSD3--- +-- -- Maintainer : conal@conal.net -- Stability : experimental--- +-- -- Linear maps ---------------------------------------------------------------------- module Data.LinearMap- ( (:-*) , linear, lapply, atBasis, idL, (*.*)- , liftMS, liftMS2, liftMS3- , liftL, liftL2, liftL3- ) where+ ( (:-*) , linear, lapply, atBasis, idL, (*.*)+ , inLMap, inLMap2, inLMap3+ , liftMS, liftMS2, liftMS3+ , liftL, liftL2, liftL3+ , exlL, exrL, forkL, firstL, secondL+ , inlL, inrL, joinL -- , leftL, rightL+ )+ where -import Control.Applicative ((<$>),Applicative,liftA2,liftA3)-import Control.Arrow (first)+#if !(MIN_VERSION_base(4,8,0))+import Control.Applicative (Applicative, liftA2)+#endif+import Control.Applicative (liftA3)+import Control.Arrow (first,second) -import Data.MemoTrie ((:->:)(..))-import Data.AdditiveGroup (Sum(..),inSum2, AdditiveGroup(..))+import Data.MemoTrie (HasTrie(..),(:->:))+import Data.AdditiveGroup (Sum(..), AdditiveGroup(..)) import Data.VectorSpace (VectorSpace(..)) import Data.Basis (HasBasis(..), linearCombo) - -- Linear maps are almost but not quite a Control.Category. The type -- class constraints interfere. They're almost an Arrow also, but for the -- constraints and the generality of arr.@@ -36,35 +42,97 @@ -- | An optional additive value type MSum a = Maybe (Sum a) --- nsum :: MSum a--- nsum = Nothing- jsum :: a -> MSum a jsum = Just . Sum +type LMap' u v = MSum (Basis u :->: v)++infixr 1 :-* -- | Linear map, represented as an optional memo-trie from basis to -- values, where 'Nothing' means the zero map (an optimization).-type u :-* v = MSum (Basis u :->: v)+newtype u :-* v = LMap { unLMap :: LMap' u v } --- TODO: Try a partial trie instead, excluding (known) zero elements.--- Then 'lapply' could be much faster for sparse situations. Make sure to--- correctly sum them. It'd be more like Jason Foutz's formulation--- <http://metavar.blogspot.com/2008/02/higher-order-multivariate-automatic.html>--- which uses in @IntMap@.+deriving instance (HasTrie (Basis u), AdditiveGroup v) => AdditiveGroup (u :-* v) +instance (HasTrie (Basis u), VectorSpace v) =>+ VectorSpace (u :-* v) where+ type Scalar (u :-* v) = Scalar v+ (*^) s = (inLMap.liftMS.fmap) (s *^) --- PROBLEM: u :-* v is a type synonym, and Basis is an associated type synonym, resulting in a subtle+-- In GHC 7.10:+-- Constraint is no smaller than the instance head+-- in the constraint: HasTrie (Basis u)+-- (Use UndecidableInstances to permit this)++exlL :: ( HasBasis a, HasTrie (Basis a), HasBasis b, HasTrie (Basis b)+ , Scalar a ~ Scalar b )+ => (a,b) :-* a+exlL = linear fst++exrL :: ( HasBasis a, HasTrie (Basis a), HasBasis b, HasTrie (Basis b)+ , Scalar a ~ Scalar b )+ => (a,b) :-* b+exrL = linear snd++forkL :: (HasTrie (Basis a), HasBasis c, HasBasis d)+ => (a :-* c) -> (a :-* d) -> (a :-* (c,d))+forkL = (inLMap2.liftL2) (,)++firstL :: ( HasBasis u, HasBasis u', HasBasis v+ , HasTrie (Basis u), HasTrie (Basis v) + , Scalar u ~ Scalar v, Scalar u ~ Scalar u'+ ) =>+ (u :-* u') -> ((u,v) :-* (u',v))+firstL = linear.first.lapply++secondL :: ( HasBasis u, HasBasis v, HasBasis v'+ , HasTrie (Basis u), HasTrie (Basis v) + , Scalar u ~ Scalar v, Scalar u ~ Scalar v'+ ) =>+ (v :-* v') -> ((u,v) :-* (u,v'))+secondL = linear.second.lapply++-- TODO: more efficient firstL++-- liftMS :: (AdditiveGroup a) => (a -> b) -> (MSum a -> MSum b)++-- (s *^) :: v -> v+-- fmap (s *^) :: (Basis u :->: v) -> (Basis u :->: v)+-- (liftMS.fmap) (s *^) :: LMap' u v -> LMap' u v+-- (inLMap.liftMS.fmap) (s *^) :: (u :-* v) -> (u :-* v)+++inlL :: (HasBasis a, HasTrie (Basis a), HasBasis b)+ => a :-* (a,b)+inlL = linear (,zeroV)++inrL :: (HasBasis a, HasBasis b, HasTrie (Basis b))+ => b :-* (a,b)+inrL = linear (zeroV,)++joinL :: ( HasBasis a, HasTrie (Basis a)+ , HasBasis b, HasTrie (Basis b)+ , Scalar a ~ Scalar b, Scalar a ~ Scalar c+ , VectorSpace c )+ => (a :-* c) -> (b :-* c) -> ((a,b) :-* c)+f `joinL` g = linear (\ (a,b) -> lapply f a ^+^ lapply g b)++-- Before version 0.7, u :-* v was a type synonym, resulting in a subtle -- ambiguity: u:-*v == u':-*v' does not imply that u==u', since Basis -- might map different types to the same basis (e.g., Float & Double).--- See <http://hackage.haskell.org/trac/ghc/ticket/1897>--- --- Work in progress. See NewLinearMap.hs+-- See <http://hackage.haskell.org/trac/ghc/ticket/1897>.+-- See also <http://thread.gmane.org/gmane.comp.lang.haskell.cafe/73271/focus=73332>. +-- TODO: Try a partial trie instead, excluding (known) zero elements.+-- Then 'lapply' could be much faster for sparse situations. Make sure to+-- correctly sum them. It'd be more like Jason Foutz's formulation+-- <http://metavar.blogspot.com/2008/02/higher-order-multivariate-automatic.html>+-- which uses in @IntMap@. -- | Function (assumed linear) as linear map. linear :: (HasBasis u, HasTrie (Basis u)) => (u -> v) -> (u :-* v)-linear f = jsum (trie (f . basisValue))+linear f = LMap (jsum (trie (f . basisValue))) atZ :: AdditiveGroup b => (a -> b) -> (MSum a -> b) atZ f = maybe zeroV (f . getSum)@@ -72,78 +140,88 @@ -- atZ :: AdditiveGroup b => (a -> b) -> (a -> b) -- atZ = id --- | Evaluate a linear map on a basis element. I've loosened the type to--- work around a typing problem in 'derivAtBasis'.--- atBasis :: (AdditiveGroup v, HasTrie (Basis u)) =>--- (u :-* v) -> Basis u -> v-atBasis :: (HasTrie a, AdditiveGroup b) => MSum (a :->: b) -> a -> b-m `atBasis` b = atZ (`untrie` b) m+inLMap :: (LMap' r s -> LMap' t u) -> ((r :-* s) -> (t :-* u))+inLMap = unLMap ~> LMap +inLMap2 :: (LMap' r s -> LMap' t u -> LMap' v w)+ -> ((r :-* s) -> (t :-* u) -> (v :-* w))+inLMap2 = unLMap ~> inLMap++inLMap3 :: (LMap' r s -> LMap' t u -> LMap' v w -> LMap' x y)+ -> ((r :-* s) -> (t :-* u) -> (v :-* w) -> (x :-* y))+inLMap3 = unLMap ~> inLMap2+ -- | Apply a linear map to a vector. lapply :: ( VectorSpace v, Scalar u ~ Scalar v , HasBasis u, HasTrie (Basis u) ) => (u :-* v) -> (u -> v)-lapply = atZ lapply'+lapply = atZ lapply' . unLMap --- Handy for 'lapply' and '(*.*)'.+-- | Evaluate a linear map on a basis element.+atBasis :: (AdditiveGroup v, HasTrie (Basis u)) =>+ (u :-* v) -> Basis u -> v+LMap m `atBasis` b = atZ (`untrie` b) m++-- | Handy for 'lapply' and '(*.*)'. lapply' :: ( VectorSpace v, Scalar u ~ Scalar v , HasBasis u, HasTrie (Basis u) ) => (Basis u :->: v) -> (u -> v) lapply' tr = linearCombo . fmap (first (untrie tr)) . decompose ----- Identity linear map-idL :: (HasBasis u, HasTrie (Basis u)) => +-- | Identity linear map+idL :: (HasBasis u, HasTrie (Basis u)) => u :-* u idL = linear id infixr 9 *.* -- | Compose linear maps-(*.*) :: ( HasBasis u, HasTrie (Basis u)+(*.*) :: ( HasTrie (Basis u) , HasBasis v, HasTrie (Basis v) , VectorSpace w , Scalar v ~ Scalar w ) => (v :-* w) -> (u :-* v) -> (u :-* w) -- Simple definition, but only optimizes out uv == zero--- --- (*.*) vw = (fmap.fmap) (lapply vw) +-- vw *.* uv = LMap ((fmap.fmap.fmap) (lapply vw) (unLMap uv))++(*.*) vw = (inLMap.fmap.fmap.fmap) (lapply vw)++-- Eep:+-- (*.*) = inLMap.fmap.fmap.fmap.lapply++ -- Instead, use Nothing/zero if /either/ map is zeroV (exploiting linearity -- when uv == zeroV.) --- Nothing *.* _ = Nothing--- _ *.* Nothing = Nothing--- Just (Sum vw) *.* Just (Sum uv) = Just (Sum (lapply' vw <$> uv))+-- LMap Nothing *.* _ = LMap Nothing+-- _ *.* LMap Nothing = LMap Nothing+-- LMap (Just (Sum vw)) *.* LMap (Just (Sum uv)) = LMap (Just (Sum (lapply' vw <$> uv))) --- (*.*) = liftA2 (\ (Sum vw) (Sum uv) -> Sum (lapply' vw <$> uv))+-- (*.*) = liftA2 (\ (LMap (Sum vw)) (LMap (Sum uv)) -> LMap (Sum (lapply' vw <$> uv))) --- (*.*) = (liftA2.inSum2) (\ vw uv -> lapply' vw <$> uv)-(*.*) = (liftA2.inSum2) (\ vw uv -> lapply' vw <$> uv)+-- (*.*) = (liftA2.inSum2.inLMap2) (\ vw uv -> lapply' vw <$> uv) --- (*.*) = (liftA2.inSum2) (\ vw -> fmap (lapply' vw))+-- (*.*) = (liftA2.inSum2.inLMap2) (\ vw -> fmap (lapply' vw)) --- (*.*) = (liftA2.inSum2) (fmap . lapply')+-- (*.*) = (liftA2.inSum2.inLMap2) (fmap . lapply') -- It may be helpful that @lapply vw@ is evaluated just once and not -- once per uv. 'untrie' can strip off all of its trie constructors. -- Less efficient definition:--- +-- -- vw `compL` uv = linear (lapply vw . lapply uv)--- +-- -- i.e., compL = inL2 (.)--- +-- -- The problem with these definitions is that basis elements get converted -- to values and then decomposed, followed by recombination of the -- results. -liftMS :: (AdditiveGroup a) =>- (a -> b)- -> (MSum a -> MSum b)+liftMS :: (a -> b) -> (MSum a -> MSum b) -- liftMS _ Nothing = Nothing -- liftMS h ma = Just (Sum (h (z ma))) @@ -168,8 +246,7 @@ -- | Apply a linear function to each element of a linear map. -- @liftL f l == linear f *.* l@, but works more efficiently.-liftL :: (Functor f, AdditiveGroup (f a)) =>- (a -> b) -> MSum (f a) -> MSum (f b)+liftL :: Functor f => (a -> b) -> MSum (f a) -> MSum (f b) liftL = liftMS . fmap -- | Apply a linear binary function (not to be confused with a bilinear@@ -186,3 +263,42 @@ (a -> b -> c -> d) -> (MSum (f a) -> MSum (f b) -> MSum (f c) -> MSum (f d)) liftL3 = liftMS3 . liftA3++{-+++infixr 9 *.*+-- | Compose linear maps+(*.*) :: ( HasBasis u, HasTrie (Basis u)+ , HasBasis v, HasTrie (Basis v)+ , VectorSpace w+ , Scalar v ~ Scalar w ) =>+ (v :-* w) -> (u :-* v) -> (u :-* w)++-- Simple definition, but only optimizes out uv == zero+--+-- (*.*) vw = (fmap.fmap) (lapply vw)++-- Instead, use Nothing/zero if /either/ map is zeroV (exploiting linearity+-- when uv == zeroV.)++-- Nothing *.* _ = Nothing+-- _ *.* Nothing = Nothing+-- Just (Sum vw) *.* Just (Sum uv) = Just (Sum (lapply' vw <$> uv))++-- (*.*) = liftA2 (\ (Sum vw) (Sum uv) -> Sum (lapply' vw <$> uv))++-- (*.*) = (liftA2.inSum2) (\ vw uv -> lapply' vw <$> uv)+(*.*) = (liftA2.inSum2) (\ vw uv -> lapply' vw <$> uv)++-- (*.*) = (liftA2.inSum2) (\ vw -> fmap (lapply' vw))++-- (*.*) = (liftA2.inSum2) (fmap . lapply')+++-}++-----++(~>) :: (a' -> a) -> (b -> b') -> ((a -> b) -> (a' -> b'))+(f ~> h) g = h . g . f
src/Data/Maclaurin.hs view
@@ -1,8 +1,7 @@ {-# LANGUAGE TypeOperators, MultiParamTypeClasses, UndecidableInstances , TypeSynonymInstances, FlexibleInstances , FlexibleContexts, TypeFamilies- , ScopedTypeVariables- #-}+ , ScopedTypeVariables, CPP #-} -- The ScopedTypeVariables is there just as a bug work-around. Without it -- I get a bogus error about context mismatch for mutually recursive@@ -25,7 +24,7 @@ -- Stability : experimental -- -- Infinite derivative towers via linear maps, using the Maclaurin--- representation. See blog posts <http://conal.net/blog/tag/derivatives/>.+-- representation. See blog posts <http://conal.net/blog/tag/derivative/>. ---------------------------------------------------------------------- module Data.Maclaurin@@ -53,6 +52,10 @@ import Data.Boolean +#if MIN_VERSION_base(4,8,0)+import Prelude hiding ((<*))+#endif+ infixr 9 `D` -- | Tower of derivatives. data a :> b = D { powVal :: b, derivative :: a :-* (a :> b) }@@ -71,11 +74,10 @@ infixl 4 <$>> -- | Map a /linear/ function over a derivative tower.-fmapD, (<$>>) :: (HasBasis a, HasTrie (Basis a), AdditiveGroup b) =>- (b -> c) -> (a :> b) -> (a :> c)+fmapD, (<$>>) :: HasTrie (Basis a) => (b -> c) -> (a :> b) -> (a :> c) fmapD f = lf where- lf (D b0 b') = D (f b0) (liftL lf b')+ lf (D b0 b') = D (f b0) ((inLMap.liftL) lf b') (<$>>) = fmapD @@ -84,7 +86,7 @@ (b -> c -> d) -> (a :> b) -> (a :> c) -> (a :> d) liftD2 f = lf where- lf (D b0 b') (D c0 c') = D (f b0 c0) (liftL2 lf b' c')+ lf (D b0 b') (D c0 c') = D (f b0 c0) ((inLMap2.liftL2) lf b' c') -- | Apply a /linear/ ternary function over derivative towers.@@ -95,7 +97,7 @@ liftD3 f = lf where lf (D b0 b') (D c0 c') (D d0 d') =- D (f b0 c0 d0) (liftL3 lf b' c' d')+ D (f b0 c0 d0) ((inLMap3.liftL3) lf b' c' d') -- TODO: Can liftD2 and liftD3 be defined in terms of a (<*>>) similar to@@ -108,9 +110,7 @@ -- | Differentiable identity function. Sometimes called "the -- derivation variable" or similar, but it's not really a variable.-idD :: ( VectorSpace u, s ~ Scalar u- , VectorSpace (u :> u), VectorSpace s- , HasBasis u, HasTrie (Basis u)) =>+idD :: (VectorSpace u , HasBasis u, HasTrie (Basis u)) => u :~> u idD = linearD id @@ -165,16 +165,14 @@ -- | Derivative tower for applying a binary function that distributes over -- addition, such as multiplication. A bit weaker assumption than -- bilinearity. Is bilinearity necessary for correctness here?-distrib :: forall a b c u.- ( HasBasis a, HasTrie (Basis a)- , AdditiveGroup b, AdditiveGroup c, AdditiveGroup u) =>+distrib :: forall a b c u. (HasBasis a, HasTrie (Basis a) , AdditiveGroup u) => (b -> c -> u) -> (a :> b) -> (a :> c) -> (a :> u) distrib op = (#) where u@(D u0 u') # v@(D v0 v') =- D (u0 `op` v0) ( liftMS (inTrie ((# v) .)) u' ^+^- liftMS (inTrie ((u #) .)) v' )+ D (u0 `op` v0) ( (inLMap.liftMS) (inTrie ((# v) .)) u' ^+^+ (inLMap.liftMS) (inTrie ((u #) .)) v' ) -- TODO: I think this distrib is exponential in increasing degree. Switch@@ -187,18 +185,19 @@ instance Show b => Show (a :> b) where show (D b0 _) = "D " ++ show b0 ++ " ..." -instance Eq b => Eq (a :> b) where (==) = noOv "(==)"+instance Eq (a :> b) where (==) = noOv "(==)" -instance (AdditiveGroup v, HasBasis u, HasTrie (Basis u), IfB b v) =>- IfB b (u :> v) where+type instance BooleanOf (a :> b) = BooleanOf b++instance (AdditiveGroup v, HasBasis u, HasTrie (Basis u), IfB v) =>+ IfB (u :> v) where ifB = liftD2 . ifB -instance (AdditiveGroup v, HasBasis u, HasTrie (Basis u), OrdB b v) =>- OrdB b (u :> v) where+instance OrdB v => OrdB (u :> v) where (<*) = (<*) `on` powVal instance ( AdditiveGroup b, HasBasis a, HasTrie (Basis a)- , OrdB bool b, IfB bool b, Ord b) => Ord (a :> b) where+ , OrdB b, IfB b, Ord b) => Ord (a :> b) where compare = compare `on` powVal min = minB max = maxB@@ -213,8 +212,7 @@ -- Less efficient: adds zero -- (^+^) = liftD2 (^+^) -instance ( HasBasis a, HasTrie (Basis a)- , VectorSpace u, AdditiveGroup (Scalar u) )+instance (HasBasis a, HasTrie (Basis a), VectorSpace u) => VectorSpace (a :> u) where type Scalar (a :> u) = (a :> Scalar u) (*^) = distrib (*^) @@ -236,11 +234,10 @@ infix 0 >-< -- | Specialized chain rule. See also '(\@.)'-(>-<) :: ( HasBasis a, HasTrie (Basis a), VectorSpace u- , AdditiveGroup (Scalar u)) =>+(>-<) :: (HasBasis a, HasTrie (Basis a), VectorSpace u) => (u -> u) -> ((a :> u) -> (a :> Scalar u)) -> (a :> u) -> (a :> u)-f >-< f' = \ u@(D u0 u') -> D (f u0) (liftMS (f' u *^) u')+f >-< f' = \ u@(D u0 u') -> D (f u0) ((inLMap.liftMS) (f' u *^) u') -- TODO: express '(>-<)' in terms of '(@.)'. If I can't, then understand why not.@@ -293,31 +290,21 @@ ---- Misc -pairD :: ( HasBasis a, HasTrie (Basis a)- , VectorSpace b, VectorSpace c- , Scalar b ~ Scalar c- ) => (a:>b,a:>c) -> a:>(b,c)+pairD :: (HasBasis a, HasTrie (Basis a), VectorSpace b, VectorSpace c)+ => (a:>b,a:>c) -> a:>(b,c) pairD (u,v) = liftD2 (,) u v -unpairD :: ( HasBasis a, HasTrie (Basis a)- , VectorSpace a, VectorSpace b, VectorSpace c- , Scalar b ~ Scalar c- ) => (a :> (b,c)) -> (a:>b, a:>c)+unpairD :: HasTrie (Basis a) => (a :> (b,c)) -> (a:>b, a:>c) unpairD d = (fst <$>> d, snd <$>> d) tripleD :: ( HasBasis a, HasTrie (Basis a) , VectorSpace b, VectorSpace c, VectorSpace d- , Scalar b ~ Scalar c, Scalar c ~ Scalar d ) => (a:>b,a:>c,a:>d) -> a:>(b,c,d) tripleD (u,v,w) = liftD3 (,,) u v w -untripleD :: ( HasBasis a, HasTrie (Basis a)- , VectorSpace a, VectorSpace b, VectorSpace c, VectorSpace d- , Scalar b ~ Scalar c, Scalar c ~ Scalar d- ) =>- (a :> (b,c,d)) -> (a:>b, a:>c, a:>d)+untripleD :: HasTrie (Basis a) => (a :> (b,c,d)) -> (a:>b, a:>c, a:>d) untripleD d = ((\ (a,_,_) -> a) <$>> d, (\ (_,b,_) -> b) <$>> d, (\ (_,_,c) -> c) <$>> d)
− src/Data/NumInstances.hs
@@ -1,174 +0,0 @@-{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}-{-# OPTIONS_GHC -fno-warn-orphans #-}-------------------------------------------------------------------------- |--- Module : Data.NumInstances--- Copyright : (c) Conal Elliott 2008--- License : BSD3--- --- Maintainer : conal@conal.net--- Stability : experimental--- --- Number class instances for functions and tuples-------------------------------------------------------------------------module Data.NumInstances () where--import Control.Applicative--noOv :: String -> String -> a-noOv ty meth = error $ meth ++ ": No overloading for " ++ ty--noFun :: String -> a-noFun = noOv "function"---- Eq & Show are prerequisites for Num, so they need to be faked here-instance Eq (a->b) where- (==) = noFun "(==)"- (/=) = noFun "(/=)"--instance Ord b => Ord (a->b) where- min = liftA2 min- max = liftA2 max--instance Show (a->b) where- show = noFun "show"- showsPrec = noFun "showsPrec"- showList = noFun "showList"--instance Num b => Num (a->b) where- negate = fmap negate- (+) = liftA2 (+)- (*) = liftA2 (*)- fromInteger = pure . fromInteger- abs = fmap abs- signum = fmap signum--instance Fractional b => Fractional (a->b) where- recip = fmap recip- fromRational = pure . fromRational--instance Floating b => Floating (a->b) where- pi = pure pi- sqrt = fmap sqrt- exp = fmap exp- log = fmap log- sin = fmap sin- cos = fmap cos- asin = fmap asin- atan = fmap atan- acos = fmap acos- sinh = fmap sinh- cosh = fmap cosh- asinh = fmap asinh- atanh = fmap atanh- acosh = fmap acosh-------- Tuples--lift2 :: (a->u) -> (b->v) -> (a,b) -> (u,v)-lift2 f g (a,b) = (f a, g b)---- Equivalently, lift2 = (***)--instance (Num a, Num b) => Num (a,b) where- fromInteger n = (fromInteger n, fromInteger n)- (a,b) + (a',b') = (a+a',b+b')- (a,b) - (a',b') = (a-a',b-b')- (a,b) * (a',b') = (a*a',b*b')- negate = lift2 negate negate- abs = lift2 abs abs- signum = lift2 signum signum--instance (Fractional a, Fractional b) => Fractional (a,b) where- fromRational x = (fromRational x, fromRational x)- recip = lift2 recip recip--instance (Floating a, Floating b) => Floating (a,b) where- pi = (pi,pi)- exp = lift2 exp exp- log = lift2 log log- sqrt = lift2 sqrt sqrt- sin = lift2 sin sin- cos = lift2 cos cos- sinh = lift2 sinh sinh- cosh = lift2 cosh cosh- asin = lift2 asin asin- acos = lift2 acos acos- atan = lift2 atan atan- asinh = lift2 asinh asinh- acosh = lift2 acosh acosh- atanh = lift2 atanh atanh--instance (Num a, Num b, Num c) => Num (a,b,c) where- fromInteger n = (fromInteger n, fromInteger n, fromInteger n)- (a,b,c) + (a',b',c') = (a+a',b+b',c+c')- (a,b,c) - (a',b',c') = (a-a',b-b',c-c')- (a,b,c) * (a',b',c') = (a*a',b*b',c*c')- negate = lift3 negate negate negate- abs = lift3 abs abs abs- signum = lift3 signum signum signum--instance (Fractional a, Fractional b, Fractional c)- => Fractional (a,b,c) where- fromRational x = (fromRational x, fromRational x, fromRational x)- recip = lift3 recip recip recip---lift3 :: (a->u) -> (b->v) -> (c->w) -> (a,b,c) -> (u,v,w)-lift3 f g h (a,b,c) = (f a, g b, h c)--instance (Floating a, Floating b, Floating c)- => Floating (a,b,c) where- pi = (pi,pi,pi)- exp = lift3 exp exp exp- log = lift3 log log log- sqrt = lift3 sqrt sqrt sqrt- sin = lift3 sin sin sin- cos = lift3 cos cos cos- sinh = lift3 sinh sinh sinh- cosh = lift3 cosh cosh cosh- asin = lift3 asin asin asin- acos = lift3 acos acos acos- atan = lift3 atan atan atan- asinh = lift3 asinh asinh asinh- acosh = lift3 acosh acosh acosh- atanh = lift3 atanh atanh atanh----lift4 :: (a->u) -> (b->v) -> (c->w) -> (d->x)- -> (a,b,c,d) -> (u,v,w,x)-lift4 f g h k (a,b,c,d) = (f a, g b, h c, k d)--instance (Num a, Num b, Num c, Num d) => Num (a,b,c,d) where- fromInteger n = (fromInteger n, fromInteger n, fromInteger n, fromInteger n)- (a,b,c,d) + (a',b',c',d') = (a+a',b+b',c+c',d+d')- (a,b,c,d) - (a',b',c',d') = (a-a',b-b',c-c',d-d')- (a,b,c,d) * (a',b',c',d') = (a*a',b*b',c*c',d*d')- negate = lift4 negate negate negate negate- abs = lift4 abs abs abs abs- signum = lift4 signum signum signum signum--instance (Fractional a, Fractional b, Fractional c, Fractional d)- => Fractional (a,b,c,d) where- fromRational x = (fromRational x, fromRational x, fromRational x, fromRational x)- recip = lift4 recip recip recip recip--instance (Floating a, Floating b, Floating c, Floating d)- => Floating (a,b,c,d) where- pi = (pi,pi,pi,pi)- exp = lift4 exp exp exp exp- log = lift4 log log log log- sqrt = lift4 sqrt sqrt sqrt sqrt- sin = lift4 sin sin sin sin- cos = lift4 cos cos cos cos- sinh = lift4 sinh sinh sinh sinh- cosh = lift4 cosh cosh cosh cosh- asin = lift4 asin asin asin asin- acos = lift4 acos acos acos acos- atan = lift4 atan atan atan atan- asinh = lift4 asinh asinh asinh asinh- acosh = lift4 acosh acosh acosh acosh- atanh = lift4 atanh atanh atanh atanh
src/Data/VectorSpace.hs view
@@ -1,18 +1,21 @@ {-# LANGUAGE MultiParamTypeClasses, TypeOperators- , TypeFamilies, UndecidableInstances- #-}+ , TypeFamilies, UndecidableInstances, CPP+ , FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE DefaultSignatures #-}+{-# LANGUAGE ScopedTypeVariables #-} {-# OPTIONS_GHC -Wall #-} ---------------------------------------------------------------------- -- | -- Module : Data.VectorSpace -- Copyright : (c) Conal Elliott and Andy J Gill 2008 -- License : BSD3--- +-- -- Maintainer : conal@conal.net, andygill@ku.edu -- Stability : experimental--- +-- -- Vector spaces--- +-- -- This version uses associated types instead of fundeps and -- requires ghc-6.10 or later ----------------------------------------------------------------------@@ -20,83 +23,130 @@ -- NB: I'm attempting to replace fundeps with associated types. See -- NewVectorSpace.hs. Ran into trouble with type equality constraints not -- getting propagated. Manuel Ch is looking into it.--- +-- -- Blocking bug: http://hackage.haskell.org/trac/ghc/ticket/2448 module Data.VectorSpace ( module Data.AdditiveGroup , VectorSpace(..), (^/), (^*) , InnerSpace(..)- , lerp, magnitudeSq, magnitude, normalized+ , lerp, linearCombo, magnitudeSq, magnitude, normalized, project ) where-+#if !(MIN_VERSION_base(4,8,0)) import Control.Applicative (liftA2)+#endif import Data.Complex hiding (magnitude)+import Foreign.C.Types (CSChar, CInt, CShort, CLong, CLLong, CIntMax, CFloat, CDouble)+import Data.Ratio+import Data.Kind import Data.AdditiveGroup import Data.MemoTrie +import Data.VectorSpace.Generic+import qualified GHC.Generics as Gnrx+import GHC.Generics (Generic, (:*:)(..))+ infixr 7 *^ -- | Vector space @v@. class AdditiveGroup v => VectorSpace v where- type Scalar v :: *+ type Scalar v :: Type+ type Scalar v = Scalar (VRep v) -- | Scale a vector (*^) :: Scalar v -> v -> v+ default (*^) :: (Generic v, VectorSpace (VRep v), Scalar (VRep v) ~ Scalar v)+ => Scalar v -> v -> v+ μ *^ v = Gnrx.to (μ *^ Gnrx.from v :: VRep v)+ {-# INLINE (*^) #-} infixr 7 <.> -- | Adds inner (dot) products.-class VectorSpace v => InnerSpace v where+class (VectorSpace v, AdditiveGroup (Scalar v)) => InnerSpace v where -- | Inner/dot product (<.>) :: v -> v -> Scalar v+ default (<.>) :: (Generic v, InnerSpace (VRep v), Scalar (VRep v) ~ Scalar v)+ => v -> v -> Scalar v+ v<.>w = (Gnrx.from v :: VRep v) <.> Gnrx.from w+ {-# INLINE (<.>) #-} infixr 7 ^/ infixl 7 ^* -- | Vector divided by scalar (^/) :: (VectorSpace v, s ~ Scalar v, Fractional s) => v -> s -> v-v ^/ s = (1/s) *^ v+v ^/ s = recip s *^ v+{-# INLINE (^/) #-} -- | Vector multiplied by scalar (^*) :: (VectorSpace v, s ~ Scalar v) => v -> s -> v (^*) = flip (*^)+{-# INLINE (^*) #-} -- | Linear interpolation between @a@ (when @t==0@) and @b@ (when @t==1@). -- lerp :: (VectorSpace v, s ~ Scalar v, Num s) => v -> v -> s -> v lerp :: VectorSpace v => v -> v -> Scalar v -> v lerp a b t = a ^+^ t *^ (b ^-^ a)+{-# INLINE lerp #-} +-- | Linear combination of vectors+linearCombo :: VectorSpace v => [(v,Scalar v)] -> v+linearCombo ps = sumV [v ^* s | (v,s) <- ps]+{-# INLINE linearCombo #-}+ -- | Square of the length of a vector. Sometimes useful for efficiency. -- See also 'magnitude'. magnitudeSq :: (InnerSpace v, s ~ Scalar v) => v -> s magnitudeSq v = v <.> v+{-# INLINE magnitudeSq #-} -- | Length of a vector. See also 'magnitudeSq'. magnitude :: (InnerSpace v, s ~ Scalar v, Floating s) => v -> s magnitude = sqrt . magnitudeSq+{-# INLINE magnitude #-} --- | Vector in same direction as given one but with length of one. If--- given the zero vector, then return it.+-- | Vector in same direction as given one but with length of one.+-- Divides by zero for the zero vector. normalized :: (InnerSpace v, s ~ Scalar v, Floating s) => v -> v normalized v = v ^/ magnitude v+{-# INLINE normalized #-} -instance VectorSpace Double where- type Scalar Double = Double- (*^) = (*)-instance InnerSpace Double where (<.>) = (*)+-- | @project u v@ computes the projection of @v@ onto @u@.+project :: (InnerSpace v, s ~ Scalar v, Fractional s) => v -> v -> v+project u v = ((v <.> u) / magnitudeSq u) *^ u+{-# INLINE project #-} -instance VectorSpace Float where- type Scalar Float = Float- (*^) = (*)-instance InnerSpace Float where (<.>) = (*)+#define ScalarType(t) \+ instance VectorSpace (t) where \+ { type Scalar (t) = (t) \+ ; (*^) = (*) } ; \+ instance InnerSpace (t) where (<.>) = (*) +ScalarType(Int)+ScalarType(Integer)+ScalarType(Double)+ScalarType(Float)+ScalarType(CSChar)+ScalarType(CInt)+ScalarType(CShort)+ScalarType(CLong)+ScalarType(CLLong)+ScalarType(CIntMax)+ScalarType(CDouble)+ScalarType(CFloat)++instance Integral a => VectorSpace (Ratio a) where+ type Scalar (Ratio a) = Ratio a+ (*^) = (*)+instance Integral a => InnerSpace (Ratio a) where (<.>) = (*)+ instance (RealFloat v, VectorSpace v) => VectorSpace (Complex v) where type Scalar (Complex v) = Scalar v s*^(u :+ v) = s*^u :+ s*^v -instance (RealFloat v, InnerSpace v, s ~ Scalar v, AdditiveGroup s)+instance (RealFloat v, InnerSpace v) => InnerSpace (Complex v) where (u :+ v) <.> (u' :+ v') = (u <.> u') ^+^ (v <.> v') @@ -116,8 +166,7 @@ s *^ (u,v) = (s*^u,s*^v) instance ( InnerSpace u, s ~ Scalar u- , InnerSpace v, s ~ Scalar v- , AdditiveGroup (Scalar v) )+ , InnerSpace v, s ~ Scalar v ) => InnerSpace (u,v) where (u,v) <.> (u',v') = (u <.> u') ^+^ (v <.> v') @@ -130,8 +179,7 @@ instance ( InnerSpace u, s ~ Scalar u , InnerSpace v, s ~ Scalar v- , InnerSpace w, s ~ Scalar w- , AdditiveGroup s )+ , InnerSpace w, s ~ Scalar w ) => InnerSpace (u,v,w) where (u,v,w) <.> (u',v',w') = u<.>u' ^+^ v<.>v' ^+^ w<.>w' @@ -146,8 +194,7 @@ instance ( InnerSpace u, s ~ Scalar u , InnerSpace v, s ~ Scalar v , InnerSpace w, s ~ Scalar w- , InnerSpace x, s ~ Scalar x- , AdditiveGroup s )+ , InnerSpace x, s ~ Scalar x ) => InnerSpace (u,v,w,x) where (u,v,w,x) <.> (u',v',w',x') = u<.>u' ^+^ v<.>v' ^+^ w<.>w' ^+^ x<.>x' @@ -185,7 +232,7 @@ -instance (InnerSpace a, AdditiveGroup (Scalar a)) => InnerSpace (Maybe a) where+instance InnerSpace a => InnerSpace (Maybe a) where -- dotting with zero (vector) yields zero (scalar) Nothing <.> _ = zeroV _ <.> Nothing = zeroV@@ -194,3 +241,30 @@ -- mu <.> mv = fromMaybe zeroV (liftA2 (<.>) mu mv) -- (<.>) = (fmap.fmap) (fromMaybe zeroV) (liftA2 (<.>))+++instance VectorSpace a => VectorSpace (Gnrx.Rec0 a s) where+ type Scalar (Gnrx.Rec0 a s) = Scalar a+ μ *^ Gnrx.K1 v = Gnrx.K1 $ μ*^v+ {-# INLINE (*^) #-}+instance VectorSpace (f p) => VectorSpace (Gnrx.M1 i c f p) where+ type Scalar (Gnrx.M1 i c f p) = Scalar (f p)+ μ *^ Gnrx.M1 v = Gnrx.M1 $ μ*^v+ {-# INLINE (*^) #-}+instance (VectorSpace (f p), VectorSpace (g p), Scalar (f p) ~ Scalar (g p))+ => VectorSpace ((f :*: g) p) where+ type Scalar ((f:*:g) p) = Scalar (f p)+ μ *^ (x:*:y) = μ*^x :*: μ*^y+ {-# INLINE (*^) #-}++instance InnerSpace a => InnerSpace (Gnrx.Rec0 a s) where+ Gnrx.K1 v <.> Gnrx.K1 w = v<.>w+ {-# INLINE (<.>) #-}+instance InnerSpace (f p) => InnerSpace (Gnrx.M1 i c f p) where+ Gnrx.M1 v <.> Gnrx.M1 w = v<.>w+ {-# INLINE (<.>) #-}+instance ( InnerSpace (f p), InnerSpace (g p)+ , Scalar (f p) ~ Scalar (g p), Num (Scalar (f p)) )+ => InnerSpace ((f :*: g) p) where+ (x:*:y) <.> (ξ:*:υ) = x<.>ξ + y<.>υ+ {-# INLINE (<.>) #-}
+ src/Data/VectorSpace/Generic.hs view
@@ -0,0 +1,20 @@+-- |+-- Module : Data.VectorSpace.Generic+-- Copyright : (c) Conal Elliott and Justus Sagemüller 2017+-- License : BSD3+-- +-- Maintainer : conal@conal.net, (@) jsagemue $ uni-koeln.de+-- Stability : experimental+-- +-- Underpinnings of the type that represents vector / affine / etc. spaces+-- with GHC generics++module Data.VectorSpace.Generic where+++import qualified GHC.Generics as Gnrx++import Data.Void+++type VRep v = Gnrx.Rep v Void
− tests/src/Perf.hs
@@ -1,203 +0,0 @@-{-# LANGUAGE TypeOperators, MultiParamTypeClasses, UndecidableInstances, FlexibleInstances- , TypeFamilies, FlexibleContexts- #-}----- 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.Basis-import Data.MemoTrie-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 ()-- let samples = samples_2d-- sequence_ [ do t <- getClockTime- loop msg fun t 0 samples- | (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 :: (HasBasis s s, Basis s ~ (), Floating s, InnerSpace s s)- => SurfPt s -> ((s,s,s),(s,s,s))-toVN3 v = ( powVal v- , powVal (normal v)- )-vector3D :: (HasBasis a s, HasTrie (Basis a), VectorSpace s s) => (a :> s,a :> s,a :> s) -> (a :> (s,s,s))-vector3D (u,v,w) = liftD3 (,,) u v w-unvector2D :: (HasBasis a s, HasTrie (Basis a), 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,7 +1,7 @@ Name: vector-space-Version: 0.5.9-Cabal-Version: >= 1.2-Synopsis: Vector & affine spaces, linear maps, and derivatives (requires ghc 6.9 or better)+Version: 0.19+Cabal-Version: >= 1.10+Synopsis: Vector & affine spaces, linear maps, and derivatives Category: math Description: /vector-space/ provides classes and generic operations for vector@@ -11,24 +11,30 @@ (scalars, vectors, matrices, ...). . /Warning/: this package depends on type families working fairly well,- and requires ghc version at least 6.9.+ requiring GHC version at least 6.9. . Project wiki page: <http://haskell.org/haskellwiki/vector-space> .- © 2008 by Conal Elliott; BSD3 license.+ © 2008-2012 by Conal Elliott; BSD3 license. Author: Conal Elliott Maintainer: conal@conal.net-Homepage: http://haskell.org/haskellwiki/vector-space-Package-Url: http://code.haskell.org/vector-space-Copyright: (c) 2008 by Conal Elliott+Copyright: (c) 2008-2012 by Conal Elliott License: BSD3+License-File: COPYING Stability: experimental build-type: Simple +source-repository head+ type: git+ location: git://github.com/conal/vector-space.git+ Library+ default-language: Haskell2010 hs-Source-Dirs: src Extensions: - Build-Depends: base, MemoTrie >= 0.4.2, Boolean+ Build-Depends: base<5, MemoTrie >= 0.5+ , Boolean >= 0.1.0+ , NumInstances >= 1.0 Exposed-Modules: Data.AdditiveGroup Data.VectorSpace@@ -39,14 +45,11 @@ Data.Derivative Data.Cross Data.AffineSpace- Data.NumInstances-+ Other-Modules: + Data.VectorSpace.Generic -- This library relies on type families working as well as in 6.10.- if impl(ghc < 6.10) {- buildable: False- }- ghc-options: -Wall -O2- ghc-prof-options: -prof -auto-all ---- For ghc-options: -ddump-simpl-stats -ddump-rules -ddump-simpl -ddump-simpl-phases+ if impl(ghc < 6.10) { buildable: False }+ if !impl(ghc >= 7.6) { Build-Depends: ghc-prim >= 0.2 }+ if !impl(ghc >= 7.9) { Build-Depends: void >= 0.4 }+ if !impl(ghc >= 8.0) { Build-Depends: semigroups >= 0.16 }