packages feed

clifford 0.1.0.13 → 0.1.0.14

raw patch · 10 files changed

+211/−102 lines, 10 filesdep +Chartdep +Chart-cairodep +colourdep ~basedep ~lensdep ~numeric-preludenew-component:exe:pendulum

Dependencies added: Chart, Chart-cairo, colour, data-default-class, gnuplot, reflection

Dependency ranges changed: base, lens, numeric-prelude, stream-fusion

Files

changelog.md view
@@ -1,4 +1,5 @@ -*-change-log-*-+	0.1.0.14 Added pendulum example 	0.1.0.13 Adding ability to create linear operators from matrices; made linear operators a category and monoid 	0.1.0.12 Adding type signatures to top level stuff 	0.1.0.11 More inlining and specialisation@@ -6,7 +7,7 @@ 	0.1.0.9 Inlined/specialised a bunch of function, hueg speed increase         0.1.0.8 Implemented algebraic/transcendental typeclasses 	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.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.13+version:             0.1.0.14  -- A short (one-line) description of the package. synopsis:            A Clifford algebra library@@ -57,11 +57,11 @@   -- other-modules:             -- LANGUAGE extensions used by modules in this package.-  other-extensions:    NoImplicitPrelude, ScopedTypeVariables, GADTs, DataKinds,KindSignatures,UnicodeSyntax,FlexibleContexts, RankNTypes,TemplateHaskell,NoMonomorphismRestriction,MultiParamTypeClasses,FlexibleInstances, TypeOperators+  other-extensions:    NoImplicitPrelude, ScopedTypeVariables, GADTs, DataKinds,KindSignatures,UnicodeSyntax,FlexibleContexts, RankNTypes,TemplateHaskell,NoMonomorphismRestriction,MultiParamTypeClasses,FlexibleInstances, TypeOperators, UnicodeSyntax, ConstraintKinds   -- 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, MemoTrie >= 0.6 && < 0.7, semigroupoids >= 4.0 && < 4.1, monoid-extras >=0.3 && <0.4+                       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, semigroupoids >= 4.0 && < 4.1, monoid-extras >=0.3 && <0.4, reflection >=1.4 && < 1.5      ghc-options: -fllvm -fexcess-precision -optlc-O=3 -O3   -- Directories containing source files.@@ -70,6 +70,15 @@   -- Base language which the package is written in.   default-language:    Haskell2010   ++executable pendulum+  build-depends: base, clifford, numeric-prelude,gnuplot,stream-fusion, Chart, Chart-cairo, colour, lens, data-default-class+  other-extensions: NoImplicitPrelude+  main-is: Pendulum.hs+  hs-source-dirs: examples+  default-language: Haskell2010+  ghc-options: -fllvm -fexcess-precision -optlc-O=3 -O3+	 test-suite spec   type: exitcode-stdio-1.0   default-extensions: DataKinds, ScopedTypeVariables@@ -86,5 +95,5 @@   main-is: benchmarks.hs   build-depends: base, clifford, criterion, numeric-prelude, stream-fusion   -- -optlo-O3-  ghc-options: -fllvm -fexcess-precision -optlc-O=3 -O3+  ghc-options: -fllvm -fexcess-precision -optlc-O=3 -O3    default-language: Haskell2010
+ examples/Pendulum.hs view
@@ -0,0 +1,53 @@+{-# LANGUAGE NoImplicitPrelude, NoMonomorphismRestriction, DataKinds #-}+import NumericPrelude (Double, fst, snd, ($), (.), seq)+import Prelude (getLine, putStrLn)+import Algebra.Transcendental+import Data.List.Stream+import Numeric.Clifford.Multivector+import Algebra.Ring +import Algebra.Additive +import Algebra.Field+import Numeric.Clifford.NumericIntegration+import Numeric.Clifford.NumericIntegration.DefaultIntegrators+import Numeric.Clifford.Blade++import Control.Lens+import Graphics.Rendering.Chart+import Data.Colour+import Data.Colour.Names+import Data.Default.Class+import Graphics.Rendering.Chart.Backend.Cairo++m = scalar 3 :: E3Vector+l = scalar 20 :: E3Vector+g = scalar 9.81 :: E3Vector++hamil _ [p,theta] =   [ (-m*g*l)* sin theta, p / (m*l*l)] ++integrator = gaussLegendreFourthOrder 0.1 hamil++++tenSeconds :: [(Double,Double,Double)]+tenSeconds = map ((\ (t, ([BladeSum [Blade a []],BladeSum [Blade b []]])) -> (t,a,b) ) )  $ take 5001 $ iterate integrator (0,[zero,one/10])+++chart = toRenderable layout+  where+    momentum = plot_lines_values .~ [ ( map (\(t,p,_) -> (t,p)) tenSeconds )]+                     $ plot_lines_style  . line_color .~ opaque blue+                     $ plot_lines_title .~ "momentum"+                     $ def++    angle = plot_lines_style . line_color .~ (opaque red)+                  $ plot_lines_values .~ [ ( map (\(t,_,theta) -> (t,theta)) tenSeconds )]+                  $ plot_lines_title .~ "angle"+                  $ def++    layout = layout_title .~ "Pendulum"+           $ layout_plots .~ [toPlot momentum,+                              toPlot angle]+           $ def++main = renderableToFile def  chart "pendulum.png"+
src/Numeric/Clifford/Blade.lhs view
@@ -63,6 +63,7 @@ \texttt{bScale} is the amplitude of the blade. \texttt{bIndices} are the indices for the basis.  \begin{code} + data Blade (p :: Nat) (q :: Nat) f where     Blade :: forall (p::Nat) (q::Nat) f . (Algebra.Field.C f, SingI p, SingI q) => {_scale :: f, _indices :: [Natural]} -> Blade p q f @@ -138,7 +139,7 @@ {-#INLINE bladeNormalForm#-} {-#SPECIALISE INLINE bladeNormalForm::E3Blade -> E3Blade #-} {-#SPECIALISE INLINE bladeNormalForm :: STBlade -> STBlade #-}-bladeNormalForm :: forall (p::Nat) (q::Nat) f.  Blade p q f -> Blade p q f+bladeNormalForm :: ∀ (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 zeroBlade else Blade scale' newIndices
src/Numeric/Clifford/ClassicalMechanics.lhs view
@@ -19,12 +19,10 @@ {-# LANGUAGE FlexibleInstances, TypeOperators #-} {-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE MultiParamTypeClasses #-}+{-# OPTIONS_HADDOCK show-extensions #-} \end{code} %if False \begin{code}-{-# OPTIONS_GHC -fllvm -fexcess-precision -optlo-O3 -O3 -optlc-O=3 -Wall #-}--- OPTIONS_GHC -Odph -fvectorise -package dph-lifted-vseg ---  LANGUAGE ParallelArrays \end{code} %endif @@ -34,7 +32,7 @@ import Numeric.Clifford.Blade import GHC.TypeLits import Data.Proxy-import NumericPrelude hiding (iterate, head, map, tail, reverse, scanl, zipWith, drop, (++), filter, null, length, foldr, foldl1, zip, foldl, concat, (!!), concatMap,any, repeat, replicate, elem, replicate, all)+import NumericPrelude hiding (iterate, head, map, tail, reverse, scanl, zipWith, drop, (++), filter, null, length, foldr, foldl1, zip, foldl, concat, (!!), concatMap,any, repeat, replicate, elem, replicate, all, (.) ) import Algebra.Absolute import Algebra.Algebraic import Algebra.Additive@@ -56,9 +54,15 @@ import Data.Maybe import Data.DeriveTH import Data.Word-import Debug.Trace---trace _ a = a+import Numeric.Clifford.Internal+import Numeric.Clifford.LinearOperators+import Control.Applicative+import Data.Monoid+import Control.Category +nonEqualFrames = "Non-equal reference frames! Insert code here to translate between them! :) Should really make reference frames as types and then have type families to convert between them :v :v :v"++ data EnergyMethod (p::Nat) (q::Nat) f = Hamiltonian{ _dqs :: [DynamicSystem p q f -> Multivector p q f], _dps :: [DynamicSystem p q f -> Multivector p q f]}  data DynamicSystem (p::Nat) (q::Nat) f = DynamicSystem {_time :: f, coordinates :: [Multivector p q f], _momenta :: [Multivector p q f], _energyFunction :: EnergyMethod p q f, _projector :: DynamicSystem p q f -> DynamicSystem p q f}@@ -80,23 +84,49 @@  Now to make a physical object. \begin{code}-data ReferenceFrame (p::Nat) (q::Nat) t = ReferenceFrame {basisFrame :: [Multivector p q t]}-psuedoScalar' :: forall f (p::Nat) (q::Nat). (Ord f, Algebra.Field.C f, SingI p, SingI q) => ReferenceFrame p q f -> Multivector p q f-psuedoScalar'  = multiplyList . basisFrame+data ReferenceFrame (p::Nat) (q::Nat) t = RelativeFrame {frameName :: String, euclideanMove :: EuclideanMove p q t, velocityRelToParentFrame :: Multivector p q t, angularVelocityRelToParentFrame :: Multivector p q t, parent :: ReferenceFrame p q t}+                                        |GlobalAbsoluteFrame deriving (Eq, Show)  +getRigidDisplacementRelToInertial :: (Algebra.Field.C t, Ord t, SingI p, SingI q) =>  ReferenceFrame p q t -> EuclideanMove   p q t+getRigidDisplacementRelToInertial GlobalAbsoluteFrame = mempty+--getRigidDisplacementRelToInertial (RelativeFrame _ displacement mother) = displacement <> (getRigidDisplacementRelToInertial mother) -a `cross` b = (negate $ one)`e`[0,1,2] * (a ∧ b)+getFrameTransformation :: forall (p::Nat) (q::Nat) t . (Algebra.Field.C t, Ord t, SingI p, SingI q) =>  ReferenceFrame p q t -> ReferenceFrame p q t -> EuclideanMove p q t+getFrameTransformation r' r = undefined++{-data InertialFrame (p::Nat) (q::Nat) f t = InertialFrame {objects :: t, changeFrame :: t -> EuclideanMove p q f -> t, frame :: ReferenceFrame p q f}+++instance Functor (InertialFrame p q f) where+    fmap func (InertialFrame objs changeFrame frame) = InertialFrame (func objs) (changeFrame . func) frame++instance (SingI p, SingI q) => Applicative (InertialFrame p q f) where+    pure a = InertialFrame a GlobalAbsoluteFrame where +    (<*>) (InertialFrame func trans1 frame1) (InertialFrame objs trans2 frame2) = if (name frame1)==(name frame2)+                                                                    then InertialFrame (func objs) frame1 +                                                                    else InertialFrame (trans2 (func objs) (getFrameTransformation frame2 frame1)) trans2 frame1+++instance (SingI p, SingI q, Algebra.Field.C f, Ord f) => Monad (InertialFrame p q f) where+    return = pure+    (>>=) (InertialFrame objA changeFrameA frameA)  func = undefined where+        (InertialFrame objB changeFrameB frameB) = func objA+        +-}++a `cross` b = (negate psuedoScalar) * (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+ RigidBody:: (Algebra.Field.C f, Algebra.Module.C f (Multivector p q f)) =>  {bodyName::String, frame::ReferenceFrame p q f, mass :: f, inertia :: Multivector p q f, position :: Multivector p q f, attitude :: Multivector p q f, velocity :: Multivector p q f, angularVelocity :: Multivector p q f                              } -> RigidBody p q f  --makeLenses ''RigidBody doesn't actually work@@ -105,8 +135,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 -} -type Vector3 f =  Multivector 3 0 f-type STVector f = Multivector 3 1 f++ \end{code} \bibliographystyle{IEEEtran} \bibliography{biblio.bib}
src/Numeric/Clifford/Internal.hs view
@@ -1,6 +1,6 @@ {-# OPTIONS_GHC -fllvm -fexcess-precision -optlo-O3 -O3 -optlc-O=3 -Wall #-}-{-# LANGUAGE TypeOperators, TypeFamilies,CPP #-}-module Numeric.Clifford.Internal (myTrace, trie, untrie, enumerate, dimension) where+{-# LANGUAGE TypeOperators, TypeFamilies,CPP, ConstraintKinds, RankNTypes, DataKinds #-}+module Numeric.Clifford.Internal (myTrace, trie, untrie, enumerate, dimension, DefaultField, AllowableCliffordType) where import Numeric.Natural import Prelude hiding (head,tail, null) import Data.MemoTrie@@ -9,12 +9,19 @@ import Data.Bits import Test.QuickCheck import Data.Word+import GHC.TypeLits+import Algebra.Field import qualified Debug.Trace as DebugTrace #ifdef DEBUG myTrace = DebugTrace.trace #else myTrace _ x = x #endif+++type AllowableCliffordType p q f = forall (p::Nat) (q::Nat) f. (Ord f, Algebra.Field.C f, SingI p, SingI q)+type DefaultField = Double+ instance HasTrie Natural where     newtype Natural :->: a = NaturalTrie ((Bool,[Bool]) :->: a)     trie f = NaturalTrie (trie (f . unbitsZ))
src/Numeric/Clifford/LinearOperators.lhs view
@@ -1,8 +1,9 @@ \begin{code}-{-# LANGUAGE NoImplicitPrelude, RankNTypes, KindSignatures, DataKinds, GADTs, FlexibleInstances, UndecidableInstances, InstanceSigs, MultiParamTypeClasses #-}+{-# LANGUAGE NoImplicitPrelude, RankNTypes, KindSignatures, DataKinds, GADTs, FlexibleInstances, UndecidableInstances, InstanceSigs, MultiParamTypeClasses, PolyKinds, ConstraintKinds, UnicodeSyntax, StandaloneDeriving #-}+{-# OPTIONS_HADDOCK show-extensions #-} module Numeric.Clifford.LinearOperators where import qualified NumericPrelude as NP ((.), id)-import NumericPrelude hiding ((.), id, (!!), zipWith, map, length)+import NumericPrelude hiding ((.), id, (!!), zipWith, map, length, zipWith3, and) import Numeric.Clifford.Multivector import Algebra.Algebraic import Algebra.Field@@ -20,6 +21,7 @@ import Data.Semigroupoid import Numeric.Natural import Data.Word+import Algebra.Additive import Numeric.Clifford.Internal import qualified Numeric.Clifford.Blade \end{code}@@ -31,38 +33,87 @@ data LinearOperator' p q f g where     LinearOperator' :: {_operator' :: Multivector p q f -> Multivector p q g} -> LinearOperator' p q f g     LinearOperator :: {_operator :: Multivector p q f -> Multivector p q f} -> LinearOperator' p q f f++getFuncFromOperator :: LinearOperator' p q f g → (Multivector p q f → Multivector p q g)+getFuncFromOperator (LinearOperator' op) = op+getFuncFromOperator (LinearOperator op) = op+ type LinearOperator p q f = LinearOperator' p q f f type LinearOperatorCreator p q f = (Algebra.Algebraic.C f, Ord f, SingI p, SingI q) => Multivector p q f -> LinearOperator p q f +instance (Show g, f ~ g) => Show (LinearOperator' p q f g) where+    show = show . getMatrixElementsFromOperator++instance (Algebra.Field.C f, Algebra.Field.C g, Ord f,  Ord g, SingI p, SingI q) => Eq (LinearOperator' p q f g) where+    a == b = and (map (\ e → (f1 e) == (f2 e)) basisVectors) where+           f1 = getFuncFromOperator a+           f2 = getFuncFromOperator b+ instance Category (LinearOperator' p q) where     id = LinearOperator' NP.id     (.) (LinearOperator' a) (LinearOperator' b)  = LinearOperator' (a NP.. b) +  instance (Algebra.Field.C f, Ord f,Algebra.Field.C g, Ord g, SingI p, SingI q, f~g) => Monoid (LinearOperator' p q f g) where     mempty = id     mappend = (.) +data EuclideanMove p q f where+    EuclideanMove :: ∀ (p::Nat) (q::Nat) f.  (Algebra.Field.C f, Ord f, SingI p, SingI q) => { _rotation :: Multivector p q f, _translation :: Multivector p q f} -> EuclideanMove p q f +deriving instance Eq(EuclideanMove p q f)+deriving instance (Show f) => Show (EuclideanMove p q f) +applyEuclideanMove (EuclideanMove r a) x = (rotate r x) + a -class LinearOperatorClass' (p::Nat) (q::Nat) f g where +++instance (Algebra.Field.C f, Ord f, SingI p, SingI q) => Monoid (EuclideanMove p q f) where+    mempty = EuclideanMove one zero+    mappend (EuclideanMove s b) (EuclideanMove r a) = EuclideanMove rot trans where+          rot = r*s+          trans = (rotate s a) + b++{-instance ∀ a b (p::Nat) (q::Nat).(Algebra.Field.C a, SingI p, SingI q, Ord a, Algebra.Field.C b, Ord b) => Category (AffineOperator' p q) where+    id:: (Algebra.Field.C c) => AffineOperator' p q c c+    id = AffineOperator id zero+    (.) = undefined -}+ {-+--GHC 7.8: The Control.Category module now has the PolyKinds extension enabled, meaning that instances of Category no longer need be of kind * -> * -> *.+class Operator (p::Nat) (q::Nat) f g where+    apply :: Operator p q f g ->  Multivector p q f -> Multivector p q g++instance forall (p::Nat) (q::Nat) . Category (Operator p q) where+    id = NP.id+    (.) a b = a (NP..) b+-}  ++{- [[f11, f12, f13],  [f21, f22, f21],  [f31, f32, f33]] -}-createFunctionalFromElements :: forall (p::Nat) (q::Nat) f . (Algebra.Field.C f, Ord f, SingI p, SingI q) => [[f]] ->(Multivector p q f -> Multivector p q f)++getMatrixElementsFromOperator :: LinearOperator' p q f' f'-> [[f']]+getMatrixElementsFromOperator operator = error "Numeric.Clifford.LinearOperator.getMatrixElementsFromOperator not implemented yet!" where+    f = getFuncFromOperator operator+    +   ++createFunctionalFromElements :: ∀  (p::Nat) (q::Nat) f . (Algebra.Field.C f, Ord f, SingI p, SingI q) => [[f]] ->(Multivector p q f -> Multivector p q f) createFunctionalFromElements elements = (\x -> f*x) where     d = (length elements) - 1     f = sumList $ map elementsForK [0..d]     column k = let transposed = transpose elements in transposed !! k        elementsForK k =sumList $   zipWith (scaleRight) basisVectors (column k)      -createLinearOperatorFromElements :: forall (p::Nat) (q::Nat) f . (Algebra.Field.C f, Ord f, SingI p, SingI q) => [[f]] -> LinearOperator p q f+createLinearOperatorFromElements :: ∀ (p::Nat) (q::Nat) f . (Algebra.Field.C f, Ord f, SingI p, SingI q) => [[f]] -> LinearOperator p q f createLinearOperatorFromElements  = LinearOperator .  createFunctionalFromElements  +reflect :: (Algebra.Algebraic.C f, Ord f, SingI p, SingI q) => Multivector p q f -> Multivector p q f -> Multivector p q f reflect u x = (-u)*x*recip u  makeReflectionOperator ::LinearOperatorCreator p q f@@ -73,10 +124,10 @@  makeRotationOperator :: LinearOperatorCreator p q f makeRotationOperator u = LinearOperator (rotate u)- makeRotationOperatorFromPlaneAngle plane angle = LinearOperator (rotatePlaneAngle plane angle)-project u x = inverse u * (u `dot` x) ++project u x = inverse u * (u `dot` x) makeProjectionOperator :: LinearOperatorCreator p q f makeProjectionOperator u = LinearOperator (project u)
src/Numeric/Clifford/Multivector.lhs view
@@ -21,16 +21,11 @@ {-# LANGUAGE NoImplicitPrelude, FlexibleContexts, RankNTypes, ScopedTypeVariables, DeriveDataTypeable #-} {-# LANGUAGE NoMonomorphismRestriction, UnicodeSyntax, GADTs#-} {-# LANGUAGE FlexibleInstances, StandaloneDeriving, KindSignatures, DataKinds #-}-{-# LANGUAGE TemplateHaskell, TypeOperators, DeriveFunctor #-}+{-# LANGUAGE TemplateHaskell, TypeOperators, DeriveFunctor, DeriveFoldable, DeriveTraversable#-} {-# LANGUAGE MultiParamTypeClasses, UndecidableInstances #-}-\end{code}-%if False-\begin{code}-{-# OPTIONS_GHC -fllvm -fexcess-precision -optlo-O3 -O3 -optlc-O=3 -Wall #-}--- OPTIONS_GHC -Odph -fvectorise -package dph-lifted-vseg ---  LANGUAGE ParallelArrays+{-# OPTIONS_HADDOCK show-extensions #-} \end{code}-%endif+ Clifford algebras are a module over a ring. They also support all the usual transcendental functions. \begin{code} module Numeric.Clifford.Multivector where@@ -77,7 +72,8 @@ import Data.Word import Control.Applicative ((<$>)) import Numeric.Clifford.Internal-+import Data.Traversable+import Data.Foldable (Foldable)  \end{code} @@ -86,8 +82,9 @@  \begin{code} data Multivector (p::Nat) (q::Nat) f where-    BladeSum :: forall (p::Nat) (q::Nat) f . (Ord f, Algebra.Field.C f,SingI p, SingI q) => { _terms :: [Blade p q f]} -> Multivector p q f+    BladeSum :: ∀ (p::Nat) (q::Nat) f . (Ord f, Algebra.Field.C f,SingI p, SingI q) => { _terms :: [Blade p q f]} → Multivector p q f + type STVector = Multivector 3 1 Double type E3Vector = Multivector 3 0 Double @@ -98,12 +95,10 @@        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) + 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)) @@ -118,10 +113,14 @@ terms = lens _terms (\bladeSum v -> bladeSum {_terms = v})  {-# INLINE mvNormalForm #-}-mvNormalForm (BladeSum terms) = BladeSum $ if null resultant then [scalarBlade Algebra.Additive.zero] else resultant  where+mvNormalForm (BladeSum terms) = BladeSum $ mvNormalForm' terms+++mvNormalForm' terms =  if null resultant then [scalarBlade Algebra.Additive.zero] else resultant  where     resultant = filter bladeNonZero $ addLikeTerms' $ Data.List.Ordered.sortBy compare $  map bladeNormalForm $ terms+ {-#INLINE mvTerms #-}-mvTerms m = m^.terms+mvTerms m = _terms m  {-# INLINE addLikeTerms' #-} addLikeTerms' = sumLikeTerms . groupLikeTerms@@ -133,7 +132,8 @@ 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) -{-#INLINE compensatedSum' #-}+{-#NOINLINE compensatedSum' #-}+{-#SPECIALISE compensatedSum' :: [Double] -> Double #-} compensatedSum' :: (Algebra.Additive.C f) => [f] -> f compensatedSum' xs = kahan zero zero xs where     kahan s _ [] = s@@ -153,20 +153,6 @@         t = s + y         newc = (t - s) - y             ---multiplyAdd a b c = a*b + c---twoProduct a b = (x,y) where---    x = a*b---z    y = multiplyAdd a b (negate x)----multiplyList [] = []---multiplyList a@(x:[])=a---multiplyList (a:b:xs) = loop a (b:xs) zero where---  loop pm [] ei = pm+ei---  loop pm1 (ai:remaining) eim1= loop pi remaining ei where---      (pi, pii) = twoProduct pm1 ai---      ei = multiplyAdd eim1 ai pii-- multiplyOutBlades :: (SingI p, SingI q, Algebra.Ring.C a) => [Blade p q a] -> [Blade p q a] -> [Blade p q a] multiplyOutBlades x y = [bladeMul l r | l <-x, r <- y] @@ -179,11 +165,9 @@ multiplyList1 l = mvNormalForm $ BladeSum listOfBlades where     expandedBlades a = foldl1 multiplyOutBlades a     listOfBlades = expandedBlades $ map mvTerms l---things to test: is 1. adding blades into a map based on indices 2. adding errything together 3. sort results quicker than---                   1. sorting by indices 2. groupBy-ing on indices 3. adding the lists of identical indices  {-#INLINE sumList #-}-sumList xs = mvNormalForm $ BladeSum $ concat $ map mvTerms xs+sumList xs = BladeSum $ mvNormalForm' $ concatMap _terms xs   {-#INLINE sumLikeTerms #-} {-#SPECIALISE INLINE sumLikeTerms :: [[STBlade]] -> [STBlade] #-}@@ -260,7 +244,6 @@      a ^ 2 = a * a     a ^ 0 = one-    a ^ 1 = a     --a ^ n  --n < 0 = Clifford.recip $ a ^ (negate n)     a ^ n  =  multiplyList (replicate (NPN.fromInteger n) a) @@ -293,7 +276,6 @@     signum (BladeSum []) = scalar Algebra.Additive.zero  instance (Algebra.Field.C f, Ord f, SingI p, SingI q) => Algebra.Module.C f (Multivector p q f) where---    (*>) zero v = Algebra.Additive.zero     {-#INLINE (*>) #-}     {-#SPECIALISE INLINE (*>) :: Double -> STVector -> STVector #-}     {-#SPECIALISE INLINE (*>) :: Double -> E3Vector -> E3Vector #-}@@ -315,8 +297,8 @@ 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 + {-# INLINE converge#-} converge [] = error "converge: empty list" converge xs = fromMaybe empty (convergeBy checkPeriodic Just xs) @@ -539,23 +521,6 @@ normalised a = (a `scaleRight` ( recip $ mag),mag) where     mag = magnitude a --\end{code}--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}----{- $(derive makeSerialize ''Blade)-$(derive makeSerialize ''Multivector)-$(derive makeData ''Blade)-$(derive makeTypeable ''Blade)-$(derive makeData ''Multivector)-$(derive makeTypeable ''Multivector)-}---- $(derive makeArbitrary ''Multivector)  \end{code} \bibliographystyle{IEEEtran}
src/Numeric/Clifford/NumericIntegration.lhs view
@@ -18,17 +18,7 @@ {-# LANGUAGE NoMonomorphismRestriction, UnicodeSyntax, GADTs, DataKinds, KindSignatures #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE TemplateHaskell #-}-{-# LANGUAGE MultiParamTypeClasses #-}-\end{code}-%if False-\begin{code}-{-# OPTIONS_GHC -fllvm -fexcess-precision -optlo-O3 -O3 -optlc-O=3 -Wall #-}--- OPTIONS_GHC -Odph -fvectorise -package dph-lifted-vseg ---  LANGUAGE ParallelArrays-\end{code}-%endif--\begin{code}+{-# LANGUAGE MultiParamTypeClasses, BangPatterns #-}  module Numeric.Clifford.NumericIntegration where import Numeric.Clifford.Multivector as MV@@ -57,6 +47,7 @@ import Numeric.Clifford.Internal import Data.DeriveTH +import qualified Debug.Trace  elementAdd = zipWith (+) elementScale = zipWith (*>) @@ -123,7 +114,7 @@ type AdaptiveStepSizeFunction f state = f -> state -> f   data RKAttribute f state = Explicit-                 | HamiltonianFunction+                 | HamiltonianFunction {totalEnergy :: state -> f}                  | AdaptiveStepSize {sigma :: AdaptiveStepSizeFunction f state}                  | ConvergenceTolerance {epsilon :: f}                  | ConvergenceFunction {converger :: ConvergerFunction f } @@ -172,7 +163,7 @@                         Just (AdaptiveStepSize sigma) -> sigma                         Nothing -> (\_ _ -> one) -    {-#INLINE rkMethodImplicitFixedPoint#-}+ --    {-#SPECIALISE rkMethodImplicitFixedPoint :: RKStepper 3 0 Double stateType #-} --    {-#SPECIALISE rkMethodImplicitFixedPoint :: RKStepper 3 0 Double [E3Vector] #-} --    {-#SPECIALISE rkMethodImplicitFixedPoint :: RKStepper 3 1 Double stateType #-}@@ -192,7 +183,8 @@             guessTime :: t             guessTime = time + h'             zkp1 :: NPN.Int -> [Multivector p q t] -> [Multivector p q t]-            zkp1 i zk = map (h*>) (sumOfJs i zk) where+            zkp1 i !zk = map (h*>) (sumOfJs i zk) where+                                 {-#INLINE sumOfJs#-}                 sumOfJs :: Int -> [Multivector p q t] -> [Multivector p q t]                 sumOfJs i zk =  sumListOfLists $ map (scaledByAij zk) (a i) where
src/Numeric/Clifford/NumericIntegration/DefaultIntegrators.hs view
@@ -36,9 +36,9 @@ import           Test.QuickCheck --trace _ a = a -gaussLegendreFourthOrderTableau = ButcherTableau [[0.25::NPN.Double, 0.25 - (1.0 NPN./6)* sqrt 3], [0.25 + (1.0 NPN./ 6) * sqrt 3, 0.25]] [0.5, 0.5] [0.5 - (1.0 NPN./6)* sqrt 3, 0.5 + (1.0 NPN./ 6)* sqrt 3]+gaussLegendreFourthOrderTableau = ButcherTableau [[0.25::NPN.Double, 0.25 - ((sqrt 3.0) NPN./6.0)],[0.25 + ((sqrt 3.0) NPN./ 6.0) , 0.25]] [0.5, 0.5] [0.5 - ((sqrt 3.0) NPN./6.0), 0.5 + ((sqrt 3.0) NPN./ 6.0)] gaussLegendreFourthOrder h f (t, state) = impl h f id id (t,state) where-    impl= genericRKMethod gaussLegendreFourthOrderTableau []+    impl= genericRKMethod gaussLegendreFourthOrderTableau [ConvergenceTolerance 1.0e-8]   rk4ClassicalFromTableau h f (t,state) = impl h f id id (t, state) where