clifford 0.1.0.4 → 0.1.0.5
raw patch · 8 files changed
+162/−42 lines, 8 filesdep +clifforddep +hspecdep ~basedep ~numeric-prelude
Dependencies added: clifford, hspec
Dependency ranges changed: base, numeric-prelude
Files
- changelog.md +1/−0
- clifford.cabal +12/−4
- src/Numeric/Clifford/Blade.lhs +28/−7
- src/Numeric/Clifford/ClassicalMechanics.lhs +10/−22
- src/Numeric/Clifford/Multivector.lhs +32/−9
- test/Numeric/Clifford/BladeSpec.lhs +34/−0
- test/Numeric/Clifford/MultivectorSpec.lhs +42/−0
- test/Spec.lhs +3/−0
changelog.md view
@@ -1,3 +1,4 @@ -*-change-log-*-+ 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.4+version: 0.1.0.5 -- A short (one-line) description of the package. synopsis: A Clifford algebra library@@ -43,8 +43,7 @@ -- Extra files to be distributed with the package, such as examples or a -- README.-extra-source-files: README.md changelog.md-+extra-source-files: README.md changelog.md test/Numeric/Clifford/BladeSpec.lhs test/Numeric/Clifford/MultivectorSpec.lhs -- Constraint on the version of Cabal needed to build this package. cabal-version: >=1.10 @@ -61,7 +60,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+ 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 -- Directories containing source files. hs-source-dirs: src@@ -69,3 +68,12 @@ -- Base language which the package is written in. default-language: Haskell2010 +test-suite spec+ type: exitcode-stdio-1.0+ default-extensions: DataKinds, ScopedTypeVariables+ default-language: Haskell2010+ ghc-options: -Wall -Werror+ hs-source-dirs: test+ main-is: Spec.lhs+ build-depends: base, clifford, hspec, numeric-prelude, QuickCheck, nats+
src/Numeric/Clifford/Blade.lhs view
@@ -21,7 +21,7 @@ {-# LANGUAGE NoImplicitPrelude, FlexibleContexts, RankNTypes, ScopedTypeVariables, DeriveDataTypeable #-} {-# LANGUAGE NoMonomorphismRestriction, UnicodeSyntax, GADTs #-} {-# LANGUAGE FlexibleInstances, UnicodeSyntax, GADTs, KindSignatures, DataKinds #-}-{-# LANGUAGE TemplateHaskell, StandaloneDeriving #-}+{-# LANGUAGE TemplateHaskell, StandaloneDeriving, TypeOperators #-} {-# LANGUAGE MultiParamTypeClasses #-} \end{code} %if False@@ -39,6 +39,7 @@ import Algebra.Laws import Algebra.Absolute import Algebra.Additive+import Control.DeepSeq import Algebra.Ring import Data.Serialize import Data.Word@@ -66,7 +67,6 @@ data Blade (p :: Nat) (q :: Nat) f where Blade :: forall p q f . (SingI p, SingI q, Algebra.Field.C f) => {_scale :: f, _indices :: [Natural]} -> Blade p q f --- makeLenses ''Blade scale :: Lens' (Blade p q f) f scale = lens _scale (\blade v -> blade {_scale = v}) indices :: Lens' (Blade p q f) [Natural]@@ -75,6 +75,7 @@ dimension _ = (toNatural ((GHC.Real.fromIntegral $ fromSing (sing :: Sing p))::Word),toNatural((GHC.Real.fromIntegral $ fromSing (sing :: Sing q))::Word)) bScale b = b^.scale bIndices b = b^.indices+instance (Control.DeepSeq.NFData f) => Control.DeepSeq.NFData (Blade p q f) instance(Show f) => Show (Blade p q f) where --TODO: Do this with HaTeX show (Blade scale indices) = pref ++ if null indices then "" else basis where@@ -117,13 +118,13 @@ \begin{code}-bladeNormalForm :: forall (p::Nat) (q::Nat) f. Blade p q f -> Blade p q f+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 trace "Blade contains vector with i >= d" zeroBlade else Blade scale' uniqueSorted p' = (fromSing (sing :: Sing p)) :: Integer q' = (fromSing (sing :: Sing q)) :: Integer- d = p' + q' + d = p' + q' numOfIndices = length indices (sorted, perm) = Data.Permute.sort numOfIndices indices scale' = if (isEven perm) /= (negated) then scale else negate scale@@ -133,7 +134,7 @@ 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 > p' + if GHC.Real.toInteger x < q' then removeDupPairs accum rest (not negated) else removeDupPairs accum rest negated | otherwise = removeDupPairs (accum++[x]) (y:rest) negated@@ -202,7 +203,10 @@ \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 = compare (bIndices a) (bIndices b)+ | otherwise = case compare ((length . bIndices) a) ((length . bIndices) b) of+ LT -> LT+ GT -> GT+ EQ -> compare (bIndices a) (bIndices b) instance Arbitrary Natural where@@ -210,6 +214,23 @@ let n' = NPN.abs n in fmap (toNatural . (\x -> (GHC.Real.fromIntegral x)::Word)) (choose (0, n')) shrink = shrinkIntegral++instance (SingI p, SingI q, Algebra.Field.C f, Arbitrary f) => Arbitrary (Blade p q f) where+ arbitrary = do+ let p'' = (fromSing (sing :: Sing p)) :: Integer+ let q'' = (fromSing (sing :: Sing q)) + let d = p'' + q''+ let maxLength = (2^d - 1) :: Integer+ scale <- arbitrary+ listSize <- choose (0, maxLength)+ indices <- vectorOf (NPN.fromIntegral listSize) (resize (NPN.fromIntegral d-1) arbitrary )+ return (Blade scale indices) + where+ + + + + -- $(derive makeArbitrary ''Blade) \end{code}
src/Numeric/Clifford/ClassicalMechanics.lhs view
@@ -86,22 +86,17 @@ -a `cross` b = (negate $ one)`e`[1,2,3] * (a ∧ b)-data PhysicalVector (p::Nat) (q::Nat) t = PhysicalVector {dimension :: Natural, r :: Multivector p q t, referenceFrame :: ReferenceFrame p q t}-{-squishToDimension (PhysicalVector d (BladeSum terms) f) = PhysicalVector d r' f where- r' = BladeSum terms' where- terms' = terms & filter (\(Blade _ ind) -> all (\k -> k <= d) ind)-squishToDimension' d (BladeSum terms) = r' where- r' = BladeSum terms' where- terms' = terms & filter (\(Blade _ ind) -> all (\k -> k <= d) ind)-}+a `cross` b = (negate $ one)`e`[0,1,2] * (a ∧ b)+data PhysicalVector (p::Nat) (q::Nat) t = PhysicalVector {r :: Multivector p q t, referenceFrame :: ReferenceFrame p q t} + data RigidBody (p::Nat) (q::Nat) f where RigidBody:: (Algebra.Field.C f, Algebra.Module.C f (Multivector p q f)) => {position :: PhysicalVector p q f,- _momentum :: PhysicalVector p q f,- _mass :: f,- _attitude :: PhysicalVector p q f,- _angularMomentum :: PhysicalVector p q f,- _inertia :: PhysicalVector p q f+ momentum :: PhysicalVector p q f,+ mass :: f,+ attitude :: PhysicalVector p q f,+ angularMomentum :: PhysicalVector p q f,+ inertia :: PhysicalVector p q f } -> RigidBody p q f --makeLenses ''RigidBody doesn't actually work@@ -110,15 +105,8 @@ 5. figure a way to take exterior product of 1 forms at a type level so i can just go like: omega = df1 ^ df2 ^ df ; omega a b c -} -{-data NDVector (n :: Nat) f where- NDVector :: (Algebra.Field.C f, Algebra.Module.C f (Multivector f)) => {value :: Multivector f} -> NDVector n f-}--{-ndVector :: forall n.(n ~ Nat) => Proxy n -> (forall f.- (Algebra.Field.C f, Algebra.Module.C f (Multivector f)) =>- Multivector f -> NDVector (n) f)-ndVector _ value = NDVector $ squishToDimension' (toNatural nummed) value where- nummed :: Word32- nummed = fromIntegral $ fromSing (sing :: Sing n)-}+type Vector3 f = Multivector 3 0 f+type STVector f = Multivector 3 1 f \end{code} \bibliographystyle{IEEEtran} \bibliography{biblio.bib}
src/Numeric/Clifford/Multivector.lhs view
@@ -21,7 +21,7 @@ {-# LANGUAGE NoImplicitPrelude, FlexibleContexts, RankNTypes, ScopedTypeVariables, DeriveDataTypeable #-} {-# LANGUAGE NoMonomorphismRestriction, UnicodeSyntax, GADTs#-} {-# LANGUAGE FlexibleInstances, StandaloneDeriving, KindSignatures, DataKinds #-}-{-# LANGUAGE TemplateHaskell, TypeOperators #-}+{-# LANGUAGE TemplateHaskell, TypeOperators, DeriveFunctor #-} {-# LANGUAGE MultiParamTypeClasses #-} \end{code} %if False@@ -73,7 +73,9 @@ import Data.Data import Data.DeriveTH import GHC.TypeLits+import Control.Lens.Lens import Data.Word+import Control.Applicative import Debug.Trace --trace _ a = a @@ -86,12 +88,21 @@ data Multivector (p::Nat) (q::Nat) f where BladeSum :: forall p q f . (Ord f, Algebra.Field.C f, SingI p, SingI q) => { _terms :: [Blade p q f]} -> Multivector p q f +instance (SingI p, SingI q, Algebra.Field.C f, Arbitrary f, Ord f) => Arbitrary (Multivector p q f) where+ arbitrary = mvNormalForm <$> BladeSum <$> (vector d) where+ p' = (fromSing (sing :: Sing p)) :: Integer+ q' = (fromSing (sing :: Sing q)) + d = fromIntegral (p' + q')+ deriving instance Eq (Multivector p q f)+--instance (SingI p, SingI q) => Functor (Multivector p q) where+-- fmap func x = func x--((terms x) & scale %~ func) deriving instance Ord (Multivector p q f) deriving instance (Show f) => Show (Multivector p q f)+--deriving instance (Read f) => Read (Multivector p q f) -dimension :: forall (p::Nat) (q::Nat) f. (SingI p, SingI q) => Multivector p q f -> (Natural,Natural)-dimension _ = (toNatural ((fromIntegral $ fromSing (sing :: Sing p))::Word),toNatural ((fromIntegral $ fromSing (sing :: Sing q))::Word))+signature :: forall (p::Nat) (q::Nat) f. (SingI p, SingI q) => Multivector p q f -> (Natural,Natural)+signature _ = (toNatural ((fromIntegral $ fromSing (sing :: Sing p))::Word),toNatural ((fromIntegral $ fromSing (sing :: Sing q))::Word)) terms :: Lens' (Multivector p q f) [Blade p q f] terms = lens _terms (\bladeSum v -> bladeSum {_terms = v})@@ -105,6 +116,9 @@ groupLikeTerms ::Eq f => [Blade p q f] -> [[Blade p q f]] groupLikeTerms = groupBy (\a b -> a^.indices == b^.indices) +compareTol :: (Algebra.Algebraic.C f, Algebra.Absolute.C f, Ord f, SingI p, SingI q) => Multivector p q f -> Multivector p q f -> f -> Bool+compareTol x y tol = ((NPN.abs $ magnitude (x-y) ) <= tol)+ compensatedSum' :: (Algebra.Additive.C f) => [f] -> f compensatedSum' xs = kahan zero zero xs where kahan s _ [] = s@@ -171,7 +185,7 @@ instance (Control.DeepSeq.NFData f) => Control.DeepSeq.NFData (Multivector p q f)-instance (Control.DeepSeq.NFData f) => Control.DeepSeq.NFData (Blade p q f)+ instance (Algebra.Field.C f, Ord f, SingI p, SingI q) => Algebra.Additive.C (Multivector p q f) where a + b = mvNormalForm $ BladeSum (mvTerms a ++ mvTerms b) a - b = mvNormalForm $ BladeSum (mvTerms a ++ map bladeNegate (mvTerms b))@@ -200,9 +214,11 @@ two = fromInteger 2 mul = (Algebra.Ring.*) -psuedoScalar :: forall (p::Nat) (q::Nat) f. (Ord f, Algebra.Field.C f, SingI p, SingI q, SingI (p+q)) => Multivector p q f-psuedoScalar = one `e` [1..(toNatural d)] where- d = fromIntegral (fromSing (sing :: Sing (p+q)) )::Word+psuedoScalar :: forall (p::Nat) (q::Nat) f. (Ord f, Algebra.Field.C f, SingI p, SingI q) => Multivector p q f+psuedoScalar = one `e` [0..(toNatural d)] where+ d = fromIntegral (p' + q' - 1 )::Word+ p'= fromSing (sing :: Sing p)+ q' = fromSing (sing :: Sing q) \end{code} @@ -338,12 +354,19 @@ \begin{code} root :: (Show f, Ord f, Algebra.Algebraic.C f, SingI p, SingI q) => NPN.Integer -> Multivector p q f -> Multivector p q f+root 0 _ = error "Cannot take 0th root"+root _ (BladeSum []) = error "Empty bladesum"+root _ (BladeSum [Blade zero []]) = error "Cannot compute a root of zero" root n (BladeSum [Blade s []]) = scalar $ Algebra.Algebraic.root n s-root n a@(BladeSum _) = converge $ rootIterationsStart n a one+root n a@(BladeSum _) = converge $ rootIterationsStart n a g where+ g = if q' <= 1 then one`e`[q',succ q'] else one + one `e` [0,1]+ (p',q') = signature a rootIterationsStart ::(Ord f, Show f, Algebra.Algebraic.C f)=> NPN.Integer -> Multivector p q f -> Multivector p q f -> [Multivector p q f] rootIterationsStart n a@(BladeSum (Blade s [] :xs)) one = rootHalleysIterations n a g where- g = if s >= NPN.zero then one else Algebra.Ring.one `e` [1,2] --BladeSum[Blade Algebra.Ring.one [1,2]]+ g = if s >= NPN.zero || q' == 1 then one else Algebra.Ring.one `e` [0,1] + (p',q') = signature a+ rootIterationsStart n a@(BladeSum _) g = rootHalleysIterations n a g
+ test/Numeric/Clifford/BladeSpec.lhs view
@@ -0,0 +1,34 @@+\begin{code}++module Numeric.Clifford.BladeSpec (main, spec) where++import Test.Hspec+import Test.QuickCheck+import Numeric.Clifford.Blade+import Algebra.Ring+main :: IO ()+main = hspec spec++spec :: Spec+spec = do+ let a = Blade 1.0 [0] :: Blade 3 1 Double+ let b = Blade 2.0 [1] :: Blade 3 1 Double+ let ab = Blade 2.0 [0,1] :: Blade 3 1 Double+ let c = Blade 3.0 [1,2] :: Blade 3 1 Double+ 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+ 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+ it "should handle duplicate blades which square to +1 due to metric signature" $ do+ c `bladeMul` d `shouldBe` Blade 12.0 [1]+ it "should handle diplicate blades which square to -1 due to metric signature" $ do+ e `bladeMul` f `shouldBe` Blade (-6.0) []+ context "leaves blades unchanged when multiplied by a scalar of 1" $ do+ it "on the left" $ property $+ \x -> (Blade one []) `bladeMul` x == (x::Blade 3 1 Double)+ it "on the right" $ property $ + \x -> x `bladeMul` (Blade one []) == (x::Blade 3 1 Double) ++\end{code}
+ test/Numeric/Clifford/MultivectorSpec.lhs view
@@ -0,0 +1,42 @@+\begin{code}++module Numeric.Clifford.MultivectorSpec (main, spec) where+import Prelude hiding ((^), (*))+import Test.Hspec+import Test.QuickCheck+import Numeric.Clifford.Multivector+import Algebra.Ring+--import Numeric.Natural+import Algebra.Additive (zero)+import Control.Exception (evaluate)+main :: IO ()+main = hspec spec++type STVector = Multivector 3 1 Double+spec :: Spec+spec = do+ let i = 1.0 `e` [1,2] :: STVector+ let fuckOffSized = (i + (scalar 3.8) + (1.1 `e` [0])) :: STVector+ let comp a b = compareTol a b 0.0000001+ describe "addition" $ do+ it "is assiocitive" $ property (\a (b::STVector) -> a + b == b + a)+ describe "multiplication" $ do+ it "should square unit bivectors to -1" $ do+ i*i `shouldBe` scalar (-1.0)+ describe "root n" $ do+ it "cannot compute the 0th root" $ do+ evaluate (root 0 i) `shouldThrow` anyErrorCall+ it "cannot compute a root of 0" $ do+ evaluate (root 1 (zero::STVector)) `shouldThrow` anyErrorCall+ it "computes the nth root of a value" $ do+ comp ((root 3 fuckOffSized)^3) fuckOffSized `shouldBe` True+ + {-it "computes the nth root of a vector. May fail to terminate." $ verbose prop where+ prop x k= (magnitude (abs ((rooted ^ n) - x))) <= 0.000001 || x == zero || n == zero where+ n :: Integer+ n = (fromIntegral (k::Natural)) `mod` 6+ rooted :: STVector+ rooted = root n x-}+++\end{code}
+ test/Spec.lhs view
@@ -0,0 +1,3 @@+\begin{code}+{-# OPTIONS_GHC -F -pgmF hspec-discover #-}+\end{code}