more work into observability

This commit is contained in:
sheaf 2023-03-13 22:09:15 +01:00
parent 52420a1169
commit 7db3cbef33
7 changed files with 724 additions and 174 deletions

View file

@ -265,7 +265,7 @@ library metabrushes
, bytestring
>= 0.10.10.0 && < 0.12
executable test-cusp-isolation
executable cusps
import:
common
@ -273,11 +273,15 @@ executable test-cusp-isolation
hs-source-dirs:
src/cusps
main-is:
Main.hs
other-modules:
Math.Interval.Abstract
build-depends:
splines
main-is:
Main.hs
executable convert-metafont

View file

@ -1,27 +1,27 @@
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE PartialTypeSignatures #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE UndecidableInstances #-}
{-# OPTIONS_GHC -Wno-orphans #-}
module Main where
-- base
import Prelude hiding ( Num(..), Fractional(..), Floating(..), (^) )
import qualified Prelude ( Num(..), Fractional(..) )
import Data.Coerce
( Coercible )
import qualified Prelude
import GHC.Generics
( Generic )
-- acts
import Data.Act
( Torsor )
( Torsor(..) )
-- containers
import Data.Map
( Map )
import qualified Data.Map as Map
( (!) )
import qualified Data.Sequence as Seq
( index )
( fromList, index )
-- generic-lens
import Data.Generics.Product.Typed
@ -36,18 +36,22 @@ import Math.Bezier.Stroke.EnvelopeEquation
( StrokeDatum(..) )
import Math.Differentiable
( I )
import Math.Interval
import Math.Linear
( (..), T(..) )
( (..), T(..)
, Fin(..)
, Representable(..) )
import Math.Module
( Module, Cross )
( Module(..), Cross(..) )
import Math.Ring
( AbelianGroup(..), Ring(..), Field(..), Transcendental(..)
, ViaPrelude(..)
)
( AbelianGroup(..), Ring(..), Transcendental(..) )
import Math.Interval.Abstract
--------------------------------------------------------------------------------
main :: IO ()
main = print ellipseTest
data PointData params
= PointData
{ pointCoords :: !( 2 )
@ -55,33 +59,14 @@ data PointData params
}
deriving stock ( Show, Generic )
{-
newtype C k u v = D { runD :: u -> D k u v }
type instance D k ( n ) = Dk𝔸n
type instance D k ( 𝕀 v ) = D k v
e.g.
data D2𝔸2 v =
D22 { _D22_v :: !v
, _D22_dx, _D22_dy :: !( T v )
, _D22_dxdx, _D22_dxdy, _D22_dydy :: !( T v )
}
-}
outlineFunction
:: forall i brushParams
:: forall {t} (i :: t) brushParams
. ( Show brushParams
, D 1 ( I i ( 2 ) ) ~ D 1 ( 2 )
, D 2 ( I i ( 2 ) ) ~ D 2 ( 2 )
, D 3 ( I i ( 1 ) ) ~ D 3 ( 1 )
, D 3 ( I i ( 2 ) ) ~ D 3 ( 2 )
, Coercible ( I i ( 1 ) ) ( I i Double )
-- , Coercible ( I i ( 1 ) ) ( I i Double )
, HasType ( 2 ) ( PointData brushParams )
, Cross ( I i Double ) ( T ( I i ( 2 ) ) )
, Module ( I i Double ) ( T ( I i brushParams ) )
@ -92,118 +77,200 @@ outlineFunction
, Traversable ( D 3 brushParams )
, Traversable ( D 3 ( I i brushParams ) )
)
=> ( forall a. a -> I i a )
=> ( I i Double -> I i ( 1 ) )
-> ( I i ( 1 ) -> I i Double )
-> ( forall a. a -> I i a )
-> C 3 ( I i brushParams ) ( Spline Closed () ( I i ( 2 ) ) )
-- ^ brush shape
-> Int -- ^ brush segment index to consider
-> PointData brushParams -> Curve Open () ( PointData brushParams )
-- ^ stroke path
-> ( I i ( 1 ) -> I i ( 1 ) -> StrokeDatum 3 i )
outlineFunction single brush brush_index sp0 crv = strokeData
outlineFunction co1 co2 single brush brush_index sp0 crv = strokeData
where
path :: C 3 ( I i ( 1 ) ) ( I i ( 2 ) )
params :: C 3 ( I i ( 1 ) ) ( I i brushParams )
(path, params) =
pathAndUsedParams @3 @i single brushParams sp0 crv
pathAndUsedParams @3 @i co2 single brushParams sp0 crv
strokeData :: I i ( 1 ) -> I i ( 1 ) -> StrokeDatum 3 i
strokeData =
fmap ( `Seq.index` brush_index ) $
brushStrokeData @3 @brushParams @i path params brush
brushStrokeData @3 @brushParams @i co1 co2 path params brush
main :: IO ()
main = return ()
boo :: ( HasType ( 2 ) ( PointData brushParams )
, Show brushParams
, Module ( AI Double ) ( T ( AI brushParams ) )
, Torsor ( T ( AI brushParams ) ) ( AI brushParams )
, HasChainRule ( AI Double ) 3 ( AI brushParams )
, Traversable ( D 3 brushParams )
)
=> C 3 ( AI brushParams ) ( Spline Closed () ( AI ( 2 ) ) )
-> Int
-> PointData brushParams -> Curve Open () ( PointData brushParams )
-> StrokeDatum 3 AI
boo x i y z
= outlineFunction @AI
Coerce1 Coerce2 (Pt . Val)
x i y z
(Coerce1 $ Pt $ Var "t") (Coerce1 $ Pt $ Var "s")
ellipseTest :: StrokeDatum 3 AI
ellipseTest =
boo ( ellipseBrush @AI ( Pt . Val ) scaleI ) 1
( mkPoint ( 2 10 -20 ) 2 1 0 )
( LineTo { curveEnd = NextPoint $ mkPoint ( 2 10 10 ) 10 5 ( pi Prelude./ 4 )
, curveData = () } )
where
mkPoint :: 2 -> Double -> Double -> Double -> PointData ( 3 )
mkPoint pt a b phi = PointData pt ( 3 a b phi )
scaleI :: Double -> AI Double -> AI Double
scaleI x iv = Scale (Pt $ Val x) iv
--------------------------------------------------------------------------------
-- EDSL for inspection
type instance I AA a = AA a
instance HasChainRule ( AI Double ) 3 ( AI ( 1 ) ) where
data AAVal a where
Val :: a -> AAVal a
Var :: String -> AAVal Double
konst x = D31 x o o o
where o = T $ fromInteger (0 :: Integer)
value ( D31 { _D31_v = x } ) = x
linearD f x = D31 (f x) (T $ f $ Pt $ Val $ 1 1) o o
where o = origin
instance Show a => Show ( AAVal a ) where
show (Val v) = show v
show (Var v) = show v
chain :: forall w. Module ( AI Double ) ( T w ) => D3𝔸1 ( AI ( 1 ) ) -> D3𝔸1 w -> D3𝔸1 w
chain !df !dg =
let !o = origin @( AI Double ) @( T w )
!p = (^+^) @( AI Double ) @( T w )
!s = (^*) @( AI Double ) @( T w )
in $$( chainRule1NQ
[|| o ||] [|| p ||] [|| s ||]
[|| df ||] [|| dg ||] )
data AA a where
Pt :: AAVal a -> AA a
Ival :: AAVal a -> AAVal a -> AA a
(:+:) :: AA Double -> AA Double -> AA Double
(:-:) :: AA Double -> AA Double -> AA Double
Negate :: AA Double -> AA Double
(:*:) :: AA Double -> AA Double -> AA Double
(:^:) :: AA Double -> Word -> AA Double
Recip :: AA Double -> AA Double
Cos :: AA Double -> AA Double
Sin :: AA Double -> AA Double
Pi :: AA Double
instance HasChainRule ( AI Double ) 3 ( AI ( 2 ) ) where
instance Show a => Show (AA a) where
showsPrec prec = \case
Pt v -> showString "[" . showsPrec 0 v . showString "]"
Ival lo hi -> showString "[" . showsPrec 0 lo . showString "," . showsPrec 0 hi . showString "]"
iv1 :+: iv2 -> showParen (prec > 6) $ showsPrec 6 iv1 . showString " + " . showsPrec 7 iv2
iv1 :-: iv2 -> showParen (prec > 6) $ showsPrec 6 iv1 . showString " - " . showsPrec 7 iv2
iv1 :*: iv2 -> showParen (prec > 7) $ showsPrec 7 iv1 . showString " - " . showsPrec 8 iv2
iv :^: n -> showParen (prec > 8) $ showsPrec 9 iv . showString " ^ " . showsPrec 8 n
Negate iv -> showParen (prec > 10) $ showString "negate " . showsPrec 11 iv
Recip iv -> showParen (prec > 10) $ showString "recip " . showsPrec 11 iv
Cos iv -> showParen (prec > 10) $ showString "cos " . showsPrec 11 iv
Sin iv -> showParen (prec > 10) $ showString "sin " . showsPrec 11 iv
Pi -> showString "pi"
konst x = D32 x o o o o o o o o o
where o = T $ fromInteger (0 :: Integer)
value ( D32 { _D32_v = x } ) = x
linearD f x =
D32 (f x)
(T $ f $ Pt $ Val $ 2 1 0)
(T $ f $ Pt $ Val $ 2 0 1)
o o o o o o o
where o = origin
infixl 6 :+:
infixl 6 :-:
infixl 7 :*:
infixr 8 :^:
chain :: forall w. Module ( AI Double ) ( T w ) => D3𝔸1 ( AI ( 2 ) ) -> D3𝔸2 w -> D3𝔸1 w
chain !df !dg =
let !o = origin @( AI Double ) @( T w )
!p = (^+^) @( AI Double ) @( T w )
!s = (^*) @( AI Double ) @( T w )
in $$( chainRule1NQ
[|| o ||] [|| p ||] [|| s ||]
[|| df ||] [|| dg ||] )
val :: Map String Double -> AAVal a -> a
val _ ( Val a ) = a
val vars ( Var v ) = vars Map.! v
instance HasChainRule ( AI Double ) 3 ( AI ( 3 ) ) where
eval :: Map String Double -> AA a -> 𝕀 a
eval vars = \ case
Pt v -> let f = val vars v in 𝕀 f f
Ival lo hi -> 𝕀 (val vars lo) (val vars hi)
iv1 :+: iv2 -> eval vars iv1 + eval vars iv2
iv1 :-: iv2 -> eval vars iv1 - eval vars iv2
Negate iv -> negate $ eval vars iv
iv1 :*: iv2 -> eval vars iv1 * eval vars iv2
iv :^: n -> eval vars iv ^ n
Recip iv -> recip $ eval vars iv
Cos iv -> cos $ eval vars iv
Sin iv -> sin $ eval vars iv
Pi -> pi
konst x = D33 x o o o o o o o o o o o o o o o o o o o
where o = T $ fromInteger (0 :: Integer)
value ( D33 { _D33_v = x } ) = x
linearD f x =
D33 (f x)
(T $ f $ Pt $ Val $ 3 1 0 0)
(T $ f $ Pt $ Val $ 3 0 1 0)
(T $ f $ Pt $ Val $ 3 0 0 1)
o o o o o o o o o o o o o o o o
where o = origin
instance Prelude.Num ( AA Double ) where
(+) = (:+:)
(-) = (:-:)
negate = Negate
(*) = (:*:)
abs = error "No abs for abstract intervals"
signum = error "No signum for abstract intervals"
fromInteger = Pt . Val . fromInteger
chain :: forall w. Module ( AI Double ) ( T w ) => D3𝔸1 ( AI ( 3 ) ) -> D3𝔸3 w -> D3𝔸1 w
chain !df !dg =
let !o = origin @( AI Double ) @( T w )
!p = (^+^) @( AI Double ) @( T w )
!s = (^*) @( AI Double ) @( T w )
in $$( chainRule1NQ
[|| o ||] [|| p ||] [|| s ||]
[|| df ||] [|| dg ||] )
instance Prelude.Fractional ( AA Double ) where
recip = Recip
(/) = error "No division for abstract intervals"
fromRational = Pt . Val . fromRational
instance HasChainRule ( AI Double ) 3 ( AI ( 4 ) ) where
instance Ring ( AA Double ) where
(*) = (:*:)
(^) = (:^:)
konst x = D34 x o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o
where o = T $ fromInteger (0 :: Integer)
value ( D34 { _D34_v = x } ) = x
linearD f x =
D34 (f x)
(T $ f $ Pt $ Val $ 4 1 0 0 0)
(T $ f $ Pt $ Val $ 4 0 1 0 0)
(T $ f $ Pt $ Val $ 4 0 0 1 0)
(T $ f $ Pt $ Val $ 4 0 0 0 1)
o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o
where o = origin
deriving via ViaPrelude ( AA Double )
instance AbelianGroup ( AA Double )
deriving via ViaPrelude ( AA Double )
instance AbelianGroup ( T ( AA Double ) )
deriving via ViaPrelude ( AA Double )
instance Field ( AA Double )
chain :: forall w. Module ( AI Double ) ( T w ) => D3𝔸1 ( AI ( 4 ) ) -> D3𝔸4 w -> D3𝔸1 w
chain !df !dg =
let !o = origin @( AI Double ) @( T w )
!p = (^+^) @( AI Double ) @( T w )
!s = (^*) @( AI Double ) @( T w )
in $$( chainRule1NQ
[|| o ||] [|| p ||] [|| s ||]
[|| df ||] [|| dg ||] )
instance Transcendental ( AA Double ) where
pi = Pi
cos = Cos
sin = Sin
--------------------------------------------------------------------------------
-- Brushes (TODO copied from MetaBrush.Asset.Brushes)
κ :: Double
κ = 0.5519150244935105707435627227925
circleSpline :: forall {t} (i :: t) 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 ) ) ()
]
lastCrv =
Bezier3To ( p κ -1 ) ( p 1 -κ ) BackToStart ()
{-# INLINE circleSpline #-}
ellipseBrush :: forall {t} (i :: t) k irec
. ( irec ~ I i ( 3 )
, Module
( D k irec ( I i Double ) )
( D k irec ( I i ( 2 ) ) )
, Module ( I i Double ) ( T ( I i Double ) )
, HasChainRule ( I i Double ) k irec
, Representable ( I i Double ) irec
, Applicative ( D k irec )
, Transcendental ( D k irec ( I i Double ) )
)
=> ( forall a. a -> I i a )
-> ( Double -> I i Double -> I i Double )
-> C k irec ( Spline 'Closed () ( I i ( 2 ) ) )
ellipseBrush mkI scaleI =
D \ params ->
let a, b, phi :: D k irec ( 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 irec ( I i ( 2 ) )
mkPt x y
= let !a' = mul x a
!b' = mul y b
!c = cos phi
!s = sin phi
in ( a' * c - b' * s ) *^ e_x
^+^ ( a' * s + b' * c ) *^ e_y
in circleSpline @i @k @( 3 ) @( 2 ) mkPt
where
e_x, e_y :: D k irec ( I i ( 2 ) )
e_x = pure $ mkI $ 2 1 0
e_y = pure $ mkI $ 2 0 1
mul :: Double -> D k irec ( I i Double ) -> D k irec ( I i Double )
mul k = fmap ( scaleI k )
{-# INLINEABLE ellipseBrush #-}

View file

@ -0,0 +1,465 @@
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE UndecidableInstances #-}
module Math.Interval.Abstract where
-- base
import Prelude hiding
( Num(..), Fractional(..), Floating(..), (^)
, unzip )
import qualified Prelude
( Num(..), Fractional(..) )
import Data.Coerce
( Coercible, coerce )
import Data.Foldable
( toList )
import Data.List
( intersperse )
import GHC.Show
( showCommaSpace )
-- acts
import Data.Act
( Act(..), Torsor(..) )
-- containers
import Data.Map
( Map )
import qualified Data.Map as Map
( (!) )
-- groups
import Data.Group
( Group(..) )
-- MetaBrush
import Math.Algebra.Dual
import Math.Bezier.Stroke.EnvelopeEquation
( StrokeDatum(..) )
import Math.Differentiable
( I )
import Math.Interval
import Math.Linear
( (..), T(..)
, Vec(..), (!), unzip
, Fin(..), RepDim, Representable(..), RepresentableQ(..)
)
import Math.Module
( Module(..), Cross(..) )
import Math.Ring
( AbelianGroup(..), Ring(..), Field(..), Transcendental(..)
, ViaPrelude(..)
)
--------------------------------------------------------------------------------
-- EDSL for inspection
type instance D k (AI a) = D k (𝕀 a)
type instance I AI a = AI a
-- | A value: a constant, or a variable name.
data Value a where
Val :: a -> Value a
Var :: String -> Value Double
-- | A scalar or a vector.
data SV a where
Scalar :: SV Double
Vector :: IsVec v => SV v
instance Show a => Show ( Value a ) where
show (Val v) = show v
show (Var v) = v
-- | An abstract interval.
data AI a where
Pt :: Value a -> AI a
Ival :: Value a -> Value a -> AI a
(:+:) :: SV a -> AI a -> AI a -> AI a
(:-:) :: SV a -> AI a -> AI a -> AI a
Negate :: SV a -> AI a -> AI a
(:*:) :: AI Double -> AI Double -> AI Double
(:^:) :: AI Double -> Word -> AI Double
Recip :: AI Double -> AI Double
Scale :: IsVec v => AI Double -> AI v -> AI v
Cross :: AI ( 2 ) -> AI ( 2 ) -> AI Double
Cos :: AI Double -> AI Double
Sin :: AI Double -> AI Double
Pi :: AI Double
Coerce1 :: AI Double -> AI ( 1 )
Coerce2 :: AI ( 1 ) -> AI Double
Tabulate :: ( Eq ( n ), Show ( n ), Representable Double ( n ) ) => Vec n ( AI Double ) -> AI ( n )
Index :: ( Eq ( n ), Show ( n ), Representable Double ( n ) ) => AI ( n ) -> Fin n -> AI Double
instance ( Show a, Eq a ) => Show ( AI a ) where
showsPrec prec ai = case fst $ normalise ai of
Pt v -> showsPrec prec v
Ival lo hi -> showString "[" . showsPrec 0 lo . showCommaSpace . showsPrec 0 hi . showString "]"
(:+:) Scalar iv1 iv2 -> showParen (prec > 6) $ showsPrec 6 iv1 . showString " + " . showsPrec 7 iv2
(:+:) Vector iv1 iv2 -> showParen (prec > 6) $ showsPrec 6 iv1 . showString " + " . showsPrec 7 iv2
(:-:) Scalar iv1 iv2 -> showParen (prec > 6) $ showsPrec 6 iv1 . showString " - " . showsPrec 7 iv2
(:-:) Vector iv1 iv2 -> showParen (prec > 6) $ showsPrec 6 iv1 . showString " - " . showsPrec 7 iv2
iv1 :*: iv2 -> showParen (prec > 7) $ showsPrec 7 iv1 . showString " * " . showsPrec 8 iv2
Cross iv1 iv2 -> showParen (prec > 7) $ showsPrec 8 iv1 . showString " × " . showsPrec 8 iv2
iv :^: n -> showParen (prec > 8) $ showsPrec 9 iv . showString " ^ " . showsPrec 8 n
Scale k v -> showParen (prec > 9) $ showsPrec 10 k . showString " * " . showsPrec 10 v
Index iv ( Fin i ) -> showParen (prec > 9) $ showsPrec 9 iv . showString " ! " . showsPrec 10 i
Negate Scalar iv -> showParen (prec > 10) $ showString " -" . showsPrec 11 iv
Negate Vector iv -> showParen (prec > 10) $ showString " -" . showsPrec 11 iv
Recip iv -> showParen (prec > 10) $ showString " 1 / " . showsPrec 11 iv
Coerce1 iv -> showsPrec prec iv
Coerce2 iv -> showsPrec prec iv
Cos iv -> showParen (prec > 10) $ showString "cos " . showsPrec 11 iv
Sin iv -> showParen (prec > 10) $ showString "sin " . showsPrec 11 iv
Pi -> showString "pi"
Tabulate vs -> showString "(" . foldr (.) id ( intersperse showCommaSpace fields ) . showString ")"
where
fields :: [ ShowS ]
fields = [ showsPrec 0 v
| v <- toList vs
]
infixl 6 :+:
infixl 6 :-:
infixl 7 :*:
infixr 8 :^:
index_maybe :: ( Show ( n ), Representable Double ( n ) )
=> AI ( n ) -> Fin n -> Maybe ( AI Double )
index_maybe expr i = case expr of
Pt ( Val v )
-> Just $ Pt $ Val $ v `index` i
Ival ( Val lo ) ( Val hi )
-> Just $ Ival ( Val $ lo `index` i ) ( Val $ hi `index` i )
(:+:) Vector iv1 iv2
| Just v1 <- iv1 `index_maybe` i
, Just v2 <- iv2 `index_maybe` i
-> Just $ v1 + v2
(:-:) Vector iv1 iv2
| Just v1 <- iv1 `index_maybe` i
, Just v2 <- iv2 `index_maybe` i
-> Just $ v1 - v2
Negate Vector iv
| Just v <- iv `index_maybe` i
-> Just $ negate v
Scale k iv
| Just v <- iv `index_maybe` i
-> Just $ k * v
Tabulate v
-> Just $ v ! i
Coerce1 v
-> Just $ v
_ -> Nothing
ival :: 𝕀 a -> AI a
ival ( 𝕀 lo hi ) = Ival ( Val lo ) ( Val hi )
normalise :: Eq a => AI a -> ( AI a, Maybe ( 𝕀 a ) )
normalise expr = case expr of
Pt ( Val v ) -> ( expr, Just $ 𝕀 v v )
Pt {} -> ( expr, Nothing )
Ival ( Val v1 ) ( Val v2 )
| v1 == v2
-> ( Pt ( Val v1 ), Just $ 𝕀 v1 v2 )
| otherwise
-> ( expr, Just $ 𝕀 v1 v2 )
Ival {} -> ( expr, Nothing )
(:+:) sv iv1 iv2
| Scalar <- sv
, Just ( 𝕀 0 0 ) <- mv1
-> ( iv2', mv2 )
| Vector <- sv
, Just o <- mv1
, T o == origin
-> ( iv2', mv2 )
| Scalar <- sv
, Just ( 𝕀 0 0 ) <- mv2
-> ( iv1', mv1 )
| Vector <- sv
, Just o <- mv2
, T o == origin
-> ( iv1', mv1 )
| Just v1 <- mv1
, Just v2 <- mv2
-> let r = case sv of { Scalar -> v1 + v2; Vector -> unT $ T v1 ^+^ T v2 }
in ( ival r, Just $ r )
| otherwise
-> ( expr', Nothing )
where
expr' = (:+:) sv iv1' iv2'
(iv1', mv1) = normalise iv1
(iv2', mv2) = normalise iv2
(:-:) sv iv1 iv2
| Scalar <- sv
, Just ( 𝕀 0 0 ) <- mv1
-> normalise $ Negate sv iv2
| Vector <- sv
, Just o <- mv1
, T o == origin
-> normalise $ Negate sv iv2
| Scalar <- sv
, Just ( 𝕀 0 0 ) <- mv2
-> ( iv1', mv1 )
| Vector <- sv
, Just o <- mv2
, T o == origin
-> ( iv1', mv1 )
| Just v1 <- mv1
, Just v2 <- mv2
-> let r = case sv of { Scalar -> v1 - v2; Vector -> unT $ T v1 ^-^ T v2 }
in ( ival r, Just $ r )
| otherwise
-> ( expr', Nothing )
where
expr' = (:-:) sv iv1' iv2'
(iv1', mv1) = normalise iv1
(iv2', mv2) = normalise iv2
Negate sv iv
| Just v <- mv
-> let r = case sv of { Scalar -> unT $ negate $ T v; Vector -> unT $ origin ^-^ T v }
in ( ival r , Just r )
| otherwise
-> ( expr', Nothing )
where
expr' = Negate sv iv'
(iv', mv) = normalise iv
iv1 :*: iv2
| Just ( 𝕀 0 0 ) <- mv1
-> zero
| Just ( 𝕀 0 0 ) <- mv2
-> zero
| Just v1 <- mv1
, Just v2 <- mv2
-> let r = v1 + v2
in ( ival r, Just $ r )
| otherwise
-> ( expr', Nothing )
where
zero = ( ival $ 𝕀 0 0, Just $ 𝕀 0 0 )
expr' = (:*:) iv1' iv2'
(iv1', mv1) = normalise iv1
(iv2', mv2) = normalise iv2
iv :^: n
| Just v <- mv
-> let r = v ^ n
in ( ival r , Just r )
| otherwise
-> ( expr', Nothing )
where
expr' = iv' :^: n
(iv', mv) = normalise iv
Recip iv
| Just v <- mv
-> let r = recip v
in ( ival r , Just r )
| otherwise
-> ( expr', Nothing )
where
expr' = Recip iv'
(iv', mv) = normalise iv
Scale k iv
| Just ( 𝕀 0 0 ) <- mk
-> ( ival $ unT $ origin, Just $ unT origin )
| Just o <- mv
, T o == origin
-> ( ival $ unT $ origin, Just $ unT origin )
| Just ( 𝕀 1 1 ) <- mk
-> ( iv', mv )
| Just val_k <- mk
, Just val_v <- mv
-> let r = unT $ val_k *^ T val_v
in ( ival r, Just $ r )
| otherwise
-> ( expr', Nothing )
where
expr' = Scale k' iv'
(k' , mk) = normalise k
(iv', mv) = normalise iv
Cross iv1 iv2
| Just v1 <- mv1
, Just v2 <- mv2
-> let r = T v1 × T v2
in ( ival r, Just $ r )
| otherwise
-> ( expr', Nothing )
where
expr' = Cross iv1' iv2'
(iv1', mv1) = normalise iv1
(iv2', mv2) = normalise iv2
Cos iv
| Just v <- mv
-> let r = cos v
in ( ival r, Just r )
| otherwise
-> ( expr', Nothing )
where
expr' = Cos iv'
(iv', mv) = normalise iv
Sin iv
| Just v <- mv
-> let r = sin v
in ( ival r, Just r )
| otherwise
-> ( expr', Nothing )
where
expr' = Sin iv'
(iv', mv) = normalise iv
Pi -> ( Pi, Just pi )
Coerce1 iv
| Just v <- mv
-> let r = coerce v
in ( ival r, Just r )
| otherwise
-> ( expr', Nothing )
where
expr' = Coerce1 iv'
(iv', mv) = normalise iv
Coerce2 iv
| Just v <- mv
-> let r = coerce v
in ( ival r, Just r )
| otherwise
-> ( expr', Nothing )
where
expr' = Coerce2 iv'
(iv', mv) = normalise iv
Tabulate @n v
| Just xs <- sequence mxs
-> let r = tabulate ( xs ! )
in ( ival r, Just r )
| otherwise
-> ( expr', Nothing )
where
expr' = Tabulate v'
v' :: Vec n ( AI Double )
mxs :: Vec n ( Maybe ( 𝕀 Double ) )
(v', mxs) = unzip $ fmap normalise v
Index iv i
| Just r <- index_maybe iv' i
-> normalise r
| Just v <- mv
-> let r = v `index` i
in ( ival r, Just r )
| otherwise
-> ( expr', Nothing )
where
expr' = Index iv' i
(iv', mv) = normalise iv
val :: Map String Double -> Value a -> a
val _ ( Val a ) = a
val vars ( Var v ) = vars Map.! v
eval :: Map String Double -> AI a -> 𝕀 a
eval vars = \ case
Pt v -> let f = val vars v in 𝕀 f f
Ival lo hi -> 𝕀 (val vars lo) (val vars hi)
(:+:) Scalar iv1 iv2 -> eval vars iv1 + eval vars iv2
(:+:) Vector iv1 iv2 -> unT $ (^+^) ( T ( eval vars iv1 ) ) ( T ( eval vars iv2 ) )
(:-:) Scalar iv1 iv2 -> eval vars iv1 - eval vars iv2
(:-:) Vector iv1 iv2 -> unT $ (^-^) ( T ( eval vars iv1 ) ) ( T ( eval vars iv2 ) )
Negate Scalar iv -> unT $ negate $ T $ eval vars iv
Negate Vector iv -> unT $ invert $ T $ eval vars iv
iv1 :*: iv2 -> eval vars iv1 * eval vars iv2
iv :^: n -> eval vars iv ^ n
Recip iv -> recip $ eval vars iv
Scale k v -> unT $ eval vars k *^ T ( eval vars v )
Cross iv1 iv2 -> T ( eval vars iv1 ) × T ( eval vars iv2 )
Cos iv -> cos $ eval vars iv
Sin iv -> sin $ eval vars iv
Pi -> pi
Coerce1 a -> coerce $ eval vars a
Coerce2 a -> coerce $ eval vars a
Tabulate v -> let !w = fmap ( eval vars ) v
in tabulate ( w ! )
Index v i -> let !w = eval vars v
in index w i
instance Prelude.Num ( AI Double ) where
(+) = (:+:) Scalar
(-) = (:-:) Scalar
negate = Negate Scalar
(*) = (:*:)
abs = error "No abs for abstract intervals"
signum = error "No signum for abstract intervals"
fromInteger = Pt . Val . fromInteger
instance Prelude.Fractional ( AI Double ) where
recip = Recip
(/) = error "No division for abstract intervals"
fromRational = Pt . Val . fromRational
instance Ring ( AI Double ) where
(*) = (:*:)
(^) = (:^:)
deriving via ViaPrelude ( AI Double )
instance AbelianGroup ( AI Double )
deriving via ViaPrelude ( AI Double )
instance AbelianGroup ( T ( AI Double ) )
deriving via ViaPrelude ( AI Double )
instance Field ( AI Double )
instance Transcendental ( AI Double ) where
pi = Pi
cos = Cos
sin = Sin
class ( Eq v, Show v, Module ( 𝕀 Double ) ( T ( 𝕀 v ) ), Group ( T ( 𝕀 v ) ) )
=> IsVec v
instance ( Eq v, Show v, Module ( 𝕀 Double ) ( T ( 𝕀 v ) ), Group ( T ( 𝕀 v ) ) )
=> IsVec v
instance IsVec v => Module ( AI Double ) ( T ( AI v ) ) where
origin = T $ Pt $ Val $ inf $ unT ( origin :: T ( 𝕀 v ) )
(^+^) = coerce $ (:+:) Vector
(^-^) = coerce $ (:-:) Vector
(*^) = coerce Scale
instance Cross ( AI Double ) ( T ( AI ( 2 ) ) ) where
(×) = coerce Cross
instance IsVec v => Semigroup ( T ( AI v ) ) where
(<>) = (^+^)
instance IsVec v => Monoid ( T ( AI v ) ) where
mempty = origin
instance IsVec v => Group ( T ( AI v ) ) where
invert = coerce $ Negate Vector
instance IsVec v => Act ( T ( AI v ) ) ( AI v ) where
u v = unT $ u ^+^ T v
instance IsVec v => Torsor ( T ( AI v ) ) ( AI v ) where
u --> v = T v ^-^ T u
deriving stock instance Show ( StrokeDatum 3 AI )
type instance RepDim ( AI u ) = RepDim u
instance Representable ( AI Double ) ( AI ( 1 ) ) where
tabulate f = Coerce1 ( f ( Fin 1 ) )
index x _ = Coerce2 x
instance Representable ( AI Double ) ( AI ( 2 ) ) where
tabulate f = Tabulate $ ( f ( Fin 1 ) ) `VS` ( f ( Fin 2 ) ) `VS` VZ
index x i = Index x i
instance Representable ( AI Double ) ( AI ( 3 ) ) where
tabulate f = Tabulate $ ( f ( Fin 1 ) ) `VS` ( f ( Fin 2 ) ) `VS` ( f ( Fin 3 ) ) `VS` VZ
index x i = Index x i
instance Representable ( AI Double ) ( AI ( 4 ) ) where
tabulate f = Tabulate $ ( f ( Fin 1 ) ) `VS` ( f ( Fin 2 ) ) `VS` ( f ( Fin 3 ) ) `VS` ( f ( Fin 4 ) ) `VS` VZ
index x i = Index x i
instance RepresentableQ ( AI Double ) ( AI ( 1 ) ) where
tabulateQ f = [|| Coerce1 $$( f ( Fin 1 ) ) ||]
indexQ x _ = [|| Coerce2 $$x ||]
instance RepresentableQ ( AI Double ) ( AI ( 2 ) ) where
tabulateQ f = [|| Tabulate $ $$( f ( Fin 1 ) ) `VS` $$( f ( Fin 2 ) ) `VS` VZ ||]
indexQ x ( Fin i ) = [|| Index $$x ( Fin i ) ||]
instance RepresentableQ ( AI Double ) ( AI ( 3 ) ) where
tabulateQ f = [|| Tabulate $ $$( f ( Fin 1 ) ) `VS` $$( f ( Fin 2 ) ) `VS` $$( f ( Fin 3 ) ) `VS` VZ ||]
indexQ x ( Fin i ) = [|| Index $$x ( Fin i ) ||]
instance RepresentableQ ( AI Double ) ( AI ( 4 ) ) where
tabulateQ f = [|| Tabulate $ $$( f ( Fin 1 ) ) `VS` $$( f ( Fin 2 ) ) `VS` $$( f ( Fin 3 ) ) `VS` $$( f ( Fin 4 ) ) `VS` VZ ||]
indexQ x ( Fin i ) = [|| Index $$x ( Fin i ) ||]

View file

@ -15,6 +15,8 @@ module Math.Algebra.Dual
, uncurryD2, uncurryD3
, linear, fun, var
, chainRuleN1Q, chainRule1NQ
, D𝔸0(..)
, D1𝔸1(..), D2𝔸1(..), D3𝔸1(..)
, D1𝔸2(..), D2𝔸2(..), D3𝔸2(..)

View file

@ -534,12 +534,13 @@ outlineFunction ptParams toBrushParams brushFromParams = \ sp0 crv ->
usedParams :: C 2 ( 1 ) usedParams
path :: C 2 ( 1 ) ( 2 )
( path, usedParams )
= pathAndUsedParams @2 @() id ptParams sp0 crv
= pathAndUsedParams @2 @() coerce id ptParams sp0 crv
curves :: 1 -- t
-> Seq ( 1 {- s -} -> StrokeDatum 2 () )
curves =
brushStrokeData @2 @brushParams
coerce coerce
path
( chainRule @Double @2
usedParams
@ -550,6 +551,7 @@ outlineFunction ptParams toBrushParams brushFromParams = \ sp0 crv ->
curvesI :: 𝕀 1 -- t
-> Seq ( 𝕀 1 {- s -} -> StrokeDatum 3 𝕀 )
curvesI = brushStrokeData @3 @brushParams
coerce coerce
pathI
( chainRule @( 𝕀 Double ) @3
usedParamsI
@ -559,7 +561,7 @@ outlineFunction ptParams toBrushParams brushFromParams = \ sp0 crv ->
usedParamsI :: C 3 ( 𝕀 1 ) ( 𝕀 usedParams )
pathI :: C 3 ( 𝕀 1 ) ( 𝕀 2 )
( pathI, usedParamsI ) = pathAndUsedParams @3 @𝕀 singleton ptParams sp0 crv
( pathI, usedParamsI ) = pathAndUsedParams @3 @𝕀 coerce singleton ptParams sp0 crv
fwdBwd :: OutlineFn
fwdBwd t
@ -601,31 +603,31 @@ pathAndUsedParams :: forall k i arr crvData ptData usedParams
, CurveOrder k
, arr ~ C k
, D k ( I i ( 1 ) ) ~ D k ( 1 )
, Coercible ( I i ( 1 ) ) ( I i Double )
, Module ( I i Double ) ( T ( I i ( 2 ) ) )
, Torsor ( T ( I i ( 2 ) ) ) ( I i ( 2 ) )
, Module ( I i Double ) ( T ( I i usedParams ) )
, Torsor ( T ( I i usedParams ) ) ( I i usedParams )
)
=> ( forall a. a -> I i a )
=> ( I i ( 1 ) -> I i Double )
-> ( forall a. a -> I i a )
-> ( ptData -> usedParams )
-> ptData
-> Curve Open crvData ptData
-> ( I i ( 1 ) `arr` I i ( 2 ), I i ( 1 ) `arr` I i usedParams )
pathAndUsedParams toI ptParams sp0 crv =
pathAndUsedParams co toI ptParams sp0 crv =
case crv of
LineTo { curveEnd = NextPoint sp1 }
| let seg = Segment sp0 sp1
-> ( line @k @i ( fmap ( toI . coords ) seg )
, line @k @i ( fmap ( toI . ptParams ) seg ) )
-> ( line @k @i co ( fmap ( toI . coords ) seg )
, line @k @i co ( fmap ( toI . ptParams ) seg ) )
Bezier2To { controlPoint = sp1, curveEnd = NextPoint sp2 }
| let bez2 = Quadratic.Bezier sp0 sp1 sp2
-> ( bezier2 @k @i ( fmap ( toI . coords ) bez2 )
, bezier2 @k @i ( fmap ( toI . ptParams ) bez2 ) )
-> ( bezier2 @k @i co ( fmap ( toI . coords ) bez2 )
, bezier2 @k @i co ( fmap ( toI . ptParams ) bez2 ) )
Bezier3To { controlPoint1 = sp1, controlPoint2 = sp2, curveEnd = NextPoint sp3 }
| let bez3 = Cubic.Bezier sp0 sp1 sp2 sp3
-> ( bezier3 @k @i ( fmap ( toI . coords ) bez3 )
, bezier3 @k @i ( fmap ( toI . ptParams ) bez3 ) )
-> ( bezier3 @k @i co ( fmap ( toI . coords ) bez3 )
, bezier3 @k @i co ( fmap ( toI . ptParams ) bez3 ) )
-----------------------------------
-- Various utility functions
@ -924,10 +926,11 @@ splineCurveFns :: forall k i
. ( CurveOrder k
, D k ( I i ( 1 ) ) ~ D k ( 1 )
, Module ( I i Double ) ( T ( I i ( 2 ) ) )
, Torsor ( T ( I i ( 2 ) ) ) ( I i ( 2 ) )
, Coercible ( I i ( 1 ) ) ( I i Double ) )
=> Spline Closed () ( I i ( 2 ) ) -> Seq ( C k ( I i ( 1 ) ) ( I i ( 2 ) ) )
splineCurveFns spls
, Torsor ( T ( I i ( 2 ) ) ) ( I i ( 2 ) ) )
=> ( I i ( 1 ) -> I i Double )
-> Spline Closed () ( I i ( 2 ) )
-> Seq ( C k ( I i ( 1 ) ) ( I i ( 2 ) ) )
splineCurveFns co spls
= runIdentity
. bifoldSpline
( \ pt crv -> Identity . Seq.singleton $ curveFn pt crv )
@ -940,11 +943,11 @@ splineCurveFns spls
-> ( C k ( I i ( 1 ) ) ( I i ( 2 ) ) )
curveFn p0 = \case
LineTo { curveEnd = NextPoint p1 }
-> line @k @i $ Segment p0 p1
-> line @k @i co $ Segment p0 p1
Bezier2To { controlPoint = p1, curveEnd = NextPoint p2 }
-> bezier2 @k @i $ Quadratic.Bezier p0 p1 p2
-> bezier2 @k @i co $ Quadratic.Bezier p0 p1 p2
Bezier3To { controlPoint1 = p1, controlPoint2 = p2, curveEnd = NextPoint p3 }
-> bezier3 @k @i $ Cubic.Bezier p0 p1 p2 p3
-> bezier3 @k @i co $ Cubic.Bezier p0 p1 p2 p3
-- | 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.
@ -1081,17 +1084,19 @@ brushStrokeData :: forall k brushParams i arr
, D k ( I i ( 2 ) ) ~ D k ( 2 )
, Cross ( I i Double ) ( T ( I i ( 2 ) ) )
, Torsor ( T ( I i ( 2 ) ) ) ( I i ( 2 ) )
, Coercible ( I i ( 1 ) ) ( I i Double )
-- , Coercible ( I i ( 1 ) ) ( I i Double )
, Show brushParams
)
=> ( I i ( 1 ) `arr` I i ( 2 ) )
=> ( I i Double -> I i ( 1 ) )
-> ( I i ( 1 ) -> I i Double )
-> ( I i ( 1 ) `arr` I i ( 2 ) )
-- ^ path
-> ( I i ( 1 ) `arr` I i brushParams )
-- ^ brush parameters
-> ( I i brushParams `arr` Spline Closed () ( I i ( 2 ) ) )
-- ^ brush from parameters
-> ( I i ( 1 ) -> Seq ( I i ( 1 ) -> StrokeDatum k i ) )
brushStrokeData path params brush =
brushStrokeData co1 co2 path params brush =
\ t ->
let
dpath_t :: D k ( I i ( 1 ) ) ( I i ( 2 ) )
@ -1101,7 +1106,7 @@ brushStrokeData path params brush =
dbrush_params :: D k ( I i brushParams ) ( Spline Closed () ( I i ( 2 ) ) )
!dbrush_params = runD brush $ value @( I i Double ) @k @( I i ( 1 ) ) dparams_t
splines :: Seq ( D k ( I i brushParams ) ( I i ( 1 ) `arr` I i ( 2 ) ) )
!splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns @k @i ) dbrush_params
!splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns @k @i co2 ) dbrush_params
dbrushes_t :: Seq ( I i ( 1 ) -> D k ( I i ( 2 ) ) ( I i ( 2 ) ) )
!dbrushes_t = force $ fmap ( uncurryD @k . ( chain @(I i Double) @k dparams_t ) ) splines
-- This is the crucial use of the chain rule.
@ -1115,7 +1120,7 @@ brushStrokeData path params brush =
mkStrokeDatum dpath_t dbrush_t s =
let dbrush_t_s = dbrush_t s
dstroke = brushStroke @k dpath_t dbrush_t_s
( ee, 𝛿E𝛿sdcdt ) = envelopeEquation @k @i dstroke
( ee, 𝛿E𝛿sdcdt ) = envelopeEquation @k @i co1 dstroke
in -- trace
-- ( unlines
-- [ "envelopeEquation:"

View file

@ -79,25 +79,25 @@ class CurveOrder k where
line :: forall i b
. ( Module ( I i Double ) ( T b ), Torsor ( T b ) b
, D k ( I i ( 1 ) ) ~ D k ( 1 )
, Coercible ( I i ( 1 ) ) ( I i Double )
)
=> Segment b -> C k ( I i ( 1 ) ) b
=> ( I i ( 1 ) -> I i Double )
-> Segment b -> C k ( I i ( 1 ) ) b
-- | A quadratic Bézier curve, as a differentiable function.
bezier2 :: forall i b
. ( Module ( I i Double ) ( T b ), Torsor ( T b ) b
, D k ( I i ( 1 ) ) ~ D k ( 1 )
, Coercible ( I i ( 1 ) ) ( I i Double )
)
=> Quadratic.Bezier b -> C k ( I i ( 1 ) ) b
=> ( I i ( 1 ) -> I i Double )
-> Quadratic.Bezier b -> C k ( I i ( 1 ) ) b
-- | A cubic Bézier curve, as a differentiable function.
bezier3 :: forall i b
. ( Module ( I i Double ) ( T b ), Torsor ( T b ) b
, D k ( I i ( 1 ) ) ~ D k ( 1 )
, Coercible ( I i ( 1 ) ) ( I i Double )
)
=> Cubic.Bezier b -> C k ( I i ( 1 ) ) b
=> ( I i ( 1 ) -> I i Double )
-> Cubic.Bezier b -> C k ( I i ( 1 ) ) b
uncurryD :: D k a ~ D k ( 1 )
=> D k ( 1 ) ( C k a b ) -> a -> D k ( 2 ) b
@ -131,27 +131,27 @@ class CurveOrder k where
, D ( k - 1 ) ( I i ( 2 ) ) ~ D ( k - 1 ) ( 2 )
, D k ( I i ( 2 ) ) ~ D k ( 2 )
, Cross ( I i Double ) ( T ( I i ( 2 ) ) )
, Coercible ( I i Double ) ( I i ( 1 ) )
)
=> D k ( I i ( 2 ) ) ( I i ( 2 ) )
=> ( I i Double -> I i ( 1 ) )
-> D k ( I i ( 2 ) ) ( I i ( 2 ) )
-> ( D ( k - 1 ) ( I i ( 2 ) ) ( I i ( 1 ) )
, D ( k - 2 ) ( I i ( 2 ) ) ( T ( I i ( 2 ) ) ) )
instance CurveOrder 2 where
line ( Segment a b :: Segment b ) =
D \ ( coerce -> t ) ->
line co ( Segment a b :: Segment b ) =
D \ ( co -> t ) ->
D21 ( lerp @( T b ) t a b )
( a --> b )
origin
bezier2 ( bez :: Quadratic.Bezier b ) =
D \ ( coerce -> t ) ->
bezier2 co ( bez :: Quadratic.Bezier b ) =
D \ ( co -> t ) ->
D21 ( Quadratic.bezier @( T b ) bez t )
( Quadratic.bezier' bez t )
( Quadratic.bezier'' bez )
bezier3 ( bez :: Cubic.Bezier b ) =
D \ ( coerce -> t ) ->
bezier3 co ( bez :: Cubic.Bezier b ) =
D \ ( co -> t ) ->
D21 ( Cubic.bezier @( T b ) bez t )
( Cubic.bezier' bez t )
( Cubic.bezier'' bez t )
@ -171,12 +171,12 @@ instance CurveOrder 2 where
-- ∂²c/∂t∂s = ∂²b/∂t∂s
-- ∂²c/∂s² = ∂²b/∂s²
envelopeEquation ( D22 _ c_t c_s c_tt c_ts c_ss ) =
envelopeEquation co ( D22 _ c_t c_s c_tt c_ts c_ss ) =
let !ee = c_t × c_s
!ee_t = c_tt × c_s + c_t × c_ts
!ee_s = c_ts × c_s + c_t × c_ss
!𝛿E𝛿sdcdt = ee_s *^ c_t ^-^ ee_t *^ c_s
in ( D12 ( coerce ee ) ( T $ coerce ee_t ) ( T $ coerce ee_s )
in ( D12 ( co ee ) ( T $ co ee_t ) ( T $ co ee_s )
, D0 𝛿E𝛿sdcdt )
-- Computation of total derivative dc/dt:
--
@ -187,22 +187,22 @@ instance CurveOrder 2 where
instance CurveOrder 3 where
line ( Segment a b :: Segment b ) =
D \ ( coerce -> t ) ->
line co ( Segment a b :: Segment b ) =
D \ ( co -> t ) ->
D31 ( lerp @( T b ) t a b )
( a --> b )
origin
origin
bezier2 ( bez :: Quadratic.Bezier b ) =
D \ ( coerce -> t ) ->
bezier2 co ( bez :: Quadratic.Bezier b ) =
D \ ( co -> t ) ->
D31 ( Quadratic.bezier @( T b ) bez t )
( Quadratic.bezier' bez t )
( Quadratic.bezier'' bez )
origin
bezier3 ( bez :: Cubic.Bezier b ) =
D \ ( coerce -> t ) ->
bezier3 co ( bez :: Cubic.Bezier b ) =
D \ ( co -> t ) ->
D31 ( Cubic.bezier @( T b ) bez t )
( Cubic.bezier' bez t )
( Cubic.bezier'' bez t )
@ -219,7 +219,7 @@ instance CurveOrder 3 where
( p_tt ^+^ b_tt ) b_ts b_ss
( p_ttt ^+^ b_ttt ) b_tts b_tss b_sss
envelopeEquation
envelopeEquation co
( D32 _ c_t c_s
c_tt c_ts c_ss
c_ttt c_tts c_tss c_sss )
@ -242,7 +242,7 @@ instance CurveOrder 3 where
!𝛿E𝛿sdcdt_s = ee_ss *^ c_t ^+^ ee_s *^ c_ts
^-^ ( ee_ts *^ c_s ^+^ ee_t *^ c_ss )
in ( D22
( coerce ee )
( T $ coerce ee_t ) ( T $ coerce ee_s )
( T $ coerce ee_tt) ( T $ coerce ee_ts ) ( T $ coerce ee_ss )
( co ee )
( T $ co ee_t ) ( T $ co ee_s )
( T $ co ee_tt) ( T $ co ee_ts ) ( T $ co ee_ss )
, D12 𝛿E𝛿sdcdt ( T 𝛿E𝛿sdcdt_t ) ( T 𝛿E𝛿sdcdt_s ) )

View file

@ -13,11 +13,12 @@ module Math.Linear
, Fin(..), MFin(..)
, RepDim, RepresentableQ(..)
, Representable(..), injection, projection
, Vec(..), (!), find, zipIndices
, Vec(..), (!), find, zipIndices, unzip
) where
-- base
import Prelude hiding ( unzip )
import Data.Coerce
( coerce )
import Data.Kind
@ -160,3 +161,9 @@ zipIndices = zipIndices_ 1
zipIndices_ :: forall i. Word -> Vec i a -> [ ( Fin n, a ) ]
zipIndices_ _ VZ = []
zipIndices_ i (a `VS` as) = ( Fin i, a ) : zipIndices_ ( i + 1 ) as
unzip :: Vec n ( a, b ) -> ( Vec n a, Vec n b )
unzip VZ = ( VZ, VZ )
unzip ( ( a, b ) `VS` ab ) =
case unzip ab of
( as, bs ) -> ( a `VS` as, b `VS` bs )