diff --git a/MetaBrush.cabal b/MetaBrush.cabal index 2ef144b..9d00843 100644 --- a/MetaBrush.cabal +++ b/MetaBrush.cabal @@ -273,6 +273,9 @@ executable cusps hs-source-dirs: src/cusps + default-language: + Haskell2010 + main-is: Main.hs diff --git a/src/app/MetaBrush/Application.hs b/src/app/MetaBrush/Application.hs index 3a2697c..beb54d8 100644 --- a/src/app/MetaBrush/Application.hs +++ b/src/app/MetaBrush/Application.hs @@ -167,9 +167,9 @@ runApplication application = do -- ] -- } Spline - { splineStart = mkPoint ( ℝ2 10 -20 ) 2 1 0 + { splineStart = mkPoint ( ℝ2 0 0 ) 1 1 0 , splineCurves = OpenCurves $ Seq.fromList - [ LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 10 10 ) 10 5 ( pi / 4 ) ), curveData = invalidateCache undefined } + [ LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 100 0 ) 10 10 (pi / 2) ), curveData = invalidateCache undefined } --, LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 -10 10 ) 8 5 ( pi / 4 ) ), curveData = invalidateCache undefined } --, LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 -10 -20 ) 10 7 ( pi / 2 ) ), curveData = invalidateCache undefined } ] diff --git a/src/cusps/Main.hs b/src/cusps/Main.hs index 278cc29..ee71b18 100644 --- a/src/cusps/Main.hs +++ b/src/cusps/Main.hs @@ -1,6 +1,7 @@ {-# LANGUAGE AllowAmbiguousTypes #-} {-# LANGUAGE PartialTypeSignatures #-} {-# LANGUAGE PolyKinds #-} +{-# LANGUAGE RebindableSyntax #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE UndecidableInstances #-} @@ -10,7 +11,11 @@ module Main where -- base -import Prelude hiding ( Num(..), Fractional(..), Floating(..), (^) ) +import Prelude + hiding + ( Num(..), Fractional(..), Floating(..), (^), (/) + , fromInteger, fromRational + ) import qualified Prelude import GHC.Generics ( Generic ) @@ -29,21 +34,29 @@ import Data.Generics.Product.Typed -- MetaBrush import Math.Algebra.Dual +import qualified Math.Bezier.Quadratic as Quadratic +import qualified Math.Bezier.Cubic as Cubic import Math.Bezier.Spline import Math.Bezier.Stroke ( brushStrokeData, pathAndUsedParams ) import Math.Bezier.Stroke.EnvelopeEquation - ( StrokeDatum(..) ) + ( HasBézier(..), StrokeDatum(..) ) import Math.Differentiable ( I ) import Math.Linear - ( ℝ(..), T(..) + ( ℝ(..), T(..), Segment(..) , Fin(..) - , Representable(..) ) + , Representable(..) + ) import Math.Module - ( Module(..), Cross(..) ) + ( Module(..), Cross(..) + , lerp + ) import Math.Ring - ( AbelianGroup(..), Ring(..), Transcendental(..) ) + ( AbelianGroup(..), Ring(..) + , Field(..), Floating(..), Transcendental(..) + , ifThenElse + ) import Math.Interval.Abstract @@ -66,16 +79,19 @@ outlineFunction , 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 ) , HasType ( ℝ 2 ) ( PointData brushParams ) , Cross ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) , Module ( I i Double ) ( T ( I i brushParams ) ) + , Module Double (T brushParams) + , Torsor (T brushParams) brushParams + , Representable Double brushParams , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) ) , Torsor ( T ( I i brushParams ) ) ( I i brushParams ) , HasChainRule ( I i Double ) 3 ( I i ( ℝ 1 ) ) , HasChainRule ( I i Double ) 3 ( I i brushParams ) , Traversable ( D 3 brushParams ) , Traversable ( D 3 ( I i brushParams ) ) + , HasBézier 3 i ) => ( I i Double -> I i ( ℝ 1 ) ) -> ( I i ( ℝ 1 ) -> I i Double ) @@ -97,11 +113,15 @@ outlineFunction co1 co2 single brush brush_index sp0 crv = strokeData strokeData :: I i ( ℝ 1 ) -> I i ( ℝ 1 ) -> StrokeDatum 3 i strokeData = fmap ( `Seq.index` brush_index ) $ - brushStrokeData @3 @brushParams @i co1 co2 path params brush + brushStrokeData @3 @brushParams co1 co2 path params brush + boo :: ( HasType ( ℝ 2 ) ( PointData brushParams ) , Show brushParams + , Representable Double brushParams + , Module Double ( T brushParams ) , Module ( AI Double ) ( T ( AI brushParams ) ) + , Torsor ( T brushParams ) brushParams , Torsor ( T ( AI brushParams ) ) ( AI brushParams ) , HasChainRule ( AI Double ) 3 ( AI brushParams ) , Traversable ( D 3 brushParams ) @@ -118,9 +138,9 @@ boo x i y z 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 ) + boo ( ellipseBrush ( Pt . Val ) scaleI ) 1 + ( mkPoint ( ℝ2 0 0 ) 1.0111 1.0222 0 ) + ( LineTo { curveEnd = NextPoint $ mkPoint ( ℝ2 100 0 ) 10.0333 10.0444 ( pi Prelude./ 2 ) , curveData = () } ) where mkPoint :: ℝ 2 -> Double -> Double -> Double -> PointData ( ℝ 3 ) @@ -214,6 +234,28 @@ instance HasChainRule ( AI Double ) 3 ( AI ( ℝ 4 ) ) where [|| o ||] [|| p ||] [|| s ||] [|| df ||] [|| dg ||] ) +instance HasBézier 3 AI where + line co ( Segment a b :: Segment b ) = + D \ ( co -> t ) -> + D31 ( lerp @( T b ) t a b ) + ( a --> b ) + origin + origin + + bezier2 co ( bez :: Quadratic.Bezier b ) = + D \ ( co -> t ) -> + D31 ( EvalBez2 bez t ) + ( Quadratic.bezier' bez t ) + ( Quadratic.bezier'' bez ) + origin + + bezier3 co ( bez :: Cubic.Bezier b ) = + D \ ( co -> t ) -> + D31 ( EvalBez3 bez t ) + ( T $ EvalBez2 ( fmap unT $ Cubic.derivative ( fmap T bez ) ) t ) + ( Cubic.bezier'' bez t ) + ( Cubic.bezier''' bez ) + -------------------------------------------------------------------------------- -- Brushes (TODO copied from MetaBrush.Asset.Brushes) @@ -237,40 +279,60 @@ circleSpline p = sequenceA $ Bezier3To ( p κ -1 ) ( p 1 -κ ) BackToStart () {-# INLINE circleSpline #-} -ellipseBrush :: forall {t} (i :: t) k irec - . ( irec ~ I i ( ℝ 3 ) +ellipseBrush :: forall k irec + . ( irec ~ AI ( ℝ 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 + ( D k irec ( AI Double ) ) + ( D k irec ( AI ( ℝ 2 ) ) ) + , Module ( AI Double ) ( T ( AI Double ) ) + , HasChainRule ( AI Double ) k irec + , Representable ( AI Double ) irec , Applicative ( D k irec ) - , Transcendental ( D k irec ( I i Double ) ) + , Transcendental ( D k irec ( AI Double ) ) ) - => ( forall a. a -> I i a ) - -> ( Double -> I i Double -> I i Double ) - -> C k irec ( Spline 'Closed () ( I i ( ℝ 2 ) ) ) + => ( forall a. a -> AI a ) + -> ( Double -> AI Double -> AI Double ) + -> C k irec ( Spline 'Closed () ( AI ( ℝ 2 ) ) ) ellipseBrush mkI scaleI = D \ params -> - let a, b, phi :: D k irec ( I i Double ) + let a, b, phi :: D k irec ( AI 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 :: Double -> Double -> D k irec ( AI ( ℝ 2 ) ) mkPt x y - = let !a' = mul x a - !b' = mul y b + = let !x' = a `scaledBy` x + !y' = b `scaledBy` y +{- !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 + in ( x' * c - y' * s ) *^ e_x + ^+^ ( x' * s + y' * c ) *^ e_y +-} +-- {- + !r = sqrt $ x' ^ 2 + y' ^ 2 + !arctgt = atan ( y' / x' ) + -- a and b are always strictly positive, so we can compute + -- the quadrant using only x and y, which are constants. + !theta + | x > 0 + = arctgt + | x < 0 + = if y >= 0 then arctgt + pi else arctgt - pi + | otherwise + = if y >= 0 then 0.5 * pi else -0.5 * pi + !phi' = phi + theta + in + ( r * cos phi' ) *^ e_x + ^+^ ( r * sin phi' ) *^ e_y +-- -} + + in circleSpline @AI @k @( ℝ 3 ) @( ℝ 2 ) mkPt where - e_x, e_y :: D k irec ( I i ( ℝ 2 ) ) + e_x, e_y :: D k irec ( AI ( ℝ 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 ) + scaledBy :: D k irec ( AI Double ) -> Double -> D k irec ( AI Double ) + scaledBy d k = fmap ( scaleI k ) d {-# INLINEABLE ellipseBrush #-} diff --git a/src/cusps/Math/Interval/Abstract.hs b/src/cusps/Math/Interval/Abstract.hs index 2024cd3..0d54cb4 100644 --- a/src/cusps/Math/Interval/Abstract.hs +++ b/src/cusps/Math/Interval/Abstract.hs @@ -13,13 +13,21 @@ import Prelude hiding import qualified Prelude ( Num(..), Fractional(..) ) import Data.Coerce - ( Coercible, coerce ) + ( coerce ) +import Data.Dynamic + ( Dynamic, toDyn, fromDynamic, dynTypeRep ) import Data.Foldable ( toList ) import Data.List ( intersperse ) +import Data.Proxy + ( Proxy(..) ) +import Data.Typeable + ( Typeable, typeRep ) import GHC.Show - ( showCommaSpace ) + ( showSpace, showCommaSpace ) +import GHC.TypeNats + ( KnownNat ) -- acts import Data.Act @@ -28,8 +36,8 @@ import Data.Act -- containers import Data.Map ( Map ) -import qualified Data.Map as Map - ( (!) ) +import qualified Data.Map.Strict as Map + ( empty, insertLookupWithKey, lookup ) -- groups import Data.Group @@ -37,6 +45,8 @@ import Data.Group -- MetaBrush import Math.Algebra.Dual +import qualified Math.Bezier.Quadratic as Quadratic +import qualified Math.Bezier.Cubic as Cubic import Math.Bezier.Stroke.EnvelopeEquation ( StrokeDatum(..) ) import Math.Differentiable @@ -50,21 +60,24 @@ import Math.Linear import Math.Module ( Module(..), Cross(..) ) import Math.Ring - ( AbelianGroup(..), Ring(..), Field(..), Transcendental(..) + ( AbelianGroup(..), Ring(..), Field(..) + , Floating(..), Transcendental(..) , ViaPrelude(..) ) -------------------------------------------------------------------------------- -- EDSL for inspection -type instance D k (AI a) = D k (𝕀 a) +type instance D k ( AI a ) = D k ( 𝕀 a ) + -- TODO: I think we want a custom representation for D k (AI a), + -- so that we can recover sharing. 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 + Var :: ( Typeable a, Eq a, Show a ) => String -> Value a -- | A scalar or a vector. data SV a where @@ -87,16 +100,20 @@ data AI a where Recip :: AI Double -> AI Double Scale :: IsVec v => AI Double -> AI v -> AI v Cross :: AI ( ℝ 2 ) -> AI ( ℝ 2 ) -> AI Double + Sqrt :: AI Double -> AI Double Cos :: AI Double -> AI Double Sin :: AI Double -> AI Double + Atan :: 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 + EvalBez2 :: Quadratic.Bezier (AI a) -> AI Double -> AI a + EvalBez3 :: Cubic.Bezier (AI a) -> AI Double -> AI a + Tabulate :: ( KnownNat n, Eq ( ℝ n ), Show ( ℝ n ), Representable Double ( ℝ n ) ) => Vec n ( AI Double ) -> AI ( ℝ n ) + Index :: ( KnownNat n, 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 +instance ( Show a, Eq a, Typeable a ) => Show ( AI a ) where + showsPrec prec ai = case fst $ normalise Map.empty 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 @@ -113,8 +130,12 @@ instance ( Show a, Eq a ) => Show ( AI a ) where Recip iv -> showParen (prec > 10) $ showString " 1 / " . showsPrec 11 iv Coerce1 iv -> showsPrec prec iv Coerce2 iv -> showsPrec prec iv + EvalBez2 bez t -> showParen (prec > 10) $ showString "bez2" . showSpace . showsPrec 11 t . showSpace . showString "[" . showBez2 bez . showString "]" + EvalBez3 bez t -> showParen (prec > 10) $ showString "bez3" . showSpace . showsPrec 11 t . showSpace . showString "[" . showBez3 bez . showString "]" + Sqrt iv -> showParen (prec > 10) $ showString "sqrt " . showsPrec 11 iv Cos iv -> showParen (prec > 10) $ showString "cos " . showsPrec 11 iv Sin iv -> showParen (prec > 10) $ showString "sin " . showsPrec 11 iv + Atan iv -> showParen (prec > 10) $ showString "atan " . showsPrec 11 iv Pi -> showString "pi" Tabulate vs -> showString "(" . foldr (.) id ( intersperse showCommaSpace fields ) . showString ")" where @@ -123,6 +144,9 @@ instance ( Show a, Eq a ) => Show ( AI a ) where | v <- toList vs ] +showBez2 (Quadratic.Bezier p0 p1 p2 ) = foldr (.) id ( intersperse showCommaSpace $ map (showsPrec 0) [p0, p1, p2] ) +showBez3 ( Cubic.Bezier p0 p1 p2 p3) = foldr (.) id ( intersperse showCommaSpace $ map (showsPrec 0) [p0, p1, p2, p3] ) + infixl 6 :+: infixl 6 :-: infixl 7 :*: @@ -153,15 +177,29 @@ index_maybe expr i = case expr of -> Just $ v ! i Coerce1 v -> Just $ v + EvalBez2 bez t + | Just bez_i <- traverse ( `index_maybe` i ) bez + -> Just $ EvalBez2 bez_i t + EvalBez3 bez t + | Just bez_i <- traverse ( `index_maybe` i ) bez + -> Just $ EvalBez3 bez_i t _ -> 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 ) +normalise :: forall a. ( Eq a, Typeable a ) + => Map String Dynamic -> AI a -> ( AI a, Maybe ( 𝕀 a ) ) +normalise vars expr = case expr of + Pt ( Val v ) + -> ( expr, Just $ 𝕀 v v ) + Pt ( Var v ) + | Just d <- Map.lookup v vars + , Just a <- fromDynamic @( AI a ) d + , Just r <- snd $ normalise vars a + -> ( ival r, Just r ) + | otherwise + -> ( expr, Nothing ) Ival ( Val v1 ) ( Val v2 ) | v1 == v2 -> ( Pt ( Val v1 ), Just $ 𝕀 v1 v2 ) @@ -191,16 +229,16 @@ normalise expr = case expr of -> ( expr', Nothing ) where expr' = (:+:) sv iv1' iv2' - (iv1', mv1) = normalise iv1 - (iv2', mv2) = normalise iv2 + (iv1', mv1) = normalise vars iv1 + (iv2', mv2) = normalise vars iv2 (:-:) sv iv1 iv2 | Scalar <- sv , Just ( 𝕀 0 0 ) <- mv1 - -> normalise $ Negate sv iv2 + -> normalise vars $ Negate sv iv2 | Vector <- sv , Just o <- mv1 , T o == origin - -> normalise $ Negate sv iv2 + -> normalise vars $ Negate sv iv2 | Scalar <- sv , Just ( 𝕀 0 0 ) <- mv2 -> ( iv1', mv1 ) @@ -216,8 +254,8 @@ normalise expr = case expr of -> ( expr', Nothing ) where expr' = (:-:) sv iv1' iv2' - (iv1', mv1) = normalise iv1 - (iv2', mv2) = normalise iv2 + (iv1', mv1) = normalise vars iv1 + (iv2', mv2) = normalise vars iv2 Negate sv iv | Just v <- mv -> let r = case sv of { Scalar -> unT $ negate $ T v; Vector -> unT $ origin ^-^ T v } @@ -226,10 +264,14 @@ normalise expr = case expr of -> ( expr', Nothing ) where expr' = Negate sv iv' - (iv', mv) = normalise iv + (iv', mv) = normalise vars iv iv1 :*: iv2 | Just ( 𝕀 0 0 ) <- mv1 -> zero + | Just ( 𝕀 1 1 ) <- mv1 + -> (iv2', mv2) + | Just ( 𝕀 1 1 ) <- mv2 + -> (iv1', mv1) | Just ( 𝕀 0 0 ) <- mv2 -> zero | Just v1 <- mv1 @@ -241,8 +283,8 @@ normalise expr = case expr of where zero = ( ival $ 𝕀 0 0, Just $ 𝕀 0 0 ) expr' = (:*:) iv1' iv2' - (iv1', mv1) = normalise iv1 - (iv2', mv2) = normalise iv2 + (iv1', mv1) = normalise vars iv1 + (iv2', mv2) = normalise vars iv2 iv :^: n | Just v <- mv -> let r = v ^ n @@ -251,7 +293,7 @@ normalise expr = case expr of -> ( expr', Nothing ) where expr' = iv' :^: n - (iv', mv) = normalise iv + (iv', mv) = normalise vars iv Recip iv | Just v <- mv -> let r = recip v @@ -260,7 +302,7 @@ normalise expr = case expr of -> ( expr', Nothing ) where expr' = Recip iv' - (iv', mv) = normalise iv + (iv', mv) = normalise vars iv Scale k iv | Just ( 𝕀 0 0 ) <- mk -> ( ival $ unT $ origin, Just $ unT origin ) @@ -277,8 +319,8 @@ normalise expr = case expr of -> ( expr', Nothing ) where expr' = Scale k' iv' - (k' , mk) = normalise k - (iv', mv) = normalise iv + (k' , mk) = normalise vars k + (iv', mv) = normalise vars iv Cross iv1 iv2 | Just v1 <- mv1 , Just v2 <- mv2 @@ -288,8 +330,17 @@ normalise expr = case expr of -> ( expr', Nothing ) where expr' = Cross iv1' iv2' - (iv1', mv1) = normalise iv1 - (iv2', mv2) = normalise iv2 + (iv1', mv1) = normalise vars iv1 + (iv2', mv2) = normalise vars iv2 + Sqrt iv + | Just v <- mv + -> let r = sqrt v + in ( ival r, Just r ) + | otherwise + -> ( expr', Nothing ) + where + expr' = Sqrt iv' + (iv', mv) = normalise vars iv Cos iv | Just v <- mv -> let r = cos v @@ -298,7 +349,7 @@ normalise expr = case expr of -> ( expr', Nothing ) where expr' = Cos iv' - (iv', mv) = normalise iv + (iv', mv) = normalise vars iv Sin iv | Just v <- mv -> let r = sin v @@ -307,7 +358,16 @@ normalise expr = case expr of -> ( expr', Nothing ) where expr' = Sin iv' - (iv', mv) = normalise iv + (iv', mv) = normalise vars iv + Atan iv + | Just v <- mv + -> let r = atan v + in ( ival r, Just r ) + | otherwise + -> ( expr', Nothing ) + where + expr' = Atan iv' + (iv', mv) = normalise vars iv Pi -> ( Pi, Just pi ) Coerce1 iv | Just v <- mv @@ -317,7 +377,7 @@ normalise expr = case expr of -> ( expr', Nothing ) where expr' = Coerce1 iv' - (iv', mv) = normalise iv + (iv', mv) = normalise vars iv Coerce2 iv | Just v <- mv -> let r = coerce v @@ -326,7 +386,32 @@ normalise expr = case expr of -> ( expr', Nothing ) where expr' = Coerce2 iv' - (iv', mv) = normalise iv + (iv', mv) = normalise vars iv + EvalBez2 @v (Quadratic.Bezier p0 p1 p2) t + -> ( EvalBez2 (Quadratic.Bezier p0' p1' p2') t' + , Nothing + -- Quadratic.bezier + -- <$> (Quadratic.Bezier <$> m0 <*> m1 <*> m2) + -- <*> mt + ) + where + (p0', m0) = normalise vars p0 + (p1', m1) = normalise vars p1 + (p2', m2) = normalise vars p2 + (t' , mt) = normalise vars t + EvalBez3 @v (Cubic.Bezier p0 p1 p2 p3) t + -> ( EvalBez3 (Cubic.Bezier p0' p1' p2' p3') t' + , Nothing + -- Cubic.bezier + -- <$> (Cubic.Bezier <$> m0 <*> m1 <*> m2 <*> m3) + -- <*> mt + ) + where + (p0', m0) = normalise vars p0 + (p1', m1) = normalise vars p1 + (p2', m2) = normalise vars p2 + (p3', m3) = normalise vars p3 + (t' , mt) = normalise vars t Tabulate @n v | Just xs <- sequence mxs -> let r = tabulate ( xs ! ) @@ -337,10 +422,10 @@ normalise expr = case expr of expr' = Tabulate v' v' :: Vec n ( AI Double ) mxs :: Vec n ( Maybe ( 𝕀 Double ) ) - (v', mxs) = unzip $ fmap normalise v + (v', mxs) = unzip $ fmap ( normalise vars ) v Index iv i | Just r <- index_maybe iv' i - -> normalise r + -> normalise vars r | Just v <- mv -> let r = v `index` i in ( ival r, Just r ) @@ -348,13 +433,21 @@ normalise expr = case expr of -> ( expr', Nothing ) where expr' = Index iv' i - (iv', mv) = normalise iv + (iv', mv) = normalise vars iv -val :: Map String Double -> Value a -> a +val :: forall a. Typeable a => Map String Dynamic -> Value a -> a val _ ( Val a ) = a -val vars ( Var v ) = vars Map.! v +val vars ( Var v ) = + case Map.lookup v vars of + Nothing -> error $ "Variable '" ++ v ++ "' not in environment." + Just d -> + case fromDynamic d of + Nothing -> error $ unlines [ "Incorrect type for variable '" ++ v ++ "'." + , " Expected: " ++ show ( typeRep $ Proxy @a ) + , " Actual: " ++ show ( dynTypeRep d ) ] + Just a -> a -eval :: Map String Double -> AI a -> 𝕀 a +eval :: Typeable a => Map String Dynamic -> 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) @@ -369,11 +462,15 @@ eval vars = \ case 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 ) + Sqrt iv -> sqrt $ eval vars iv Cos iv -> cos $ eval vars iv Sin iv -> sin $ eval vars iv + Atan iv -> atan $ eval vars iv Pi -> pi Coerce1 a -> coerce $ eval vars a Coerce2 a -> coerce $ eval vars a + EvalBez2 bez t -> error "TODO" -- Quadratic.bezier (fmap (eval vars) bez) (eval vars t) + EvalBez3 bez t -> error "TODO" -- Cubic.bezier (fmap (eval vars) bez) (eval vars t) Tabulate v -> let !w = fmap ( eval vars ) v in tabulate ( w ! ) Index v i -> let !w = eval vars v @@ -390,7 +487,7 @@ instance Prelude.Num ( AI Double ) where instance Prelude.Fractional ( AI Double ) where recip = Recip - (/) = error "No division for abstract intervals" + a / b = a :*: Recip b fromRational = Pt . Val . fromRational instance Ring ( AI Double ) where @@ -404,10 +501,14 @@ deriving via ViaPrelude ( AI Double ) deriving via ViaPrelude ( AI Double ) instance Field ( AI Double ) +instance Floating ( AI Double ) where + sqrt = Sqrt + instance Transcendental ( AI Double ) where pi = Pi cos = Cos sin = Sin + atan = Atan class ( Eq v, Show v, Module ( 𝕀 Double ) ( T ( 𝕀 v ) ), Group ( T ( 𝕀 v ) ) ) => IsVec v diff --git a/src/metabrushes/MetaBrush/Asset/Brushes.hs b/src/metabrushes/MetaBrush/Asset/Brushes.hs index b10e50f..70a052a 100644 --- a/src/metabrushes/MetaBrush/Asset/Brushes.hs +++ b/src/metabrushes/MetaBrush/Asset/Brushes.hs @@ -1,6 +1,7 @@ {-# LANGUAGE AllowAmbiguousTypes #-} {-# LANGUAGE OverloadedStrings #-} {-# LANGUAGE PolyKinds #-} +{-# LANGUAGE RebindableSyntax #-} {-# LANGUAGE ScopedTypeVariables #-} module MetaBrush.Asset.Brushes @@ -12,9 +13,9 @@ module MetaBrush.Asset.Brushes -- base import Prelude hiding - ( Num(..), Floating(..) ) + ( Num(..), Floating(..), (^), (/), fromInteger, fromRational ) import GHC.Exts - ( Proxy# ) + ( Proxy#, fromString ) -- containers import qualified Data.Sequence as Seq @@ -117,16 +118,16 @@ circleBrush _ mkI = let r :: D k irec ( I i Double ) r = runD ( var @_ @k ( Fin 1 ) ) params mkPt :: Double -> Double -> D k irec ( I i ( ℝ 2 ) ) - mkPt ( kon -> x ) ( kon -> y ) - = ( x * r ) *^ e_x - ^+^ ( y * r ) *^ e_y + mkPt x y + = ( r `scaledBy` x ) *^ e_x + ^+^ ( r `scaledBy` y ) *^ e_y in circleSpline @i @k @( Record CircleBrushFields ) @( ℝ 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 - kon = konst @( I i Double ) @k @irec . mkI + scaledBy d x = fmap ( mkI x * ) d {-# INLINEABLE circleBrush #-} ellipseBrush :: forall {t} (i :: t) k irec @@ -152,14 +153,37 @@ ellipseBrush _ mkI = b = runD ( var @_ @k ( Fin 2 ) ) params phi = runD ( var @_ @k ( Fin 3 ) ) params mkPt :: Double -> Double -> D k irec ( 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 + mkPt x y + = let !x' = a `scaledBy` x + !y' = b `scaledBy` y +-- {- + in + ( x' * cos phi - y' * sin phi ) *^ e_x + ^+^ ( y' * cos phi + x' * sin phi ) *^ e_y +-- -} +{- + r = sqrt ( x' ^ 2 + y' ^ 2 ) + arctgt = atan ( y' / x' ) + -- a and b are always strictly positive, so we can compute + -- the quadrant using only x and y, which are constants. + !theta + | x > 0 + = arctgt + | x < 0 + = if y >= 0 then arctgt + pi else arctgt - pi + | otherwise + = if y >= 0 then 0.5 * pi else -0.5 * pi + !phi' = phi + theta + in + ( r * cos phi' ) *^ e_x + ^+^ ( r * sin phi' ) *^ e_y +-} + in circleSpline @i @k @( Record EllipseBrushFields ) @( ℝ 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 - kon = konst @( I i Double ) @k @irec . mkI + scaledBy d x = fmap ( mkI x * ) d {-# INLINEABLE ellipseBrush #-} diff --git a/src/splines/Math/Algebra/Dual.hs b/src/splines/Math/Algebra/Dual.hs index c577809..7406ce3 100644 --- a/src/splines/Math/Algebra/Dual.hs +++ b/src/splines/Math/Algebra/Dual.hs @@ -25,15 +25,14 @@ module Math.Algebra.Dual ) where -- base -import Prelude hiding ( Num(..), Floating(..), (^) ) -import Control.Applicative - ( liftA2 ) +import Prelude + hiding + ( Num(..), Floating(..) + , (^), recip, fromInteger, fromRational ) import Data.Coerce ( coerce ) import Data.Kind ( Type ) -import Data.Monoid - ( Ap(..) ) import GHC.TypeNats ( Nat ) @@ -238,24 +237,37 @@ deriving via ApAp2 3 ( ℝ 4 ) r -------------------------------------------------------------------------------- -- Ring instances. --- TODO: should also implement a power operation specifically, --- as this is important for interval arithmethic. +d1pow :: Ring r => Word -> r -> D1𝔸1 r +d1pow i x = + let !j = fromInteger $ fromIntegral i + in D11 + ( x ^ i ) + ( T $ j * x ^ (i - 1) ) +{-# INLINE d1pow #-} +d2pow :: Ring r => Word -> r -> D2𝔸1 r +d2pow i x = + let !j = fromInteger $ fromIntegral i + in D21 + ( x ^ i ) + ( T $ j * x ^ (i - 1) ) + ( T $ j * (j - 1) * x ^ ( i - 2 ) ) +{-# INLINE d2pow #-} +d3pow :: Ring r => Word -> r -> D3𝔸1 r +d3pow i x = + let !j = fromInteger $ fromIntegral i + in D31 + ( x ^ i ) + ( T $ j * x ^ (i - 1) ) + ( T $ j * (j - 1) * x ^ ( i - 2 ) ) + ( T $ j * (j - 1) * (j - 2) * x ^ ( i - 3 ) ) +{-# INLINE d3pow #-} 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 ( D1𝔸2 r ) where +--instance Ring r => Ring ( D1𝔸3 r ) where +--instance Ring r => Ring ( D1𝔸4 r ) where instance Ring r => Ring ( D2𝔸1 r ) where !dr1 * !dr2 = @@ -270,39 +282,18 @@ instance Ring r => Ring ( D2𝔸1 r ) where $$( prodRuleQ [|| o ||] [|| p ||] [|| m ||] [|| dr1 ||] [|| dr2 ||] ) + df ^ i = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d2pow @r i ( _D21_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Ring ( D2𝔸1 Double ) #-} {-# SPECIALISE instance Ring ( D2𝔸1 ( 𝕀 Double ) ) #-} -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 ||] ) - {-# SPECIALISE instance Ring ( D3𝔸1 Double ) #-} - {-# SPECIALISE instance Ring ( D3𝔸1 ( 𝕀 Double ) ) #-} - ---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 @@ -316,39 +307,18 @@ instance Ring r => Ring ( D2𝔸2 r ) where $$( prodRuleQ [|| o ||] [|| p ||] [|| m ||] [|| dr1 ||] [|| dr2 ||] ) + df ^ i = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d2pow @r i ( _D22_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Ring ( D2𝔸2 Double ) #-} {-# SPECIALISE instance Ring ( D2𝔸2 ( 𝕀 Double ) ) #-} -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 ||] ) - {-# SPECIALISE instance Ring ( D3𝔸2 Double ) #-} - {-# SPECIALISE instance Ring ( D3𝔸2 ( 𝕀 Double ) ) #-} - ---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 @@ -362,39 +332,18 @@ instance Ring r => Ring ( D2𝔸3 r ) where $$( prodRuleQ [|| o ||] [|| p ||] [|| m ||] [|| dr1 ||] [|| dr2 ||] ) + df ^ i = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d2pow @r i ( _D23_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Ring ( D2𝔸3 Double ) #-} {-# SPECIALISE instance Ring ( D2𝔸3 ( 𝕀 Double ) ) #-} -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 ||] ) - {-# SPECIALISE instance Ring ( D3𝔸3 Double ) #-} - {-# SPECIALISE instance Ring ( D3𝔸3 ( 𝕀 Double ) ) #-} - ---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 @@ -408,9 +357,93 @@ instance Ring r => Ring ( D2𝔸4 r ) where $$( prodRuleQ [|| o ||] [|| p ||] [|| m ||] [|| dr1 ||] [|| dr2 ||] ) + df ^ i = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d2pow @r i ( _D24_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Ring ( D2𝔸4 Double ) #-} {-# SPECIALISE instance Ring ( D2𝔸4 ( 𝕀 Double ) ) #-} +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 ||] ) + df ^ i = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d3pow @r i ( _D31_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Ring ( D3𝔸1 Double ) #-} + {-# SPECIALISE instance Ring ( D3𝔸1 ( 𝕀 Double ) ) #-} + +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 ||] ) + df ^ i = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d3pow @r i ( _D32_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Ring ( D3𝔸2 Double ) #-} + {-# SPECIALISE instance Ring ( D3𝔸2 ( 𝕀 Double ) ) #-} + +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 ||] ) + df ^ i = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d3pow @r i ( _D33_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Ring ( D3𝔸3 Double ) #-} + {-# SPECIALISE instance Ring ( D3𝔸3 ( 𝕀 Double ) ) #-} + instance Ring r => Ring ( D3𝔸4 r ) where !dr1 * !dr2 = let @@ -424,11 +457,36 @@ instance Ring r => Ring ( D3𝔸4 r ) where $$( prodRuleQ [|| o ||] [|| p ||] [|| m ||] [|| dr1 ||] [|| dr2 ||] ) + df ^ i = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d2pow @r i ( _D34_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Ring ( D3𝔸4 Double ) #-} {-# SPECIALISE instance Ring ( D3𝔸4 ( 𝕀 Double ) ) #-} -------------------------------------------------------------------------------- --- Field & transcendental instances +-- Field, floating & transcendental instances + +d1recip :: Field r => r -> D1𝔸1 r +d1recip x = + let !d = recip x + in D11 d ( T $ -1 * d ^ 2 ) +{-# INLINE d1recip #-} +d2recip :: Field r => r -> D2𝔸1 r +d2recip x = + let !d = recip x + in D21 d ( T $ -1 * d ^ 2 ) ( T $ 2 * d ^ 3 ) +{-# INLINE d2recip #-} +d3recip :: Field r => r -> D3𝔸1 r +d3recip x = + let !d = recip x + in D31 d ( T $ -1 * d ^ 2 ) ( T $ 2 * d ^ 3 ) ( T $ -6 * d ^ 4 ) +{-# INLINE d3recip #-} deriving newtype instance Field r => Field ( D𝔸0 r ) @@ -437,16 +495,232 @@ deriving newtype instance Field r => Field ( D𝔸0 r ) --instance Field r => Field ( D1𝔸3 r ) where --instance Field r => Field ( D1𝔸4 r ) where instance Field r => Field ( D2𝔸1 r ) where + fromRational q = konst @Double @2 @( ℝ 1 ) ( fromRational q ) + recip df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d2recip @r ( _D21_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Field ( D2𝔸1 Double ) #-} + {-# SPECIALISE instance Field ( D2𝔸1 ( 𝕀 Double ) ) #-} instance Field r => Field ( D2𝔸2 r ) where + fromRational q = konst @Double @2 @( ℝ 2 ) ( fromRational q ) + recip df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d2recip @r ( _D22_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Field ( D2𝔸2 Double ) #-} + {-# SPECIALISE instance Field ( D2𝔸2 ( 𝕀 Double ) ) #-} instance Field r => Field ( D2𝔸3 r ) where + fromRational q = konst @Double @2 @( ℝ 3 ) ( fromRational q ) + recip df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d2recip @r ( _D23_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Field ( D2𝔸3 Double ) #-} + {-# SPECIALISE instance Field ( D2𝔸3 ( 𝕀 Double ) ) #-} instance Field r => Field ( D2𝔸4 r ) where + fromRational q = konst @Double @2 @( ℝ 4 ) ( fromRational q ) + recip df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d2recip @r ( _D24_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Field ( D2𝔸4 Double ) #-} + {-# SPECIALISE instance Field ( D2𝔸4 ( 𝕀 Double ) ) #-} instance Field r => Field ( D3𝔸1 r ) where + fromRational q = konst @Double @3 @( ℝ 1 ) ( fromRational q ) + recip df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d3recip @r ( _D31_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Field ( D3𝔸1 Double ) #-} + {-# SPECIALISE instance Field ( D3𝔸1 ( 𝕀 Double ) ) #-} instance Field r => Field ( D3𝔸2 r ) where + fromRational q = konst @Double @3 @( ℝ 2 ) ( fromRational q ) + recip df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d3recip @r ( _D32_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Field ( D3𝔸2 Double ) #-} + {-# SPECIALISE instance Field ( D3𝔸2 ( 𝕀 Double ) ) #-} instance Field r => Field ( D3𝔸3 r ) where + fromRational q = konst @Double @3 @( ℝ 3 ) ( fromRational q ) + recip df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d3recip @r ( _D33_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Field ( D3𝔸3 Double ) #-} + {-# SPECIALISE instance Field ( D3𝔸3 ( 𝕀 Double ) ) #-} instance Field r => Field ( D3𝔸4 r ) where --- TODO + fromRational q = konst @Double @3 @( ℝ 4 ) ( fromRational q ) + recip df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d3recip @r ( _D34_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Field ( D3𝔸4 Double ) #-} + {-# SPECIALISE instance Field ( D3𝔸4 ( 𝕀 Double ) ) #-} -d1sin, d1cos :: Transcendental r => r -> D1𝔸1 r +d1sqrt :: Floating r => r -> D1𝔸1 r +d1sqrt x = + let !v = sqrt x + !d = recip v + in D11 v ( T $ 0.5 * d ) +{-# INLINE d1sqrt #-} +d2sqrt :: Floating r => r -> D2𝔸1 r +d2sqrt x = + let !v = sqrt x + !d = recip v + in D21 v ( T $ 0.5 * d ) ( T $ -0.25 * d ^ 2 ) +{-# INLINE d2sqrt #-} +d3sqrt :: Floating r => r -> D3𝔸1 r +d3sqrt x = + let !v = sqrt x + !d = recip v + in D31 v ( T $ 0.5 * d ) ( T $ -0.25 * d ^ 2 ) ( T $ 0.375 * d ^ 3 ) +{-# INLINE d3sqrt #-} + +deriving newtype instance Floating r => Floating ( D𝔸0 r ) +--instance Floating r => Floating ( D1𝔸1 r ) where +--instance Floating r => Floating ( D1𝔸2 r ) where +--instance Floating r => Floating ( D1𝔸3 r ) where +--instance Floating r => Floating ( D1𝔸4 r ) where +instance Floating r => Floating ( D2𝔸1 r ) where + sqrt df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d2sqrt @r ( _D21_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Floating ( D2𝔸1 Double ) #-} + {-# SPECIALISE instance Floating ( D2𝔸1 ( 𝕀 Double ) ) #-} +instance Floating r => Floating ( D2𝔸2 r ) where + sqrt df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d2sqrt @r ( _D22_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Floating ( D2𝔸2 Double ) #-} + {-# SPECIALISE instance Floating ( D2𝔸2 ( 𝕀 Double ) ) #-} +instance Floating r => Floating ( D2𝔸3 r ) where + sqrt df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d2sqrt @r ( _D23_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Floating ( D2𝔸3 Double ) #-} + {-# SPECIALISE instance Floating ( D2𝔸3 ( 𝕀 Double ) ) #-} +instance Floating r => Floating ( D2𝔸4 r ) where + sqrt df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d2sqrt @r ( _D24_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Floating ( D2𝔸4 Double ) #-} + {-# SPECIALISE instance Floating ( D2𝔸4 ( 𝕀 Double ) ) #-} +instance Floating r => Floating ( D3𝔸1 r ) where + sqrt df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d3sqrt @r ( _D31_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Floating ( D3𝔸1 Double ) #-} + {-# SPECIALISE instance Floating ( D3𝔸1 ( 𝕀 Double ) ) #-} +instance Floating r => Floating ( D3𝔸2 r ) where + sqrt df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d3sqrt @r ( _D32_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Floating ( D3𝔸2 Double ) #-} + {-# SPECIALISE instance Floating ( D3𝔸2 ( 𝕀 Double ) ) #-} +instance Floating r => Floating ( D3𝔸3 r ) where + sqrt df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d3sqrt @r ( _D33_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Floating ( D3𝔸3 Double ) #-} + {-# SPECIALISE instance Floating ( D3𝔸3 ( 𝕀 Double ) ) #-} +instance Floating r => Floating ( D3𝔸4 r ) where + sqrt df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d3sqrt @r ( _D34_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Floating ( D3𝔸4 Double ) #-} + {-# SPECIALISE instance Floating ( D3𝔸4 ( 𝕀 Double ) ) #-} + +d1sin, d1cos, d1atan :: Transcendental r => r -> D1𝔸1 r d1sin x = let !s = sin x !c = cos x @@ -455,8 +729,16 @@ d1cos x = let !s = sin x !c = cos x in D11 c ( T -s ) +d1atan x + = let !d = recip $ 1 + x ^ 2 + in D11 + ( atan x ) + ( T $ d ) +{-# INLINE d1sin #-} +{-# INLINE d1cos #-} +{-# INLINE d1atan #-} -d2sin, d2cos :: Transcendental r => r -> D2𝔸1 r +d2sin, d2cos, d2atan :: Transcendental r => r -> D2𝔸1 r d2sin x = let !s = sin x !c = cos x @@ -465,8 +747,17 @@ d2cos x = let !s = sin x !c = cos x in D21 c ( T -s ) ( T -c ) +d2atan x + = let !d = recip $ 1 + x ^ 2 + in D21 + ( atan x ) + ( T $ d ) + ( T $ -2 * x * d ^ 2 ) +{-# INLINE d2sin #-} +{-# INLINE d2cos #-} +{-# INLINE d2atan #-} -d3sin, d3cos :: Transcendental r => r -> D3𝔸1 r +d3sin, d3cos, d3atan :: Transcendental r => r -> D3𝔸1 r d3sin x = let !s = sin x !c = cos x @@ -475,6 +766,16 @@ d3cos x = let !s = sin x !c = cos x in D31 c ( T -s ) ( T -c ) ( T s ) +d3atan x + = let !d = recip $ 1 + x ^ 2 + in D31 + ( atan x ) + ( T $ d ) + ( T $ -2 * x * d ^ 2 ) + ( T $ ( -2 + 6 * x ^ 2 ) * d ^ 3 ) +{-# INLINE d3sin #-} +{-# INLINE d3cos #-} +{-# INLINE d3atan #-} deriving newtype instance Transcendental r => Transcendental ( D𝔸0 r ) --instance Transcendental r => Transcendental ( D1𝔸1 r ) where @@ -501,6 +802,17 @@ instance Transcendental r => Transcendental ( D2𝔸1 r ) where dg = d2cos @r ( _D21_v df ) in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) + atan df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d2atan @r ( _D21_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Transcendental ( D2𝔸1 Double ) #-} + {-# SPECIALISE instance Transcendental ( D2𝔸1 ( 𝕀 Double ) ) #-} instance Transcendental r => Transcendental ( D2𝔸2 r ) where pi = konst @Double @2 @( ℝ 2 ) pi sin df = @@ -521,6 +833,17 @@ instance Transcendental r => Transcendental ( D2𝔸2 r ) where dg = d2cos @r ( _D22_v df ) in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) + atan df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d2atan @r ( _D22_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Transcendental ( D2𝔸2 Double ) #-} + {-# SPECIALISE instance Transcendental ( D2𝔸2 ( 𝕀 Double ) ) #-} instance Transcendental r => Transcendental ( D2𝔸3 r ) where pi = konst @Double @2 @( ℝ 3 ) pi sin df = @@ -541,6 +864,17 @@ instance Transcendental r => Transcendental ( D2𝔸3 r ) where dg = d2cos @r ( _D23_v df ) in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) + atan df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d2atan @r ( _D23_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Transcendental ( D2𝔸3 Double ) #-} + {-# SPECIALISE instance Transcendental ( D2𝔸3 ( 𝕀 Double ) ) #-} instance Transcendental r => Transcendental ( D2𝔸4 r ) where pi = konst @Double @2 @( ℝ 4 ) pi @@ -562,6 +896,17 @@ instance Transcendental r => Transcendental ( D2𝔸4 r ) where dg = d2cos @r ( _D24_v df ) in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) + atan df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d2atan @r ( _D24_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Transcendental ( D2𝔸4 Double ) #-} + {-# SPECIALISE instance Transcendental ( D2𝔸4 ( 𝕀 Double ) ) #-} instance Transcendental r => Transcendental ( D3𝔸1 r ) where pi = konst @Double @3 @( ℝ 1 ) pi @@ -583,6 +928,17 @@ instance Transcendental r => Transcendental ( D3𝔸1 r ) where dg = d3cos @r ( _D31_v df ) in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) + atan df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d3atan @r ( _D31_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Transcendental ( D3𝔸1 Double ) #-} + {-# SPECIALISE instance Transcendental ( D3𝔸1 ( 𝕀 Double ) ) #-} instance Transcendental r => Transcendental ( D3𝔸2 r ) where pi = konst @Double @3 @( ℝ 2 ) pi @@ -604,6 +960,17 @@ instance Transcendental r => Transcendental ( D3𝔸2 r ) where dg = d3cos @r ( _D32_v df ) in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) + atan df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d3atan @r ( _D32_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Transcendental ( D3𝔸2 Double ) #-} + {-# SPECIALISE instance Transcendental ( D3𝔸2 ( 𝕀 Double ) ) #-} instance Transcendental r => Transcendental ( D3𝔸3 r ) where pi = konst @Double @3 @( ℝ 3 ) pi @@ -625,6 +992,17 @@ instance Transcendental r => Transcendental ( D3𝔸3 r ) where dg = d3cos @r ( _D33_v df ) in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) + atan df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d3atan @r ( _D33_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Transcendental ( D3𝔸3 Double ) #-} + {-# SPECIALISE instance Transcendental ( D3𝔸3 ( 𝕀 Double ) ) #-} instance Transcendental r => Transcendental ( D3𝔸4 r ) where pi = konst @Double @3 @( ℝ 4 ) pi @@ -646,6 +1024,17 @@ instance Transcendental r => Transcendental ( D3𝔸4 r ) where dg = d3cos @r ( _D34_v df ) in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) + atan df = + let + fromInt = fromInteger @r + add = (+) @r + times = (*) @r + pow = (^) @r + dg = d3atan @r ( _D34_v df ) + in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] + [|| df ||] [|| dg ||] ) + {-# SPECIALISE instance Transcendental ( D3𝔸4 Double ) #-} + {-# SPECIALISE instance Transcendental ( D3𝔸4 ( 𝕀 Double ) ) #-} -------------------------------------------------------------------------------- -- HasChainRule instances. diff --git a/src/splines/Math/Bezier/Stroke.hs b/src/splines/Math/Bezier/Stroke.hs index 27de833..21e3512 100644 --- a/src/splines/Math/Bezier/Stroke.hs +++ b/src/splines/Math/Bezier/Stroke.hs @@ -22,8 +22,6 @@ module Math.Bezier.Stroke -- base import Control.Arrow ( first, (***) ) -import Control.Applicative - ( Applicative(..) ) import Control.Monad ( unless ) import Control.Monad.ST @@ -31,7 +29,7 @@ import Control.Monad.ST import Data.Bifunctor ( Bifunctor(bimap) ) import Data.Coerce - ( Coercible, coerce ) + ( coerce ) import Data.Foldable ( for_ ) import Data.Functor.Identity @@ -1273,7 +1271,7 @@ intervalNewtonGS precondMethod minWidth eqs = , 𝛿E𝛿sdcdt = D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) } <- ( eqs t `Seq.index` i ) s - , StrokeDatum { 𝛿E𝛿sdcdt = D12 ( T f_mid ) ( T ( T f_t_mid ) ) ( T ( T f_s_mid ) ) } + , StrokeDatum { 𝛿E𝛿sdcdt = D12 ( T f_mid ) ( T ( T _f_t_mid ) ) ( T ( T _f_s_mid ) ) } <- ( eqs i_t_mid `Seq.index` i ) i_s_mid = if | inf ee < ℝ1 0 , sup ee > ℝ1 0 diff --git a/src/splines/Math/Interval.hs b/src/splines/Math/Interval.hs index a904273..f5cdfd4 100644 --- a/src/splines/Math/Interval.hs +++ b/src/splines/Math/Interval.hs @@ -19,8 +19,6 @@ module Math.Interval -- base import Prelude hiding ( Num(..) ) -import Data.List.NonEmpty - ( NonEmpty(..) ) -- acts import Data.Act diff --git a/src/splines/Math/Interval/Internal.hs b/src/splines/Math/Interval/Internal.hs index a8d692a..31a135d 100644 --- a/src/splines/Math/Interval/Internal.hs +++ b/src/splines/Math/Interval/Internal.hs @@ -116,12 +116,24 @@ instance Prelude.Fractional ( 𝕀 Double ) where | otherwise = error "BAD interval recip; should use extendedRecip instead" -- #endif - _ / _ = error "TODO: interval division is not implemented" + 𝕀 x_lo x_hi / 𝕀 y_lo y_hi +-- #ifdef ASSERTS + | y_lo >= 0 || y_hi <= 0 +-- #endif + = 𝕀 ( fst $ divI x_lo y_hi ) ( snd $ divI x_hi y_lo ) +-- #ifdef ASSERTS + | otherwise + = error "BAD interval division; should use extendedRecip instead" +-- #endif + +instance Floating ( 𝕀 Double ) where + sqrt = withHW Prelude.sqrt instance Transcendental ( 𝕀 Double ) where pi = 𝕀 3.141592653589793 3.1415926535897936 cos = withHW Prelude.cos sin = withHW Prelude.sin + atan = withHW Prelude.atan {-# INLINE withHW #-} withHW :: (Interval.Interval a -> Interval.Interval b) -> 𝕀 a -> 𝕀 b diff --git a/src/splines/Math/Linear.hs b/src/splines/Math/Linear.hs index 9deaf6f..e6a6863 100644 --- a/src/splines/Math/Linear.hs +++ b/src/splines/Math/Linear.hs @@ -96,7 +96,6 @@ instance {-# INCOHERENT #-} Show v => Show ( T v ) where showsPrec p ( T v ) = showParen ( p >= 11 ) $ showString "T" - . showSpace . showsPrec 0 v instance Applicative T where diff --git a/src/splines/Math/Module.hs b/src/splines/Math/Module.hs index f342512..b525286 100644 --- a/src/splines/Math/Module.hs +++ b/src/splines/Math/Module.hs @@ -19,8 +19,6 @@ module Math.Module where -- base -import Control.Applicative - ( liftA2 ) import Control.Monad ( guard ) import Data.Coerce diff --git a/src/splines/Math/Ring.hs b/src/splines/Math/Ring.hs index 1a8a8bb..ae3773a 100644 --- a/src/splines/Math/Ring.hs +++ b/src/splines/Math/Ring.hs @@ -14,8 +14,6 @@ module Math.Ring import Prelude ( Num, Fractional ) import Prelude hiding ( Num(..), Fractional(..), Floating(..) ) import qualified Prelude -import Control.Applicative - ( liftA2 ) import Data.Coerce ( coerce ) import Data.Monoid @@ -64,10 +62,11 @@ class Ring a => Field a where class Field a => Floating a where sqrt :: a -> a -class Field a => Transcendental a where +class Floating a => Transcendental a where pi :: a cos :: a -> a sin :: a -> a + atan :: a -> a -------------------------------------------------------------------------------- @@ -136,6 +135,7 @@ instance Prelude.Floating a => Transcendental ( ViaPrelude a ) where pi = coerce $ Prelude.pi @a sin = coerce $ Prelude.sin @a cos = coerce $ Prelude.cos @a + atan = coerce $ Prelude.atan @a --------------------------------------------------------------------------------