From d1b3765335ced46ffd4bdd45534dedad5a94526f Mon Sep 17 00:00:00 2001 From: sheaf Date: Fri, 23 Feb 2024 17:03:28 +0100 Subject: [PATCH] Split out benchmark for cusp finding --- brush-strokes/bench/Main.hs | 272 ++++++++++++++++++++ brush-strokes/brush-strokes.cabal | 78 ++++-- brush-strokes/cabal.project | 4 + brush-strokes/src/Calligraphy/Brushes.hs | 215 ++++++++++++++++ brush-strokes/src/Math/Algebra/Dual.hs | 40 +++ brush-strokes/src/Math/Bezier/Stroke.hs | 176 ++++++++----- brush-strokes/src/Math/Differentiable.hs | 2 + brush-strokes/src/Math/Interval.hs | 3 +- brush-strokes/src/Math/Interval/Internal.hs | 2 + brush-strokes/src/Math/Linear/Internal.hs | 12 + src/metabrushes/MetaBrush/Asset/Brushes.hs | 194 +------------- src/metabrushes/MetaBrush/Records.hs | 9 +- 12 files changed, 723 insertions(+), 284 deletions(-) create mode 100644 brush-strokes/bench/Main.hs create mode 100644 brush-strokes/src/Calligraphy/Brushes.hs diff --git a/brush-strokes/bench/Main.hs b/brush-strokes/bench/Main.hs new file mode 100644 index 0000000..3df3cb3 --- /dev/null +++ b/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 ) diff --git a/brush-strokes/brush-strokes.cabal b/brush-strokes/brush-strokes.cabal index d8e06aa..e825794 100644 --- a/brush-strokes/brush-strokes.cabal +++ b/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 - >= 2.18 && < 2.22 + , primitive + ^>= 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 diff --git a/brush-strokes/cabal.project b/brush-strokes/cabal.project index 8f7e859..51f427c 100644 --- a/brush-strokes/cabal.project +++ b/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 -- ------------- diff --git a/brush-strokes/src/Calligraphy/Brushes.hs b/brush-strokes/src/Calligraphy/Brushes.hs new file mode 100644 index 0000000..774c2b1 --- /dev/null +++ b/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 #-} diff --git a/brush-strokes/src/Math/Algebra/Dual.hs b/brush-strokes/src/Math/Algebra/Dual.hs index 9296c1e..4f73d74 100644 --- a/brush-strokes/src/Math/Algebra/Dual.hs +++ b/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 #-} diff --git a/brush-strokes/src/Math/Bezier/Stroke.hs b/brush-strokes/src/Math/Bezier/Stroke.hs index c4a343a..0c138cc 100644 --- a/brush-strokes/src/Math/Bezier/Stroke.hs +++ b/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 area (IntervalNewtonLeaf l) = Node (showArea area ++ " " ++ show l) [] -showIntervalNewtonTree area (IntervalNewtonStep s ts) - = Node (showArea area ++ " " ++ show s) $ map (\ (d,t) -> showIntervalNewtonTree d t) ts +showIntervalNewtonTree :: Box -> IntervalNewtonTree Box -> Tree String +showIntervalNewtonTree cand (IntervalNewtonLeaf l) = Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ show l) [] +showIntervalNewtonTree cand (IntervalNewtonStep s 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 - -> Double - -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) - -> ( [ IntervalNewtonTree ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ], ( [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ], [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] ) ) -intervalNewtonGS precondMethod minWidth eqs = - first concat $ runWriter $ - traverse go - [ ( 𝕀 ( ℝ1 zero ) ( ℝ1 one ), i, 𝕀 ( ℝ1 zero ) ( ℝ1 one ) ) - | i <- [ 0 .. length ( eqs ( 𝕀 ( ℝ1 zero ) ( ℝ1 one ) ) ) - 1 ] - ] +intervalNewtonGSFrom + :: Preconditioner + -> Double + -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) + -> Box + -> ( [ ( Box, IntervalNewtonTree Box ) ], ( [ Box ], [ Box ] ) ) +intervalNewtonGSFrom precondMethod minWidth eqs initBox = + runWriter $ + 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 - boxArea ( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ), _, 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) ) - = abs ( t_hi - t_lo ) * abs ( s_hi - s_lo ) - - go :: ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) -- box to work on - -> Writer ( [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ], [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] ) - [ IntervalNewtonTree ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] + go :: Box -- box to work on + -> Writer ( [ Box ], [ Box ] ) + [ IntervalNewtonTree Box ] 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 - , sup ee >= ℝ1 0 + ee_potential_zero -- Check the 𝛿E𝛿sdcdt interval contains zero. - , cmpℝ2 (<=) ( inf f ) ( ℝ2 0 0 ) - , cmpℝ2 (>=) ( sup f ) ( ℝ2 0 0 ) + , 𝛿E𝛿sdcdt_potential_zero -> 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_s_mid = 𝕀 ( ℝ1 s_mid ) ( ℝ1 s_mid ) + i_t_mid = singleton ( ℝ1 t_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 ) diff --git a/brush-strokes/src/Math/Differentiable.hs b/brush-strokes/src/Math/Differentiable.hs index 8ecbe97..521bd9a 100644 --- a/brush-strokes/src/Math/Differentiable.hs +++ b/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 diff --git a/brush-strokes/src/Math/Interval.hs b/brush-strokes/src/Math/Interval.hs index 1a656d7..f3f2d4a 100644 --- a/brush-strokes/src/Math/Interval.hs +++ b/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 diff --git a/brush-strokes/src/Math/Interval/Internal.hs b/brush-strokes/src/Math/Interval/Internal.hs index 294c11c..aa649fe 100644 --- a/brush-strokes/src/Math/Interval/Internal.hs +++ b/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 ) ) diff --git a/brush-strokes/src/Math/Linear/Internal.hs b/brush-strokes/src/Math/Linear/Internal.hs index 9d360fc..8830247 100644 --- a/brush-strokes/src/Math/Linear/Internal.hs +++ b/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 #-} diff --git a/src/metabrushes/MetaBrush/Asset/Brushes.hs b/src/metabrushes/MetaBrush/Asset/Brushes.hs index 52cbeb9..7f6b160 100644 --- a/src/metabrushes/MetaBrush/Asset/Brushes.hs +++ b/src/metabrushes/MetaBrush/Asset/Brushes.hs @@ -16,11 +16,7 @@ import Prelude hiding ( Num(..), Floating(..), (^), (/), fromInteger, fromRational ) import GHC.Exts - ( Proxy#, fromString ) - --- containers -import qualified Data.Sequence as Seq - ( fromList, empty ) + ( fromString ) -- text import Data.Text @@ -32,15 +28,12 @@ import Data.HashMap.Strict import qualified Data.HashMap.Strict as HashMap ( fromList, lookup ) --- MetaBrush -import Math.Algebra.Dual -import Math.Bezier.Spline -import Math.Differentiable - ( I ) +-- brush-strokes +import Calligraphy.Brushes 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 #-} diff --git a/src/metabrushes/MetaBrush/Records.hs b/src/metabrushes/MetaBrush/Records.hs index 6017317..0d86207 100644 --- a/src/metabrushes/MetaBrush/Records.hs +++ b/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 ) ) - => Representable r ( Record ks ) +instance Representable r ( ℝ ( Length ks ) ) + => Representable r ( Record ks ) where + 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 ) )