clifford 0.1.0.6 → 0.1.0.7
raw patch · 8 files changed
+148/−16 lines, 8 filesdep ~basedep ~criteriondep ~numeric-prelude
Dependency ranges changed: base, criterion, numeric-prelude, stream-fusion
Files
- bench/benchmarks.hs +39/−0
- changelog.md +1/−0
- clifford.cabal +12/−2
- src/Numeric/Clifford/Internal.hs +42/−0
- src/Numeric/Clifford/LinearOperators.lhs +26/−0
- src/Numeric/Clifford/Multivector.lhs +24/−10
- src/Numeric/Clifford/NumericIntegration.lhs +1/−1
- src/Numeric/Clifford/NumericIntegration/DefaultIntegrators.hs +3/−3
+ bench/benchmarks.hs view
@@ -0,0 +1,39 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE NoMonomorphismRestriction #-}+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE DataKinds #-}+{-# OPTIONS_GHC -fllvm -fexcess-precision -optlo-O3 -optlc-O=3 -O3 #-}++import Numeric.Clifford.Multivector+import Numeric.Clifford.NumericIntegration.DefaultIntegrators+import Criterion.Main+import Data.List.Stream+import NumericPrelude hiding (iterate, last, map, take, log)+import Prelude hiding (iterate, last, map, negate, take,log, (*),+ (+))+ +type STVector = Multivector 3 1 Double+scalar2 = scalar (2.0::NumericPrelude.Double) :: STVector+ij2 = (2.0::NumericPrelude.Double) `e` [1,2] :: STVector +ik3 = (3::NumericPrelude.Double) `e` [1,3] :: STVector +ijk4 = (4::NumericPrelude.Double) `e` [1,2,3] :: STVector +ijl5 = (5::NumericPrelude.Double) `e` [1,2,4] :: STVector+a = ij2 + ik3 + ijk4 + ijl5 + (scalar 1.5)+enormousThing = a*a*a*a*a*a*a + scalar2+expDecay _ x = map negate $ map ((*) (1.3 `e` [] :: STVector)) x+thelambda init = lobattoIIIAFourthOrder 0.01 expDecay init+main = defaultMain [+ bgroup "log" [ bench "scalar 2.0" $ nf log scalar2+ , bench "2ij" $ nf log ij2+ , bench "3ik" $ nf log ik3+ , bench "4ijk" $ nf log ijk4+ , bench "5ijl" $ nf log ijl5+ , bench "sum" $ nf log a+ , bench "enormous thing" $ nf log enormousThing+ ],+ bgroup "lobatto IIIA 4th order RK solver"+ [+ bench "200 iterations exponential decay" $ nf (\x -> last $ take 200 (iterate thelambda x)) (0.0,[scalar 1.0])+ ]+ ]
changelog.md view
@@ -1,4 +1,5 @@ -*-change-log-*-+ 0.1.0.7 Adding basic linear operators; made multivector a field 0.1.0.6 Memoising the blade index comparision function for a 20% speed increase 0.1.0.5 Adding hspec tests, fixed blade comparison to order blades in terms of grade first
clifford.cabal view
@@ -10,7 +10,7 @@ -- PVP summary: +-+------- breaking API changes -- | | +----- non-breaking API additions -- | | | +--- code changes with no API change-version: 0.1.0.6+version: 0.1.0.7 -- A short (one-line) description of the package. synopsis: A Clifford algebra library@@ -51,7 +51,7 @@ library -- Modules exported by the library.- exposed-modules: Numeric.Clifford.Blade, Numeric.Clifford.Multivector, Numeric.Clifford.NumericIntegration, Numeric.Clifford.NumericIntegration.DefaultIntegrators, Numeric.Clifford.ClassicalMechanics, Numeric.Clifford.LinearOperators+ exposed-modules: Numeric.Clifford.Blade, Numeric.Clifford.Multivector, Numeric.Clifford.NumericIntegration, Numeric.Clifford.NumericIntegration.DefaultIntegrators, Numeric.Clifford.ClassicalMechanics, Numeric.Clifford.LinearOperators, Numeric.Clifford.Internal -- Modules included in this library but not exported. -- other-modules: @@ -63,6 +63,7 @@ data-ordlist >= 0.4.5 && < 0.5, converge >= 0.1.0.1 && < 0.2, lens >= 4.0.3 && < 4.1, deepseq >= 1.3.0.1 && < 1.4, vector >= 0.10.0.1 && < 0.11, stream-fusion >= 0.1 && < 0.2, criterion >= 0.8.0.0 && < 0.9, derive, QuickCheck, nats, tagged, cereal,hspec, MemoTrie >= 0.6 && < 0.7 + ghc-options: -fllvm -fexcess-precision -optlc-O=3 -O3 -- Directories containing source files. hs-source-dirs: src @@ -78,3 +79,12 @@ main-is: Spec.lhs build-depends: base, clifford, hspec, numeric-prelude, QuickCheck, nats ++benchmark basic-ops+ type: exitcode-stdio-1.0+ hs-source-dirs: bench+ main-is: benchmarks.hs+ build-depends: base, clifford, criterion, numeric-prelude, stream-fusion+ -- -optlo-O3+ ghc-options: -fllvm -fexcess-precision -optlc-O=3 -O3+ default-language: Haskell2010
+ src/Numeric/Clifford/Internal.hs view
@@ -0,0 +1,42 @@+{-# OPTIONS_GHC -fllvm -fexcess-precision -optlo-O3 -O3 -optlc-O=3 -Wall #-}+{-# LANGUAGE TypeOperators, TypeFamilies,CPP #-}+module Numeric.Clifford.Internal (myTrace, trie, untrie, enumerate) where+import Numeric.Natural+import Prelude hiding (head,tail, null)+import Data.MemoTrie+import Data.List.Stream+import Control.Arrow+import Data.Bits+import qualified Debug.Trace as DebugTrace+#ifdef DEBUG+myTrace = DebugTrace.trace+#else+myTrace _ x = x+#endif+instance HasTrie Natural where+ newtype Natural :->: a = NaturalTrie ((Bool,[Bool]) :->: a)+ trie f = NaturalTrie (trie (f . unbitsZ)) + untrie (NaturalTrie t) = untrie t . bitsZ+ enumerate (NaturalTrie t) = enum' unbitsZ t+++unbitsZ :: (Prelude.Num n, Bits n) => (Bool,[Bool]) -> n+unbitsZ (headder,bs) = (unbits (headder:bs))++bitsZ :: (Prelude.Num n, Ord n, Bits n) => n -> (Bool,[Bool])+bitsZ i = (h, t ) where + theBits = bits i+ (h,t) = if null theBits+ then (False,[])+ else (head theBits, tail theBits)+bits :: (Prelude.Num t, Bits t) => t -> [Bool]+bits 0 = []+bits x = testBit x 0 : bits (shiftR x 1)+unbits :: (Prelude.Num t, Bits t) => [Bool] -> t+unbits [] = 0+unbits (x:xs) = unbit x .|. shiftL (unbits xs) 1+unbit :: Prelude.Num t => Bool -> t+unbit False = 0+unbit True = 1+enum' :: (HasTrie a) => (a -> a') -> (a :->: b) -> [(a', b)]+enum' f = (fmap.first) f . enumerate
src/Numeric/Clifford/LinearOperators.lhs view
@@ -1,4 +1,30 @@ \begin{code}+{-# LANGUAGE NoImplicitPrelude, RankNTypes #-} module Numeric.Clifford.LinearOperators where+import NumericPrelude+import Numeric.Clifford.Multivector+import Algebra.Algebraic+import GHC.TypeLits+\end{code}+What is a linear operator? Just a Vector -> Vector!++\begin{code}+type LinearOperator p q f = Multivector p q f -> Multivector p q f+type LinearOperatorCreator p q f = (Algebra.Algebraic.C f, Ord f, SingI p, SingI q) => Multivector p q f -> LinearOperator p q f+reflect u x = (-u)*x*recip u++makeReflectionOperator ::LinearOperatorCreator p q f+makeReflectionOperator u = reflect u++rotate spinor x = (reverseMultivector spinor) * x * spinor+rotatePlaneAngle plane angle = rotate (Numeric.Clifford.Multivector.exp ((normalised plane) * (angle/2)))++makeRotationOperator :: LinearOperatorCreator p q f+makeRotationOperator u = rotate u++project u x = inverse u * (u `dot` x)++makeProjectionOperator :: LinearOperatorCreator p q f+makeProjectionOperator u = project u \end{code}
src/Numeric/Clifford/Multivector.lhs view
@@ -63,7 +63,7 @@ import Test.QuickCheck import Math.Sequence.Converge (convergeBy) import Control.DeepSeq -import Number.Ratio hiding (scale)+import Number.Ratio hiding (scale, recip) import Algebra.ToRational import qualified GHC.Num as PNum import Control.Lens hiding (indices)@@ -75,7 +75,7 @@ import GHC.TypeLits import Control.Lens.Lens import Data.Word-import Control.Applicative+import Control.Applicative ((<$>)) import Numeric.Clifford.Internal @@ -240,14 +240,20 @@ -(/) :: (Algebra.Field.C f, Ord f, SingI p, SingI q) => Multivector p q f -> f -> Multivector p q f-(/) v d = BladeSum $ map (bladeScaleLeft (NPN.recip d)) $ mvTerms v --Algebra.Field.recip d *> v+--(/) :: (Algebra.Field.C f, Ord f, SingI p, SingI q) => Multivector p q f -> f -> Multivector p q f+--(/) v d = BladeSum $ map (bladeScaleLeft (NPN.recip d)) $ mvTerms v --Algebra.Field.recip d *> v (</) n d = Numeric.Clifford.Multivector.inverse d * n (/>) n d = n * Numeric.Clifford.Multivector.inverse d (</>) n d = n /> d -integratePoly c x = c : zipWith (Numeric.Clifford.Multivector./) x progression+{-#INLINE scaleLeft #-}+scaleLeft s v = BladeSum $ map (bladeScaleLeft s) $ mvTerms v+{-#INLINE scaleRight #-}+scaleRight v s = BladeSum $ map (bladeScaleRight s) $ mvTerms v+{-#INLINE divideRight #-}+divideRight v s = scaleRight v (recip s)+--integratePoly c x = c : zipWith (Numeric.Clifford.Multivector./) x progression --converge :: (Eq f, Show f) => [f] -> f converge [] = error "converge: empty list"@@ -291,7 +297,7 @@ expScaled = converge $ shanksTransformation.shanksTransformation . compensatedRunningSum $ expTerms scaled convergeTerms terms = converge $ shanksTransformation.shanksTransformation.compensatedRunningSum $ expTerms terms mag = myTrace ("In exponential, magnitude is " ++ show ( magnitude x)) magnitude x- scaled = let val = (Numeric.Clifford.Multivector./) x mag in myTrace ("In exponential, scaled is" ++ show val) val+ scaled = let val = (recip mag) *> x in myTrace ("In exponential, scaled is" ++ show val) val @@ -314,7 +320,8 @@ cos x = converge $ shanksTransformation $ compensatedRunningSum (Algebra.Ring.one : cosTerms x) cosTerms x = seriesMinusPlus $ takeEvery 2 $ tail $ expTerms x -expTerms x = map snd $ iterate (\(n,b) -> (n + 1, (x*b) Numeric.Clifford.Multivector./ fromInteger n )) (1::NPN.Integer,one)+expTerms :: (Algebra.Algebraic.C f, SingI p, SingI q, Ord f) => Multivector p q f -> [Multivector p q f]+expTerms x = map snd $ iterate (\(n,b) -> (n + 1, (recip $ fromInteger n ) `scaleLeft` (x*b) )) (1::NPN.Integer,one) dot a@(BladeSum _) b@(BladeSum _) = mvNormalForm $ BladeSum [x `bDot` y | x <- mvTerms a, y <- mvTerms b] wedge a@(BladeSum _) b@(BladeSum _) = mvNormalForm $ BladeSum [x `bWedge` y | x <- mvTerms a, y <- mvTerms b]@@ -325,9 +332,14 @@ reverseBlade b = bladeNormalForm $ b & indices %~ reverse reverseMultivector v = mvNormalForm $ v & terms.traverse%~ reverseBlade -inverse a = assert (a /= zero) $ reverseMultivector a Numeric.Clifford.Multivector./ bScale (head $ mvTerms (a * reverseMultivector a))-recip=Numeric.Clifford.Multivector.inverse +inverse a@(BladeSum _) = assert (a /= zero) $ (recip scalarComponent) *> (reverseMultivector a) where+ scalarComponent = bScale (head $ mvTerms (a * reverseMultivector a))+++instance (Algebra.Field.C f, Ord f, SingI p, SingI q) => Algebra.Field.C (Multivector p q f) where+ recip = inverse+ instance (Algebra.Field.C f, Ord f, SingI p, SingI q) => Algebra.OccasionallyScalar.C f (Multivector p q f) where toScalar = bScale . bladeGetGrade 0 . head . mvTerms toMaybeScalar (BladeSum [Blade s []]) = Just s@@ -412,8 +424,10 @@ \begin{code} -normalised a = a * (scalar $ NPN.recip $ magnitude a)+normalised :: (Ord f, Algebra.Algebraic.C f, SingI p, SingI q) => Multivector p q f -> Multivector p q f+normalised a = a `scaleRight` ( NPN.recip $ magnitude a) +log :: (Algebra.Transcendental.C f, Ord f, SingI p, SingI q, Show f) => Multivector p q f -> Multivector p q f log (BladeSum [Blade s []]) = scalar $ NPN.log s log a = scalar (NPN.log mag) + log' scaled where scaled = normalised a
src/Numeric/Clifford/NumericIntegration.lhs view
@@ -169,7 +169,7 @@ (time + (stepSizeAdapter time state)*h*(c s), newState) where zi :: Int -> [Multivector p q t] zi i = (\out -> myTrace ("initialGuess is " ++ show initialGuess++" whereas the final one is " ++ show out) out) $- assert (i <= s && i>= 1) $ converger $ iterate (zkp1 i) initialGuess where+ converger $ iterate (zkp1 i) initialGuess where initialGuess :: [Multivector p q t] initialGuess = if i == 1 || null (zi (i-1)) then map (h'*>) $ unproject $ f guessTime state else zi (i-1) adaptiveStepSizeFraction :: t
src/Numeric/Clifford/NumericIntegration/DefaultIntegrators.hs view
@@ -62,13 +62,13 @@ rk4Classical :: (Ord a, Algebra.Algebraic.C a, SingI p, SingI q) => stateType -> a -> (stateType->stateType) -> ([Multivector p q a] -> stateType) -> (stateType -> [Multivector p q a]) -> stateType rk4Classical state h f project unproject = project newState where- update = map (\(k1', k2', k3', k4') -> sumList [k1',2*k2',2*k3',k4'] MV./ Algebra.Ring.fromInteger 6) $ zip4 k1 k2 k3 k4+ update = map (\(k1', k2', k3', k4') -> sumList [k1',2*k2',2*k3',k4'] `divideRight` Algebra.Ring.fromInteger 6) $ zip4 k1 k2 k3 k4 newState = zipWith (+) state' update state' = unproject state evalDerivatives x = unproject $ f $ project x k1 = map (h*>) $ evalDerivatives state'- k2 = map (h*>) $ evalDerivatives . map (uncurry (+)) $ zip state' (map (MV./ two) k1)- k3 = map (h*>) $ evalDerivatives . map (uncurry (+)) $ zip state' (map (MV./ two) k2)+ k2 = map (h*>) $ evalDerivatives . map (uncurry (+)) $ zip state' (map (`divideRight` two) k1)+ k3 = map (h*>) $ evalDerivatives . map (uncurry (+)) $ zip state' (map (`divideRight` two) k2) k4 = map (h*>) $ evalDerivatives . map (uncurry (+)) $ zip state' k3 rk4ClassicalList state h f = rk4Classical state h f id id