diff --git a/brush-strokes/brush-strokes.cabal b/brush-strokes/brush-strokes.cabal index 36703eb..5beaaf9 100644 --- a/brush-strokes/brush-strokes.cabal +++ b/brush-strokes/brush-strokes.cabal @@ -12,10 +12,15 @@ description: extra-source-files: cbits/**/*.{cpp, hpp} +flag use-simd + description: Use SIMD instructions to implement interval arithmetic. + default: True + manual: True + flag use-fma description: Use fused-muliply add instructions to implement interval arithmetic. - default: True - manual: False + default: False + manual: True flag asserts description: Enable debugging assertions. @@ -97,6 +102,18 @@ common common base >= 4.19 + if flag(use-simd) + cpp-options: + -DUSE_SIMD + ghc-options: + -mavx2 + build-depends: + base + >= 4.20 + -- 4.21 + , ghc-prim + >= 0.11 + if flag(asserts) cpp-options: -DASSERTS @@ -160,9 +177,17 @@ library , Math.Module.Internal , TH.Utils - if flag(use-fma) + if flag(use-simd) other-modules: - Math.Interval.FMA + Math.Interval.Internal.SIMD + + if flag(use-fma) && !flag(use-simd) + other-modules: + Math.Interval.Internal.FMA + + if !flag(use-fma) && !flag(use-simd) + other-modules: + Math.Interval.Internal.RoundedHW build-depends: template-haskell diff --git a/brush-strokes/cabal.project b/brush-strokes/cabal.project index 269675e..d4e54c1 100644 --- a/brush-strokes/cabal.project +++ b/brush-strokes/cabal.project @@ -22,3 +22,18 @@ source-repository-package type: git location: https://github.com/sheaf/eigen tag: 3bc1d4bc53edb75e2ee3893763361aa0d5c7714a + +-------------- +-- GHC HEAD -- +-------------- + +repository head.hackage.ghc.haskell.org + url: https://ghc.gitlab.haskell.org/head.hackage/ + secure: True + key-threshold: 3 + root-keys: + 7541f32a4ccca4f97aea3b22f5e593ba2c0267546016b992dfadcd2fe944e55d + 26021a13b401500c8eb2761ca95c61f2d625bfef951b939a8124ed12ecf07329 + f76d08be13e9a61a377a85e2fb63f4c5435d40f8feb3e12eb05905edb8cdea89 + +active-repositories: hackage.haskell.org, head.hackage.ghc.haskell.org:override diff --git a/brush-strokes/src/cusps/bench/Bench/Cases.hs b/brush-strokes/src/cusps/bench/Bench/Cases.hs index 3d206f5..8aae46c 100644 --- a/brush-strokes/src/cusps/bench/Bench/Cases.hs +++ b/brush-strokes/src/cusps/bench/Bench/Cases.hs @@ -69,7 +69,7 @@ trickyCusp2BrushStroke = defaultStartBoxes :: [ Int ] -> [ ( Int, [ Box 2 ] ) ] defaultStartBoxes is = - [ ( i, [ 𝕀 ( ℝ2 zero zero ) ( ℝ2 one one ) ] ) | i <- is ] + [ ( i, [ 𝕀ℝ2 ( 𝕀 zero one ) ( 𝕀 zero one ) ] ) | i <- is ] zero, one :: Double zero = 1e-6 diff --git a/brush-strokes/src/cusps/bench/Bench/Types.hs b/brush-strokes/src/cusps/bench/Bench/Types.hs index bc67c07..2971df6 100644 --- a/brush-strokes/src/cusps/bench/Bench/Types.hs +++ b/brush-strokes/src/cusps/bench/Bench/Types.hs @@ -43,7 +43,7 @@ import Math.Root.Isolation data BrushStroke = forall nbParams. ParamsCt nbParams => BrushStroke - { brush :: !( Brush ( ℝ nbParams ) ) + { brush :: !( Brush nbParams ) , stroke :: !( Point nbParams, Curve Open () ( Point nbParams ) ) } @@ -60,18 +60,19 @@ data TestCase = type ParamsCt nbParams = ( Show ( ℝ nbParams ) , HasChainRule Double 2 ( ℝ nbParams ) - , HasChainRule ( 𝕀 Double ) 3 ( 𝕀 ( ℝ nbParams ) ) + , HasChainRule 𝕀 3 ( 𝕀ℝ nbParams ) , Applicative ( D 2 ( ℝ nbParams ) ) , Applicative ( D 3 ( ℝ nbParams ) ) , Traversable ( D 2 ( ℝ nbParams ) ) , Traversable ( D 3 ( ℝ nbParams ) ) , Representable Double ( ℝ nbParams ) + , Representable 𝕀 ( 𝕀ℝ nbParams ) , Module Double ( T ( ℝ nbParams ) ) - , Module ( 𝕀 Double ) ( T ( 𝕀 ( ℝ nbParams ) ) ) + , Module 𝕀 ( T ( 𝕀ℝ nbParams ) ) , Module ( D 2 ( ℝ nbParams ) Double ) ( D 2 ( ℝ nbParams ) ( ℝ 2 ) ) - , Module ( D 3 ( ℝ nbParams ) ( 𝕀 Double ) ) ( D 3 ( ℝ nbParams ) ( 𝕀 ( ℝ 2 ) ) ) + , Module ( D 3 ( ℝ nbParams ) 𝕀 ) ( D 3 ( ℝ nbParams ) ( 𝕀ℝ 2 ) ) , Transcendental ( D 2 ( ℝ nbParams ) Double ) - , Transcendental ( D 3 ( ℝ nbParams ) ( 𝕀 Double ) ) + , Transcendental ( D 3 ( ℝ nbParams ) 𝕀 ) ) newtype Params nbParams = Params { getParams :: ( ℝ nbParams ) } @@ -87,70 +88,70 @@ deriving stock instance Show ( ℝ nbParams ) => Show ( Point nbParams ) getStrokeFunctions :: forall nbParams . ParamsCt nbParams - => Brush ( ℝ nbParams ) + => Brush nbParams -- ^ brush shape -> Point nbParams -- ^ start point -> Curve Open () ( Point nbParams ) -- ^ curve points - -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) + -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 ℝ ) , 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) getStrokeFunctions ( Brush brushShape brushShapeI mbRot ) sp0 crv = let usedParams :: C 2 ( ℝ 1 ) ( ℝ nbParams ) path :: C 2 ( ℝ 1 ) ( ℝ 2 ) ( path, usedParams ) = - pathAndUsedParams @2 @() coerce id ( getParams . pointParams ) + pathAndUsedParams @2 @ℝ coerce id id ( getParams . pointParams ) sp0 crv usedParamsI :: C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ nbParams ) pathI :: C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 ) ( pathI, usedParamsI ) = - pathAndUsedParams @3 @𝕀 coerce singleton ( getParams . pointParams ) + pathAndUsedParams @3 @𝕀 coerce point point ( getParams . pointParams ) sp0 crv - in ( brushStrokeData @2 @( ℝ nbParams ) coerce coerce + in ( brushStrokeData @2 @nbParams coerce coerce path usedParams brushShape mbRot - , brushStrokeData @3 @( ℝ nbParams ) coerce coerce + , brushStrokeData @3 @nbParams coerce coerce pathI usedParamsI brushShapeI - ( fmap nonDecreasing mbRot ) + ( fmap ( \ rot -> un𝕀ℝ1 . nonDecreasing ( ℝ1 . rot ) ) mbRot ) ) {-# INLINEABLE getStrokeFunctions #-} {-# SPECIALISE getStrokeFunctions - :: Brush ( ℝ 1 ) + :: Brush 1 -> Point 1 -> Curve Open () ( Point 1 ) - -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) + -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 ℝ ) , 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) #-} {-# SPECIALISE getStrokeFunctions - :: Brush ( ℝ 2 ) + :: Brush 2 -> Point 2 -> Curve Open () ( Point 2 ) - -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) + -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 ℝ ) , 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) #-} {-# SPECIALISE getStrokeFunctions - :: Brush ( ℝ 3 ) + :: Brush 3 -> Point 3 -> Curve Open () ( Point 3 ) - -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) + -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 ℝ ) , 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) #-} {-# SPECIALISE getStrokeFunctions - :: Brush ( ℝ 4 ) + :: Brush 4 -> Point 4 -> Curve Open () ( Point 4 ) - -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) + -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 ℝ ) , 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) #-} brushStrokeFunctions :: BrushStroke - -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) + -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 ℝ ) , 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) -brushStrokeFunctions ( BrushStroke { stroke = ( sp0, crv ), brush = brush :: Brush ( ℝ nbParams ) } ) +brushStrokeFunctions ( BrushStroke { stroke = ( sp0, crv ), brush = brush :: Brush nbParams } ) -- Trick for triggering specialisation... TODO improve this. | Just Refl <- sameNat @nbParams @1 Proxy Proxy = getStrokeFunctions @1 brush sp0 crv diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs index 4742526..d9fd97d 100644 --- a/brush-strokes/src/cusps/bench/Main.hs +++ b/brush-strokes/src/cusps/bench/Main.hs @@ -167,10 +167,8 @@ instance Semigroup BigBoxes where = j1 | otherwise = j2 - widthX ( 𝕀 ( ℝ2 x_lo _y_lo ) ( ℝ2 x_hi _y_hi ) ) - = x_hi - x_lo - widthY ( 𝕀 ( ℝ2 _x_lo y_lo ) ( ℝ2 _x_hi y_hi ) ) - = y_hi - y_lo + widthX ( 𝕀ℝ2 x _y ) = width x + widthY ( 𝕀ℝ2 _x y ) = width y instance Monoid BigBoxes where mempty = BigBoxes Nothing Nothing Nothing @@ -206,8 +204,8 @@ benchGroups = ( defaultStartBoxes [ 0 .. 3 ] ) | let Ξ΅_bis = 5e-3 -- <- [ 1e-6, 5e-6, 1e-5, 5e-5, 1e-4, 5e-4, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 0.1, 0.2, 0.3 ] , (newtMeth, newtOpts) - <- [ ( "LP", \ hist -> NewtonLP ) - , ( "GS_Complete", defaultNewtonOptions @2 @3 ) + <- [ --( "LP", \ hist -> NewtonLP ) + ( "GS_Complete", defaultNewtonOptions @2 @3 ) , ( "GS_Partial", \hist -> NewtonGaussSeidel $ ( defaultGaussSeidelOptions @2 @3 hist ) { gsUpdate = GS_Partial } ) ] , let opts = @@ -230,15 +228,15 @@ benchGroups = getR1 (ℝ1 u) = u eval - :: ( I i ( ℝ 1 ) -> Seq ( I i ( ℝ 1 ) -> StrokeDatum k i ) ) - -> ( I i ( ℝ 1 ), Int, I i ( ℝ 1 ) ) + :: ( I i 1 -> Seq ( I i 1 -> StrokeDatum k i ) ) + -> ( I i 1, Int, I i 1 ) -> StrokeDatum k i eval f ( t, i, s ) = ( f t `Seq.index` i ) s mkVal :: Double -> Int -> Double -> ( ℝ 1, Int, ℝ 1 ) mkVal t i s = ( ℝ1 t, i, ℝ1 s ) -mkI :: ( Double, Double ) -> 𝕀 Double +mkI :: ( Double, Double ) -> 𝕀 mkI ( lo, hi ) = 𝕀 lo hi mkBox :: ( Double, Double ) -> ( Double, Double ) -> Box 2 @@ -255,7 +253,7 @@ potentialCusp && vx_min <= 0 && vx_max >= 0 && vy_min <= 0 && vy_max >= 0 -dEdsdcdt :: StrokeDatum k i -> D ( k - 2 ) ( I i ( ℝ 2 ) ) ( T ( I i ( ℝ 2 ) ) ) +dEdsdcdt :: StrokeDatum k i -> D ( k - 2 ) ( I i 2 ) ( T ( I i 2 ) ) dEdsdcdt ( StrokeDatum { 𝛿E𝛿sdcdt = v } ) = v width :: 𝕀ℝ 1 -> Double diff --git a/brush-strokes/src/cusps/inspect/Main.hs b/brush-strokes/src/cusps/inspect/Main.hs index 485330f..5e0f265 100644 --- a/brush-strokes/src/cusps/inspect/Main.hs +++ b/brush-strokes/src/cusps/inspect/Main.hs @@ -75,42 +75,42 @@ data PointData params outlineFunction :: 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 ) + , 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 ) , HasType ( ℝ 2 ) ( PointData brushParams ) - , Cross ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) + , 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 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 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 ) + => ( 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 ) ) ) + -> 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 ) + -> ( I i 1 -> I i 1 -> StrokeDatum 3 i ) 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 :: C 3 ( I i 1 ) ( I i 2 ) + params :: C 3 ( I i 1 ) ( I i brushParams ) (path, params) = pathAndUsedParams @3 @i co2 single brushParams sp0 crv - strokeData :: I i ( ℝ 1 ) -> I i ( ℝ 1 ) -> StrokeDatum 3 i + strokeData :: I i 1 -> I i 1 -> StrokeDatum 3 i strokeData = fmap ( `Seq.index` brush_index ) $ brushStrokeData @3 @brushParams co1 co2 path params brush diff --git a/brush-strokes/src/cusps/inspect/Math/Interval/Abstract.hs b/brush-strokes/src/cusps/inspect/Math/Interval/Abstract.hs index 2d4baa6..4be86c6 100644 --- a/brush-strokes/src/cusps/inspect/Math/Interval/Abstract.hs +++ b/brush-strokes/src/cusps/inspect/Math/Interval/Abstract.hs @@ -421,7 +421,7 @@ normalise vars expr = case expr of where expr' = Tabulate v' v' :: Vec n ( AI Double ) - mxs :: Vec n ( Maybe ( 𝕀 Double ) ) + mxs :: Vec n ( Maybe 𝕀 ) (v', mxs) = unzip $ fmap ( normalise vars ) v Index iv i | Just r <- index_maybe iv' i @@ -510,9 +510,9 @@ instance Transcendental ( AI Double ) where sin = Sin atan = Atan -class ( Eq v, Show v, Module ( 𝕀 Double ) ( T ( 𝕀 v ) ), Group ( T ( 𝕀 v ) ) ) +class ( Eq v, Show v, Module 𝕀 ( T ( 𝕀 v ) ), Group ( T ( 𝕀 v ) ) ) => IsVec v -instance ( Eq v, Show v, Module ( 𝕀 Double ) ( T ( 𝕀 v ) ), Group ( T ( 𝕀 v ) ) ) +instance ( Eq v, Show v, Module 𝕀 ( T ( 𝕀 v ) ), Group ( T ( 𝕀 v ) ) ) => IsVec v instance IsVec v => Module ( AI Double ) ( T ( AI v ) ) where diff --git a/brush-strokes/src/lib/Calligraphy/Brushes.hs b/brush-strokes/src/lib/Calligraphy/Brushes.hs index 402eaae..151ab70 100644 --- a/brush-strokes/src/lib/Calligraphy/Brushes.hs +++ b/brush-strokes/src/lib/Calligraphy/Brushes.hs @@ -17,7 +17,7 @@ import Prelude hiding ( Num(..), Floating(..), (^), (/), fromInteger, fromRational ) import Data.Kind - ( Type ) + ( Type, Constraint ) import GHC.Exts ( Proxy#, proxy# ) import GHC.TypeNats @@ -33,7 +33,7 @@ import Math.Bezier.Spline import Math.Differentiable ( I ) import Math.Interval - ( 𝕀, singleton ) + ( 𝕀, singleton, point ) import Math.Linear import Math.Module ( Module((^+^), (*^)) ) @@ -42,32 +42,32 @@ import Math.Ring -------------------------------------------------------------------------------- -- | The shape of a brush (before applying any rotation). -type BrushFn :: forall {kd}. kd -> Nat -> Type -> Type -type BrushFn i k brushParams - = C k ( I i brushParams ) - ( Spline Closed () ( I i ( ℝ 2 ) ) ) +type BrushFn :: forall {kd}. kd -> Nat -> Nat -> Type +type BrushFn i k nbBrushParams + = C k ( I i nbBrushParams ) + ( Spline Closed () ( I i 2 ) ) -- | A brush, described as a base shape + an optional rotation. -data Brush brushParams +data Brush nbBrushParams = Brush { -- | Base brush shape, before applying any rotation (if any). - brushBaseShape :: BrushFn () 2 brushParams + brushBaseShape :: BrushFn ℝ 2 nbBrushParams -- | Base brush shape, before applying any rotation (if any). - , brushBaseShapeI :: BrushFn 𝕀 3 brushParams + , brushBaseShapeI :: BrushFn 𝕀 3 nbBrushParams -- | Optional rotation angle function -- (assumed to be a linear function). - , mbRotation :: Maybe ( brushParams -> Double ) + , mbRotation :: Maybe ( ℝ nbBrushParams -> Double ) } -------------------------------------------------------------------------------- -- Brushes -- Some convenience type synonyms for brush types... a bit horrible -type ParamsCt rec = ( ParamsICt 2 () rec, ParamsICt 3 𝕀 rec ) +type ParamsICt :: Nat -> k -> Nat -> Constraint type ParamsICt k i rec = ( Module ( D k ( I i rec ) ( I i Double ) ) - ( D k ( I i rec ) ( I i ( ℝ 2 ) ) ) + ( D k ( I i rec ) ( I i 2 ) ) , Module ( I i Double ) ( T ( I i Double ) ) , HasChainRule ( I i Double ) k ( I i rec ) , Representable ( I i Double ) ( I i rec ) @@ -75,29 +75,29 @@ type ParamsICt k i rec = ) {-# INLINEABLE circleBrush #-} -circleBrush :: ( 1 <= RepDim params, ParamsCt params ) => Brush params +circleBrush :: Brush 1 circleBrush = Brush - { brushBaseShape = circleBrushFn @() @2 proxy# id - , brushBaseShapeI = circleBrushFn @𝕀 @3 proxy# singleton + { brushBaseShape = circleBrushFn @ℝ @2 @1 proxy# id id + , brushBaseShapeI = circleBrushFn @𝕀 @3 @1 proxy# singleton point , mbRotation = Nothing } {-# INLINEABLE ellipseBrush #-} -ellipseBrush :: ( 3 <= RepDim params, ParamsCt params ) => Brush params +ellipseBrush :: Brush 3 ellipseBrush = Brush - { brushBaseShape = ellipseBrushFn @() @2 proxy# id - , brushBaseShapeI = ellipseBrushFn @𝕀 @3 proxy# singleton + { brushBaseShape = ellipseBrushFn @ℝ @2 @3 proxy# id id + , brushBaseShapeI = ellipseBrushFn @𝕀 @3 @3 proxy# singleton point , mbRotation = Just ( `index` ( Fin 3 ) ) } {-# INLINEABLE tearDropBrush #-} -tearDropBrush :: ( 3 <= RepDim params, ParamsCt params ) => Brush params +tearDropBrush :: Brush 3 tearDropBrush = Brush - { brushBaseShape = tearDropBrushFn @() @2 proxy# id - , brushBaseShapeI = tearDropBrushFn @𝕀 @3 proxy# singleton + { brushBaseShape = tearDropBrushFn @ℝ @2 @3 proxy# id id + , brushBaseShapeI = tearDropBrushFn @𝕀 @3 @3 proxy# singleton point , mbRotation = Just ( `index` ( Fin 3 ) ) } @@ -127,54 +127,59 @@ circleSpline p = sequenceA $ Bezier3To ( p ΞΊ -1 ) ( p 1 -ΞΊ ) BackToStart () {-# INLINE circleSpline #-} -circleBrushFn :: forall {t} (i :: t) k rec - . ( 1 <= RepDim ( I i rec ) - , ParamsICt k i rec +circleBrushFn :: forall {t} (i :: t) k nbParams + . ( 1 <= nbParams + , ParamsICt k i nbParams ) => Proxy# i - -> ( forall a. a -> I i a ) - -> C k ( I i rec ) ( Spline 'Closed () ( I i ( ℝ 2 ) ) ) -circleBrushFn _ mkI = + -> ( Double -> I i Double ) + -> ( I ℝ 2 -> I i 2 ) + -> C k ( I i nbParams ) ( Spline 'Closed () ( I i 2 ) ) +circleBrushFn _ mkI1 mkI2 = D \ params -> - let r :: D k ( I i rec ) ( I i Double ) + let r :: D k ( I i nbParams ) ( I i Double ) r = runD ( var @_ @k $ Fin 1 ) params - mkPt :: Double -> Double -> D k ( I i rec ) ( I i ( ℝ 2 ) ) + mkPt :: Double -> Double -> D k ( I i nbParams ) ( I i 2 ) mkPt x y = ( r `scaledBy` x ) *^ e_x ^+^ ( r `scaledBy` y ) *^ e_y in circleSpline mkPt where - e_x, e_y :: D k ( I i rec ) ( I i ( ℝ 2 ) ) - e_x = pure $ mkI $ ℝ2 1 0 - e_y = pure $ mkI $ ℝ2 0 1 + e_x, e_y :: D k ( I i nbParams ) ( I i 2 ) + e_x = pure $ mkI2 $ ℝ2 1 0 + e_y = pure $ mkI2 $ ℝ2 0 1 - scaledBy d x = fmap ( mkI x * ) d + scaledBy :: D k ( I i nbParams ) ( I i Double ) -> Double -> D k ( I i nbParams ) ( I i Double ) + scaledBy d x = fmap ( mkI1 x * ) d {-# INLINEABLE circleBrushFn #-} -ellipseBrushFn :: forall {t} (i :: t) k rec - . ( 3 <= RepDim ( I i rec ) - , ParamsICt k i rec +ellipseBrushFn :: forall {t} (i :: t) k nbParams + . ( 3 <= nbParams + , ParamsICt k i nbParams ) => Proxy# i - -> ( forall a. a -> I i a ) - -> C k ( I i rec ) ( Spline 'Closed () ( I i ( ℝ 2 ) ) ) -ellipseBrushFn _ mkI = + -> ( Double -> I i Double ) + -> ( I ℝ 2 -> I i 2 ) + + -> C k ( I i nbParams ) ( Spline 'Closed () ( I i 2 ) ) +ellipseBrushFn _ mkI1 mkI2 = D \ params -> - let a, b :: D k ( I i rec ) ( I i Double ) + let a, b :: D k ( I i nbParams ) ( I i Double ) a = runD ( var @_ @k $ Fin 1 ) params b = runD ( var @_ @k $ Fin 2 ) params - mkPt :: Double -> Double -> D k ( I i rec ) ( I i ( ℝ 2 ) ) + mkPt :: Double -> Double -> D k ( I i nbParams ) ( I i 2 ) mkPt x y = let !x' = a `scaledBy` x !y' = b `scaledBy` y in x' *^ e_x ^+^ y' *^ e_y in circleSpline mkPt where - e_x, e_y :: D k ( I i rec ) ( I i ( ℝ 2 ) ) - e_x = pure $ mkI $ ℝ2 1 0 - e_y = pure $ mkI $ ℝ2 0 1 + e_x, e_y :: D k ( I i nbParams ) ( I i 2 ) + e_x = pure $ mkI2 $ ℝ2 1 0 + e_y = pure $ mkI2 $ ℝ2 0 1 - scaledBy d x = fmap ( mkI x * ) d + scaledBy :: D k ( I i nbParams ) ( I i Double ) -> Double -> D k ( I i nbParams ) ( I i Double ) + scaledBy d x = fmap ( mkI1 x * ) d {-# INLINEABLE ellipseBrushFn #-} -------------------------------------------------------------------------------- @@ -197,19 +202,20 @@ tearHeight = 3 * sqrt 3 / 8 sqrt3_over_2 :: Double sqrt3_over_2 = 0.5 * sqrt 3 -tearDropBrushFn :: forall {t} (i :: t) k rec - . ( 3 <= RepDim ( I i rec ) - , ParamsICt k i rec - ) +tearDropBrushFn :: forall {t} (i :: t) k nbParams + . ( 3 <= nbParams + , ParamsICt k i nbParams + ) => Proxy# i - -> ( forall a. a -> I i a ) - -> C k ( I i rec ) ( Spline 'Closed () ( I i ( ℝ 2 ) ) ) -tearDropBrushFn _ mkI = + -> ( Double -> I i Double ) + -> ( I ℝ 2 -> I i 2 ) + -> C k ( I i nbParams ) ( Spline 'Closed () ( I i 2 ) ) +tearDropBrushFn _ mkI1 mkI2 = D \ params -> - let w, h :: D k ( I i rec ) ( I i Double ) + let w, h :: D k ( I i nbParams ) ( I i Double ) w = runD ( var @_ @k ( Fin 1 ) ) params h = runD ( var @_ @k ( Fin 2 ) ) params - mkPt :: Double -> Double -> D k ( I i rec ) ( I i ( ℝ 2 ) ) + mkPt :: Double -> Double -> D k ( I i nbParams ) ( I i 2 ) mkPt x y -- 1. translate the teardrop so that the centre of mass is at the origin -- 2. scale the teardrop so that it has the requested width/height @@ -226,9 +232,10 @@ tearDropBrushFn _ mkI = ( mkPt -0.5 sqrt3_over_2 ) BackToStart () } where - e_x, e_y :: D k ( I i rec ) ( I i ( ℝ 2 ) ) - e_x = pure $ mkI $ ℝ2 1 0 - e_y = pure $ mkI $ ℝ2 0 1 + e_x, e_y :: D k ( I i nbParams ) ( I i 2 ) + e_x = pure $ mkI2 $ ℝ2 1 0 + e_y = pure $ mkI2 $ ℝ2 0 1 - scaledBy d x = fmap ( mkI x * ) d + scaledBy :: D k ( I i nbParams ) ( I i Double ) -> Double -> D k ( I i nbParams ) ( I i Double ) + scaledBy d x = fmap ( mkI1 x * ) d {-# INLINEABLE tearDropBrushFn #-} diff --git a/brush-strokes/src/lib/Math/Algebra/Dual.hs b/brush-strokes/src/lib/Math/Algebra/Dual.hs index 470538f..6150f44 100644 --- a/brush-strokes/src/lib/Math/Algebra/Dual.hs +++ b/brush-strokes/src/lib/Math/Algebra/Dual.hs @@ -289,7 +289,7 @@ instance Ring r => Ring ( D2𝔸1 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Ring ( D2𝔸1 Double ) #-} - {-# SPECIALISE instance Ring ( D2𝔸1 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Ring ( D2𝔸1 𝕀 ) #-} instance Ring r => Ring ( D2𝔸2 r ) where !dr1 * !dr2 = @@ -314,7 +314,7 @@ instance Ring r => Ring ( D2𝔸2 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Ring ( D2𝔸2 Double ) #-} - {-# SPECIALISE instance Ring ( D2𝔸2 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Ring ( D2𝔸2 𝕀 ) #-} instance Ring r => Ring ( D2𝔸3 r ) where !dr1 * !dr2 = @@ -339,7 +339,7 @@ instance Ring r => Ring ( D2𝔸3 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Ring ( D2𝔸3 Double ) #-} - {-# SPECIALISE instance Ring ( D2𝔸3 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Ring ( D2𝔸3 𝕀 ) #-} instance Ring r => Ring ( D2𝔸4 r ) where !dr1 * !dr2 = @@ -364,7 +364,7 @@ instance Ring r => Ring ( D2𝔸4 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Ring ( D2𝔸4 Double ) #-} - {-# SPECIALISE instance Ring ( D2𝔸4 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Ring ( D2𝔸4 𝕀 ) #-} instance Ring r => Ring ( D3𝔸1 r ) where !dr1 * !dr2 = @@ -389,7 +389,7 @@ instance Ring r => Ring ( D3𝔸1 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Ring ( D3𝔸1 Double ) #-} - {-# SPECIALISE instance Ring ( D3𝔸1 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Ring ( D3𝔸1 𝕀 ) #-} instance Ring r => Ring ( D3𝔸2 r ) where !dr1 * !dr2 = @@ -414,7 +414,7 @@ instance Ring r => Ring ( D3𝔸2 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Ring ( D3𝔸2 Double ) #-} - {-# SPECIALISE instance Ring ( D3𝔸2 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Ring ( D3𝔸2 𝕀 ) #-} instance Ring r => Ring ( D3𝔸3 r ) where !dr1 * !dr2 = @@ -439,7 +439,7 @@ instance Ring r => Ring ( D3𝔸3 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Ring ( D3𝔸3 Double ) #-} - {-# SPECIALISE instance Ring ( D3𝔸3 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Ring ( D3𝔸3 𝕀 ) #-} instance Ring r => Ring ( D3𝔸4 r ) where !dr1 * !dr2 = @@ -464,7 +464,7 @@ instance Ring r => Ring ( D3𝔸4 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Ring ( D3𝔸4 Double ) #-} - {-# SPECIALISE instance Ring ( D3𝔸4 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Ring ( D3𝔸4 𝕀 ) #-} -------------------------------------------------------------------------------- -- Field, floating & transcendental instances @@ -503,7 +503,7 @@ instance Field r => Field ( D2𝔸1 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Field ( D2𝔸1 Double ) #-} - {-# SPECIALISE instance Field ( D2𝔸1 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Field ( D2𝔸1 𝕀 ) #-} instance Field r => Field ( D2𝔸2 r ) where fromRational q = konst @Double @2 @( ℝ 2 ) ( fromRational q ) recip df = @@ -516,7 +516,7 @@ instance Field r => Field ( D2𝔸2 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Field ( D2𝔸2 Double ) #-} - {-# SPECIALISE instance Field ( D2𝔸2 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Field ( D2𝔸2 𝕀 ) #-} instance Field r => Field ( D2𝔸3 r ) where fromRational q = konst @Double @2 @( ℝ 3 ) ( fromRational q ) recip df = @@ -529,7 +529,7 @@ instance Field r => Field ( D2𝔸3 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Field ( D2𝔸3 Double ) #-} - {-# SPECIALISE instance Field ( D2𝔸3 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Field ( D2𝔸3 𝕀 ) #-} instance Field r => Field ( D2𝔸4 r ) where fromRational q = konst @Double @2 @( ℝ 4 ) ( fromRational q ) recip df = @@ -542,7 +542,7 @@ instance Field r => Field ( D2𝔸4 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Field ( D2𝔸4 Double ) #-} - {-# SPECIALISE instance Field ( D2𝔸4 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Field ( D2𝔸4 𝕀 ) #-} instance Field r => Field ( D3𝔸1 r ) where fromRational q = konst @Double @3 @( ℝ 1 ) ( fromRational q ) recip df = @@ -555,7 +555,7 @@ instance Field r => Field ( D3𝔸1 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Field ( D3𝔸1 Double ) #-} - {-# SPECIALISE instance Field ( D3𝔸1 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Field ( D3𝔸1 𝕀 ) #-} instance Field r => Field ( D3𝔸2 r ) where fromRational q = konst @Double @3 @( ℝ 2 ) ( fromRational q ) recip df = @@ -568,7 +568,7 @@ instance Field r => Field ( D3𝔸2 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Field ( D3𝔸2 Double ) #-} - {-# SPECIALISE instance Field ( D3𝔸2 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Field ( D3𝔸2 𝕀 ) #-} instance Field r => Field ( D3𝔸3 r ) where fromRational q = konst @Double @3 @( ℝ 3 ) ( fromRational q ) recip df = @@ -581,7 +581,7 @@ instance Field r => Field ( D3𝔸3 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Field ( D3𝔸3 Double ) #-} - {-# SPECIALISE instance Field ( D3𝔸3 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Field ( D3𝔸3 𝕀 ) #-} instance Field r => Field ( D3𝔸4 r ) where fromRational q = konst @Double @3 @( ℝ 4 ) ( fromRational q ) recip df = @@ -594,7 +594,7 @@ instance Field r => Field ( D3𝔸4 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Field ( D3𝔸4 Double ) #-} - {-# SPECIALISE instance Field ( D3𝔸4 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Field ( D3𝔸4 𝕀 ) #-} d1sqrt :: Floating r => r -> D1𝔸1 r d1sqrt x = @@ -631,7 +631,7 @@ instance Floating r => Floating ( D2𝔸1 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Floating ( D2𝔸1 Double ) #-} - {-# SPECIALISE instance Floating ( D2𝔸1 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Floating ( D2𝔸1 𝕀 ) #-} instance Floating r => Floating ( D2𝔸2 r ) where sqrt df = let @@ -643,7 +643,7 @@ instance Floating r => Floating ( D2𝔸2 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Floating ( D2𝔸2 Double ) #-} - {-# SPECIALISE instance Floating ( D2𝔸2 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Floating ( D2𝔸2 𝕀 ) #-} instance Floating r => Floating ( D2𝔸3 r ) where sqrt df = let @@ -655,7 +655,7 @@ instance Floating r => Floating ( D2𝔸3 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Floating ( D2𝔸3 Double ) #-} - {-# SPECIALISE instance Floating ( D2𝔸3 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Floating ( D2𝔸3 𝕀 ) #-} instance Floating r => Floating ( D2𝔸4 r ) where sqrt df = let @@ -667,7 +667,7 @@ instance Floating r => Floating ( D2𝔸4 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Floating ( D2𝔸4 Double ) #-} - {-# SPECIALISE instance Floating ( D2𝔸4 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Floating ( D2𝔸4 𝕀 ) #-} instance Floating r => Floating ( D3𝔸1 r ) where sqrt df = let @@ -679,7 +679,7 @@ instance Floating r => Floating ( D3𝔸1 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Floating ( D3𝔸1 Double ) #-} - {-# SPECIALISE instance Floating ( D3𝔸1 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Floating ( D3𝔸1 𝕀 ) #-} instance Floating r => Floating ( D3𝔸2 r ) where sqrt df = let @@ -691,7 +691,7 @@ instance Floating r => Floating ( D3𝔸2 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Floating ( D3𝔸2 Double ) #-} - {-# SPECIALISE instance Floating ( D3𝔸2 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Floating ( D3𝔸2 𝕀 ) #-} instance Floating r => Floating ( D3𝔸3 r ) where sqrt df = let @@ -703,7 +703,7 @@ instance Floating r => Floating ( D3𝔸3 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Floating ( D3𝔸3 Double ) #-} - {-# SPECIALISE instance Floating ( D3𝔸3 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Floating ( D3𝔸3 𝕀 ) #-} instance Floating r => Floating ( D3𝔸4 r ) where sqrt df = let @@ -715,7 +715,7 @@ instance Floating r => Floating ( D3𝔸4 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Floating ( D3𝔸4 Double ) #-} - {-# SPECIALISE instance Floating ( D3𝔸4 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Floating ( D3𝔸4 𝕀 ) #-} d1sin, d1cos, d1atan :: Transcendental r => r -> D1𝔸1 r d1sin x = @@ -809,7 +809,7 @@ instance Transcendental r => Transcendental ( D2𝔸1 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Transcendental ( D2𝔸1 Double ) #-} - {-# SPECIALISE instance Transcendental ( D2𝔸1 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Transcendental ( D2𝔸1 𝕀 ) #-} instance Transcendental r => Transcendental ( D2𝔸2 r ) where pi = konst @Double @2 @( ℝ 2 ) pi sin df = @@ -840,7 +840,7 @@ instance Transcendental r => Transcendental ( D2𝔸2 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Transcendental ( D2𝔸2 Double ) #-} - {-# SPECIALISE instance Transcendental ( D2𝔸2 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Transcendental ( D2𝔸2 𝕀 ) #-} instance Transcendental r => Transcendental ( D2𝔸3 r ) where pi = konst @Double @2 @( ℝ 3 ) pi sin df = @@ -871,7 +871,7 @@ instance Transcendental r => Transcendental ( D2𝔸3 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Transcendental ( D2𝔸3 Double ) #-} - {-# SPECIALISE instance Transcendental ( D2𝔸3 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Transcendental ( D2𝔸3 𝕀 ) #-} instance Transcendental r => Transcendental ( D2𝔸4 r ) where pi = konst @Double @2 @( ℝ 4 ) pi @@ -903,7 +903,7 @@ instance Transcendental r => Transcendental ( D2𝔸4 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Transcendental ( D2𝔸4 Double ) #-} - {-# SPECIALISE instance Transcendental ( D2𝔸4 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Transcendental ( D2𝔸4 𝕀 ) #-} instance Transcendental r => Transcendental ( D3𝔸1 r ) where pi = konst @Double @3 @( ℝ 1 ) pi @@ -935,7 +935,7 @@ instance Transcendental r => Transcendental ( D3𝔸1 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Transcendental ( D3𝔸1 Double ) #-} - {-# SPECIALISE instance Transcendental ( D3𝔸1 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Transcendental ( D3𝔸1 𝕀 ) #-} instance Transcendental r => Transcendental ( D3𝔸2 r ) where pi = konst @Double @3 @( ℝ 2 ) pi @@ -967,7 +967,7 @@ instance Transcendental r => Transcendental ( D3𝔸2 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Transcendental ( D3𝔸2 Double ) #-} - {-# SPECIALISE instance Transcendental ( D3𝔸2 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Transcendental ( D3𝔸2 𝕀 ) #-} instance Transcendental r => Transcendental ( D3𝔸3 r ) where pi = konst @Double @3 @( ℝ 3 ) pi @@ -999,7 +999,7 @@ instance Transcendental r => Transcendental ( D3𝔸3 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Transcendental ( D3𝔸3 Double ) #-} - {-# SPECIALISE instance Transcendental ( D3𝔸3 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Transcendental ( D3𝔸3 𝕀 ) #-} instance Transcendental r => Transcendental ( D3𝔸4 r ) where pi = konst @Double @3 @( ℝ 4 ) pi @@ -1031,7 +1031,7 @@ instance Transcendental r => Transcendental ( D3𝔸4 r ) where in $$( chainRuleN1Q [|| fromInt ||] [|| add ||] [|| times ||] [|| pow ||] [|| df ||] [|| dg ||] ) {-# SPECIALISE instance Transcendental ( D3𝔸4 Double ) #-} - {-# SPECIALISE instance Transcendental ( D3𝔸4 ( 𝕀 Double ) ) #-} + {-# SPECIALISE instance Transcendental ( D3𝔸4 𝕀 ) #-} -------------------------------------------------------------------------------- -- HasChainRule instances. diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke.hs b/brush-strokes/src/lib/Math/Bezier/Stroke.hs index acfc766..a7e0f20 100644 --- a/brush-strokes/src/lib/Math/Bezier/Stroke.hs +++ b/brush-strokes/src/lib/Math/Bezier/Stroke.hs @@ -62,7 +62,7 @@ import GHC.STRef import GHC.Generics ( Generic, Generic1, Generically(..) ) import GHC.TypeNats - ( type (-) ) + ( Nat, type (-) ) -- acts import Data.Act @@ -229,34 +229,35 @@ data OutlineInfo = type N = 2 computeStrokeOutline :: - forall ( clo :: SplineType ) usedParams brushParams crvData ptData s + forall ( clo :: SplineType ) ( nbUsedParams :: Nat ) ( nbBrushParams :: Nat ) crvData ptData s . ( KnownSplineType clo , HasType ( ℝ 2 ) ptData , HasType ( CachedStroke s ) crvData , NFData ptData, NFData crvData -- Differentiability. - , Interpolatable Double usedParams - , Interpolatable ( 𝕀 Double ) ( 𝕀 usedParams ) - , DiffInterp 2 () brushParams - , DiffInterp 3 𝕀 brushParams - , HasChainRule Double 2 usedParams - , HasChainRule ( 𝕀 Double ) 3 ( 𝕀 usedParams ) - , HasChainRule Double 2 brushParams - , HasChainRule ( 𝕀 Double ) 3 ( 𝕀 brushParams ) - , Traversable ( D 2 brushParams ) - , Representable Double usedParams + , Interpolatable Double ( ℝ nbUsedParams ) + , DiffInterp 2 ℝ nbBrushParams + , DiffInterp 3 𝕀 nbBrushParams + , HasChainRule Double 2 ( ℝ nbUsedParams ) + , HasChainRule 𝕀 3 ( 𝕀ℝ nbUsedParams ) + , HasChainRule Double 2 ( ℝ nbBrushParams ) + , HasChainRule 𝕀 3 ( 𝕀ℝ nbBrushParams ) + , Traversable ( D 2 ( ℝ nbBrushParams ) ) + , Representable Double ( ℝ nbUsedParams ) + , Representable 𝕀 ( 𝕀ℝ nbUsedParams ) + , Module 𝕀 (T (𝕀ℝ nbUsedParams)) -- Debugging. - , Show ptData, Show brushParams + , Show ptData, Show ( ℝ nbBrushParams ) ) => RootSolvingAlgorithm -> Maybe ( RootIsolationOptions N 3 ) -> FitParameters - -> ( ptData -> usedParams ) - -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing - -> Brush brushParams + -> ( ptData -> ℝ nbUsedParams ) + -> ( ℝ nbUsedParams -> ℝ nbBrushParams ) -- ^ assumed to be linear and non-decreasing + -> Brush nbBrushParams -> Spline clo crvData ptData -> ST s ( Either ( SplinePts Closed ) ( SplinePts Closed, SplinePts Closed ) @@ -514,31 +515,32 @@ computeStrokeOutline rootAlgo mbCuspOptions fitParams ptParams toBrushParams bru -- | Computes the forward and backward stroke outline functions for a single curve. outlineFunction - :: forall usedParams brushParams crvData ptData + :: forall (nbUsedParams :: Nat) (nbBrushParams :: Nat) crvData ptData . ( HasType ( ℝ 2 ) ptData -- Differentiability. - , Interpolatable Double usedParams - , Interpolatable ( 𝕀 Double ) ( 𝕀 usedParams ) - , DiffInterp 2 () brushParams - , DiffInterp 3 𝕀 brushParams - , HasChainRule Double 2 usedParams - , HasChainRule ( 𝕀 Double ) 3 ( 𝕀 usedParams ) - , HasChainRule Double 2 brushParams - , HasChainRule ( 𝕀 Double ) 3 ( 𝕀 brushParams ) - , Traversable ( D 2 brushParams ) + , Interpolatable Double ( ℝ nbUsedParams ) + , DiffInterp 2 ℝ nbBrushParams + , DiffInterp 3 𝕀 nbBrushParams + , HasChainRule Double 2 ( ℝ nbUsedParams ) + , HasChainRule 𝕀 3 ( 𝕀ℝ nbUsedParams ) + , HasChainRule Double 2 ( ℝ nbBrushParams ) + , HasChainRule 𝕀 3 ( 𝕀ℝ nbBrushParams ) + , Traversable ( D 2 ( ℝ nbBrushParams ) ) + , Module 𝕀 (T (𝕀ℝ nbUsedParams)) -- Computing AABBs - , Representable Double usedParams + , Representable Double ( ℝ nbUsedParams ) + , Representable 𝕀 ( 𝕀ℝ nbUsedParams ) -- Debugging. - , Show ptData, Show brushParams + , Show ptData, Show ( ℝ nbBrushParams ) ) => RootSolvingAlgorithm -> Maybe ( RootIsolationOptions N 3 ) - -> ( ptData -> usedParams ) - -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing - -> Brush brushParams + -> ( ptData -> ℝ nbUsedParams ) + -> ( ℝ nbUsedParams -> ℝ nbBrushParams ) -- ^ assumed to be linear and non-decreasing + -> Brush nbBrushParams -> ptData -> Curve Open crvData ptData -> OutlineInfo @@ -546,15 +548,15 @@ outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams ( Brush { brushBaseShape, brushBaseShapeI, mbRotation } ) = \ sp0 crv -> let - usedParams :: C 2 ( ℝ 1 ) usedParams + usedParams :: C 2 ( ℝ 1 ) ( ℝ nbUsedParams ) path :: C 2 ( ℝ 1 ) ( ℝ 2 ) ( path, usedParams ) - = pathAndUsedParams @2 @() coerce id ptParams sp0 crv + = pathAndUsedParams @2 @ℝ @nbUsedParams coerce id id ptParams sp0 crv curves :: ℝ 1 -- t - -> Seq ( ℝ 1 {- s -} -> StrokeDatum 2 () ) + -> Seq ( ℝ 1 {- s -} -> StrokeDatum 2 ℝ ) curves = - brushStrokeData @2 @brushParams + brushStrokeData @2 @nbBrushParams coerce coerce path ( fmap toBrushParams usedParams ) @@ -564,16 +566,16 @@ outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams curvesI :: 𝕀ℝ 1 -- t -> Seq ( 𝕀ℝ 1 {- s -} -> StrokeDatum 3 𝕀 ) curvesI = - brushStrokeData @3 @brushParams + brushStrokeData @3 @nbBrushParams coerce coerce pathI ( fmap ( nonDecreasing toBrushParams ) usedParamsI ) brushBaseShapeI - ( fmap nonDecreasing mbRotation ) + ( fmap ( \ rot -> un𝕀ℝ1 . nonDecreasing ( ℝ1 . rot ) ) mbRotation ) - usedParamsI :: C 3 ( 𝕀ℝ 1 ) ( 𝕀 usedParams ) + usedParamsI :: C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ nbUsedParams ) pathI :: C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 ) - ( pathI, usedParamsI ) = pathAndUsedParams @3 @𝕀 coerce singleton ptParams sp0 crv + ( pathI, usedParamsI ) = pathAndUsedParams @3 @𝕀 @nbUsedParams coerce point point ptParams sp0 crv fwdBwd :: OutlineFn fwdBwd t @@ -597,7 +599,7 @@ outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams in fmap ( unT . rotate cosΞΈ sinΞΈ . T ) in applyRot $ - value @Double @2 @brushParams $ + value @Double @2 @( ℝ nbBrushParams ) $ runD brushBaseShape brushParams ( potentialCusps, definiteCusps ) = @@ -618,39 +620,41 @@ outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams } {-# INLINEABLE outlineFunction #-} -pathAndUsedParams :: forall k i arr crvData ptData usedParams +pathAndUsedParams :: forall k i (nbUsedParams :: Nat) arr crvData ptData . ( HasType ( ℝ 2 ) ptData , HasBΓ©zier k i , arr ~ C 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 ) ) - , Module ( I i Double ) ( T ( I i usedParams ) ) - , Torsor ( T ( I i usedParams ) ) ( I i usedParams ) - , Module Double ( T usedParams ) - , Representable Double usedParams - , Torsor ( T usedParams ) usedParams + , D k ( I i 1 ) ~ D k ( ℝ 1 ) + , Module ( I i Double ) ( T ( I i 2 ) ) + , Torsor ( T ( I i 2 ) ) ( I i 2 ) + , Module ( I i Double ) ( T ( I i nbUsedParams ) ) + , Torsor ( T ( I i nbUsedParams ) ) ( I i nbUsedParams ) + , Module Double ( T ( I ℝ nbUsedParams ) ) + , Torsor ( T ( ( I ℝ nbUsedParams ) ) ) ( I ℝ nbUsedParams ) + , Representable Double ( ℝ nbUsedParams ) + , Representable 𝕀 ( 𝕀ℝ nbUsedParams ) ) - => ( I i ( ℝ 1 ) -> I i Double ) - -> ( forall a. a -> I i a ) - -> ( ptData -> usedParams ) + => ( I i 1 -> I i Double ) + -> ( I ℝ 2 -> I i 2 ) + -> ( I ℝ nbUsedParams -> I i nbUsedParams ) + -> ( ptData -> I ℝ nbUsedParams ) -> ptData -> Curve Open crvData ptData - -> ( I i ( ℝ 1 ) `arr` I i ( ℝ 2 ), I i ( ℝ 1 ) `arr` I i usedParams ) -pathAndUsedParams co toI ptParams sp0 crv = + -> ( I i 1 `arr` I i 2, I i 1 `arr` I i nbUsedParams ) +pathAndUsedParams co toI1 toI2 ptParams sp0 crv = case crv of LineTo { curveEnd = NextPoint sp1 } | let seg = Segment sp0 sp1 - -> ( line @k @i @( ℝ 2 ) co ( fmap ( toI . coords ) seg ) - , line @k @i @usedParams co ( fmap ( toI . ptParams ) seg ) ) + -> ( line @k @i @2 co ( fmap ( toI1 . coords ) seg ) + , line @k @i @nbUsedParams co ( fmap ( toI2 . ptParams ) seg ) ) Bezier2To { controlPoint = sp1, curveEnd = NextPoint sp2 } | let bez2 = Quadratic.Bezier sp0 sp1 sp2 - -> ( bezier2 @k @i @( ℝ 2 ) co ( fmap ( toI . coords ) bez2 ) - , bezier2 @k @i @usedParams co ( fmap ( toI . ptParams ) bez2 ) ) + -> ( bezier2 @k @i @2 co ( fmap ( toI1 . coords ) bez2 ) + , bezier2 @k @i @nbUsedParams co ( fmap ( toI2 . ptParams ) bez2 ) ) Bezier3To { controlPoint1 = sp1, controlPoint2 = sp2, curveEnd = NextPoint sp3 } | let bez3 = Cubic.Bezier sp0 sp1 sp2 sp3 - -> ( bezier3 @k @i @( ℝ 2 ) co ( fmap ( toI . coords ) bez3 ) - , bezier3 @k @i @usedParams co ( fmap ( toI . ptParams ) bez3 ) ) + -> ( bezier3 @k @i @2 co ( fmap ( toI1 . coords ) bez3 ) + , bezier3 @k @i @nbUsedParams co ( fmap ( toI2 . ptParams ) bez3 ) ) {-# INLINEABLE pathAndUsedParams #-} ----------------------------------- @@ -951,12 +955,12 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) splineCurveFns :: forall k i . ( HasBΓ©zier k i - , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) - , Module ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) - , 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 ) ) ) + , D k ( I i 1 ) ~ D k ( ℝ 1 ) + , Module ( I i Double ) ( T ( I i 2 ) ) + , 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 @@ -965,16 +969,16 @@ splineCurveFns co spls . adjustSplineType @Open $ spls where - curveFn :: I i ( ℝ 2 ) - -> Curve Open () ( I i ( ℝ 2 ) ) - -> C k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) + curveFn :: I i 2 + -> Curve Open () ( I i 2 ) + -> C k ( I i 1 ) ( I i 2 ) curveFn p0 = \case LineTo { curveEnd = NextPoint p1 } - -> line @k @i @( ℝ 2 ) co $ Segment p0 p1 + -> line @k @i @2 co $ Segment p0 p1 Bezier2To { controlPoint = p1, curveEnd = NextPoint p2 } - -> bezier2 @k @i @( ℝ 2 ) co $ Quadratic.Bezier p0 p1 p2 + -> bezier2 @k @i @2 co $ Quadratic.Bezier p0 p1 p2 Bezier3To { controlPoint1 = p1, controlPoint2 = p2, curveEnd = NextPoint p3 } - -> bezier3 @k @i @( ℝ 2 ) co $ Cubic.Bezier p0 p1 p2 p3 + -> bezier3 @k @i @2 co $ Cubic.Bezier p0 p1 p2 p3 newtype ZipSeq a = ZipSeq { getZipSeq :: Seq a } deriving stock Functor @@ -983,99 +987,102 @@ instance Applicative ZipSeq where liftA2 f ( ZipSeq xs ) ( ZipSeq ys ) = ZipSeq ( Seq.zipWith f xs ys ) {-# INLINE liftA2 #-} -brushStrokeData :: forall k brushParams i arr +brushStrokeData :: forall {kd} (k :: Nat) (nbBrushParams :: Nat) (i :: kd) arr . ( arr ~ C k , HasBΓ©zier k i, HasEnvelopeEquation k - , Differentiable k i brushParams - , HasChainRule ( I i Double ) k ( I i ( ℝ 1 ) ) + , Differentiable k i nbBrushParams + , HasChainRule ( I i Double ) k ( I i 1 ) , Applicative ( D k ( ℝ 1 ) ) - , D ( k - 2 ) ( I i ( ℝ 2 ) ) ~ D ( k - 2 ) ( ℝ 2 ) - , D ( k - 1 ) ( I i ( ℝ 2 ) ) ~ D ( k - 1 ) ( ℝ 2 ) - , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) - , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) - , D k ( I i ( ℝ 2 ) ) ~ D k ( ℝ 2 ) + , D ( k - 2 ) ( I i 2 ) ~ D ( k - 2 ) ( ℝ 2 ) + , D ( k - 1 ) ( I i 2 ) ~ D ( k - 1 ) ( ℝ 2 ) + , D k ( I i 1 ) ~ D k ( ℝ 1 ) + , D k ( I i 1 ) ~ D k ( ℝ 1 ) + , D k ( I i 2 ) ~ D k ( ℝ 2 ) , Transcendental ( I i Double ) - , Module ( I i Double ) ( T ( I i ( ℝ 1 ) ) ) - , Cross ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) - , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) ) - , Show brushParams - , Representable ( I i Double ) ( I i ( ℝ 2 ) ), RepDim ( I i ( ℝ 2 ) ) ~ 2 + , Module ( I i Double ) ( T ( I i 1 ) ) + , Cross ( I i Double ) ( T ( I i 2 ) ) + , Torsor ( T ( I i 2 ) ) ( I i 2 ) + , Show ( ℝ nbBrushParams ) + , Representable ( I i Double ) ( I i 2 ), RepDim ( I i 2 ) ~ 2 ) - => ( I i Double -> I i ( ℝ 1 ) ) - -> ( I i ( ℝ 1 ) -> I i Double ) - -> ( 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 ) + -> ( I i 1 `arr` I i nbBrushParams ) -- ^ brush parameters - -> ( I i brushParams `arr` Spline Closed () ( I i ( ℝ 2 ) ) ) + -> ( I i nbBrushParams `arr` Spline Closed () ( I i 2 ) ) -- ^ brush from parameters - -> ( Maybe ( I i brushParams -> I i Double ) ) + -> ( Maybe ( I i nbBrushParams -> I i Double ) ) -- ^ rotation parameter - -> ( I i ( ℝ 1 ) -> Seq ( I i ( ℝ 1 ) -> StrokeDatum k i ) ) + -> ( I i 1 -> Seq ( I i 1 -> StrokeDatum k i ) ) brushStrokeData co1 co2 path params brush mbBrushRotation = \ t -> let - dpath_t :: D k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) + dpath_t :: D k ( I i 1 ) ( I i 2 ) !dpath_t = runD path t - dparams_t :: D k ( I i ( ℝ 1 ) ) ( I i brushParams ) + dparams_t :: D k ( I i 1 ) ( I i nbBrushParams ) !dparams_t = runD params t - 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 ) ) ) + dbrush_params :: D k ( I i nbBrushParams ) ( 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 nbBrushParams ) ( I i 1 `arr` I i 2 ) ) !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 :: 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. in fmap ( mkStrokeDatum dpath_t dparams_t ) dbrushes_t where - mkStrokeDatum :: D k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) - -> D k ( I i ( ℝ 1 ) ) ( I i brushParams ) - -> ( I i ( ℝ 1 ) -> D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) ) - -> ( I i ( ℝ 1 ) -> StrokeDatum k i ) + mkStrokeDatum :: D k ( I i 1 ) ( I i 2 ) + -> D k ( I i 1 ) ( I i nbBrushParams ) + -> ( I i 1 -> D k ( I i 2 ) ( I i 2 ) ) + -> ( I i 1 -> StrokeDatum k i ) mkStrokeDatum dpath_t dparams_t dbrush_t s = let dbrush_t_s = dbrush_t s mbRotation = mbBrushRotation <&> \ getTheta -> fmap getTheta dparams_t in envelopeEquation @k @i Proxy co1 dpath_t dbrush_t_s mbRotation {-# INLINEABLE brushStrokeData #-} + +{- {-# SPECIALISE brushStrokeData - :: ( I 𝕀 Double -> I 𝕀 ( ℝ 1 ) ) - -> ( I 𝕀 ( ℝ 1 ) -> I 𝕀 Double ) - -> ( C 3 ( I 𝕀 ( ℝ 1 ) ) ( I 𝕀 ( ℝ 2 ) ) ) - -> ( C 3 ( I 𝕀 ( ℝ 1 ) ) ( I 𝕀 ( ℝ 1 ) ) ) - -> ( C 3 ( I 𝕀 ( ℝ 1 ) ) ( Spline Closed () ( I 𝕀 ( ℝ 2 ) ) ) ) - -> ( Maybe ( I 𝕀 ( ℝ 1 ) -> I 𝕀 Double ) ) - -> ( I 𝕀 ( ℝ 1 ) -> Seq ( I 𝕀 ( ℝ 1 ) -> StrokeDatum 3 𝕀 ) ) + :: ( 𝕀 -> ( 𝕀ℝ 1 ) ) + -> ( 𝕀ℝ 1 -> 𝕀 ) + -> ( C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 ) ) + -> ( C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 1 ) ) + -> ( C 3 ( 𝕀ℝ 1 ) ( Spline Closed () ( 𝕀ℝ 2 ) ) ) + -> ( Maybe ( 𝕀ℝ 1 -> 𝕀 ) ) + -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) #-} {-# SPECIALISE brushStrokeData - :: ( I 𝕀 Double -> I 𝕀 ( ℝ 1 ) ) - -> ( I 𝕀 ( ℝ 1 ) -> I 𝕀 Double ) - -> ( C 3 ( I 𝕀 ( ℝ 1 ) ) ( I 𝕀 ( ℝ 2 ) ) ) - -> ( C 3 ( I 𝕀 ( ℝ 1 ) ) ( I 𝕀 ( ℝ 2 ) ) ) - -> ( C 3 ( I 𝕀 ( ℝ 2 ) ) ( Spline Closed () ( I 𝕀 ( ℝ 2 ) ) ) ) - -> ( Maybe ( I 𝕀 ( ℝ 2 ) -> I 𝕀 Double ) ) - -> ( I 𝕀 ( ℝ 1 ) -> Seq ( I 𝕀 ( ℝ 1 ) -> StrokeDatum 3 𝕀 ) ) + :: ( 𝕀 -> ( 𝕀ℝ 1 ) ) + -> ( 𝕀ℝ 1 -> 𝕀 ) + -> ( C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 ) ) + -> ( C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 ) ) + -> ( C 3 ( 𝕀ℝ 2 ) ( Spline Closed () ( 𝕀ℝ 2 ) ) ) + -> ( Maybe ( 𝕀ℝ 2 -> 𝕀) ) + -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) #-} {-# SPECIALISE brushStrokeData - :: ( I 𝕀 Double -> I 𝕀 ( ℝ 1 ) ) - -> ( I 𝕀 ( ℝ 1 ) -> I 𝕀 Double ) - -> ( C 3 ( I 𝕀 ( ℝ 1 ) ) ( I 𝕀 ( ℝ 2 ) ) ) - -> ( C 3 ( I 𝕀 ( ℝ 1 ) ) ( I 𝕀 ( ℝ 3 ) ) ) - -> ( C 3 ( I 𝕀 ( ℝ 3 ) ) ( Spline Closed () ( I 𝕀 ( ℝ 2 ) ) ) ) - -> ( Maybe ( I 𝕀 ( ℝ 3 ) -> I 𝕀 Double ) ) - -> ( I 𝕀 ( ℝ 1 ) -> Seq ( I 𝕀 ( ℝ 1 ) -> StrokeDatum 3 𝕀 ) ) + :: ( 𝕀 -> ( 𝕀ℝ 1 ) ) + -> ( 𝕀ℝ 1 -> 𝕀 ) + -> ( C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 ) ) + -> ( C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 3 ) ) + -> ( C 3 ( 𝕀ℝ 3 ) ( Spline Closed () ( 𝕀ℝ 2 ) ) ) + -> ( Maybe ( 𝕀ℝ 3 -> 𝕀 ) ) + -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) #-} {-# SPECIALISE brushStrokeData - :: ( I 𝕀 Double -> I 𝕀 ( ℝ 1 ) ) - -> ( I 𝕀 ( ℝ 1 ) -> I 𝕀 Double ) - -> ( C 3 ( I 𝕀 ( ℝ 1 ) ) ( I 𝕀 ( ℝ 2 ) ) ) - -> ( C 3 ( I 𝕀 ( ℝ 1 ) ) ( I 𝕀 ( ℝ 4 ) ) ) - -> ( C 3 ( I 𝕀 ( ℝ 4 ) ) ( Spline Closed () ( I 𝕀 ( ℝ 2 ) ) ) ) - -> ( Maybe ( I 𝕀 ( ℝ 4 ) -> I 𝕀 Double ) ) - -> ( I 𝕀 ( ℝ 1 ) -> Seq ( I 𝕀 ( ℝ 1 ) -> StrokeDatum 3 𝕀 ) ) + :: ( 𝕀 -> ( 𝕀ℝ 1 ) ) + -> ( 𝕀ℝ 1 -> 𝕀 ) + -> ( C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 ) ) + -> ( C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 4 ) ) + -> ( C 3 ( 𝕀ℝ 4 ) ( Spline Closed () ( 𝕀ℝ 2 ) ) ) + -> ( Maybe ( 𝕀ℝ 4 -> 𝕀 ) ) + -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) #-} +-} -- TODO: these specialisations fire in the benchmarking code because -- we instantiate brushParams with ( ℝ nbBrushParams ), but they won't fire -- in the main app code because we are using types such as "Record EllipseBrushFields". @@ -1097,7 +1104,7 @@ solveEnvelopeEquations :: RootSolvingAlgorithm -> ℝ 2 -> T ( ℝ 2 ) -> ( Offset, Offset ) - -> Seq ( ℝ 1 -> StrokeDatum 2 () ) + -> Seq ( ℝ 1 -> StrokeDatum 2 ℝ ) -> ( ( ℝ 2, T ( ℝ 2 ) ), ( ℝ 2, T ( ℝ 2 ) ) ) solveEnvelopeEquations rootAlgo ( ℝ1 _t ) path_t path'_t ( fwdOffset, bwdOffset ) strokeData = ( fwdSol, ( bwdPt, -1 *^ bwdTgt ) ) @@ -1121,7 +1128,7 @@ solveEnvelopeEquations rootAlgo ( ℝ1 _t ) path_t path'_t ( fwdOffset, bwdOffse | otherwise -> ( off β€’ path_t, path'_t ) - sol :: String -> Seq ( ℝ 1 -> StrokeDatum 2 () ) -> Double -> ( Bool, ℝ 2, T ( ℝ 2 ) ) + sol :: String -> Seq ( ℝ 1 -> StrokeDatum 2 ℝ ) -> Double -> ( Bool, ℝ 2, T ( ℝ 2 ) ) sol desc f is0 = let ( solRes, _solSteps ) = runSolveMethod ( eqn f ) is0 ( good, is ) = @@ -1141,7 +1148,7 @@ solveEnvelopeEquations rootAlgo ( ℝ1 _t ) path_t path'_t ( fwdOffset, bwdOffse NewtonRaphson { maxIters, precision } -> newtonRaphson maxIters precision domain - finish :: Seq ( ℝ 1 -> StrokeDatum 2 () ) -> Double -> ( Double, ℝ 2, T ( ℝ 2 ) ) + finish :: Seq ( ℝ 1 -> StrokeDatum 2 ℝ ) -> Double -> ( Double, ℝ 2, T ( ℝ 2 ) ) finish fs is = let (i, s) = fromDomain is in case evalStrokeDatum fs is of -- TODO: a bit redundant to have to compute this again... @@ -1174,12 +1181,12 @@ solveEnvelopeEquations rootAlgo ( ℝ1 _t ) path_t path'_t ( fwdOffset, bwdOffse -- The sign of this quantity determines which side of the envelope -- we are on. - evalStrokeDatum :: Seq ( ℝ 1 -> StrokeDatum 2 () ) -> ( Double -> StrokeDatum 2 () ) + evalStrokeDatum :: Seq ( ℝ 1 -> StrokeDatum 2 ℝ ) -> ( Double -> StrokeDatum 2 ℝ ) evalStrokeDatum fs is = let (i, s) = fromDomain is in ( fs `Seq.index` i ) ( ℝ1 $ max 1e-6 $ min (1 - 1e-6) $ s ) - eqn :: Seq ( ℝ 1 -> StrokeDatum 2 () ) -> ( Double -> (# Double, Double #) ) + eqn :: Seq ( ℝ 1 -> StrokeDatum 2 ℝ ) -> ( Double -> (# Double, Double #) ) eqn fs is = case evalStrokeDatum fs is of StrokeDatum { ee = D12 ee _ ee_s } -> coerce (# ee, ee_s #) @@ -1202,7 +1209,7 @@ solveEnvelopeEquations rootAlgo ( ℝ1 _t ) path_t path'_t ( fwdOffset, bwdOffse -- -- TODO: use Newton's method starting at the midpoint of the box, -- instead of just taking the midpoint of the box. -cuspCoords :: ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) ) +cuspCoords :: ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 ℝ ) ) -> ( Int, Box N ) -> Cusp cuspCoords eqs ( i, box ) @@ -1231,8 +1238,8 @@ findCusps findCusps opts boxStrokeData = findCuspsIn opts boxStrokeData $ IntMap.fromList - [ ( i, [ 𝕀 ( ℝ2 zero zero ) ( ℝ2 one one ) ] ) --[ 𝕀 ( ℝ3 zero zero -1e6 ) ( ℝ3 one one 1e6 ) ] ) - | i <- [ 0 .. length ( boxStrokeData ( 𝕀 ( ℝ1 zero ) ( ℝ1 one ) ) ) - 1 ] + [ ( i, [ 𝕀ℝ2 ( 𝕀 zero one ) ( 𝕀 zero one ) ] ) + | i <- [ 0 .. length ( boxStrokeData ( 𝕀ℝ1 ( 𝕀 zero one ) ) ) - 1 ] ] where zero, one :: Double @@ -1251,23 +1258,20 @@ findCuspsIn findCuspsIn opts boxStrokeData initBoxes = IntMap.mapWithKey ( \ i -> foldMap ( isolateRootsIn opts ( eqnPiece i ) ) ) initBoxes where - eqnPiece i ( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) = - let t = 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) - s = 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) - -- This conversion is quite painful, but hey... - StrokeDatum + eqnPiece i ( 𝕀ℝ2 t s ) = + let StrokeDatum { ee = - D22 { _D22_v = 𝕀 ( ℝ1 ee_lo ) ( ℝ1 ee_hi ) - , _D22_dx = T ( 𝕀 ( ℝ1 ee_t_lo ) ( ℝ1 ee_t_hi ) ) - , _D22_dy = T ( 𝕀 ( ℝ1 ee_s_lo ) ( ℝ1 ee_s_hi ) ) } + D22 { _D22_v = 𝕀ℝ1 ee + , _D22_dx = T ( 𝕀ℝ1 ee_t ) + , _D22_dy = T ( 𝕀ℝ1 ee_s ) } , 𝛿E𝛿sdcdt = - D12 { _D12_v = T ( 𝕀 ( ℝ2 f_lo g_lo ) ( ℝ2 f_hi g_hi ) ) - , _D12_dx = T ( T ( 𝕀 ( ℝ2 f_t_lo g_t_lo ) ( ℝ2 f_t_hi g_t_hi ) ) ) - , _D12_dy = T ( T ( 𝕀 ( ℝ2 f_s_lo g_s_lo ) ( ℝ2 f_s_hi g_s_hi ) ) ) } - } = ( boxStrokeData t `Seq.index` i ) s - in D12 ( 𝕀 ( ℝ3 ee_lo f_lo g_lo ) ( ℝ3 ee_hi f_hi g_hi ) ) - ( T $ 𝕀 ( ℝ3 ee_t_lo f_t_lo g_t_lo ) ( ℝ3 ee_t_hi f_t_hi g_t_hi ) ) - ( T $ 𝕀 ( ℝ3 ee_s_lo f_s_lo g_s_lo ) ( ℝ3 ee_s_hi f_s_hi g_s_hi ) ) + D12 { _D12_v = T ( 𝕀ℝ2 f g ) + , _D12_dx = T ( T ( 𝕀ℝ2 f_t g_t ) ) + , _D12_dy = T ( T ( 𝕀ℝ2 f_s g_s ) ) } + } = ( boxStrokeData ( 𝕀ℝ1 t ) `Seq.index` i ) ( 𝕀ℝ1 s ) + in D12 ( 𝕀ℝ3 ee f g ) + ( T $ 𝕀ℝ3 ee_t f_t g_t ) + ( T $ 𝕀ℝ3 ee_s f_s g_s ) {- findCuspsIn @@ -1344,14 +1348,14 @@ findCuspsIn opts boxStrokeData initBoxes = ( T $ 𝕀 ( ℝ3 f1_s_lo f2_s_lo f3_s_lo ) ( ℝ3 f1_s_hi f2_s_hi f3_s_hi ) ) ( T $ 𝕀 ( ℝ3 f1_Ξ»_lo f2_Ξ»_lo f3_Ξ»_lo ) ( ℝ3 f1_Ξ»_hi f2_Ξ»_hi f3_Ξ»_hi ) ) -intersectMany :: [𝕀 Double] -> [𝕀 Double] -> [𝕀 Double] +intersectMany :: [𝕀] -> [𝕀] -> [𝕀] intersectMany _ [] = [] intersectMany is (j : js) = intersectOne is j ++ intersectMany is js -intersectOne :: [ 𝕀 Double ] -> 𝕀 Double -> [ 𝕀 Double ] +intersectOne :: [ 𝕀 ] -> 𝕀 -> [ 𝕀 ] intersectOne is i = concatMap ( intersect i ) is -intersect :: 𝕀 Double -> 𝕀 Double -> [ 𝕀 Double ] +intersect :: 𝕀 -> 𝕀 -> [ 𝕀 ] intersect ( 𝕀 lo1 hi1 ) ( 𝕀 lo2 hi2 ) | lo > hi = [ ] diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs b/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs index 60f2710..16b8f99 100644 --- a/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs +++ b/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs @@ -45,35 +45,35 @@ type StrokeDatum :: Nat -> k -> Type data StrokeDatum k i = StrokeDatum { -- | Path \( p(t) \). - dpath :: D k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) + dpath :: D k ( I i 1 ) ( I i 2 ) -- | Brush shape \( b(t, s) \). - , dbrush :: D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) + , dbrush :: D k ( I i 2 ) ( I i 2 ) -- | (Optional) rotation angle \( \theta(t) \). - , mbRotation :: Maybe ( D k ( I i ( ℝ 1 ) ) ( I i Double ) ) + , mbRotation :: Maybe ( D k ( I i 1 ) ( I i Double ) ) -- Everything below is computed in terms of the first three fields. -- | Stroke shape \( c(t,s) = p(t) + R(\theta(t)) b(t,s) \). - , stroke :: I i ( ℝ 2 ) + , stroke :: I i 2 -- | \( u(t,s) = R(-\theta(t)) \frac{\partial c}{\partial t} \), -- \( v(t,s) = R(-\theta(t)) \frac{\partial c}{\partial s} \) - , du, dv :: D ( k - 1 ) ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) + , du, dv :: D ( k - 1 ) ( I i 2 ) ( I i 2 ) -- | Envelope function -- -- \[ E(t_0,s_0) = \left ( \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} \right )_{(t_0,s_0)}. \] - , ee :: D ( k - 1 ) ( I i ( ℝ 2 ) ) ( I i ( ℝ 1 ) ) + , ee :: D ( k - 1 ) ( I i 2 ) ( I i 1 ) -- \[ \frac{\partial E}{\partial s} \frac{\mathrm{d} c}{\mathrm{d} t}, \] -- -- where \( \frac{\mathrm{d} c}{\mathrm{d} t} \) -- -- denotes a total derivative. - , 𝛿E𝛿sdcdt :: D ( k - 2 ) ( I i ( ℝ 2 ) ) ( T ( I i ( ℝ 2 ) ) ) + , 𝛿E𝛿sdcdt :: D ( k - 2 ) ( I i 2 ) ( T ( I i 2 ) ) } -deriving stock instance Show ( StrokeDatum 2 () ) +deriving stock instance Show ( StrokeDatum 2 ℝ ) deriving stock instance Show ( StrokeDatum 3 𝕀 ) -------------------------------------------------------------------------------- @@ -82,33 +82,35 @@ type HasBΓ©zier :: forall {t}. Nat -> t -> Constraint class HasBΓ©zier k i where -- | Linear interpolation, as a differentiable function. - line :: forall b - . ( Module Double ( T b ), Torsor ( T b ) b - , Module ( I i Double ) ( T (I i b ) ), Torsor ( T ( I i b ) ) ( I i b ) - , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) + line :: forall ( n :: Nat ) + . ( Module Double ( T ( ℝ n ) ), Torsor ( T ( ℝ n ) ) ( ℝ n ) + , Module ( I i Double ) ( T ( I i n ) ), Torsor ( T ( I i n ) ) ( I i n ) + , D k ( I i 1 ) ~ D k ( ℝ 1 ) ) - => ( I i ( ℝ 1 ) -> I i Double ) - -> Segment ( I i b ) -> C k ( I i ( ℝ 1 ) ) ( I i b ) + => ( I i 1 -> I i Double ) + -> Segment ( I i n ) -> C k ( I i 1 ) ( I i n ) -- | A quadratic BΓ©zier curve, as a differentiable function. - bezier2 :: forall b - . ( Module Double ( T b ), Torsor ( T b ) b - , Module ( I i Double ) ( T (I i b ) ), Torsor ( T ( I i b ) ) ( I i b ) - , Representable Double b - , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) + bezier2 :: forall ( n :: Nat ) + . ( Module Double ( T ( ℝ n ) ), Torsor ( T ( ℝ n ) ) ( ℝ n ) + , Module ( I i Double ) ( T ( I i n ) ), Torsor ( T ( I i n ) ) ( I i n ) + , Representable Double ( ℝ n ) + , Representable 𝕀 ( 𝕀ℝ n ) + , D k ( I i 1 ) ~ D k ( ℝ 1 ) ) - => ( I i ( ℝ 1 ) -> I i Double ) - -> Quadratic.Bezier ( I i b ) -> C k ( I i ( ℝ 1 ) ) ( I i b ) + => ( I i 1 -> I i Double ) + -> Quadratic.Bezier ( I i n ) -> C k ( I i 1 ) ( I i n ) -- | A cubic BΓ©zier curve, as a differentiable function. - bezier3 :: forall b - . ( Module Double ( T b ), Torsor ( T b ) b - , Module ( I i Double ) ( T (I i b ) ), Torsor ( T ( I i b ) ) ( I i b ) - , Representable Double b - , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) + bezier3 :: forall ( n :: Nat ) + . ( Module Double ( T ( ℝ n ) ), Torsor ( T ( ℝ n ) ) ( ℝ n ) + , Module ( I i Double ) ( T ( I i n ) ), Torsor ( T ( I i n ) ) ( I i n ) + , Representable Double ( ℝ n ) + , Representable 𝕀 ( 𝕀ℝ n ) + , D k ( I i 1 ) ~ D k ( ℝ 1 ) ) - => ( I i ( ℝ 1 ) -> I i Double ) - -> Cubic.Bezier ( I i b ) -> C k ( I i ( ℝ 1 ) ) ( I i b ) + => ( I i 1 -> I i Double ) + -> Cubic.Bezier ( I i n ) -> C k ( I i 1 ) ( I i n ) type HasEnvelopeEquation :: Nat -> Constraint class HasEnvelopeEquation k where @@ -128,24 +130,24 @@ class HasEnvelopeEquation k where -- -- denotes a total derivative. envelopeEquation :: forall i - . ( D ( k - 2 ) ( I i ( ℝ 2 ) ) ~ D ( k - 2 ) ( ℝ 2 ) - , D ( k - 1 ) ( I i ( ℝ 2 ) ) ~ D ( k - 1 ) ( ℝ 2 ) - , D k ( I i ( ℝ 2 ) ) ~ D k ( ℝ 2 ) - , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) - , Module ( I i Double ) ( T ( I i ( ℝ 1 ) ) ) - , Cross ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) + . ( D ( k - 2 ) ( I i 2 ) ~ D ( k - 2 ) ( ℝ 2 ) + , D ( k - 1 ) ( I i 2 ) ~ D ( k - 1 ) ( ℝ 2 ) + , D k ( I i 2 ) ~ D k ( ℝ 2 ) + , D k ( I i 1 ) ~ D k ( ℝ 1 ) + , Module ( I i Double ) ( T ( I i 1 ) ) + , Cross ( I i Double ) ( T ( I i 2 ) ) , Transcendental ( I i Double ) - , Representable ( I i Double ) ( I i ( ℝ 2 ) ) - , RepDim ( I i ( ℝ 2 ) ) ~ 2 + , Representable ( I i Double ) ( I i 2 ) + , RepDim ( I i 2 ) ~ 2 ) => Proxy i - -> ( I i Double -> I i ( ℝ 1 ) ) - -> D k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) - -> D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) - -> Maybe ( D k ( I i ( ℝ 1 ) ) ( I i Double ) ) + -> ( I i Double -> I i 1 ) + -> D k ( I i 1 ) ( I i 2 ) + -> D k ( I i 2 ) ( I i 2 ) + -> Maybe ( D k ( I i 1 ) ( I i Double ) ) -> StrokeDatum k i -instance HasBΓ©zier 2 () where +instance HasBΓ©zier 2 ℝ where line co ( Segment a b :: Segment b ) = D \ ( co -> t ) -> D21 ( lerp @( T b ) t a b ) @@ -217,14 +219,14 @@ instance HasEnvelopeEquation 2 where -- = ( R(-ΞΈ(t)) βˆ‚c/βˆ‚t ) Γ— ( R(-ΞΈ(t)) βˆ‚c/ds ) -- = ( R(-ΞΈ(t)) p'(t) + ΞΈ'(t) S b(t,s) + βˆ‚b/βˆ‚t ) Γ— βˆ‚b/ds -- - let rot, rot' :: T ( I i ( ℝ 2 ) ) -> T ( I i ( ℝ 2 ) ) + let rot, rot' :: T ( I i 2 ) -> T ( I i 2 ) cosΞΈ = cos ΞΈ sinΞΈ = sin ΞΈ -- rot = R(-ΞΈ), rot' = R'(-ΞΈ) -- NB: rot' is not the derivative of f(ΞΈ) = R(-ΞΈ) rot = rotate cosΞΈ -sinΞΈ rot' = rotate sinΞΈ cosΞΈ - swap :: T ( I i ( ℝ 2 ) ) -> T ( I i ( ℝ 2 ) ) + swap :: T ( I i 2 ) -> T ( I i 2 ) swap ( T xy ) = let x = xy `index` Fin 1 y = xy `index` Fin 2 @@ -232,7 +234,7 @@ instance HasEnvelopeEquation 2 where Fin 1 -> -y _ -> x - u, v, u_t, u_s, v_t, v_s :: T ( I i ( ℝ 2 ) ) + u, v, u_t, u_s, v_t, v_s :: T ( I i 2 ) u = rot p_t ^+^ ΞΈ_t *^ swap b ^+^ b_t v = b_s u_t = ( -ΞΈ_t *^ rot' p_t @@ -253,7 +255,7 @@ instance HasEnvelopeEquation 2 where ) {-# INLINEABLE envelopeEquation #-} -instance HasBΓ©zier 3 () where +instance HasBΓ©zier 3 ℝ where line co ( Segment a b :: Segment b ) = D \ ( co -> t ) -> @@ -347,7 +349,7 @@ instance HasEnvelopeEquation 3 where -- = ( R(-ΞΈ(t)) βˆ‚c/βˆ‚t ) Γ— ( R(-ΞΈ(t)) βˆ‚c/ds ) -- = ( R(-ΞΈ(t)) p'(t) + ΞΈ'(t) S b(t,s) + βˆ‚b/βˆ‚t ) Γ— βˆ‚b/ds -- - let rot, rot', rot'' :: T ( I i ( ℝ 2 ) ) -> T ( I i ( ℝ 2 ) ) + let rot, rot', rot'' :: T ( I i 2 ) -> T ( I i 2 ) cosΞΈ = cos ΞΈ sinΞΈ = sin ΞΈ -- rot = R(-ΞΈ), rot' = R'(-ΞΈ), rot'' = R''(-ΞΈ) @@ -355,7 +357,7 @@ instance HasEnvelopeEquation 3 where rot = rotate cosΞΈ -sinΞΈ rot' = rotate sinΞΈ cosΞΈ rot'' z = -1 *^ rot z - swap :: T ( I i ( ℝ 2 ) ) -> T ( I i ( ℝ 2 ) ) + swap :: T ( I i 2 ) -> T ( I i 2 ) swap ( T xy ) = let x = xy `index` Fin 1 y = xy `index` Fin 2 @@ -363,7 +365,7 @@ instance HasEnvelopeEquation 3 where Fin 1 -> -y _ -> x - u, v, u_t, u_s, u_tt, u_ts, u_ss, v_t, v_s, v_tt, v_ts, v_ss :: T ( I i ( ℝ 2 ) ) + u, v, u_t, u_s, u_tt, u_ts, u_ss, v_t, v_s, v_tt, v_ts, v_ss :: T ( I i 2 ) u = rot p_t ^+^ ΞΈ_t *^ swap b ^+^ b_t v = b_s @@ -458,7 +460,7 @@ with respect to t and root-finding. -- | Evaluate a cubic BΓ©zier curve, when both the coefficients and the -- parameter are intervals. -evaluateCubic :: Cubic.Bezier ( 𝕀 Double ) -> 𝕀 Double -> 𝕀 Double +evaluateCubic :: Cubic.Bezier 𝕀 -> 𝕀 -> 𝕀 evaluateCubic bez t = -- assert (inf t >= 0 && sup t <= 1) "evaluateCubic: t ⊊ [0,1]" $ -- Requires t βŠ† [0,1] let inf_bez = fmap inf bez @@ -472,7 +474,7 @@ evaluateCubic bez t = -- | Evaluate a quadratic BΓ©zier curve, when both the coefficients and the -- parameter are intervals. -evaluateQuadratic :: Quadratic.Bezier ( 𝕀 Double ) -> 𝕀 Double -> 𝕀 Double +evaluateQuadratic :: Quadratic.Bezier 𝕀 -> 𝕀 -> 𝕀 evaluateQuadratic bez t = -- assert (inf t >= 0 && sup t <= 1) "evaluateCubic: t ⊊ [0,1]" $ -- Requires t βŠ† [0,1] let inf_bez = fmap inf bez diff --git a/brush-strokes/src/lib/Math/Differentiable.hs b/brush-strokes/src/lib/Math/Differentiable.hs index 521bd9a..b6c102a 100644 --- a/brush-strokes/src/lib/Math/Differentiable.hs +++ b/brush-strokes/src/lib/Math/Differentiable.hs @@ -15,7 +15,7 @@ import GHC.TypeNats import Math.Algebra.Dual ( D, HasChainRule ) import Math.Interval - ( 𝕀 ) + ( 𝕀, 𝕀ℝ ) import Math.Linear import Math.Module import Math.Ring @@ -25,13 +25,15 @@ import Math.Ring -- | Type family to parametrise over both pointwise and interval computations. -- --- Use '()' parameter for points, and '𝕀' parameter for intervals. -type I :: k -> Type -> Type -type family I i a -type instance I () a = a -type instance I 𝕀 a = 𝕀 a +-- Use 'ℝ' parameter for points, and '𝕀' parameter for intervals. +type I :: k -> l -> Type +type family I i t +type instance I @(Nat -> Type) @Type ℝ Double = Double +type instance I @Type @Type 𝕀 Double = 𝕀 +type instance I @(Nat -> Type) @Nat ℝ n = ℝ n +type instance I @Type @Nat 𝕀 n = 𝕀ℝ n -type Differentiable :: Nat -> k -> Type -> Constraint +type Differentiable :: Nat -> k -> Nat -> Constraint class ( Module ( I i Double ) ( T ( I i u ) ) , HasChainRule ( I i Double ) k ( I i u ) @@ -43,13 +45,13 @@ instance , Traversable ( D k ( I i u ) ) ) => Differentiable k i u -type DiffInterp :: Nat -> k -> Type -> Constraint +type DiffInterp :: Nat -> k -> Nat -> Constraint class ( Differentiable k i u , Interpolatable ( I i Double ) ( I i u ) , Module ( I i Double ) ( T ( I i Double ) ) , Module ( D k ( I i u ) ( I i Double ) ) - ( D k ( I i u ) ( I i ( ℝ 2 ) ) ) + ( D k ( I i u ) ( I i 2 ) ) , Transcendental ( D k ( I i u ) ( I i Double ) ) , Applicative ( D k ( I i u ) ) , Representable ( I i Double ) ( I i u ) @@ -60,7 +62,7 @@ instance , Interpolatable ( I i Double ) ( I i u ) , Module ( I i Double ) ( T ( I i Double ) ) , Module ( D k ( I i u ) ( I i Double ) ) - ( D k ( I i u ) ( I i ( ℝ 2 ) ) ) + ( D k ( I i u ) ( I i 2 ) ) , Transcendental ( D k ( I i u ) ( I i Double ) ) , Applicative ( D k ( I i u ) ) , Representable ( I i Double ) ( I i u ) diff --git a/brush-strokes/src/lib/Math/Interval.hs b/brush-strokes/src/lib/Math/Interval.hs index a97408d..203497e 100644 --- a/brush-strokes/src/lib/Math/Interval.hs +++ b/brush-strokes/src/lib/Math/Interval.hs @@ -8,14 +8,17 @@ module Math.Interval ( 𝕀(𝕀), inf, sup - , width, midpoint - , scaleInterval - , 𝕀ℝ + , midpoint , scaleInterval, width + , isCanonical - , singleton, nonDecreasing - , (∈), (βŠ–), (∩), (⊘) + , (∩), (⊘), (∈), (βŠ–) , extendedRecip + + , singleton, point + , nonDecreasing , aabb, bisect + + , 𝕀ℝ(..) ) where @@ -40,7 +43,10 @@ import Math.Algebra.Dual import Math.Float.Utils ( prevFP ) import Math.Interval.Internal - ( 𝕀(𝕀), inf, sup, scaleInterval ) + ( 𝕀(𝕀), inf, sup + , scaleInterval, width, singleton + , 𝕀ℝ(..) + ) import Math.Linear ( ℝ(..), T(..) , RepresentableQ(..), Representable(..) @@ -52,46 +58,53 @@ import Math.Ring -------------------------------------------------------------------------------- -- Interval arithmetic. -type 𝕀ℝ n = 𝕀 ( ℝ n ) -type instance D k ( 𝕀 v ) = D k v - --- | An interval reduced to a single point. -singleton :: a -> 𝕀 a -singleton a = 𝕀 a a - --- | The width of an interval. --- --- NB: assumes the endpoints are neither @NaN@ nor infinity. -width :: AbelianGroup a => 𝕀 a -> a -width ( 𝕀 lo hi ) = hi - lo -{-# INLINEABLE width #-} +type instance D k 𝕀 = D k Double +type instance D k ( 𝕀ℝ n ) = D k ( ℝ n ) -- | Turn a non-decreasing function into a function on intervals. -nonDecreasing :: ( a -> b ) -> 𝕀 a -> 𝕀 b -nonDecreasing f ( 𝕀 lo hi ) = 𝕀 ( f lo ) ( f hi ) +nonDecreasing :: forall n m + . ( Representable Double ( ℝ n ) + , Representable 𝕀 ( 𝕀ℝ n ) + , Representable Double ( ℝ m ) + , Representable 𝕀 ( 𝕀ℝ m ) ) + => ( ℝ n -> ℝ m ) -> 𝕀ℝ n -> 𝕀ℝ m +nonDecreasing f u = + let + lo, hi :: ℝ n + lo = tabulate \ i -> inf ( u `index` i ) + hi = tabulate \ i -> sup ( u `index` i ) + lo', hi' :: ℝ m + lo' = f lo + hi' = f hi + in tabulate \ j -> 𝕀 ( lo' `index` j ) ( hi' `index` j ) +{-# INLINEABLE nonDecreasing #-} -- | Does the given value lie inside the specified interval? -(∈) :: Ord a => a -> 𝕀 a -> Bool +(∈) :: Double -> 𝕀 -> Bool x ∈ 𝕀 lo hi = x >= lo && x <= hi -{-# INLINEABLE (∈) #-} infix 4 ∈ -- | Is this interval canonical, i.e. it consists of either 1 or 2 floating point -- values only? -isCanonical :: 𝕀 Double -> Bool +isCanonical :: 𝕀 -> Bool isCanonical ( 𝕀 x_inf x_sup ) = x_inf >= prevFP x_sup -- | The midpoint of an interval. -- -- NB: assumes the endpoints are neither @NaN@ nor infinity. -midpoint :: 𝕀 Double -> Double +midpoint :: 𝕀 -> Double midpoint ( 𝕀 x_inf x_sup ) = 0.5 * ( x_inf + x_sup ) +point :: ( Representable Double ( ℝ n ), Representable 𝕀 ( 𝕀ℝ n ) ) + => ℝ n -> 𝕀ℝ n +point x = tabulate \ i -> singleton ( x `index` i ) +{-# INLINEABLE point #-} + -- | Compute the intersection of two intervals. -- -- Returns whether the first interval is a strict subset of the second interval -- (or the intersection is a single point). -(∩) :: 𝕀 Double -> 𝕀 Double -> [ ( 𝕀 Double, Bool ) ] +(∩) :: 𝕀 -> 𝕀 -> [ ( 𝕀, Bool ) ] 𝕀 lo1 hi1 ∩ 𝕀 lo2 hi2 | lo > hi = [ ] @@ -102,11 +115,19 @@ midpoint ( 𝕀 x_inf x_sup ) = 0.5 * ( x_inf + x_sup ) hi = min hi1 hi2 infix 3 ∩ +infixl 6 βŠ– +(βŠ–) :: 𝕀 -> 𝕀 -> 𝕀 +(βŠ–) a@( 𝕀 lo1 hi1 ) b@( 𝕀 lo2 hi2 ) + | width a >= width b + = 𝕀 ( lo1 - lo2 ) ( hi1 - hi2 ) + | otherwise + = 𝕀 ( hi1 - hi2 ) ( lo1 - lo2 ) + -- | Bisect an interval. -- -- Generally returns two sub-intervals, but returns the input interval if -- it is canonical (and thus bisection would be pointless). -bisect :: 𝕀 Double -> NE.NonEmpty ( 𝕀 Double ) +bisect :: 𝕀 -> NE.NonEmpty 𝕀 bisect x@( 𝕀 x_inf x_sup ) | isCanonical x = NE.singleton x @@ -114,25 +135,16 @@ bisect x@( 𝕀 x_inf x_sup ) = 𝕀 x_inf x_mid NE.:| [ 𝕀 x_mid x_sup ] where x_mid = midpoint x -infixl 6 βŠ– -(βŠ–) :: ( AbelianGroup a, Ord a ) => 𝕀 a -> 𝕀 a -> 𝕀 a -(βŠ–) a@( 𝕀 lo1 hi1 ) b@( 𝕀 lo2 hi2 ) - | width a >= width b - = 𝕀 ( lo1 - lo2 ) ( hi1 - hi2 ) - | otherwise - = 𝕀 ( hi1 - hi2 ) ( lo1 - lo2 ) -{-# INLINEABLE (βŠ–) #-} +deriving via ViaAbelianGroup ( T 𝕀 ) + instance Semigroup ( T 𝕀 ) +deriving via ViaAbelianGroup ( T 𝕀 ) + instance Monoid ( T 𝕀 ) +deriving via ViaAbelianGroup ( T 𝕀 ) + instance Group ( T 𝕀 ) -deriving via ViaAbelianGroup ( T ( 𝕀 Double ) ) - instance Semigroup ( T ( 𝕀 Double ) ) -deriving via ViaAbelianGroup ( T ( 𝕀 Double ) ) - instance Monoid ( T ( 𝕀 Double ) ) -deriving via ViaAbelianGroup ( T ( 𝕀 Double ) ) - instance Group ( T ( 𝕀 Double ) ) - -instance Act ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where +instance Act ( T 𝕀 ) 𝕀 where T g β€’ a = g + a -instance Torsor ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where +instance Torsor ( T 𝕀 ) 𝕀 where a --> b = T $ b - a ------------------------------------------------------------------------------- @@ -142,7 +154,7 @@ instance Torsor ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where -- -- NB: this function can return overlapping intervals if both arguments contain 0. -- Otherwise, the returned intervals are guaranteed to be disjoint. -(⊘) :: 𝕀 Double -> 𝕀 Double -> [ 𝕀 Double ] +(⊘) :: 𝕀 -> 𝕀 -> [ 𝕀 ] x ⊘ y #ifdef ASSERTS | 0 ∈ x && 0 ∈ y @@ -155,7 +167,7 @@ x ⊘ y infixl 7 ⊘ -- | Extended reciprocal, returning either 1 or 2 intervals. -extendedRecip :: 𝕀 Double -> [ 𝕀 Double ] +extendedRecip :: 𝕀 -> [ 𝕀 ] extendedRecip x@( 𝕀 lo hi ) | lo == 0 && hi == 0 = [ 𝕀 negInf posInf ] @@ -172,85 +184,77 @@ posInf = 1 / 0 ------------------------------------------------------------------------------- -- Lattices. -aabb :: ( Representable r v, Ord r, Functor f ) - => f ( 𝕀 v ) -> ( f ( 𝕀 r ) -> 𝕀 r ) -> 𝕀 v +aabb :: ( Representable 𝕀 ( 𝕀ℝ n ), Functor f ) + => f ( 𝕀ℝ n ) -> ( f 𝕀 -> 𝕀 ) -> 𝕀ℝ n aabb fv extrema = tabulate \ i -> extrema ( fmap ( `index` i ) fv ) {-# INLINEABLE aabb #-} -------------------------------------------------------------------------------- -instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ 0 ) ) where - origin = T $ singleton ℝ0 - _ ^+^ _ = T $ singleton ℝ0 - _ ^-^ _ = T $ singleton ℝ0 - _ *^ _ = T $ singleton ℝ0 +instance Module 𝕀 ( T ( 𝕀ℝ 0 ) ) where + origin = T 𝕀ℝ0 + _ ^+^ _ = T 𝕀ℝ0 + _ ^-^ _ = T 𝕀ℝ0 + _ *^ _ = T 𝕀ℝ0 -instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ 1 ) ) where - origin = T $$( tabulateQ \ _ -> [|| unT $ origin ||] ) - T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] ) - T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] ) - k *^ T a = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] ) +instance Module 𝕀 ( T ( 𝕀ℝ 1 ) ) where + origin = T $ 𝕀ℝ1 ( unT origin ) + T ( 𝕀ℝ1 x1 ) ^+^ T ( 𝕀ℝ1 x2 ) = T $ 𝕀ℝ1 ( x1 + x2 ) + T ( 𝕀ℝ1 x1 ) ^-^ T ( 𝕀ℝ1 x2 ) = T $ 𝕀ℝ1 ( x1 - x2 ) + k *^ T ( 𝕀ℝ1 x1 ) = T $ 𝕀ℝ1 ( k * x1 ) -instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ 2 ) ) where - origin = T $$( tabulateQ \ _ -> [|| unT $ origin ||] ) - T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] ) - T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] ) - k *^ T a = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] ) +instance Module 𝕀 ( T ( 𝕀ℝ 2 ) ) where + origin = T $ 𝕀ℝ2 ( unT origin ) ( unT origin ) + T ( 𝕀ℝ2 x1 y1 ) ^+^ T ( 𝕀ℝ2 x2 y2 ) = T $ 𝕀ℝ2 ( x1 + x2 ) ( y1 + y2 ) + T ( 𝕀ℝ2 x1 y1 ) ^-^ T ( 𝕀ℝ2 x2 y2 ) = T $ 𝕀ℝ2 ( x1 - x2 ) ( y1 - y2 ) + k *^ T ( 𝕀ℝ2 x1 y1 ) = T $ 𝕀ℝ2 ( k * x1 ) ( k * y1 ) -instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ 3 ) ) where - origin = T $$( tabulateQ \ _ -> [|| unT $ origin ||] ) - T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] ) - T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] ) - k *^ T a = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] ) +instance Module 𝕀 ( T ( 𝕀ℝ 3 ) ) where + origin = T $ 𝕀ℝ3 ( unT origin ) ( unT origin ) ( unT origin ) + T ( 𝕀ℝ3 x1 y1 z1 ) ^+^ T ( 𝕀ℝ3 x2 y2 z2 ) = T $ 𝕀ℝ3 ( x1 + x2 ) ( y1 + y2 ) ( z1 + z2 ) + T ( 𝕀ℝ3 x1 y1 z1 ) ^-^ T ( 𝕀ℝ3 x2 y2 z2 ) = T $ 𝕀ℝ3 ( x1 - x2 ) ( y1 - y2 ) ( z1 - z2 ) + k *^ T ( 𝕀ℝ3 x1 y1 z1 ) = T $ 𝕀ℝ3 ( k * x1 ) ( k * y1 ) ( k * z1 ) -instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ 4 ) ) where - origin = T $$( tabulateQ \ _ -> [|| unT $ origin ||] ) - T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] ) - T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] ) - k *^ T a = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] ) +instance Module 𝕀 ( T ( 𝕀ℝ 4 ) ) where + origin = T $ 𝕀ℝ4 ( unT origin ) ( unT origin ) ( unT origin ) ( unT origin ) + T ( 𝕀ℝ4 x1 y1 z1 w1 ) ^+^ T ( 𝕀ℝ4 x2 y2 z2 w2 ) = T $ 𝕀ℝ4 ( x1 + x2 ) ( y1 + y2 ) ( z1 + z2 ) ( w1 + w2 ) + T ( 𝕀ℝ4 x1 y1 z1 w1 ) ^-^ T ( 𝕀ℝ4 x2 y2 z2 w2 ) = T $ 𝕀ℝ4 ( x1 - x2 ) ( y1 - y2 ) ( z1 - z2 ) ( w1 - w2 ) + k *^ T ( 𝕀ℝ4 x1 y1 z1 w1 ) = T $ 𝕀ℝ4 ( k * x1 ) ( k * y1 ) ( k * z1 ) ( k * w1 ) -instance Inner ( 𝕀 Double ) ( T ( 𝕀ℝ 2 ) ) where - T ( 𝕀 ( ℝ2 x1_lo y1_lo ) ( ℝ2 x1_hi y1_hi ) ) ^.^ - T ( 𝕀 ( ℝ2 x2_lo y2_lo ) ( ℝ2 x2_hi y2_hi ) ) - = let !x1x2 = 𝕀 x1_lo x1_hi * 𝕀 x2_lo x2_hi - !y1y2 = 𝕀 y1_lo y1_hi * 𝕀 y2_lo y2_hi - in x1x2 + y1y2 +instance Inner 𝕀 ( T ( 𝕀ℝ 2 ) ) where + T ( 𝕀ℝ2 x1 y1 ) ^.^ T ( 𝕀ℝ2 x2 y2 ) = x1 * x2 + y1 * y2 -instance Cross ( 𝕀 Double ) ( T ( 𝕀ℝ 2 ) ) where - T ( 𝕀 ( ℝ2 x1_lo y1_lo ) ( ℝ2 x1_hi y1_hi ) ) Γ— - T ( 𝕀 ( ℝ2 x2_lo y2_lo ) ( ℝ2 x2_hi y2_hi ) ) - = let !x1y2 = 𝕀 x1_lo x1_hi * 𝕀 y2_lo y2_hi - !y2x1 = 𝕀 x2_lo x2_hi * 𝕀 y1_lo y1_hi - in x1y2 - y2x1 +instance Cross 𝕀 ( T ( 𝕀ℝ 2 ) ) where + T ( 𝕀ℝ2 x1 y1 ) Γ— T ( 𝕀ℝ2 x2 y2 ) = x1 * y2 - x2 * y1 -deriving via ViaModule ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) - instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) => Semigroup ( T ( 𝕀ℝ n ) ) -deriving via ViaModule ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) - instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) => Monoid ( T ( 𝕀ℝ n ) ) -deriving via ViaModule ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) - instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) => Group ( T ( 𝕀ℝ n ) ) +deriving via ViaModule 𝕀 ( T ( 𝕀ℝ n ) ) + instance Module 𝕀 ( T ( 𝕀ℝ n ) ) => Semigroup ( T ( 𝕀ℝ n ) ) +deriving via ViaModule 𝕀 ( T ( 𝕀ℝ n ) ) + instance Module 𝕀 ( T ( 𝕀ℝ n ) ) => Monoid ( T ( 𝕀ℝ n ) ) +deriving via ViaModule 𝕀 ( T ( 𝕀ℝ n ) ) + instance Module 𝕀 ( T ( 𝕀ℝ n ) ) => Group ( T ( 𝕀ℝ n ) ) -deriving via ViaModule ( 𝕀 Double ) ( 𝕀ℝ n ) - instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) => Act ( T ( 𝕀ℝ n ) ) ( 𝕀ℝ n ) -deriving via ( ViaModule ( 𝕀 Double ) ( 𝕀ℝ n ) ) - instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) => Torsor ( T ( 𝕀ℝ n ) ) ( 𝕀ℝ n ) +deriving via ViaModule 𝕀 ( 𝕀ℝ n ) + instance Module 𝕀 ( T ( 𝕀ℝ n ) ) => Act ( T ( 𝕀ℝ n ) ) ( 𝕀ℝ n ) +deriving via ( ViaModule 𝕀 ( 𝕀ℝ n ) ) + instance Module 𝕀 ( T ( 𝕀ℝ n ) ) => Torsor ( T ( 𝕀ℝ n ) ) ( 𝕀ℝ n ) -------------------------------------------------------------------------------- -- HasChainRule instances. -instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 0 ) where +instance HasChainRule 𝕀 2 ( 𝕀ℝ 0 ) where konst w = D0 w linearD f v = D0 ( f v ) value ( D0 v ) = v chain _ ( D0 gfx ) = D21 gfx origin origin -instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 0 ) where +instance HasChainRule 𝕀 3 ( 𝕀ℝ 0 ) where konst w = D0 w linearD f v = D0 ( f v ) value ( D0 v ) = v chain _ ( D0 gfx ) = D31 gfx origin origin origin -instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 1 ) where +instance HasChainRule 𝕀 2 ( 𝕀ℝ 1 ) where konst :: forall w. AbelianGroup w => w -> D2𝔸1 w konst w = @@ -259,9 +263,9 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 1 ) where value df = $$( monIndexQ [|| df ||] zeroMonomial ) - linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 1 -> w ) -> 𝕀ℝ 1 -> D2𝔸1 w + linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 1 -> w ) -> 𝕀ℝ 1 -> D2𝔸1 w linearD f v = - let !o = origin @( 𝕀 Double ) @( T w ) + let !o = origin @𝕀 @( T w ) in $$( monTabulateQ \ mon -> if | isZeroMonomial mon -> [|| f v ||] @@ -276,16 +280,16 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 1 ) where -> [|| unT o ||] ) - chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D2𝔸1 ( 𝕀ℝ 1 ) -> D2𝔸1 w -> D2𝔸1 w + chain :: forall w. Module 𝕀 ( T w ) => D2𝔸1 ( 𝕀ℝ 1 ) -> D2𝔸1 w -> D2𝔸1 w chain !df !dg = - let !o = origin @( 𝕀 Double ) @( T w ) - !p = (^+^) @( 𝕀 Double ) @( T w ) - !s = (^*) @( 𝕀 Double ) @( T w ) + let !o = origin @𝕀 @( T w ) + !p = (^+^) @𝕀 @( T w ) + !s = (^*) @𝕀 @( T w ) in $$( chainRule1NQ [|| o ||] [|| p ||] [|| s ||] [|| df ||] [|| dg ||] ) -instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 1 ) where +instance HasChainRule 𝕀 3 ( 𝕀ℝ 1 ) where konst :: forall w. AbelianGroup w => w -> D3𝔸1 w konst w = @@ -294,9 +298,9 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 1 ) where value df = $$( monIndexQ [|| df ||] zeroMonomial ) - linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 1 -> w ) -> 𝕀ℝ 1 -> D3𝔸1 w + linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 1 -> w ) -> 𝕀ℝ 1 -> D3𝔸1 w linearD f v = - let !o = origin @( 𝕀 Double ) @( T w ) + let !o = origin @𝕀 @( T w ) in $$( monTabulateQ \ mon -> if | isZeroMonomial mon -> [|| f v ||] @@ -311,16 +315,16 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 1 ) where -> [|| unT o ||] ) - chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D3𝔸1 ( 𝕀ℝ 1 ) -> D3𝔸1 w -> D3𝔸1 w + chain :: forall w. Module 𝕀 ( T w ) => D3𝔸1 ( 𝕀ℝ 1 ) -> D3𝔸1 w -> D3𝔸1 w chain !df !dg = - let !o = origin @( 𝕀 Double ) @( T w ) - !p = (^+^) @( 𝕀 Double ) @( T w ) - !s = (^*) @( 𝕀 Double ) @( T w ) + let !o = origin @𝕀 @( T w ) + !p = (^+^) @𝕀 @( T w ) + !s = (^*) @𝕀 @( T w ) in $$( chainRule1NQ [|| o ||] [|| p ||] [|| s ||] [|| df ||] [|| dg ||] ) -instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 2 ) where +instance HasChainRule 𝕀 2 ( 𝕀ℝ 2 ) where konst :: forall w. AbelianGroup w => w -> D2𝔸2 w konst w = @@ -329,9 +333,9 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 2 ) where value df = $$( monIndexQ [|| df ||] zeroMonomial ) - linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 2 -> w ) -> 𝕀ℝ 2 -> D2𝔸2 w + linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 2 -> w ) -> 𝕀ℝ 2 -> D2𝔸2 w linearD f v = - let !o = origin @( 𝕀 Double ) @( T w ) + let !o = origin @𝕀 @( T w ) in $$( monTabulateQ \ mon -> if | isZeroMonomial mon -> [|| f v ||] @@ -346,16 +350,16 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 2 ) where -> [|| unT o ||] ) - chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D2𝔸1 ( 𝕀ℝ 2 ) -> D2𝔸2 w -> D2𝔸1 w + chain :: forall w. Module 𝕀 ( T w ) => D2𝔸1 ( 𝕀ℝ 2 ) -> D2𝔸2 w -> D2𝔸1 w chain !df !dg = - let !o = origin @( 𝕀 Double ) @( T w ) - !p = (^+^) @( 𝕀 Double ) @( T w ) - !s = (^*) @( 𝕀 Double ) @( T w ) + let !o = origin @𝕀 @( T w ) + !p = (^+^) @𝕀 @( T w ) + !s = (^*) @𝕀 @( T w ) in $$( chainRule1NQ [|| o ||] [|| p ||] [|| s ||] [|| df ||] [|| dg ||] ) -instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 2 ) where +instance HasChainRule 𝕀 3 ( 𝕀ℝ 2 ) where konst :: forall w. AbelianGroup w => w -> D3𝔸2 w konst w = @@ -364,9 +368,9 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 2 ) where value df = $$( monIndexQ [|| df ||] zeroMonomial ) - linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 2 -> w ) -> 𝕀ℝ 2 -> D3𝔸2 w + linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 2 -> w ) -> 𝕀ℝ 2 -> D3𝔸2 w linearD f v = - let !o = origin @( 𝕀 Double ) @( T w ) + let !o = origin @𝕀 @( T w ) in $$( monTabulateQ \ mon -> if | isZeroMonomial mon -> [|| f v ||] @@ -381,16 +385,16 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 2 ) where -> [|| unT o ||] ) - chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D3𝔸1 ( 𝕀ℝ 2 ) -> D3𝔸2 w -> D3𝔸1 w + chain :: forall w. Module 𝕀 ( T w ) => D3𝔸1 ( 𝕀ℝ 2 ) -> D3𝔸2 w -> D3𝔸1 w chain !df !dg = - let !o = origin @( 𝕀 Double ) @( T w ) - !p = (^+^) @( 𝕀 Double ) @( T w ) - !s = (^*) @( 𝕀 Double ) @( T w ) + let !o = origin @𝕀 @( T w ) + !p = (^+^) @𝕀 @( T w ) + !s = (^*) @𝕀 @( T w ) in $$( chainRule1NQ [|| o ||] [|| p ||] [|| s ||] [|| df ||] [|| dg ||] ) -instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 3 ) where +instance HasChainRule 𝕀 2 ( 𝕀ℝ 3 ) where konst :: forall w. AbelianGroup w => w -> D2𝔸3 w konst w = @@ -399,9 +403,9 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 3 ) where value df = $$( monIndexQ [|| df ||] zeroMonomial ) - linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 3 -> w ) -> 𝕀ℝ 3 -> D2𝔸3 w + linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 3 -> w ) -> 𝕀ℝ 3 -> D2𝔸3 w linearD f v = - let !o = origin @( 𝕀 Double ) @( T w ) + let !o = origin @𝕀 @( T w ) in $$( monTabulateQ \ mon -> if | isZeroMonomial mon -> [|| f v ||] @@ -416,16 +420,16 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 3 ) where -> [|| unT o ||] ) - chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D2𝔸1 ( 𝕀ℝ 3 ) -> D2𝔸3 w -> D2𝔸1 w + chain :: forall w. Module 𝕀 ( T w ) => D2𝔸1 ( 𝕀ℝ 3 ) -> D2𝔸3 w -> D2𝔸1 w chain !df !dg = - let !o = origin @( 𝕀 Double ) @( T w ) - !p = (^+^) @( 𝕀 Double ) @( T w ) - !s = (^*) @( 𝕀 Double ) @( T w ) + let !o = origin @𝕀 @( T w ) + !p = (^+^) @𝕀 @( T w ) + !s = (^*) @𝕀 @( T w ) in $$( chainRule1NQ [|| o ||] [|| p ||] [|| s ||] [|| df ||] [|| dg ||] ) -instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 3 ) where +instance HasChainRule 𝕀 3 ( 𝕀ℝ 3 ) where konst :: forall w. AbelianGroup w => w -> D3𝔸3 w konst w = @@ -434,9 +438,9 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 3 ) where value df = $$( monIndexQ [|| df ||] zeroMonomial ) - linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 3 -> w ) -> 𝕀ℝ 3 -> D3𝔸3 w + linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 3 -> w ) -> 𝕀ℝ 3 -> D3𝔸3 w linearD f v = - let !o = origin @( 𝕀 Double ) @( T w ) + let !o = origin @𝕀 @( T w ) in $$( monTabulateQ \ mon -> if | isZeroMonomial mon -> [|| f v ||] @@ -451,16 +455,16 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 3 ) where -> [|| unT o ||] ) - chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D3𝔸1 ( 𝕀ℝ 3 ) -> D3𝔸3 w -> D3𝔸1 w + chain :: forall w. Module 𝕀 ( T w ) => D3𝔸1 ( 𝕀ℝ 3 ) -> D3𝔸3 w -> D3𝔸1 w chain !df !dg = - let !o = origin @( 𝕀 Double ) @( T w ) - !p = (^+^) @( 𝕀 Double ) @( T w ) - !s = (^*) @( 𝕀 Double ) @( T w ) + let !o = origin @𝕀 @( T w ) + !p = (^+^) @𝕀 @( T w ) + !s = (^*) @𝕀 @( T w ) in $$( chainRule1NQ [|| o ||] [|| p ||] [|| s ||] [|| df ||] [|| dg ||] ) -instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 4 ) where +instance HasChainRule 𝕀 2 ( 𝕀ℝ 4 ) where konst :: forall w. AbelianGroup w => w -> D2𝔸4 w konst w = @@ -469,9 +473,9 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 4 ) where value df = $$( monIndexQ [|| df ||] zeroMonomial ) - linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 4 -> w ) -> 𝕀ℝ 4 -> D2𝔸4 w + linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 4 -> w ) -> 𝕀ℝ 4 -> D2𝔸4 w linearD f v = - let !o = origin @( 𝕀 Double ) @( T w ) + let !o = origin @𝕀 @( T w ) in $$( monTabulateQ \ mon -> if | isZeroMonomial mon -> [|| f v ||] @@ -486,16 +490,16 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ 4 ) where -> [|| unT o ||] ) - chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D2𝔸1 ( 𝕀ℝ 4 ) -> D2𝔸4 w -> D2𝔸1 w + chain :: forall w. Module 𝕀 ( T w ) => D2𝔸1 ( 𝕀ℝ 4 ) -> D2𝔸4 w -> D2𝔸1 w chain !df !dg = - let !o = origin @( 𝕀 Double ) @( T w ) - !p = (^+^) @( 𝕀 Double ) @( T w ) - !s = (^*) @( 𝕀 Double ) @( T w ) + let !o = origin @𝕀 @( T w ) + !p = (^+^) @𝕀 @( T w ) + !s = (^*) @𝕀 @( T w ) in $$( chainRule1NQ [|| o ||] [|| p ||] [|| s ||] [|| df ||] [|| dg ||] ) -instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 4 ) where +instance HasChainRule 𝕀 3 ( 𝕀ℝ 4 ) where konst :: forall w. AbelianGroup w => w -> D3𝔸4 w konst w = @@ -504,9 +508,9 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 4 ) where value df = $$( monIndexQ [|| df ||] zeroMonomial ) - linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀ℝ 4 -> w ) -> 𝕀ℝ 4 -> D3𝔸4 w + linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 4 -> w ) -> 𝕀ℝ 4 -> D3𝔸4 w linearD f v = - let !o = origin @( 𝕀 Double ) @( T w ) + let !o = origin @𝕀 @( T w ) in $$( monTabulateQ \ mon -> if | isZeroMonomial mon -> [|| f v ||] @@ -521,11 +525,11 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ 4 ) where -> [|| unT o ||] ) - chain :: forall w. Module ( 𝕀 Double ) ( T w ) => D3𝔸1 ( 𝕀ℝ 4 ) -> D3𝔸4 w -> D3𝔸1 w + chain :: forall w. Module 𝕀 ( T w ) => D3𝔸1 ( 𝕀ℝ 4 ) -> D3𝔸4 w -> D3𝔸1 w chain !df !dg = - let !o = origin @( 𝕀 Double ) @( T w ) - !p = (^+^) @( 𝕀 Double ) @( T w ) - !s = (^*) @( 𝕀 Double ) @( T w ) + let !o = origin @𝕀 @( T w ) + !p = (^+^) @𝕀 @( T w ) + !s = (^*) @𝕀 @( T w ) in $$( chainRule1NQ [|| o ||] [|| p ||] [|| s ||] [|| df ||] [|| dg ||] ) diff --git a/brush-strokes/src/lib/Math/Interval/FMA.hs b/brush-strokes/src/lib/Math/Interval/FMA.hs deleted file mode 100644 index d5ca39d..0000000 --- a/brush-strokes/src/lib/Math/Interval/FMA.hs +++ /dev/null @@ -1,92 +0,0 @@ -{-# LANGUAGE AllowAmbiguousTypes #-} -{-# LANGUAGE ScopedTypeVariables #-} - -module Math.Interval.FMA - ( addI, subI, prodI, divI, posPowI ) - where - --- base -import GHC.Exts - ( Double(D#), fmsubDouble#, fnmaddDouble# ) - --- brush-strokes -import Math.Float.Utils - ( prevFP, succFP ) - --------------------------------------------------------------------------------- - -{-# INLINE withError #-} -withError :: Double -> Double -> ( Double, Double ) -withError f e = - case compare e 0 of - LT -> ( prevFP f, f ) - EQ -> ( f, f ) - GT -> ( f, succFP f ) - -addI :: Double -> Double -> ( Double, Double ) -addI a b = - let !s = a + b - !a' = s - b - !b' = s - a' - !da = a - a' - !db = b - b' - !e = da + db - in s `withError` e - -subI :: Double -> Double -> ( Double, Double ) -subI a b = - let !s = a - b - !a' = s - b - !b' = s - a' - !da = a - a' - !db = b - b' - !e = da + db - in s `withError` e - -{-# INLINE fmsubDouble #-} -fmsubDouble :: Double -> Double -> Double -> Double -fmsubDouble ( D# x ) ( D# y ) ( D# z ) = D# ( fmsubDouble# x y z ) -{-# INLINE fnmaddDouble #-} -fnmaddDouble :: Double -> Double -> Double -> Double -fnmaddDouble ( D# x ) ( D# y ) ( D# z ) = D# ( fnmaddDouble# x y z ) - -prodI :: Double -> Double -> ( Double, Double ) -prodI a b = - let !p = a * b - !e = fmsubDouble a b p - in p `withError` e - -divI :: Double -> Double -> ( Double, Double ) -divI a b = - let !r = a / b - !e = fnmaddDouble r b a - in r `withError` e - --- | Power of a **non-negative** number to a natural power. -posPowI :: Double -- ^ Assumed to be non-negative! - -> Word -> ( Double, Double ) -posPowI _ 0 = ( 1, 1 ) -posPowI f 1 = ( f, f ) -posPowI f 2 = prodI f f -posPowI f n - | even n - , let m = n `quot` 2 - ( fΒ²_lo, fΒ²_hi ) = prodI f f - = ( fst $ posPowI fΒ²_lo m, snd $ posPowI fΒ²_hi m ) - | otherwise - , let m = n `quot` 2 - ( fΒ²_lo, fΒ²_hi ) = prodI f f - = ( fst $ posPowAcc fΒ²_lo m f, snd $ posPowAcc fΒ²_hi m f ) - -posPowAcc :: Double -> Word -> Double -> ( Double, Double ) -posPowAcc f 1 x = prodI f x -posPowAcc f n x - | even n - , let m = n `quot` 2 - ( fΒ²_lo, fΒ²_hi ) = prodI f f - = ( fst $ posPowAcc fΒ²_lo m x, snd $ posPowAcc fΒ²_hi m x ) - | otherwise - , let m = n `quot` 2 - ( fΒ²_lo, fΒ²_hi ) = prodI f f - ( y_lo, y_hi ) = prodI f x - = ( fst $ posPowAcc fΒ²_lo m y_lo, snd $ posPowAcc fΒ²_hi m y_hi ) diff --git a/brush-strokes/src/lib/Math/Interval/Internal.hs b/brush-strokes/src/lib/Math/Interval/Internal.hs index 3a81335..271bd98 100644 --- a/brush-strokes/src/lib/Math/Interval/Internal.hs +++ b/brush-strokes/src/lib/Math/Interval/Internal.hs @@ -6,190 +6,68 @@ {-# OPTIONS_GHC -Wno-orphans #-} module Math.Interval.Internal - ( 𝕀(𝕀, inf, sup), scaleInterval ) + ( 𝕀(𝕀, inf, sup) + , scaleInterval, width, singleton + , 𝕀ℝ(..) + ) where -- base import Prelude hiding ( Num(..), Fractional(..), Floating(..), (^) ) -import qualified Prelude +import Data.Kind + ( Type ) import Data.Monoid ( Sum(..) ) +import GHC.Generics + ( Generic ) import GHC.Show ( showCommaSpace ) +import GHC.TypeNats + ( Nat ) -- deepseq import Control.DeepSeq ( NFData(..) ) --- rounded-hw -import Numeric.Rounded.Hardware - ( Rounded(..) ) -import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval - ( Interval(..) ) - -#ifdef USE_FMA -- brush-strokes -import Math.Interval.FMA - ( addI, subI, prodI, divI, posPowI ) - +#if defined(USE_SIMD) +import Math.Interval.Internal.SIMD +#elif defined(USE_FMA) +import Math.Interval.Internal.FMA #else --- rounded-hw -import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval - ( powInt ) +import Math.Interval.Internal.RoundedHW #endif -- brush-strokes import Math.Linear ( T(..) , RepDim, RepresentableQ(..), Representable(..) + , Fin(..) ) import Math.Module ( Module(..) ) import Math.Ring -------------------------------------------------------------------------------- --- Intervals. -#ifdef USE_FMA -data 𝕀 a = 𝕀 { inf, sup :: !a } +-- | An interval reduced to a single point. +singleton :: Double -> 𝕀 +singleton x = 𝕀 x x -instance NFData a => NFData ( 𝕀 a ) where - rnf ( 𝕀 lo hi ) = rnf lo `seq` rnf hi - -instance Prelude.Num ( 𝕀 Double ) where - 𝕀 x_lo x_hi + 𝕀 y_lo y_hi - | let !z_lo = fst $ addI x_lo y_lo - !z_hi = snd $ addI x_hi y_hi - = 𝕀 z_lo z_hi - 𝕀 x_lo x_hi - 𝕀 y_lo y_hi - | let !z_lo = fst $ subI x_lo y_hi - !z_hi = snd $ subI x_hi y_lo - = 𝕀 z_lo z_hi - negate (𝕀 lo hi) = 𝕀 (Prelude.negate hi) (Prelude.negate lo) - (*) = (*) - fromInteger i = - let !j = Prelude.fromInteger i - in 𝕀 j j - abs (𝕀 lo hi) - | 0 <= lo - = 𝕀 lo hi - | hi <= 0 - = 𝕀 (Prelude.negate hi) (Prelude.negate lo) - | otherwise - = 𝕀 0 (max (Prelude.negate lo) hi) - signum _ = error "No implementation of signum for intervals" - -instance Ring ( 𝕀 Double ) where - 𝕀 lo1 hi1 * 𝕀 lo2 hi2 - | let !( x_min, x_max ) = prodI lo1 lo2 - !( y_min, y_max ) = prodI lo1 hi2 - !( z_min, z_max ) = prodI hi1 lo2 - !( w_min, w_max ) = prodI hi1 hi2 - = 𝕀 ( min ( min x_min y_min ) ( min z_min w_min ) ) - ( max ( max x_max y_max ) ( max z_max w_max ) ) - _ ^ 0 = 𝕀 1 1 - iv ^ 1 = iv - 𝕀 lo hi ^ n - | odd n || 0 <= lo - , let !lo' = fst $ posPowI lo n - !hi' = snd $ posPowI hi n - = 𝕀 lo' hi' - | hi <= 0 - , let !lo' = fst $ posPowI (negate hi) n - !hi' = snd $ posPowI (negate lo) n - = 𝕀 lo' hi' - | otherwise - , let !hi1 = snd $ posPowI (negate lo) n - , let !hi2 = snd $ posPowI hi n - = 𝕀 0 ( max hi1 hi2 ) - -instance Prelude.Fractional ( 𝕀 Double ) where - fromRational r = - let q = Prelude.fromRational r - in 𝕀 q q - recip ( 𝕀 lo hi ) --- #ifdef ASSERTS - | lo == 0 - = 𝕀 ( fst $ divI 1 hi ) ( 1 Prelude./ 0 ) - | hi == 0 - = 𝕀 ( -1 Prelude./ 0 ) ( snd $ divI 1 lo ) - | lo > 0 || hi < 0 --- #endif - = 𝕀 ( fst $ divI 1 hi ) ( snd $ divI 1 lo ) --- #ifdef ASSERTS - | otherwise - = error "BAD interval recip; should use extendedRecip instead" --- #endif - p / q = p * Prelude.recip q - -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 - -{- -TODO: consider using https://github.com/JishinMaster/simd_utils/blob/160c50f07e08d2077ae4368f0aed2f10f7173c67/simd_utils_sse_double.h#L530 -for sin/cos? Or Intel SVML? Or metalibm? - --} - -{-# INLINE withHW #-} --- | Internal function: use @rounded-hw@ to define a function on intervals. -withHW :: (Interval.Interval a -> Interval.Interval b) -> 𝕀 a -> 𝕀 b -withHW f = \ ( 𝕀 lo hi ) -> - case f ( Interval.I ( Rounded lo ) ( Rounded hi ) ) of - Interval.I ( Rounded x ) ( Rounded y ) -> 𝕀 x y - -scaleInterval :: Double -> 𝕀 Double -> 𝕀 Double -scaleInterval s ( 𝕀 lo hi ) = - case compare s 0 of - LT -> 𝕀 ( fst $ prodI s hi ) ( snd $ prodI s lo ) - EQ -> 𝕀 0 0 - GT -> 𝕀 ( fst $ prodI s lo ) ( snd $ prodI s hi ) -#else - -newtype 𝕀 a = MkI { ival :: Interval.Interval a } - deriving newtype ( Prelude.Num, Prelude.Fractional, Prelude.Floating ) - deriving newtype NFData - -{-# COMPLETE 𝕀 #-} -pattern 𝕀 :: a -> a -> 𝕀 a -pattern 𝕀 { inf, sup } = MkI ( Interval.I ( Rounded inf ) ( Rounded sup ) ) - -instance Ring ( 𝕀 Double ) where - MkI i1 * MkI i2 = MkI $ i1 Prelude.* i2 - MkI x ^ n = MkI { ival = Interval.powInt x ( Prelude.fromIntegral n ) } - -- This is very important, as x^2 is not the same as x * x - -- in interval arithmetic. This ensures we don't - -- accidentally use (^) from Prelude. - -deriving via ViaPrelude ( 𝕀 Double ) - instance Floating ( 𝕀 Double ) -deriving via ViaPrelude ( 𝕀 Double ) - instance Transcendental ( 𝕀 Double ) - -scaleInterval :: Double -> 𝕀 Double -> 𝕀 Double -scaleInterval s iv = 𝕀 s s * iv -- TODO: could be better -#endif - -deriving via ViaPrelude ( 𝕀 Double ) - instance AbelianGroup ( 𝕀 Double ) -deriving via ViaPrelude ( 𝕀 Double ) - instance AbelianGroup ( T ( 𝕀 Double ) ) -deriving via ViaPrelude ( 𝕀 Double ) - instance Field ( 𝕀 Double ) +-- | The width of an interval. +-- +-- NB: assumes the endpoints are neither @NaN@ nor infinity. +width :: 𝕀 -> Double +width ( 𝕀 lo hi ) = hi - lo -------------------------------------------------------------------------------- +-- Instances for (1D) intervals. -instance Eq a => Eq ( 𝕀 a ) where +instance Eq 𝕀 where 𝕀 a b == 𝕀 c d = a == c && b == d -instance Show a => Show ( 𝕀 a ) where +instance Show 𝕀 where showsPrec _ ( 𝕀 x y ) = showString "[" . showsPrec 0 x @@ -197,35 +75,104 @@ instance Show a => Show ( 𝕀 a ) where . showsPrec 0 y . showString "]" +deriving via ViaPrelude 𝕀 + instance AbelianGroup ( T 𝕀 ) +deriving via Sum 𝕀 + instance Module 𝕀 ( T 𝕀 ) + -------------------------------------------------------------------------------- -type instance RepDim ( 𝕀 u ) = RepDim u -instance RepresentableQ r u => RepresentableQ ( 𝕀 r ) ( 𝕀 u ) where - tabulateQ f = - let !lo = tabulateQ @r @u ( \ i -> [|| inf $ $$( f i ) ||] ) - !hi = tabulateQ @r @u ( \ i -> [|| sup $ $$( f i ) ||] ) - in [|| 𝕀 $$lo $$hi ||] +type 𝕀ℝ :: Nat -> Type +data family 𝕀ℝ n - indexQ d i = - [|| case $$d of - 𝕀 lo hi -> - let !lo_i = $$( indexQ @r @u [|| lo ||] i ) - !hi_i = $$( indexQ @r @u [|| hi ||] i ) - in 𝕀 lo_i hi_i - ||] -instance Representable r u => Representable ( 𝕀 r ) ( 𝕀 u ) where - tabulate f = - let !lo = tabulate @r @u ( \ i -> inf $ f i ) - !hi = tabulate @r @u ( \ i -> sup $ f i ) - in 𝕀 lo hi +data instance 𝕀ℝ 0 = 𝕀ℝ0 + deriving stock ( Show, Eq, Ord, Generic ) + deriving anyclass NFData +newtype instance 𝕀ℝ 1 = 𝕀ℝ1 { un𝕀ℝ1 :: 𝕀 } + deriving stock ( Show, Generic ) + deriving newtype ( Eq, NFData ) +data instance 𝕀ℝ 2 = 𝕀ℝ2 { _𝕀ℝ2_x, _𝕀ℝ2_y :: !𝕀 } + deriving stock Generic + deriving anyclass NFData + deriving stock ( Show, Eq ) +data instance 𝕀ℝ 3 = 𝕀ℝ3 { _𝕀ℝ3_x, _𝕀ℝ3_y, _𝕀ℝ3_z :: !𝕀 } + deriving stock Generic + deriving anyclass NFData + deriving stock ( Show, Eq ) +data instance 𝕀ℝ 4 = 𝕀ℝ4 { _𝕀ℝ4_x, _𝕀ℝ4_y, _𝕀ℝ4_z, _𝕀ℝ4_w :: !𝕀 } + deriving stock Generic + deriving anyclass NFData + deriving stock ( Show, Eq ) + +type instance RepDim ( 𝕀ℝ n ) = n + +instance RepresentableQ 𝕀 ( 𝕀ℝ 0 ) where + tabulateQ _ = [|| 𝕀ℝ0 ||] + indexQ _ _ = [|| 0 ||] + +instance RepresentableQ 𝕀 ( 𝕀ℝ 1 ) where + tabulateQ f = [|| 𝕀ℝ1 $$( f ( Fin 1 ) ) ||] + indexQ p = \ case + _ -> [|| un𝕀ℝ1 $$p ||] + +instance RepresentableQ 𝕀 ( 𝕀ℝ 2 ) where + tabulateQ f = [|| 𝕀ℝ2 $$( f ( Fin 1 ) ) $$( f ( Fin 2 ) ) ||] + indexQ p = \ case + Fin 1 -> [|| _𝕀ℝ2_x $$p ||] + _ -> [|| _𝕀ℝ2_y $$p ||] + +instance RepresentableQ 𝕀 ( 𝕀ℝ 3 ) where + tabulateQ f = [|| 𝕀ℝ3 $$( f ( Fin 1 ) ) $$( f ( Fin 2 ) ) $$( f ( Fin 3 ) ) ||] + indexQ p = \ case + Fin 1 -> [|| _𝕀ℝ3_x $$p ||] + Fin 2 -> [|| _𝕀ℝ3_y $$p ||] + _ -> [|| _𝕀ℝ3_z $$p ||] + +instance RepresentableQ 𝕀 ( 𝕀ℝ 4 ) where + tabulateQ f = [|| 𝕀ℝ4 $$( f ( Fin 1 ) ) $$( f ( Fin 2 ) ) $$( f ( Fin 3 ) ) $$( f ( Fin 4 ) ) ||] + indexQ p = \ case + Fin 1 -> [|| _𝕀ℝ4_x $$p ||] + Fin 2 -> [|| _𝕀ℝ4_y $$p ||] + Fin 3 -> [|| _𝕀ℝ4_z $$p ||] + _ -> [|| _𝕀ℝ4_w $$p ||] + + +instance Representable 𝕀 ( 𝕀ℝ 0 ) where + tabulate _ = 𝕀ℝ0 {-# INLINE tabulate #-} - - index d i = - case d of - 𝕀 lo hi -> - let !lo_i = index @r @u lo i - !hi_i = index @r @u hi i - in 𝕀 lo_i hi_i + index _ _ = 0 {-# INLINE index #-} -deriving via Sum ( 𝕀 Double ) instance Module ( 𝕀 Double ) ( T ( 𝕀 Double ) ) +instance Representable 𝕀 ( 𝕀ℝ 1 ) where + tabulate f = 𝕀ℝ1 ( f ( Fin 1 ) ) + {-# INLINE tabulate #-} + index p = \ case + _ -> un𝕀ℝ1 p + {-# INLINE index #-} + +instance Representable 𝕀 ( 𝕀ℝ 2 ) where + tabulate f = 𝕀ℝ2 ( f ( Fin 1 ) ) ( f ( Fin 2 ) ) + {-# INLINE tabulate #-} + index p = \ case + Fin 1 -> _𝕀ℝ2_x p + _ -> _𝕀ℝ2_y p + {-# INLINE index #-} + +instance Representable 𝕀 ( 𝕀ℝ 3 ) where + tabulate f = 𝕀ℝ3 ( f ( Fin 1 ) ) ( f ( Fin 2 ) ) ( f ( Fin 3 ) ) + {-# INLINE tabulate #-} + index p = \ case + Fin 1 -> _𝕀ℝ3_x p + Fin 2 -> _𝕀ℝ3_y p + _ -> _𝕀ℝ3_z p + {-# INLINE index #-} + +instance Representable 𝕀 ( 𝕀ℝ 4 ) where + tabulate f = 𝕀ℝ4 ( f ( Fin 1 ) ) ( f ( Fin 2 ) ) ( f ( Fin 3 ) ) ( f ( Fin 4 ) ) + {-# INLINE tabulate #-} + index p = \ case + Fin 1 -> _𝕀ℝ4_x p + Fin 2 -> _𝕀ℝ4_y p + Fin 3 -> _𝕀ℝ4_z p + _ -> _𝕀ℝ4_w p + {-# INLINE index #-} diff --git a/brush-strokes/src/lib/Math/Interval/Internal/FMA.hs b/brush-strokes/src/lib/Math/Interval/Internal/FMA.hs new file mode 100644 index 0000000..2779350 --- /dev/null +++ b/brush-strokes/src/lib/Math/Interval/Internal/FMA.hs @@ -0,0 +1,214 @@ +{-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} + +module Math.Interval.Internal.FMA + ( 𝕀(..) + , scaleInterval + ) where + +-- base +import Prelude hiding ( Num(..), Fractional(..), Floating(..) ) +import qualified Prelude +import GHC.Exts + ( Double(D#), fmsubDouble#, fnmaddDouble# ) + +-- brush-strokes +import Math.Float.Utils + ( prevFP, succFP ) +import Math.Ring + +-- deepseq +import Control.DeepSeq + ( NFData(..) ) + +-- rounded-hw +import Numeric.Rounded.Hardware + ( Rounded(..) ) +import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval + ( Interval(..) ) + +-------------------------------------------------------------------------------- + +{-# INLINE withError #-} +withError :: Double -> Double -> ( Double, Double ) +withError f e = + case compare e 0 of + LT -> ( prevFP f, f ) + EQ -> ( f, f ) + GT -> ( f, succFP f ) + +addI :: Double -> Double -> ( Double, Double ) +addI a b = + let !s = a + b + !a' = s - b + !b' = s - a' + !da = a - a' + !db = b - b' + !e = da + db + in s `withError` e + +subI :: Double -> Double -> ( Double, Double ) +subI a b = + let !s = a - b + !a' = s - b + !b' = s - a' + !da = a - a' + !db = b - b' + !e = da + db + in s `withError` e + +{-# INLINE fmsubDouble #-} +fmsubDouble :: Double -> Double -> Double -> Double +fmsubDouble ( D# x ) ( D# y ) ( D# z ) = D# ( fmsubDouble# x y z ) +{-# INLINE fnmaddDouble #-} +fnmaddDouble :: Double -> Double -> Double -> Double +fnmaddDouble ( D# x ) ( D# y ) ( D# z ) = D# ( fnmaddDouble# x y z ) + +prodI :: Double -> Double -> ( Double, Double ) +prodI a b = + let !p = a * b + !e = fmsubDouble a b p + in p `withError` e + +divI :: Double -> Double -> ( Double, Double ) +divI a b = + let !r = a / b + !e = fnmaddDouble r b a + in r `withError` e + +-- | Power of a **non-negative** number to a natural power. +posPowI :: Double -- ^ Assumed to be non-negative! + -> Word -> ( Double, Double ) +posPowI _ 0 = ( 1, 1 ) +posPowI f 1 = ( f, f ) +posPowI f 2 = prodI f f +posPowI f n + | even n + , let m = n `quot` 2 + ( fΒ²_lo, fΒ²_hi ) = prodI f f + = ( fst $ posPowI fΒ²_lo m, snd $ posPowI fΒ²_hi m ) + | otherwise + , let m = n `quot` 2 + ( fΒ²_lo, fΒ²_hi ) = prodI f f + = ( fst $ posPowAcc fΒ²_lo m f, snd $ posPowAcc fΒ²_hi m f ) + +posPowAcc :: Double -> Word -> Double -> ( Double, Double ) +posPowAcc f 1 x = prodI f x +posPowAcc f n x + | even n + , let m = n `quot` 2 + ( fΒ²_lo, fΒ²_hi ) = prodI f f + = ( fst $ posPowAcc fΒ²_lo m x, snd $ posPowAcc fΒ²_hi m x ) + | otherwise + , let m = n `quot` 2 + ( fΒ²_lo, fΒ²_hi ) = prodI f f + ( y_lo, y_hi ) = prodI f x + = ( fst $ posPowAcc fΒ²_lo m y_lo, snd $ posPowAcc fΒ²_hi m y_hi ) + +-------------------------------------------------------------------------------- + +-- | A non-empty interval of real numbers (possibly unbounded). +data 𝕀 = 𝕀 { inf, sup :: !Double } + +instance NFData 𝕀 where + rnf ( 𝕀 lo hi ) = rnf lo `seq` rnf hi + +instance Prelude.Num 𝕀 where + 𝕀 x_lo x_hi + 𝕀 y_lo y_hi + | let !z_lo = fst $ addI x_lo y_lo + !z_hi = snd $ addI x_hi y_hi + = 𝕀 z_lo z_hi + 𝕀 x_lo x_hi - 𝕀 y_lo y_hi + | let !z_lo = fst $ subI x_lo y_hi + !z_hi = snd $ subI x_hi y_lo + = 𝕀 z_lo z_hi + negate (𝕀 lo hi) = 𝕀 (Prelude.negate hi) (Prelude.negate lo) + (*) = (*) + fromInteger i = + let !j = Prelude.fromInteger i + in 𝕀 j j + abs (𝕀 lo hi) + | 0 <= lo + = 𝕀 lo hi + | hi <= 0 + = 𝕀 (Prelude.negate hi) (Prelude.negate lo) + | otherwise + = 𝕀 0 (max (Prelude.negate lo) hi) + signum _ = error "No implementation of signum for intervals" + +instance Ring 𝕀 where + 𝕀 lo1 hi1 * 𝕀 lo2 hi2 + | let !( x_min, x_max ) = prodI lo1 lo2 + !( y_min, y_max ) = prodI lo1 hi2 + !( z_min, z_max ) = prodI hi1 lo2 + !( w_min, w_max ) = prodI hi1 hi2 + = 𝕀 ( min ( min x_min y_min ) ( min z_min w_min ) ) + ( max ( max x_max y_max ) ( max z_max w_max ) ) + _ ^ 0 = 𝕀 1 1 + iv ^ 1 = iv + 𝕀 lo hi ^ n + | odd n || 0 <= lo + , let !lo' = fst $ posPowI lo n + !hi' = snd $ posPowI hi n + = 𝕀 lo' hi' + | hi <= 0 + , let !lo' = fst $ posPowI (negate hi) n + !hi' = snd $ posPowI (negate lo) n + = 𝕀 lo' hi' + | otherwise + , let !hi1 = snd $ posPowI (negate lo) n + , let !hi2 = snd $ posPowI hi n + = 𝕀 0 ( max hi1 hi2 ) + +instance Prelude.Fractional 𝕀 where + fromRational r = + let q = Prelude.fromRational r + in 𝕀 q q + recip ( 𝕀 lo hi ) +-- #ifdef ASSERTS + | lo == 0 + = 𝕀 ( fst $ divI 1 hi ) ( 1 Prelude./ 0 ) + | hi == 0 + = 𝕀 ( -1 Prelude./ 0 ) ( snd $ divI 1 lo ) + | lo > 0 || hi < 0 +-- #endif + = 𝕀 ( fst $ divI 1 hi ) ( snd $ divI 1 lo ) +-- #ifdef ASSERTS + | otherwise + = error "BAD interval recip; should use extendedRecip instead" +-- #endif + p / q = p * Prelude.recip q + +instance Floating 𝕀 where + sqrt = withHW Prelude.sqrt + +instance Transcendental 𝕀 where + pi = 𝕀 3.141592653589793 3.1415926535897936 + cos = withHW Prelude.cos + sin = withHW Prelude.sin + atan = withHW Prelude.atan + +deriving via ViaPrelude 𝕀 + instance AbelianGroup 𝕀 +deriving via ViaPrelude 𝕀 + instance Field 𝕀 + +{-# INLINE withHW #-} +-- | Internal function: use @rounded-hw@ to define a function on intervals. +withHW :: (Interval.Interval Double -> Interval.Interval Double) -> 𝕀 -> 𝕀 +withHW f = \ ( 𝕀 lo hi ) -> + case f ( Interval.I ( Rounded lo ) ( Rounded hi ) ) of + Interval.I ( Rounded x ) ( Rounded y ) -> 𝕀 x y + +-- | The width of an interval. +-- +-- NB: assumes the endpoints are neither @NaN@ nor infinity. +width :: 𝕀 -> Double +width ( 𝕀 lo hi ) = hi - lo + +scaleInterval :: Double -> 𝕀 -> 𝕀 +scaleInterval s ( 𝕀 lo hi ) = + case compare s 0 of + LT -> 𝕀 ( fst $ prodI s hi ) ( snd $ prodI s lo ) + EQ -> 𝕀 0 0 + GT -> 𝕀 ( fst $ prodI s lo ) ( snd $ prodI s hi ) diff --git a/brush-strokes/src/lib/Math/Interval/Internal/RoundedHW.hs b/brush-strokes/src/lib/Math/Interval/Internal/RoundedHW.hs new file mode 100644 index 0000000..55e97a2 --- /dev/null +++ b/brush-strokes/src/lib/Math/Interval/Internal/RoundedHW.hs @@ -0,0 +1,53 @@ + +module Math.Interval.Internal.RoundedHW + ( 𝕀(𝕀, inf, sup) + , scaleInterval + ) + where + +-- base +import Prelude hiding ( Num(..), Fractional(..), Floating(..) ) +import qualified Prelude + +-- brush-strokes +import Math.Ring + +-- deepseq +import Control.DeepSeq + ( NFData(..) ) + +-- rounded-hw +import Numeric.Rounded.Hardware + ( Rounded(..) ) +import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval + ( Interval(..), powInt ) + +-------------------------------------------------------------------------------- + +-- | A non-empty interval of real numbers (possibly unbounded). +newtype 𝕀 = MkI { ival :: Interval.Interval Double } + deriving newtype ( Prelude.Num, Prelude.Fractional, Prelude.Floating ) + deriving newtype NFData + +{-# COMPLETE 𝕀 #-} +pattern 𝕀 :: Double -> Double -> 𝕀 +pattern 𝕀 { inf, sup } = MkI ( Interval.I ( Rounded inf ) ( Rounded sup ) ) + +instance Ring 𝕀 where + MkI i1 * MkI i2 = MkI $ i1 Prelude.* i2 + MkI x ^ n = MkI { ival = Interval.powInt x ( Prelude.fromIntegral n ) } + -- This is very important, as x^2 is not the same as x * x + -- in interval arithmetic. This ensures we don't + -- accidentally use (^) from Prelude. + +deriving via ViaPrelude 𝕀 + instance AbelianGroup 𝕀 +deriving via ViaPrelude 𝕀 + instance Field 𝕀 +deriving via ViaPrelude 𝕀 + instance Floating 𝕀 +deriving via ViaPrelude 𝕀 + instance Transcendental 𝕀 + +scaleInterval :: Double -> 𝕀 -> 𝕀 +scaleInterval s iv = 𝕀 s s * iv -- TODO: could be better diff --git a/brush-strokes/src/lib/Math/Interval/Internal/SIMD.hs b/brush-strokes/src/lib/Math/Interval/Internal/SIMD.hs new file mode 100644 index 0000000..15f3267 --- /dev/null +++ b/brush-strokes/src/lib/Math/Interval/Internal/SIMD.hs @@ -0,0 +1,173 @@ +{-# LANGUAGE CPP #-} + +module Math.Interval.Internal.SIMD + ( 𝕀(𝕀, inf, sup) + , scaleInterval + ) + where + +-- base +import Prelude hiding ( Num(..), Fractional(..), Floating(..) ) +import qualified Prelude +import GHC.Exts + ( isTrue# + , Double(D#), Double# + , (<=##), (*##), negateDouble# + , DoubleX2# + , packDoubleX2#, unpackDoubleX2# + , broadcastDoubleX2# + , plusDoubleX2#, timesDoubleX2# + ) + +-- ghc-prim +import GHC.Prim + ( maxDouble#, shuffleDoubleX2# ) + +-- deepseq +import Control.DeepSeq + ( NFData(..) ) + +-- rounded-hw +import Numeric.Rounded.Hardware + ( Rounded(..) ) +import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval + ( Interval(..) ) + +-- brush-strokes +import Math.Ring + +-------------------------------------------------------------------------------- + +-- | A non-empty interval of real numbers (possibly unbounded). +data 𝕀 = MkI DoubleX2# + -- Interval [lo..hi] represented as (hi, -lo) + +instance NFData 𝕀 where + rnf ( MkI !_ ) = () + +{-# COMPLETE PackI #-} +pattern PackI :: Double# -> Double# -> 𝕀 +pattern PackI hi m_lo <- ( ( \ ( MkI x ) -> unpackDoubleX2# x ) -> (# hi, m_lo #) ) + where + PackI hi m_lo = MkI ( packDoubleX2# (# hi, m_lo #) ) + + +{-# COMPLETE 𝕀 #-} +pattern 𝕀 :: Double -> Double -> 𝕀 +pattern 𝕀 { inf, sup } <- ( ( \ ( PackI hi m_lo ) -> (# D# ( negateDouble# m_lo ), D# hi #) ) + -> (# inf, sup #) ) + where + 𝕀 (D# lo) (D# hi) = + let m_lo = negateDouble# lo + in PackI hi m_lo + +deriving via ViaPrelude 𝕀 + instance AbelianGroup 𝕀 +deriving via ViaPrelude 𝕀 + instance Field 𝕀 + +instance Prelude.Num 𝕀 where + MkI x + MkI y = MkI ( x `plusDoubleX2#` y ) + negate ( MkI x ) = MkI ( shuffleDoubleX2# x x (# 1#, 0# #) ) + (*) = (*) + fromInteger i = + let !j = Prelude.fromInteger i + in 𝕀 j j + abs x@(PackI hi m_lo) + | D# m_lo <= 0 + = PackI hi m_lo + | D# hi <= 0 + = Prelude.negate x + | otherwise + = 𝕀 0 (max (D# m_lo) (D# hi)) + signum _ = error "No implementation of signum for intervals" + + +instance Ring 𝕀 where + ( PackI b ma ) * ( PackI d mc ) = + case ( zeroPos b ma, zeroPos d mc ) of + -- From "Fast and Correct SIMD Algorithms for Interval Arithmetic" + -- FrΓ©dΓ©ric Goualard, 2008 + ( LT, LT ) -> PackI (ma *## mc) ( ( negateDouble# b ) *## d) + ( LT, EQ ) -> PackI (ma *## mc) ( ma *## d ) + ( LT, GT ) -> PackI ( b *## ( negateDouble# mc ) ) ( ma *## d ) + ( EQ, LT ) -> PackI (ma *## mc) ( b *## mc ) + ( EQ, EQ ) -> PackI + ( maxDouble# ( ma *## d ) ( b *## mc ) ) + ( maxDouble# ( ma *## mc ) ( b *## d ) ) + ( EQ, GT ) -> PackI ( b *## d ) ( ma *## d ) + ( GT, LT ) -> PackI ( ( negateDouble# ma ) *## d ) ( b *## mc ) + ( GT, EQ ) -> PackI ( b *## d ) ( b *## mc ) + ( GT, GT ) -> PackI ( b *## d ) ( ( negateDouble# ma ) *## mc ) + where + zeroPos :: Double# -> Double# -> Ordering + zeroPos hi m_lo + | isTrue# (hi <=## 0.0##) + = LT + | isTrue# (m_lo <=## 0.0##) + = GT + | otherwise + = EQ + iv ^ 1 = iv + PackI hi m_lo ^ n + | odd n || isTrue# (m_lo <=## 0.0##) + , let !(D# m_lo') = D# m_lo Prelude.^ n + !(D# hi' ) = D# hi Prelude.^ n + = PackI hi' m_lo' + | isTrue# (hi <=## 0.0##) + , let !(D# m_lo') = (D# hi ) Prelude.^ n + !(D# hi') = (D# m_lo) Prelude.^ n + = PackI hi' m_lo' + | otherwise + , let !(D# hi1) = (D# m_lo) Prelude.^ n + !(D# hi2) = (D# hi ) Prelude.^ n + hi' = maxDouble# hi1 hi2 + = PackI hi' 0.0## + +instance Prelude.Fractional 𝕀 where + fromRational r = + let q = Prelude.fromRational r + in 𝕀 q q + recip ( 𝕀 lo hi ) +#ifdef ASSERTS + | lo == 0 + = 𝕀 ( Prelude.recip hi ) ( 1 Prelude./ 0 ) + | hi == 0 + = 𝕀 ( -1 Prelude./ 0 ) ( Prelude.recip lo ) + | lo > 0 || hi < 0 +#endif + = 𝕀 ( Prelude.recip hi ) ( Prelude.recip lo ) +#ifdef ASSERTS + | otherwise + = error "BAD interval recip; should use extendedRecip instead" +#endif + p / q = p * Prelude.recip q + +instance Floating 𝕀 where + sqrt = withHW Prelude.sqrt + +instance Transcendental 𝕀 where + pi = 𝕀 3.141592653589793 3.1415926535897936 + cos = withHW Prelude.cos + sin = withHW Prelude.sin + atan = withHW Prelude.atan + +{-# INLINE withHW #-} +-- | Internal function: use @rounded-hw@ to define a function on intervals. +withHW :: (Interval.Interval Double -> Interval.Interval Double) -> 𝕀 -> 𝕀 +withHW f = \ ( 𝕀 lo hi ) -> + case f ( Interval.I ( Rounded lo ) ( Rounded hi ) ) of + Interval.I ( Rounded x ) ( Rounded y ) -> 𝕀 x y + +scaleInterval :: Double -> 𝕀 -> 𝕀 +scaleInterval s@(D# s#) i = + case compare s 0 of + LT -> case negate i of { MkI x -> MkI ( broadcastDoubleX2# ( negateDouble# s# ) `timesDoubleX2#` x ) } + EQ -> 𝕀 0 0 + GT -> case i of { MkI x -> MkI ( broadcastDoubleX2# s# `timesDoubleX2#` x ) } + +{- +TODO: consider using https://github.com/JishinMaster/simd_utils/blob/160c50f07e08d2077ae4368f0aed2f10f7173c67/simd_utils_sse_double.h#L530 +for sin/cos? Or Intel SVML? Or metalibm? + +-} diff --git a/brush-strokes/src/lib/Math/Root/Isolation.hs b/brush-strokes/src/lib/Math/Root/Isolation.hs index 16d6b2c..51a9c76 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation.hs @@ -204,7 +204,7 @@ isolateRootsIn ( RootIsolationOptions { rootIsolationAlgorithms } ) [ RootIsolationTree ( Box n ) ] go history cand | -- Check the range of the equations contains zero. - not $ ( unT ( origin @Double ) ∈ iRange ) + not $ ( all ( unT ( origin @Double ) ∈ ) $ coordinates iRange ) -- Box doesn't contain a solution: discard it. = return [ RootIsolationLeaf "rangeTest" cand ] | otherwise diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Bisection.hs b/brush-strokes/src/lib/Math/Root/Isolation/Bisection.hs index fc2db0a..de4bc67 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation/Bisection.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation/Bisection.hs @@ -120,7 +120,7 @@ defaultBisectionOptions minWidth _Ξ΅_eq box = -- box(0)-consistency let iRange' :: Box d iRange' = eqs box' `monIndex` zeroMonomial - in unT ( origin @Double ) ∈ iRange' + in all ( unT ( origin @Double ) ∈ ) ( coordinates iRange' ) -- box(1)-consistency --let box1Options = Box1Options _Ξ΅_eq ( toList $ universe @n ) ( toList $ universe @d ) @@ -173,33 +173,33 @@ type CoordBisectionData :: Nat -> Nat -> Type -> Type data CoordBisectionData n d r = CoordBisectionData { coordIndex :: !( Fin n ) - , coordInterval :: !( 𝕀 Double ) + , coordInterval :: !𝕀 , coordJacobianColumn :: !( 𝕀ℝ d ) , coordBisectionData :: !r } -deriving stock instance ( Show ( ℝ d ), Show r ) +deriving stock instance ( Show ( 𝕀ℝ d ), Show r ) => Show ( CoordBisectionData n d r ) -spread :: ( BoxCt n d, Representable Double ( ℝ d ) ) +spread :: ( BoxCt n d, Representable 𝕀 ( 𝕀ℝ d ) ) => CoordBisectionData n d r -> Double spread ( CoordBisectionData { coordInterval = cd, coordJacobianColumn = j_cd } ) = width cd * normVI j_cd {-# INLINEABLE spread #-} -normVI :: ( Applicative ( Vec d ), Representable Double ( ℝ d ) ) => 𝕀ℝ d -> Double -normVI ( 𝕀 los his ) = - sqrt $ sum ( nm1 <$> coordinates los <*> coordinates his ) +normVI :: ( Applicative ( Vec d ), Representable 𝕀 ( 𝕀ℝ d ) ) => 𝕀ℝ d -> Double +normVI box = + sqrt $ sum ( nm1 <$> coordinates box ) where - nm1 :: Double -> Double -> Double - nm1 lo hi = max ( abs lo ) ( abs hi ) Ring.^ 2 + nm1 :: 𝕀 -> Double + nm1 ( 𝕀 lo hi ) = max ( abs lo ) ( abs hi ) Ring.^ 2 {-# INLINEABLE normVI #-} -maxVI :: ( Applicative ( Vec d ), Representable Double ( ℝ d ) ) => 𝕀ℝ d -> Double -maxVI ( 𝕀 los his ) = - maximum ( maxAbs <$> coordinates los <*> coordinates his ) +maxVI :: ( Applicative ( Vec d ), Representable 𝕀 ( 𝕀ℝ d ) ) => 𝕀ℝ d -> Double +maxVI box = + maximum ( maxAbs <$> coordinates box ) where - maxAbs :: Double -> Double -> Double - maxAbs lo hi = max ( abs lo ) ( abs hi ) + maxAbs :: 𝕀 -> Double + maxAbs ( 𝕀 lo hi ) = max ( abs lo ) ( abs hi ) {-# INLINEABLE maxVI #-} -------------------------------------------------------------------------------- @@ -210,7 +210,7 @@ maxVI ( 𝕀 los his ) = -- dimension to bisect.) bisection :: forall n - . ( 1 <= n, KnownNat n, Representable Double ( ℝ n ) ) + . ( 1 <= n, KnownNat n, Representable 𝕀 ( 𝕀ℝ n ) ) => ( Box n -> Bool ) -- ^ how to check whether a box contains solutions -> ( forall r. NE.NonEmpty ( Fin n, r ) -> ( r, String ) ) @@ -250,7 +250,7 @@ bisection canHaveSols pickFallbackBisCoord box = -- | Bisect a box in the given coordinate dimension. bisectInCoord - :: Representable Double ( ℝ n ) + :: Representable 𝕀 ( 𝕀ℝ n ) => Box n -> Fin n -> ( Double, NE.NonEmpty ( Box n ) ) bisectInCoord box i = let z = box `index` i diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Core.hs b/brush-strokes/src/lib/Math/Root/Isolation/Core.hs index d436724..16d5195 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation/Core.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation/Core.hs @@ -97,17 +97,21 @@ type BoxCt n d = , Show ( 𝕀ℝ n ), Show ( ℝ n ) , Eq ( ℝ n ) , Representable Double ( ℝ n ) + , Representable 𝕀 ( 𝕀ℝ n ) , MonomialBasis ( D 1 ( ℝ n ) ) , Deg ( D 1 ( ℝ n ) ) ~ 1 , Vars ( D 1 ( ℝ n ) ) ~ n , Module Double ( T ( ℝ n ) ) - , Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) + , Module 𝕀 ( T ( 𝕀ℝ n ) ) + , NFData ( ℝ n ) , NFData ( 𝕀ℝ n ) , Ord ( ℝ d ) , Module Double ( T ( ℝ d ) ) , Representable Double ( ℝ d ) + , Representable 𝕀 ( 𝕀ℝ d ) , NFData ( ℝ d ) + , NFData ( 𝕀ℝ d ) ) -- | Boxes we are done with and will not continue processing. data DoneBoxes n = @@ -235,15 +239,15 @@ data RootIsolationTree d | RootIsolationStep RootIsolationStep [ ( d, RootIsolationTree d ) ] showRootIsolationTree - :: ( Representable Double ( ℝ n ), Show ( Box n ) ) + :: ( Representable 𝕀 ( 𝕀ℝ n ), Show ( Box n ) ) => Box n -> RootIsolationTree ( Box n ) -> Tree String showRootIsolationTree cand (RootIsolationLeaf why l) = Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ why ++ " " ++ show l) [] showRootIsolationTree cand (RootIsolationStep s ts) = Node (show cand ++ " abc " ++ showArea (boxArea cand) ++ " " ++ show s) $ map (\ (c,t) -> showRootIsolationTree c t) ts -boxArea :: Representable Double ( ℝ n ) => Box n -> Double -boxArea ( 𝕀 lo hi ) = - product ( ( \ l h -> abs ( h - l ) ) <$> coordinates lo <*> coordinates hi ) +{-# INLINEABLE boxArea #-} +boxArea :: Representable 𝕀 ( 𝕀ℝ n ) => Box n -> Double +boxArea box = product ( width <$> coordinates box ) showArea :: Double -> String showArea area = "(area " ++ showFFloat (Just 6) area "" ++ ")" diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Degree.hs b/brush-strokes/src/lib/Math/Root/Isolation/Degree.hs index 8546cc5..a440e28 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation/Degree.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation/Degree.hs @@ -31,14 +31,14 @@ topologicalDegree -> 𝕀ℝ 2 -- ^ box -> Maybe Int -topologicalDegree Ξ΅_bis f ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = +topologicalDegree Ξ΅_bis f ( 𝕀ℝ2 ( 𝕀 x_lo x_hi ) ( 𝕀 y_lo y_hi ) ) = topologicalDegreeFromContour $ smallArrayFromList $ concatMap ( tagFacet Ξ΅_bis f ) - [ ( Tag ( Fin 2 ) P, 𝕀 ( ℝ2 x_hi y_lo ) ( ℝ2 x_hi y_hi ) ) - , ( Tag ( Fin 1 ) N, 𝕀 ( ℝ2 x_lo y_hi ) ( ℝ2 x_hi y_hi ) ) - , ( Tag ( Fin 2 ) N, 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_lo y_hi ) ) - , ( Tag ( Fin 1 ) P, 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_lo ) ) ] + [ ( Tag ( Fin 2 ) P, 𝕀ℝ2 ( 𝕀 x_hi x_hi ) ( 𝕀 y_lo y_hi ) ) + , ( Tag ( Fin 1 ) N, 𝕀ℝ2 ( 𝕀 x_lo x_hi ) ( 𝕀 y_hi y_hi ) ) + , ( Tag ( Fin 2 ) N, 𝕀ℝ2 ( 𝕀 x_lo x_lo ) ( 𝕀 y_lo y_hi ) ) + , ( Tag ( Fin 1 ) P, 𝕀ℝ2 ( 𝕀 x_lo x_hi ) ( 𝕀 y_lo y_lo ) ) ] -- | Tag a facet with the sign of one of the equations that is non-zero -- on that facet. @@ -77,11 +77,13 @@ bisectFacet ( Tag i s ) x = bisectFacet NoTag _ = error "impossible: input always tagged" bisectInCoord :: Fin n -> 𝕀ℝ 2 -> [ 𝕀ℝ 2 ] -bisectInCoord ( Fin 1 ) ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = - [ 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_mid y_hi ), 𝕀 ( ℝ2 x_mid y_lo ) ( ℝ2 x_hi y_hi ) ] +bisectInCoord ( Fin 1 ) ( 𝕀ℝ2 ( 𝕀 x_lo x_hi ) ( 𝕀 y_lo y_hi ) ) = + [ 𝕀ℝ2 ( 𝕀 x_lo x_mid ) ( 𝕀 y_lo y_hi ) + , 𝕀ℝ2 ( 𝕀 x_mid x_hi ) ( 𝕀 y_lo y_hi ) ] where x_mid = 0.5 * ( x_lo + x_hi ) -bisectInCoord _ ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = - [ 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_mid ), 𝕀 ( ℝ2 x_lo y_mid ) ( ℝ2 x_hi y_hi ) ] +bisectInCoord _ ( 𝕀ℝ2 ( 𝕀 x_lo x_hi ) ( 𝕀 y_lo y_hi ) ) = + [ 𝕀ℝ2 ( 𝕀 x_lo x_hi ) ( 𝕀 y_lo y_mid ) + , 𝕀ℝ2 ( 𝕀 x_lo x_hi ) ( 𝕀 y_mid y_hi ) ] where y_mid = 0.5 * ( y_lo + y_hi ) -- | Compute the topological degree of \( f \colon \mathbb{IR}^n \to \mathbb{IR}^n \) diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Narrowing.hs b/brush-strokes/src/lib/Math/Root/Isolation/Narrowing.hs index 05c04fb..5611e47 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation/Narrowing.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation/Narrowing.hs @@ -156,11 +156,14 @@ defaultBox2Options minWidth Ξ΅_eq = -- See also -- "Presentation of a highly tuned multithreaded interval solver for underdetermined and well-determined nonlinear systems" makeBox1Consistent - :: ( KnownNat n, Representable Double ( ℝ n ) + :: ( KnownNat n + , Representable Double ( ℝ n ) + , Representable 𝕀 ( 𝕀ℝ n ) , MonomialBasis ( D 1 ( ℝ n ) ) , Deg ( D 1 ( ℝ n ) ) ~ 1 , Vars ( D 1 ( ℝ n ) ) ~ n , Representable Double ( ℝ d ) + , Representable 𝕀 ( 𝕀ℝ d ) ) => Box1Options n d -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) @@ -174,11 +177,14 @@ makeBox1Consistent box1Options eqs x = -- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems" makeBox2Consistent :: forall n d - . ( KnownNat n, Representable Double ( ℝ n ) + . ( KnownNat n + , Representable Double ( ℝ n ) + , Representable 𝕀 ( 𝕀ℝ n ) , MonomialBasis ( D 1 ( ℝ n ) ) , Deg ( D 1 ( ℝ n ) ) ~ 1 , Vars ( D 1 ( ℝ n ) ) ~ n , Representable Double ( ℝ d ) + , Representable 𝕀 ( 𝕀ℝ d ) ) => Box2Options n d -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) @@ -240,7 +246,7 @@ data NarrowingMethod narrowingMethods :: Double -> Double -> NarrowingMethod - -> [ ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) ) -> 𝕀 Double -> [ 𝕀 Double ] ] + -> [ ( 𝕀 -> ( 𝕀, 𝕀 ) ) -> 𝕀 -> [ 𝕀 ] ] narrowingMethods Ξ΅_eq Ξ΅_bis Kubica = [ leftNarrow Ξ΅_eq Ξ΅_bis, rightNarrow Ξ΅_eq Ξ΅_bis ] narrowingMethods Ξ΅_eq Ξ΅_bis ( AdaptiveShaving opts ) @@ -252,11 +258,14 @@ narrowingMethods _Ξ΅_eq Ξ΅_bis TwoSidedShaving allNarrowingOperators :: forall n d - . ( KnownNat n, Representable Double ( ℝ n ) + . ( KnownNat n + , Representable Double ( ℝ n ) + , Representable 𝕀 ( 𝕀ℝ n ) , MonomialBasis ( D 1 ( ℝ n ) ) , Deg ( D 1 ( ℝ n ) ) ~ 1 , Vars ( D 1 ( ℝ n ) ) ~ n , Representable Double ( ℝ d ) + , Representable 𝕀 ( 𝕀ℝ d ) ) => Box1Options n d -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) @@ -281,7 +290,7 @@ allNarrowingOperators ( Box1Options Ξ΅_eq Ξ΅_bis coordsToNarrow eqsToUse narrowi , eqnIndex <- eqsToUse ] where - ff' :: Fin n -> Fin d -> 𝕀ℝ n -> ( 𝕀 Double, 𝕀 Double ) + ff' :: Fin n -> Fin d -> 𝕀ℝ n -> ( 𝕀, 𝕀 ) ff' i d ts = let df = eqs ts f, f' :: 𝕀ℝ d @@ -301,9 +310,9 @@ allNarrowingOperators ( Box1Options Ξ΅_eq Ξ΅_bis coordsToNarrow eqsToUse narrowi -- by BartΕ‚omiej Jacek Kubica, 2017 leftNarrow :: Double -> Double - -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) ) - -> 𝕀 Double - -> [ 𝕀 Double ] + -> ( 𝕀 -> ( 𝕀, 𝕀 ) ) + -> 𝕀 + -> [ 𝕀 ] leftNarrow Ξ΅_eq Ξ΅_bis ff' = left_narrow where left_narrow ( 𝕀 x_inf x_sup ) = @@ -334,9 +343,9 @@ leftNarrow Ξ΅_eq Ξ΅_bis ff' = left_narrow -- TODO: de-duplicate with 'leftNarrow'? rightNarrow :: Double -> Double - -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) ) - -> 𝕀 Double - -> [ 𝕀 Double ] + -> ( 𝕀 -> ( 𝕀, 𝕀 ) ) + -> 𝕀 + -> [ 𝕀 ] rightNarrow Ξ΅_eq Ξ΅_bis ff' = right_narrow where right_narrow ( 𝕀 x_inf x_sup ) = @@ -388,23 +397,23 @@ defaultAdaptiveShavingOptions = leftShave :: Double -> Double -> AdaptiveShavingOptions - -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) ) - -> 𝕀 Double - -> [ 𝕀 Double ] + -> ( 𝕀 -> ( 𝕀, 𝕀 ) ) + -> 𝕀 + -> [ 𝕀 ] leftShave Ξ΅_eq Ξ΅_bis ( AdaptiveShavingOptions { Ξ³_init, Οƒ_good, Οƒ_bad, Ξ²_good, Ξ²_bad } ) ff' i0 = left_narrow Ξ³_init i0 where w0 = width i0 - left_narrow :: Double -> 𝕀 Double -> [ 𝕀 Double ] + left_narrow :: Double -> 𝕀 -> [ 𝕀 ] left_narrow Ξ³ i@( 𝕀 x_inf x_sup ) -- Stop if the box is too small. | width i < Ξ΅_bis = [ i ] | otherwise = go Ξ³ x_sup ( 𝕀 x_inf ( if x_inf == x_sup then x_inf else succFP x_inf ) ) - go :: Double -> Double -> 𝕀 Double -> [ 𝕀 Double ] + go :: Double -> Double -> 𝕀 -> [ 𝕀 ] go Ξ³ x_sup x_left = let ( f_x_left, _f'_x_left ) = ff' x_left in @@ -429,7 +438,7 @@ leftShave Ξ΅_eq Ξ΅_bis -- Do a Newton step to try to reduce the guess interval. -- Starting from "guess", we get a collection of sub-intervals -- "guesses'" refining where the function can be zero. - let guesses' :: [ ( 𝕀 Double ) ] + let guesses' :: [ 𝕀 ] guesses' = do Ξ΄ <- f_x_left ⊘ f'_guess let guess' = singleton ( inf guess ) - Ξ΄ @@ -463,11 +472,11 @@ leftShave Ξ΅_eq Ξ΅_bis -- "A Data-Parallel Algorithm to Reliably Solve Systems of Nonlinear Equations", -- (FrΓ©dΓ©ric Goualard, Alexandre Goldsztejn, 2008). sbc :: Double - -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) ) - -> 𝕀 Double -> [ 𝕀 Double ] + -> ( 𝕀 -> ( 𝕀, 𝕀 ) ) + -> 𝕀 -> [ 𝕀 ] sbc Ξ΅ ff' = go where - go :: 𝕀 Double -> [ 𝕀 Double ] + go :: 𝕀 -> [ 𝕀 ] go x@( 𝕀 x_l x_r ) | width x <= Ξ΅ = [ x ] diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Newton.hs b/brush-strokes/src/lib/Math/Root/Isolation/Newton.hs index acd6760..b14501a 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation/Newton.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation/Newton.hs @@ -86,6 +86,8 @@ defaultNewtonOptions , 1 <= n, 1 <= d, n <= d , Representable Double ( ℝ n ) , Representable Double ( ℝ d ) + , Representable 𝕀 ( 𝕀ℝ n ) + , Representable 𝕀 ( 𝕀ℝ d ) ) => BoxHistory n -> NewtonOptions n d @@ -110,7 +112,7 @@ intervalNewton opts eqs x = case opts of , gsPickEqs = pickEqs , gsUpdate } ) -> - let x_mid = singleton $ boxMidpoint x + let x_mid = point $ boxMidpoint x f :: 𝕀ℝ n -> 𝕀ℝ n f = \ x_0 -> pickEqs $ eqs x_0 `monIndex` zeroMonomial f'_x :: Vec n ( 𝕀ℝ n ) @@ -123,7 +125,7 @@ intervalNewton opts eqs x = case opts of -- Precondition the above linear system into A ( x - x_mid ) = B. !( !a, !b ) = force $ precondition precondMeth - f'_x ( singleton minus_f_x_mid ) + f'_x ( point minus_f_x_mid ) !x'_0 = force ( T x ^-^ T x_mid ) -- NB: we have to change coordinates, putting the midpoint of the box @@ -144,14 +146,14 @@ intervalNewton opts eqs x = case opts of return ( dt, todo ) NewtonLP -> -- TODO: reduce duplication with the above. - let x_mid = singleton $ boxMidpoint x + let x_mid = point $ boxMidpoint x f :: 𝕀ℝ 2 -> 𝕀ℝ d f = \ x_0 -> eqs x_0 `monIndex` zeroMonomial f'_x :: Vec 2 ( 𝕀ℝ d ) f'_x = fmap ( \ i -> eqs x `monIndex` linearMonomial i ) ( universe @2 ) minus_f_x_mid = unT $ -1 *^ T ( boxMidpoint $ f x_mid ) - !( !a, !b ) = force ( f'_x, singleton minus_f_x_mid ) + !( !a, !b ) = force ( f'_x, point minus_f_x_mid ) !x'_0 = force ( T x ^-^ T x_mid ) ( x's, dt ) = timeInterval diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Newton/GaussSeidel.hs b/brush-strokes/src/lib/Math/Root/Isolation/Newton/GaussSeidel.hs index 99e8dd3..194d2d5 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation/Newton/GaussSeidel.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation/Newton/GaussSeidel.hs @@ -68,8 +68,8 @@ defaultGaussSeidelOptions :: forall n d . ( KnownNat n, KnownNat d , 1 <= n, 1 <= d, n <= d - , Representable Double ( ℝ n ) - , Representable Double ( ℝ d ) + , Representable Double ( ℝ n ), Representable 𝕀 ( 𝕀ℝ n ) + , Representable Double ( ℝ d ), Representable 𝕀 ( 𝕀ℝ d ) ) => BoxHistory n -> GaussSeidelOptions n d @@ -102,7 +102,10 @@ data Preconditioner -- refining the initial guess box for \( X \) into up to \( 2^n \) (disjoint) new boxes. gaussSeidelUpdate :: forall n - . ( Representable Double ( ℝ n ), n ~ RepDim ( ℝ n ), Eq ( ℝ n ) ) + . ( Representable Double ( ℝ n ) + , Representable 𝕀 ( 𝕀ℝ n ) + , Eq ( ℝ n ) + ) => GaussSeidelUpdateMethod -- ^ which step method to use -> Vec n ( 𝕀ℝ n ) -- ^ columns of \( A \) -> 𝕀ℝ n -- ^ \( B \) @@ -120,7 +123,10 @@ gaussSeidelUpdate upd as b x = -- The boolean indicates whether the Gauss–Seidel step was a contraction. gaussSeidelStep :: forall n - . ( Representable Double ( ℝ n ), n ~ RepDim ( ℝ n ), Eq ( ℝ n ) ) + . ( Representable Double ( ℝ n ) + , Representable 𝕀 ( 𝕀ℝ n ) + , Eq ( ℝ n ) + ) => Vec n ( 𝕀ℝ n ) -- ^ columns of \( A \) -> 𝕀ℝ n -- ^ \( B \) -> T ( 𝕀ℝ n ) -- ^ initial box \( X \) @@ -160,7 +166,9 @@ gaussSeidelStep as b ( T x0 ) = coerce $ -- (Montanher, Domes, Schichl, Neumaier) (2017) gaussSeidelStep_Complete :: forall n - . ( Representable Double ( ℝ n ), n ~ RepDim ( ℝ n ) ) + . ( Representable Double ( ℝ n ) + , Representable 𝕀 ( 𝕀ℝ n ) + ) => Vec n ( 𝕀ℝ n ) -- ^ columns of \( A \) -> 𝕀ℝ n -- ^ \( B \) -> T ( 𝕀ℝ n ) -- ^ initial box \( X \) @@ -193,7 +201,10 @@ gaussSeidelStep_Complete as b ( T x0 ) = coerce $ do -- | Pre-condition the system \( AX = B \). precondition :: forall n - . ( KnownNat n, Representable Double ( ℝ n ) ) + . ( KnownNat n + , Representable Double ( ℝ n ) + , Representable 𝕀 ( 𝕀ℝ n ) + ) => Preconditioner -- ^ pre-conditioning method to use -> Vec n ( 𝕀ℝ n ) -- ^ columns of \( A \) -> 𝕀ℝ n -- ^ \( B \) diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Newton/LP.hs b/brush-strokes/src/lib/Math/Root/Isolation/Newton/LP.hs index 61419a1..76c0289 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation/Newton/LP.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation/Newton/LP.hs @@ -42,7 +42,7 @@ import Math.Linear -- | Solve the system of linear equations \( A X = B \) -- using linear programming. solveIntervalLinearEquations - :: ( KnownNat d, Representable Double ( ℝ d ) ) + :: ( KnownNat d, Representable 𝕀 ( 𝕀ℝ d ) ) => Vec 2 ( 𝕀ℝ d ) -- ^ columns of \( A \) -> 𝕀ℝ d -- ^ \( B \) -> T ( 𝕀ℝ 2 ) -- ^ initial box \( X \) @@ -61,18 +61,18 @@ solveIntervalLinearEquations a b x = -- is in fact a strict subset. isSubBox :: T ( 𝕀ℝ 2 ) -> T ( 𝕀ℝ 2 ) -> Bool isSubBox - ( T ( 𝕀 ( ℝ2 x1_min y1_min ) ( ℝ2 x1_max y1_max ) ) ) - ( T ( 𝕀 ( ℝ2 x2_min y2_min ) ( ℝ2 x2_max y2_max ) ) ) + ( T ( 𝕀ℝ2 ( 𝕀 x1_min x1_max ) ( 𝕀 y1_min y1_max ) ) ) + ( T ( 𝕀ℝ2 ( 𝕀 x2_min x2_max ) ( 𝕀 y2_min y2_max ) ) ) = ( ( x2_min > x1_min && x2_max < x1_max ) || x2_min == x2_max ) && ( ( y2_min > y1_min && y2_max < y1_max ) || y2_min == y2_max ) hasNaNs :: T (𝕀ℝ 2) -> Bool -hasNaNs ( T ( 𝕀 ( ℝ2 x_min y_min ) ( ℝ2 x_max y_max ) ) ) = +hasNaNs ( T ( 𝕀ℝ2 ( 𝕀 x_min x_max ) ( 𝕀 y_min y_max ) ) ) = any ( \ x -> isNaN x || isInfinite x ) [ x_min, y_min, x_max, y_max ] intervalSystem2D_LP :: forall d - . ( KnownNat d, Representable Double ( ℝ d ) ) + . ( KnownNat d, Representable 𝕀 ( 𝕀ℝ d ) ) => Vec 2 ( 𝕀ℝ d ) -> 𝕀ℝ d -> T ( 𝕀ℝ 2 ) @@ -94,7 +94,7 @@ intervalSystem2D_LP a b x = d = natVal' @d proxy# mkEquationArray - :: ( KnownNat d, Representable Double ( ℝ d ) ) + :: ( KnownNat d, Representable 𝕀 ( 𝕀ℝ d ) ) => Vec 2 ( 𝕀ℝ d ) -> 𝕀ℝ d -> [ CEqn ] mkEquationArray ( Vec [ a_x, a_y ] ) b = [ CEqn a_x_i a_y_i b_i @@ -107,7 +107,7 @@ mkEquationArray _ _ = error "impossible" foreign import ccall "interval_system_2d" interval_system_2d :: Ptr CBox -> Ptr CBox -> Ptr CEqn -> CUInt -> IO CInt -data CEqn = CEqn !( 𝕀 Double ) !( 𝕀 Double ) !( 𝕀 Double ) +data CEqn = CEqn !𝕀 !𝕀 !𝕀 instance Storable CEqn where sizeOf _ = 6 * sizeOf @Double undefined alignment _ = 4 * alignment @Double undefined @@ -127,6 +127,6 @@ instance Storable CBox where peek ptr = do [ CDouble x_min, CDouble x_max, CDouble y_min, CDouble y_max ] <- peekArray 4 ( castPtr ptr :: Ptr CDouble ) return $ - CBox ( 𝕀 ( ℝ2 x_min y_min ) ( ℝ2 x_max y_max ) ) - poke ptr ( CBox ( 𝕀 ( ℝ2 x_min y_min ) ( ℝ2 x_max y_max ) ) ) = + CBox ( 𝕀ℝ2 ( 𝕀 x_min x_max ) ( 𝕀 y_min y_max ) ) + poke ptr ( CBox ( 𝕀ℝ2 ( 𝕀 x_min x_max ) ( 𝕀 y_min y_max ) ) ) = pokeArray ( castPtr ptr ) [ CDouble x_min, CDouble x_max, CDouble y_min, CDouble y_max ] diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Utils.hs b/brush-strokes/src/lib/Math/Root/Isolation/Utils.hs index f273bcf..4a6e931 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation/Utils.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation/Utils.hs @@ -20,8 +20,8 @@ import Math.Linear fromComponents :: forall n - . ( Representable Double ( ℝ n ), n ~ RepDim ( ℝ n ) ) - => ( Fin n -> [ ( 𝕀 Double, Bool ) ] ) -> [ ( 𝕀ℝ n, Vec n Bool ) ] + . ( Representable 𝕀 ( 𝕀ℝ n ), n ~ RepDim ( 𝕀ℝ n ) ) + => ( Fin n -> [ ( 𝕀, Bool ) ] ) -> [ ( 𝕀ℝ n, Vec n Bool ) ] fromComponents f = do ( xs, bs ) <- unzip <$> traverse f ( universe @n ) return $ ( tabulate $ \ i -> xs ! i, bs ) @@ -29,7 +29,7 @@ fromComponents f = do {-# INLINEABLE fromComponents #-} -- | The midpoint of a box. -boxMidpoint :: Representable Double ( ℝ n ) => 𝕀ℝ n -> ℝ n +boxMidpoint :: ( Representable 𝕀 ( 𝕀ℝ n ), Representable Double ( ℝ n ) ) => 𝕀ℝ n -> ℝ n boxMidpoint box = tabulate $ \ i -> midpoint ( box `index` i ) {-# INLINEABLE boxMidpoint #-} @@ -37,7 +37,9 @@ boxMidpoint box = -- | Matrix multiplication \( A v \). matMulVec :: forall n m - . ( Representable Double ( ℝ n ), Representable Double ( ℝ m ) ) + . ( Representable 𝕀 ( 𝕀ℝ n ), Representable 𝕀 ( 𝕀ℝ m ) + , Representable Double ( ℝ n ), Representable Double ( ℝ m ) + ) => Vec m ( ℝ n ) -- ^ columns of the matrix \( A ) -> 𝕀ℝ m -- ^ vector \( v \) -> 𝕀ℝ n diff --git a/cabal.project b/cabal.project index 38c14c5..6cfca0a 100644 --- a/cabal.project +++ b/cabal.project @@ -50,3 +50,18 @@ source-repository-package type: git location: https://github.com/sheaf/eigen tag: 3bc1d4bc53edb75e2ee3893763361aa0d5c7714a + +-------------- +-- GHC HEAD -- +-------------- + +repository head.hackage.ghc.haskell.org + url: https://ghc.gitlab.haskell.org/head.hackage/ + secure: True + key-threshold: 3 + root-keys: + f76d08be13e9a61a377a85e2fb63f4c5435d40f8feb3e12eb05905edb8cdea89 + 26021a13b401500c8eb2761ca95c61f2d625bfef951b939a8124ed12ecf07329 + 7541f32a4ccca4f97aea3b22f5e593ba2c0267546016b992dfadcd2fe944e55d + +active-repositories: hackage.haskell.org, head.hackage.ghc.haskell.org:override