diff --git a/Makefile b/Makefile deleted file mode 100644 --- a/Makefile +++ /dev/null @@ -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 diff --git a/src/Data/AdditiveGroup.hs b/src/Data/AdditiveGroup.hs --- a/src/Data/AdditiveGroup.hs +++ b/src/Data/AdditiveGroup.hs @@ -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 = () diff --git a/src/Data/Horner.hs b/src/Data/Horner.hs deleted file mode 100644 --- a/src/Data/Horner.hs +++ /dev/null @@ -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 . ----------------------------------------------------------------------- - -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) - diff --git a/src/Data/OldLinearMap.hs b/src/Data/OldLinearMap.hs deleted file mode 100644 --- a/src/Data/OldLinearMap.hs +++ /dev/null @@ -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 --- --- 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 --- --- 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 diff --git a/tests/src/Perf.hs b/tests/src/Perf.hs deleted file mode 100644 --- a/tests/src/Perf.hs +++ /dev/null @@ -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 diff --git a/vector-space.cabal b/vector-space.cabal --- a/vector-space.cabal +++ b/vector-space.cabal @@ -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