diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -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 diff --git a/clifford.cabal b/clifford.cabal --- a/clifford.cabal +++ b/clifford.cabal @@ -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 diff --git a/src/Numeric/Clifford/Blade.lhs b/src/Numeric/Clifford/Blade.lhs --- a/src/Numeric/Clifford/Blade.lhs +++ b/src/Numeric/Clifford/Blade.lhs @@ -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 -> diff --git a/src/Numeric/Clifford/LinearOperators.lhs b/src/Numeric/Clifford/LinearOperators.lhs new file mode 100644 --- /dev/null +++ b/src/Numeric/Clifford/LinearOperators.lhs @@ -0,0 +1,4 @@ +\begin{code} +module Numeric.Clifford.LinearOperators where + +\end{code} diff --git a/src/Numeric/Clifford/Multivector.lhs b/src/Numeric/Clifford/Multivector.lhs --- a/src/Numeric/Clifford/Multivector.lhs +++ b/src/Numeric/Clifford/Multivector.lhs @@ -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} - diff --git a/src/Numeric/Clifford/NumericIntegration.lhs b/src/Numeric/Clifford/NumericIntegration.lhs --- a/src/Numeric/Clifford/NumericIntegration.lhs +++ b/src/Numeric/Clifford/NumericIntegration.lhs @@ -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) diff --git a/test/Numeric/Clifford/BladeSpec.lhs b/test/Numeric/Clifford/BladeSpec.lhs --- a/test/Numeric/Clifford/BladeSpec.lhs +++ b/test/Numeric/Clifford/BladeSpec.lhs @@ -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