From 25d738b25253b168c1cea5e2584130a1fcbe972a Mon Sep 17 00:00:00 2001 From: sheaf Date: Fri, 20 Jan 2023 16:34:04 +0100 Subject: [PATCH] more WIP with TH now --- MetaBrush.cabal | 19 +- cabal.project | 8 +- src/app/MetaBrush/Application.hs | 6 +- src/app/MetaBrush/Render/Document.hs | 10 +- src/metabrushes/MetaBrush/Asset/Brushes.hs | 66 +- src/metabrushes/MetaBrush/Brush.hs | 47 +- src/metabrushes/MetaBrush/Records.hs | 58 +- src/metabrushes/MetaBrush/Serialisable.hs | 8 +- src/splines/Math/Algebra/Dual.hs | 761 ++++++++++++++++++ src/splines/Math/Algebra/Dual/Internal.hs | 588 ++++++++++++++ src/splines/Math/Bezier/Cubic.hs | 5 +- src/splines/Math/Bezier/Quadratic.hs | 5 +- src/splines/Math/Bezier/Stroke.hs | 157 ++-- src/splines/Math/Differentiable.hs | 66 ++ src/splines/Math/Interval.hs | 441 +++++++++++ src/splines/Math/Interval/Internal.hs | 82 ++ src/splines/Math/Linear.hs | 254 +----- src/splines/Math/Linear/Dual.hs | 863 --------------------- src/splines/Math/Linear/Internal.hs | 154 ++++ src/splines/Math/Module.hs | 149 ++-- src/splines/Math/Module/Internal.hs | 90 +++ src/splines/Math/Monomial.hs | 174 +++++ src/splines/Math/Ring.hs | 177 +++++ src/splines/TH/Utils.hs | 23 + 24 files changed, 2885 insertions(+), 1326 deletions(-) create mode 100644 src/splines/Math/Algebra/Dual.hs create mode 100644 src/splines/Math/Algebra/Dual/Internal.hs create mode 100644 src/splines/Math/Differentiable.hs create mode 100644 src/splines/Math/Interval.hs create mode 100644 src/splines/Math/Interval/Internal.hs delete mode 100644 src/splines/Math/Linear/Dual.hs create mode 100644 src/splines/Math/Linear/Internal.hs create mode 100644 src/splines/Math/Module/Internal.hs create mode 100644 src/splines/Math/Monomial.hs create mode 100644 src/splines/Math/Ring.hs create mode 100644 src/splines/TH/Utils.hs diff --git a/MetaBrush.cabal b/MetaBrush.cabal index 84ba0ee..62b376f 100644 --- a/MetaBrush.cabal +++ b/MetaBrush.cabal @@ -75,6 +75,7 @@ common common MagicHash MultiWayIf NamedFieldPuns + NoStarIsType PatternSynonyms RankNTypes RecordWildCards @@ -92,7 +93,6 @@ common common -O1 -fexpose-all-unfoldings -funfolding-use-threshold=16 - -fexcess-precision -fspecialise-aggressively -optc-O3 -Wall @@ -173,20 +173,31 @@ library splines Haskell2010 exposed-modules: - Math.Bezier.Cubic + Math.Algebra.Dual + , Math.Bezier.Cubic , Math.Bezier.Cubic.Fit , Math.Bezier.Quadratic , Math.Bezier.Spline , Math.Bezier.Stroke + , Math.Differentiable , Math.Epsilon + , Math.Interval , Math.Linear - , Math.Linear.Dual , Math.Linear.Solve , Math.Module + , Math.Monomial , Math.Orientation + , Math.Ring , Math.Roots , Debug.Utils + other-modules: + Math.Algebra.Dual.Internal + , Math.Interval.Internal + , Math.Linear.Internal + , Math.Module.Internal + , TH.Utils + build-depends: bifunctors >= 5.5.4 && < 5.6 @@ -198,6 +209,8 @@ library splines ^>= 0.2 , vector >= 0.12.1.2 && < 0.14 + , template-haskell + >= 2.18 && < 2.20 library metabrushes diff --git a/cabal.project b/cabal.project index ef5dd78..b162dbc 100644 --- a/cabal.project +++ b/cabal.project @@ -28,10 +28,10 @@ allow-newer: waargonaut:vector, waargonaut:witherable, -- eigen -source-repository-package - type: git - location: https://github.com/chessai/eigen - tag: 8fff32a43df743c8c83428a86dd566a0936a4fba +--source-repository-package +-- type: git +-- location: https://github.com/chessai/eigen +-- tag: 8fff32a43df743c8c83428a86dd566a0936a4fba -- records-sop source-repository-package diff --git a/src/app/MetaBrush/Application.hs b/src/app/MetaBrush/Application.hs index a59b501..b3349a4 100644 --- a/src/app/MetaBrush/Application.hs +++ b/src/app/MetaBrush/Application.hs @@ -200,11 +200,11 @@ runApplication application = do maxHistorySizeTVar <- STM.newTVarIO @Int 1000 fitParametersTVar <- STM.newTVarIO @FitParameters ( FitParameters - { maxSubdiv = 3 --2 --3 -- 6 - , nbSegments = 40 --3 --6 -- 12 + { maxSubdiv = 1 --2 --3 -- 6 + , nbSegments = 2 --3 --6 -- 12 , dist_tol = 0.1 -- 5e-3 , t_tol = 0.1 -- 1e-4 - , maxIters = 2 -- 100 + , maxIters = 1 -- 100 } ) diff --git a/src/app/MetaBrush/Render/Document.hs b/src/app/MetaBrush/Render/Document.hs index b826c14..561f6e9 100644 --- a/src/app/MetaBrush/Render/Document.hs +++ b/src/app/MetaBrush/Render/Document.hs @@ -56,6 +56,8 @@ import Control.Monad.Trans.State.Strict ( StateT, evalStateT, get, put ) -- MetaBrush +import Math.Algebra.Dual + ( fun ) import qualified Math.Bezier.Cubic as Cubic ( Bezier(..), fromQuadratic ) import Math.Bezier.Cubic.Fit @@ -73,10 +75,10 @@ import Math.Bezier.Stroke ( CachedStroke(..), invalidateCache , computeStrokeOutline ) +import Math.Interval + ( Extent(Point) ) import Math.Linear - ( ℝ(..), T(..), Extent(Point) ) -import Math.Linear.Dual - ( fun ) + ( ℝ(..), T(..) ) import MetaBrush.Asset.Colours ( Colours, ColourRecord(..) ) import MetaBrush.Brush @@ -305,7 +307,7 @@ strokeRenderData fitParams StrokeWithOutlineRenderData { strokeDataSpline = spline , strokeOutlineData = ( outline, fitPts ) - , strokeBrushFunction = fun ( brushFn @Point proxy# id ) + , strokeBrushFunction = fun @Double ( brushFn @Point proxy# id ) . embedUsedParams . toUsedParams } diff --git a/src/metabrushes/MetaBrush/Asset/Brushes.hs b/src/metabrushes/MetaBrush/Asset/Brushes.hs index 886755a..b2e6637 100644 --- a/src/metabrushes/MetaBrush/Asset/Brushes.hs +++ b/src/metabrushes/MetaBrush/Asset/Brushes.hs @@ -1,10 +1,14 @@ {-# LANGUAGE AllowAmbiguousTypes #-} {-# LANGUAGE OverloadedStrings #-} {-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TemplateHaskell #-} module MetaBrush.Asset.Brushes where -- base +import Prelude + hiding + ( Num(..), Floating(..) ) import GHC.Exts ( Proxy# ) @@ -23,12 +27,16 @@ import qualified Data.HashMap.Strict as HashMap ( fromList, lookup ) -- MetaBrush +import Math.Algebra.Dual import Math.Bezier.Spline +import Math.Differentiable + ( DiffInterp ) +import Math.Interval + ( type I ) import Math.Linear -import Math.Linear.Dual - ( D, type (~>)(..), Differentiable, Diffy(konst), var ) import Math.Module ( Module((^+^), (*^)) ) +import Math.Ring import MetaBrush.Brush ( Brush(..), SomeBrush(..), WithParams(..) ) import MetaBrush.Records @@ -70,62 +78,62 @@ ellipse = BrushData "ellipse" ( WithParams deflts ellipseBrush ) -------------------------------------------------------------------------------- -- Differentiable brushes. -circleSpline :: forall i u v - . Applicative ( D ( I i u ) ) - => ( Double -> Double -> D ( I i u ) ( I i v ) ) - -> D ( I i u ) ( Spline 'Closed () ( I i v ) ) +circleSpline :: forall i k u v + . Applicative ( D k ( I i u ) ) + => ( Double -> Double -> D k ( I i u ) ( I i v ) ) + -> D k ( I i u ) ( Spline 'Closed () ( I i v ) ) circleSpline p = sequenceA $ Spline { splineStart = p 1 0 , splineCurves = ClosedCurves crvs lastCrv } where crvs = Seq.fromList - [ Bezier3To ( p 1 κ ) ( p κ 1 ) ( NextPoint (p 0 1) ) () - , Bezier3To ( p -κ 1 ) ( p -1 κ ) ( NextPoint (p -1 0) ) () - , Bezier3To ( p -1 -κ ) ( p -κ -1 ) ( NextPoint (p 0 -1) ) () + [ Bezier3To ( p 1 κ ) ( p κ 1 ) ( NextPoint ( p 0 1 ) ) () + , Bezier3To ( p -κ 1 ) ( p -1 κ ) ( NextPoint ( p -1 0 ) ) () + , Bezier3To ( p -1 -κ ) ( p -κ -1 ) ( NextPoint ( p 0 -1 ) ) () ] lastCrv = Bezier3To ( p κ -1 ) ( p 1 -κ ) BackToStart () -circleBrush :: forall i - . ( Differentiable i ( Record CircleBrushFields ) ) +circleBrush :: forall i k + . ( k ~ 2, DiffInterp i ( Record CircleBrushFields ) ) => Proxy# i -> ( forall a. a -> I i a ) - -> I i ( Record CircleBrushFields ) ~> Spline 'Closed () ( I i ( ℝ 2 ) ) + -> C k ( I i ( Record CircleBrushFields ) ) ( Spline 'Closed () ( I i ( ℝ 2 ) ) ) circleBrush _ mkI = D \ params -> - let r :: D ( I i ( Record CircleBrushFields ) ) ( I i Double ) - r = runD ( var ( Fin 1## ) ) params - mkPt :: Double -> Double -> D ( I i ( Record CircleBrushFields ) ) ( I i ( ℝ 2 ) ) + let r :: D k ( I i ( Record CircleBrushFields ) ) ( I i Double ) + r = runD ( var @_ @k ( Fin 1 ) ) params + mkPt :: Double -> Double -> D k ( I i ( Record CircleBrushFields ) ) ( I i ( ℝ 2 ) ) mkPt ( kon -> x ) ( kon -> y ) = ( x * r ) *^ e_x ^+^ ( y * r ) *^ e_y - in circleSpline @i @( Record CircleBrushFields ) @( ℝ 2 ) mkPt + in circleSpline @i @k @( Record CircleBrushFields ) @( ℝ 2 ) mkPt where - e_x, e_y :: D ( I i ( Record CircleBrushFields ) ) ( I i ( ℝ 2 ) ) + e_x, e_y :: D k ( I i ( Record CircleBrushFields ) ) ( I i ( ℝ 2 ) ) e_x = pure $ mkI $ ℝ2 1 0 e_y = pure $ mkI $ ℝ2 0 1 - kon = konst @( I i Double ) @( I i ( Record CircleBrushFields ) ) . mkI + kon = konst @( I i Double ) @k @( I i ( Record CircleBrushFields ) ) . mkI -ellipseBrush :: forall i - . ( Differentiable i ( Record EllipseBrushFields ) ) +ellipseBrush :: forall i k + . ( k ~ 2, DiffInterp i ( Record EllipseBrushFields ) ) => Proxy# i -> ( forall a. a -> I i a ) - -> I i ( Record EllipseBrushFields ) ~> Spline 'Closed () ( I i ( ℝ 2 ) ) + -> C k ( I i ( Record EllipseBrushFields ) ) ( Spline 'Closed () ( I i ( ℝ 2 ) ) ) ellipseBrush _ mkI = D \ params -> - let a, b, phi :: D ( I i ( Record EllipseBrushFields ) ) ( I i Double ) - a = runD ( var ( Fin 1## ) ) params - b = runD ( var ( Fin 2## ) ) params - phi = runD ( var ( Fin 3## ) ) params - mkPt :: Double -> Double -> D ( I i ( Record EllipseBrushFields ) ) ( I i ( ℝ 2 ) ) + let a, b, phi :: D k ( I i ( Record EllipseBrushFields ) ) ( I i Double ) + a = runD ( var @_ @k ( Fin 1 ) ) params + b = runD ( var @_ @k ( Fin 2 ) ) params + phi = runD ( var @_ @k ( Fin 3 ) ) params + mkPt :: Double -> Double -> D k ( I i ( Record EllipseBrushFields ) ) ( I i ( ℝ 2 ) ) mkPt ( kon -> x ) ( kon -> y ) = ( x * a * cos phi - y * b * sin phi ) *^ e_x ^+^ ( y * b * cos phi + x * a * sin phi ) *^ e_y - in circleSpline @i @( Record EllipseBrushFields ) @( ℝ 2 ) mkPt + in circleSpline @i @k @( Record EllipseBrushFields ) @( ℝ 2 ) mkPt where - e_x, e_y :: D ( I i ( Record EllipseBrushFields ) ) ( I i ( ℝ 2 ) ) + e_x, e_y :: D k ( I i ( Record EllipseBrushFields ) ) ( I i ( ℝ 2 ) ) e_x = pure $ mkI $ ℝ2 1 0 e_y = pure $ mkI $ ℝ2 0 1 - kon = konst @( I i Double ) @( I i ( Record EllipseBrushFields ) ) . mkI + kon = konst @( I i Double ) @k @( I i ( Record EllipseBrushFields ) ) . mkI diff --git a/src/metabrushes/MetaBrush/Brush.hs b/src/metabrushes/MetaBrush/Brush.hs index 229d7c5..e4b3996 100644 --- a/src/metabrushes/MetaBrush/Brush.hs +++ b/src/metabrushes/MetaBrush/Brush.hs @@ -39,12 +39,16 @@ import qualified Data.Text as Text ( unpack ) -- MetaBrush +import Math.Algebra.Dual + ( type (~>) ) import Math.Linear - ( ℝ, type I, Extent(Point, Interval) ) -import Math.Linear.Dual - ( type (~>), Differentiable ) +import Math.Interval + ( type I, Extent(Point, Interval) ) + import Math.Bezier.Spline ( SplineType(Closed), Spline ) +import Math.Differentiable + ( DiffInterp ) import MetaBrush.Records ( KnownSymbols, Length, Record ) import MetaBrush.Serialisable @@ -54,30 +58,33 @@ import MetaBrush.Serialisable -- | A differentiable function from a given record type, -- with provided default values that can be overridden. -type WithParams :: [ Symbol ] -> ( Type -> Type ) -> Type +type WithParams :: Type -> ( Type -> Type ) -> Type data WithParams params f = WithParams - { defaultParams :: Record params - , withParams :: forall i - . ( Differentiable i ( Record params ) ) - => Proxy# i - -> ( forall a. a -> I i a ) - -> I i ( Record params ) ~> f ( I i ( ℝ 2 ) ) + { defaultParams :: params + , withParams + :: forall i + . ( DiffInterp i params ) + => Proxy# i + -> ( forall a. a -> I i a ) + -> I i params ~> f ( I i ( ℝ 2 ) ) } -------------------------------------------------------------------------------- -- | A brush function: a function from a record of parameters to a closed spline. type BrushFunction :: [ Symbol ] -> Type -type BrushFunction brushFields = WithParams brushFields ( Spline Closed () ) +type BrushFunction brushFields = + WithParams ( Record brushFields ) ( Spline Closed () ) type Brush :: [ Symbol ] -> Type data Brush brushFields where BrushData :: forall brushFields . ( KnownSymbols brushFields, Typeable brushFields - , Differentiable Point ( ℝ ( Length brushFields ) ) - , Differentiable Interval ( ℝ ( Length brushFields ) ) + , Representable Double ( ℝ ( Length brushFields ) ) + , DiffInterp Point ( ℝ ( Length brushFields ) ) + , DiffInterp Interval ( ℝ ( Length brushFields ) ) ) => { brushName :: !Text , brushFunction :: BrushFunction brushFields @@ -111,21 +118,25 @@ class ( KnownSymbols pointFields, Typeable pointFields , Serialisable ( Record pointFields ) , Show ( Record pointFields ) , NFData ( Record pointFields ) - , Differentiable Point ( ℝ ( Length pointFields ) ) - , Differentiable Interval ( ℝ ( Length pointFields ) ) + , Representable Double ( ℝ ( Length pointFields ) ) + , DiffInterp Point ( ℝ ( Length pointFields ) ) + , DiffInterp Interval ( ℝ ( Length pointFields ) ) ) => PointFields pointFields where { } instance ( KnownSymbols pointFields, Typeable pointFields , Serialisable ( Record pointFields ) , Show ( Record pointFields ) , NFData ( Record pointFields ) - , Differentiable Point ( ℝ ( Length pointFields ) ) - , Differentiable Interval ( ℝ ( Length pointFields ) ) + , Representable Double ( ℝ ( Length pointFields ) ) + , DiffInterp Point ( ℝ ( Length pointFields ) ) + , DiffInterp Interval ( ℝ ( Length pointFields ) ) ) => PointFields pointFields where { } -- | Assumes the input has no duplicates (doesn't check.) -provePointFields :: [ Text ] -> ( forall pointFields. PointFields pointFields => Proxy# pointFields -> r ) -> r +provePointFields :: [ Text ] + -> ( forall pointFields. PointFields pointFields => Proxy# pointFields -> r ) + -> r provePointFields fieldNames k = case fieldNames of [] diff --git a/src/metabrushes/MetaBrush/Records.hs b/src/metabrushes/MetaBrush/Records.hs index f5cadcc..4609e1a 100644 --- a/src/metabrushes/MetaBrush/Records.hs +++ b/src/metabrushes/MetaBrush/Records.hs @@ -17,7 +17,7 @@ import Data.Typeable import Data.Type.Equality ( (:~:)(Refl) ) import GHC.Exts - ( Word(W#), Proxy#, proxy# ) + ( Proxy#, proxy# ) import GHC.Show ( showCommaSpace ) import GHC.TypeLits @@ -54,9 +54,10 @@ import qualified Data.Text as Text ( pack, unpack ) -- MetaBrush +import Math.Algebra.Dual +import Math.Differentiable +import Math.Interval import Math.Linear -import Math.Linear.Dual - ( type (~>)(..), D, Diffy, Differentiable ) import Math.Module ( Module ) @@ -87,7 +88,7 @@ instance ( KnownSymbols ks, Representable Double ( ℝ ( Length ks ) ) ) where fields :: [ ShowS ] fields = - zip [ 1.. ] ( knownSymbols @ks ) <&> \ ( W# i, fld ) -> + zip [ 1.. ] ( knownSymbols @ks ) <&> \ ( i, fld ) -> let v = index r ( Fin i ) in showString ( Text.unpack fld ) . showString " = " . showsPrec 0 v @@ -143,17 +144,18 @@ instance ( Torsor ( T ( 𝕀ℝ ( Length ks ) ) ) ( 𝕀ℝ ( Length ks ) ) T ( I ( Rounded c_lo ) ( Rounded c_hi ) ) -> T ( I ( Rounded ( MkR c_lo ) ) ( Rounded ( MkR c_hi ) ) ) +type instance RepDim ( Record ks ) = Length ks deriving newtype instance Representable r ( ℝ ( Length ks ) ) => Representable r ( Record ks ) -type instance D ( Record ks ) = D ( ℝ ( Length ks ) ) -deriving newtype instance Diffy Double ( ℝ ( Length ks ) ) - => Diffy Double ( Record ks ) +type instance D k ( Record ks ) = D k ( ℝ ( Length ks ) ) +deriving newtype instance HasChainRule Double 2 ( ℝ ( Length ks ) ) + => HasChainRule Double 2 ( Record ks ) deriving via 𝕀ℝ ( Length ks ) - instance Diffy ( 𝕀 Double ) ( 𝕀ℝ ( Length ks ) ) - => Diffy ( 𝕀 Double ) ( 𝕀 ( Record ks ) ) + instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ ( Length ks ) ) + => HasChainRule ( 𝕀 Double ) 2 ( 𝕀 ( Record ks ) ) -------------------------------------------------------------------------------- @@ -181,6 +183,7 @@ intersect :: forall r1 r2 l1 l2 , KnownSymbols r1, KnownSymbols r2 , l1 ~ Length r1, l2 ~ Length r2 , Representable Double ( ℝ l1 ) + , Representable Double ( ℝ l2 ) , Differentiable 'Point ( ℝ l2 ) , Differentiable 'Interval ( ℝ l2 ) ) @@ -197,7 +200,7 @@ intersect project = \ ( MkR r1 ) -> MkR $ projection ( (!) r1_idxs ) r1 inject :: Record r2 -> Record r1r2 -> Record r2 - inject = \ ( MkR r2 ) -> \ ( MkR r1r2 ) -> MkR $ injection ( find eqFin r2_idxs ) r1r2 r2 + inject = \ ( MkR r2 ) -> \ ( MkR r1r2 ) -> MkR $ injection ( \ i -> find ( == i ) r2_idxs ) r1r2 r2 in Intersection { project, inject } data Intersection r1 r2 where @@ -205,6 +208,7 @@ data Intersection r1 r2 where :: forall r1r2 r1 r2 l12 . ( l12 ~ Length r1r2 , KnownSymbols r1r2 + , Representable Double ( ℝ l12 ) , Differentiable 'Point ( ℝ l12 ) , Differentiable 'Interval ( ℝ l12 ) ) @@ -221,13 +225,13 @@ doIntersection :: forall r1 r2 l1 l2 kont . ( KnownSymbols r1, KnownSymbols r2 , l1 ~ Length r1, l2 ~ Length r2 - , Representable Double ( ℝ l1 ) - , Representable Double ( ℝ l2 ) ) => ( forall r1r2 l12. - ( r1r2 ~ Intersect r1 r2, KnownSymbols r1r2, l12 ~ Length r1r2 + ( r1r2 ~ Intersect r1 r2 + , KnownSymbols r1r2, l12 ~ Length r1r2 , Differentiable 'Point ( ℝ l12 ) , Differentiable 'Interval ( ℝ l12 ) + , Representable Double ( ℝ l12 ) ) => Proxy# r1r2 -> Vec l12 ( Fin l1 ) -> Vec l12 ( Fin l2 ) -> kont ) -> kont @@ -237,11 +241,11 @@ doIntersection k = [ ] | ( _ :: Proxy# r1r2 ) <- proxy# @'[ ] , Refl <- ( unsafeCoerce Refl :: r1r2 :~: Intersect r1 r2 ) - -> k @'[] proxy# + -> k @r1r2 proxy# VZ VZ - [ ( f1, W# r1_i1, W# r2_i1 ) ] + [ ( f1, r1_i1, r2_i1 ) ] | SomeSymbol @f1 _ <- someSymbolVal ( Text.unpack f1 ) , ( _ :: Proxy# r1r2 ) <- proxy# @'[ f1 ] , Refl <- ( unsafeCoerce Refl :: r1r2 :~: Intersect r1 r2 ) @@ -249,8 +253,8 @@ doIntersection k = ( VS ( Fin r1_i1 ) VZ ) ( VS ( Fin r2_i1 ) VZ ) - [ ( f1, W# r1_i1, W# r2_i1 ) - , ( f2, W# r1_i2, W# r2_i2 ) ] + [ ( f1, r1_i1, r2_i1 ) + , ( f2, r1_i2, r2_i2 ) ] | SomeSymbol @f1 _ <- someSymbolVal ( Text.unpack f1 ) , SomeSymbol @f2 _ <- someSymbolVal ( Text.unpack f2 ) , ( _ :: Proxy# r1r2 ) <- proxy# @'[ f1, f2 ] @@ -259,9 +263,9 @@ doIntersection k = ( VS ( Fin r1_i1 ) $ VS ( Fin r1_i2 ) VZ ) ( VS ( Fin r2_i1 ) $ VS ( Fin r2_i2 ) VZ ) - [ ( f1, W# r1_i1, W# r2_i1 ) - , ( f2, W# r1_i2, W# r2_i2 ) - , ( f3, W# r1_i3, W# r2_i3 ) ] + [ ( f1, r1_i1, r2_i1 ) + , ( f2, r1_i2, r2_i2 ) + , ( f3, r1_i3, r2_i3 ) ] | SomeSymbol @f1 _ <- someSymbolVal ( Text.unpack f1 ) , SomeSymbol @f2 _ <- someSymbolVal ( Text.unpack f2 ) , SomeSymbol @f3 _ <- someSymbolVal ( Text.unpack f3 ) @@ -271,6 +275,20 @@ doIntersection k = ( VS ( Fin r1_i1 ) $ VS ( Fin r1_i2 ) $ VS ( Fin r1_i3 ) VZ ) ( VS ( Fin r2_i1 ) $ VS ( Fin r2_i2 ) $ VS ( Fin r2_i3 ) VZ ) + [ ( f1, r1_i1, r2_i1 ) + , ( f2, r1_i2, r2_i2 ) + , ( f3, r1_i3, r2_i3 ) + , ( f4, r1_i4, r2_i4 ) ] + | SomeSymbol @f1 _ <- someSymbolVal ( Text.unpack f1 ) + , SomeSymbol @f2 _ <- someSymbolVal ( Text.unpack f2 ) + , SomeSymbol @f3 _ <- someSymbolVal ( Text.unpack f3 ) + , SomeSymbol @f4 _ <- someSymbolVal ( Text.unpack f4 ) + , ( _ :: Proxy# r1r2 ) <- proxy# @'[ f1, f2, f3, f4 ] + , Refl <- ( unsafeCoerce Refl :: r1r2 :~: Intersect r1 r2 ) + -> k @r1r2 proxy# + ( VS ( Fin r1_i1 ) $ VS ( Fin r1_i2 ) $ VS ( Fin r1_i3 ) $ VS ( Fin r1_i4 ) VZ ) + ( VS ( Fin r2_i1 ) $ VS ( Fin r2_i2 ) $ VS ( Fin r2_i3 ) $ VS ( Fin r2_i4 ) VZ ) + other -> error $ "Intersection not defined in dimension " ++ show ( length other ) ------ diff --git a/src/metabrushes/MetaBrush/Serialisable.hs b/src/metabrushes/MetaBrush/Serialisable.hs index dd953a2..d7c75a8 100644 --- a/src/metabrushes/MetaBrush/Serialisable.hs +++ b/src/metabrushes/MetaBrush/Serialisable.hs @@ -29,8 +29,6 @@ import Data.STRef ( newSTRef ) import Data.Traversable ( for ) -import GHC.Exts - ( Word(W#) ) -- containers import Data.Map.Strict @@ -116,14 +114,14 @@ instance ( KnownSymbols ks, Representable Double ( ℝ ( Length ks ) ) ) where encodeFields :: Record ks -> [ ( Text, Double ) ] encodeFields ( MkR r ) = - zip [1..] ( knownSymbols @ks ) <&> \ ( W# i, fld ) -> + zip [1..] ( knownSymbols @ks ) <&> \ ( i, fld ) -> ( fld, index r ( Fin i ) ) decoder = fmap decodeFields $ for ( knownSymbols @ks ) \ k -> JSON.Decoder.atKey k ( decoder @Double ) where decodeFields :: [ Double ] -> Record ks - decodeFields coords = MkR $ tabulate \ ( Fin i# ) -> - coords !! ( fromIntegral ( W# i# ) - 1 ) + decodeFields coords = MkR $ tabulate \ ( Fin i ) -> + coords !! ( fromIntegral i - 1 ) -------------------------------------------------------------------------------- diff --git a/src/splines/Math/Algebra/Dual.hs b/src/splines/Math/Algebra/Dual.hs new file mode 100644 index 0000000..cc5e841 --- /dev/null +++ b/src/splines/Math/Algebra/Dual.hs @@ -0,0 +1,761 @@ +{-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE RebindableSyntax #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE UndecidableInstances #-} + +{-# OPTIONS_GHC -Wno-orphans -O2 #-} + +{-# OPTIONS_GHC -ddump-splices -ddump-simpl -ddump-to-file -dno-typeable-binds + -dsuppress-unfoldings -dsuppress-coercions #-} + +module Math.Algebra.Dual + ( C(..), D, type (~>), type (~~>) + , HasChainRule(..), chainRule + , uncurryD2, uncurryD3 + , linear, fun, var + + , D𝔸0(..) + , D1𝔸1(..), D2𝔸1(..), D3𝔸1(..) + , D1𝔸2(..), D2𝔸2(..), D3𝔸2(..) + , D1𝔸3(..), D2𝔸3(..), D3𝔸3(..) + , D1𝔸4(..), D2𝔸4(..), D3𝔸4(..) + ) where + +-- base +import Prelude hiding ( Num(..), Floating(..), (^) ) +import Control.Applicative + ( liftA2 ) +import Data.Coerce + ( coerce ) +import Data.Kind + ( Type ) +import Data.Monoid + ( Ap(..) ) +import GHC.TypeNats + ( Nat ) + +-- MetaBrush +import Math.Algebra.Dual.Internal +import Math.Linear +import Math.Module +import Math.Monomial +import Math.Ring + +-------------------------------------------------------------------------------- + +-- | @C n u v@ is the space of @C^k@-differentiable maps from @u@ to @v@. +type C :: Nat -> Type -> Type -> Type +newtype C k u v = D { runD :: u -> D k u v } +deriving stock instance Functor ( D k u ) => Functor ( C k u ) + +-- | \( C^2 \)-differentiable mappings. +type (~>) = C 2 +-- | \( C^3 \)-differentiable mappings. +type (~~>) = C 3 + +-- | @D k u v@ is the space of @k@-th order germs of functions from @u@ to @v@, +-- represented by the algebra: +-- +-- \[ \mathbb{Z}[x_1, \ldots, x_n]/(x_1, \ldots, x_n)^(k+1) \otimes_\mathbb{Z} v \] +-- +-- when @u@ is of dimension @n@. +type D :: Nat -> Type -> Type -> Type +type family D k u + +type instance D k ( ℝ 0 ) = D𝔸0 + +type instance D 1 ( ℝ 1 ) = D1𝔸1 +type instance D 1 ( ℝ 2 ) = D1𝔸2 +type instance D 1 ( ℝ 3 ) = D1𝔸3 +type instance D 1 ( ℝ 4 ) = D1𝔸4 + +type instance D 2 ( ℝ 1 ) = D2𝔸1 +type instance D 2 ( ℝ 2 ) = D2𝔸2 +type instance D 2 ( ℝ 3 ) = D2𝔸3 +type instance D 2 ( ℝ 4 ) = D2𝔸4 + +type instance D 3 ( ℝ 1 ) = D3𝔸1 +type instance D 3 ( ℝ 2 ) = D3𝔸2 +type instance D 3 ( ℝ 3 ) = D3𝔸3 +type instance D 3 ( ℝ 4 ) = D3𝔸4 + +-------------------------------------------------------------------------------- + +-- Weird instance needed in just one place; +-- see use of chain in 'Math.Bezier.Stroke.brushStrokeData'. +instance ( Applicative ( D k u ), Module r ( T v ) ) + => Module r ( T ( C k u v ) ) where + origin = T $ D \ _ -> pure $ coerce $ origin @r @( T v ) + T ( D f ) ^+^ T ( D g ) = T $ D \ t -> liftA2 ( coerce $ (^+^) @r @( T v ) ) ( f t ) ( g t ) + T ( D f ) ^-^ T ( D g ) = T $ D \ t -> liftA2 ( coerce $ (^-^) @r @( T v ) ) ( f t ) ( g t ) + a *^ T ( D f ) = T $ D \ t -> fmap ( coerce $ (*^) @r @( T v ) a ) $ f t + +-------------------------------------------------------------------------------- + +-- | @HasChainRule r k v@ means we have a chain rule +-- with @D k v w@ in the middle, for any @r@-module @w@. +class HasChainRule r k v where + chain :: Module r ( T w ) + => D k ( ℝ 1 ) v -> D k v w -> D k ( ℝ 1 ) w + konst :: AbelianGroup w => w -> D k v w + value :: D k v w -> w + linearD :: Module r ( T w ) => ( v -> w ) -> v -> D k v w + +linear :: forall k r v w + . ( HasChainRule r k v, Module r ( T w ) ) + => ( v -> w ) -> C k v w +linear f = D \ x -> linearD @r @k @v @w f x + +chainRule :: forall r k u v w + . ( HasChainRule r k v, Module r ( T w ) + , D k u ~ D k ( ℝ 1 ), HasChainRule r k u + ) + => C k u v -> C k v w -> C k u w +chainRule ( D df ) ( D dg ) = + D \ x -> + case df x of + df_x -> + chain @r @k @v df_x ( dg $ value @r @k @u df_x ) + + +uncurryD2 :: D 2 a ~ D 2 ( ℝ 1 ) + => D 2 ( ℝ 1 ) ( C 2 a b ) -> a -> D 2 ( ℝ 2 ) b +uncurryD2 ( D21 ( D b_t0 ) ( T ( D dbdt_t0 ) ) ( T ( D d2bdt2_t0 ) ) ) s0 = + let !( D21 b_t0s0 dbds_t0s0 d2bds2_t0s0 ) = b_t0 s0 + !( D21 dbdt_t0s0 d2bdtds_t0s0 _ ) = dbdt_t0 s0 + !( D21 d2bdt2_t0s0 _ _ ) = d2bdt2_t0 s0 + in D22 + b_t0s0 + ( T dbdt_t0s0 ) dbds_t0s0 + ( T d2bdt2_t0s0 ) d2bdtds_t0s0 d2bds2_t0s0 + +uncurryD3 :: D 3 a ~ D 3 ( ℝ 1 ) + => D 3 ( ℝ 1 ) ( C 3 a b ) -> a -> D 3 ( ℝ 2 ) b +uncurryD3 ( D31 ( D b_t0 ) ( T ( D dbdt_t0 ) ) ( T ( D d2bdt2_t0 ) ) ( T ( D d3bdt3_t0 ) ) ) s0 = + let !( D31 b_t0s0 dbds_t0s0 d2bds2_t0s0 d3bds3_t0s0 ) = b_t0 s0 + !( D31 dbdt_t0s0 d2bdtds_t0s0 d3bdtds2_t0s0 _ ) = dbdt_t0 s0 + !( D31 d2bdt2_t0s0 d3bdt2ds_t0s0 _ _ ) = d2bdt2_t0 s0 + !( D31 d3bdt3_t0s0 _ _ _ ) = d3bdt3_t0 s0 + in D32 + b_t0s0 + ( T dbdt_t0s0 ) dbds_t0s0 + ( T d2bdt2_t0s0 ) d2bdtds_t0s0 d2bds2_t0s0 + ( T d3bdt3_t0s0 ) d3bdt2ds_t0s0 d3bdtds2_t0s0 d3bds3_t0s0 + +-- | Recover the underlying function, discarding all infinitesimal information. +fun :: forall r k v w. HasChainRule r k v => C k v w -> ( v -> w ) +fun ( D df ) = value @r @k @v . df +{-# INLINE fun #-} + +-- | The differentiable germ of a coordinate variable. +var :: forall r k v + . ( Module r ( T r ), Representable r v, HasChainRule r k v ) + => Fin ( RepDim v ) -> C k v r +var i = D $ linearD @r @k @v ( `index` i ) +{-# INLINE var #-} + +-------------------------------------------------------------------------------- + +-- | Newtype for the module instance @Module r v => Module ( dr r ) ( dr v )@. +newtype ApAp r dr v = ApAp { unApAp :: dr v } + +instance ( Ring ( dr r ), Module r ( T v ), Applicative dr ) + => Module ( dr r ) ( ApAp r dr v ) where + ApAp !u ^+^ ApAp !v = ApAp $ liftA2 ( coerce $ (^+^) @r @( T v ) ) u v + ApAp !u ^-^ ApAp !v = ApAp $ liftA2 ( coerce $ (^-^) @r @( T v ) ) u v + origin = ApAp $ pure $ coerce $ origin @r @( T v ) + !k *^ ApAp !u = ApAp $ liftA2 ( coerce $ (*^) @r @( T v ) ) k u + + +deriving via ApAp r D𝔸0 v + instance Module r ( T v ) => Module ( D𝔸0 r ) ( D𝔸0 v ) +--deriving via ApAp r D1𝔸1 v +-- instance Module r ( T v ) => Module ( D1𝔸1 r ) ( D1𝔸1 v ) +deriving via ApAp r D2𝔸1 v + instance Module r ( T v ) => Module ( D2𝔸1 r ) ( D2𝔸1 v ) +deriving via ApAp r D3𝔸1 v + instance Module r ( T v ) => Module ( D3𝔸1 r ) ( D3𝔸1 v ) +--deriving via ApAp r D1𝔸2 v +-- instance Module r ( T v ) => Module ( D1𝔸2 r ) ( D1𝔸2 v ) +deriving via ApAp r D2𝔸2 v + instance Module r ( T v ) => Module ( D2𝔸2 r ) ( D2𝔸2 v ) +deriving via ApAp r D3𝔸2 v + instance Module r ( T v ) => Module ( D3𝔸2 r ) ( D3𝔸2 v ) +--deriving via ApAp r D1𝔸3 v +-- instance Module r ( T v ) => Module ( D1𝔸3 r ) ( D1𝔸3 v ) +deriving via ApAp r D2𝔸3 v + instance Module r ( T v ) => Module ( D2𝔸3 r ) ( D2𝔸3 v ) +deriving via ApAp r D3𝔸3 v + instance Module r ( T v ) => Module ( D3𝔸3 r ) ( D3𝔸3 v ) +--deriving via ApAp r D1𝔸4 v +-- instance Module r ( T v ) => Module ( D1𝔸4 r ) ( D1𝔸4 v ) +deriving via ApAp r D2𝔸4 v + instance Module r ( T v ) => Module ( D2𝔸4 r ) ( D2𝔸4 v ) +deriving via ApAp r D3𝔸4 v + instance Module r ( T v ) => Module ( D3𝔸4 r ) ( D3𝔸4 v ) + +-------------------------------------------------------------------------------- +-- AbelianGroup instances + +newtype ApAp2 k u r = ApAp2 { unApAp2 :: D k u r } + +instance ( Applicative ( D k u ) + , AbelianGroup r + , HasChainRule Double k u + ) => AbelianGroup ( ApAp2 k u r ) where + ApAp2 !x + ApAp2 !y = ApAp2 $ liftA2 ( (+) @r ) x y + ApAp2 !x - ApAp2 !y = ApAp2 $ liftA2 ( (-) @r ) x y + negate ( ApAp2 !x ) = ApAp2 $ fmap ( negate @r ) x + + -- DO NOT USE PURE!! + fromInteger !i = ApAp2 $ konst @Double @k @u ( fromInteger @r i ) + +deriving newtype instance AbelianGroup r => AbelianGroup ( D𝔸0 r ) + +deriving via ApAp2 2 ( ℝ 1 ) r + instance AbelianGroup r => AbelianGroup ( D2𝔸1 r ) +deriving via ApAp2 3 ( ℝ 1 ) r + instance AbelianGroup r => AbelianGroup ( D3𝔸1 r ) +deriving via ApAp2 2 ( ℝ 2 ) r + instance AbelianGroup r => AbelianGroup ( D2𝔸2 r ) +deriving via ApAp2 3 ( ℝ 2 ) r + instance AbelianGroup r => AbelianGroup ( D3𝔸2 r ) +deriving via ApAp2 2 ( ℝ 3 ) r + instance AbelianGroup r => AbelianGroup ( D2𝔸3 r ) +deriving via ApAp2 3 ( ℝ 3 ) r + instance AbelianGroup r => AbelianGroup ( D3𝔸3 r ) +deriving via ApAp2 2 ( ℝ 4 ) r + instance AbelianGroup r => AbelianGroup ( D2𝔸4 r ) +deriving via ApAp2 3 ( ℝ 4 ) r + instance AbelianGroup r => AbelianGroup ( D3𝔸4 r ) + +-------------------------------------------------------------------------------- +-- Ring instances. + +-- TODO: should also implement a power operation specifically, +-- as this is important for interval arithmethic. + +deriving newtype instance Ring r => Ring ( D𝔸0 r ) + +--instance Ring r => Ring ( D1𝔸1 r ) where +-- !dr1 * !dr2 = +-- let +-- o :: r +-- o = fromInteger 0 +-- p :: r -> r -> r +-- p = (+) @r +-- m :: r -> r -> r +-- m = (*) @r +-- in +-- $$( prodRuleQ +-- [|| o ||] [|| p ||] [|| m ||] +-- [|| dr1 ||] [|| dr2 ||] ) + +instance Ring r => Ring ( D2𝔸1 r ) where + !dr1 * !dr2 = + let + o :: r + o = fromInteger 0 + p :: r -> r -> r + p = (+) @r + m :: r -> r -> r + m = (*) @r + in + $$( prodRuleQ + [|| o ||] [|| p ||] [|| m ||] + [|| dr1 ||] [|| dr2 ||] ) + +instance Ring r => Ring ( D3𝔸1 r ) where + !dr1 * !dr2 = + let + o :: r + o = fromInteger 0 + p :: r -> r -> r + p = (+) @r + m :: r -> r -> r + m = (*) @r + in + $$( prodRuleQ + [|| o ||] [|| p ||] [|| m ||] + [|| dr1 ||] [|| dr2 ||] ) + +--instance Ring r => Ring ( D1𝔸2 r ) where +-- !dr1 * !dr2 = +-- let +-- o :: r +-- o = fromInteger 0 +-- p :: r -> r -> r +-- p = (+) @r +-- m :: r -> r -> r +-- m = (*) @r +-- in +-- $$( prodRuleQ +-- [|| o ||] [|| p ||] [|| m ||] +-- [|| dr1 ||] [|| dr2 ||] ) + +instance Ring r => Ring ( D2𝔸2 r ) where + !dr1 * !dr2 = + let + o :: r + o = fromInteger 0 + p :: r -> r -> r + p = (+) @r + m :: r -> r -> r + m = (*) @r + in + $$( prodRuleQ + [|| o ||] [|| p ||] [|| m ||] + [|| dr1 ||] [|| dr2 ||] ) + +instance Ring r => Ring ( D3𝔸2 r ) where + !dr1 * !dr2 = + let + o :: r + o = fromInteger 0 + p :: r -> r -> r + p = (+) @r + m :: r -> r -> r + m = (*) @r + in + $$( prodRuleQ + [|| o ||] [|| p ||] [|| m ||] + [|| dr1 ||] [|| dr2 ||] ) + +--instance Ring r => Ring ( D1𝔸3 r ) where +-- !dr1 * !dr2 = +-- let +-- o :: r +-- o = fromInteger 0 +-- p :: r -> r -> r +-- p = (+) @r +-- m :: r -> r -> r +-- m = (*) @r +-- in +-- $$( prodRuleQ +-- [|| o ||] [|| p ||] [|| m ||] +-- [|| dr1 ||] [|| dr2 ||] ) + +instance Ring r => Ring ( D2𝔸3 r ) where + !dr1 * !dr2 = + let + o :: r + o = fromInteger 0 + p :: r -> r -> r + p = (+) @r + m :: r -> r -> r + m = (*) @r + in + $$( prodRuleQ + [|| o ||] [|| p ||] [|| m ||] + [|| dr1 ||] [|| dr2 ||] ) + +instance Ring r => Ring ( D3𝔸3 r ) where + !dr1 * !dr2 = + let + o :: r + o = fromInteger 0 + p :: r -> r -> r + p = (+) @r + m :: r -> r -> r + m = (*) @r + in + $$( prodRuleQ + [|| o ||] [|| p ||] [|| m ||] + [|| dr1 ||] [|| dr2 ||] ) + +--instance Ring r => Ring ( D1𝔸4 r ) where +-- !dr1 * !dr2 = +-- let +-- o :: r +-- o = fromInteger 0 +-- p :: r -> r -> r +-- p = (+) @r +-- m :: r -> r -> r +-- m = (*) @r +-- in +-- $$( prodRuleQ +-- [|| o ||] [|| p ||] [|| m ||] +-- [|| dr1 ||] [|| dr2 ||] ) + +instance Ring r => Ring ( D2𝔸4 r ) where + !dr1 * !dr2 = + let + o :: r + o = fromInteger 0 + p :: r -> r -> r + p = (+) @r + m :: r -> r -> r + m = (*) @r + in + $$( prodRuleQ + [|| o ||] [|| p ||] [|| m ||] + [|| dr1 ||] [|| dr2 ||] ) + +instance Ring r => Ring ( D3𝔸4 r ) where + !dr1 * !dr2 = + let + o :: r + o = fromInteger 0 + p :: r -> r -> r + p = (+) @r + m :: r -> r -> r + m = (*) @r + in + $$( prodRuleQ + [|| o ||] [|| p ||] [|| m ||] + [|| dr1 ||] [|| dr2 ||] ) + +-------------------------------------------------------------------------------- +-- Field & transcendental instances + +-- TODO!! + +deriving newtype instance Field r => Field ( D𝔸0 r ) +--instance Field r => Field ( D1𝔸1 r ) where +--instance Field r => Field ( D1𝔸2 r ) where +--instance Field r => Field ( D1𝔸3 r ) where +--instance Field r => Field ( D1𝔸4 r ) where +instance Field r => Field ( D2𝔸1 r ) where +instance Field r => Field ( D2𝔸2 r ) where +instance Field r => Field ( D2𝔸3 r ) where +instance Field r => Field ( D2𝔸4 r ) where +instance Field r => Field ( D3𝔸1 r ) where +instance Field r => Field ( D3𝔸2 r ) where +instance Field r => Field ( D3𝔸3 r ) where +instance Field r => Field ( D3𝔸4 r ) where + + +deriving newtype instance Transcendental r => Transcendental ( D𝔸0 r ) +--instance Transcendental r => Transcendental ( D1𝔸1 r ) where +--instance Transcendental r => Transcendental ( D1𝔸2 r ) where +--instance Transcendental r => Transcendental ( D1𝔸3 r ) where +--instance Transcendental r => Transcendental ( D1𝔸4 r ) where +instance Transcendental r => Transcendental ( D2𝔸1 r ) where +instance Transcendental r => Transcendental ( D2𝔸2 r ) where +instance Transcendental r => Transcendental ( D2𝔸3 r ) where + pi = konst @Double @2 @( ℝ 3 ) pi + sin ( D23 v ( T dx ) ( T dy ) ( T dz ) ( T ddx ) ( T dxdy ) ( T ddy ) ( T dxdz ) ( T dydz ) ( T ddz ) ) + = let !s = sin v + !c = cos v + in D23 s + ( T $ c * dx ) ( T $ c * dy ) ( T $ c * dz ) + ( T $ 2 * c * ddx - s * ( dx ^ 2 ) ) + ( T $ 2 * c * dxdy - 2 * s * dx * dy ) + ( T $ 2 * c * ddy - s * ( dy ^ 2 ) ) + ( T $ 2 * c * dxdz - 2 * s * dx * dz ) + ( T $ 2 * c * dydz - 2 * s * dy * dz ) + ( T $ 2 * c * ddz - s * ( dz ^ 2 ) ) + + cos ( D23 v ( T dx ) ( T dy ) ( T dz ) ( T ddx ) ( T dxdy ) ( T ddy ) ( T dxdz ) ( T dydz ) ( T ddz ) ) + = let !s = sin v + !c = cos v + in D23 c + ( T $ -s * dx ) ( T $ -s * dy ) ( T $ -s * dz ) + ( T $ -2 * s * ddx - c * ( dx ^ 2 ) ) + ( T $ -2 * s * dxdy - 2 * c * dx * dy ) + ( T $ -2 * s * ddy - c * ( dy ^ 2 ) ) + ( T $ -2 * s * dxdz - 2 * c * dx * dz ) + ( T $ -2 * s * dydz - 2 * c * dy * dz ) + ( T $ -2 * s * ddz - c * ( dz ^ 2 ) ) + +instance Transcendental r => Transcendental ( D2𝔸4 r ) where +instance Transcendental r => Transcendental ( D3𝔸1 r ) where +instance Transcendental r => Transcendental ( D3𝔸2 r ) where +instance Transcendental r => Transcendental ( D3𝔸3 r ) where +instance Transcendental r => Transcendental ( D3𝔸4 r ) where + +-------------------------------------------------------------------------------- +-- HasChainRule instances. + +instance HasChainRule Double 2 ( ℝ 0 ) where + konst w = D0 w + linearD f v = D0 ( f v ) + value ( D0 v ) = v + chain _ ( D0 gfx ) = D21 gfx origin origin + +instance HasChainRule Double 3 ( ℝ 0 ) where + konst w = D0 w + linearD f v = D0 ( f v ) + value ( D0 v ) = v + chain _ ( D0 gfx ) = D31 gfx origin origin origin + +instance HasChainRule Double 2 ( ℝ 1 ) where + + konst :: forall w. AbelianGroup w => w -> D2𝔸1 w + konst w = + let !o = fromInteger @w 0 + in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] ) + + value df = $$( monIndex [|| df ||] zeroMonomial ) + + linearD :: forall w. Module Double ( T w ) => ( ℝ 1 -> w ) -> ℝ 1 -> D2𝔸1 w + linearD f v = + let !o = origin @Double @( T w ) + in $$( monTabulate \ mon -> + if | isZeroMonomial mon + -> [|| f v ||] + | Just i <- isLinear mon + -> [|| f $$( tabulateQ \ j -> + if | j == i + -> [|| 1 ||] + | otherwise + -> [|| 0 ||] + ) ||] + | otherwise + -> [|| unT o ||] + ) + + chain :: forall w. Module Double ( T w ) => D2𝔸1 ( ℝ 1 ) -> D2𝔸1 w -> D2𝔸1 w + chain !df !dg = + let !o = origin @Double @( T w ) + !p = (^+^) @Double @( T w ) + !s = (^*) @Double @( T w ) + in $$( chainRuleQ + [|| o ||] [|| p ||] [|| s ||] + [|| df ||] [|| dg ||] ) + +instance HasChainRule Double 3 ( ℝ 1 ) where + + konst :: forall w. AbelianGroup w => w -> D3𝔸1 w + konst w = + let !o = fromInteger @w 0 + in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] ) + + value df = $$( monIndex [|| df ||] zeroMonomial ) + + linearD :: forall w. Module Double ( T w ) => ( ℝ 1 -> w ) -> ℝ 1 -> D3𝔸1 w + linearD f v = + let !o = origin @Double @( T w ) + in $$( monTabulate \ mon -> + if | isZeroMonomial mon + -> [|| f v ||] + | Just i <- isLinear mon + -> [|| f $$( tabulateQ \ j -> + if | j == i + -> [|| 1 ||] + | otherwise + -> [|| 0 ||] + ) ||] + | otherwise + -> [|| unT o ||] + ) + + chain :: forall w. Module Double ( T w ) => D3𝔸1 ( ℝ 1 ) -> D3𝔸1 w -> D3𝔸1 w + chain !df !dg = + let !o = origin @Double @( T w ) + !p = (^+^) @Double @( T w ) + !s = (^*) @Double @( T w ) + in $$( chainRuleQ + [|| o ||] [|| p ||] [|| s ||] + [|| df ||] [|| dg ||] ) + +instance HasChainRule Double 2 ( ℝ 2 ) where + + konst :: forall w. AbelianGroup w => w -> D2𝔸2 w + konst w = + let !o = fromInteger @w 0 + in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] ) + + value df = $$( monIndex [|| df ||] zeroMonomial ) + + linearD :: forall w. Module Double ( T w ) => ( ℝ 2 -> w ) -> ℝ 2 -> D2𝔸2 w + linearD f v = + let !o = origin @Double @( T w ) + in $$( monTabulate \ mon -> + if | isZeroMonomial mon + -> [|| f v ||] + | Just i <- isLinear mon + -> [|| f $$( tabulateQ \ j -> + if | j == i + -> [|| 1 ||] + | otherwise + -> [|| 0 ||] + ) ||] + | otherwise + -> [|| unT o ||] + ) + + chain :: forall w. Module Double ( T w ) => D2𝔸1 ( ℝ 2 ) -> D2𝔸2 w -> D2𝔸1 w + chain !df !dg = + let !o = origin @Double @( T w ) + !p = (^+^) @Double @( T w ) + !s = (^*) @Double @( T w ) + in $$( chainRuleQ + [|| o ||] [|| p ||] [|| s ||] + [|| df ||] [|| dg ||] ) + +instance HasChainRule Double 3 ( ℝ 2 ) where + + konst :: forall w. AbelianGroup w => w -> D3𝔸2 w + konst w = + let !o = fromInteger @w 0 + in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] ) + + value df = $$( monIndex [|| df ||] zeroMonomial ) + + linearD :: forall w. Module Double ( T w ) => ( ℝ 2 -> w ) -> ℝ 2 -> D3𝔸2 w + linearD f v = + let !o = origin @Double @( T w ) + in $$( monTabulate \ mon -> + if | isZeroMonomial mon + -> [|| f v ||] + | Just i <- isLinear mon + -> [|| f $$( tabulateQ \ j -> + if | j == i + -> [|| 1 ||] + | otherwise + -> [|| 0 ||] + ) ||] + | otherwise + -> [|| unT o ||] + ) + + chain :: forall w. Module Double ( T w ) => D3𝔸1 ( ℝ 2 ) -> D3𝔸2 w -> D3𝔸1 w + chain !df !dg = + let !o = origin @Double @( T w ) + !p = (^+^) @Double @( T w ) + !s = (^*) @Double @( T w ) + in $$( chainRuleQ + [|| o ||] [|| p ||] [|| s ||] + [|| df ||] [|| dg ||] ) + +instance HasChainRule Double 2 ( ℝ 3 ) where + + konst :: forall w. AbelianGroup w => w -> D2𝔸3 w + konst w = + let !o = fromInteger @w 0 + in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] ) + + value df = $$( monIndex [|| df ||] zeroMonomial ) + + linearD :: forall w. Module Double ( T w ) => ( ℝ 3 -> w ) -> ℝ 3 -> D2𝔸3 w + linearD f v = + let !o = origin @Double @( T w ) + in $$( monTabulate \ mon -> + if | isZeroMonomial mon + -> [|| f v ||] + | Just i <- isLinear mon + -> [|| f $$( tabulateQ \ j -> + if | j == i + -> [|| 1 ||] + | otherwise + -> [|| 0 ||] + ) ||] + | otherwise + -> [|| unT o ||] + ) + + chain :: forall w. Module Double ( T w ) => D2𝔸1 ( ℝ 3 ) -> D2𝔸3 w -> D2𝔸1 w + chain !df !dg = + let !o = origin @Double @( T w ) + !p = (^+^) @Double @( T w ) + !s = (^*) @Double @( T w ) + in $$( chainRuleQ + [|| o ||] [|| p ||] [|| s ||] + [|| df ||] [|| dg ||] ) + +instance HasChainRule Double 3 ( ℝ 3 ) where + + konst :: forall w. AbelianGroup w => w -> D3𝔸3 w + konst w = + let !o = fromInteger @w 0 + in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] ) + + value df = $$( monIndex [|| df ||] zeroMonomial ) + + linearD :: forall w. Module Double ( T w ) => ( ℝ 3 -> w ) -> ℝ 3 -> D3𝔸3 w + linearD f v = + let !o = origin @Double @( T w ) + in $$( monTabulate \ mon -> + if | isZeroMonomial mon + -> [|| f v ||] + | Just i <- isLinear mon + -> [|| f $$( tabulateQ \ j -> + if | j == i + -> [|| 1 ||] + | otherwise + -> [|| 0 ||] + ) ||] + | otherwise + -> [|| unT o ||] + ) + + chain :: forall w. Module Double ( T w ) => D3𝔸1 ( ℝ 3 ) -> D3𝔸3 w -> D3𝔸1 w + chain !df !dg = + let !o = origin @Double @( T w ) + !p = (^+^) @Double @( T w ) + !s = (^*) @Double @( T w ) + in $$( chainRuleQ + [|| o ||] [|| p ||] [|| s ||] + [|| df ||] [|| dg ||] ) + +instance HasChainRule Double 2 ( ℝ 4 ) where + + konst :: forall w. AbelianGroup w => w -> D2𝔸4 w + konst w = + let !o = fromInteger @w 0 + in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] ) + + value df = $$( monIndex [|| df ||] zeroMonomial ) + + linearD :: forall w. Module Double ( T w ) => ( ℝ 4 -> w ) -> ℝ 4 -> D2𝔸4 w + linearD f v = + let !o = origin @Double @( T w ) + in $$( monTabulate \ mon -> + if | isZeroMonomial mon + -> [|| f v ||] + | Just i <- isLinear mon + -> [|| f $$( tabulateQ \ j -> + if | j == i + -> [|| 1 ||] + | otherwise + -> [|| 0 ||] + ) ||] + | otherwise + -> [|| unT o ||] + ) + + chain :: forall w. Module Double ( T w ) => D2𝔸1 ( ℝ 4 ) -> D2𝔸4 w -> D2𝔸1 w + chain !df !dg = + let !o = origin @Double @( T w ) + !p = (^+^) @Double @( T w ) + !s = (^*) @Double @( T w ) + in $$( chainRuleQ + [|| o ||] [|| p ||] [|| s ||] + [|| df ||] [|| dg ||] ) + +instance HasChainRule Double 3 ( ℝ 4 ) where + + konst :: forall w. AbelianGroup w => w -> D3𝔸4 w + konst w = + let !o = fromInteger @w 0 + in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] ) + + value df = $$( monIndex [|| df ||] zeroMonomial ) + + linearD :: forall w. Module Double ( T w ) => ( ℝ 4 -> w ) -> ℝ 4 -> D3𝔸4 w + linearD f v = + let !o = origin @Double @( T w ) + in $$( monTabulate \ mon -> + if | isZeroMonomial mon + -> [|| f v ||] + | Just i <- isLinear mon + -> [|| f $$( tabulateQ \ j -> + if | j == i + -> [|| 1 ||] + | otherwise + -> [|| 0 ||] + ) ||] + | otherwise + -> [|| unT o ||] + ) + + chain :: forall w. Module Double ( T w ) => D3𝔸1 ( ℝ 4 ) -> D3𝔸4 w -> D3𝔸1 w + chain !df !dg = + let !o = origin @Double @( T w ) + !p = (^+^) @Double @( T w ) + !s = (^*) @Double @( T w ) + in $$( chainRuleQ + [|| o ||] [|| p ||] [|| s ||] + [|| df ||] [|| dg ||] ) diff --git a/src/splines/Math/Algebra/Dual/Internal.hs b/src/splines/Math/Algebra/Dual/Internal.hs new file mode 100644 index 0000000..f70b978 --- /dev/null +++ b/src/splines/Math/Algebra/Dual/Internal.hs @@ -0,0 +1,588 @@ + +{-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE QuantifiedConstraints #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE UndecidableInstances #-} + +module Math.Algebra.Dual.Internal + ( D𝔸0(..) + , D1𝔸1(..), D2𝔸1(..), D3𝔸1(..) + , D1𝔸2(..), D2𝔸2(..), D3𝔸2(..) + , D1𝔸3(..), D2𝔸3(..), D3𝔸3(..) + , D1𝔸4(..), D2𝔸4(..), D3𝔸4(..) + + , chainRuleQ + ) where + +-- base +import GHC.Generics + ( Generic, Generic1, Generically1(..) ) +import GHC.TypeNats + ( KnownNat ) + +-- template-haskell +import Language.Haskell.TH + ( CodeQ ) +import Language.Haskell.TH.Syntax + ( Lift(..) ) + +-- MetaBrush +import Math.Linear + ( Vec(..), T(..) + , RepresentableQ(..), RepDim + , zipIndices + ) +import Math.Monomial + ( Mon(..), MonomialBasis(..), Vars, Deg + , mons, faà, multiSubsetsSum, zeroMonomial + ) +import Math.Ring + ( Ring ) +import qualified Math.Ring as Ring +import TH.Utils + ( foldQ, powQ ) + +-------------------------------------------------------------------------------- + +-- | \( \mathbb{Z} \). +newtype D𝔸0 v = D0 { _D0_v :: v } + deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) + deriving Applicative + via Generically1 D𝔸0 + +-- | \( \mathbb{Z}[x]/(x)^2 \). +data D1𝔸1 v = + D11 { _D11_v :: !v + , _D11_dx :: !( T v ) + } + deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) + deriving Applicative + via Generically1 D1𝔸1 + +-- | \( \mathbb{Z}[x]/(x)^3 \). +data D2𝔸1 v = + D21 { _D21_v :: !v + , _D21_dx :: !( T v ) + , _D21_dxdx :: !( T v ) + } + deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) + deriving Applicative + via Generically1 D2𝔸1 + +-- | \( \mathbb{Z}[x]/(x)^4 \). +data D3𝔸1 v = + D31 { _D31_v :: !v + , _D31_dx :: !( T v ) + , _D31_dxdx :: !( T v ) + , _D31_dxdxdx :: !( T v ) + } + deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) + deriving Applicative + via Generically1 D3𝔸1 + +-- | \( \mathbb{Z}[x, y]/(x, y)^2 \). +data D1𝔸2 v = + D12 { _D12_v :: !v + , _D12_dx, _D12_dy :: !( T v ) + } + deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) + deriving Applicative + via Generically1 D1𝔸2 + +-- | \( \mathbb{Z}[x, y]/(x, y)^3 \). +data D2𝔸2 v = + D22 { _D22_v :: !v + , _D22_dx, _D22_dy :: !( T v ) + , _D22_dxdx, _D22_dxdy, _D22_dydy :: !( T v ) + } + deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) + deriving Applicative + via Generically1 D2𝔸2 + +-- | \( \mathbb{Z}[x, y]/(x, y)^4 \). +data D3𝔸2 v = + D32 { _D32_v :: !v + , _D32_dx, _D32_dy :: !( T v ) + , _D32_dxdx, _D32_dxdy, _D32_dydy :: !( T v ) + , _D32_dxdxdx, _D32_dxdxdy, _D32_dxdydy, _D32_dydydy :: !( T v ) + } + deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) + deriving Applicative + via Generically1 D3𝔸2 + +-- | \( \mathbb{Z}[x, y, z]/(x, y, z)^2 \). +data D1𝔸3 v = + D13 { _D13_v :: !v + , _D13_dx, _D13_dy, _D13_dz :: !( T v ) + } + deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) + deriving Applicative + via Generically1 D1𝔸3 + +-- | \( \mathbb{Z}[x, y, z]/(x, y, z)^3 \). +data D2𝔸3 v = + D23 { _D23_v :: !v + , _D23_dx, _D23_dy, _D23_dz :: !( T v ) + , _D23_dxdx, _D23_dxdy, _D23_dydy, _D23_dxdz, _D23_dydz, _D23_dzdz :: !( T v ) + } + deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) + deriving Applicative + via Generically1 D2𝔸3 + +-- | \( \mathbb{Z}[x, y, z]/(x, y, z)^4 \). +data D3𝔸3 v = + D33 { _D33_v :: !v + , _D33_dx, _D33_dy, _D33_dz :: !( T v ) + , _D33_dxdx, _D33_dxdy, _D33_dydy, _D33_dxdz, _D33_dydz, _D33_dzdz :: !( T v ) + , _D33_dxdxdx, _D33_dxdxdy, _D33_dxdydy, _D33_dydydy + , _D33_dxdxdz, _D33_dxdydz, _D33_dxdzdz, _D33_dydydz, _D33_dydzdz, _D33_dzdzdz :: !( T v ) + } + deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) + deriving Applicative + via Generically1 D3𝔸3 + +-- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^2 \). +data D1𝔸4 v = + D14 { _D14_v :: !v + , _D14_dx, _D14_dy, _D14_dz, _D14_dw :: !( T v ) + } + deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) + deriving Applicative + via Generically1 D1𝔸4 + +-- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^3 \). +data D2𝔸4 v = + D24 { _D24_v :: !v + , _D24_dx, _D24_dy, _D24_dz, _D24_dw :: !( T v ) + , _D24_dxdx, _D24_dxdy, _D24_dydy, _D24_dxdz + , _D24_dydz, _D24_dzdz, _D24_dxdw, _D24_dydw, _D24_dzdw, _D24_dwdw :: !( T v ) + } + deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) + deriving Applicative + via Generically1 D2𝔸4 + +-- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^3 \). +data D3𝔸4 v = + D34 { _D34_v :: !v + , _D34_dx, _D34_dy, _D34_dz, _D34_dw :: !( T v ) + , _D34_dxdx, _D34_dxdy, _D34_dydy, _D34_dxdz, _D34_dydz, _D34_dzdz + , _D34_dxdw, _D34_dydw, _D34_dzdw, _D34_dwdw :: !( T v ) + , _D34_dxdxdx, _D34_dxdxdy, _D34_dxdydy, _D34_dydydy, + _D34_dxdxdz, _D34_dxdydz, _D34_dxdzdz, _D34_dydzdz, _D34_dydydz, _D34_dzdzdz, + _D34_dxdxdw, _D34_dxdydw, _D34_dydydw, _D34_dxdzdw, _D34_dydzdw, _D34_dzdzdw, + _D34_dxdwdw, _D34_dydwdw, _D34_dzdwdw, _D34_dwdwdw :: !( T v ) + } + deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) + deriving Applicative + via Generically1 D3𝔸4 + +-------------------------------------------------------------------------------- + +-- | The chain rule, to be spliced in using Template Haskell. +chainRuleQ :: forall dr1 dv v r w d + . ( Ring r, RepresentableQ r v + , MonomialBasis dr1, Vars dr1 ~ 1 + , MonomialBasis dv , Vars dv ~ RepDim v + , Deg dr1 ~ Deg dv + , d ~ Vars dv, KnownNat d + ) + => CodeQ ( T w ) -- Module r ( T w ) + -> CodeQ ( T w -> T w -> T w ) -- + -> CodeQ ( T w -> r -> T w ) -- (circumvent TH constraint issue) + -> CodeQ ( dr1 v ) + -> CodeQ ( dv w ) + -> CodeQ ( dr1 w ) +chainRuleQ zero_w sum_w scale_w df dg = + monTabulate @dr1 \ ( Mon ( k `VS` _ ) ) -> + case k of + -- Set the value of the composition separately, + -- as that isn't handled by the Faà di Bruno formula. + 0 -> monIndex @dv dg zeroMonomial + _ -> + [|| unT $ $$( foldQ sum_w zero_w + [ [|| $$scale_w ( T $$( monIndex @dv dg m_g ) ) + $$( foldQ [|| (Ring.+) ||] [|| Ring.fromInteger ( 0 :: Integer ) ||] + [ [|| Ring.fromInteger $$( liftTyped $ fromIntegral $ faà k is ) Ring.* + $$( foldQ [|| (Ring.*) ||] [|| Ring.fromInteger ( 1 :: Integer ) ||] + [ foldQ [|| (Ring.*) ||] [|| Ring.fromInteger ( 1 :: Integer ) ||] + [ ( indexQ @r @v ( monIndex @dr1 df ( Mon $ f_deg `VS` VZ ) ) v_index ) + `powQ` pow + | ( f_deg, pow ) <- pows + ] + | ( v_index, pows ) <- zipIndices is + ] + ) ||] + | is <- mss + ] + ) + ||] + | m_g <- mons @_ @d k + , let mss = multiSubsetsSum [ 1 .. k ] k ( monDegs m_g ) + , not ( null mss ) -- avoid terms of the form x * 0 + ] + ) ||] + +-------------------------------------------------------------------------------- +-- MonomialBasis instances follow (nothing else). + +type instance Deg D𝔸0 = 0 +type instance Vars D𝔸0 = 0 +instance MonomialBasis D𝔸0 where + monTabulate f = + [|| let + !_D0_v = $$( f $ Mon VZ ) + in D0 { .. } + ||] + + monIndex d _ = [|| _D0_v $$d ||] + +type instance Deg D1𝔸1 = 1 +type instance Vars D1𝔸1 = 1 +instance MonomialBasis D1𝔸1 where + monTabulate f = + [|| let + !_D11_v = $$( f $ Mon ( 0 `VS` VZ ) ) + !_D11_dx = T $$( f $ Mon ( 1 `VS` VZ ) ) + in D11 { .. } + ||] + + monIndex d = \ case + Mon ( 1 `VS` VZ ) -> [|| unT $ _D11_dx $$d ||] + _ -> [|| _D11_v $$d ||] + +type instance Deg D2𝔸1 = 2 +type instance Vars D2𝔸1 = 1 +instance MonomialBasis D2𝔸1 where + monTabulate f = + [|| let + !_D21_v = $$( f $ Mon ( 0 `VS` VZ ) ) + !_D21_dx = T $$( f $ Mon ( 1 `VS` VZ ) ) + !_D21_dxdx = T $$( f $ Mon ( 2 `VS` VZ ) ) + in D21 { .. } + ||] + + monIndex d = \ case + Mon ( 1 `VS` VZ ) -> [|| unT $ _D21_dx $$d ||] + Mon ( 2 `VS` VZ ) -> [|| unT $ _D21_dxdx $$d ||] + _ -> [|| _D21_v $$d ||] + +type instance Deg D3𝔸1 = 3 +type instance Vars D3𝔸1 = 1 +instance MonomialBasis D3𝔸1 where + monTabulate f = + [|| let + !_D31_v = $$( f $ Mon ( 0 `VS` VZ ) ) + !_D31_dx = T $$( f $ Mon ( 1 `VS` VZ ) ) + !_D31_dxdx = T $$( f $ Mon ( 2 `VS` VZ ) ) + !_D31_dxdxdx = T $$( f $ Mon ( 3 `VS` VZ ) ) + in D31 { .. } + ||] + + monIndex d = \ case + Mon ( 1 `VS` VZ ) -> [|| unT $ _D31_dx $$d ||] + Mon ( 2 `VS` VZ ) -> [|| unT $ _D31_dxdx $$d ||] + Mon ( 3 `VS` VZ ) -> [|| unT $ _D31_dxdxdx $$d ||] + _ -> [|| _D31_v $$d ||] + +type instance Deg D1𝔸2 = 1 +type instance Vars D1𝔸2 = 2 +instance MonomialBasis D1𝔸2 where + monTabulate f = + [|| let + !_D12_v = $$( f $ Mon ( 0 `VS` 0 `VS` VZ ) ) + !_D12_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` VZ ) ) + !_D12_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` VZ ) ) + in D12 { .. } + ||] + + monIndex d = \ case + Mon ( 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D12_dx $$d ||] + Mon ( 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D12_dy $$d ||] + _ -> [|| _D12_v $$d ||] + +type instance Deg D2𝔸2 = 2 +type instance Vars D2𝔸2 = 2 +instance MonomialBasis D2𝔸2 where + monTabulate f = + [|| let + !_D22_v = $$( f $ Mon ( 0 `VS` 0 `VS` VZ ) ) + !_D22_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` VZ ) ) + !_D22_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` VZ ) ) + !_D22_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` VZ ) ) + !_D22_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` VZ ) ) + !_D22_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` VZ ) ) + in D22 { .. } + ||] + + monIndex d = \ case + Mon ( 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D22_dx $$d ||] + Mon ( 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D22_dy $$d ||] + Mon ( 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D22_dxdx $$d ||] + Mon ( 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D22_dxdy $$d ||] + Mon ( 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D22_dydy $$d ||] + _ -> [|| _D22_v $$d ||] + +type instance Deg D3𝔸2 = 3 +type instance Vars D3𝔸2 = 2 +instance MonomialBasis D3𝔸2 where + monTabulate f = + [|| let + !_D32_v = $$( f $ Mon ( 0 `VS` 0 `VS` VZ ) ) + !_D32_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` VZ ) ) + !_D32_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` VZ ) ) + !_D32_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` VZ ) ) + !_D32_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` VZ ) ) + !_D32_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` VZ ) ) + !_D32_dxdxdx = T $$( f $ Mon ( 3 `VS` 0 `VS` VZ ) ) + !_D32_dxdxdy = T $$( f $ Mon ( 2 `VS` 1 `VS` VZ ) ) + !_D32_dxdydy = T $$( f $ Mon ( 1 `VS` 2 `VS` VZ ) ) + !_D32_dydydy = T $$( f $ Mon ( 0 `VS` 3 `VS` VZ ) ) + in D32 { .. } ||] + + monIndex d = \ case + Mon ( 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D32_dx $$d ||] + Mon ( 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D32_dy $$d ||] + Mon ( 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D32_dxdx $$d ||] + Mon ( 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D32_dxdy $$d ||] + Mon ( 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D32_dydy $$d ||] + Mon ( 3 `VS` 0 `VS` VZ ) -> [|| unT $ _D32_dxdxdx $$d ||] + Mon ( 2 `VS` 1 `VS` VZ ) -> [|| unT $ _D32_dxdxdy $$d ||] + Mon ( 1 `VS` 2 `VS` VZ ) -> [|| unT $ _D32_dxdydy $$d ||] + Mon ( 0 `VS` 3 `VS` VZ ) -> [|| unT $ _D32_dydydy $$d ||] + _ -> [|| _D32_v $$d ||] + +type instance Deg D1𝔸3 = 1 +type instance Vars D1𝔸3 = 3 +instance MonomialBasis D1𝔸3 where + monTabulate f = + [|| let + !_D13_v = $$( f ( Mon ( 0 `VS` 0 `VS` 0 `VS` VZ ) ) ) + !_D13_dx = T $$( f ( Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) ) ) + !_D13_dy = T $$( f ( Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) ) ) + !_D13_dz = T $$( f ( Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) ) ) + in D13 { .. } + ||] + + monIndex d = \ case + Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D13_dx $$d ||] + Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D13_dy $$d ||] + Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D13_dz $$d ||] + _ -> [|| _D13_v $$d ||] + +type instance Deg D2𝔸3 = 2 +type instance Vars D2𝔸3 = 3 +instance MonomialBasis D2𝔸3 where + monTabulate f = + [|| let + !_D23_v = $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D23_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D23_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) ) + !_D23_dz = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) ) + !_D23_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D23_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) ) + !_D23_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) ) + !_D23_dxdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) ) + !_D23_dydz = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) ) + !_D23_dzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) ) + in D23 { .. } + ||] + + monIndex d = \ case + Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D23_dx $$d ||] + Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D23_dy $$d ||] + Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D23_dz $$d ||] + Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D23_dxdx $$d ||] + Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D23_dxdy $$d ||] + Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D23_dydy $$d ||] + Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D23_dxdz $$d ||] + Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D23_dydz $$d ||] + Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D23_dzdz $$d ||] + _ -> [|| _D23_v $$d ||] + + +type instance Deg D3𝔸3 = 3 +type instance Vars D3𝔸3 = 3 +instance MonomialBasis D3𝔸3 where + monTabulate f = + [|| let + !_D33_v = $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D33_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D33_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) ) + !_D33_dz = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) ) + !_D33_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D33_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) ) + !_D33_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) ) + !_D33_dxdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) ) + !_D33_dydz = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) ) + !_D33_dzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) ) + !_D33_dxdxdx = T $$( f $ Mon ( 3 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D33_dxdxdy = T $$( f $ Mon ( 2 `VS` 1 `VS` 0 `VS` VZ ) ) + !_D33_dxdydy = T $$( f $ Mon ( 1 `VS` 2 `VS` 0 `VS` VZ ) ) + !_D33_dydydy = T $$( f $ Mon ( 0 `VS` 3 `VS` 0 `VS` VZ ) ) + !_D33_dxdxdz = T $$( f $ Mon ( 2 `VS` 0 `VS` 1 `VS` VZ ) ) + !_D33_dxdydz = T $$( f $ Mon ( 1 `VS` 1 `VS` 1 `VS` VZ ) ) + !_D33_dxdzdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 2 `VS` VZ ) ) + !_D33_dydydz = T $$( f $ Mon ( 0 `VS` 2 `VS` 1 `VS` VZ ) ) + !_D33_dydzdz = T $$( f $ Mon ( 0 `VS` 1 `VS` 2 `VS` VZ ) ) + !_D33_dzdzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 3 `VS` VZ ) ) + in D33 { .. } ||] + + monIndex d = \ case + Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dx $$d ||] + Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dy $$d ||] + Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dz $$d ||] + Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dxdx $$d ||] + Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dxdy $$d ||] + Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dydy $$d ||] + Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dxdz $$d ||] + Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dydz $$d ||] + Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D33_dzdz $$d ||] + Mon ( 3 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dxdxdx $$d ||] + Mon ( 2 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dxdxdy $$d ||] + Mon ( 1 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dxdydy $$d ||] + Mon ( 0 `VS` 3 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dydydy $$d ||] + Mon ( 2 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dxdxdz $$d ||] + Mon ( 1 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dxdydz $$d ||] + Mon ( 1 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D33_dxdzdz $$d ||] + Mon ( 0 `VS` 1 `VS` 2 `VS` VZ ) -> [|| unT $ _D33_dydzdz $$d ||] + Mon ( 0 `VS` 0 `VS` 3 `VS` VZ ) -> [|| unT $ _D33_dzdzdz $$d ||] + _ -> [|| _D33_v $$d ||] + +type instance Deg D1𝔸4 = 1 +type instance Vars D1𝔸4 = 4 +instance MonomialBasis D1𝔸4 where + monTabulate f = + [|| let + !_D14_v = $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D14_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D14_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D14_dz = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) ) + !_D14_dw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) ) + in D14 { .. } ||] + + monIndex d = \ case + Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D14_dx $$d ||] + Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D14_dy $$d ||] + Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D14_dz $$d ||] + Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D14_dw $$d ||] + _ -> [|| _D14_v $$d ||] + +type instance Deg D2𝔸4 = 2 +type instance Vars D2𝔸4 = 4 +instance MonomialBasis D2𝔸4 where + monTabulate f = + [|| let + !_D24_v = $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D24_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D24_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D24_dz = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) ) + !_D24_dw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) ) + !_D24_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D24_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D24_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D24_dxdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) ) + !_D24_dydz = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) ) + !_D24_dzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) ) + !_D24_dxdw = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) ) + !_D24_dydw = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) ) + !_D24_dzdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) ) + !_D24_dwdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) ) + in D24 { .. } ||] + + monIndex d = \ case + Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dx $$d ||] + Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dy $$d ||] + Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dz $$d ||] + Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D24_dw $$d ||] + Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dxdx $$d ||] + Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dxdy $$d ||] + Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dydy $$d ||] + Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dxdz $$d ||] + Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dydz $$d ||] + Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dzdz $$d ||] + Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D24_dxdw $$d ||] + Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D24_dydw $$d ||] + Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D24_dzdw $$d ||] + Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D24_dwdw $$d ||] + _ -> [|| _D24_v $$d ||] + + +type instance Deg D3𝔸4 = 3 +type instance Vars D3𝔸4 = 4 +instance MonomialBasis D3𝔸4 where + monTabulate f = + [|| let + !_D34_v = $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D34_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D34_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D34_dz = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) ) + !_D34_dw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) ) + !_D34_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D34_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D34_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D34_dxdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) ) + !_D34_dydz = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) ) + !_D34_dzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) ) + !_D34_dxdw = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) ) + !_D34_dydw = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) ) + !_D34_dzdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) ) + !_D34_dwdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) ) + !_D34_dxdxdx = T $$( f $ Mon ( 3 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D34_dxdxdy = T $$( f $ Mon ( 2 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D34_dxdydy = T $$( f $ Mon ( 1 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D34_dydydy = T $$( f $ Mon ( 0 `VS` 3 `VS` 0 `VS` 0 `VS` VZ ) ) + !_D34_dxdxdz = T $$( f $ Mon ( 2 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) ) + !_D34_dxdydz = T $$( f $ Mon ( 1 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) ) + !_D34_dxdzdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) ) + !_D34_dydydz = T $$( f $ Mon ( 0 `VS` 2 `VS` 1 `VS` 0 `VS` VZ ) ) + !_D34_dydzdz = T $$( f $ Mon ( 0 `VS` 1 `VS` 2 `VS` 0 `VS` VZ ) ) + !_D34_dzdzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 3 `VS` 0 `VS` VZ ) ) + !_D34_dxdxdw = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) ) + !_D34_dxdydw = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) ) + !_D34_dydydw = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 1 `VS` VZ ) ) + !_D34_dxdzdw = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) ) + !_D34_dydzdw = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 1 `VS` VZ ) ) + !_D34_dzdzdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 1 `VS` VZ ) ) + !_D34_dxdwdw = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) ) + !_D34_dydwdw = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 2 `VS` VZ ) ) + !_D34_dzdwdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 2 `VS` VZ ) ) + !_D34_dwdwdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 3 `VS` VZ ) ) + in D34 { .. } ||] + + monIndex d = \ case + Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dx $$d ||] + Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dy $$d ||] + Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dz $$d ||] + Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dw $$d ||] + Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdx $$d ||] + Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdy $$d ||] + Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dydy $$d ||] + Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdz $$d ||] + Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dydz $$d ||] + Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dzdz $$d ||] + Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dxdw $$d ||] + Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dydw $$d ||] + Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dzdw $$d ||] + Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D34_dwdw $$d ||] + Mon ( 3 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdxdx $$d ||] + Mon ( 2 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdxdy $$d ||] + Mon ( 1 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdydy $$d ||] + Mon ( 0 `VS` 3 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dydydy $$d ||] + Mon ( 2 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdxdz $$d ||] + Mon ( 1 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdydz $$d ||] + Mon ( 1 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdzdz $$d ||] + Mon ( 0 `VS` 2 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dydydz $$d ||] + Mon ( 0 `VS` 1 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dydzdz $$d ||] + Mon ( 0 `VS` 0 `VS` 3 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dzdzdz $$d ||] + Mon ( 2 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dxdxdw $$d ||] + Mon ( 1 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dxdydw $$d ||] + Mon ( 0 `VS` 2 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dydydw $$d ||] + Mon ( 1 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dxdzdw $$d ||] + Mon ( 0 `VS` 1 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dydzdw $$d ||] + Mon ( 0 `VS` 0 `VS` 2 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dzdzdw $$d ||] + Mon ( 1 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D34_dxdwdw $$d ||] + Mon ( 0 `VS` 1 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D34_dydwdw $$d ||] + Mon ( 0 `VS` 0 `VS` 1 `VS` 2 `VS` VZ ) -> [|| unT $ _D34_dzdwdw $$d ||] + Mon ( 0 `VS` 0 `VS` 0 `VS` 3 `VS` VZ ) -> [|| unT $ _D34_dwdwdw $$d ||] + _ -> [|| _D34_v $$d ||] diff --git a/src/splines/Math/Bezier/Cubic.hs b/src/splines/Math/Bezier/Cubic.hs index 784cf35..b40a76f 100644 --- a/src/splines/Math/Bezier/Cubic.hs +++ b/src/splines/Math/Bezier/Cubic.hs @@ -63,6 +63,7 @@ import Math.Roots ( realRoots, solveQuadratic ) import Math.Linear ( ℝ(..), T(..) ) +import qualified Math.Ring as Ring -------------------------------------------------------------------------------- @@ -107,13 +108,13 @@ bezier ( Bezier {..} ) t = -- | Derivative of a cubic Bézier curve. bezier' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> v bezier' ( Bezier {..} ) - = ( 3 *^ ) + = ( Ring.fromInteger 3 *^ ) . Quadratic.bezier @v ( Quadratic.Bezier ( p0 --> p1 ) ( p1 --> p2 ) ( p2 --> p3 ) ) -- | Second derivative of a cubic Bézier curve. bezier'' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> v bezier'' ( Bezier {..} ) t - = ( 6 *^ ) + = ( Ring.fromInteger 6 *^ ) $ lerp @v t ( p1 --> p0 ^+^ p1 --> p2 ) ( p2 --> p1 ^+^ p2 --> p3 ) diff --git a/src/splines/Math/Bezier/Quadratic.hs b/src/splines/Math/Bezier/Quadratic.hs index 601e2f7..c8c5c72 100644 --- a/src/splines/Math/Bezier/Quadratic.hs +++ b/src/splines/Math/Bezier/Quadratic.hs @@ -60,6 +60,7 @@ import Math.Roots ( realRoots ) import Math.Linear ( ℝ(..), T(..) ) +import qualified Math.Ring as Ring -------------------------------------------------------------------------------- @@ -91,11 +92,11 @@ bezier ( Bezier {..} ) t = lerp @v t ( lerp @v t p0 p1 ) ( lerp @v t p1 p2 ) -- | Derivative of a quadratic Bézier curve. bezier' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> v -bezier' ( Bezier {..} ) t = 2 *^ lerp @v t ( p0 --> p1 ) ( p1 --> p2 ) +bezier' ( Bezier {..} ) t = Ring.fromInteger 2 *^ lerp @v t ( p0 --> p1 ) ( p1 --> p2 ) -- | Second derivative of a quadratic Bézier curve. bezier'' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> v -bezier'' ( Bezier {..} ) = 2 *^ ( p1 --> p0 ^+^ p1 --> p2 ) +bezier'' ( Bezier {..} ) = Ring.fromInteger 2 *^ ( p1 --> p0 ^+^ p1 --> p2 ) -- | Curvature of a quadratic Bézier curve. curvature :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> r -> r diff --git a/src/splines/Math/Bezier/Stroke.hs b/src/splines/Math/Bezier/Stroke.hs index 6b6faa1..da0f098 100644 --- a/src/splines/Math/Bezier/Stroke.hs +++ b/src/splines/Math/Bezier/Stroke.hs @@ -100,6 +100,7 @@ import Control.Monad.Trans.Writer.CPS ( WriterT, execWriterT, runWriter, tell ) -- MetaBrush +import Math.Algebra.Dual import qualified Math.Bezier.Cubic as Cubic import Math.Bezier.Cubic.Fit ( FitPoint, FitParameters, fitSpline ) @@ -114,8 +115,12 @@ import Math.Bezier.Spline , showSplinePoints ) import qualified Math.Bezier.Quadratic as Quadratic +import Math.Differentiable + ( Differentiable, DiffInterp ) import Math.Epsilon ( epsilon ) +import Math.Interval +import Math.Linear import Math.Module ( Module(..), Inner((^.^)), Cross(cross), Interpolatable , lerp, convexCombination, strictlyParallel @@ -126,8 +131,7 @@ import Math.Orientation ) import Math.Roots ( solveQuadratic, newtonRaphson ) -import Math.Linear -import Math.Linear.Dual + import Debug.Utils @@ -195,12 +199,15 @@ computeStrokeOutline :: , NFData ptData, NFData crvData -- Differentiability. - , Differentiable 'Point brushParams - , Differentiable 'Interval brushParams + , DiffInterp 'Point brushParams + , DiffInterp 'Interval brushParams , Interpolatable Double usedParams , Interpolatable ( 𝕀 Double ) ( 𝕀 usedParams ) - , Diffy Double usedParams - , Diffy ( 𝕀 Double ) ( 𝕀 usedParams ) + , HasChainRule Double 2 brushParams + , HasChainRule ( 𝕀 Double ) 2 ( 𝕀 brushParams ) + , HasChainRule Double 2 usedParams + , HasChainRule ( 𝕀 Double ) 2 ( 𝕀 usedParams ) + , Traversable ( D 2 brushParams ) -- Debugging. , Show ptData, Show brushParams @@ -209,7 +216,8 @@ computeStrokeOutline :: => FitParameters -> ( ptData -> usedParams ) -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing - -> ( forall i. Differentiable i brushParams + -> ( forall i + . DiffInterp i brushParams => Proxy# i -> ( forall a. a -> I i a ) -> I i brushParams ~> Spline Closed () ( I i ( ℝ 2 ) ) @@ -323,7 +331,7 @@ computeStrokeOutline fitParams ptParams toBrushParams brushFn spline@( Spline { outlineFunction ptParams toBrushParams brushFn p0 crv :<| go ( openCurveEnd crv ) crvs brushShape :: ptData -> SplinePts Closed - brushShape pt = fun ( brushFn @Point proxy# id ) $ toBrushParams $ ptParams pt + brushShape pt = fun @Double ( brushFn @Point proxy# id ) $ toBrushParams $ ptParams pt updateSpline :: ( T ( ℝ 2 ), T ( ℝ 2 ), T ( ℝ 2 ) ) -> ST s OutlineData updateSpline ( lastTgt, lastTgtFwd, lastTgtBwd ) @@ -438,19 +446,26 @@ outlineFunction . ( HasType ( ℝ 2 ) ptData -- Differentiability. - , Differentiable 'Point brushParams - , Differentiable 'Interval brushParams , Interpolatable Double usedParams , Interpolatable ( 𝕀 Double ) ( 𝕀 usedParams ) - , Diffy Double usedParams - , Diffy ( 𝕀 Double ) ( 𝕀 usedParams ) + , DiffInterp 'Point brushParams + , DiffInterp 'Interval brushParams + , HasChainRule Double 2 usedParams + , HasChainRule ( 𝕀 Double ) 2 ( 𝕀 usedParams ) + , HasChainRule Double 2 brushParams + , HasChainRule ( 𝕀 Double ) 2 ( 𝕀 brushParams ) + , Traversable ( D 2 brushParams ) + + -- , Diffy Double usedParams + -- , Diffy ( 𝕀 Double ) ( 𝕀 usedParams ) -- Debugging. , Show ptData, Show brushParams ) => ( ptData -> usedParams ) -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing - -> ( forall i. Differentiable i brushParams + -> ( forall i + . DiffInterp i brushParams => Proxy# i -> ( forall a. a -> I i a ) -> I i brushParams ~> Spline Closed () ( I i ( ℝ 2 ) ) @@ -461,7 +476,7 @@ outlineFunction outlineFunction ptParams toBrushParams brushFromParams sp0 crv = let pathAndUsedParams :: forall i - . ( D ( I i ( ℝ 1 ) ) ~ D ( ℝ 1 ) + . ( D 2 ( I i ( ℝ 1 ) ) ~ D 2 ( ℝ 1 ) , Coercible ( I i ( ℝ 1 ) ) ( I i Double ) , Module ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) ) @@ -492,7 +507,7 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv = curvesI :: 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 'Interval ) curvesI = brushStrokeData @'Interval @brushParams pathI - ( usedParamsI `chainRule` linear ( nonDecreasing toBrushParams ) ) + ( chainRule @( 𝕀 Double ) @2 usedParamsI $ linear ( nonDecreasing toBrushParams ) ) ( brushFromParams @'Interval proxy# singleton ) usedParamsI :: 𝕀ℝ 1 ~> 𝕀 usedParams @@ -501,7 +516,7 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv = fwdBwd :: OutlineFn fwdBwd t - = solveEnvelopeEquations path_t path'_t ( fwdOffset, bwdOffset ) curves + = solveEnvelopeEquations t path_t path'_t ( fwdOffset, bwdOffset ) curves -- = ( ( offset fwdOffset • path_t, path'_t ) -- , ( offset bwdOffset • path_t, -1 *^ path'_t ) ) where @@ -509,26 +524,26 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv = curves :: Seq ( ℝ 1 -> StrokeDatum Point ) curves = brushStrokeData @Point @brushParams path - ( usedParams `chainRule` linear toBrushParams ) + ( chainRule @Double usedParams $ linear toBrushParams ) ( brushFromParams @Point proxy# id ) t fwdOffset = withTangent path'_t brush_t bwdOffset = withTangent ( -1 *^ path'_t ) brush_t - D1 path_t path'_t _ = runD path t - D1 params_t _ _ = runD usedParams t - brush_t = value @Double @brushParams - $ runD ( brushFromParams @Point proxy# id ) - $ toBrushParams params_t + D21 path_t path'_t _ = runD path t + D21 params_t _ _ = runD usedParams t + brush_t = value @Double @2 @brushParams + $ runD ( brushFromParams @Point proxy# id ) + $ toBrushParams params_t bisSols = bisection 0.0001 curvesI - in trace - ( unlines $ - ( "bisectionMethod: #(possible zeroes) = " ++ show ( length bisSols ) ) : - "" : - map show bisSols ) + in --trace + -- ( unlines $ + -- ( "bisectionMethod: #(possible zeroes) = " ++ show ( length bisSols ) ) : + -- "" : + -- map show bisSols ) fwdBwd ----------------------------------- @@ -821,11 +836,11 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) -- - \( p(t) \) is the path that the brush follows, and -- - \( b(t,s) \) is the brush shape, as it varies along the path. brushStroke :: Module r ( T v ) - => D ( ℝ 1 ) v -- ^ stroke path \( p(t) \) - -> D ( ℝ 2 ) v -- ^ brush \( b(t,s) \) - -> D ( ℝ 2 ) v -brushStroke ( D1 p dpdt d2pdt2 ) ( D2 b dbdt dbds d2bdt2 d2bdtds d2bds2 ) = - D2 ( unT $ T p ^+^ T b ) + => D 2 ( ℝ 1 ) v -- ^ stroke path \( p(t) \) + -> D 2 ( ℝ 2 ) v -- ^ brush \( b(t,s) \) + -> D 2 ( ℝ 2 ) v +brushStroke ( D21 p dpdt d2pdt2 ) ( D22 b dbdt dbds d2bdt2 d2bdtds d2bds2 ) = + D22 ( unT $ T p ^+^ T b ) -- c = p + b ( dpdt ^+^ dbdt ) dbds @@ -851,13 +866,13 @@ brushStroke ( D1 p dpdt d2pdt2 ) ( D2 b dbdt dbds d2bdt2 d2bdtds d2bds2 ) = -- -- NB: if \( \frac{\partial E}{\partial s} \) is zero, the total derivative is ill-defined. envelopeEquation :: forall i - . ( D ( I i ( ℝ 2 ) ) ~ D ( ℝ 2 ) + . ( D 2 ( I i ( ℝ 2 ) ) ~ D 2 ( ℝ 2 ) , Cross ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) , Fractional ( I i Double ) ) - => D ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) + => D 2 ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) -> ( I i Double, T ( I i ( ℝ 2 ) ), T ( I i ( ℝ 2 ) ), I i Double, I i Double ) -envelopeEquation ( D2 _ dcdt dcds d2cdt2 d2cdtds d2cds2 ) = +envelopeEquation ( D22 _ dcdt dcds d2cdt2 d2cdtds d2cds2 ) = let ee = dcdt `cross` dcds dEdt = d2cdt2 `cross` dcds + dcdt `cross` d2cdtds dEds = d2cdtds `cross` dcds + dcdt `cross` d2cds2 @@ -874,41 +889,41 @@ envelopeEquation ( D2 _ dcdt dcds d2cdt2 d2cdtds d2cds2 ) = -- | Linear interpolation, as a differentiable function. line :: forall i b . ( Module ( I i Double ) ( T b ), Torsor ( T b ) b - , D ( I i ( ℝ 1 ) ) ~ D ( ℝ 1 ) + , D 2 ( I i ( ℝ 1 ) ) ~ D 2 ( ℝ 1 ) , Coercible ( I i ( ℝ 1 ) ) ( I i Double ) ) => Segment b -> I i ( ℝ 1 ) ~> b line ( Segment a b ) = D \ ( coerce -> t ) -> - D1 ( lerp @( T b ) t a b ) - ( a --> b ) - origin + D21 ( lerp @( T b ) t a b ) + ( a --> b ) + origin -- | A quadratic Bézier curve, as a differentiable function. bezier2 :: forall i b . ( Module ( I i Double ) ( T b ), Torsor ( T b ) b - , D ( I i ( ℝ 1 ) ) ~ D ( ℝ 1 ) + , D 2 ( I i ( ℝ 1 ) ) ~ D 2 ( ℝ 1 ) , Coercible ( I i ( ℝ 1 ) ) ( I i Double ) ) => Quadratic.Bezier b -> I i ( ℝ 1 ) ~> b bezier2 bez = D \ ( coerce -> t ) -> - D1 ( Quadratic.bezier @( T b ) bez t ) - ( Quadratic.bezier' bez t ) - ( Quadratic.bezier'' bez ) + D21 ( Quadratic.bezier @( T b ) bez t ) + ( Quadratic.bezier' bez t ) + ( Quadratic.bezier'' bez ) -- | A cubic Bézier curve, as a differentiable function. bezier3 :: forall i b . ( Module ( I i Double ) ( T b ), Torsor ( T b ) b - , D ( I i ( ℝ 1 ) ) ~ D ( ℝ 1 ) + , D 2 ( I i ( ℝ 1 ) ) ~ D 2 ( ℝ 1 ) , Coercible ( I i ( ℝ 1 ) ) ( I i Double ) ) => Cubic.Bezier b -> I i ( ℝ 1 ) ~> b bezier3 bez = D \ ( coerce -> t ) -> - D1 ( Cubic.bezier @( T b ) bez t ) - ( Cubic.bezier' bez t ) - ( Cubic.bezier'' bez t ) + D21 ( Cubic.bezier @( T b ) bez t ) + ( Cubic.bezier' bez t ) + ( Cubic.bezier'' bez t ) splineCurveFns :: forall i - . ( D ( I i ( ℝ 1 ) ) ~ D ( ℝ 1 ) + . ( D 2 ( I i ( ℝ 1 ) ) ~ D 2 ( ℝ 1 ) , Module ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) ) , Coercible ( I i ( ℝ 1 ) ) ( I i Double ) ) @@ -934,12 +949,13 @@ splineCurveFns spls -- | Solve the envelope equations at a given point \( t = t_0 \), to find -- \( s_0 \) such that \( c(t_0, s_0) \) is on the envelope of the brush stroke. -solveEnvelopeEquations :: ℝ 2 +solveEnvelopeEquations :: ℝ 1 -- ^ @t@ (for debugging only) + -> ℝ 2 -> T ( ℝ 2 ) -> ( Offset, Offset ) -> Seq ( ℝ 1 -> StrokeDatum Point ) -> ( ( ℝ 2, T ( ℝ 2 ) ), ( ℝ 2, T ( ℝ 2 ) ) ) -solveEnvelopeEquations path_t path'_t ( fwdOffset, bwdOffset ) strokeData +solveEnvelopeEquations _t path_t path'_t ( fwdOffset, bwdOffset ) strokeData = ( fwdSol, ( bwdPt, -1 *^ bwdTgt ) ) where @@ -947,6 +963,7 @@ solveEnvelopeEquations path_t path'_t ( fwdOffset, bwdOffset ) strokeData -- !_ = trace -- ( unlines -- [ "solveEnvelopeEquation" + -- , " t: " ++ show _t -- , " pt: " ++ show path_t -- , " tgt: " ++ show path'_t -- , " fwdOffset: " ++ show fwdOffset @@ -1016,7 +1033,7 @@ solveEnvelopeEquations path_t path'_t ( fwdOffset, bwdOffset ) strokeData in --trace -- ( unlines -- [ "solveEnvelopeEquations" - -- , " t = " ++ show t + -- , " t = " ++ show _t -- , " s = " ++ show s -- , " c = " ++ show dstroke -- , " E = " ++ show _ee @@ -1024,7 +1041,8 @@ solveEnvelopeEquations path_t path'_t ( fwdOffset, bwdOffset ) strokeData -- , " ∂E/∂s = " ++ show 𝛿E𝛿s -- , " dc/dt = " ++ show totDeriv -- ] ) - ( good, ℝ1 s, value @Double @( ℝ 2 ) dstroke, totDeriv ) + ( good, ℝ1 s, value @Double @2 @( ℝ 2 ) dstroke + , totDeriv ) eqn :: ( ℝ 1 -> StrokeDatum Point ) -> ( Double -> ( Double, Double ) ) eqn f s = @@ -1046,10 +1064,10 @@ instance Applicative ZipSeq where liftA2 f ( ZipSeq xs ) ( ZipSeq ys ) = ZipSeq ( Seq.zipWith f xs ys ) brushStrokeData :: forall i brushParams - . ( Diffy ( I i Double ) ( I i brushParams ) + . ( Differentiable i brushParams , Fractional ( I i Double ) - , D ( I i ( ℝ 1 ) ) ~ D ( ℝ 1 ) - , D ( I i ( ℝ 2 ) ) ~ D ( ℝ 2 ) + , D 2 ( I i ( ℝ 1 ) ) ~ D 2 ( ℝ 1 ) + , D 2 ( I i ( ℝ 2 ) ) ~ D 2 ( ℝ 2 ) , Cross ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) ) , Coercible ( I i ( ℝ 1 ) ) ( I i Double ) @@ -1065,26 +1083,27 @@ brushStrokeData :: forall i brushParams brushStrokeData path params brush = \ t -> let - dpath_t :: D ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) + dpath_t :: D 2 ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) !dpath_t = runD path t - dparams_t :: D ( I i ( ℝ 1 ) ) ( I i brushParams ) - !dparams_t@( D1 { v = params_t } ) = runD params t - dbrush_params :: D ( I i brushParams ) ( Spline Closed () ( I i ( ℝ 2 ) ) ) + dparams_t :: D 2 ( I i ( ℝ 1 ) ) ( I i brushParams ) + !dparams_t@( D21 { _D21_v = params_t } ) = runD params t + dbrush_params :: D 2 ( I i brushParams ) ( Spline Closed () ( I i ( ℝ 2 ) ) ) !dbrush_params = runD brush params_t - splines :: Seq ( D ( I i brushParams ) ( I i ( ℝ 1 ) ~> I i ( ℝ 2 ) ) ) + splines :: Seq ( D 2 ( I i brushParams ) ( I i ( ℝ 1 ) ~> I i ( ℝ 2 ) ) ) !splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns @i ) dbrush_params - dbrushes_t :: Seq ( I i ( ℝ 1 ) -> D ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) ) - !dbrushes_t = force $ fmap ( uncurryD . ( dparams_t `chain` ) ) splines + dbrushes_t :: Seq ( I i ( ℝ 1 ) -> D 2 ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) ) + !dbrushes_t = force $ fmap ( uncurryD2 . ( chain @(I i Double) @2 dparams_t ) ) splines + -- This is the crucial use of the chain rule. in fmap ( mkStrokeDatum dpath_t ) dbrushes_t where - mkStrokeDatum :: D ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) - -> ( I i ( ℝ 1 ) -> D ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) ) + mkStrokeDatum :: D 2 ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) + -> ( I i ( ℝ 1 ) -> D 2 ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) ) -> ( I i ( ℝ 1 ) -> StrokeDatum i ) mkStrokeDatum dpath_t dbrush_t s = let dbrush_t_s = dbrush_t s - dstroke@( D2 _c _𝛿c𝛿t _𝛿c𝛿s _ _ _ ) = brushStroke dpath_t dbrush_t_s + dstroke@( D22 _c _𝛿c𝛿t _𝛿c𝛿s _ _ _ ) = brushStroke dpath_t dbrush_t_s ( ee, dcdt, 𝛿E𝛿sdcdt, 𝛿E𝛿t, 𝛿E𝛿s ) = envelopeEquation @i dstroke in -- trace -- ( unlines @@ -1111,14 +1130,14 @@ brushStrokeData path params brush = data StrokeDatum i = StrokeDatum { -- | Path \( p(t_0) \) (with its 1st and 2nd derivatives). - dpath :: D ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) + dpath :: D 2 ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) -- | Brush shape \( b(t_0, s_0) \) (with its 1st and 2nd derivatives). - , dbrush :: D ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) + , dbrush :: D 2 ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) -- Everything below can be computed in terms of the first two fields. -- | Stroke \( c(t_0,s_0) = p(t_0) + b(t_0,s_0) \) (with its 1st and 2nd derivatives). - , dstroke :: D ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) + , dstroke :: D 2 ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) -- | Envelope -- -- \[ E(t_0,s_0) = \left ( \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} \right )_{(t_0,s_0)}. \] @@ -1209,6 +1228,6 @@ bisection minWidth eqs = bisect initialCands [] [] && cmpℝ2 (>) ( getRounded ( Interval.sup 𝛿E𝛿sdcdt ) ) ( ℝ2 0 0 ) return ( ee, 𝛿E𝛿sdcdt ) -cmpℝ2 :: (Double -> Double -> Bool) -> ℝ 2 -> ℝ 2 -> Bool +cmpℝ2 :: ( Double -> Double -> Bool ) -> ℝ 2 -> ℝ 2 -> Bool cmpℝ2 cmp ( ℝ2 x1 y1 ) ( ℝ2 x2 y2 ) = cmp x1 x2 && cmp y1 y2 diff --git a/src/splines/Math/Differentiable.hs b/src/splines/Math/Differentiable.hs new file mode 100644 index 0000000..9c77026 --- /dev/null +++ b/src/splines/Math/Differentiable.hs @@ -0,0 +1,66 @@ +{-# LANGUAGE UndecidableInstances #-} + +module Math.Differentiable + ( Differentiable, DiffInterp ) + where + +-- base +import Data.Kind + ( Type, Constraint ) +import GHC.TypeNats + ( Nat ) + +-- MetaBrush +import Math.Algebra.Dual + ( D, HasChainRule ) +import Math.Linear +import Math.Module +import Math.Interval + ( Extent(..), type I ) +import Math.Ring + ( Transcendental ) + +-------------------------------------------------------------------------------- + +type ExtentOrder :: Extent -> Nat +type family ExtentOrder e where + ExtentOrder i = 2 + --ExtentOrder 'Point = 2 + --ExtentOrder 'Interval = 2 +-- Currently we're doing order 2 derivatives for the brush stroke fitting, +-- but order 3 derivatives for the interval Newton method to find cusps. +-- TODO: using 2 for both until migration finishes. + +type Differentiable :: Extent -> Type -> Constraint +class + ( Module ( I i Double ) ( T ( I i u ) ) + , HasChainRule ( I i Double ) ( ExtentOrder i ) ( I i u ) + , Traversable ( D ( ExtentOrder i ) ( I i u ) ) + ) => Differentiable i u +instance + ( Module ( I i Double ) ( T ( I i u ) ) + , HasChainRule ( I i Double ) ( ExtentOrder i ) ( I i u ) + , Traversable ( D ( ExtentOrder i ) ( I i u ) ) + ) => Differentiable i u + +type DiffInterp :: Extent -> Type -> Constraint +class + ( Differentiable i u + , Interpolatable ( I i Double ) ( I i u ) + , Module ( I i Double ) ( T ( I i Double ) ) + , Module ( D ( ExtentOrder i ) ( I i u ) ( I i Double ) ) + ( D ( ExtentOrder i ) ( I i u ) ( I i ( ℝ 2 ) ) ) + , Transcendental ( D ( ExtentOrder i ) ( I i u ) ( I i Double ) ) + , Applicative ( D ( ExtentOrder i ) ( I i u ) ) + , Representable ( I i Double ) ( I i u ) + ) => DiffInterp i u +instance + ( Differentiable i u + , Interpolatable ( I i Double ) ( I i u ) + , Module ( I i Double ) ( T ( I i Double ) ) + , Module ( D ( ExtentOrder i ) ( I i u ) ( I i Double ) ) + ( D ( ExtentOrder i ) ( I i u ) ( I i ( ℝ 2 ) ) ) + , Transcendental ( D ( ExtentOrder i ) ( I i u ) ( I i Double ) ) + , Applicative ( D ( ExtentOrder i ) ( I i u ) ) + , Representable ( I i Double ) ( I i u ) + ) => DiffInterp i u diff --git a/src/splines/Math/Interval.hs b/src/splines/Math/Interval.hs new file mode 100644 index 0000000..8265088 --- /dev/null +++ b/src/splines/Math/Interval.hs @@ -0,0 +1,441 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE UndecidableInstances #-} + +{-# OPTIONS_GHC -Wno-orphans #-} + +module Math.Interval + ( 𝕀, 𝕀ℝ + , Extent(..), type I + , singleton, nonDecreasing + ) + where + +-- base +import Prelude hiding ( Num(..) ) +import Data.Coerce + ( coerce ) +import Data.Kind + ( Type ) +import Data.Monoid + ( Sum(..) ) + +-- acts +import Data.Act + ( Act((•)), Torsor((-->)) ) + +-- groups +import Data.Group + ( Group(..) ) + +-- groups-generic +import Data.Group.Generics + ( ) + +-- rounded-hw +import Numeric.Rounded.Hardware + ( Rounded(..) ) +import Numeric.Rounded.Hardware.Interval.NonEmpty + ( Interval(..) ) + +-- splines +import Math.Algebra.Dual +import Math.Algebra.Dual.Internal + ( chainRuleQ ) +import Math.Interval.Internal + ( type 𝕀 ) +import Math.Linear + ( ℝ(..), T(..) + , RepresentableQ(..) + ) +import Math.Module +import Math.Monomial +import Math.Ring + +-------------------------------------------------------------------------------- +-- Interval arithmetic using rounded-hw library. + +type 𝕀ℝ n = 𝕀 ( ℝ n ) +type instance D k ( 𝕀 v ) = D k v + +-- Handling points and intervals uniformly. +data Extent = Point | Interval + +type I :: Extent -> Type -> Type +type family I i a where + I 'Point a = a + I 'Interval a = 𝕀 a + +singleton :: a -> 𝕀 a +singleton a = I ( Rounded a ) ( Rounded a ) + +-- | Turn a non-decreasing function into a function on intervals. +nonDecreasing :: ( a -> b ) -> 𝕀 a -> 𝕀 b +nonDecreasing f ( I ( Rounded lo ) ( Rounded hi ) ) = + I ( Rounded $ f lo ) ( Rounded $ f hi ) + + +deriving via ViaAbelianGroup ( T ( 𝕀 Double ) ) + instance Semigroup ( T ( 𝕀 Double ) ) +deriving via ViaAbelianGroup ( T ( 𝕀 Double ) ) + instance Monoid ( T ( 𝕀 Double ) ) +deriving via ViaAbelianGroup ( T ( 𝕀 Double ) ) + instance Group ( T ( 𝕀 Double ) ) + +instance Act ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where + T g • a = coerce ( Sum g • a ) +instance Torsor ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where + a --> b = T $ getSum ( a --> b ) + +------------------------------------------------------------------------------- + +instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ 0 ) ) where + origin = T $$( tabulateQ \ _ -> [|| unT $ origin ||] ) + T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] ) + T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] ) + k *^ T a = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] ) + +instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ 1 ) ) where + origin = T $$( tabulateQ \ _ -> [|| unT $ origin ||] ) + T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] ) + T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] ) + k *^ T a = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] ) + +instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ 2 ) ) where + origin = T $$( tabulateQ \ _ -> [|| unT $ origin ||] ) + T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] ) + T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] ) + k *^ T a = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] ) + +instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ 3 ) ) where + origin = T $$( tabulateQ \ _ -> [|| unT $ origin ||] ) + T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] ) + T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] ) + k *^ T a = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] ) + +instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ 4 ) ) where + origin = T $$( tabulateQ \ _ -> [|| unT $ origin ||] ) + T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] ) + T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] ) + k *^ T a = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] ) + +instance Inner ( 𝕀 Double ) ( T ( 𝕀ℝ 2 ) ) where + T ( I ( Rounded ( ℝ2 x1_lo y1_lo ) ) ( Rounded ( ℝ2 x1_hi y1_hi ) ) ) ^.^ + T ( I ( Rounded ( ℝ2 x2_lo y2_lo ) ) ( Rounded ( ℝ2 x2_hi y2_hi ) ) ) + = let !x1x2 = I ( Rounded x1_lo ) ( Rounded x1_hi ) * I ( Rounded x2_lo ) ( Rounded x2_hi ) + !y1y2 = I ( Rounded y1_lo ) ( Rounded y1_hi ) * I ( Rounded y2_lo ) ( Rounded y2_hi ) + in x1x2 + y1y2 + +instance Cross ( 𝕀 Double ) ( T ( 𝕀ℝ 2 ) ) where + T ( I ( Rounded ( ℝ2 x1_lo y1_lo ) ) ( Rounded ( ℝ2 x1_hi y1_hi ) ) ) `cross` + T ( I ( Rounded ( ℝ2 x2_lo y2_lo ) ) ( Rounded ( ℝ2 x2_hi y2_hi ) ) ) + = let !x1y2 = I ( Rounded x1_lo ) ( Rounded x1_hi ) * I ( Rounded y2_lo ) ( Rounded y2_hi ) + !y2x1 = I ( Rounded x2_lo ) ( Rounded x2_hi ) * I ( Rounded y1_lo ) ( Rounded y1_hi ) + in x1y2 - y2x1 + +deriving via ViaModule ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) + instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) => Semigroup ( T ( 𝕀ℝ n ) ) +deriving via ViaModule ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) + instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) => Monoid ( T ( 𝕀ℝ n ) ) +deriving via ViaModule ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) + instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) => Group ( T ( 𝕀ℝ n ) ) + +deriving via ViaModule ( 𝕀 Double ) ( 𝕀ℝ n ) + instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) => Act ( T ( 𝕀ℝ n ) ) ( 𝕀ℝ n ) +deriving via ( ViaModule ( 𝕀 Double ) ( 𝕀ℝ n ) ) + instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) => Torsor ( T ( 𝕀ℝ n ) ) ( 𝕀ℝ n ) + +-------------------------------------------------------------------------------- +-- HasChainRule instances. + +instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 0 ) where + konst w = D0 w + linearD f v = D0 ( f v ) + value ( D0 v ) = v + chain _ ( D0 gfx ) = D21 gfx origin origin + +instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 0 ) where + konst w = D0 w + linearD f v = D0 ( f v ) + value ( D0 v ) = v + chain _ ( D0 gfx ) = D31 gfx origin origin origin + +instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 1 ) where + + konst :: forall w. AbelianGroup w => w -> D2𝔸1 w + konst w = + let !o = fromInteger @w 0 + in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] ) + + value df = $$( monIndex [|| df ||] zeroMonomial ) + + linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 1 -> w ) -> 𝕀ℝ 1 -> D2𝔸1 w + linearD f v = + let !o = origin @( 𝕀 Double ) @( T w ) + in $$( monTabulate \ mon -> + if | isZeroMonomial mon + -> [|| f v ||] + | Just i <- isLinear mon + -> [|| f $$( tabulateQ \ j -> + if | j == i + -> [|| 1 ||] + | otherwise + -> [|| 0 ||] + ) ||] + | otherwise + -> [|| unT o ||] + ) + + chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D2𝔸1 ( 𝕀ℝ 1 ) -> D2𝔸1 w -> D2𝔸1 w + chain !df !dg = + let !o = origin @( 𝕀 Double ) @( T w ) + !p = (^+^) @( 𝕀 Double ) @( T w ) + !s = (^*) @( 𝕀 Double ) @( T w ) + in $$( chainRuleQ + [|| o ||] [|| p ||] [|| s ||] + [|| df ||] [|| dg ||] ) + +instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 1 ) where + + konst :: forall w. AbelianGroup w => w -> D3𝔸1 w + konst w = + let !o = fromInteger @w 0 + in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] ) + + value df = $$( monIndex [|| df ||] zeroMonomial ) + + linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 1 -> w ) -> 𝕀ℝ 1 -> D3𝔸1 w + linearD f v = + let !o = origin @( 𝕀 Double ) @( T w ) + in $$( monTabulate \ mon -> + if | isZeroMonomial mon + -> [|| f v ||] + | Just i <- isLinear mon + -> [|| f $$( tabulateQ \ j -> + if | j == i + -> [|| 1 ||] + | otherwise + -> [|| 0 ||] + ) ||] + | otherwise + -> [|| unT o ||] + ) + + chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D3𝔸1 ( 𝕀ℝ 1 ) -> D3𝔸1 w -> D3𝔸1 w + chain !df !dg = + let !o = origin @( 𝕀 Double ) @( T w ) + !p = (^+^) @( 𝕀 Double ) @( T w ) + !s = (^*) @( 𝕀 Double ) @( T w ) + in $$( chainRuleQ + [|| o ||] [|| p ||] [|| s ||] + [|| df ||] [|| dg ||] ) + +instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 2 ) where + + konst :: forall w. AbelianGroup w => w -> D2𝔸2 w + konst w = + let !o = fromInteger @w 0 + in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] ) + + value df = $$( monIndex [|| df ||] zeroMonomial ) + + linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 2 -> w ) -> 𝕀ℝ 2 -> D2𝔸2 w + linearD f v = + let !o = origin @( 𝕀 Double ) @( T w ) + in $$( monTabulate \ mon -> + if | isZeroMonomial mon + -> [|| f v ||] + | Just i <- isLinear mon + -> [|| f $$( tabulateQ \ j -> + if | j == i + -> [|| 1 ||] + | otherwise + -> [|| 0 ||] + ) ||] + | otherwise + -> [|| unT o ||] + ) + + chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D2𝔸1 ( 𝕀ℝ 2 ) -> D2𝔸2 w -> D2𝔸1 w + chain !df !dg = + let !o = origin @( 𝕀 Double ) @( T w ) + !p = (^+^) @( 𝕀 Double ) @( T w ) + !s = (^*) @( 𝕀 Double ) @( T w ) + in $$( chainRuleQ + [|| o ||] [|| p ||] [|| s ||] + [|| df ||] [|| dg ||] ) + +instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 2 ) where + + konst :: forall w. AbelianGroup w => w -> D3𝔸2 w + konst w = + let !o = fromInteger @w 0 + in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] ) + + value df = $$( monIndex [|| df ||] zeroMonomial ) + + linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 2 -> w ) -> 𝕀ℝ 2 -> D3𝔸2 w + linearD f v = + let !o = origin @( 𝕀 Double ) @( T w ) + in $$( monTabulate \ mon -> + if | isZeroMonomial mon + -> [|| f v ||] + | Just i <- isLinear mon + -> [|| f $$( tabulateQ \ j -> + if | j == i + -> [|| 1 ||] + | otherwise + -> [|| 0 ||] + ) ||] + | otherwise + -> [|| unT o ||] + ) + + chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D3𝔸1 ( 𝕀ℝ 2 ) -> D3𝔸2 w -> D3𝔸1 w + chain !df !dg = + let !o = origin @( 𝕀 Double ) @( T w ) + !p = (^+^) @( 𝕀 Double ) @( T w ) + !s = (^*) @( 𝕀 Double ) @( T w ) + in $$( chainRuleQ + [|| o ||] [|| p ||] [|| s ||] + [|| df ||] [|| dg ||] ) + +instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 3 ) where + + konst :: forall w. AbelianGroup w => w -> D2𝔸3 w + konst w = + let !o = fromInteger @w 0 + in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] ) + + value df = $$( monIndex [|| df ||] zeroMonomial ) + + linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 3 -> w ) -> 𝕀ℝ 3 -> D2𝔸3 w + linearD f v = + let !o = origin @( 𝕀 Double ) @( T w ) + in $$( monTabulate \ mon -> + if | isZeroMonomial mon + -> [|| f v ||] + | Just i <- isLinear mon + -> [|| f $$( tabulateQ \ j -> + if | j == i + -> [|| 1 ||] + | otherwise + -> [|| 0 ||] + ) ||] + | otherwise + -> [|| unT o ||] + ) + + chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D2𝔸1 ( 𝕀ℝ 3 ) -> D2𝔸3 w -> D2𝔸1 w + chain !df !dg = + let !o = origin @( 𝕀 Double ) @( T w ) + !p = (^+^) @( 𝕀 Double ) @( T w ) + !s = (^*) @( 𝕀 Double ) @( T w ) + in $$( chainRuleQ + [|| o ||] [|| p ||] [|| s ||] + [|| df ||] [|| dg ||] ) + +instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 3 ) where + + konst :: forall w. AbelianGroup w => w -> D3𝔸3 w + konst w = + let !o = fromInteger @w 0 + in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] ) + + value df = $$( monIndex [|| df ||] zeroMonomial ) + + linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 3 -> w ) -> 𝕀ℝ 3 -> D3𝔸3 w + linearD f v = + let !o = origin @( 𝕀 Double ) @( T w ) + in $$( monTabulate \ mon -> + if | isZeroMonomial mon + -> [|| f v ||] + | Just i <- isLinear mon + -> [|| f $$( tabulateQ \ j -> + if | j == i + -> [|| 1 ||] + | otherwise + -> [|| 0 ||] + ) ||] + | otherwise + -> [|| unT o ||] + ) + + chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D3𝔸1 ( 𝕀ℝ 3 ) -> D3𝔸3 w -> D3𝔸1 w + chain !df !dg = + let !o = origin @( 𝕀 Double ) @( T w ) + !p = (^+^) @( 𝕀 Double ) @( T w ) + !s = (^*) @( 𝕀 Double ) @( T w ) + in $$( chainRuleQ + [|| o ||] [|| p ||] [|| s ||] + [|| df ||] [|| dg ||] ) + +instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 4 ) where + + konst :: forall w. AbelianGroup w => w -> D2𝔸4 w + konst w = + let !o = fromInteger @w 0 + in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] ) + + value df = $$( monIndex [|| df ||] zeroMonomial ) + + linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 4 -> w ) -> 𝕀ℝ 4 -> D2𝔸4 w + linearD f v = + let !o = origin @( 𝕀 Double ) @( T w ) + in $$( monTabulate \ mon -> + if | isZeroMonomial mon + -> [|| f v ||] + | Just i <- isLinear mon + -> [|| f $$( tabulateQ \ j -> + if | j == i + -> [|| 1 ||] + | otherwise + -> [|| 0 ||] + ) ||] + | otherwise + -> [|| unT o ||] + ) + + chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D2𝔸1 ( 𝕀ℝ 4 ) -> D2𝔸4 w -> D2𝔸1 w + chain !df !dg = + let !o = origin @( 𝕀 Double ) @( T w ) + !p = (^+^) @( 𝕀 Double ) @( T w ) + !s = (^*) @( 𝕀 Double ) @( T w ) + in $$( chainRuleQ + [|| o ||] [|| p ||] [|| s ||] + [|| df ||] [|| dg ||] ) + +instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 4 ) where + + konst :: forall w. AbelianGroup w => w -> D3𝔸4 w + konst w = + let !o = fromInteger @w 0 + in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] ) + + value df = $$( monIndex [|| df ||] zeroMonomial ) + + linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 4 -> w ) -> 𝕀ℝ 4 -> D3𝔸4 w + linearD f v = + let !o = origin @( 𝕀 Double ) @( T w ) + in $$( monTabulate \ mon -> + if | isZeroMonomial mon + -> [|| f v ||] + | Just i <- isLinear mon + -> [|| f $$( tabulateQ \ j -> + if | j == i + -> [|| 1 ||] + | otherwise + -> [|| 0 ||] + ) ||] + | otherwise + -> [|| unT o ||] + ) + + chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D3𝔸1 ( 𝕀ℝ 4 ) -> D3𝔸4 w -> D3𝔸1 w + chain !df !dg = + let !o = origin @( 𝕀 Double ) @( T w ) + !p = (^+^) @( 𝕀 Double ) @( T w ) + !s = (^*) @( 𝕀 Double ) @( T w ) + in $$( chainRuleQ + [|| o ||] [|| p ||] [|| s ||] + [|| df ||] [|| dg ||] ) diff --git a/src/splines/Math/Interval/Internal.hs b/src/splines/Math/Interval/Internal.hs new file mode 100644 index 0000000..6b62fd3 --- /dev/null +++ b/src/splines/Math/Interval/Internal.hs @@ -0,0 +1,82 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE UndecidableInstances #-} + +{-# OPTIONS_GHC -Wno-orphans #-} + +module Math.Interval.Internal + ( type 𝕀 ) + where + +-- base +import Data.Monoid + ( Sum(..) ) + +-- rounded-hw +import Numeric.Rounded.Hardware + ( Rounded(..) ) +import Numeric.Rounded.Hardware.Interval.NonEmpty + ( Interval(..) ) +import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval + ( sup, inf, powInt ) + +-- MetaBrush +import Math.Linear + ( T(..) + , RepDim, RepresentableQ(..), Representable(..) + ) +import Math.Module + ( Module(..) ) +import Math.Ring + +-------------------------------------------------------------------------------- +-- Intervals. + +type 𝕀 = Interval + +deriving via ViaPrelude ( 𝕀 Double ) + instance AbelianGroup ( 𝕀 Double ) +deriving via ViaPrelude ( 𝕀 Double ) + instance AbelianGroup ( T ( 𝕀 Double ) ) + +instance Ring ( 𝕀 Double ) where + (*) = (Prelude.*) + x ^ n = Interval.powInt x ( Prelude.fromIntegral n ) + -- This is very important, as x^2 is not the same as x * x + -- in interval arithmetic. This ensures we don't + -- accidentally use (^) from Prelude. + +deriving via ViaPrelude ( 𝕀 Double ) + instance Field ( 𝕀 Double ) + +deriving via ViaPrelude ( 𝕀 Double ) + instance Transcendental ( 𝕀 Double ) + +type instance RepDim ( 𝕀 u ) = RepDim u +instance RepresentableQ r u => RepresentableQ ( 𝕀 r ) ( 𝕀 u ) where + tabulateQ f = + let !lo = tabulateQ @r @u ( \ i -> [|| getRounded $ Interval.inf $$( f i ) ||] ) + !hi = tabulateQ @r @u ( \ i -> [|| getRounded $ Interval.sup $$( f i ) ||] ) + in [|| I ( Rounded $$lo ) ( Rounded $$hi ) ||] + + indexQ d i = + [|| case $$d of + I ( Rounded lo ) ( Rounded hi ) -> + let !lo_i = $$( indexQ @r @u [|| lo ||] i ) + !hi_i = $$( indexQ @r @u [|| hi ||] i ) + in I ( Rounded lo_i ) ( Rounded hi_i ) + ||] +instance Representable r u => Representable ( 𝕀 r ) ( 𝕀 u ) where + tabulate f = + let !lo = tabulate @r @u ( \ i -> getRounded $ Interval.inf ( f i ) ) + !hi = tabulate @r @u ( \ i -> getRounded $ Interval.sup ( f i ) ) + in I ( Rounded lo ) ( Rounded hi ) + + index d i = + case d of + I ( Rounded lo ) ( Rounded hi ) -> + let !lo_i = index @r @u lo i + !hi_i = index @r @u hi i + in I ( Rounded lo_i ) ( Rounded hi_i ) + +deriving via Sum ( 𝕀 Double ) instance Module ( 𝕀 Double ) ( T ( 𝕀 Double ) ) diff --git a/src/splines/Math/Linear.hs b/src/splines/Math/Linear.hs index a78bd73..3b6d1bb 100644 --- a/src/splines/Math/Linear.hs +++ b/src/splines/Math/Linear.hs @@ -1,8 +1,8 @@ {-# LANGUAGE AllowAmbiguousTypes #-} {-# LANGUAGE PolyKinds #-} {-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE UndecidableInstances #-} -{-# LANGUAGE UnliftedNewtypes #-} module Math.Linear ( -- * Points and vectors @@ -11,22 +11,17 @@ module Math.Linear -- * Points and vectors (second version) , ℝ(..), T(.., V2, V3) , Fin(..), MFin(..) - , Dim, Representable(..), ApRep(..) - , injection, projection + , RepDim, RepresentableQ(..) + , Representable(..), injection, projection , Vec(..), (!), find, zipIndices - -- * Intervals - , 𝕀, 𝕀ℝ, singleton, nonDecreasing - , Extent(..), I ) where -- base import Data.Coerce - ( Coercible, coerce ) + ( coerce ) import Data.Kind - ( Type, Constraint ) -import Data.Monoid - ( Sum(..) ) + ( Type ) import GHC.Generics ( Generic, Generic1, Generically(..), Generically1(..) ) import GHC.TypeNats @@ -48,18 +43,12 @@ import Data.Group import Data.Group.Generics ( ) --- rounded-hw -import Numeric.Rounded.Hardware - ( Rounded(..) ) -import Numeric.Rounded.Hardware.Interval.NonEmpty - ( Interval(..) ) -import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval - ( sup, inf ) +-- MetaBrush +import Math.Linear.Internal -------------------------------------------------------------------------------- data Mat22 = Mat22 !Double !Double !Double !Double - data Segment p = Segment { segmentStart :: !p @@ -72,45 +61,11 @@ data Segment p = via Generically1 Segment deriving anyclass ( NFData, NFData1 ) -instance Show p => Show (Segment p) where - show (Segment s e) = show s ++ " -> " ++ show e +instance Show p => Show ( Segment p ) where + show ( Segment s e ) = show s ++ " -> " ++ show e -------------------------------------------------------------------------------- --- | Euclidean space \( \mathbb{R}^n \). -type ℝ :: Nat -> Type -data family ℝ n -data instance ℝ 0 = ℝ0 - deriving stock ( Show, Eq, Ord, Generic ) - deriving anyclass NFData -newtype instance ℝ 1 = ℝ1 { unℝ1 :: Double } - deriving stock ( Generic ) - deriving newtype ( Show, Eq, Ord, NFData ) -data instance ℝ 2 = ℝ2 {-# UNPACK #-} !Double {-# UNPACK #-} !Double - deriving stock Generic - deriving anyclass NFData - deriving stock ( Show, Eq, Ord ) -data instance ℝ 3 = ℝ3 {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-# UNPACK #-} !Double - deriving stock Generic - deriving anyclass NFData - deriving stock ( Show, Eq, Ord ) - -data instance ℝ 4 = ℝ4 {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-# UNPACK #-} !Double - deriving stock Generic - deriving anyclass NFData - deriving stock ( Show, Eq, Ord ) - -deriving via ApRep ( Sum Double ) ( ℝ n ) - instance Representable Double ( ℝ n ) => Semigroup ( T ( ℝ n ) ) -deriving via ApRep ( Sum Double ) ( ℝ n ) - instance Representable Double ( ℝ n ) => Monoid ( T ( ℝ n ) ) -deriving via ApRep ( Sum Double ) ( ℝ n ) - instance Representable Double ( ℝ n ) => Group ( T ( ℝ n ) ) -deriving via ApRep Double ( ℝ n ) - instance Representable Double ( ℝ n ) => Act ( T ( ℝ n ) ) ( ℝ n ) -deriving via ApRep Double ( ℝ n ) - instance Representable Double ( ℝ n ) => Torsor ( T ( ℝ n ) ) ( ℝ n ) - -- | Tangent space to Euclidean space. type T :: Type -> Type newtype T e = T { unT :: e } @@ -130,9 +85,8 @@ instance Act ( T Double ) Double where instance Torsor ( T Double ) Double where a --> b = T ( b - a ) - -instance {-# OVERLAPPING #-} Show ( ℝ n ) => Show ( T ( ℝ n ) ) where - show ( T p ) = "V" ++ drop 1 ( show p ) +--instance {-# OVERLAPPING #-} Show ( ℝ n ) => Show ( T ( ℝ n ) ) where +-- show ( T p ) = "V" ++ drop 1 ( show p ) deriving stock instance Show v => Show ( T v ) instance Applicative T where @@ -151,85 +105,15 @@ pattern V3 x y z = T ( ℝ3 x y z ) -------------------------------------------------------------------------------- --- | 1, ..., n -type Fin :: Nat -> Type -newtype Fin n = Fin Word - deriving stock Eq - --- | 0, ..., n -type MFin :: Nat -> Type -newtype MFin n = MFin Word - -type Dim :: k -> Nat -type family Dim v - -type instance Dim ( ℝ n ) = n - -type Representable :: Type -> Type -> Constraint -class Representable r v | v -> r where - tabulate :: ( Fin ( Dim v ) -> r ) -> v - index :: v -> Fin ( Dim v ) -> r - -instance Representable Double ( ℝ 0 ) where - {-# INLINE tabulate #-} - tabulate _ = ℝ0 - {-# INLINE index #-} - index _ _ = 0 - -instance Representable Double ( ℝ 1 ) where - {-# INLINE tabulate #-} - tabulate f = ℝ1 ( f ( Fin 1 ) ) - {-# INLINE index #-} - index ( ℝ1 x ) _ = x - -instance Representable Double ( ℝ 2 ) where - {-# INLINE tabulate #-} - tabulate f = ℝ2 ( f ( Fin 1 ) ) ( f ( Fin 2 ) ) - {-# INLINE index #-} - index ( ℝ2 x y ) = \ case - Fin 1 -> x - _ -> y - -instance Representable Double ( ℝ 3 ) where - {-# INLINE tabulate #-} - tabulate f = ℝ3 ( f ( Fin 1 ) ) ( f ( Fin 2 ) ) ( f ( Fin 3 ) ) - {-# INLINE index #-} - index ( ℝ3 x y z ) = \ case - Fin 1 -> x - Fin 2 -> y - _ -> z - -instance Representable Double ( ℝ 4 ) where - {-# INLINE tabulate #-} - tabulate f = ℝ4 ( f ( Fin 1 ) ) ( f ( Fin 2 ) ) ( f ( Fin 3 ) ) ( f ( Fin 4 ) ) - {-# INLINE index #-} - index ( ℝ4 x y z w ) = \ case - Fin 1 -> x - Fin 2 -> y - Fin 3 -> z - _ -> w - -{-# INLINE projection #-} -projection :: ( Representable r u, Representable r v ) - => ( Fin ( Dim v ) -> Fin ( Dim u ) ) - -> u -> v -projection f = \ u -> - tabulate \ i -> index u ( f i ) - -{-# INLINE injection #-} -injection :: ( Representable r u, Representable r v ) - => ( Fin ( Dim v ) -> MFin ( Dim u ) ) - -> u -> v -> v -injection f = \ u v -> - tabulate \ i -> case f i of - MFin 0 -> index v i - MFin j -> index u ( Fin j ) - infixr 5 `VS` type Vec :: Nat -> Type -> Type data Vec n a where VZ :: Vec 0 a VS :: a -> Vec n a -> Vec ( 1 + n ) a + -- can't be strict, otherwise we can't conveniently + -- unsafeCoerce from lists + +deriving stock instance Show a => Show ( Vec n a ) deriving stock instance Functor ( Vec n ) deriving stock instance Foldable ( Vec n ) @@ -241,106 +125,20 @@ VS a _ ! Fin 1 = a VS _ a ! Fin i = a ! Fin ( i - 1 ) _ ! _ = error "impossible: Fin 0 is uninhabited" -find :: forall l a. ( a -> a -> Bool ) -> Vec l a -> a -> MFin l -find eq v b = MFin ( go 1 v ) +find :: forall l a. ( a -> Bool ) -> Vec l a -> MFin l +find f v = MFin ( find_ 1 v ) where - go :: Word -> Vec n a -> Word - go j ( VS a as ) - | a `eq` b + find_ :: Word -> Vec n a -> Word + find_ j ( VS a as ) + | f a = j | otherwise - = go ( j + 1 ) as - go _ VZ = 0 + = find_ ( j + 1 ) as + find_ _ VZ = 0 zipIndices :: forall n a. Vec n a -> [ ( Fin n, a ) ] -zipIndices = go 0 +zipIndices = zipIndices_ 1 where - go :: forall i. Word -> Vec i a -> [ ( Fin n, a ) ] - go _ VZ = [] - go i (a `VS` as) = ( Fin i, a ) : go ( i + 1 ) as - --------------------------------------------------------------------------------- --- Instances in terms of representable. - --- | A newtype to hang off instances for representable functors. -newtype ApRep r u = ApRep { unApRep :: u } - -instance ( Representable r u, Coercible r m, Semigroup m ) => Semigroup ( ApRep m u ) where - ApRep a <> ApRep b = ApRep $ tabulate @r @u \ i -> - coerce $ (<>) @m ( coerce ( index @r @u a i ) ) ( coerce ( index @r @u b i ) ) - {-# INLINE (<>) #-} -instance ( Representable r u, Coercible r m, Monoid m ) => Monoid ( ApRep m u ) where - mempty = ApRep $ tabulate @r @u \ _ -> coerce $ mempty @m - {-# INLINE mempty #-} -instance ( Representable r u, Coercible r m, Group m ) => Group ( ApRep m u ) where - invert ( ApRep a ) = ApRep $ tabulate @r @u \ i -> - coerce $ invert @m $ coerce ( index @r @u a i ) - {-# INLINE invert #-} -instance ( Act ( T r ) r , Semigroup ( T u ), Representable r u ) => Act ( T u ) ( ApRep r u ) where - T g • ApRep a = ApRep $ tabulate @r @u \ i -> - coerce $ (•) @(T r) @r ( coerce $ index @r @u g i ) ( coerce ( index @r @u a i ) ) - {-# INLINE (•) #-} -instance ( Torsor ( T r ) r , Group ( T u ), Representable r u ) => Torsor ( T u ) ( ApRep r u ) where - ApRep a --> ApRep b = T $ tabulate @r @u \ i -> - coerce $ (-->) @(T r) @r ( coerce $ index @r @u a i ) ( coerce ( index @r @u b i ) ) - {-# INLINE (-->) #-} - --------------------------------------------------------------------------------- --- Intervals. - -type 𝕀 = Interval -type 𝕀ℝ n = 𝕀 ( ℝ n ) - --- Handling points and intervals uniformly. -data Extent = Point | Interval - -type I :: Extent -> Type -> Type -type family I i a where - I 'Point a = a - I 'Interval a = 𝕀 a - -singleton :: a -> 𝕀 a -singleton a = I ( Rounded a ) ( Rounded a ) - --- | Turn a non-decreasing function into a function on intervals. -nonDecreasing :: ( a -> b ) -> 𝕀 a -> 𝕀 b -nonDecreasing f ( I ( Rounded lo ) ( Rounded hi ) ) = - I ( Rounded $ f lo ) ( Rounded $ f hi ) - -type instance Dim ( 𝕀 u ) = Dim u - -instance Representable r u => Representable ( 𝕀 r ) ( 𝕀 u ) where - tabulate f = - let !lo = tabulate @r @u ( \ i -> getRounded $ Interval.inf ( f i ) ) - !hi = tabulate @r @u ( \ i -> getRounded $ Interval.sup ( f i ) ) - in I ( Rounded lo ) ( Rounded hi ) - {-# INLINE tabulate #-} - index ( I ( Rounded lo ) ( Rounded hi ) ) i = - let !lo_i = index @r @u lo i - !hi_i = index @r @u hi i - in I ( Rounded lo_i ) ( Rounded hi_i ) - {-# INLINE index #-} - -deriving via ApRep ( Sum ( 𝕀 Double ) ) ( 𝕀ℝ n ) - instance Representable ( 𝕀 Double ) ( 𝕀ℝ n ) => Semigroup ( T ( 𝕀ℝ n ) ) -deriving via ApRep ( Sum ( 𝕀 Double ) ) ( 𝕀ℝ n ) - instance Representable ( 𝕀 Double ) ( 𝕀ℝ n ) => Monoid ( T ( 𝕀ℝ n ) ) -deriving via ApRep ( Sum ( 𝕀 Double ) ) ( 𝕀ℝ n ) - instance Representable ( 𝕀 Double ) ( 𝕀ℝ n ) => Group ( T ( 𝕀ℝ n ) ) - -deriving via Sum ( 𝕀 Double ) - instance Semigroup ( T ( 𝕀 Double ) ) -deriving via Sum ( 𝕀 Double ) - instance Monoid ( T ( 𝕀 Double ) ) -deriving via Sum ( 𝕀 Double ) - instance Group ( T ( 𝕀 Double ) ) - -instance Act ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where - T g • a = coerce ( Sum g • a ) -instance Torsor ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where - a --> b = T $ getSum ( a --> b ) - -deriving via ApRep ( 𝕀 Double ) ( 𝕀ℝ n ) - instance Representable ( 𝕀 Double ) ( 𝕀ℝ n ) => Act ( T ( 𝕀ℝ n ) ) ( 𝕀ℝ n ) -deriving via ApRep ( 𝕀 Double ) ( 𝕀ℝ n ) - instance Representable ( 𝕀 Double ) ( 𝕀ℝ n ) => Torsor ( T ( 𝕀ℝ n ) ) ( 𝕀ℝ n ) + zipIndices_ :: forall i. Word -> Vec i a -> [ ( Fin n, a ) ] + zipIndices_ _ VZ = [] + zipIndices_ i (a `VS` as) = ( Fin i, a ) : zipIndices_ ( i + 1 ) as diff --git a/src/splines/Math/Linear/Dual.hs b/src/splines/Math/Linear/Dual.hs deleted file mode 100644 index 6afba83..0000000 --- a/src/splines/Math/Linear/Dual.hs +++ /dev/null @@ -1,863 +0,0 @@ -{-# LANGUAGE AllowAmbiguousTypes #-} -{-# LANGUAGE DuplicateRecordFields #-} -{-# LANGUAGE QuantifiedConstraints #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE UndecidableInstances #-} - -module Math.Linear.Dual where - --- base -import Control.Applicative - ( liftA2 ) -import Data.Coerce - ( Coercible, coerce ) -import Data.Foldable - ( toList ) -import Data.Kind - ( Type, Constraint ) -import GHC.Exts - ( Any, Word#, proxy#, plusWord# ) -import GHC.Generics - ( Generic, Generic1(..), Generically1(..) ) -import GHC.TypeNats -import Unsafe.Coerce - ( unsafeCoerce ) - --- acts -import Data.Act - ( Torsor ) - --- rounded-hw -import Numeric.Rounded.Hardware - ( Rounded(..) ) -import Numeric.Rounded.Hardware.Interval.NonEmpty - ( Interval(..) ) - --- MetaBrush -import Math.Module - ( Module(..) ) -import Math.Linear - --------------------------------------------------------------------------------- - --- | @C n u v@ is the space of @C^k@-differentiable maps from @u@ to @v@. -type C :: Nat -> Type -> Type -> Type -newtype C k u v = D { runD :: u -> D k u v } -deriving stock instance Functor ( D k u ) => Functor ( C k u ) - --- | \( C^2 \)-differentiable mappings. -type (~>) = C 2 --- | \( C^3 \)-differentiable mappings. -type (~~>) = C 3 - --- | @D k u v@ is the space of @k@-th order germs of functions from @u@ to @v@, --- represented by the algebra: --- --- \[ \mathbb{Z}[x_1, \ldots, x_n]/(x_1, \ldots, x_n)^(k+1) \otimes_\mathbb{Z} v \] --- --- when @u@ is of dimension @n@. -type D :: Nat -> Type -> Type -> Type -type family D k u - -type instance D n ( ℝ 0 ) = D𝔸0 - -type instance D 2 ( ℝ 1 ) = D2𝔸1 -type instance D 2 ( ℝ 2 ) = D2𝔸2 -type instance D 2 ( ℝ 3 ) = D2𝔸3 -type instance D 2 ( ℝ 4 ) = D2𝔸4 - -type instance D 3 ( ℝ 1 ) = D3𝔸1 -type instance D 3 ( ℝ 2 ) = D3𝔸2 -type instance D 3 ( ℝ 3 ) = D3𝔸3 -type instance D 3 ( ℝ 4 ) = D3𝔸3 - --- | \( \mathbb{Z} \otimes_\mathbb{Z} v \). -newtype D𝔸0 v = D0 { v :: v } - deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) - deriving Applicative - via Generically1 D𝔸0 -deriving via MonApRep Double D𝔸0 - instance HasChainRule 2 Double ( ℝ 0 ) D𝔸0 -deriving via MonApRep ( 𝕀 Double ) D𝔸0 - instance HasChainRule 2 ( 𝕀 Double ) ( 𝕀ℝ 0 ) D𝔸0 -deriving via MonApRep r D𝔸0 r - instance Num r => Num ( D𝔸0 r ) -deriving via ( ApAp r D𝔸0 v ) - instance ( Num r, Module r ( T v ) ) => Module ( D𝔸0 r ) ( D𝔸0 v ) - --- | \( \mathbb{Z}[x]/(x)^3 \otimes_\mathbb{Z} v \). -data D2𝔸1 v = - D21 { v :: !v - , dx :: !( T v ) - , dxdx :: !( T v ) - } - deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 ) - deriving Applicative - via Generically1 D2𝔸1 -deriving stock instance ( Show v, Show ( T v ) ) => Show ( D2𝔸1 v ) -deriving via MonApRep Double D2𝔸1 - instance HasChainRule 2 Double ( ℝ 1 ) D2𝔸1 -deriving via MonApRep ( 𝕀 Double ) D2𝔸1 - instance HasChainRule 2 ( 𝕀 Double ) ( 𝕀ℝ 1 ) D2𝔸1 -deriving via MonApRep r D2𝔸1 r - instance Num r => Num ( D2𝔸1 r ) -deriving via ( ApAp r D2𝔸1 v ) - instance ( Num r, Module r ( T v ) ) => Module ( D2𝔸1 r ) ( D2𝔸1 v ) - --- | \( \mathbb{Z}[x]/(x)^4 \otimes_\mathbb{Z} v \). -data D3𝔸1 v = - D31 { v :: !v - , dx :: !( T v ) - , dxdx :: !( T v ) - , dxdxdx :: !( T v ) - } - deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 ) - deriving Applicative - via Generically1 D3𝔸1 -deriving stock instance ( Show v, Show ( T v ) ) => Show ( D3𝔸1 v ) -deriving via MonApRep Double D3𝔸1 - instance HasChainRule 3 Double ( ℝ 1 ) D3𝔸1 -deriving via MonApRep ( 𝕀 Double ) D3𝔸1 - instance HasChainRule 3 ( 𝕀 Double ) ( 𝕀ℝ 1 ) D3𝔸1 -deriving via MonApRep r D3𝔸1 r - instance Num r => Num ( D3𝔸1 r ) -deriving via ( ApAp r D3𝔸1 v ) - instance ( Num r, Module r ( T v ) ) => Module ( D3𝔸1 r ) ( D3𝔸1 v ) - --- | \( \mathbb{Z}[x, y]/(x, y)^3 \otimes_\mathbb{Z} v \). -data D2𝔸2 v = - D22 { v :: !v - , dx, dy :: !( T v ) - , dxdx, dxdy, dydy :: !( T v ) - } - deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 ) - deriving Applicative - via Generically1 D2𝔸2 -deriving stock instance ( Show v, Show ( T v ) ) => Show ( D2𝔸2 v ) -deriving via MonApRep Double D2𝔸2 - instance HasChainRule 2 Double ( ℝ 2 ) D2𝔸2 -deriving via MonApRep ( 𝕀 Double ) D2𝔸2 - instance HasChainRule 2 ( 𝕀 Double ) ( 𝕀ℝ 2 ) D2𝔸2 -deriving via MonApRep r D2𝔸2 r - instance Num r => Num ( D2𝔸2 r ) -deriving via ( ApAp r D2𝔸2 v ) - instance ( Num r, Module r ( T v ) ) => Module ( D2𝔸2 r ) ( D2𝔸2 v ) - --- | \( \mathbb{Z}[x, y]/(x, y)^4 \otimes_\mathbb{Z} v \). -data D3𝔸2 v = - D32 { v :: !v - , dx, dy :: !( T v ) - , dxdx, dxdy, dydy :: !( T v ) - , dxdxdx, dxdxdy, dxdydy, dydydy :: !( T v ) - } - deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 ) - deriving Applicative - via Generically1 D3𝔸2 -deriving stock instance ( Show v, Show ( T v ) ) => Show ( D3𝔸2 v ) -deriving via MonApRep Double D3𝔸2 - instance HasChainRule 3 Double ( ℝ 2 ) D3𝔸2 -deriving via MonApRep ( 𝕀 Double ) D3𝔸2 - instance HasChainRule 3 ( 𝕀 Double ) ( 𝕀ℝ 2 ) D3𝔸2 -deriving via MonApRep r D3𝔸2 r - instance Num r => Num ( D3𝔸2 r ) -deriving via ( ApAp r D3𝔸2 v ) - instance ( Num r, Module r ( T v ) ) => Module ( D3𝔸2 r ) ( D3𝔸2 v ) - --- | \( \mathbb{Z}[x, y, z]/(x, y, z)^3 \otimes_\mathbb{Z} v \). -data D2𝔸3 v = - D23 { v :: !v - , dx, dy, dz :: !( T v ) - , dxdx, dxdy, dydy, dxdz, dydz, dzdz :: !( T v ) - } - deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 ) - deriving Applicative - via Generically1 D2𝔸3 -deriving stock instance ( Show v, Show ( T v ) ) => Show ( D2𝔸3 v ) -deriving via MonApRep Double D2𝔸3 - instance HasChainRule 2 Double ( ℝ 3 ) D2𝔸3 -deriving via MonApRep ( 𝕀 Double ) D2𝔸3 - instance HasChainRule 2 ( 𝕀 Double ) ( 𝕀ℝ 3 ) D2𝔸3 -deriving via MonApRep r D2𝔸3 r - instance Num r => Num ( D2𝔸3 r ) -deriving via ( ApAp r D2𝔸3 v ) - instance ( Num r, Module r ( T v ) ) => Module ( D2𝔸3 r ) ( D2𝔸3 v ) - --- | \( \mathbb{Z}[x, y, z]/(x, y, z)^4 \otimes_\mathbb{Z} v \). -data D3𝔸3 v = - D33 { v :: !v - , dx, dy, dz :: !( T v ) - , dxdx, dxdy, dydy, dxdz, dydz, dzdz :: !( T v ) - , dxdxdx, dxdxdy, dxdydy, dydydy, - dxdxdz, dxdydz, dxdzdz, dydydz, dydzdz, dzdzdz :: !( T v ) - } - deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 ) - deriving Applicative - via Generically1 D3𝔸3 -deriving stock instance ( Show v, Show ( T v ) ) => Show ( D3𝔸3 v ) -deriving via MonApRep Double D3𝔸3 - instance HasChainRule 3 Double ( ℝ 3 ) D3𝔸3 -deriving via MonApRep ( 𝕀 Double ) D3𝔸3 - instance HasChainRule 3 ( 𝕀 Double ) ( 𝕀ℝ 3 ) D3𝔸3 -deriving via MonApRep r D3𝔸3 r - instance Num r => Num ( D3𝔸3 r ) -deriving via ( ApAp r D3𝔸3 v ) - instance ( Num r, Module r ( T v ) ) => Module ( D3𝔸3 r ) ( D3𝔸3 v ) - --- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^3 \otimes_\mathbb{Z} v \). -data D2𝔸4 v = - D24 { v :: !v - , dx, dy, dz, dw :: !( T v ) - , dxdx, dxdy, dydy, dxdz, dydz, dzdz, dxdw, dydw, dzdw, dwdw :: !( T v ) - } - deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 ) - deriving Applicative - via Generically1 D2𝔸4 -deriving stock instance ( Show v, Show ( T v ) ) => Show ( D2𝔸4 v ) -deriving via MonApRep Double D2𝔸4 - instance HasChainRule 2 Double ( ℝ 4 ) D2𝔸4 -deriving via MonApRep ( 𝕀 Double ) D2𝔸4 - instance HasChainRule 2 ( 𝕀 Double ) ( 𝕀ℝ 4 ) D2𝔸4 -deriving via MonApRep r D2𝔸4 r - instance Num r => Num ( D2𝔸4 r ) -deriving via ( ApAp r D2𝔸4 v ) - instance ( Num r, Module r ( T v ) ) => Module ( D2𝔸4 r ) ( D2𝔸4 v ) - --- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^3 \otimes_\mathbb{Z} v \). -data D3𝔸4 v = - D34 { v :: !v - , dx, dy, dz, dw :: !( T v ) - , dxdx, dxdy, dydy, dxdz, dydz, dzdz, dxdw, dydw, dzdw, dwdw :: !( T v ) - , dxdxdx, dxdxdy, dxdydy, dydydy, - dxdxdz, dxdydz, dxdzdz, dydzdz, dydydz, dzdzdz, - dxdxdw, dxdydw, dydydw, dxdzdw, dydzdw, dzdzdw - , dxdwdw, dydwdw, dzdwdw, dwdwdw :: !( T v ) - } - deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 ) - deriving Applicative - via Generically1 D3𝔸4 -deriving stock instance ( Show v, Show ( T v ) ) => Show ( D3𝔸4 v ) -deriving via MonApRep Double D3𝔸4 - instance HasChainRule 3 Double ( ℝ 4 ) D3𝔸4 -deriving via MonApRep ( 𝕀 Double ) D3𝔸4 - instance HasChainRule 3 ( 𝕀 Double ) ( 𝕀ℝ 4 ) D3𝔸4 -deriving via MonApRep r D3𝔸4 r - instance Num r => Num ( D3𝔸4 r ) -deriving via ( ApAp r D3𝔸4 v ) - instance ( Num r, Module r ( T v ) ) => Module ( D3𝔸4 r ) ( D3𝔸4 v ) - --------------------------------------------------------------------------------- - -{- -type Differentiable :: Extent -> Type -> Constraint -class - ( Module ( I i Double ) ( T ( I i Double ) ) - , Torsor ( T ( I i u ) ) ( I i u ) - , Module ( D ( I i u ) ( I i Double ) ) ( D ( I i u ) ( I i ( ℝ 2 ) ) ) - , Representable ( I i Double ) ( I i u ) - , Floating ( D ( I i u ) ( I i Double ) ) - , Applicative ( D ( I i u ) ) - ) => Differentiable i u - -instance - ( Module ( I i Double ) ( T ( I i Double ) ) - , Torsor ( T ( I i u ) ) ( I i u ) - , Module ( D ( I i u ) ( I i Double ) ) ( D ( I i u ) ( I i ( ℝ 2 ) ) ) - , Representable ( I i Double ) ( I i u ) - , Floating ( D ( I i u ) ( I i Double ) ) - , Applicative ( D ( I i u ) ) - ) => Differentiable i u --} - --------------------------------------------------------------------------------- - --- | @Mon k n@ is the set of monomials in @n@ variables of degree less than or equal to @k@. -type Mon :: Nat -> Nat -> Type -newtype Mon k n = Mon { monDegs :: Vec n Word } -- sum <= k - -zeroMonomial :: forall n k. KnownNat n => Mon k n -zeroMonomial = Mon - $ unsafeCoerce @[ Word ] @( Vec n Word ) - $ replicate ( fromIntegral $ natVal' @n proxy# ) 0 - -isZeroMonomial :: Mon k n -> Bool -isZeroMonomial ( Mon ds ) = isZeroVector ds - - -isZeroVector :: forall i. Vec i Word -> Bool -isZeroVector VZ = True -isZeroVector ( 0 `VS` ds ) = isZeroVector ds -isZeroVector _ = False - -isLinear :: Mon k n -> Maybe Word -isLinear = go 1 . monDegs - where - go :: forall i. Word -> Vec i Word -> Maybe Word - go _ VZ = Nothing - go i ( 1 `VS` ds ) - | isZeroVector ds - = Just i - | otherwise - = Nothing - go i ( 0 `VS` ds ) - = go ( i + 1 ) ds - go _ _ = Nothing - -type Deg :: ( Type -> Type ) -> Nat -type Vars :: ( Type -> Type ) -> Nat -type family Deg f -type family Vars f - --- | @'MonomialBasis' f@ exhibits @f u@ as a free @r@-module with basis the --- monomials in @Vars u@ variables, of degree up to (and including) @Deg u@. --- --- This is used as an accessory class to derive all the algebra, but calls to --- @monTabulate@ and @monIndex@ should not remain in the optimised Core. -type MonomialBasis :: ( Type -> Type ) -> Constraint -class MonomialBasis f where - monTabulate :: ( Mon ( Deg f ) ( Vars f ) -> u ) -> f u - monIndex :: f u -> Mon ( Deg f ) ( Vars f ) -> u - --- | A newtype to hang off instances for functors representable --- by a monomial basis (see 'MonomialBasis'). -newtype MonApRep r dv w = MonApRep { unMonApRep :: dv w } - - -instance ( Num v, Applicative dr, MonomialBasis dr ) => Num ( MonApRep r dr v ) where - (+) = coerce $ liftA2 @dr ( (+) @v ) - (-) = coerce $ liftA2 @dr ( (-) @v ) - fromInteger i = MonApRep $ pure @dr $ fromInteger @v i - - MonApRep df1 * MonApRep df2 = - MonApRep $ - monTabulate \ mon -> - sum [ monIndex df1 m1 * monIndex df2 m2 | (m1, m2) <- split mon ] - - abs = error "no abs please" - signum = error "no signum please" - --- | Newtype for the module instance @Module r v => Module ( dr r ) ( dr v )@. -newtype ApAp r dr v = ApAp { unApAp :: dr v } - -instance ( Num ( dr r ), Module r ( T v ), Applicative dr ) => Module ( dr r ) ( ApAp r dr v ) where - ApAp u ^+^ ApAp v = ApAp $ liftA2 ( coerce $ (^+^) @r @( T v ) ) u v - ApAp u ^-^ ApAp v = ApAp $ liftA2 ( coerce $ (^-^) @r @( T v ) ) u v - origin = ApAp $ pure $ coerce $ origin @r @( T v ) - k *^ ApAp u = ApAp $ liftA2 ( coerce $ (*^) @r @( T v ) ) k u - --------------------------------------------------------------------------------- --- Chain rule and Faà di Bruno formula. - -class HasChainRule n r v dn_v where - chain :: Module r ( T w ) - => D n ( ℝ 1 ) v -> dn_v w -> D n ( ℝ 1 ) w - konst :: Module s ( T w ) => w -> dn_v w - value :: dn_v w -> w - linearD :: Module r ( T w ) => ( v -> w ) -> v -> dn_v w - -chainRule2 :: forall r v w - . ( HasChainRule 2 r v ( D 2 v ), Module r ( T w ) ) - => C 2 ( ℝ 1 ) v -> C 2 v w -> C 2 ( ℝ 1 ) w -chainRule2 ( D df ) ( D dg ) = - D \ x -> - case df x of - df_x@( D21 { v = f_x } ) -> - chain @2 @r @v @( D 2 v ) df_x ( dg f_x ) - - -uncurryD2 :: D 2 a ~ D 2 ( ℝ 1 ) - => D 2 ( ℝ 1 ) ( C 2 a b ) -> a -> D 2 ( ℝ 2 ) b -uncurryD2 ( D21 ( D b_t0 ) ( T ( D dbdt_t0 ) ) ( T ( D d2bdt2_t0 ) ) ) s0 = - let !( D21 b_t0s0 dbds_t0s0 d2bds2_t0s0 ) = b_t0 s0 - !( D21 dbdt_t0s0 d2bdtds_t0s0 _ ) = dbdt_t0 s0 - !( D21 d2bdt2_t0s0 _ _ ) = d2bdt2_t0 s0 - in D22 b_t0s0 ( T dbdt_t0s0 ) dbds_t0s0 ( T d2bdt2_t0s0 ) d2bdtds_t0s0 d2bds2_t0s0 - - - --- | Recover the underlying function, discarding all infinitesimal information. -fun :: forall n r v w. HasChainRule n r v ( D n v ) => C n v w -> ( v -> w ) -fun ( D df ) = value @n @r @v @( D n v ) . df -{-# INLINE fun #-} - --- | The differentiable germ of a coordinate variable. -var :: forall n r v - . ( Module r ( T r ), Representable r v, HasChainRule n r v ( D n v ) ) - => Fin ( Dim v ) -> C n v r -var i = D $ linearD @n @r @v @( D n v ) ( `index` i ) -{-# INLINE var #-} - - -instance forall r dn_v n v - . ( Representable r v - , MonomialBasis dn_v - , Deg dn_v ~ n - , MonomialBasis ( D n ( ℝ 1 ) ) - , Deg ( D n ( ℝ 1 ) ) ~ n - , Vars ( D n ( ℝ 1 ) ) ~ 1 - , KnownNat n, KnownNat ( Dim v ) - , Vars dn_v ~ Dim v - ) => HasChainRule n r v ( MonApRep r dn_v ) where - chain :: forall w - . Module r ( T w ) - => D n ( ℝ 1 ) v -> MonApRep r dn_v w -> D n ( ℝ 1 ) w - chain = coerce $ monChain @r @v @w @( D n ( ℝ 1 ) ) @dn_v - value ( MonApRep u ) = u `monIndex` zeroMonomial - {-# INLINE value #-} - konst :: forall s w . Module s ( T w ) => w -> MonApRep r dn_v w - konst k = MonApRep $ monTabulate \ mon -> - if | isZeroMonomial mon - -> k - | otherwise - -> unT $ origin @s @( T w ) - {-# INLINE konst #-} - linearD :: forall w. Module r ( T w ) => ( v -> w ) -> v -> MonApRep r dn_v w - linearD f v = MonApRep $ monTabulate \ mon -> - if | isZeroMonomial mon - -> f v - | Just i <- isLinear mon - -> f $ tabulate @r @v \ ( Fin j ) -> - if | j == i - -> 1 :: r - | otherwise - -> 0 :: r - | otherwise - -> unT $ origin @r @( T w ) - {-# INLINE linearD #-} - - -monChain :: forall r v w dr1 dv - . ( Module r ( T w ), Representable r v - , MonomialBasis dr1, Vars dr1 ~ 1 - , MonomialBasis dv , Vars dv ~ Dim v - , Deg dr1 ~ Deg dv - , KnownNat ( Dim v ) - , KnownNat ( Deg dr1 ) - ) - => dr1 v -> dv w -> dr1 w -monChain df dg = monTabulate @dr1 - $ doChain @r @w ( index @r @v ) ( monIndex @dr1 df ) ( monIndex @dv dg ) -{-# INLINE monChain #-} - --- | Compute the chain rule for the composition \( g(f_1(t), .., f_d(t)) \), --- for derivatives up to order @k@. -doChain :: forall r w k d v - . ( KnownNat k, KnownNat d, Module r ( T w ) ) - => ( v -> Fin d -> r ) - -> ( Mon k 1 -> v ) - -> ( Mon k d -> w ) - -> ( Mon k 1 -> w ) -doChain index_v df dg = \ ( Mon ( deg `VS` _ ) ) -> - unT $ moduleSum - [ T ( dg m_g ) ^* - ( sum [ fromIntegral ( faà k is ) - * product - [ product - [ ( df ( Mon $ f_deg `VS` VZ ) `index_v` v_index ) ^ pow - | ( f_deg, pow ) <- pows ] - | ( v_index, pows ) <- zipIndices is ] - | is <- multiSubsetsSum [ 1 .. k ] k ( monDegs m_g ) ] - ) - | m <- monomials @_ @d deg - , m_g <- subs m - ] - where - k = word @k -{-# INLINE doChain #-} - -moduleSum :: Module r u => [ u ] -> u -moduleSum = ( foldr (^+^) ) origin - --- | Faà di Bruno coefficient. -faà :: Word -> Vec n [ ( Word, Word ) ] -> Word -faà k multisubsets = - factorial k `div` - product [ factorial p * ( factorial i ) ^ p - | multisubset <- toList multisubsets - , ( i, p ) <- multisubset ] - -factorial :: Word -> Word -factorial i = product [ 1 .. i ] - --- | All monomials of degree **exactly** @k@ in @n@ variables, --- in lexicographic order. -monomials :: forall k n - . ( KnownNat n ) - => Word -- ^ degree - -> [ Mon k n ] -monomials k = unsafeCoerce @[ [ Word ] ] @[ Mon k n ] - $ monomials' k ( word @n ) - -monomials' :: Word -> Word -> [ [ Word ] ] -monomials' k n | k < 0 || n <= 0 = [] -monomials' 0 n = [ replicate ( fromIntegral n ) 0 ] -monomials' k n = addOne ( monomials' ( k - 1 ) n ) - ++ map ( 0 : ) ( monomials' k ( n - 1 ) ) - where - addOne :: [ [ Word ] ] -> [ [ Word ] ] - addOne [] = [] - addOne ( [] : _ ) = [] - addOne ( ( d : ds ) : dss ) = ( ( d + 1 ) : ds ) : addOne dss - --- | All monomials less than or equal to a given monomial. -subs :: Mon k n -> [ Mon k n ] -subs = unsafeCoerce subs' - -subs' :: [ Word ] -> [ [ Word ] ] -subs' [] = [ [] ] -subs' ( d : ds ) = - [ i : as - | i <- [ 0 .. d ] - , as <- subs' ds ] - -split :: Mon k n -> [ ( Mon k n, Mon k n ) ] -split = unsafeCoerce split' - -split' :: [ Word ] -> [ ( [ Word ], [ Word ] ) ] -split' [] = [ ( [], [] ) ] -split' ( d : ds ) = - [ ( i : as, d - i : bs ) - | i <- [ 0 .. d ] - , (as, bs) <- split' ds ] - --- | @multiSubsetsSum is s ns@ computes all collection of multisubsets of @is@, --- with sizes specifeid by @ns@, such that the total sum is @s@. -multiSubsetsSum :: forall n - . [ Word ] -- ^ set to pick from - -> Word -- ^ desired total sum - -> Vec n Word -- ^ sizes of each multisets - -> [ Vec n [ ( Word, Word ) ] ] -multiSubsetsSum is = go - where - go :: forall i. Word -> Vec i Word -> [ Vec i [ ( Word, Word ) ] ] - go 0 VZ = [ VZ ] - go _ VZ = [ ] - go s (n `VS` ns) = - [ multi `VS` rest - | s_i <- [ n * i_min .. s ] - , multi <- multiSubsetSum n s_i is - , rest <- go ( s - s_i ) ns ] - i_min = case is of - [] -> 0 - _ -> max 0 $ minimum is - --- | Computes the multisubsets of the given set which have the specified sum --- and number of elements. -multiSubsetSum :: Word -- ^ size of multisubset - -> Word -- ^ desired sum - -> [ Word ] -- ^ set to pick from - -> [ [ ( Word, Word ) ] ] -multiSubsetSum 0 0 _ = [ [] ] -multiSubsetSum 0 _ _ = [] -multiSubsetSum _ _ [] = [] -multiSubsetSum n s ( i : is ) = - [ if j == 0 then js else ( i, j ) : js - | j <- [ 0 .. n ] - , js <- multiSubsetSum ( n - j ) ( s - i * j ) is - ] - -word :: forall n. KnownNat n => Word -word = fromIntegral $ natVal' @n proxy# - - --------------------------------------------------------------------------------- --- MonomialBasis instances. - -type instance Deg D𝔸0 = 2 -type instance Vars D𝔸0 = 0 -instance MonomialBasis D𝔸0 where - monTabulate f = D0 { v } - where - v = f $ Mon VZ - {-# INLINE monTabulate #-} - monIndex ( D0 { v } ) = \ _ -> v - {-# INLINE monIndex #-} - -type instance Deg D2𝔸1 = 2 -type instance Vars D2𝔸1 = 1 -instance MonomialBasis D2𝔸1 where - monTabulate f = - D21 { v, dx, dxdx } - where - v = f $ Mon ( 0 `VS` VZ ) - dx = T $ f $ Mon ( 1 `VS` VZ ) - dxdx = T $ f $ Mon ( 2 `VS` VZ ) - {-# INLINE monTabulate #-} - - monIndex ( D21 { .. } ) = \ case - Mon ( 1 `VS` VZ ) -> unT dx - Mon ( 2 `VS` VZ ) -> unT dxdx - _ -> v - {-# INLINE monIndex #-} - -type instance Deg D3𝔸1 = 3 -type instance Vars D3𝔸1 = 1 -instance MonomialBasis D3𝔸1 where - monTabulate f = - D31 { v, dx, dxdx, dxdxdx } - where - v = f $ Mon ( 0 `VS` VZ ) - dx = T $ f $ Mon ( 1 `VS` VZ ) - dxdx = T $ f $ Mon ( 2 `VS` VZ ) - dxdxdx = T $ f $ Mon ( 3 `VS` VZ ) - {-# INLINE monTabulate #-} - - monIndex ( D31 { .. } ) = \ case - Mon ( 1 `VS` VZ ) -> unT dx - Mon ( 2 `VS` VZ ) -> unT dxdx - Mon ( 3 `VS` VZ ) -> unT dxdxdx - _ -> v - {-# INLINE monIndex #-} - - -type instance Deg D2𝔸2 = 2 -type instance Vars D2𝔸2 = 2 -instance MonomialBasis D2𝔸2 where - monTabulate f = D22 { .. } - where - v = f $ Mon ( 0 `VS` 0 `VS` VZ ) - dx = T $ f $ Mon ( 1 `VS` 0 `VS` VZ ) - dy = T $ f $ Mon ( 0 `VS` 1 `VS` VZ ) - dxdx = T $ f $ Mon ( 2 `VS` 0 `VS` VZ ) - dxdy = T $ f $ Mon ( 1 `VS` 1 `VS` VZ ) - dydy = T $ f $ Mon ( 0 `VS` 2 `VS` VZ ) - {-# INLINE monTabulate #-} - - monIndex ( D22 { .. } ) = \ case - Mon ( 1 `VS` 0 `VS` VZ ) -> unT dx - Mon ( 0 `VS` 1 `VS` VZ ) -> unT dy - Mon ( 2 `VS` 0 `VS` VZ ) -> unT dxdx - Mon ( 1 `VS` 1 `VS` VZ ) -> unT dxdy - Mon ( 0 `VS` 2 `VS` VZ ) -> unT dydy - _ -> v - {-# INLINE monIndex #-} - -type instance Deg D3𝔸2 = 3 -type instance Vars D3𝔸2 = 2 -instance MonomialBasis D3𝔸2 where - monTabulate f = D32 { .. } - where - v = f $ Mon ( 0 `VS` 0 `VS` VZ ) - dx = T $ f $ Mon ( 1 `VS` 0 `VS` VZ ) - dy = T $ f $ Mon ( 0 `VS` 1 `VS` VZ ) - dxdx = T $ f $ Mon ( 2 `VS` 0 `VS` VZ ) - dxdy = T $ f $ Mon ( 1 `VS` 1 `VS` VZ ) - dydy = T $ f $ Mon ( 0 `VS` 2 `VS` VZ ) - dxdxdx = T $ f $ Mon ( 3 `VS` 0 `VS` VZ ) - dxdxdy = T $ f $ Mon ( 2 `VS` 1 `VS` VZ ) - dxdydy = T $ f $ Mon ( 1 `VS` 2 `VS` VZ ) - dydydy = T $ f $ Mon ( 0 `VS` 3 `VS` VZ ) - {-# INLINE monTabulate #-} - - monIndex ( D32 { .. } ) = \ case - Mon ( 1 `VS` 0 `VS` VZ ) -> unT dx - Mon ( 0 `VS` 1 `VS` VZ ) -> unT dy - Mon ( 2 `VS` 0 `VS` VZ ) -> unT dxdx - Mon ( 1 `VS` 1 `VS` VZ ) -> unT dxdy - Mon ( 0 `VS` 2 `VS` VZ ) -> unT dydy - Mon ( 3 `VS` 0 `VS` VZ ) -> unT dxdxdx - Mon ( 2 `VS` 1 `VS` VZ ) -> unT dxdxdy - Mon ( 1 `VS` 2 `VS` VZ ) -> unT dxdydy - Mon ( 0 `VS` 3 `VS` VZ ) -> unT dydydy - _ -> v - {-# INLINE monIndex #-} - -type instance Deg D2𝔸3 = 2 -type instance Vars D2𝔸3 = 3 -instance MonomialBasis D2𝔸3 where - monTabulate f = D23 { .. } - where - v = f $ Mon ( 0 `VS` 0 `VS` 0 `VS` VZ ) - dx = T $ f $ Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) - dy = T $ f $ Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) - dz = T $ f $ Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) - dxdx = T $ f $ Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) - dxdy = T $ f $ Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) - dydy = T $ f $ Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) - dxdz = T $ f $ Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) - dydz = T $ f $ Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) - dzdz = T $ f $ Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) - {-# INLINE monTabulate #-} - - - monIndex ( D23 { .. } ) = \ case - Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) -> unT dx - Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) -> unT dy - Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) -> unT dz - Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) -> unT dxdx - Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) -> unT dxdy - Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) -> unT dydy - Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) -> unT dxdz - Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) -> unT dydz - Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) -> unT dzdz - _ -> v - {-# INLINE monIndex #-} - -type instance Deg D3𝔸3 = 3 -type instance Vars D3𝔸3 = 3 -instance MonomialBasis D3𝔸3 where - monTabulate f = D33 { .. } - where - v = f $ Mon ( 0 `VS` 0 `VS` 0 `VS` VZ ) - dx = T $ f $ Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) - dy = T $ f $ Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) - dz = T $ f $ Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) - dxdx = T $ f $ Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) - dxdy = T $ f $ Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) - dydy = T $ f $ Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) - dxdz = T $ f $ Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) - dydz = T $ f $ Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) - dzdz = T $ f $ Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) - dxdxdx = T $ f $ Mon ( 3 `VS` 0 `VS` 0 `VS` VZ ) - dxdxdy = T $ f $ Mon ( 2 `VS` 1 `VS` 0 `VS` VZ ) - dxdydy = T $ f $ Mon ( 1 `VS` 2 `VS` 0 `VS` VZ ) - dydydy = T $ f $ Mon ( 0 `VS` 3 `VS` 0 `VS` VZ ) - dxdxdz = T $ f $ Mon ( 2 `VS` 0 `VS` 1 `VS` VZ ) - dxdydz = T $ f $ Mon ( 1 `VS` 1 `VS` 1 `VS` VZ ) - dxdzdz = T $ f $ Mon ( 1 `VS` 0 `VS` 2 `VS` VZ ) - dydydz = T $ f $ Mon ( 0 `VS` 2 `VS` 1 `VS` VZ ) - dydzdz = T $ f $ Mon ( 0 `VS` 1 `VS` 2 `VS` VZ ) - dzdzdz = T $ f $ Mon ( 0 `VS` 0 `VS` 3 `VS` VZ ) - {-# INLINE monTabulate #-} - - monIndex ( D33 { .. } ) = \ case - Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) -> unT dx - Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) -> unT dy - Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) -> unT dz - Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) -> unT dxdx - Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) -> unT dxdy - Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) -> unT dydy - Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) -> unT dxdz - Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) -> unT dydz - Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) -> unT dzdz - Mon ( 3 `VS` 0 `VS` 0 `VS` VZ ) -> unT dxdxdx - Mon ( 2 `VS` 1 `VS` 0 `VS` VZ ) -> unT dxdxdy - Mon ( 1 `VS` 2 `VS` 0 `VS` VZ ) -> unT dxdydy - Mon ( 0 `VS` 3 `VS` 0 `VS` VZ ) -> unT dydydy - Mon ( 2 `VS` 0 `VS` 1 `VS` VZ ) -> unT dxdxdz - Mon ( 1 `VS` 1 `VS` 1 `VS` VZ ) -> unT dxdydz - Mon ( 1 `VS` 0 `VS` 2 `VS` VZ ) -> unT dxdzdz - Mon ( 0 `VS` 1 `VS` 2 `VS` VZ ) -> unT dydzdz - Mon ( 0 `VS` 0 `VS` 3 `VS` VZ ) -> unT dzdzdz - _ -> v - {-# INLINE monIndex #-} - -type instance Deg D2𝔸4 = 2 -type instance Vars D2𝔸4 = 4 -instance MonomialBasis D2𝔸4 where - monTabulate f = D24 { .. } - where - v = f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) - dx = T $ f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) - dy = T $ f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) - dz = T $ f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) - dw = T $ f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) - dxdx = T $ f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) - dxdy = T $ f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) - dydy = T $ f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) - dxdz = T $ f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) - dydz = T $ f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) - dzdz = T $ f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) - dxdw = T $ f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) - dydw = T $ f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) - dzdw = T $ f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) - dwdw = T $ f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) - {-# INLINE monTabulate #-} - - monIndex ( D24 { .. } ) = \ case - Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> unT dx - Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> unT dy - Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> unT dz - Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> unT dw - Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> unT dxdx - Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> unT dxdy - Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) -> unT dydy - Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> unT dxdz - Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) -> unT dydz - Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) -> unT dzdz - Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> unT dxdw - Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) -> unT dydw - Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) -> unT dzdw - Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) -> unT dwdw - _ -> v - {-# INLINE monIndex #-} - -type instance Deg D3𝔸4 = 3 -type instance Vars D3𝔸4 = 4 -instance MonomialBasis D3𝔸4 where - monTabulate f = D34 { .. } - where - v = f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) - dx = T $ f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) - dy = T $ f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) - dz = T $ f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) - dw = T $ f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) - dxdx = T $ f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) - dxdy = T $ f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) - dydy = T $ f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) - dxdz = T $ f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) - dydz = T $ f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) - dzdz = T $ f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) - dxdw = T $ f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) - dydw = T $ f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) - dzdw = T $ f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) - dwdw = T $ f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) - dxdxdx = T $ f $ Mon ( 3 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) - dxdxdy = T $ f $ Mon ( 2 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) - dxdydy = T $ f $ Mon ( 1 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) - dydydy = T $ f $ Mon ( 0 `VS` 3 `VS` 0 `VS` 0 `VS` VZ ) - dxdxdz = T $ f $ Mon ( 2 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) - dxdydz = T $ f $ Mon ( 1 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) - dxdzdz = T $ f $ Mon ( 1 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) - dydydz = T $ f $ Mon ( 0 `VS` 2 `VS` 1 `VS` 0 `VS` VZ ) - dydzdz = T $ f $ Mon ( 0 `VS` 1 `VS` 2 `VS` 0 `VS` VZ ) - dzdzdz = T $ f $ Mon ( 0 `VS` 0 `VS` 3 `VS` 0 `VS` VZ ) - dxdxdw = T $ f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) - dxdydw = T $ f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) - dydydw = T $ f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 1 `VS` VZ ) - dxdzdw = T $ f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) - dydzdw = T $ f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 1 `VS` VZ ) - dzdzdw = T $ f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 1 `VS` VZ ) - dxdwdw = T $ f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) - dydwdw = T $ f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 2 `VS` VZ ) - dzdwdw = T $ f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 2 `VS` VZ ) - dwdwdw = T $ f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 3 `VS` VZ ) - - {-# INLINE monTabulate #-} - - monIndex ( D34 { .. } ) = \ case - Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> unT dx - Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> unT dy - Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> unT dz - Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> unT dw - Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> unT dxdx - Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> unT dxdy - Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) -> unT dydy - Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> unT dxdz - Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) -> unT dydz - Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) -> unT dzdz - Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> unT dxdw - Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) -> unT dydw - Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) -> unT dzdw - Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) -> unT dwdw - Mon ( 3 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> unT dxdxdx - Mon ( 2 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> unT dxdxdy - Mon ( 1 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) -> unT dxdydy - Mon ( 0 `VS` 3 `VS` 0 `VS` 0 `VS` VZ ) -> unT dydydy - Mon ( 2 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> unT dxdxdz - Mon ( 1 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) -> unT dxdydz - Mon ( 1 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) -> unT dxdzdz - Mon ( 0 `VS` 2 `VS` 1 `VS` 0 `VS` VZ ) -> unT dydydz - Mon ( 0 `VS` 1 `VS` 2 `VS` 0 `VS` VZ ) -> unT dydzdz - Mon ( 0 `VS` 0 `VS` 3 `VS` 0 `VS` VZ ) -> unT dzdzdz - Mon ( 2 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> unT dxdxdw - Mon ( 1 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) -> unT dxdydw - Mon ( 0 `VS` 2 `VS` 0 `VS` 1 `VS` VZ ) -> unT dydydw - Mon ( 1 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) -> unT dxdzdw - Mon ( 0 `VS` 1 `VS` 1 `VS` 1 `VS` VZ ) -> unT dydzdw - Mon ( 0 `VS` 0 `VS` 2 `VS` 1 `VS` VZ ) -> unT dzdzdw - Mon ( 1 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) -> unT dxdwdw - Mon ( 0 `VS` 1 `VS` 0 `VS` 2 `VS` VZ ) -> unT dydwdw - Mon ( 0 `VS` 0 `VS` 1 `VS` 2 `VS` VZ ) -> unT dzdwdw - Mon ( 0 `VS` 0 `VS` 0 `VS` 3 `VS` VZ ) -> unT dwdwdw - _ -> v - {-# INLINE monIndex #-} diff --git a/src/splines/Math/Linear/Internal.hs b/src/splines/Math/Linear/Internal.hs new file mode 100644 index 0000000..fab9641 --- /dev/null +++ b/src/splines/Math/Linear/Internal.hs @@ -0,0 +1,154 @@ +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE TemplateHaskell #-} + + +module Math.Linear.Internal + ( ℝ(..) + , Fin(..), MFin(..) + , RepDim + , RepresentableQ(..) + , Representable(..), projection, injection + ) + where + +-- base +import Data.Kind + ( Type, Constraint ) +import GHC.Generics + ( Generic ) +import GHC.TypeNats + ( Nat ) + +-- template-haskell +import Language.Haskell.TH + ( CodeQ ) + +-- deepseq +import Control.DeepSeq + ( NFData) + +-------------------------------------------------------------------------------- + +-- | Euclidean space \( \mathbb{R}^n \). +type ℝ :: Nat -> Type +data family ℝ n + +data instance ℝ 0 = ℝ0 + deriving stock ( Show, Eq, Ord, Generic ) + deriving anyclass NFData +newtype instance ℝ 1 = ℝ1 { unℝ1 :: Double } + deriving stock ( Generic ) + deriving newtype ( Show, Eq, Ord, NFData ) +data instance ℝ 2 = ℝ2 { _ℝ2_x, _ℝ2_y :: {-# UNPACK #-} !Double } + deriving stock Generic + deriving anyclass NFData + deriving stock ( Show, Eq, Ord ) +data instance ℝ 3 = ℝ3 { _ℝ3_x, _ℝ3_y, _ℝ3_z :: {-# UNPACK #-} !Double } + deriving stock Generic + deriving anyclass NFData + deriving stock ( Show, Eq, Ord ) +data instance ℝ 4 = ℝ4 { _ℝ4_x, _ℝ4_y, _ℝ4_z, _ℝ4_w :: {-# UNPACK #-} !Double } + deriving stock Generic + deriving anyclass NFData + deriving stock ( Show, Eq, Ord ) + +-------------------------------------------------------------------------------- + +-- | 1, ..., n +type Fin :: Nat -> Type +newtype Fin n = Fin Word + deriving stock Eq + +-- | 0, ..., n +type MFin :: Nat -> Type +newtype MFin n = MFin Word + +type RepDim :: k -> Nat +type family RepDim v + +type RepresentableQ :: Type -> Type -> Constraint +class RepresentableQ r v | v -> r where + tabulateQ :: ( Fin ( RepDim v ) -> CodeQ r ) -> CodeQ v + indexQ :: CodeQ v -> Fin ( RepDim v ) -> CodeQ r + +type Representable :: Type -> Type -> Constraint +class Representable r v | v -> r where + tabulate :: ( Fin ( RepDim v ) -> r ) -> v + index :: v -> Fin ( RepDim v ) -> r + + +projection :: ( Representable r u, Representable r v ) + => ( Fin ( RepDim v ) -> Fin ( RepDim u ) ) + -> u -> v +projection f = \ u -> + tabulate \ i -> index u ( f i ) + +injection :: ( Representable r u, Representable r v ) + => ( Fin ( RepDim v ) -> MFin ( RepDim u ) ) + -> u -> v -> v +injection f = \ u v -> + tabulate \ i -> case f i of + MFin 0 -> index v i + MFin j -> index u ( Fin j ) + +-------------------------------------------------------------------------------- + +type instance RepDim ( ℝ n ) = n + +instance RepresentableQ Double ( ℝ 0 ) where + tabulateQ _ = [|| ℝ0 ||] + indexQ _ _ = [|| 0 ||] + +instance RepresentableQ Double ( ℝ 1 ) where + tabulateQ f = [|| ℝ1 $$( f ( Fin 1 ) ) ||] + indexQ p _ = [|| unℝ1 $$p ||] + +instance RepresentableQ Double ( ℝ 2 ) where + tabulateQ f = [|| ℝ2 $$( f ( Fin 1 ) ) $$( f ( Fin 2 ) ) ||] + indexQ p = \ case + Fin 1 -> [|| _ℝ2_x $$p ||] + _ -> [|| _ℝ2_y $$p ||] + +instance RepresentableQ Double ( ℝ 3 ) where + tabulateQ f = [|| ℝ3 $$( f ( Fin 1 ) ) $$( f ( Fin 2 ) ) $$( f ( Fin 3 ) ) ||] + indexQ p = \ case + Fin 1 -> [|| _ℝ3_x $$p ||] + Fin 2 -> [|| _ℝ3_y $$p ||] + _ -> [|| _ℝ3_z $$p ||] + +instance RepresentableQ Double ( ℝ 4 ) where + tabulateQ f = [|| ℝ4 $$( f ( Fin 1 ) ) $$( f ( Fin 2 ) ) $$( f ( Fin 3 ) ) $$( f ( Fin 4 ) ) ||] + indexQ p = \ case + Fin 1 -> [|| _ℝ4_x $$p ||] + Fin 2 -> [|| _ℝ4_y $$p ||] + Fin 3 -> [|| _ℝ4_z $$p ||] + _ -> [|| _ℝ4_w $$p ||] + +instance Representable Double ( ℝ 0 ) where + tabulate _ = ℝ0 + index _ _ = 0 + +instance Representable Double ( ℝ 1 ) where + tabulate f = ℝ1 ( f ( Fin 1 ) ) + index p _ = unℝ1 p + +instance Representable Double ( ℝ 2 ) where + tabulate f = ℝ2 ( f ( Fin 1 ) ) ( f ( Fin 2 ) ) + index p = \ case + Fin 1 -> _ℝ2_x p + _ -> _ℝ2_y p + +instance Representable Double ( ℝ 3 ) where + tabulate f = ℝ3 ( f ( Fin 1 ) ) ( f ( Fin 2 ) ) ( f ( Fin 3 ) ) + index p = \ case + Fin 1 -> _ℝ3_x p + Fin 2 -> _ℝ3_y p + _ -> _ℝ3_z p + +instance Representable Double ( ℝ 4 ) where + tabulate f = ℝ4 ( f ( Fin 1 ) ) ( f ( Fin 2 ) ) ( f ( Fin 3 ) ) ( f ( Fin 4 ) ) + index p = \ case + Fin 1 -> _ℝ4_x p + Fin 2 -> _ℝ4_y p + Fin 3 -> _ℝ4_z p + _ -> _ℝ4_w p diff --git a/src/splines/Math/Module.hs b/src/splines/Math/Module.hs index 7570d34..dea92d3 100644 --- a/src/splines/Math/Module.hs +++ b/src/splines/Math/Module.hs @@ -1,7 +1,10 @@ {-# LANGUAGE AllowAmbiguousTypes #-} {-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE UndecidableInstances #-} +{-# OPTIONS_GHC -Wno-orphans #-} + module Math.Module ( Module(..), lerp , Inner(..), Cross(..) @@ -9,6 +12,8 @@ module Math.Module , norm, squaredNorm, quadrance, distance , proj, projC, closestPointOnSegment , strictlyParallel, convexCombination + + , ViaModule(..) ) where @@ -18,7 +23,7 @@ import Control.Applicative import Control.Monad ( guard ) import Data.Coerce - ( Coercible, coerce ) + ( coerce ) import Data.Kind ( Type, Constraint ) import Data.Monoid @@ -32,36 +37,22 @@ import Data.Act ( (-->) ) ) --- rounded-hw -import Numeric.Rounded.Hardware - ( Rounded(..) ) -import Numeric.Rounded.Hardware.Interval.NonEmpty - ( Interval(..) ) +-- groups +import Data.Group + ( Group(..) ) -- MetaBrush import Math.Epsilon ( epsilon ) import Math.Linear +import Math.Module.Internal + ( Module(..), Inner(..) ) +import Math.Ring + ( Ring ) +import qualified Math.Ring as Ring -------------------------------------------------------------------------------- -infixl 6 ^+^, ^-^ -infix 9 ^*, *^ - -class Num r => Module r m | m -> r where - - {-# MINIMAL origin, (^+^), ( (^*) | (*^) ) #-} - - origin :: m - (^+^) :: m -> m -> m - (^-^) :: m -> m -> m - (*^) :: r -> m -> m - (^*) :: m -> r -> m - - (*^) = flip (^*) - (^*) = flip (*^) - m ^-^ n = m ^+^ -1 *^ n - instance ( Applicative f, Module r m ) => Module r ( Ap f m ) where origin = pure origin (^+^) = liftA2 (^+^) @@ -71,11 +62,6 @@ instance ( Applicative f, Module r m ) => Module r ( Ap f m ) where lerp :: forall v r p. ( Module r v, Torsor v p ) => r -> p -> p -> p lerp t p0 p1 = ( t *^ ( p0 --> p1 :: v ) ) • p0 -infixl 8 ^.^ - -class Module r m => Inner r m where - (^.^) :: m -> m -> r - class Module r m => Cross r m where cross :: m -> m -> r @@ -129,42 +115,14 @@ instance ( Torsor ( T u ) u, Module r ( T u ) ) => Interpolatable r u -------------------------------------------------------------------------------- -instance ( Representable r u, Coercible r m, Module r m ) => Module r ( ApRep m u ) where - origin = ApRep $ tabulate @r @u \ _ -> coerce $ origin @r @m - {-# INLINE origin #-} - ApRep a ^+^ ApRep b = ApRep $ tabulate @r @u \ i -> - coerce $ (^+^) @r @m ( coerce ( index @r @u a i ) ) ( coerce ( index @r @u b i ) ) - {-# INLINE (^+^) #-} - ApRep a ^-^ ApRep b = ApRep $ tabulate @r @u \ i -> - coerce $ (^-^) @r @m ( coerce ( index @r @u a i ) ) ( coerce ( index @r @u b i ) ) - {-# INLINE (^-^) #-} - k *^ ApRep a = ApRep $ tabulate @r @u \ i -> - coerce $ (*^) @r @m k ( coerce ( index @r @u a i ) ) - {-# INLINE (*^) #-} - -deriving via ( ApRep ( Sum Double ) ( ℝ n ) ) - instance Representable Double ( ℝ n ) => Module Double ( T ( ℝ n ) ) - -instance Num a => Module a ( Sum a ) where - - origin = Sum 0 - - (^+^) = (<>) - Sum x ^-^ Sum y = Sum ( x - y ) - - c *^ Sum x = Sum ( c * x ) - Sum x ^* c = Sum ( x * c ) - -instance Num a => Inner a ( Sum a ) where - Sum a ^.^ Sum b = a * b - -deriving via Sum Double instance Module Double ( T Double ) +instance Ring a => Inner a ( Sum a ) where + Sum a ^.^ Sum b = a Ring.* b instance Inner Double ( T ( ℝ 2 ) ) where - V2 x1 y1 ^.^ V2 x2 y2 = x1 * x2 + y1 * y2 + V2 x1 y1 ^.^ V2 x2 y2 = x1 Ring.* x2 + y1 Ring.* y2 instance Cross Double ( T ( ℝ 2 ) ) where - cross ( V2 x1 y1 ) ( V2 x2 y2 ) = x1 * y2 - x2 * y1 + cross ( V2 x1 y1 ) ( V2 x2 y2 ) = x1 Ring.* y2 Ring.- x2 Ring.* y1 -- | Compute whether two vectors point in the same direction, -- that is, whether each vector is a (strictly) positive multiple of the other. @@ -205,23 +163,62 @@ convexCombination v0 v1 u c10 = ( v0 ^-^ v1 ) `cross` u -------------------------------------------------------------------------------- --- Interval arithmetic using rounded-hw library. +-- Not sure how to set things up to automate the following... -deriving via Sum ( 𝕀 Double ) instance Module ( 𝕀 Double ) ( T ( 𝕀 Double ) ) +instance Module Double ( T ( ℝ 0 ) ) where + origin = T $$( tabulateQ \ _ -> [|| unT $ origin ||] ) + T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] ) + T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] ) + k *^ T a = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] ) -deriving via ApRep ( Sum ( 𝕀 Double ) ) ( 𝕀ℝ n ) - instance Representable ( 𝕀 Double ) ( 𝕀ℝ n ) => Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) +instance Module Double ( T ( ℝ 1 ) ) where + origin = T $$( tabulateQ \ _ -> [|| unT $ origin ||] ) + T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] ) + T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] ) + k *^ T a = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] ) -instance Inner ( 𝕀 Double ) ( T ( 𝕀ℝ 2 ) ) where - T ( I ( Rounded ( ℝ2 x1_lo y1_lo ) ) ( Rounded ( ℝ2 x1_hi y1_hi ) ) ) ^.^ - T ( I ( Rounded ( ℝ2 x2_lo y2_lo ) ) ( Rounded ( ℝ2 x2_hi y2_hi ) ) ) - = let !x1x2 = I ( Rounded x1_lo ) ( Rounded x1_hi ) * I ( Rounded x2_lo ) ( Rounded x2_hi ) - !y1y2 = I ( Rounded y1_lo ) ( Rounded y1_hi ) * I ( Rounded y2_lo ) ( Rounded y2_hi ) - in x1x2 + y1y2 +instance Module Double ( T ( ℝ 2 ) ) where + origin = T $$( tabulateQ \ _ -> [|| unT $ origin ||] ) + T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] ) + T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] ) + k *^ T a = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] ) -instance Cross ( 𝕀 Double ) ( T ( 𝕀ℝ 2 ) ) where - T ( I ( Rounded ( ℝ2 x1_lo y1_lo ) ) ( Rounded ( ℝ2 x1_hi y1_hi ) ) ) `cross` - T ( I ( Rounded ( ℝ2 x2_lo y2_lo ) ) ( Rounded ( ℝ2 x2_hi y2_hi ) ) ) - = let !x1y2 = I ( Rounded x1_lo ) ( Rounded x1_hi ) * I ( Rounded y2_lo ) ( Rounded y2_hi ) - !y2x1 = I ( Rounded x2_lo ) ( Rounded x2_hi ) * I ( Rounded y1_lo ) ( Rounded y1_hi ) - in x1y2 - y2x1 +instance Module Double ( T ( ℝ 3 ) ) where + origin = T $$( tabulateQ \ _ -> [|| unT $ origin ||] ) + T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] ) + T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] ) + k *^ T a = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] ) + +instance Module Double ( T ( ℝ 4 ) ) where + origin = T $$( tabulateQ \ _ -> [|| unT $ origin ||] ) + T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] ) + T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] ) + k *^ T a = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] ) + +deriving via ViaModule Double ( T ( ℝ n ) ) + instance Module Double ( T ( ℝ n ) ) => Semigroup ( T ( ℝ n ) ) +deriving via ViaModule Double ( T ( ℝ n ) ) + instance Module Double ( T ( ℝ n ) ) => Monoid ( T ( ℝ n ) ) +deriving via ViaModule Double ( T ( ℝ n ) ) + instance Module Double ( T ( ℝ n ) ) => Group ( T ( ℝ n ) ) + +deriving via ViaModule Double ( ℝ n ) + instance Module Double ( T ( ℝ n ) ) => Act ( T ( ℝ n ) ) ( ℝ n ) +deriving via ( ViaModule Double ( ℝ n ) ) + instance Module Double ( T ( ℝ n ) ) => Torsor ( T ( ℝ n ) ) ( ℝ n ) + +-------------------------------------------------------------------------------- + +newtype ViaModule r m = ViaModule { unViaModule :: m } + +instance Module r m => Semigroup ( ViaModule r m ) where + (<>) = coerce $ (^+^) @r @m +instance Module r m => Monoid ( ViaModule r m ) where + mempty = coerce $ origin @r @m +instance Module r m => Group ( ViaModule r m ) where + invert = coerce $ ( \ x -> (^-^) @r @m ( origin @r @m ) x ) + +instance ( Semigroup ( T m ), Module r ( T m ) ) => Act ( T m ) ( ViaModule r m ) where + g • ViaModule m = ViaModule ( unT $ g ^+^ T m ) +instance ( Group ( T m ), Module r ( T m ) ) => Torsor ( T m ) ( ViaModule r m ) where + ViaModule a --> ViaModule b = T ( unT $ T b ^-^ T a ) diff --git a/src/splines/Math/Module/Internal.hs b/src/splines/Math/Module/Internal.hs new file mode 100644 index 0000000..6c30f0a --- /dev/null +++ b/src/splines/Math/Module/Internal.hs @@ -0,0 +1,90 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE UndecidableInstances #-} + +module Math.Module.Internal where + +-- base +import Data.Coerce + ( coerce ) +import Data.Monoid + ( Sum(..) ) + +-- MetaBrush +import Math.Ring + ( Ring ) +import qualified Math.Ring as Ring +import Math.Linear + +-------------------------------------------------------------------------------- + +infixl 6 ^+^, ^-^ +infix 9 ^*, *^ + +class Ring r => Module r m | m -> r where + + {-# MINIMAL origin, (^+^), ( (^*) | (*^) ) #-} + + origin :: m + (^+^) :: m -> m -> m + (^-^) :: m -> m -> m + (*^) :: r -> m -> m + (^*) :: m -> r -> m + + (*^) = flip (^*) + (^*) = flip (*^) + m ^-^ n = m ^+^ Ring.fromInteger ( -1 :: Integer ) *^ n + +infixl 8 ^.^ + +class Module r m => Inner r m where + (^.^) :: m -> m -> r + +-------------------------------------------------------------------------------- + +instance Ring a => Module a ( Sum a ) where + + origin = Sum $ Ring.fromInteger ( 0 :: Integer ) + + (^+^) = coerce $ (Ring.+) @a + (^-^) = coerce $ (Ring.-) @a + + c *^ Sum x = Sum ( c Ring.* x ) + Sum x ^* c = Sum ( x Ring.* c ) + +deriving via Sum Double instance Module Double ( T Double ) + +{- +moduleViaRepresentable :: Q TH.Type -> Q TH.Type -> TH.DecsQ +moduleViaRepresentable r m = do + i <- TH.newName "i" + let tab b = ( ( TH.varE 'tabulateQ `TH.appTypeE` r ) `TH.appTypeE` m ) `TH.appE` TH.lamE [ TH.varP i ] b + idx a = ( ( TH.varE 'indexQ `TH.appTypeE` r ) `TH.appTypeE` m ) `TH.appE` TH.varE i + + boo1 <- tab ( TH.varE 'origin `TH.appTypeE` r `TH.appTypeE` m ) + + + let + + meths = + [ mkMeth 'origin $ pure boo1 -- $ tab ( TH.varE 'origin `TH.appTypeE` r `TH.appTypeE` m ) + --, mkMeth '(^+^) $ TH.varE '(^+^) `TH.appE` ( TH ) `TH.appE` () + ] + + i1 <- TH.instanceD ( pure [] ) ( TH.conT ''Module `TH.appT` r `TH.appT` m ) meths + return [i1] + + where + mkMeth nm body = TH.funD nm [ TH.clause [] ( TH.normalB body ) [] ] + + +moduleViaRepresentable :: ( ( i -> a ) -> TH.CodeQ b ) -> ( b -> i -> TH.CodeQ a ) -> Q TH.Type -> Q TH.Type -> TH.DecsQ +moduleViaRepresentable tab idx r m = do + [d| + instance Module $r $m where + origin = $( TH.unTypeCode $ tab \ _ -> origin ) + a ^+^ b = $( TH.unTypeCode $ tab \ i -> idx [|| a ||] i ^+^ idx [|| b ||] i ) + a ^-^ b = $( TH.unTypeCode $ tab \ i -> idx [|| a ||] i ^-^ idx [|| b ||] i ) + k *^ a = $( TH.unTypeCode $ tab \ i -> [|| {- k *^ -} $$( idx [|| a ||] i ) ||] ) + |] +-} diff --git a/src/splines/Math/Monomial.hs b/src/splines/Math/Monomial.hs new file mode 100644 index 0000000..d0eff63 --- /dev/null +++ b/src/splines/Math/Monomial.hs @@ -0,0 +1,174 @@ +{-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE QuantifiedConstraints #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE UndecidableInstances #-} + +{-# OPTIONS_GHC -ddump-splices -ddump-to-file #-} + +module Math.Monomial + ( Mon(..) + , MonomialBasis(..), Deg, Vars + , zeroMonomial, isZeroMonomial, isLinear + + , split, mons + + , faà, multiSubsetsSum, multiSubsetSum + + , prodRuleQ + + ) where + +-- base +import Data.Foldable + ( toList ) +import Data.Kind + ( Type, Constraint ) +import GHC.Exts + ( proxy# ) +import GHC.TypeNats + ( KnownNat, Nat, natVal' ) +import Unsafe.Coerce + ( unsafeCoerce ) + +-- template-haskell +import Language.Haskell.TH + ( CodeQ ) + +-- MetaBrush +import Math.Linear + ( Vec(..), Fin(..) ) +import TH.Utils + ( foldQ ) + +-------------------------------------------------------------------------------- + +-- | @Mon k n@ is the set of monomials in @n@ variables of degree less than or equal to @k@. +type Mon :: Nat -> Nat -> Type +newtype Mon k n = Mon { monDegs :: Vec n Word } -- sum <= k + deriving stock Show + +type Deg :: ( Type -> Type ) -> Nat +type Vars :: ( Type -> Type ) -> Nat +type family Deg f +type family Vars f + +zeroMonomial :: forall k n. KnownNat n => Mon k n +zeroMonomial = Mon $ unsafeCoerce @[ Word ] @( Vec n Word ) + $ replicate ( fromIntegral $ word @n ) 0 + +isZeroMonomial :: Mon k n -> Bool +isZeroMonomial ( Mon VZ ) = True +isZeroMonomial ( Mon ( i `VS` is ) ) + | i == 0 + = isZeroMonomial ( Mon is ) + | otherwise + = False + +isLinear :: Mon k n -> Maybe ( Fin n ) +isLinear = fmap Fin . go 1 + where + go :: forall k' n'. Word -> Mon k' n' -> Maybe Word + go _ ( Mon VZ ) + = Nothing + go j ( Mon ( i `VS` is ) ) + | i == 0 + = go ( j + 1 ) ( Mon is ) + | i == 1 && isZeroMonomial ( Mon is ) + = Just j + | otherwise + = Nothing + +-- | Comultiplication of monomials. +split :: Mon k n -> [ ( Mon k n, Mon k n ) ] +split ( Mon VZ ) = [ ( Mon VZ, Mon VZ ) ] +split ( Mon ( d `VS` ds ) ) = + [ ( Mon ( i `VS` as ), Mon ( ( d - i ) `VS` bs ) ) + | i <- [ 0 .. d ] + , ( Mon as, Mon bs ) <- ( split ( Mon ds ) ) + ] + +-- | All monomials of degree less than or equal to @k@ in @n@ variables, +-- in lexicographic order. +mons :: forall k n. ( KnownNat n ) => Word -> [ Mon k n ] +mons k = unsafeCoerce ( mons' k ( word @n ) ) + +mons' :: Word -> Word -> [ [ Word ] ] +mons' k _ | k < 0 = [] +mons' _ 0 = [ [] ] +mons' 0 n = [ replicate ( fromIntegral n ) 0 ] +mons' k n = [ i : is | i <- reverse [ 0 .. k ], is <- mons' ( k - i ) ( n - 1 ) ] + +word :: forall n. KnownNat n => Word +word = fromIntegral $ natVal' @n proxy# + +-- | Faà di Bruno coefficient (naive implementation). +faà :: Word -> Vec n [ ( Word, Word ) ] -> Word +faà k multisubsets = + factorial k `div` + product [ factorial p * ( factorial i ) ^ p + | multisubset <- toList multisubsets + , ( i, p ) <- multisubset ] + +-- | Computes the multisubsets of the given set which have the specified sum +-- and number of elements. +multiSubsetSum :: Word -- ^ size of multisubset + -> Word -- ^ desired sum + -> [ Word ] -- ^ set to pick from + -> [ [ ( Word, Word ) ] ] +multiSubsetSum 0 0 _ = [ [] ] +multiSubsetSum 0 _ _ = [] +multiSubsetSum _ _ [] = [] +multiSubsetSum n s ( i : is ) = + [ if j == 0 then js else ( i, j ) : js + | j <- [ 0 .. n ] + , js <- multiSubsetSum ( n - j ) ( s - i * j ) is + ] + +-- | @multiSubsetsSum is s ns@ computes all collection of multisubsets of @is@, +-- with sizes specified by @ns@, such that the total sum is @s@. +multiSubsetsSum :: forall n + . [ Word ] -- ^ set to pick from + -> Word -- ^ desired total sum + -> Vec n Word -- ^ sizes of each multisets + -> [ Vec n [ ( Word, Word ) ] ] +multiSubsetsSum is = goMSS + where + goMSS :: forall i. Word -> Vec i Word -> [ Vec i [ ( Word, Word ) ] ] + goMSS 0 VZ = [ VZ ] + goMSS _ VZ = [ ] + goMSS s (n `VS` ns) = + [ multi `VS` rest + | s_i <- [ n * i_min .. s ] + , multi <- multiSubsetSum n s_i is + , rest <- goMSS ( s - s_i ) ns ] + i_min = case is of + [] -> 0 + _ -> max 0 $ minimum is + +-- | The factorial function \( n! = n \cdot (n-1) \cdot `ldots` \cdot 2 `cdot` 1 \). +factorial :: Word -> Word +factorial i = product [ 1 .. i ] + +-------------------------------------------------------------------------------- + +-- | @'MonomialBasis' f@ exhibits @f u@ as a free @r@-module with basis the +-- monomials in @Vars u@ variables, of degree up to (and including) @Deg u@. +type MonomialBasis :: ( Type -> Type ) -> Constraint +class MonomialBasis f where + monTabulate :: ( ( Mon ( Deg f ) ( Vars f ) ) -> CodeQ u ) -> CodeQ ( f u ) + monIndex :: CodeQ ( f u ) -> Mon ( Deg f ) ( Vars f ) -> CodeQ u + +-------------------------------------------------------------------------------- + +prodRuleQ :: forall f r. MonomialBasis f + => CodeQ r -> CodeQ ( r -> r -> r ) -> CodeQ ( r -> r -> r ) + -- Ring r constraint (circumvent TH constraint problem) + -> CodeQ ( f r ) -> CodeQ ( f r ) -> CodeQ ( f r ) +prodRuleQ zero plus times df1 df2 = + monTabulate @f \ mon -> + [|| $$( foldQ plus zero + [ [|| $$times $$( monIndex @f df1 m1 ) $$( monIndex @f df2 m2 ) ||] + | (m1, m2) <- split mon + ] ) + ||] diff --git a/src/splines/Math/Ring.hs b/src/splines/Math/Ring.hs new file mode 100644 index 0000000..e1d088a --- /dev/null +++ b/src/splines/Math/Ring.hs @@ -0,0 +1,177 @@ +{-# LANGUAGE ScopedTypeVariables #-} + +module Math.Ring + ( AbelianGroup(..), Signed(..), Ring(..), Field(..), Transcendental(..) + + , ViaPrelude(..), ViaAbelianGroup(..) + + , ifThenElse + ) + where + +-- base +import Prelude ( Num, Fractional ) +import Prelude hiding ( Num(..), Fractional(..) ) +import qualified Prelude +import Control.Applicative + ( liftA2 ) +import Data.Coerce + ( coerce ) +import Data.Monoid + ( Ap(..) ) + +-- groups +import Data.Group + ( Group(invert) ) + +-------------------------------------------------------------------------------- + +infixl 6 +, - +infixl 7 *, / +infixr 8 ^ + +class AbelianGroup a where + {-# MINIMAL (+), fromInteger, ( (-) | negate ) #-} + + (+) :: a -> a -> a + (-) :: a -> a -> a + negate :: a -> a + fromInteger :: Integer -> a + + a - b = a + negate b + negate b = fromInteger 0 - b + +class AbelianGroup a => Signed a where + abs :: a -> a + signum :: a -> a + +class AbelianGroup a => Ring a where + (*) :: a -> a -> a + (^) :: a -> Word -> a + (^) = defaultPow + +class Ring a => Field a where + {-# MINIMAL fromRational, ( recip | (/) ) #-} + + fromRational :: Rational -> a + (/) :: a -> a -> a + recip :: a -> a + + recip x = fromInteger 1 / x + x / y = x * recip y + +class Field a => Transcendental a where + pi :: a + cos :: a -> a + sin :: a -> a +-- sqrt :: a -> a + +-------------------------------------------------------------------------------- + +ifThenElse :: Bool -> a -> a -> a +ifThenElse b x y = if b then x else y + +-------------------------------------------------------------------------------- +-- A default power function, taken from base (GHC.Real). +-- +-- DO NOT use with interval arithmetic! + +{-# RULES +"defaultPow 2" forall x. defaultPow x 2 = x*x +"defaultPow 3" forall x. defaultPow x 3 = x*x*x +"defaultPow 4" forall x. defaultPow x 4 = let u = x*x in u*u +"defaultPow 5" forall x. defaultPow x 5 = let u = x*x in u*u*x + #-} + +{-# INLINE [1] defaultPow #-} +defaultPow :: Ring a => a -> Word -> a +defaultPow x0 y0 + | y0 < 0 = errorWithoutStackTrace "Negative exponent" + | y0 == 0 = fromInteger 1 + | otherwise = powImpl x0 y0 + +powImpl :: Ring a => a -> Word -> a +-- powImpl : x0 ^ y0 = (x ^ y) +powImpl x y + | even y = powImpl (x * x) (y `quot` 2) + | y == 1 = x + | otherwise = powImplAcc (x * x) (y `quot` 2) x + +{-# INLINABLE powImplAcc #-} +powImplAcc :: Ring a => a -> Word -> a -> a +-- powImplAcc : x0 ^ y0 = (x ^ y) * z +powImplAcc x y z + | even y = powImplAcc (x * x) (y `quot` 2) z + | y == 1 = x * z + | otherwise = powImplAcc (x * x) (y `quot` 2) (x * z) + +-------------------------------------------------------------------------------- + +newtype ViaPrelude a = ViaPrelude { viaPrelude :: a } + +instance Num a => AbelianGroup ( ViaPrelude a ) where + (+) = coerce $ (Prelude.+) @a + (-) = coerce $ (Prelude.-) @a + fromInteger = coerce $ Prelude.fromInteger @a + +instance Num a => Signed ( ViaPrelude a ) where + abs = coerce $ Prelude.abs @a + signum = coerce $ Prelude.signum @a + +instance Num a => Ring ( ViaPrelude a ) where + (*) = coerce $ (Prelude.*) @a + (^) = coerce $ (Prelude.^) @a @Word + +instance Fractional a => Field ( ViaPrelude a ) where + fromRational = coerce $ Prelude.fromRational @a + (/) = coerce $ (Prelude./) @a + +instance Floating a => Transcendental ( ViaPrelude a ) where + pi = coerce $ Prelude.pi @a + sin = coerce $ Prelude.sin @a + cos = coerce $ Prelude.cos @a +-- sqrt = coerce $ Prelude.sqrt @a + +-------------------------------------------------------------------------------- + +newtype ViaAbelianGroup a = ViaAbelianGroup { viaAbelianGroup :: a } + +instance AbelianGroup a => Semigroup ( ViaAbelianGroup a ) where + (<>) = coerce $ (+) @a +instance AbelianGroup a => Monoid ( ViaAbelianGroup a ) where + mempty = ViaAbelianGroup $ fromInteger 0 +instance AbelianGroup a => Group ( ViaAbelianGroup a ) where + invert = coerce $ negate @a + +-------------------------------------------------------------------------------- + +deriving via ViaPrelude Integer instance AbelianGroup Integer +deriving via ViaPrelude Integer instance Ring Integer +deriving via ViaPrelude Integer instance Signed Integer + +deriving via ViaPrelude Int instance AbelianGroup Int +deriving via ViaPrelude Int instance Ring Int +deriving via ViaPrelude Int instance Signed Int + +deriving via ViaPrelude Word instance AbelianGroup Word +deriving via ViaPrelude Word instance Ring Word + +deriving via ViaPrelude Float instance AbelianGroup Float +deriving via ViaPrelude Float instance Ring Float +deriving via ViaPrelude Float instance Signed Float +deriving via ViaPrelude Float instance Field Float +deriving via ViaPrelude Float instance Transcendental Float + +deriving via ViaPrelude Double instance AbelianGroup Double +deriving via ViaPrelude Double instance Ring Double +deriving via ViaPrelude Double instance Signed Double +deriving via ViaPrelude Double instance Field Double +deriving via ViaPrelude Double instance Transcendental Double + +-------------------------------------------------------------------------------- + +instance ( AbelianGroup r, Applicative f ) => AbelianGroup ( Ap f r ) where + (+) = liftA2 $ (+) @r + (-) = liftA2 $ (-) @r + negate = fmap $ negate @r + fromInteger = pure . fromInteger @r diff --git a/src/splines/TH/Utils.hs b/src/splines/TH/Utils.hs new file mode 100644 index 0000000..23af12b --- /dev/null +++ b/src/splines/TH/Utils.hs @@ -0,0 +1,23 @@ +{-# LANGUAGE TemplateHaskell #-} + +module TH.Utils where + +-- template-haskell +import Language.Haskell.TH + ( CodeQ ) + +-- MetaBrush +import Math.Ring ( Ring ) +import qualified Math.Ring as Ring + +-------------------------------------------------------------------------------- + +foldQ :: CodeQ ( a -> a -> a ) -> CodeQ a -> [ CodeQ a ] -> CodeQ a +foldQ _ a0 [] = a0 +foldQ _ _ [a] = a +foldQ f a0 (a:as) = [|| $$f $$a $$( foldQ f a0 as ) ||] + +powQ :: Ring a => CodeQ a -> Word -> CodeQ a +powQ _ 0 = [|| Ring.fromInteger ( 1 :: Integer ) ||] +powQ x 1 = x +powQ x n = [|| $$x Ring.^ ( n :: Word ) ||]