vector-space 0.7.2 → 0.7.3

raw patch · 6 files changed

+5/−623 lines, 6 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• Makefile +0/−9
• src/Data/AdditiveGroup.hs +4/−2
• src/Data/Horner.hs +0/−220
• src/Data/OldLinearMap.hs +0/−188
• tests/src/Perf.hs +0/−203
• vector-space.cabal +1/−1

@@ -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

@@ -17,8 +17,11 @@   , Sum(..), inSum, inSum2   ) where +import Prelude hiding (foldr)+ import Control.Applicative import Data.Monoid (Monoid(..))+import Data.Foldable (Foldable,foldr) import Data.Complex hiding (magnitude)  import Data.MemoTrie@@ -39,9 +42,8 @@ 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-  instance AdditiveGroup () where   zeroV     = ()

@@ -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)-

@@ -1,188 +0,0 @@-{-# LANGUAGE TypeOperators, FlexibleContexts, TypeFamilies #-}-{-# 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--- License     :  BSD3--- --- Maintainer  :  conal@conal.net--- Stability   :  experimental--- --- Linear maps-------------------------------------------------------------------------module Data.LinearMap-  ( (:-*) , linear, lapply, atBasis, idL, (*.*)-  , liftMS, liftMS2, liftMS3-  , liftL, liftL2, liftL3-  ) where--import Control.Applicative ((<$>),Applicative,liftA2,liftA3)-import Control.Arrow       (first)--import Data.MemoTrie      ((:->:)(..))-import Data.AdditiveGroup (Sum(..),inSum2, 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.---- | An optional additive value-type MSum a = Maybe (Sum a)---- nsum :: MSum a--- nsum = Nothing--jsum :: a -> MSum a-jsum = Just . Sum---- | 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)---- 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@.----- PROBLEM: u :-* v is a type synonym, and Basis is an associated 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----- | Function (assumed linear) as linear map.-linear :: (HasBasis u, HasTrie (Basis u)) =>-          (u -> v) -> (u :-* v)-linear f = jsum (trie (f . basisValue))--atZ :: AdditiveGroup b => (a -> b) -> (MSum a -> b)-atZ f = maybe zeroV (f . getSum)---- 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---- | 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'---- 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)) => -       u :-* u-idL = linear id---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')----- 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 _ Nothing = Nothing--- liftMS h ma = Just (Sum (h (z ma)))--liftMS = fmap.fmap--liftMS2 :: (AdditiveGroup a, AdditiveGroup b) =>-           (a -> b -> c) ->-           (MSum a -> MSum b -> MSum c)-liftMS2 _ Nothing Nothing = Nothing-liftMS2 h ma mb = Just (Sum (h (fromMS ma) (fromMS mb)))--liftMS3 :: (AdditiveGroup a, AdditiveGroup b, AdditiveGroup c) =>-           (a -> b -> c -> d) ->-           (MSum a -> MSum b -> MSum c -> MSum d)-liftMS3 _ Nothing Nothing Nothing = Nothing-liftMS3 h ma mb mc = Just (Sum (h (fromMS ma) (fromMS mb) (fromMS mc)))--fromMS :: AdditiveGroup u => MSum u -> u-fromMS Nothing        = zeroV-fromMS (Just (Sum u)) = u----- | 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 = liftMS . fmap---- | Apply a linear binary function (not to be confused with a bilinear--- function) to each element of a linear map.-liftL2 :: (Applicative f, AdditiveGroup (f a), AdditiveGroup (f b)) =>-          (a -> b -> c)-       -> (MSum (f a) -> MSum (f b) -> MSum (f c))-liftL2 = liftMS2 . liftA2---- | Apply a linear ternary function (not to be confused with a trilinear--- function) to each element of a linear map.-liftL3 :: ( Applicative f-          , AdditiveGroup (f a), AdditiveGroup (f b), AdditiveGroup (f c)) =>-          (a -> b -> c -> d)-       -> (MSum (f a) -> MSum (f b) -> MSum (f c) -> MSum (f d))-liftL3 = liftMS3 . liftA3

@@ -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

@@ -1,5 +1,5 @@ Name:                vector-space-Version:             0.7.2+Version:             0.7.3 Cabal-Version:       >= 1.2 Synopsis:            Vector & affine spaces, linear maps, and derivatives (requires ghc 6.9 or better) Category:            math