clifford 0.1.0.5 → 0.1.0.6
raw patch · 7 files changed
+62/−47 lines, 7 filesdep +MemoTrie
Dependencies added: MemoTrie
Files
- changelog.md +1/−0
- clifford.cabal +6/−5
- src/Numeric/Clifford/Blade.lhs +31/−22
- src/Numeric/Clifford/LinearOperators.lhs +4/−0
- src/Numeric/Clifford/Multivector.lhs +13/−13
- src/Numeric/Clifford/NumericIntegration.lhs +6/−7
- test/Numeric/Clifford/BladeSpec.lhs +1/−0
changelog.md view
@@ -1,4 +1,5 @@ -*-change-log-*-+ 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 0.1.0.4 Made multivectors have a (p,q) metric signature at the type level
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.5+version: 0.1.0.6 -- A short (one-line) description of the package. synopsis: A Clifford algebra library@@ -43,14 +43,15 @@ -- Extra files to be distributed with the package, such as examples or a -- README.-extra-source-files: README.md changelog.md test/Numeric/Clifford/BladeSpec.lhs test/Numeric/Clifford/MultivectorSpec.lhs+extra-source-files: README.md changelog.md test/Numeric/Clifford/BladeSpec.lhs test/Numeric/Clifford/MultivectorSpec.lhs +--test/Numeric/Clifford/LinearOperatorsSpec.lhs -- Constraint on the version of Cabal needed to build this package. cabal-version: >=1.10 library -- Modules exported by the library.- exposed-modules: Numeric.Clifford.Blade, Numeric.Clifford.Multivector, Numeric.Clifford.NumericIntegration, Numeric.Clifford.NumericIntegration.DefaultIntegrators, Numeric.Clifford.ClassicalMechanics+ exposed-modules: Numeric.Clifford.Blade, Numeric.Clifford.Multivector, Numeric.Clifford.NumericIntegration, Numeric.Clifford.NumericIntegration.DefaultIntegrators, Numeric.Clifford.ClassicalMechanics, Numeric.Clifford.LinearOperators -- Modules included in this library but not exported. -- other-modules: @@ -60,7 +61,7 @@ -- Other library packages from which modules are imported. build-depends: base >=4.6 && <4.9, numeric-prelude >= 0.4.0.1 && < 0.5.0, permutation >= 0.4.1 && < 0.5, 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+ 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 -- Directories containing source files. hs-source-dirs: src@@ -72,7 +73,7 @@ type: exitcode-stdio-1.0 default-extensions: DataKinds, ScopedTypeVariables default-language: Haskell2010- ghc-options: -Wall -Werror+ ghc-options: -Wall hs-source-dirs: test main-is: Spec.lhs build-depends: base, clifford, hspec, numeric-prelude, QuickCheck, nats
src/Numeric/Clifford/Blade.lhs view
@@ -22,7 +22,7 @@ {-# LANGUAGE NoMonomorphismRestriction, UnicodeSyntax, GADTs #-} {-# LANGUAGE FlexibleInstances, UnicodeSyntax, GADTs, KindSignatures, DataKinds #-} {-# LANGUAGE TemplateHaskell, StandaloneDeriving, TypeOperators #-}-{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE MultiParamTypeClasses, TypeFamilies #-} \end{code} %if False \begin{code}@@ -53,9 +53,8 @@ import GHC.TypeLits hiding (isEven, isOdd) import GHC.Real (fromIntegral, toInteger) import Algebra.Field-import Debug.Trace---trace _ a = a-+import Data.MemoTrie+import Numeric.Clifford.Internal \end{code} @@ -118,26 +117,35 @@ \begin{code}+ bladeNormalForm :: forall (p::Nat) (q::Nat) f. Blade p q f -> Blade p q f bladeNormalForm (Blade scale indices) = result where- result = if (any (\i -> (GHC.Real.toInteger i) >= d) indices) then trace "Blade contains vector with i >= d" zeroBlade else Blade scale' uniqueSorted+ result = if (any (\i -> (GHC.Real.toInteger i) >= d) indices) then zeroBlade else Blade scale' newIndices p' = (fromSing (sing :: Sing p)) :: Integer q' = (fromSing (sing :: Sing q)) :: Integer d = p' + q'- numOfIndices = length indices- (sorted, perm) = Data.Permute.sort numOfIndices indices- scale' = if (isEven perm) /= (negated) then scale else negate scale- (uniqueSorted,negated) = removeDupPairs [] sorted False- where- removeDupPairs :: [Natural] -> [Natural] -> Bool -> ([Natural],Bool)- removeDupPairs accum [] negated = (accum,negated)- removeDupPairs accum [x] negated = (accum++[x],negated)- removeDupPairs accum (x:y:rest) negated | x == y = - if GHC.Real.toInteger x < q' - then removeDupPairs accum rest (not negated)- else removeDupPairs accum rest negated- | otherwise = removeDupPairs (accum++[x]) (y:rest) negated+ + + scale' = if doNotNegate then scale else negate scale+ (newIndices, doNotNegate) = sortIndices (indices,q')++sortIndices = memo sortIndices' where+sortIndices' :: ([Natural],Integer) -> ([Natural],Bool) +sortIndices' (indices,q') = (uniqueSorted, doNotNegate) where+ (sorted, perm) = Data.Permute.sort numOfIndices indices+ numOfIndices = length indices+ doNotNegate = (isEven perm) /= (negated)+ (uniqueSorted,negated) = removeDupPairs [] sorted False+ where+ removeDupPairs :: [Natural] -> [Natural] -> Bool -> ([Natural],Bool)+ removeDupPairs accum [] negated = (accum,negated)+ removeDupPairs accum [x] negated = (accum++[x],negated)+ removeDupPairs accum (x:y:rest) negated | x == y = + if GHC.Real.toInteger x < q' + then removeDupPairs accum rest (not negated)+ else removeDupPairs accum rest negated+ | otherwise = removeDupPairs (accum++[x]) (y:rest) negated \end{code} What is the grade of a blade? Just the number of indices.@@ -160,7 +168,7 @@ \begin{code} bladeMul :: Blade p q f -> Blade p q f-> Blade p q f-bladeMul x@(Blade _ _) y@(Blade _ _)= bladeNormalForm $ Blade (bScale x Algebra.Ring.* bScale y) (bIndices x ++ bIndices y) -- HAVE TO MAKE Q INDICES SQUARE NEGATIVE+bladeMul x@(Blade _ _) y@(Blade _ _)= bladeNormalForm $ Blade (bScale x Algebra.Ring.* bScale y) (bIndices x ++ bIndices y) multiplyBladeList :: (SingI p, SingI q, Algebra.Field.C f) => [Blade p q f] -> Blade p q f multiplyBladeList [] = error "Empty blade list!" multiplyBladeList (a:[]) = a@@ -203,11 +211,12 @@ \begin{code} instance (Algebra.Additive.C f, Ord f) => Ord (Blade p q f) where compare a b | bIndices a == bIndices b = compare (bScale a) (bScale b)- | otherwise = case compare ((length . bIndices) a) ((length . bIndices) b) of+ | otherwise = compareIndices (bIndices a) (bIndices b)+compareIndices = memo compareIndices' where+ compareIndices' a b = case compare (length a) (length b) of LT -> LT GT -> GT- EQ -> compare (bIndices a) (bIndices b)-+ EQ -> compare a b instance Arbitrary Natural where arbitrary = sized $ \n ->
+ src/Numeric/Clifford/LinearOperators.lhs view
@@ -0,0 +1,4 @@+\begin{code}+module Numeric.Clifford.LinearOperators where++\end{code}
src/Numeric/Clifford/Multivector.lhs view
@@ -76,9 +76,9 @@ import Control.Lens.Lens import Data.Word import Control.Applicative-import Debug.Trace---trace _ a = a+import Numeric.Clifford.Internal + \end{code} @@ -255,7 +255,7 @@ where empty = error "converge: error in implmentation" checkPeriodic (a:b:c:_)- | (trace ("Converging at " ++ show a) a) == b = Just a+ | (myTrace ("Converging at " ++ show a) a) == b = Just a | a == c = Just a checkPeriodic _ = Nothing @@ -277,21 +277,21 @@ shanksTransformation (xnm1:xn:xnp1:xs) | xnm1 == xn = [xn] | xnm1 == xnp1 = [xnm1] | denominator == zero = [xnp1]- | otherwise = trace ("Shanks transformation input = " ++ show xn ++ "\nShanks transformation output = " ++ show out) out:shanksTransformation (xn:xnp1:xs) where+ | otherwise = myTrace ("Shanks transformation input = " ++ show xn ++ "\nShanks transformation output = " ++ show out) out:shanksTransformation (xn:xnp1:xs) where out = numerator /> denominator numerator = sumList [xnp1*xnm1, negate (xn^2)] denominator = sumList [xnp1, (-2)*xn, xnm1] --exp ::(Ord f, Show f, Algebra.Transcendental.C f)=> Multivector f -> Multivector f-exp (BladeSum [ Blade s []]) = trace ("scalar exponential of " ++ show s) scalar $ Algebra.Transcendental.exp s-exp x = trace ("Computing exponential of " ++ show x) convergeTerms x where --(expMag ^ expScaled) where+exp (BladeSum [ Blade s []]) = myTrace ("scalar exponential of " ++ show s) scalar $ Algebra.Transcendental.exp s+exp x = myTrace ("Computing exponential of " ++ show x) convergeTerms x where --(expMag ^ expScaled) where --todo: compute a ^ p via a^n where n = floor p then multiply remaining power expMag = Algebra.Transcendental.exp mag expScaled = converge $ shanksTransformation.shanksTransformation . compensatedRunningSum $ expTerms scaled convergeTerms terms = converge $ shanksTransformation.shanksTransformation.compensatedRunningSum $ expTerms terms- mag = trace ("In exponential, magnitude is " ++ show ( magnitude x)) magnitude x- scaled = let val = (Numeric.Clifford.Multivector./) x mag in trace ("In exponential, scaled is" ++ show val) val+ 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 @@ -316,11 +316,12 @@ expTerms x = map snd $ iterate (\(n,b) -> (n + 1, (x*b) Numeric.Clifford.Multivector./ fromInteger n )) (1::NPN.Integer,one) -dot a b = mvNormalForm $ BladeSum [x `bDot` y | x <- mvTerms a, y <- mvTerms b]-wedge a b = mvNormalForm $ BladeSum [x `bWedge` y | x <- mvTerms a, y <- mvTerms b]-(∧) = wedge-(⋅) = dot+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] +(∧) = wedge :: Multivector p q f -> Multivector p q f -> Multivector p q f+(⋅) = dot :: Multivector p q f -> Multivector p q f -> Multivector p q f+ reverseBlade b = bladeNormalForm $ b & indices %~ reverse reverseMultivector v = mvNormalForm $ v & terms.traverse%~ reverseBlade @@ -426,7 +427,6 @@ Now let's do (slow as fuck probably) numerical integration! :D~! Since this is gonna be used for physical applications, it's we're gonna start off with a Hamiltonian structure and then a symplectic integrator. \begin{code}-
src/Numeric/Clifford/NumericIntegration.lhs view
@@ -54,9 +54,8 @@ import Control.Exception (assert) import Data.Maybe import GHC.TypeLits+import Numeric.Clifford.Internal import Data.DeriveTH-import Debug.Trace---trace _ a = a elementAdd = zipWith (+)@@ -95,7 +94,7 @@ convergeList = converge -showOutput name x = trace ("output of " ++ name ++" is " ++ show x) x+showOutput name x = myTrace ("output of " ++ name ++" is " ++ show x) x convergeTolLists :: (Ord f, Algebra.Absolute.C f, Algebra.Algebraic.C f, Show f, SingI p, SingI q) => f -> [[Multivector p q f]] -> [Multivector p q f]@@ -104,7 +103,7 @@ where empty = error "converge: error in impl" check (a:b:c:_)- | (trace ("Converging at " ++ show a) a) == b = Just b+ | (myTrace ("Converging at " ++ show a) a) == b = Just b | a == c = Just c | ((showOutput ("convergence check with tolerance " ++ show t) $ magnitude (sumList $ (zipWith (\x y -> NPN.abs (x-y)) b c))) <= t) = showOutput ("convergence with tolerance "++ show t )$ Just c@@ -157,8 +156,8 @@ converger :: [[Multivector p q t]] -> [Multivector p q t] converger = case find (\x -> isConvergenceTolerance x || isConvergenceFunction x) attributes of Just (ConvergenceFunction conv) -> conv- Just (ConvergenceTolerance tol) -> convergeTolLists (trace ("Convergence tolerance set to " ++ show tol)tol)- Nothing -> trace "No convergence tolerance specified, defaulting to equality" convergeList + Just (ConvergenceTolerance tol) -> convergeTolLists (myTrace ("Convergence tolerance set to " ++ show tol)tol)+ Nothing -> myTrace "No convergence tolerance specified, defaulting to equality" convergeList stepSizeAdapter :: AdaptiveStepSizeFunction t stateType stepSizeAdapter = case find isAdaptiveStepSize attributes of@@ -169,7 +168,7 @@ rkMethodImplicitFixedPoint h f project unproject (time, state) = (time + (stepSizeAdapter time state)*h*(c s), newState) where zi :: Int -> [Multivector p q t]- zi i = (\out -> trace ("initialGuess is " ++ show initialGuess++" whereas the final one is " ++ show out) out) $+ 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 initialGuess :: [Multivector p q t] initialGuess = if i == 1 || null (zi (i-1)) then map (h'*>) $ unproject $ f guessTime state else zi (i-1)
test/Numeric/Clifford/BladeSpec.lhs view
@@ -18,6 +18,7 @@ let d = Blade 4.0 [2] :: Blade 3 1 Double let e = Blade 2.0 [0] :: Blade 3 1 Double let f = Blade 3.0 [0] :: Blade 3 1 Double+ --todo: Memoise the blade index sorting function. describe "bladeMul" $ do it "multiplies an n-blade with an m-blade to give an n+m blade if each index is unique" $ do (a `bladeMul` b) `shouldBe` ab