From 7db3cbef332d5fd636aa866055fcda91e7d6a490 Mon Sep 17 00:00:00 2001 From: sheaf Date: Mon, 13 Mar 2023 22:09:15 +0100 Subject: [PATCH] more work into observability --- MetaBrush.cabal | 10 +- src/cusps/Main.hs | 309 +++++++----- src/cusps/Math/Interval/Abstract.hs | 465 ++++++++++++++++++ src/splines/Math/Algebra/Dual.hs | 2 + src/splines/Math/Bezier/Stroke.hs | 51 +- .../Math/Bezier/Stroke/EnvelopeEquation.hs | 52 +- src/splines/Math/Linear.hs | 9 +- 7 files changed, 724 insertions(+), 174 deletions(-) create mode 100644 src/cusps/Math/Interval/Abstract.hs diff --git a/MetaBrush.cabal b/MetaBrush.cabal index dbe18e8..2ef144b 100644 --- a/MetaBrush.cabal +++ b/MetaBrush.cabal @@ -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 diff --git a/src/cusps/Main.hs b/src/cusps/Main.hs index fb72764..278cc29 100644 --- a/src/cusps/Main.hs +++ b/src/cusps/Main.hs @@ -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 #-} diff --git a/src/cusps/Math/Interval/Abstract.hs b/src/cusps/Math/Interval/Abstract.hs new file mode 100644 index 0000000..2024cd3 --- /dev/null +++ b/src/cusps/Math/Interval/Abstract.hs @@ -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 ) ||] diff --git a/src/splines/Math/Algebra/Dual.hs b/src/splines/Math/Algebra/Dual.hs index 64a2402..c577809 100644 --- a/src/splines/Math/Algebra/Dual.hs +++ b/src/splines/Math/Algebra/Dual.hs @@ -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(..) diff --git a/src/splines/Math/Bezier/Stroke.hs b/src/splines/Math/Bezier/Stroke.hs index 9f20ba0..87750ed 100644 --- a/src/splines/Math/Bezier/Stroke.hs +++ b/src/splines/Math/Bezier/Stroke.hs @@ -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:" diff --git a/src/splines/Math/Bezier/Stroke/EnvelopeEquation.hs b/src/splines/Math/Bezier/Stroke/EnvelopeEquation.hs index 41c4b3c..594eb9c 100644 --- a/src/splines/Math/Bezier/Stroke/EnvelopeEquation.hs +++ b/src/splines/Math/Bezier/Stroke/EnvelopeEquation.hs @@ -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 ) ) diff --git a/src/splines/Math/Linear.hs b/src/splines/Math/Linear.hs index 3f89801..9deaf6f 100644 --- a/src/splines/Math/Linear.hs +++ b/src/splines/Math/Linear.hs @@ -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 )