Split out benchmark for cusp finding

brush-strokes/bench/Main.hs

@@ -0,0 +1,272 @@
+{-# LANGUAGE PolyKinds #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+module Main
+  ( main
+
+    -- Testing
+  , TestCase(..)
+  , testCases
+  , testCaseStrokeFunctions
+  , eval
+  , mkVal, mkBox
+  , potentialCusp
+  , dEdsdcdt
+  )
+  where
+
+-- base
+import Data.Coerce
+  ( coerce )
+import Data.Foldable
+  ( for_ )
+import GHC.Exts
+  ( Proxy#, proxy# )
+import GHC.Generics
+  ( Generic )
+import GHC.TypeNats
+  ( type (-) )
+
+-- containers
+import Data.Sequence
+  ( Seq )
+import qualified Data.Sequence as Seq
+  ( index )
+import Data.Tree
+  ( foldTree )
+
+-- brush-strokes
+import Calligraphy.Brushes
+import Math.Algebra.Dual
+import Math.Bezier.Spline
+import Math.Bezier.Stroke
+import Math.Bezier.Stroke.EnvelopeEquation
+import Math.Differentiable
+import Math.Interval
+import Math.Linear
+import Math.Module
+import Math.Ring
+  ( Transcendental )
+
+--------------------------------------------------------------------------------
+
+main :: IO ()
+main = for_ testCases $ \ testCase@( TestCase { testName, testAlgorithmParams } ) -> do
+  let ( _, testStrokeFnI ) = testCaseStrokeFunctions testCase
+      ( newtTrees, ( dunno, sols ) ) = computeCusps testAlgorithmParams testStrokeFnI
+      showedTrees = map ( uncurry showIntervalNewtonTree ) newtTrees
+  putStrLn $ unlines $
+    [ ""
+    , "Test case '" ++ testName ++ "':" ] ++
+    map ( "    " ++ )
+      [ " #sols: " ++ show (length sols)
+      , "#dunno: " ++ show (length dunno)
+      , "#trees: " ++ show @Int (sum @_ @Int $ map (foldTree ( \ _ bs -> 1 + sum bs )) showedTrees)
+      , " dunno: " ++ show dunno
+      , "  sols: " ++ show sols
+      ]
+
+testCases :: [ TestCase ]
+testCases = [ ellipse , trickyCusp2 ]
+
+--------------------------------------------------------------------------------
+
+data TestCase =
+  forall nbParams. ParamsCt nbParams =>
+  TestCase
+  { testName   :: !String
+  , testBrush  :: !( Brush nbParams )
+  , testStroke :: !( Point nbParams, Curve Open () ( Point nbParams ))
+  , testAlgorithmParams :: !CuspAlgorithmParams
+  }
+
+testCaseStrokeFunctions
+  :: TestCase
+  -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () )
+     , 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
+testCaseStrokeFunctions ( TestCase { testStroke = ( sp0, crv ), testBrush } ) =
+  getStrokeFunctions testBrush sp0 crv
+
+-- Utilities to use in GHCi to help debugging.
+
+eval
+  :: ( 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 )
+
+mkBox :: ( Double, Double ) -> Int -> ( Double, Double ) -> Box
+mkBox ( t_min, t_max ) i ( s_min, s_max ) =
+  ( 𝕀 ( ℝ1 t_min ) ( ℝ1 t_max ) , i, 𝕀 ( ℝ1 s_min ) ( ℝ1 s_max ) )
+
+potentialCusp :: StrokeDatum 3 𝕀 -> Bool
+potentialCusp
+  ( StrokeDatum
+    { ee = D22 { _D22_v = 𝕀 ( ℝ1 ee_min ) ( ℝ1 ee_max ) }
+    , 𝛿E𝛿sdcdt = D12 { _D12_v = T ( 𝕀 ( ℝ2 vx_min vy_min ) ( ℝ2 vx_max vy_max ) )}
+    }
+  ) = ee_min <= 0 && ee_max >= 0
+    && 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 { 𝛿E𝛿sdcdt = v } ) = v
+
+{-
+let (f, fI) = testCaseStrokeFunctions trickyCusp2
+
+take 10 $ Data.List.sortOn ( \ ( _, ℝ1 e, v) -> abs e + norm v ) [ let { v = mkVal x 3 y; d = eval f v } in ( v, _D12_v $ ee d, _D0_v $ dEdsdcdt d ) | x <- [0.57,0.5701 .. 0.58], y <- [0.29,0.291..0.3] ]
+> ((ℝ1 0.5798800000000057,3,ℝ1 0.267980000000008),ℝ1 -2.8596965543670194e-4,V2 7.79559474412963e-2 2.0389671921293484e-2)
+
+potentialCusp $ eval fI $ mkBox (0.5798, 0.5799) 3 (0.26798, 0.26799)
+> True
+
+let nbPotentialSols b = let ( _newtTrees, ( dunno, sols ) ) = intervalNewtonGSFrom NoPreconditioning 1e-7 fI b in length dunno + length sols
+
+nbPotentialSols $ mkBox (0.5798, 0.5799) 3 (0.26798, 0.26799)
+1
+
+nbPotentialSols $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799)
+0
+
+let showTrees b = map ( uncurry showIntervalNewtonTree ) $ fst $ intervalNewtonGSFrom NoPreconditioning 1e-7 fI b
+
+putStrLn $ unlines $ map Data.Tree.View.showTree $ showTrees $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799)
+
+([ℝ1 0.5798, ℝ1 0.675],3,[ℝ1 0.26798, ℝ1 0.26799]) (area 0.000001) N []
+ └─ ([ℝ1 0.5973000285624527, ℝ1 0.6750000000000002],3,[ℝ1 0.26798, ℝ1 0.26799000000000006]) (area 0.000001) NoSolution "ee" ([ℝ1 0.5973000285624527, ℝ1 0.6750000000000002],3,[ℝ1 0.26798, ℝ1 0.26799000000000006])
+
+
+-}
+
+--------------------------------------------------------------------------------
+
+ellipse :: TestCase
+ellipse =
+  TestCase
+    { testName = "ellipse"
+    , testBrush = ellipseBrush
+    , testStroke = ( p0, LineTo ( NextPoint p1 ) () )
+    , testAlgorithmParams =
+        CuspAlgorithmParams
+         { preconditioning = NoPreconditioning
+         , maxWidth = 1e-7
+         }
+    }
+  where
+    mkPt x y w h phi =
+      Point
+        { pointCoords = ℝ2 x y
+        , pointParams = Params $ ℝ3 w h phi
+        }
+    p0 = mkPt 0 0 10 25 0
+    p1 = mkPt 100 0 15 40 pi
+
+trickyCusp2 :: TestCase
+trickyCusp2 =
+  TestCase
+    { testName = "trickyCusp2"
+    , testBrush = circleBrush
+    , testStroke = ( p0, Bezier3To p1 p2 ( NextPoint p3 ) () )
+    , testAlgorithmParams =
+        CuspAlgorithmParams
+         { preconditioning = NoPreconditioning
+         , maxWidth = 1e-7
+         }
+    }
+  where
+    mkPt x y =
+      Point
+        { pointCoords = ℝ2 x y
+        , pointParams = Params $ ℝ1 5.0
+        }
+    p0 = mkPt 5e+1 -5e+1
+    p1 = mkPt -7.72994362904069e+1 -3.124468786098509e+1
+    p2 = mkPt -5.1505430313958364e+1 -3.9826386521527986e+1
+    p3 = mkPt -5e+1 -5e+1
+
+--------------------------------------------------------------------------------
+
+type ParamsCt nbParams
+  = ( Show ( ℝ nbParams )
+    , HasChainRule Double 2 ( ℝ nbParams )
+    , HasChainRule ( 𝕀 Double ) 3 ( 𝕀 ( ℝ nbParams ) )
+    , Applicative ( D 2 ( ℝ nbParams ) )
+    , Applicative ( D 3 ( ℝ nbParams ) )
+    , Traversable ( D 2 ( ℝ nbParams ) )
+    , Traversable ( D 3 ( ℝ nbParams ) )
+    , Representable Double ( ℝ nbParams )
+    , Module Double ( T ( ℝ nbParams ) )
+    , Module ( 𝕀 Double ) ( T ( 𝕀 ( ℝ nbParams ) ) )
+    , Module ( D 2 ( ℝ nbParams ) Double ) ( D 2 ( ℝ nbParams ) ( ℝ 2 ) )
+    , Module ( D 3 ( ℝ nbParams ) ( 𝕀 Double ) ) ( D 3 ( ℝ nbParams ) ( 𝕀 ( ℝ 2 ) ) )
+    , Transcendental ( D 2 ( ℝ nbParams ) Double )
+    , Transcendental ( D 3 ( ℝ nbParams ) ( 𝕀 Double ) )
+    )
+
+newtype Params nbParams = Params { getParams :: ( ℝ nbParams ) }
+deriving newtype instance Show ( ℝ nbParams ) => Show ( Params nbParams )
+
+data Point nbParams =
+  Point
+    { pointCoords :: !( ℝ 2 )
+    , pointParams :: !( Params nbParams ) }
+  deriving stock Generic
+deriving stock instance Show ( ℝ nbParams ) => Show ( Point nbParams )
+
+data CuspAlgorithmParams =
+  CuspAlgorithmParams
+    { preconditioning :: !Preconditioner
+    , maxWidth :: !Double
+    }
+  deriving stock Show
+
+type Brush nbParams
+  = forall {t} k (i :: t)
+  .  DiffInterp k i ( ℝ nbParams )
+  => Proxy# i
+  -> ( forall a. a -> I i a )
+  -> C k ( I i ( ℝ nbParams ) )
+         ( Spline Closed () ( I i ( ℝ 2 ) ) )
+
+getStrokeFunctions
+  :: forall nbParams
+  .  ParamsCt nbParams
+  => Brush nbParams
+      -- ^ brush shape
+  -> Point nbParams
+      -- ^ start point
+  -> Curve Open () ( Point nbParams )
+      -- ^ curve points
+  -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () )
+     , 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
+getStrokeFunctions brush sp0 crv =
+  let
+    usedParams :: C 2 ( ℝ 1 ) ( ℝ nbParams )
+    path :: C 2 ( ℝ 1 ) ( ℝ 2 )
+    ( path, usedParams ) =
+      pathAndUsedParams @2 @() coerce id ( getParams . pointParams )
+        sp0 crv
+    usedParamsI :: C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ nbParams )
+    pathI :: C 3 ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 )
+    ( pathI, usedParamsI ) =
+      pathAndUsedParams @3 @𝕀 coerce singleton ( getParams . pointParams )
+        sp0 crv
+  in ( brushStrokeData @2 @( ℝ nbParams ) coerce coerce
+        path usedParams $
+        brush @2 @() proxy# id
+     , brushStrokeData @3 @( ℝ nbParams ) coerce coerce
+        pathI usedParamsI $
+        brush @3 @𝕀 proxy# singleton )
+{-# INLINEABLE getStrokeFunctions #-}
+
+computeCusps
+  :: CuspAlgorithmParams
+  -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
+  -> ( [ ( Box, IntervalNewtonTree Box ) ], ( [ Box ], [ Box ] ) )
+computeCusps params =
+  intervalNewtonGS ( preconditioning params ) ( maxWidth params )

brush-strokes/brush-strokes.cabal

@@ -19,32 +19,8 @@ common common
     build-depends:
         base
            >= 4.17    && < 4.20
-      , acts
-          ^>= 0.3.1.0
-      , code-page
-          ^>= 0.2.1
       , containers
            >= 0.6.0.1 && < 0.8
-      , deepseq
-           >= 1.4.4.0 && < 1.6
-      , directory
-           >= 1.3.7.1 && < 1.4
-      , filepath
-           >= 1.4     && < 1.6
-      , generic-lens
-           >= 2.2     && < 2.3
-      , groups
-          ^>= 0.5.3
-      , groups-generic
-          ^>= 0.3.1.0
-      , primitive
-          ^>= 0.9.0.0
-      , rounded-hw
-          ^>= 0.4
-      , time
-          ^>= 1.12.2
-      , transformers
-           >= 0.5.6.2 && < 0.7
       , tree-view
           ^>= 0.5
 
@@ -119,7 +95,8 @@ library
         Haskell2010
 
     exposed-modules:
-        Math.Algebra.Dual
+        Calligraphy.Brushes
+      , Math.Algebra.Dual
       , Math.Bezier.Cubic
       , Math.Bezier.Cubic.Fit
       , Math.Bezier.Quadratic
@@ -150,11 +127,56 @@ library
         Math.Interval.FMA
 
     build-depends:
-        bifunctors
+        template-haskell
+           >= 2.18       && < 2.22
+
+      , acts
+          ^>= 0.3.1.0
+      , bifunctors
            >= 5.5.4      && < 5.7
+      , code-page
+          ^>= 0.2.1
+      , deepseq
+           >= 1.4.4.0 && < 1.6
+      , directory
+           >= 1.3.7.1 && < 1.4
       , eigen
           ^>= 3.3.7.0
+      , filepath
+           >= 1.4     && < 1.6
+      , generic-lens
+           >= 2.2     && < 2.3
+      , groups
+          ^>= 0.5.3
+      , groups-generic
+          ^>= 0.3.1.0
       , parallel
           ^>= 3.2.2.0
-      , template-haskell
+      , primitive
-           >= 2.18       && < 2.22
+          ^>= 0.9.0.0
+      , rounded-hw
+          ^>= 0.4
+      , time
+          ^>= 1.12.2
+      , transformers
+           >= 0.5.6.2 && < 0.7
+
+benchmark cusps
+
+  import:
+    common
+
+  hs-source-dirs:
+    bench
+
+  main-is:
+    Main.hs
+
+  build-depends:
+    brush-strokes
+
+  type:
+    exitcode-stdio-1.0
+
+  default-language:
+    Haskell2010

brush-strokes/cabal.project

@@ -9,6 +9,10 @@ allow-newer:
   groups-generic:base,
   eigen:primitive,
 
+--package brush-strokes
+profiling: True
+profiling-detail: late
+
 -------------
 -- GHC 9.4 --
 -------------

brush-strokes/src/Calligraphy/Brushes.hs

@@ -0,0 +1,215 @@
+
+{-# LANGUAGE AllowAmbiguousTypes #-}
+{-# LANGUAGE OverloadedStrings   #-}
+{-# LANGUAGE PolyKinds           #-}
+{-# LANGUAGE RebindableSyntax    #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+module Calligraphy.Brushes
+  ( circleBrush
+  , ellipseBrush
+  , tearDropBrush
+  ) where
+
+-- base
+import Prelude
+  hiding
+    ( Num(..), Floating(..), (^), (/), fromInteger, fromRational )
+import GHC.Exts
+  ( Proxy# )
+import GHC.TypeNats
+  ( type (<=) )
+
+-- containers
+import qualified Data.Sequence as Seq
+  ( empty, fromList )
+
+-- brush-strokes
+import Math.Algebra.Dual
+import Math.Bezier.Spline
+import Math.Differentiable
+  ( I )
+import Math.Linear
+import Math.Module
+  ( Module((^+^), (*^)) )
+import Math.Ring
+
+--------------------------------------------------------------------------------
+-- Circle & ellipse
+
+-- | Root of @(Sqrt[2] (4 + 3 κ) - 16) (2 - 3 κ)^2 - 8 (1 - 3 κ) Sqrt[8 - 24 κ + 12 κ^2 + 8 κ^3 + 3 κ^4]@.
+--
+-- Used to approximate circles and ellipses with Bézier curves.
+κ :: Double
+κ = 0.5519150244935105707435627227925
+
+circleSpline :: forall d v
+             .  Applicative d
+             => ( Double -> Double -> d v )
+             -> d ( Spline 'Closed () v )
+circleSpline p = sequenceA $
+  Spline { splineStart  = p 1 0
+         , splineCurves = ClosedCurves crvs lastCrv }
+  where
+    crvs = Seq.fromList
+      [ Bezier3To ( p  1  κ ) ( p  κ  1 ) ( NextPoint ( p  0  1 ) ) ()
+      , Bezier3To ( p -κ  1 ) ( p -1  κ ) ( NextPoint ( p -1  0 ) ) ()
+      , Bezier3To ( p -1 -κ ) ( p -κ -1 ) ( NextPoint ( p  0 -1 ) ) ()
+      ]
+    lastCrv =
+        Bezier3To ( p  κ -1 ) ( p  1 -κ ) BackToStart ()
+{-# INLINE circleSpline #-}
+
+circleBrush :: forall {t} (i :: t) k irec
+            . ( 1 <= RepDim irec
+              , Module
+                 ( D k irec ( I i Double ) )
+                 ( D k irec ( I i ( ℝ 2 ) ) )
+              , Module ( I i Double ) ( T ( I i Double ) )
+              , HasChainRule ( I i Double ) k irec
+              , Representable ( I i Double ) irec
+              , Applicative ( D k irec )
+              )
+            => Proxy# i
+            -> ( forall a. a -> I i a )
+            -> C k irec ( Spline 'Closed () ( I i ( ℝ 2 ) ) )
+circleBrush _ mkI =
+  D \ params ->
+    let r :: D k irec ( I i Double )
+        r = runD ( var @_ @k $ Fin 1 ) params
+        mkPt :: Double -> Double -> D k irec ( I i ( ℝ 2 ) )
+        mkPt x y
+          =   ( r `scaledBy` x ) *^ e_x
+          ^+^ ( r `scaledBy` y ) *^ e_y
+    in circleSpline mkPt
+  where
+    e_x, e_y :: D k irec ( I i ( ℝ 2 ) )
+    e_x = pure $ mkI $ ℝ2 1 0
+    e_y = pure $ mkI $ ℝ2 0 1
+
+    scaledBy d x = fmap ( mkI x * ) d
+{-# INLINEABLE circleBrush #-}
+
+ellipseBrush :: forall {t} (i :: t) k irec
+             . ( 3 <= RepDim irec
+               , Module
+                  ( D k irec ( I i Double ) )
+                  ( D k irec ( I i ( ℝ 2 ) ) )
+               , Module ( I i Double ) ( T ( I i Double ) )
+               , HasChainRule ( I i Double ) k irec
+               , Representable ( I i Double ) irec
+               , Applicative ( D k irec )
+               , Transcendental ( D k irec ( I i Double ) )
+                 -- TODO: make a synonym for the above...
+                 -- it seems DiffInterp isn't exactly right
+               )
+              => Proxy# i
+              -> ( forall a. a -> I i a )
+              -> C k irec ( Spline 'Closed () ( I i ( ℝ 2 ) ) )
+ellipseBrush _ mkI =
+  D \ params ->
+    let a, b, phi :: D k irec ( I i Double )
+        a   = runD ( var @_ @k $ Fin 1 ) params
+        b   = runD ( var @_ @k $ Fin 2 ) params
+        phi = runD ( var @_ @k $ Fin 3 ) params
+        mkPt :: Double -> Double -> D k irec ( I i ( ℝ 2 ) )
+        mkPt x y
+          = let !x' = a `scaledBy` x
+                !y' = b `scaledBy` y
+-- {-
+            in
+                ( x' * cos phi - y' * sin phi ) *^ e_x
+            ^+^ ( y' * cos phi + x' * sin phi ) *^ e_y
+-- -}
+{-
+                r  = sqrt ( x' ^ 2 + y' ^ 2 )
+                arctgt = atan ( y' / x' )
+                -- a and b are always strictly positive, so we can compute
+                -- the quadrant using only x and y, which are constants.
+                !theta
+                  | x > 0
+                  = arctgt
+                  | x < 0
+                  = if y >= 0 then arctgt + pi else arctgt - pi
+                  | otherwise
+                  = if y >= 0 then 0.5 * pi else -0.5 * pi
+                !phi' = phi + theta
+            in
+                ( r * cos phi' ) *^ e_x
+            ^+^ ( r * sin phi' ) *^ e_y
+-}
+
+    in circleSpline mkPt
+  where
+    e_x, e_y :: D k irec ( I i ( ℝ 2 ) )
+    e_x = pure $ mkI $ ℝ2 1 0
+    e_y = pure $ mkI $ ℝ2 0 1
+
+    scaledBy d x = fmap ( mkI x * ) d
+{-# INLINEABLE ellipseBrush #-}
+
+--------------------------------------------------------------------------------
+-- Tear drop
+
+-- | y-coordinate of the center of mass of the cubic Bézier teardrop
+-- with control points at (0,0), (±0.5,√3/2).
+tearCenter :: Double
+tearCenter = 3 * sqrt 3 / 14
+
+-- | Width of the cubic Bézier teardrop with control points at (0,0), (±0.5,√3/2).
+tearWidth :: Double
+tearWidth = 1 / ( 2 * sqrt 3 )
+
+-- | Height of the cubic Bézier teardrop with control points at (0,0), (±0.5,√3/2).
+tearHeight :: Double
+tearHeight = 3 * sqrt 3 / 8
+
+-- √3/2
+sqrt3_over_2 :: Double
+sqrt3_over_2 = 0.5 * sqrt 3
+
+tearDropBrush :: forall {t} (i :: t) k irec
+             . ( Module
+                  ( D k irec ( I i Double ) )
+                  ( D k irec ( I i ( ℝ 2 ) ) )
+               , Module ( I i Double ) ( T ( I i Double ) )
+               , HasChainRule ( I i Double ) k irec
+               , Representable ( I i Double ) irec
+               , Applicative ( D k irec )
+               , Transcendental ( D k irec ( I i Double ) )
+               )
+              => Proxy# i
+              -> ( forall a. a -> I i a )
+              -> C k irec ( Spline 'Closed () ( I i ( ℝ 2 ) ) )
+tearDropBrush _ mkI =
+  D \ params ->
+    let w, h, phi :: D k irec ( I i Double )
+        w   = runD ( var @_ @k ( Fin 1 ) ) params
+        h   = runD ( var @_ @k ( Fin 2 ) ) params
+        phi = runD ( var @_ @k ( Fin 3 ) ) params
+
+        mkPt :: Double -> Double -> D k irec ( I i ( ℝ 2 ) )
+        mkPt x y
+          -- 1. translate the teardrop so that the center of mass is at the origin
+          -- 2. scale the teardrop so that it has the requested width/height
+          -- 3. rotate
+          = let !x' =   w `scaledBy` (x / tearWidth)
+                !y' = ( h `scaledBy` ( ( y - tearCenter ) / tearHeight) )
+            in
+                ( x' * cos phi - y' * sin phi ) *^ e_x
+            ^+^ ( y' * cos phi + x' * sin phi ) *^ e_y
+
+    in sequenceA $
+         Spline { splineStart  = mkPt 0 0
+                , splineCurves = ClosedCurves Seq.empty $
+                  Bezier3To
+                    ( mkPt  0.5 sqrt3_over_2 )
+                    ( mkPt -0.5 sqrt3_over_2 )
+                    BackToStart () }
+  where
+    e_x, e_y :: D k irec ( I i ( ℝ 2 ) )
+    e_x = pure $ mkI $ ℝ2 1 0
+    e_y = pure $ mkI $ ℝ2 0 1
+
+    scaledBy d x = fmap ( mkI x * ) d
+{-# INLINEABLE tearDropBrush #-}

brush-strokes/src/Math/Algebra/Dual.hs

@@ -1041,12 +1041,20 @@ instance HasChainRule Double 2 ( ℝ 0 ) where
   linearD f v = D0 ( f v )
   value ( D0 v ) = v
   chain _ ( D0 gfx ) = D21 gfx origin origin
+  {-# INLINEABLE konst #-}
+  {-# INLINEABLE linearD #-}
+  {-# INLINEABLE value #-}
+  {-# INLINEABLE chain #-}
 
 instance HasChainRule Double 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
+  {-# INLINEABLE konst #-}
+  {-# INLINEABLE linearD #-}
+  {-# INLINEABLE value #-}
+  {-# INLINEABLE chain #-}
 
 instance HasChainRule Double 2 ( ℝ 1 ) where
 
@@ -1082,6 +1090,10 @@ instance HasChainRule Double 2 ( ℝ 1 ) where
     in $$( chainRule1NQ
          [|| o ||] [|| p ||] [|| s ||]
          [|| df ||] [|| dg ||] )
+  {-# INLINEABLE konst #-}
+  {-# INLINEABLE linearD #-}
+  {-# INLINEABLE value #-}
+  {-# INLINEABLE chain #-}
 
 instance HasChainRule Double 3 ( ℝ 1 ) where
 
@@ -1117,6 +1129,10 @@ instance HasChainRule Double 3 ( ℝ 1 ) where
     in $$( chainRule1NQ
          [|| o ||] [|| p ||] [|| s ||]
          [|| df ||] [|| dg ||] )
+  {-# INLINEABLE konst #-}
+  {-# INLINEABLE linearD #-}
+  {-# INLINEABLE value #-}
+  {-# INLINEABLE chain #-}
 
 instance HasChainRule Double 2 ( ℝ 2 ) where
 
@@ -1152,6 +1168,10 @@ instance HasChainRule Double 2 ( ℝ 2 ) where
     in $$( chainRule1NQ
          [|| o ||] [|| p ||] [|| s ||]
          [|| df ||] [|| dg ||] )
+  {-# INLINEABLE konst #-}
+  {-# INLINEABLE linearD #-}
+  {-# INLINEABLE value #-}
+  {-# INLINEABLE chain #-}
 
 instance HasChainRule Double 3 ( ℝ 2 ) where
 
@@ -1187,6 +1207,10 @@ instance HasChainRule Double 3 ( ℝ 2 ) where
     in $$( chainRule1NQ
          [|| o ||] [|| p ||] [|| s ||]
          [|| df ||] [|| dg ||] )
+  {-# INLINEABLE konst #-}
+  {-# INLINEABLE linearD #-}
+  {-# INLINEABLE value #-}
+  {-# INLINEABLE chain #-}
 
 instance HasChainRule Double 2 ( ℝ 3 ) where
 
@@ -1222,6 +1246,10 @@ instance HasChainRule Double 2 ( ℝ 3 ) where
     in $$( chainRule1NQ
          [|| o ||] [|| p ||] [|| s ||]
          [|| df ||] [|| dg ||] )
+  {-# INLINEABLE konst #-}
+  {-# INLINEABLE linearD #-}
+  {-# INLINEABLE value #-}
+  {-# INLINEABLE chain #-}
 
 instance HasChainRule Double 3 ( ℝ 3 ) where
 
@@ -1257,6 +1285,10 @@ instance HasChainRule Double 3 ( ℝ 3 ) where
     in $$( chainRule1NQ
          [|| o ||] [|| p ||] [|| s ||]
          [|| df ||] [|| dg ||] )
+  {-# INLINEABLE konst #-}
+  {-# INLINEABLE linearD #-}
+  {-# INLINEABLE value #-}
+  {-# INLINEABLE chain #-}
 
 instance HasChainRule Double 2 ( ℝ 4 ) where
 
@@ -1292,6 +1324,10 @@ instance HasChainRule Double 2 ( ℝ 4 ) where
     in $$( chainRule1NQ
          [|| o ||] [|| p ||] [|| s ||]
          [|| df ||] [|| dg ||] )
+  {-# INLINEABLE konst #-}
+  {-# INLINEABLE linearD #-}
+  {-# INLINEABLE value #-}
+  {-# INLINEABLE chain #-}
 
 instance HasChainRule Double 3 ( ℝ 4 ) where
 
@@ -1327,3 +1363,7 @@ instance HasChainRule Double 3 ( ℝ 4 ) where
     in $$( chainRule1NQ
          [|| o ||] [|| p ||] [|| s ||]
          [|| df ||] [|| dg ||] )
+  {-# INLINEABLE konst #-}
+  {-# INLINEABLE linearD #-}
+  {-# INLINEABLE value #-}
+  {-# INLINEABLE chain #-}

brush-strokes/src/Math/Bezier/Stroke.hs

@@ -18,6 +18,14 @@ module Math.Bezier.Stroke
   , line, bezier2, bezier3
   , brushStrokeData, pathAndUsedParams
 
+  -- ** Cusp finding
+  , Preconditioner(..)
+  , IntervalNewtonTree(..), showIntervalNewtonTree
+  , IntervalNewtonStep(..)
+  , IntervalNewtonLeaf(..)
+  , Box
+  , intervalNewtonGS, intervalNewtonGSFrom
+
   )
   where
 
@@ -79,7 +87,7 @@ import Data.Sequence
 import qualified Data.Sequence as Seq
   ( empty, index, length, reverse, singleton, zipWith )
 import Data.Tree
-  ( Tree(..) )
+  ( Tree(..), foldTree )
 
 -- deepseq
 import Control.DeepSeq
@@ -614,11 +622,19 @@ outlineFunction rootAlgo ptParams toBrushParams brushFromParams = \ sp0 crv ->
          ]
         ) ++ "}"
       ]
-    newtTreeLines = map (showTree . showIntervalNewtonTree 1) newtTrees
+    showedTrees = map ( uncurry showIntervalNewtonTree ) newtTrees
+    solLines =
+      [ " #sols: " ++ show (length newtSols)
+      , "#dunno: " ++ show (length newtDunno)
+      , "Tree size: " ++ show (sum @_ @Int $ map (foldTree ( \ _ bs -> 1 + sum bs )) showedTrees)
+      ]
+    treeLines =
+      solLines ++ map showTree showedTrees
 
-    logContents = unlines $ functionDataLines ++ newtTreeLines
+    logContents = unlines $ functionDataLines ++ treeLines
 
-  in logToFile cuspFindingMVar logContents `seq`
+  in trace (unlines solLines) $
+     logToFile cuspFindingMVar logContents `seq`
      OutlineInfo
        { outlineFn = fwdBwd
        , outlineDefiniteCusps  = map ( cuspCoords curves ) newtSols
@@ -1190,7 +1206,7 @@ brushStrokeData co1 co2 path params brush =
         splines :: Seq ( D k ( I i brushParams ) ( 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 = force $ fmap ( uncurryD @k . ( chain @(I i Double) @k dparams_t ) ) splines
+        !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 ) dbrushes_t
@@ -1248,22 +1264,26 @@ gaussSeidel
           !x2  = 𝕀  x2_lo  x2_hi
       in nub $ do
 
-        -- x1' = (b1 - a12 * x2) / a11
+        -- x1' = ( ( b1 - a12 * x2 ) / a11 ) ∩ x1
         x1'_0 <- ( b1 - a12 * x2 ) `extendedDivide` a11
         ( x1'@( 𝕀 x1'_lo x1'_hi ), sub_x1 ) <- x1'_0 `intersect` x1
-        -- x2' = (b2 - a21 * x1') / a22
+        -- x2' = ( ( b2 - a21 * x1' ) / a22 ) ∩ x2
         x2'_0 <- ( b2 - a21 * x1' ) `extendedDivide` a22
         ( 𝕀 x2'_lo x2'_hi, sub_x2 ) <- x2'_0 `intersect` x2
 
         return ( ( 𝕀 ( ℝ1 x1'_lo ) ( ℝ1 x1'_hi ), 𝕀 ( ℝ1 x2'_lo ) ( ℝ1 x2'_hi ) )
                , sub_x1 && sub_x2 )
 
+-- | 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.
 intersect :: 𝕀 Double -> 𝕀 Double -> [ ( 𝕀 Double, Bool ) ]
 intersect ( 𝕀 lo1 hi1 ) ( 𝕀 lo2 hi2 )
   | lo > hi
   = [ ]
   | otherwise
-  = [ ( 𝕀 lo hi, lo == lo1 && hi == hi1 ) ]
+  = [ ( 𝕀 lo hi, ( lo1 > lo2 && hi1 < hi2 ) || ( lo == hi ) ) ]
   where
     lo = max lo1 lo2
     hi = min hi1 hi2
@@ -1283,7 +1303,7 @@ extendedRecip x@( 𝕀 lo hi )
 -- | Computes the brush stroke coordinates of a cusp from
 -- the @(t,s)@ parameter values.
 cuspCoords :: ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) )
-           -> ( 𝕀ℝ 1, Int, 𝕀ℝ 1 )
+           -> Box
            -> Cusp
 cuspCoords eqs ( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ), i, 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
   | StrokeDatum
@@ -1308,18 +1328,16 @@ data Preconditioner
 -- | A tree recording the steps taken with the interval Newton method.
 data IntervalNewtonTree d
   = IntervalNewtonLeaf (IntervalNewtonLeaf d)
-  | IntervalNewtonStep (IntervalNewtonStep d) [(Double, IntervalNewtonTree d)]
+  | IntervalNewtonStep (IntervalNewtonStep d) [(d, IntervalNewtonTree d)]
 
 data IntervalNewtonStep d
-  = IntervalNewtonContraction [d] [d]
+  = IntervalNewtonContraction [d]
   | IntervalNewtonBisection (String, Double)
 
 instance Show d => Show (IntervalNewtonStep d) where
-  showsPrec _ ( IntervalNewtonContraction v w )
+  showsPrec _ ( IntervalNewtonContraction v )
     = showString "N "
     . showsPrec 0 v
-    . showString " "
-    . showsPrec 0 w
   showsPrec _ ( IntervalNewtonBisection (coord, w) )
     = showString "B "
     . showParen True
@@ -1330,17 +1348,38 @@ instance Show d => Show (IntervalNewtonStep d) where
 
 data IntervalNewtonLeaf d
   = TooSmall d
-  | NoSolution d
+  | NoSolution String d
   deriving stock Show
 
-showIntervalNewtonTree :: Show d => Double -> IntervalNewtonTree d -> Tree String
+showIntervalNewtonTree :: Box -> IntervalNewtonTree Box -> Tree String
-showIntervalNewtonTree area (IntervalNewtonLeaf l) = Node (showArea area ++ " " ++ show l) []
+showIntervalNewtonTree cand (IntervalNewtonLeaf l) = Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ show l) []
-showIntervalNewtonTree area (IntervalNewtonStep s ts)
+showIntervalNewtonTree cand (IntervalNewtonStep s ts)
-  = Node (showArea area ++ " " ++ show s) $ map (\ (d,t) -> showIntervalNewtonTree d t) ts
+  = Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ show s) $ map (\ (c,t) -> showIntervalNewtonTree c t) ts
+
+boxArea :: Box -> Double
+boxArea ( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ), _, 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
+  = abs ( t_hi - t_lo ) * abs ( s_hi - s_lo )
 
 showArea :: Double -> [Char]
 showArea area = "(area " ++ showFFloat (Just 6) area "" ++ ")"
 
+type Box = ( 𝕀ℝ 1, Int, 𝕀ℝ 1 )
+
+intervalNewtonGS
+  :: Preconditioner
+  -> Double
+  -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
+  -> ( [ ( Box, IntervalNewtonTree Box ) ], ( [ Box ], [ Box ] ) )
+intervalNewtonGS precondMethod minWidth eqs =
+  foldMap
+    ( intervalNewtonGSFrom precondMethod minWidth eqs )
+    initBoxes
+  where
+    initBoxes =
+      [ ( 𝕀 ( ℝ1 zero ) ( ℝ1 one ), i, 𝕀 ( ℝ1 zero ) ( ℝ1 one ) )
+      | i <- [ 0 .. length ( eqs ( 𝕀 ( ℝ1 zero ) ( ℝ1 one ) ) ) - 1 ]
+      ]
+
 -- | Interval Newton method with Gauss–Seidel step for inversion
 -- of the interval Jacobian.
 --
@@ -1348,35 +1387,25 @@ showArea area = "(area " ++ showFFloat (Just 6) area "" ++ ")"
 -- @sols@ are boxes that contain a unique solution (and to which Newton's method
 -- will converge starting from anywhere inside the box), and @dunno@ are small
 -- boxes which might or might not contain solutions.
-intervalNewtonGS :: Preconditioner
+intervalNewtonGSFrom
-                 -> Double
+  :: Preconditioner
-                 -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
+  -> Double
-                 -> ( [ IntervalNewtonTree ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ], ( [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ], [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] ) )
+  -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-intervalNewtonGS precondMethod minWidth eqs =
+  -> Box
-  first concat $ runWriter $
+  -> ( [ ( Box, IntervalNewtonTree Box ) ], ( [ Box ], [ Box ] ) )
-    traverse go
+intervalNewtonGSFrom precondMethod minWidth eqs initBox =
-      [ ( 𝕀 ( ℝ1 zero ) ( ℝ1 one ), i, 𝕀 ( ℝ1 zero ) ( ℝ1 one ) )
+  runWriter $
-      | i <- [ 0 .. length ( eqs ( 𝕀 ( ℝ1 zero ) ( ℝ1 one ) ) ) - 1 ]
+    map ( initBox , ) <$> go initBox
-      ]
 
   where
-    zero, one :: Double
-    zero = 1e-6
-    one = 1 - zero
-    {-# INLINE zero #-}
-    {-# INLINE one #-}
 
     recur f cands = do
-      rest <- traverse ( \ c -> do { trees <- go c; return [ (boxArea c, tree) | tree <- trees ] } ) cands
+      rest <- traverse ( \ c -> do { trees <- go c; return [ (c, tree) | tree <- trees ] } ) cands
       return [ f $ concat rest ]
 
-    boxArea :: ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) -> Double
+    go :: Box -- box to work on
-    boxArea ( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ), _, 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
+       -> Writer ( [ Box ], [ Box ] )
-      = abs ( t_hi - t_lo ) * abs ( s_hi - s_lo )
+                 [ IntervalNewtonTree Box ]
-
-    go :: ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) -- box to work on
-       -> Writer ( [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ], [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] )
-                 [ IntervalNewtonTree ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ]
     go cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) )
             , i
             , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
@@ -1387,23 +1416,37 @@ intervalNewtonGS precondMethod minWidth eqs =
            tell ( [ cand ], [] )
            return [ IntervalNewtonLeaf res ]
 
-      | StrokeDatum { ee = D22 ee _ _ _ _ _
+      | StrokeDatum { ee = D22 ee ( T _ee_t ) ( T _ee_s ) _ _ _
                     , 𝛿E𝛿sdcdt = D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) }
         <- ( eqs t `Seq.index` i ) s
 
-      , StrokeDatum { 𝛿E𝛿sdcdt = D12 ( T f_mid ) ( T ( T _f_t_mid ) ) ( T ( T _f_s_mid ) ) }
+      , StrokeDatum { ee = D22 _ee_mid _ _ _ _ _
+                    , 𝛿E𝛿sdcdt = D12 ( T f_mid ) ( T ( T f_t_mid ) ) ( T ( T f_s_mid ) ) }
         <- ( eqs i_t_mid `Seq.index` i ) i_s_mid
+      , let ee_potential_zero = inf ee <= ℝ1 0 && sup ee >= ℝ1 0
+            𝛿E𝛿sdcdt_potential_zero = cmpℝ2 (<=) ( inf f ) ( ℝ2 0 0 ) && cmpℝ2 (>=) ( sup f ) ( ℝ2 0 0 )
       = if | -- Check the envelope equation interval contains zero.
-             inf ee <= ℝ1 0
+             ee_potential_zero
-           , sup ee >= ℝ1 0
              -- Check the 𝛿E𝛿sdcdt interval contains zero.
-           , cmpℝ2 (<=) ( inf f ) ( ℝ2 0 0 )
+           , 𝛿E𝛿sdcdt_potential_zero
-           , cmpℝ2 (>=) ( sup f ) ( ℝ2 0 0 )
            -> let -- Interval Newton method: take one Gauss–Seidel step
-                  -- for the equation f'(X) v = - f(x_mid).
+                  -- for the equation f'(x) v = - f(x_mid),
+                  -- where f = 𝛿E/𝛿s * dc/dt
                   !( a, b ) = precondition precondMethod
-                                ( f_t, f_s )
+                                ( f_t_mid, f_s_mid )
                                 ( f_t, f_s ) ( neg f_mid )
+                  --(a, b)
+                  --  | 𝕀 (ℝ1 ee_lo) (ℝ1 ee_hi) <- ee_mid
+                  --  , 𝕀 (ℝ1 ee_t_lo) (ℝ1 ee_t_hi) <- ee_t
+                  --  , 𝕀 (ℝ1 ee_s_lo) (ℝ1 ee_s_hi) <- ee_s
+                  --  , 𝕀 (ℝ2 fx_lo fy_lo) (ℝ2 fx_hi fy_hi) <- f_mid
+                  --  , 𝕀 (ℝ2 f_tx_lo f_ty_lo) (ℝ2 f_tx_hi f_ty_hi) <- f_t
+                  --  , 𝕀 (ℝ2 f_sx_lo f_sy_lo) (ℝ2 f_sx_hi f_sy_hi) <- f_s
+                  --  = ( ( 𝕀 (ℝ2 f_tx_lo ee_t_lo) (ℝ2 f_tx_hi ee_t_hi)
+                  --      , 𝕀 (ℝ2 f_sx_lo ee_s_lo) (ℝ2 f_sx_hi ee_s_hi)
+                  --      )
+                  --    , neg $ 𝕀 (ℝ2 fx_lo ee_lo) (ℝ2 fx_hi ee_hi)
+                  --    )
 
                   !gsGuesses = gaussSeidel a b
                                  ( coerce ( (-) @( 𝕀 Double ) ) t i_t_mid
@@ -1418,7 +1461,9 @@ intervalNewtonGS precondMethod minWidth eqs =
                    let !(done, todo) = bimap ( map ( mkGuess . fst ) ) ( map ( mkGuess . fst ) )
                                      $ partition snd gsGuesses
                    in do tell ([], done)
-                         recur (IntervalNewtonStep (IntervalNewtonContraction done todo)) todo
+                         case todo of
+                           [] -> return [ IntervalNewtonLeaf $ NoSolution "GaussSeidel" cand ]
+                           _ -> recur (IntervalNewtonStep (IntervalNewtonContraction done)) todo
                  else
                     -- Gauss–Seidel failed to shrink the boxes.
                     -- Bisect along the widest dimension instead.
@@ -1435,12 +1480,12 @@ intervalNewtonGS precondMethod minWidth eqs =
 
            -- Box doesn't contain a solution: discard it.
            | otherwise
-           -> return [ IntervalNewtonLeaf $ NoSolution cand ]
+           -> return [ IntervalNewtonLeaf $ NoSolution (if ee_potential_zero then "dc/dt" else "ee") cand ]
       where
         t_mid = 0.5 * ( t_lo + t_hi )
         s_mid = 0.5 * ( s_lo + s_hi )
-        i_t_mid = 𝕀 ( ℝ1 t_mid ) ( ℝ1 t_mid )
+        i_t_mid = singleton ( ℝ1 t_mid )
-        i_s_mid = 𝕀 ( ℝ1 s_mid ) ( ℝ1 s_mid )
+        i_s_mid = singleton ( ℝ1 s_mid )
         mkGuess ( t0, s0 ) = ( coerce ( (+) @( 𝕀 Double ) ) t0 i_t_mid
                              , i
                              , coerce ( (+) @( 𝕀 Double ) ) s0 i_s_mid )
@@ -1454,6 +1499,12 @@ intervalNewtonGS precondMethod minWidth eqs =
                 !( 𝕀 y'_lo y'_hi ) = negate $ 𝕀 y_lo y_hi
             in 𝕀 ( ℝ2 x'_lo y'_lo ) ( ℝ2 x'_hi y'_hi )
 
+zero, one :: Double
+zero = 1e-6
+one = 1 - zero
+{-# INLINE zero #-}
+{-# INLINE one #-}
+
 precondition :: Preconditioner
              -> ( 𝕀ℝ 2, 𝕀ℝ 2 )
              -> ( 𝕀ℝ 2, 𝕀ℝ 2 )
@@ -1481,20 +1532,19 @@ precondition meth jac_mid a@( a1, a2 ) b =
 
 scale :: Double -> 𝕀ℝ 2 -> 𝕀ℝ 2
 scale s ( 𝕀 ( ℝ2 a1_lo a2_lo ) ( ℝ2 a1_hi a2_hi ) )
-  | 𝕀 b1_lo b1_hi <- scaleInterval s ( 𝕀 a1_lo a1_hi )
-  , 𝕀 b2_lo b2_hi <- scaleInterval s ( 𝕀 a2_lo a2_hi )
   = 𝕀 ( ℝ2 b1_lo b2_lo ) ( ℝ2 b1_hi b2_hi )
+  where
+    𝕀 b1_lo b1_hi = scaleInterval s ( 𝕀 a1_lo a1_hi )
+    𝕀 b2_lo b2_hi = scaleInterval s ( 𝕀 a2_lo a2_hi )
 
 matMulVec :: ( ℝ 2, ℝ 2 ) -> 𝕀ℝ 2 -> 𝕀ℝ 2
 matMulVec ( ℝ2 a11 a21, ℝ2 a12 a22 ) ( 𝕀 ( ℝ2 u_lo v_lo ) ( ℝ2 u_hi v_hi ) )
-  | 𝕀 u'_lo u'_hi <-
-         𝕀 a11 a11 * 𝕀 u_lo u_hi
-       + 𝕀 a12 a12 * 𝕀 v_lo v_hi
-  , 𝕀 v'_lo v'_hi <-
-         𝕀 a21 a21 * 𝕀 u_lo u_hi
-       + 𝕀 a22 a22 * 𝕀 v_lo v_hi
   = 𝕀 ( ℝ2 u'_lo v'_lo ) ( ℝ2 u'_hi v'_hi )
-
+  where
+    u = 𝕀 u_lo u_hi
+    v = 𝕀 v_lo v_hi
+    𝕀 u'_lo u'_hi = scaleInterval a11 u + scaleInterval a12 v
+    𝕀 v'_lo v'_hi = scaleInterval a21 u + scaleInterval a22 v
 
 cmpℝ2 :: ( Double -> Double -> Bool ) -> ℝ 2 -> ℝ 2 -> Bool
 cmpℝ2 cmp ( ℝ2 x1 y1 ) ( ℝ2 x2 y2 )

brush-strokes/src/Math/Differentiable.hs

@@ -53,6 +53,7 @@ class
   , Transcendental ( D k ( I i u ) ( I i Double ) )
   , Applicative ( D k ( I i u ) )
   , Representable ( I i Double ) ( I i u )
+  , RepDim ( I i u ) ~ RepDim u
   ) => DiffInterp k i u
 instance
   ( Differentiable k i u
@@ -63,4 +64,5 @@ instance
   , Transcendental ( D k ( I i u ) ( I i Double ) )
   , Applicative ( D k ( I i u ) )
   , Representable ( I i Double ) ( I i u )
+  , RepDim ( I i u ) ~ RepDim u
   ) => DiffInterp k i u

brush-strokes/src/Math/Interval.hs

@@ -75,6 +75,7 @@ instance Torsor ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where
 aabb :: ( Representable r v, Ord r, Functor f )
      => f ( 𝕀 v ) -> ( f ( 𝕀 r ) -> 𝕀 r ) -> 𝕀 v
 aabb fv extrema = tabulate \ i -> extrema ( fmap ( `index` i ) fv )
+{-# INLINEABLE aabb #-}
 
 --------------------------------------------------------------------------------
 
@@ -109,7 +110,7 @@ instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ 4 ) ) where
   k *^ T a    = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] )
 
 instance Inner  ( 𝕀 Double ) ( T ( 𝕀ℝ 2 ) ) where
-  T ( 𝕀 ( ℝ2 x1_lo y1_lo ) ( ℝ2 x1_hi y1_hi ) )  ^.^
+  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

brush-strokes/src/Math/Interval/Internal.hs

@@ -216,6 +216,7 @@ instance Representable r u => Representable ( 𝕀 r ) ( 𝕀 u ) where
     let !lo = tabulate @r @u ( \ i -> inf $ f i )
         !hi = tabulate @r @u ( \ i -> sup $ f i )
     in 𝕀 lo hi
+  {-# INLINE tabulate #-}
 
   index d i =
     case d of
@@ -223,5 +224,6 @@ instance Representable r u => Representable ( 𝕀 r ) ( 𝕀 u ) where
         let !lo_i = index @r @u lo i
             !hi_i = index @r @u hi i
         in 𝕀 lo_i hi_i
+  {-# INLINE index #-}
 
 deriving via Sum ( 𝕀 Double ) instance Module ( 𝕀 Double ) ( T ( 𝕀 Double ) )

brush-strokes/src/Math/Linear/Internal.hs

@@ -117,6 +117,7 @@ projection :: ( Representable r u, Representable r v )
            -> u -> v
 projection f = \ u ->
   tabulate \ i -> index u ( f i )
+{-# INLINEABLE projection #-}
 
 injection :: ( Representable r u, Representable r v )
           => ( Fin ( RepDim v ) -> MFin ( RepDim u ) )
@@ -125,6 +126,7 @@ injection f = \ u v ->
   tabulate \ i -> case f i of
     MFin 0 -> index v i
     MFin j -> index u ( Fin j )
+{-# INLINEABLE injection #-}
 
 --------------------------------------------------------------------------------
 
@@ -161,29 +163,39 @@ instance RepresentableQ Double ( ℝ 4 ) where
 
 instance Representable Double ( ℝ 0 ) where
   tabulate _ = ℝ0
+  {-# INLINE tabulate #-}
   index _ _ = 0
+  {-# INLINE index #-}
 
 instance Representable Double ( ℝ 1 ) where
   tabulate f = ℝ1 ( f ( Fin 1 ) )
+  {-# INLINE tabulate #-}
   index p _ = unℝ1 p
+  {-# INLINE index #-}
 
 instance Representable Double ( ℝ 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 Double ( ℝ 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 Double ( ℝ 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 #-}

src/metabrushes/MetaBrush/Asset/Brushes.hs

@@ -16,11 +16,7 @@ import Prelude
   hiding
     ( Num(..), Floating(..), (^), (/), fromInteger, fromRational )
 import GHC.Exts
-  ( Proxy#, fromString )
+  ( fromString )
-
--- containers
-import qualified Data.Sequence as Seq
-  ( fromList, empty )
 
 -- text
 import Data.Text
@@ -32,15 +28,12 @@ import Data.HashMap.Strict
 import qualified Data.HashMap.Strict as HashMap
   ( fromList, lookup )
 
--- MetaBrush
+-- brush-strokes
-import Math.Algebra.Dual
+import Calligraphy.Brushes
-import Math.Bezier.Spline
-import Math.Differentiable
-  ( I )
 import Math.Linear
-import Math.Module
-  ( Module((^+^), (*^)) )
 import Math.Ring
+
+-- MetaBrush
 import MetaBrush.Brush
   ( Brush(..), SomeBrush(..), WithParams(..) )
 import MetaBrush.Records
@@ -61,12 +54,6 @@ brushes = HashMap.fromList
 
 --------------------------------------------------------------------------------
 
--- | Root of @(Sqrt[2] (4 + 3 κ) - 16) (2 - 3 κ)^2 - 8 (1 - 3 κ) Sqrt[8 - 24 κ + 12 κ^2 + 8 κ^3 + 3 κ^4]@.
---
--- Used to approximate circles and ellipses with Bézier curves.
-κ :: Double
-κ = 0.5519150244935105707435627227925
-
 type CircleBrushFields = '[ "r" ]
 -- | A circular brush with the given radius.
 circle :: Brush CircleBrushFields
@@ -86,22 +73,6 @@ ellipse = BrushData "ellipse" ( WithParams deflts ellipseBrush )
     deflts = MkR ( ℝ3 1 1 0 )
 {-# INLINE ellipse #-}
 
--- | y-coordinate of the center of mass of the cubic Bézier teardrop
--- with control points at (0,0), (±0.5,√3/2).
-tearCenter :: Double
-tearCenter = 3 * sqrt 3 / 14
-
--- | Width of the cubic Bézier teardrop with control points at (0,0), (±0.5,√3/2).
-tearWidth :: Double
-tearWidth = 1 / ( 2 * sqrt 3 )
-
--- | Height of the cubic Bézier teardrop with control points at (0,0), (±0.5,√3/2).
-tearHeight :: Double
-tearHeight = 3 * sqrt 3 / 8
-
-sqrt3_over_2 :: Double
-sqrt3_over_2 = 0.5 * sqrt 3
-
 type TearDropBrushFields = '[ "w", "h", "phi" ]
 -- | A tear-drop shape with the given width, height and angle of rotation.
 tearDrop :: Brush TearDropBrushFields
@@ -110,158 +81,3 @@ tearDrop = BrushData "tear-drop" ( WithParams deflts tearDropBrush )
     deflts :: Record TearDropBrushFields
     deflts = MkR ( ℝ3 1 2.25 0 )
 {-# INLINE tearDrop #-}
-
---------------------------------------------------------------------------------
--- Differentiable brushes.
-
-circleSpline :: forall {t} (i :: t) k u v
-             .  Applicative ( D k ( I i u ) )
-             => ( Double -> Double -> D k ( I i u ) ( I i v ) )
-             -> D k ( I i u ) ( Spline 'Closed () ( I i v ) )
-circleSpline p = sequenceA $
-  Spline { splineStart  = p 1 0
-         , splineCurves = ClosedCurves crvs lastCrv }
-  where
-    crvs = Seq.fromList
-      [ Bezier3To ( p  1  κ ) ( p  κ  1 ) ( NextPoint ( p  0  1 ) ) ()
-      , Bezier3To ( p -κ  1 ) ( p -1  κ ) ( NextPoint ( p -1  0 ) ) ()
-      , Bezier3To ( p -1 -κ ) ( p -κ -1 ) ( NextPoint ( p  0 -1 ) ) ()
-      ]
-    lastCrv =
-        Bezier3To ( p  κ -1 ) ( p  1 -κ ) BackToStart ()
-{-# INLINE circleSpline #-}
-
-circleBrush :: forall {t} (i :: t) k irec
-            . ( irec ~ I i ( Record CircleBrushFields )
-              , Module
-                 ( D k irec ( I i Double ) )
-                 ( D k irec ( I i ( ℝ 2 ) ) )
-              , Module ( I i Double ) ( T ( I i Double ) )
-              , HasChainRule ( I i Double ) k irec
-              , Representable ( I i Double ) irec
-              , Applicative ( D k irec )
-              )
-            => Proxy# i
-            -> ( forall a. a -> I i a )
-            -> C k irec ( Spline 'Closed () ( I i ( ℝ 2 ) ) )
-circleBrush _ mkI =
-  D \ params ->
-    let r :: D k irec ( I i Double )
-        r = runD ( var @_ @k ( Fin 1 ) ) params
-        mkPt :: Double -> Double -> D k irec ( I i ( ℝ 2 ) )
-        mkPt x y
-          =   ( r `scaledBy` x ) *^ e_x
-          ^+^ ( r `scaledBy` y ) *^ e_y
-    in circleSpline @i @k @( Record CircleBrushFields ) @( ℝ 2 ) mkPt
-  where
-    e_x, e_y :: D k irec ( I i ( ℝ 2 ) )
-    e_x = pure $ mkI $ ℝ2 1 0
-    e_y = pure $ mkI $ ℝ2 0 1
-
-    scaledBy d x = fmap ( mkI x * ) d
-{-# INLINEABLE circleBrush #-}
-
-ellipseBrush :: forall {t} (i :: t) k irec
-             . ( irec ~ I i ( Record EllipseBrushFields )
-               , Module
-                  ( D k irec ( I i Double ) )
-                  ( D k irec ( I i ( ℝ 2 ) ) )
-               , Module ( I i Double ) ( T ( I i Double ) )
-               , HasChainRule ( I i Double ) k irec
-               , Representable ( I i Double ) irec
-               , Applicative ( D k irec )
-               , Transcendental ( D k irec ( I i Double ) )
-                 -- TODO: make a synonym for the above...
-                 -- it seems DiffInterp isn't exactly right
-               )
-              => Proxy# i
-              -> ( forall a. a -> I i a )
-              -> C k irec ( Spline 'Closed () ( I i ( ℝ 2 ) ) )
-ellipseBrush _ mkI =
-  D \ params ->
-    let a, b, phi :: D k irec ( I i Double )
-        a   = runD ( var @_ @k ( Fin 1 ) ) params
-        b   = runD ( var @_ @k ( Fin 2 ) ) params
-        phi = runD ( var @_ @k ( Fin 3 ) ) params
-        mkPt :: Double -> Double -> D k irec ( I i ( ℝ 2 ) )
-        mkPt x y
-          = let !x' = a `scaledBy` x
-                !y' = b `scaledBy` y
--- {-
-            in
-                ( x' * cos phi - y' * sin phi ) *^ e_x
-            ^+^ ( y' * cos phi + x' * sin phi ) *^ e_y
--- -}
-{-
-                r  = sqrt ( x' ^ 2 + y' ^ 2 )
-                arctgt = atan ( y' / x' )
-                -- a and b are always strictly positive, so we can compute
-                -- the quadrant using only x and y, which are constants.
-                !theta
-                  | x > 0
-                  = arctgt
-                  | x < 0
-                  = if y >= 0 then arctgt + pi else arctgt - pi
-                  | otherwise
-                  = if y >= 0 then 0.5 * pi else -0.5 * pi
-                !phi' = phi + theta
-            in
-                ( r * cos phi' ) *^ e_x
-            ^+^ ( r * sin phi' ) *^ e_y
--}
-
-    in circleSpline @i @k @( Record EllipseBrushFields ) @( ℝ 2 ) mkPt
-  where
-    e_x, e_y :: D k irec ( I i ( ℝ 2 ) )
-    e_x = pure $ mkI $ ℝ2 1 0
-    e_y = pure $ mkI $ ℝ2 0 1
-
-    scaledBy d x = fmap ( mkI x * ) d
-{-# INLINEABLE ellipseBrush #-}
-
-tearDropBrush :: forall {t} (i :: t) k irec
-             . ( irec ~ I i ( Record TearDropBrushFields )
-               , Module
-                  ( D k irec ( I i Double ) )
-                  ( D k irec ( I i ( ℝ 2 ) ) )
-               , Module ( I i Double ) ( T ( I i Double ) )
-               , HasChainRule ( I i Double ) k irec
-               , Representable ( I i Double ) irec
-               , Applicative ( D k irec )
-               , Transcendental ( D k irec ( I i Double ) )
-               )
-              => Proxy# i
-              -> ( forall a. a -> I i a )
-              -> C k irec ( Spline 'Closed () ( I i ( ℝ 2 ) ) )
-tearDropBrush _ mkI =
-  D \ params ->
-    let w, h, phi :: D k irec ( I i Double )
-        w   = runD ( var @_ @k ( Fin 1 ) ) params
-        h   = runD ( var @_ @k ( Fin 2 ) ) params
-        phi = runD ( var @_ @k ( Fin 3 ) ) params
-
-        mkPt :: Double -> Double -> D k irec ( I i ( ℝ 2 ) )
-        mkPt x y
-          -- 1. translate the teardrop so that the center of mass is at the origin
-          -- 2. scale the teardrop so that it has the requested width/height
-          -- 3. rotate
-          = let !x' =   w `scaledBy` (x / tearWidth)
-                !y' = ( h `scaledBy` ( ( y - tearCenter ) / tearHeight) )
-            in
-                ( x' * cos phi - y' * sin phi ) *^ e_x
-            ^+^ ( y' * cos phi + x' * sin phi ) *^ e_y
-
-    in sequenceA $
-         Spline { splineStart  = mkPt 0 0
-                , splineCurves = ClosedCurves Seq.empty $
-                  Bezier3To
-                    ( mkPt  0.5 sqrt3_over_2 )
-                    ( mkPt -0.5 sqrt3_over_2 )
-                    BackToStart () }
-  where
-    e_x, e_y :: D k irec ( I i ( ℝ 2 ) )
-    e_x = pure $ mkI $ ℝ2 1 0
-    e_y = pure $ mkI $ ℝ2 0 1
-
-    scaledBy d x = fmap ( mkI x * ) d
-{-# INLINEABLE tearDropBrush #-}

src/metabrushes/MetaBrush/Records.hs

@@ -139,9 +139,12 @@ instance ( Torsor ( T ( 𝕀ℝ ( Length ks ) ) ) ( 𝕀ℝ ( Length ks ) )
             T ( 𝕀 ( MkR c_lo ) ( MkR c_hi ) )
 
 type instance RepDim ( Record ks ) = Length ks
-deriving newtype
+instance Representable r ( ℝ ( Length ks ) )
-  instance Representable r ( ℝ ( Length ks ) )
+      => Representable r ( Record ks ) where
-        => Representable r ( Record ks )
+  tabulate f = Record $ tabulate f
+  {-# INLINE tabulate #-}
+  index f (Record r) = index f r
+  {-# INLINE index #-}
 
 type instance D k ( Record ks ) = D k ( ℝ ( Length ks ) )
 deriving newtype instance HasChainRule Double 2 ( ℝ ( Length ks ) )