hashirama/metabrush
1
0
Fork 0
mirror of https://gitlab.com/sheaf/metabrush.git synced 2024-11-23 15:34:06 +00:00
Code Issues Actions Projects Releases Wiki Activity

WIP: SIMD experiments

Browse source
This commit is contained in:
sheaf 2024-08-01 21:36:44 +02:00
parent e64ce375d2
commit 5f98165276
28 changed files with 1252 additions and 854 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

33
brush-strokes/brush-strokes.cabal
View file

@ -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

15
brush-strokes/cabal.project
View file

@ -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

2
brush-strokes/src/cusps/bench/Bench/Cases.hs
View file

@ -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

45
brush-strokes/src/cusps/bench/Bench/Types.hs
View file

@ -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

18
brush-strokes/src/cusps/bench/Main.hs
View file

@ -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

28
brush-strokes/src/cusps/inspect/Main.hs
View file

@ -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

6
brush-strokes/src/cusps/inspect/Math/Interval/Abstract.hs
View file

@ -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

123
brush-strokes/src/lib/Calligraphy/Brushes.hs
View file

@ -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 #-}

64
brush-strokes/src/lib/Math/Algebra/Dual.hs
View file

@ -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.

316
brush-strokes/src/lib/Math/Bezier/Stroke.hs
View file

@ -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
= [ ]

102
brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs
View file

@ -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

22
brush-strokes/src/lib/Math/Differentiable.hs
View file

@ -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 )

310
brush-strokes/src/lib/Math/Interval.hs
View file

@ -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 ||] )

92
brush-strokes/src/lib/Math/Interval/FMA.hs
View file

@ -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 )

297
brush-strokes/src/lib/Math/Interval/Internal.hs
View file

@ -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 #-}

214
brush-strokes/src/lib/Math/Interval/Internal/FMA.hs Normal file
View file

@ -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 )

53
brush-strokes/src/lib/Math/Interval/Internal/RoundedHW.hs Normal file
View file

@ -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

173
brush-strokes/src/lib/Math/Interval/Internal/SIMD.hs Normal file
View file

@ -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?
-}

2
brush-strokes/src/lib/Math/Root/Isolation.hs
View file

@ -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

32
brush-strokes/src/lib/Math/Root/Isolation/Bisection.hs
View file

@ -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

14
brush-strokes/src/lib/Math/Root/Isolation/Core.hs
View file

@ -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 "" ++ ")"

20
brush-strokes/src/lib/Math/Root/Isolation/Degree.hs
View file

@ -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 \)

49
brush-strokes/src/lib/Math/Root/Isolation/Narrowing.hs
View file

@ -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 ]

10
brush-strokes/src/lib/Math/Root/Isolation/Newton.hs
View file

@ -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

23
brush-strokes/src/lib/Math/Root/Isolation/Newton/GaussSeidel.hs
View file

@ -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 \)

18
brush-strokes/src/lib/Math/Root/Isolation/Newton/LP.hs
View file

@ -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 ]

10
brush-strokes/src/lib/Math/Root/Isolation/Utils.hs
View file

@ -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

15
cabal.project
View file

@ -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