From 2289468a8464308ac3dab8d70d904aa069634c11 Mon Sep 17 00:00:00 2001 From: sheaf Date: Wed, 28 Feb 2024 17:20:34 +0100 Subject: [PATCH] Fixes and restructuring --- MetaBrush.cabal | 43 -- brush-strokes/bench/Main.hs | 370 ----------- brush-strokes/brush-strokes.cabal | 68 +- brush-strokes/src/cusps/bench/Main.hs | 616 ++++++++++++++++++ .../src/cusps/inspect}/Main.hs | 3 +- .../cusps/inspect}/Math/Interval/Abstract.hs | 2 +- .../src/{ => lib}/Calligraphy/Brushes.hs | 0 brush-strokes/src/{ => lib}/Debug/Utils.hs | 2 + .../src/{ => lib}/Math/Algebra/Dual.hs | 0 .../{ => lib}/Math/Algebra/Dual/Internal.hs | 66 +- .../src/{ => lib}/Math/Bezier/Cubic.hs | 13 +- .../src/{ => lib}/Math/Bezier/Cubic/Fit.hs | 0 .../src/{ => lib}/Math/Bezier/Quadratic.hs | 13 +- .../src/{ => lib}/Math/Bezier/Spline.hs | 0 .../src/{ => lib}/Math/Bezier/Stroke.hs | 320 +++++++-- .../Math/Bezier/Stroke/EnvelopeEquation.hs | 28 + .../src/{ => lib}/Math/Differentiable.hs | 0 brush-strokes/src/{ => lib}/Math/Epsilon.hs | 0 brush-strokes/src/{ => lib}/Math/Interval.hs | 0 .../src/{ => lib}/Math/Interval/FMA.hs | 0 .../src/{ => lib}/Math/Interval/Internal.hs | 14 +- brush-strokes/src/{ => lib}/Math/Linear.hs | 0 .../src/{ => lib}/Math/Linear/Internal.hs | 0 .../src/{ => lib}/Math/Linear/Solve.hs | 0 brush-strokes/src/{ => lib}/Math/Module.hs | 0 .../src/{ => lib}/Math/Module/Internal.hs | 0 brush-strokes/src/{ => lib}/Math/Monomial.hs | 0 .../src/{ => lib}/Math/Orientation.hs | 0 brush-strokes/src/{ => lib}/Math/Ring.hs | 0 brush-strokes/src/{ => lib}/Math/Roots.hs | 0 brush-strokes/src/{ => lib}/TH/Utils.hs | 0 .../src/metafont}/Main.hs | 0 .../metafont}/MetaBrush/MetaFont/Convert.hs | 0 src/app/MetaBrush/Application.hs | 12 +- src/metabrushes/MetaBrush/Records.hs | 8 +- 35 files changed, 1021 insertions(+), 557 deletions(-) delete mode 100644 brush-strokes/bench/Main.hs create mode 100644 brush-strokes/src/cusps/bench/Main.hs rename {src/cusps => brush-strokes/src/cusps/inspect}/Main.hs (99%) rename {src/cusps => brush-strokes/src/cusps/inspect}/Math/Interval/Abstract.hs (99%) rename brush-strokes/src/{ => lib}/Calligraphy/Brushes.hs (100%) rename brush-strokes/src/{ => lib}/Debug/Utils.hs (98%) rename brush-strokes/src/{ => lib}/Math/Algebra/Dual.hs (100%) rename brush-strokes/src/{ => lib}/Math/Algebra/Dual/Internal.hs (95%) rename brush-strokes/src/{ => lib}/Math/Bezier/Cubic.hs (93%) rename brush-strokes/src/{ => lib}/Math/Bezier/Cubic/Fit.hs (100%) rename brush-strokes/src/{ => lib}/Math/Bezier/Quadratic.hs (90%) rename brush-strokes/src/{ => lib}/Math/Bezier/Spline.hs (100%) rename brush-strokes/src/{ => lib}/Math/Bezier/Stroke.hs (82%) rename brush-strokes/src/{ => lib}/Math/Bezier/Stroke/EnvelopeEquation.hs (89%) rename brush-strokes/src/{ => lib}/Math/Differentiable.hs (100%) rename brush-strokes/src/{ => lib}/Math/Epsilon.hs (100%) rename brush-strokes/src/{ => lib}/Math/Interval.hs (100%) rename brush-strokes/src/{ => lib}/Math/Interval/FMA.hs (100%) rename brush-strokes/src/{ => lib}/Math/Interval/Internal.hs (94%) rename brush-strokes/src/{ => lib}/Math/Linear.hs (100%) rename brush-strokes/src/{ => lib}/Math/Linear/Internal.hs (100%) rename brush-strokes/src/{ => lib}/Math/Linear/Solve.hs (100%) rename brush-strokes/src/{ => lib}/Math/Module.hs (100%) rename brush-strokes/src/{ => lib}/Math/Module/Internal.hs (100%) rename brush-strokes/src/{ => lib}/Math/Monomial.hs (100%) rename brush-strokes/src/{ => lib}/Math/Orientation.hs (100%) rename brush-strokes/src/{ => lib}/Math/Ring.hs (100%) rename brush-strokes/src/{ => lib}/Math/Roots.hs (100%) rename brush-strokes/src/{ => lib}/TH/Utils.hs (100%) rename {src/convert => brush-strokes/src/metafont}/Main.hs (100%) rename {src/convert => brush-strokes/src/metafont}/MetaBrush/MetaFont/Convert.hs (100%) diff --git a/MetaBrush.cabal b/MetaBrush.cabal index 3f9edf9..641b871 100644 --- a/MetaBrush.cabal +++ b/MetaBrush.cabal @@ -187,49 +187,6 @@ library metabrushes , bytestring >= 0.10.10.0 && < 0.12 -executable cusps - - import: - common - - hs-source-dirs: - src/cusps - - default-language: - Haskell2010 - - main-is: - Main.hs - - other-modules: - Math.Interval.Abstract - - -executable convert-metafont - - import: - common, extras - - hs-source-dirs: - src/convert - - default-language: - Haskell2010 - - main-is: - Main.hs - - other-modules: - MetaBrush.MetaFont.Convert - - build-depends: - metabrushes, - diagrams-contrib, - diagrams-lib, - linear, - parsec - - executable MetaBrush import: diff --git a/brush-strokes/bench/Main.hs b/brush-strokes/bench/Main.hs deleted file mode 100644 index 929fb61..0000000 --- a/brush-strokes/bench/Main.hs +++ /dev/null @@ -1,370 +0,0 @@ -{-# LANGUAGE PolyKinds #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE UndecidableInstances #-} - -module Main - ( main - - -- Testing - , TestCase(..) - , testCases - , testCaseStrokeFunctions - , eval - , mkVal, mkBox - , potentialCusp - , dEdsdcdt - ) - where - --- base -import Control.Concurrent.MVar - ( newMVar ) -import Data.Coerce - ( coerce ) -import Data.Foldable - ( for_ ) -import Data.List - ( intercalate ) -import GHC.Exts - ( Proxy#, proxy# ) -import GHC.Generics - ( Generic ) -import GHC.TypeNats - ( type (-) ) -import Numeric - ( showFFloat ) - --- containers -import Data.Sequence - ( Seq ) -import qualified Data.Sequence as Seq - ( index ) -import Data.Tree - ( foldTree ) - --- brush-strokes -import Calligraphy.Brushes -import Debug.Utils - ( logToFile ) -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 - ] - --logFileMVar <- newMVar "logs/trickyCusp.log" - --logToFile logFileMVar (unlines logLines) - -- `seq` return () - -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]) - -eval fI $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799) -> D12 {_D12_v = T[ℝ2 -10088.6674944889 -3281.3820867312834, ℝ2 4124.668381545453 4524.807156085763], _D12_dx = TT[ℝ2 -173746.97965005718 -33281.18494907289, ℝ2 298.2609121556852 23639.772884799597], _D12_dy = TT[ℝ2 -18454.27716258352 -28337.509817580823, ℝ2 1163.6949532017436 -13936.383137525536]}} -i.e. -> f = [ℝ2 -10088.6674944889 -3281.3820867312834, ℝ2 4124.668381545453 4524.807156085763] -> f_t = [ℝ2 -173746.97965005718 -33281.18494907289, ℝ2 298.2609121556852 23639.772884799597] -> f_s = [ℝ2 -18454.27716258352 -28337.509817580823, ℝ2 1163.6949532017436 -13936.383137525536] - -(f, fI) = testCaseStrokeFunctions trickyCusp2 -t = 𝕀 (ℝ1 0.5798) (ℝ1 0.675) -s = 𝕀 (ℝ1 0.26798) (ℝ1 0.26799) -t_mid = 0.5 * ( 0.5798 + 0.675 ) -s_mid = 0.5 * ( 0.26798 + 0.26799 ) -D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) = dEdsdcdt $ eval fI (t, 3, s) -t' = coerce ( (-) @( 𝕀 Double ) ) t ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1 -s' = coerce ( (-) @( 𝕀 Double ) ) s ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1 -a = ( f_t, f_s ) -b = negV2 $ singleton $ midV2 f -[((t2', s2'), isContr)] = gaussSeidel a b (t', s') -t2 = coerce ( (+) @( 𝕀 Double ) ) t2' ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1 -s2 = coerce ( (+) @( 𝕀 Double ) ) s2' ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1 - -t2 -> [ℝ1 0.6102365832093095, ℝ1 0.6750000000000002] -s2 -> [ℝ1 0.26798, ℝ1 0.26799000000000006] - - -let ( 𝕀 ( ℝ2 a11_lo a21_lo ) ( ℝ2 a11_hi a21_hi ), 𝕀 ( ℝ2 a12_lo a22_lo ) ( ℝ2 a12_hi a22_hi ) ) = a -let ( 𝕀 ( ℝ2 b1_lo b2_lo ) ( ℝ2 b1_hi b2_hi ) ) = b -let ( 𝕀 ( ℝ1 x1_lo ) ( ℝ1 x1_hi ), 𝕀 ( ℝ1 x2_lo ) ( ℝ1 x2_hi ) ) = ( t', s' ) - -a11 = 𝕀 a11_lo a11_hi -a12 = 𝕀 a12_lo a12_hi -a21 = 𝕀 a21_lo a21_hi -a22 = 𝕀 a22_lo a22_hi -b1 = 𝕀 b1_lo b1_hi -b2 = 𝕀 b2_lo b2_hi -x1 = 𝕀 x1_lo x1_hi -x2 = 𝕀 x2_lo x2_hi - -( b1 - a12 * x2 ) -> [2981.90728508591, 2982.0918278575364] - -extendedRecip a11 - --} - -negV2 :: 𝕀ℝ 2 -> 𝕀ℝ 2 -negV2 ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = - let !( 𝕀 x'_lo x'_hi ) = negate $ 𝕀 x_lo x_hi - !( 𝕀 y'_lo y'_hi ) = negate $ 𝕀 y_lo y_hi - in 𝕀 ( ℝ2 x'_lo y'_lo ) ( ℝ2 x'_hi y'_hi ) - -midV2 :: 𝕀ℝ 2 -> ℝ 2 -midV2 ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = - ℝ2 ( 0.5 * ( x_lo + x_hi ) ) ( 0.5 * ( y_lo + y_hi ) ) - -logLines :: [ String ] -logLines = - [ "f = dE/ds * dc/dt: f, df/dt, df/ds" - , "{" ++ - (intercalate "," - [ "{" ++ showD (midPoint t) ++ "," ++ showD (midPoint s) ++ ",{" ++ intercalate "," vals ++ "}}" - | t <- map ( around 0.5798 ) [-0.05, -0.049.. 0.05] - , let i = 3 - , s <- map ( around 0.26798 ) [-0.05, -0.049.. 0.05] - , let StrokeDatum - { 𝛿E𝛿sdcdt = D12 (T f) (T (T f_t)) (T (T f_s)) - } = (curvesI t `Seq.index` i) s - ℝ2 vx vy = midPoint2 f - ℝ2 vx_t vy_t = midPoint2 f_t - ℝ2 vx_s vy_s = midPoint2 f_s - vals = [ "{" ++ showD vx ++ "," ++ showD vy ++ "}" - , "{" ++ showD vx_t ++ "," ++ showD vy_t ++ "}" - , "{" ++ showD vx_s ++ "," ++ showD vy_s ++ "}" - ] - ] - ) ++ "}" - ] - where - around :: Double -> Double -> 𝕀ℝ 1 - around z0 z = 𝕀 ( ℝ1 ( z + z0 - 1e-6 ) ) ( ℝ1 ( z + z0 + 1e-6 ) ) - ( _, curvesI ) = testCaseStrokeFunctions trickyCusp2 - midPoint (𝕀 (ℝ1 lo) (ℝ1 hi)) = 0.5 * ( lo + hi ) - midPoint2 (𝕀 (ℝ2 lo_x lo_y) (ℝ2 hi_x hi_y)) - = ℝ2 ( 0.5 * ( lo_x + hi_x ) ) ( 0.5 * ( lo_y + hi_y ) ) - -showD :: Double -> String -showD float = showFFloat (Just 6) float "" - --------------------------------------------------------------------------------- - -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 e825794..9ff3b55 100644 --- a/brush-strokes/brush-strokes.cabal +++ b/brush-strokes/brush-strokes.cabal @@ -83,13 +83,23 @@ common common base >= 4.19 +common extra + + build-depends: + acts + ^>= 0.3.1.0 + , generic-lens + >= 2.2 && < 2.3 + , groups + ^>= 0.5.3 + library import: - common + common, extra hs-source-dirs: - src + src/lib default-language: Haskell2010 @@ -129,9 +139,6 @@ library build-depends: template-haskell >= 2.18 && < 2.22 - - , acts - ^>= 0.3.1.0 , bifunctors >= 5.5.4 && < 5.7 , code-page @@ -144,10 +151,6 @@ library ^>= 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 @@ -161,13 +164,58 @@ library , transformers >= 0.5.6.2 && < 0.7 +--executable convert-metafont +-- +-- import: +-- common +-- +-- hs-source-dirs: +-- src/metafont +-- +-- default-language: +-- Haskell2010 +-- +-- main-is: +-- Main.hs +-- +-- other-modules: +-- Calligraphy.MetaFont.Convert +-- +-- build-depends: +-- diagrams-contrib, +-- diagrams-lib, +-- linear, +-- parsec + +executable inspect + + import: + common, extra + + hs-source-dirs: + src/cusps/inspect + + default-language: + Haskell2010 + + main-is: + Main.hs + + other-modules: + Math.Interval.Abstract + + build-depends: + brush-strokes, + data-reify + ^>= 0.6.3 + benchmark cusps import: common hs-source-dirs: - bench + src/cusps/bench main-is: Main.hs diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs new file mode 100644 index 0000000..36e793c --- /dev/null +++ b/brush-strokes/src/cusps/bench/Main.hs @@ -0,0 +1,616 @@ +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE UndecidableInstances #-} + +module Main + ( main + + -- Testing + , TestCase(..) + , testCases + , BrushStroke(..) + , brushStrokeFunctions + , eval + , mkVal, mkBox + , potentialCusp + , dEdsdcdt + ) + where + +-- base +import Control.Concurrent.MVar + ( newMVar ) +import Data.Coerce + ( coerce ) +import Data.Foldable + ( for_ ) +import Data.List + ( intercalate ) +import GHC.Exts + ( Proxy#, proxy# ) +import GHC.Generics + ( Generic ) +import GHC.TypeNats + ( type (-) ) +import Numeric + ( showFFloat ) + +-- containers +import Data.Sequence + ( Seq ) +import qualified Data.Sequence as Seq + ( index ) +import Data.Tree + ( foldTree ) + +-- tree-view +import Data.Tree.View + ( showTree ) + +-- brush-strokes +import Calligraphy.Brushes +import Debug.Utils + ( logToFile ) +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 = do + for_ testCases $ \ testCase@( TestCase { testName, testBrushStroke, testAlgorithmParams, testStartBoxes } ) -> do + let ( _, testStrokeFnI ) = brushStrokeFunctions testBrushStroke + ( newtTrees, ( dunno, sols ) ) = computeCusps testAlgorithmParams testStrokeFnI testStartBoxes + showedTrees = map ( uncurry showIntervalNewtonTree ) newtTrees + testHeader = + [ "", "Test case '" ++ testName ++ "':" ] + putStrLn $ unlines $ + testHeader ++ + 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 + ] + -- logFileMVar <- newMVar "logs/fnData.log" + -- logToFile logFileMVar (unlines logLines) + -- `seq` return () + +testCases :: [ TestCase ] +testCases = benchCases + -- [ --trickyCusp2 + -- ellipse "full" (0,1) pi $ defaultStartBoxes [ 2 ] + -- ] + -- ++ +{- + [ ellipse ( "(k1, k2) = " ++ show (k1, k2) ) (k1, k2) pi $ defaultStartBoxes [ 2 ] + | (k1, k2) <- + [(0.5,0.6), (0.55, 0.56)] + ] ++ + [ ellipse ( "'(k1, k2) = " ++ show (k1, k2) ) (0,1) pi [ mkBox (k1 + zero, k2 + zero) 2 (zero, one) ] + | (k1, k2) <- + [(0.5,0.6), (0.55, 0.56)] + ] +-} + +benchCases :: [ TestCase ] +benchCases = [ ellipseTestCase "full" ( 0, 1 ) pi $ defaultStartBoxes [ 0 .. 3 ] ] + +-------------------------------------------------------------------------------- + +data BrushStroke = + forall nbParams. ParamsCt nbParams => + BrushStroke + { brush :: !( Brush nbParams ) + , stroke :: !( Point nbParams, Curve Open () ( Point nbParams ) ) + } + + +data TestCase = + TestCase + { testName :: String + , testBrushStroke :: BrushStroke + , testAlgorithmParams :: CuspAlgorithmParams + , testStartBoxes :: [ Box ] + } + +brushStrokeFunctions + :: BrushStroke + -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) + , 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) +brushStrokeFunctions ( BrushStroke { stroke = ( sp0, crv ), brush } ) = + getStrokeFunctions brush 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 ) ) + +zero, one :: Double +zero = 1e-6 +one = 1 - zero +{-# INLINE zero #-} +{-# INLINE one #-} + +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]) + +eval fI $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799) +> D12 {_D12_v = T[ℝ2 -10088.6674944889 -3281.3820867312834, ℝ2 4124.668381545453 4524.807156085763], _D12_dx = TT[ℝ2 -173746.97965005718 -33281.18494907289, ℝ2 298.2609121556852 23639.772884799597], _D12_dy = TT[ℝ2 -18454.27716258352 -28337.509817580823, ℝ2 1163.6949532017436 -13936.383137525536]}} +i.e. +> f = [ℝ2 -10088.6674944889 -3281.3820867312834, ℝ2 4124.668381545453 4524.807156085763] +> f_t = [ℝ2 -173746.97965005718 -33281.18494907289, ℝ2 298.2609121556852 23639.772884799597] +> f_s = [ℝ2 -18454.27716258352 -28337.509817580823, ℝ2 1163.6949532017436 -13936.383137525536] + +(f, fI) = testCaseStrokeFunctions trickyCusp2 +t = 𝕀 (ℝ1 0.5798) (ℝ1 0.675) +s = 𝕀 (ℝ1 0.26798) (ℝ1 0.26799) +t_mid = 0.5 * ( 0.5798 + 0.675 ) +s_mid = 0.5 * ( 0.26798 + 0.26799 ) +D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) = dEdsdcdt $ eval fI (t, 3, s) +t' = coerce ( (-) @( 𝕀 Double ) ) t ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1 +s' = coerce ( (-) @( 𝕀 Double ) ) s ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1 +a = ( f_t, f_s ) +b = negV2 $ singleton $ midV2 f +[((t2', s2'), isContr)] = gaussSeidel a b (t', s') +t2 = coerce ( (+) @( 𝕀 Double ) ) t2' ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1 +s2 = coerce ( (+) @( 𝕀 Double ) ) s2' ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1 + +t2 +> [ℝ1 0.6102365832093095, ℝ1 0.6750000000000002] +s2 +> [ℝ1 0.26798, ℝ1 0.26799000000000006] + + +let ( 𝕀 ( ℝ2 a11_lo a21_lo ) ( ℝ2 a11_hi a21_hi ), 𝕀 ( ℝ2 a12_lo a22_lo ) ( ℝ2 a12_hi a22_hi ) ) = a +let ( 𝕀 ( ℝ2 b1_lo b2_lo ) ( ℝ2 b1_hi b2_hi ) ) = b +let ( 𝕀 ( ℝ1 x1_lo ) ( ℝ1 x1_hi ), 𝕀 ( ℝ1 x2_lo ) ( ℝ1 x2_hi ) ) = ( t', s' ) + +a11 = 𝕀 a11_lo a11_hi +a12 = 𝕀 a12_lo a12_hi +a21 = 𝕀 a21_lo a21_hi +a22 = 𝕀 a22_lo a22_hi +b1 = 𝕀 b1_lo b1_hi +b2 = 𝕀 b2_lo b2_hi +x1 = 𝕀 x1_lo x1_hi +x2 = 𝕀 x2_lo x2_hi + +( b1 - a12 * x2 ) +> [2981.90728508591, 2982.0918278575364] + +extendedRecip a11 + +----------------------- + +fI k rot = snd $ testCaseStrokeFunctions $ ellipse k rot +d k rot = width $ _D22_v $ ee $ eval (fI k rot) $ mkBox (1e-6, 1-1e-6) 3 (1e-6, 1-1e-6) + +----------------------- + +(f, fI) = testCaseStrokeFunctions $ ellipse 1 pi +t = 𝕀 (ℝ1 0.001) (ℝ1 0.099) +s = 𝕀 (ℝ1 0.001) (ℝ1 0.999) +t_mid = 0.5 * ( 0.001 + 0.099 ) +s_mid = 0.5 * ( 0.001 + 0.999 ) +D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) = dEdsdcdt $ eval fI (t, 3, s) +t' = coerce ( (-) @( 𝕀 Double ) ) t ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1 +s' = coerce ( (-) @( 𝕀 Double ) ) s ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1 +a = ( f_t, f_s ) +b = negV2 $ singleton $ midV2 f +[((t2', s2'), isContr)] = gaussSeidel a b (t', s') +t2 = coerce ( (+) @( 𝕀 Double ) ) t2' ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1 +s2 = coerce ( (+) @( 𝕀 Double ) ) s2' ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1 + +> t = [0.001, 0.099] +> t' = [-0.049, 0.049] + +(f, fI) = testCaseStrokeFunctions $ ellipse 0.1 pi +t = 𝕀 (ℝ1 0.001) (ℝ1 0.999) +s = 𝕀 (ℝ1 0.001) (ℝ1 0.999) +t_mid = 0.5 * ( 0.001 + 0.999 ) +s_mid = 0.5 * ( 0.001 + 0.999 ) +D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) = dEdsdcdt $ eval fI (t, 3, s) +t' = coerce ( (-) @( 𝕀 Double ) ) t ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1 +s' = coerce ( (-) @( 𝕀 Double ) ) s ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1 +a = ( f_t, f_s ) +b = negV2 $ singleton $ midV2 f +[((t2', s2'), isContr)] = gaussSeidel a b (t', s') +t2 = coerce ( (+) @( 𝕀 Double ) ) t2' ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1 +s2 = coerce ( (+) @( 𝕀 Double ) ) s2' ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1 + +> t = [0.001, 0.999] +> t' = [-0.499, 0.499] + +-} + +width :: 𝕀ℝ 1 -> Double +width (𝕀 (ℝ1 lo) (ℝ1 hi)) = hi - lo + +negV2 :: 𝕀ℝ 2 -> 𝕀ℝ 2 +negV2 ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = + let !( 𝕀 x'_lo x'_hi ) = negate $ 𝕀 x_lo x_hi + !( 𝕀 y'_lo y'_hi ) = negate $ 𝕀 y_lo y_hi + in 𝕀 ( ℝ2 x'_lo y'_lo ) ( ℝ2 x'_hi y'_hi ) + +midV2 :: 𝕀ℝ 2 -> ℝ 2 +midV2 ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = + ℝ2 ( 0.5 * ( x_lo + x_hi ) ) ( 0.5 * ( y_lo + y_hi ) ) + + +logLines :: [ String ] +logLines = + [ "E, dE/ds * dc/dt" + , "{" ++ + (intercalate "," + [ "{" ++ showD t ++ "," ++ showD s ++ ",{" ++ intercalate "," vals ++ "}}" + | t <- [ 0.5484, 0.5484 + 0.00001 .. 0.5488 ] + , s <- [ 0.5479, 0.5479 + 0.00001 .. 0.5483 ] + , let StrokeDatum + { ee = D22 ee _ _ _ _ _ + , 𝛿E𝛿sdcdt = D12 (T f) (T (T f_t)) (T (T f_s)) + } = (curvesI (singleton (ℝ1 t)) `Seq.index` i) (singleton (ℝ1 s)) + ℝ2 vx vy = midPoint2 f + --ℝ2 vx_t vy_t = midPoint2 f_t + --ℝ2 vx_s vy_s = midPoint2 f_s + vals = [ showD ( midPoint ee ) + , "{" ++ showD vx ++ "," ++ showD vy ++ "}" + -- , "{" ++ showD vx_t ++ "," ++ showD vy_t ++ "}" + -- , "{" ++ showD vx_s ++ "," ++ showD vy_s ++ "}" + ] + ] + ) ++ "}" + ] + where + i = 2 + ( curves, curvesI ) = brushStrokeFunctions $ ellipseBrushStroke ( 0, 1 ) pi + midPoint (𝕀 (ℝ1 lo) (ℝ1 hi)) = 0.5 * ( lo + hi ) + midPoint2 (𝕀 (ℝ2 lo_x lo_y) (ℝ2 hi_x hi_y)) + = ℝ2 ( 0.5 * ( lo_x + hi_x ) ) ( 0.5 * ( lo_y + hi_y ) ) + + +-- t = 0.5486102, s = 0.5480951 + + +bloo = + [ ( e * e + vx * vx + vy * vy, ( t, s ) ) + | t <- [ 0.548609, 0.548609 + 0.0000001 .. 0.54862 ] + , s <- [ 0.548094, 0.548094 + 0.0000001 .. 0.548096 ] + , let StrokeDatum + { ee = D22 ee _ _ _ _ _ + , 𝛿E𝛿sdcdt = D12 (T f) _ _ + } = (curvesI (singleton (ℝ1 t)) `Seq.index` i) (singleton (ℝ1 s)) + e = midPoint ee + ℝ2 vx vy = midPoint2 f + vals = [ showD ( midPoint ee ) + , "{" ++ showD vx ++ "," ++ showD vy ++ "}" + ] + ] + where + i = 2 + ( curves, curvesI ) = brushStrokeFunctions $ ellipseBrushStroke ( 0, 1 ) pi + midPoint (𝕀 (ℝ1 lo) (ℝ1 hi)) = 0.5 * ( lo + hi ) + midPoint2 (𝕀 (ℝ2 lo_x lo_y) (ℝ2 hi_x hi_y)) + = ℝ2 ( 0.5 * ( lo_x + hi_x ) ) ( 0.5 * ( lo_y + hi_y ) ) + +{- + +(f, fI) = brushStrokeFunctions $ ellipseBrushStroke ( 0, 1 ) pi + + + + +_D22_v $ ee $ eval fI $ mkBox (0.5, 0.6) 2 (0,1) +> [ℝ1 -9531.427315889887, ℝ1 10135.074304695485] + +minimum $ map inf $ [ _D22_v $ ee $ eval fI $ mkBox (t, t + 0.01) 2 (0,1) | t <- [ 0.5, 0.51 .. 0.59 ] ] +> ℝ1 -5718.905635365308 + +maximum $ map sup $ [ _D22_v $ ee $ eval fI $ mkBox (t, t + 0.01) 2 (0,1) | t <- [ 0.5, 0.51 .. 0.59 ] ] +> ℝ1 5099.008191092755 + +minimum $ map inf $ [ _D22_v $ ee $ eval fI $ mkBox (t, t) 2 (s, s) | t <- [ 0.5, 0.501 .. 0.6 ], s <- [ zero, zero + 0.001 .. one ] ] +> ℝ1 -675.9595496147458 + +maximum $ map sup $ [ _D22_v $ ee $ eval fI $ mkBox (t, t) 2 (s, s) | t <- [ 0.5, 0.501 .. 0.6 ], s <- [ zero, zero + 0.001 .. one ] ] +> ℝ1 2401.9644509525997 + +_D12_v $ dEdsdcdt $ eval fI $ mkBox (0.5, 0.6) 2 (0,1) +> T[ℝ2 -1.7300637136531524e7 -1.262824151868635e7, ℝ2 1.632868898735965e7 1.1869759856947478e7] + +minimum [ _x $ inf $ unT $ _D12_v $ dEdsdcdt $ eval fI $ mkBox (t, t + 0.01) 2 (0,1) | t <- [ 0.5, 0.51 .. 0.59 ] ] +-5606615.948203902 + +maximum [ _x $ sup $ unT $ _D12_v $ dEdsdcdt $ eval fI $ mkBox (t, t + 0.01) 2 (0,1) | t <- [ 0.5, 0.51 .. 0.59 ] ] +4340858.832347277 + +minimum [ _x $ inf $ unT $ _D12_v $ dEdsdcdt $ eval fI $ mkBox (t, t) 2 (s, s) | t <- [ 0.5, 0.501 .. 0.6 ], s <- [ zero, zero + 0.001 .. one ] ] +-1785730.2396688666 + +maximum [ _x $ sup $ unT $ _D12_v $ dEdsdcdt $ eval fI $ mkBox (t, t) 2 (s, s) | t <- [ 0.5, 0.501 .. 0.6 ], s <- [ zero, zero + 0.001 .. one ] ] +974842.6547409865 + +maximum [ _y $ sup $ unT $ _D12_v $ dEdsdcdt $ eval fI $ mkBox (t, t) 2 (s, s) | t <- [ 0.5, 0.501 .. 0.6 ], s <- [ zero, zero + 0.001 .. one ] ] +845211.4833711373 + + + +let showTrees b = map ( uncurry showIntervalNewtonTree ) $ fst $ intervalNewtonGSFrom NoPreconditioning 1e-7 fI b + +putStrLn $ unlines $ map Data.Tree.View.showTree $ showTrees $ mkBox (0.5486101933245248, 0.5486102071493622) 2 (0.548095036738487, 0.5480952) + + + +([ℝ1 0.5486101960904595, ℝ1 0.5486102071493623],2,[ℝ1 0.5480950771755867, ℝ1 0.5480952000000001]) + +-} + +_x ( ℝ2 x _ ) = x +_y ( ℝ2 _ y ) = y + +showD :: Double -> String +showD float = showFFloat (Just 6) float "" + +-------------------------------------------------------------------------------- + +ellipseTestCase :: String -> ( Double, Double ) -> Double -> [ Box ] -> TestCase +ellipseTestCase str k0k1 rot startBoxes = + TestCase + { testName = "ellipse (" ++ str ++ ")" + , testBrushStroke = ellipseBrushStroke k0k1 rot + , testAlgorithmParams = + CuspAlgorithmParams + { preconditioning = InverseMidJacobian + , maxWidth = 1e-7 + } + , testStartBoxes = startBoxes + } + +ellipseBrushStroke :: ( Double, Double ) -> Double -> BrushStroke +ellipseBrushStroke ( k0, k1 ) rot = + BrushStroke + { brush = ellipseBrush + , stroke = ( p0, LineTo ( NextPoint p1 ) () ) } + where + mkPt x y w h phi = + Point + { pointCoords = ℝ2 x y + , pointParams = Params $ ℝ3 w h phi + } + l k = lerp @( T Double ) k + p k = mkPt ( l k 0 100 ) 0 ( l k 10 15 ) ( l k 25 40 ) ( l k 0 rot ) + p0 = p k0 + p1 = p k1 + +trickyCusp2TestCase :: TestCase +trickyCusp2TestCase = + TestCase + { testName = "trickyCusp2" + , testBrushStroke = trickyCusp2BrushStroke + , testAlgorithmParams = + CuspAlgorithmParams + { preconditioning = InverseMidJacobian + , maxWidth = 1e-7 + } + , testStartBoxes = defaultStartBoxes [ 0 .. 3 ] + } + +trickyCusp2BrushStroke :: BrushStroke +trickyCusp2BrushStroke = + BrushStroke + { brush = circleBrush + , stroke = ( p0, Bezier3To p1 p2 ( NextPoint p3 ) () ) + } + 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 ] + -> ( [ ( Box, IntervalNewtonTree Box ) ], ( [ Box ], [ Box ] ) ) +computeCusps params eqs startBoxes = + foldMap + ( intervalNewtonGSFrom ( preconditioning params ) ( maxWidth params ) eqs ) + startBoxes + +defaultStartBoxes :: [ Int ] -> [ Box ] +defaultStartBoxes is = + [ mkBox (zero, one) i (zero, one) | i <- is ] + +getR1 (ℝ1 u) = u + +{- + +(f, fI) = brushStrokeFunctions $ ellipseBrushStroke (0,1) pi +nbPotentialSols box = let ( _newtTrees, ( dunno, sols ) ) = intervalNewtonGSFrom NoPreconditioning 1e-7 fI box in length dunno + length sols +showTrees box = map ( uncurry showIntervalNewtonTree ) $ fst $ intervalNewtonGSFrom NoPreconditioning 1e-7 fI box + +(t, i, s) = mkBox (0.548610200176363, 0.5486102071493623) 2 (0.5480950215354709, 0.5480952) +putStrLn $ unlines $ map Data.Tree.View.showTree $ showTrees (t,i,s) + + +t_mid = 0.5 * ( getR1 ( inf t ) + getR1 ( sup t ) ) +s_mid = 0.5 * ( getR1 ( inf s ) + getR1 ( sup s ) ) +D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) = dEdsdcdt $ eval fI (t, i, s) +t' = coerce ( (-) @( 𝕀 Double ) ) t ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1 +s' = coerce ( (-) @( 𝕀 Double ) ) s ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1 +a = ( f_t, f_s ) +b = negV2 $ singleton $ midV2 f +[((t2', s2'), isContr)] = gaussSeidel a b (t', s') +t2 = coerce ( (+) @( 𝕀 Double ) ) t2' ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1 +s2 = coerce ( (+) @( 𝕀 Double ) ) s2' ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1 + +t2 +> [ℝ1 0.548610200176363, ℝ1 0.5486102071493624] +s2 +> [ℝ1 0.5480950911334656, ℝ1 0.5480952000000001] + +mkBox (0.548610200176363, 0.5486102071493624) i (0.5480950911334656, 0.5480952000000001) + +t inf (no change) + +t sup (no change) + +s_inf: +0.5480950215354709 +0.5480950911334656 + +s_sup (no change) + +ghci> potentialCusp $ eval fI $ mkBox (0.548610200176363, 0.5486102071493623) 2 (0.54809502, 0.5480952) +True +ghci> potentialCusp $ eval fI $ mkBox (0.548610200176363, 0.5486102071493623) 2 (0.54809503, 0.5480952) +False + + + +let ( 𝕀 ( ℝ2 a11_lo a21_lo ) ( ℝ2 a11_hi a21_hi ), 𝕀 ( ℝ2 a12_lo a22_lo ) ( ℝ2 a12_hi a22_hi ) ) = a +let ( 𝕀 ( ℝ2 b1_lo b2_lo ) ( ℝ2 b1_hi b2_hi ) ) = b +let ( 𝕀 ( ℝ1 x1_lo ) ( ℝ1 x1_hi ), 𝕀 ( ℝ1 x2_lo ) ( ℝ1 x2_hi ) ) = ( t', s' ) + +a11 = 𝕀 a11_lo a11_hi +a12 = 𝕀 a12_lo a12_hi +a21 = 𝕀 a21_lo a21_hi +a22 = 𝕀 a22_lo a22_hi +b1 = 𝕀 b1_lo b1_hi +b2 = 𝕀 b2_lo b2_hi +x1 = 𝕀 x1_lo x1_hi +x2 = 𝕀 x2_lo x2_hi + +-} \ No newline at end of file diff --git a/src/cusps/Main.hs b/brush-strokes/src/cusps/inspect/Main.hs similarity index 99% rename from src/cusps/Main.hs rename to brush-strokes/src/cusps/inspect/Main.hs index ee71b18..a876c12 100644 --- a/src/cusps/Main.hs +++ b/brush-strokes/src/cusps/inspect/Main.hs @@ -32,7 +32,7 @@ import qualified Data.Sequence as Seq import Data.Generics.Product.Typed ( HasType ) --- MetaBrush +-- brush-strokes import Math.Algebra.Dual import qualified Math.Bezier.Quadratic as Quadratic import qualified Math.Bezier.Cubic as Cubic @@ -58,6 +58,7 @@ import Math.Ring , ifThenElse ) +-- brush-strokes:inspect-cusps import Math.Interval.Abstract -------------------------------------------------------------------------------- diff --git a/src/cusps/Math/Interval/Abstract.hs b/brush-strokes/src/cusps/inspect/Math/Interval/Abstract.hs similarity index 99% rename from src/cusps/Math/Interval/Abstract.hs rename to brush-strokes/src/cusps/inspect/Math/Interval/Abstract.hs index 0d54cb4..2d4baa6 100644 --- a/src/cusps/Math/Interval/Abstract.hs +++ b/brush-strokes/src/cusps/inspect/Math/Interval/Abstract.hs @@ -43,7 +43,7 @@ import qualified Data.Map.Strict as Map import Data.Group ( Group(..) ) --- MetaBrush +-- brush-strokes import Math.Algebra.Dual import qualified Math.Bezier.Quadratic as Quadratic import qualified Math.Bezier.Cubic as Cubic diff --git a/brush-strokes/src/Calligraphy/Brushes.hs b/brush-strokes/src/lib/Calligraphy/Brushes.hs similarity index 100% rename from brush-strokes/src/Calligraphy/Brushes.hs rename to brush-strokes/src/lib/Calligraphy/Brushes.hs diff --git a/brush-strokes/src/Debug/Utils.hs b/brush-strokes/src/lib/Debug/Utils.hs similarity index 98% rename from brush-strokes/src/Debug/Utils.hs rename to brush-strokes/src/lib/Debug/Utils.hs index 1687a61..0b32b91 100644 --- a/brush-strokes/src/Debug/Utils.hs +++ b/brush-strokes/src/lib/Debug/Utils.hs @@ -66,3 +66,5 @@ logToFile logFileMVar logContents = FilePath.takeDirectory logFile appendFile logFile logContentsWithHeader return logFile +{-# NOINLINE logToFile #-} + diff --git a/brush-strokes/src/Math/Algebra/Dual.hs b/brush-strokes/src/lib/Math/Algebra/Dual.hs similarity index 100% rename from brush-strokes/src/Math/Algebra/Dual.hs rename to brush-strokes/src/lib/Math/Algebra/Dual.hs diff --git a/brush-strokes/src/Math/Algebra/Dual/Internal.hs b/brush-strokes/src/lib/Math/Algebra/Dual/Internal.hs similarity index 95% rename from brush-strokes/src/Math/Algebra/Dual/Internal.hs rename to brush-strokes/src/lib/Math/Algebra/Dual/Internal.hs index 5b2a1cd..181ec9f 100644 --- a/brush-strokes/src/Math/Algebra/Dual/Internal.hs +++ b/brush-strokes/src/lib/Math/Algebra/Dual/Internal.hs @@ -60,8 +60,8 @@ newtype D𝔸0 v = D0 { _D0_v :: v } -- | \( \mathbb{Z}[x]/(x)^2 \). data D1𝔸1 v = - D11 { _D11_v :: !v - , _D11_dx :: !( T v ) + D11 { _D11_v :: v + , _D11_dx :: ( T v ) } deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) deriving anyclass NFData @@ -70,9 +70,9 @@ data D1𝔸1 v = -- | \( \mathbb{Z}[x]/(x)^3 \). data D2𝔸1 v = - D21 { _D21_v :: !v - , _D21_dx :: !( T v ) - , _D21_dxdx :: !( T v ) + D21 { _D21_v :: v + , _D21_dx :: ( T v ) + , _D21_dxdx :: ( T v ) } deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) deriving anyclass NFData @@ -81,10 +81,10 @@ data D2𝔸1 v = -- | \( \mathbb{Z}[x]/(x)^4 \). data D3𝔸1 v = - D31 { _D31_v :: !v - , _D31_dx :: !( T v ) - , _D31_dxdx :: !( T v ) - , _D31_dxdxdx :: !( T v ) + D31 { _D31_v :: v + , _D31_dx :: ( T v ) + , _D31_dxdx :: ( T v ) + , _D31_dxdxdx :: ( T v ) } deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) deriving anyclass NFData @@ -93,8 +93,8 @@ data D3𝔸1 v = -- | \( \mathbb{Z}[x, y]/(x, y)^2 \). data D1𝔸2 v = - D12 { _D12_v :: !v - , _D12_dx, _D12_dy :: !( T v ) + D12 { _D12_v :: v + , _D12_dx, _D12_dy :: ( T v ) } deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) deriving anyclass NFData @@ -103,9 +103,9 @@ data D1𝔸2 v = -- | \( \mathbb{Z}[x, y]/(x, y)^3 \). data D2𝔸2 v = - D22 { _D22_v :: !v - , _D22_dx, _D22_dy :: !( T v ) - , _D22_dxdx, _D22_dxdy, _D22_dydy :: !( T v ) + D22 { _D22_v :: v + , _D22_dx, _D22_dy :: ( T v ) + , _D22_dxdx, _D22_dxdy, _D22_dydy :: ( T v ) } deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) deriving anyclass NFData @@ -114,10 +114,10 @@ data D2𝔸2 v = -- | \( \mathbb{Z}[x, y]/(x, y)^4 \). data D3𝔸2 v = - D32 { _D32_v :: !v - , _D32_dx, _D32_dy :: !( T v ) - , _D32_dxdx, _D32_dxdy, _D32_dydy :: !( T v ) - , _D32_dxdxdx, _D32_dxdxdy, _D32_dxdydy, _D32_dydydy :: !( T v ) + D32 { _D32_v :: v + , _D32_dx, _D32_dy :: ( T v ) + , _D32_dxdx, _D32_dxdy, _D32_dydy :: ( T v ) + , _D32_dxdxdx, _D32_dxdxdy, _D32_dxdydy, _D32_dydydy :: ( T v ) } deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) deriving anyclass NFData @@ -126,8 +126,8 @@ data D3𝔸2 v = -- | \( \mathbb{Z}[x, y, z]/(x, y, z)^2 \). data D1𝔸3 v = - D13 { _D13_v :: !v - , _D13_dx, _D13_dy, _D13_dz :: !( T v ) + D13 { _D13_v :: v + , _D13_dx, _D13_dy, _D13_dz :: ( T v ) } deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) deriving anyclass NFData @@ -136,8 +136,8 @@ data D1𝔸3 v = -- | \( \mathbb{Z}[x, y, z]/(x, y, z)^3 \). data D2𝔸3 v = - D23 { _D23_v :: !v - , _D23_dx, _D23_dy, _D23_dz :: !( T v ) + D23 { _D23_v :: v + , _D23_dx, _D23_dy, _D23_dz :: ( T v ) , _D23_dxdx, _D23_dxdy, _D23_dydy, _D23_dxdz, _D23_dydz, _D23_dzdz :: !( T v ) } deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) @@ -147,8 +147,8 @@ data D2𝔸3 v = -- | \( \mathbb{Z}[x, y, z]/(x, y, z)^4 \). data D3𝔸3 v = - D33 { _D33_v :: !v - , _D33_dx, _D33_dy, _D33_dz :: !( T v ) + D33 { _D33_v :: v + , _D33_dx, _D33_dy, _D33_dz :: ( T v ) , _D33_dxdx, _D33_dxdy, _D33_dydy, _D33_dxdz, _D33_dydz, _D33_dzdz :: !( T v ) , _D33_dxdxdx, _D33_dxdxdy, _D33_dxdydy, _D33_dydydy , _D33_dxdxdz, _D33_dxdydz, _D33_dxdzdz, _D33_dydydz, _D33_dydzdz, _D33_dzdzdz :: !( T v ) @@ -160,8 +160,8 @@ data D3𝔸3 v = -- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^2 \). data D1𝔸4 v = - D14 { _D14_v :: !v - , _D14_dx, _D14_dy, _D14_dz, _D14_dw :: !( T v ) + D14 { _D14_v :: v + , _D14_dx, _D14_dy, _D14_dz, _D14_dw :: ( T v ) } deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) deriving anyclass NFData @@ -170,10 +170,10 @@ data D1𝔸4 v = -- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^3 \). data D2𝔸4 v = - D24 { _D24_v :: !v - , _D24_dx, _D24_dy, _D24_dz, _D24_dw :: !( T v ) + D24 { _D24_v :: v + , _D24_dx, _D24_dy, _D24_dz, _D24_dw :: ( T v ) , _D24_dxdx, _D24_dxdy, _D24_dydy, _D24_dxdz - , _D24_dydz, _D24_dzdz, _D24_dxdw, _D24_dydw, _D24_dzdw, _D24_dwdw :: !( T v ) + , _D24_dydz, _D24_dzdz, _D24_dxdw, _D24_dydw, _D24_dzdw, _D24_dwdw :: ( T v ) } deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) deriving anyclass NFData @@ -182,14 +182,14 @@ data D2𝔸4 v = -- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^4 \). data D3𝔸4 v = - D34 { _D34_v :: !v - , _D34_dx, _D34_dy, _D34_dz, _D34_dw :: !( T v ) + D34 { _D34_v :: v + , _D34_dx, _D34_dy, _D34_dz, _D34_dw :: ( T v ) , _D34_dxdx, _D34_dxdy, _D34_dydy, _D34_dxdz, _D34_dydz, _D34_dzdz - , _D34_dxdw, _D34_dydw, _D34_dzdw, _D34_dwdw :: !( T v ) + , _D34_dxdw, _D34_dydw, _D34_dzdw, _D34_dwdw :: ( T v ) , _D34_dxdxdx, _D34_dxdxdy, _D34_dxdydy, _D34_dydydy, _D34_dxdxdz, _D34_dxdydz, _D34_dxdzdz, _D34_dydzdz, _D34_dydydz, _D34_dzdzdz, _D34_dxdxdw, _D34_dxdydw, _D34_dydydw, _D34_dxdzdw, _D34_dydzdw, _D34_dzdzdw, - _D34_dxdwdw, _D34_dydwdw, _D34_dzdwdw, _D34_dwdwdw :: !( T v ) + _D34_dxdwdw, _D34_dydwdw, _D34_dzdwdw, _D34_dwdwdw :: ( T v ) } deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) deriving anyclass NFData diff --git a/brush-strokes/src/Math/Bezier/Cubic.hs b/brush-strokes/src/lib/Math/Bezier/Cubic.hs similarity index 93% rename from brush-strokes/src/Math/Bezier/Cubic.hs rename to brush-strokes/src/lib/Math/Bezier/Cubic.hs index 173a01e..ba4fc6d 100644 --- a/brush-strokes/src/Math/Bezier/Cubic.hs +++ b/brush-strokes/src/lib/Math/Bezier/Cubic.hs @@ -8,7 +8,7 @@ module Math.Bezier.Cubic , bezier, bezier', bezier'', bezier''' , derivative , curvature, squaredCurvature, signedCurvature - , subdivide + , subdivide, restrict , ddist, closestPoint , drag, selfIntersectionParameters , extrema @@ -180,6 +180,17 @@ subdivide ( Bezier {..} ) t = ( Bezier p0 q1 q2 pt, Bezier pt r1 r2 p3 ) pt = lerp @v t q2 r1 {-# INLINEABLE subdivide #-} +-- | Restrict a cubic BΓ©zier curve to a sub-interval, re-parametrising +-- to \( [0,1] \). +restrict :: forall v r p. ( Torsor v p, Ring.Field r, Module r v ) => Bezier p -> ( r , r ) -> Bezier p +restrict bez ( a, b ) = fst $ ( flip ( subdivide @v ) b' ) $ snd $ subdivide @v bez a + where + b' = ( b Ring.- a ) Ring./ ( Ring.fromInteger 1 Ring.- a ) + -- TODO: this could be made more efficient. + -- See e.g. "https://math.stackexchange.com/questions/4172835/cubic-b%C3%A9zier-spline-multiple-split" + -- or the paper "On the numerical condition of Bernstein-BΓ©zier subdivision process". +{-# INLINEABLE restrict #-} + -- | Polynomial coefficients of the derivative of the distance to a cubic BΓ©zier curve. ddist :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> p -> [ r ] ddist ( Bezier {..} ) c = [ a5, a4, a3, a2, a1, a0 ] diff --git a/brush-strokes/src/Math/Bezier/Cubic/Fit.hs b/brush-strokes/src/lib/Math/Bezier/Cubic/Fit.hs similarity index 100% rename from brush-strokes/src/Math/Bezier/Cubic/Fit.hs rename to brush-strokes/src/lib/Math/Bezier/Cubic/Fit.hs diff --git a/brush-strokes/src/Math/Bezier/Quadratic.hs b/brush-strokes/src/lib/Math/Bezier/Quadratic.hs similarity index 90% rename from brush-strokes/src/Math/Bezier/Quadratic.hs rename to brush-strokes/src/lib/Math/Bezier/Quadratic.hs index 8f5c8d5..4846036 100644 --- a/brush-strokes/src/Math/Bezier/Quadratic.hs +++ b/brush-strokes/src/lib/Math/Bezier/Quadratic.hs @@ -6,7 +6,7 @@ module Math.Bezier.Quadratic ( Bezier(..) , bezier, bezier', bezier'' , curvature, squaredCurvature, signedCurvature - , subdivide + , subdivide, restrict , ddist, closestPoint , interpolate , extrema @@ -141,6 +141,17 @@ subdivide ( Bezier {..} ) t = ( Bezier p0 q1 pt, Bezier pt r1 p2 ) pt = lerp @v t q1 r1 {-# INLINEABLE subdivide #-} +-- | Restrict a quadratic BΓ©zier curve to a sub-interval, re-parametrising +-- to \( [0,1] \). +restrict :: forall v r p. ( Torsor v p, Ring.Field r, Module r v ) => Bezier p -> ( r , r ) -> Bezier p +restrict bez ( a, b ) = fst $ ( flip ( subdivide @v ) b' ) $ snd $ subdivide @v bez a + where + b' = ( b Ring.- a ) Ring./ ( Ring.fromInteger 1 Ring.- a ) + -- TODO: this could be made more efficient. + -- See e.g. "https://math.stackexchange.com/questions/4172835/cubic-b%C3%A9zier-spline-multiple-split" + -- or the paper "On the numerical condition of Bernstein-BΓ©zier subdivision process". +{-# INLINEABLE restrict #-} + -- | Polynomial coefficients of the derivative of the distance to a quadratic BΓ©zier curve. ddist :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> p -> [ r ] ddist ( Bezier {..} ) c = [ a3, a2, a1, a0 ] diff --git a/brush-strokes/src/Math/Bezier/Spline.hs b/brush-strokes/src/lib/Math/Bezier/Spline.hs similarity index 100% rename from brush-strokes/src/Math/Bezier/Spline.hs rename to brush-strokes/src/lib/Math/Bezier/Spline.hs diff --git a/brush-strokes/src/Math/Bezier/Stroke.hs b/brush-strokes/src/lib/Math/Bezier/Stroke.hs similarity index 82% rename from brush-strokes/src/Math/Bezier/Stroke.hs rename to brush-strokes/src/lib/Math/Bezier/Stroke.hs index 0db8ed5..6113b74 100644 --- a/brush-strokes/src/Math/Bezier/Stroke.hs +++ b/brush-strokes/src/lib/Math/Bezier/Stroke.hs @@ -148,6 +148,7 @@ import Math.Orientation ( Orientation(..), splineOrientation , between ) +import qualified Math.Ring as Ring import Math.Roots import Debug.Utils @@ -600,7 +601,7 @@ outlineFunction rootAlgo ptParams toBrushParams brushFromParams = \ sp0 crv -> ( newtTrees, ( newtDunno, newtSols ) ) = intervalNewtonGS - NoPreconditioning --InverseMidJacobian + InverseMidJacobian 1e-7 curvesI @@ -634,7 +635,7 @@ outlineFunction rootAlgo ptParams toBrushParams brushFromParams = \ sp0 crv -> logContents = unlines $ functionDataLines ++ treeLines in trace (unlines solLines) $ - logToFile cuspFindingMVar logContents `seq` + --logToFile cuspFindingMVar logContents `seq` OutlineInfo { outlineFn = fwdBwd , outlineDefiniteCusps = map ( cuspCoords curves ) newtSols @@ -1112,7 +1113,7 @@ solveEnvelopeEquations rootAlgo _t path_t path'_t ( fwdOffset, bwdOffset ) strok ++ "}" ] - in logToFile rootSolvingMVar logContents `seq` + in --logToFile rootSolvingMVar logContents `seq` ( good, ds, dcdt ) (runSolveMethod, methodName) = case rootAlgo of @@ -1206,7 +1207,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 @@ -1274,6 +1275,115 @@ gaussSeidel return ( ( 𝕀 ( ℝ1 x1'_lo ) ( ℝ1 x1'_hi ), 𝕀 ( ℝ1 x2'_lo ) ( ℝ1 x2'_hi ) ) , sub_x1 && sub_x2 ) +{- +gaussSeidel2 :: Int + -> Double + -> Double + -> ( 𝕀ℝ 2, 𝕀ℝ 2 ) -- ^ columns of \( A \) + -> 𝕀ℝ 2 -- ^ \( B \) + -> ( 𝕀ℝ 1, 𝕀ℝ 1 ) -- ^ initial box \( X \) + -> [ ( ( 𝕀ℝ 1, 𝕀ℝ 1 ), Bool ) ] +gaussSeidel2 maxIters eps_abs eps_rel + ( 𝕀 ( ℝ2 a11_lo a21_lo ) ( ℝ2 a11_hi a21_hi ) + , 𝕀 ( ℝ2 a12_lo a22_lo ) ( ℝ2 a12_hi a22_hi ) ) + ( 𝕀 ( ℝ2 b1_lo b2_lo ) ( ℝ2 b1_hi b2_hi ) ) + x0 + = let !a11 = 𝕀 a11_lo a11_hi + !a12 = 𝕀 a12_lo a12_hi + !a21 = 𝕀 a21_lo a21_hi + !a22 = 𝕀 a22_lo a22_hi + !b1 = 𝕀 b1_lo b1_hi + !b2 = 𝕀 b2_lo b2_hi + + -- See "Algorithm 2" in + -- "Using interval unions to solve linear systems of equations with uncertainties" + iter ( 𝕀 ( ℝ1 x1_lo ) ( ℝ1 x1_hi ), 𝕀 ( ℝ1 x2_lo ) ( ℝ1 x2_hi ) ) = do + let + !x1 = 𝕀 x1_lo x1_hi + !x2 = 𝕀 x2_lo x2_hi + blah1 = do + let s = b1 - a12 * x2 + let s1 = s `monus` ( a11 * x1 ) + y1 <- + if not $ containsZero ( s1 - a11 * x1 ) + then [] + else + if containsZero s1 && containsZero a11 + then [ ( x1, False ) ] + else do + x1'' <- s1 `extendedDivide` a11 + x1'' `intersect` x1 + + let s2 = s `monus` ( a12 * x2 ) + y2 <- + if not $ containsZero ( s2 - a11 * x1 ) + then [] + else + if containsZero s2 && containsZero a11 + then [ ( x1, False ) ] + else do + x1'' <- s2 `extendedDivide` a12 + x1'' `intersect` x1 + return ( y1 `cart` y2 ) + + blah2 = do + let s = b2 - a21 * x1 + let s1 = s `monus` ( a21 * x1 ) + y1 <- + if not $ containsZero ( s1 - a22 * x2 ) + then [] + else + if containsZero s1 && containsZero a22 + then [ ( x1, False ) ] + else do + x1'' <- s1 `extendedDivide` a21 + x1'' `intersect` x1 + + let s2 = s `monus` ( a12 * x2 ) + y2 <- + if not $ containsZero ( s2 - a22 * x2 ) + then [] + else + if containsZero s2 && containsZero a22 + then [ ( x1, False ) ] + else do + x1'' <- s2 `extendedDivide` a22 + x1'' `intersect` x1 + return ( y1 `cart` y2 ) + + blah1 ++ blah2 + + go :: Int -> ( 𝕀ℝ 1, 𝕀ℝ 1 ) -> [ ( ( 𝕀ℝ 1, 𝕀ℝ 1 ), Bool ) ] + go !i x + = do { nxt@( x', sub ) <- iter x + ; let dw_abs = maxWidth x - maxWidth x' + dw_rel = 1 - ( maxWidth x' / maxWidth x ) + ; if sub + || i >= maxIters + || dw_abs < eps_abs + || dw_rel < eps_rel + then return nxt + else go ( i + 1 ) x' + } + + in go 1 x0 + where + maxWidth :: ( 𝕀ℝ 1, 𝕀ℝ 1 ) -> Double + maxWidth ( 𝕀 ( ℝ1 x1_lo ) ( ℝ1 x1_hi ), 𝕀 ( ℝ1 x2_lo ) ( ℝ1 x2_hi ) ) + = max ( x1_hi - x1_lo ) ( x2_hi - x2_lo ) + containsZero :: 𝕀 Double -> Bool + containsZero ( 𝕀 lo hi ) = lo <= 0 && hi >= 0 + monus :: 𝕀 Double -> 𝕀 Double -> 𝕀 Double + monus ( 𝕀 lo1 hi1 ) ( 𝕀 lo2 hi2 ) + | hi1 - lo1 >= hi2 - lo2 + = 𝕀 ( lo1 - lo2 ) ( hi1 - hi2 ) + | otherwise + = 𝕀 ( hi1 - hi2 ) ( lo1 - lo2 ) + cart :: ( 𝕀 Double, Bool ) -> ( 𝕀 Double, Bool ) -> ( ( 𝕀ℝ 1, 𝕀ℝ 1 ), Bool ) + cart ( 𝕀 lo1 hi1, sub1 ) ( 𝕀 lo2 hi2, sub2 ) = + ( ( 𝕀 ( ℝ1 lo1 ) ( ℝ1 hi1 ), 𝕀 ( ℝ1 lo2 ) ( ℝ1 hi2 ) ), sub1 && sub2 ) +-} + -- | Compute the intersection of two intervals. -- -- Returns whether the first interval is a strict subset of the second interval, @@ -1342,7 +1452,7 @@ data IntervalNewtonLeaf d 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 + = Node (show cand ++ " abc " ++ 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 ) ) @@ -1383,96 +1493,152 @@ intervalNewtonGSFrom -> ( [ ( Box, IntervalNewtonTree Box ) ], ( [ Box ], [ Box ] ) ) intervalNewtonGSFrom precondMethod minWidth eqs initBox = runWriter $ - map ( initBox , ) <$> go initBox + map ( initBox , ) <$> evalStrokeDataAndGo initBox where - recur f cands = do - rest <- traverse ( \ c -> do { trees <- go c; return [ (c, tree) | tree <- trees ] } ) cands + recur :: ( cand -> Writer ( [ Box ], [ Box ] ) [ IntervalNewtonTree Box ] ) + -> ( [ ( cand, IntervalNewtonTree Box ) ] -> IntervalNewtonTree Box ) + -> [ cand ] + -> Writer ( [Box], [Box] ) [ IntervalNewtonTree Box ] + recur r f cands = do + rest <- traverse ( \ c -> do { trees <- r c; return [ (c, tree) | tree <- trees ] } ) cands return [ f $ concat rest ] - go :: Box -- box to work on + evalStrokeDataAndGo :: Box -> Writer ( [Box], [Box] ) [ IntervalNewtonTree Box ] + evalStrokeDataAndGo box@( t, i, s ) = go ( box, ( eqs t `Seq.index` i ) s ) + + go :: ( Box, StrokeDatum 3 𝕀 ) -- 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 ) ) - ) + go ( cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) ) + , i + , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) ) + ) + , sd@( StrokeDatum { ee = D22 _ ( T ee_t ) (T ee_s ) _ _ _ + , 𝛿E𝛿sdcdt = D12 _ ( T ( T f_t ) ) ( T ( T f_s ) ) }) + ) -- Box is small: stop processing it. | t_hi - t_lo < minWidth && s_hi - s_lo < minWidth = do let res = TooSmall cand tell ( [ cand ], [] ) return [ IntervalNewtonLeaf res ] - - | 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 - + | -- Check the envelope equation interval contains zero. + ee_potential_zero sd + -- Check the 𝛿E𝛿sdcdt interval contains zero. + , 𝛿E𝛿sdcdt_potential_zero sd , 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. - ee_potential_zero - -- Check the 𝛿E𝛿sdcdt interval contains zero. - , 𝛿E𝛿sdcdt_potential_zero - -> let -- Interval Newton method: take one Gauss–Seidel step - -- for the equation f'(x) v = - f(x_mid), - -- where f = 𝛿E/𝛿s * dc/dt - !( a, b ) = precondition precondMethod - ( midI f_t, midI f_s ) - ( 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) - -- ) + = let -- Interval Newton method: take one Gauss–Seidel step + -- for the equation f'(x) v = - f(x_mid), + -- where f = 𝛿E/𝛿s * dc/dt + !( a, b ) = precondition precondMethod + ( midI f_t, midI f_s ) + ( 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 - , coerce ( (-) @( 𝕀 Double ) ) s i_s_mid ) - in if all ( smaller . fst ) gsGuesses - then - -- If the Gauss–Seidel step was a contraction, then the box - -- contains a unique solution (by the Banach fixed point theorem). - -- - -- These boxes can thus be directly added to the solution set: - -- Newton's method is guaranteed to converge to the unique solution. - let !(done, todo) = bimap ( map ( mkGuess . fst ) ) ( map ( mkGuess . fst ) ) - $ partition snd gsGuesses - in do tell ([], done) - 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. - let (bisGuesses, whatBis) - | t_hi - t_lo > s_hi - s_lo - = ( [ ( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_mid ), i, s ) - , ( 𝕀 ( ℝ1 t_mid ) ( ℝ1 t_hi ), i, s ) ] - , ("t", t_mid) ) - | otherwise - = ( [ ( t, i, 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_mid ) ) - , ( t, i, 𝕀 ( ℝ1 s_mid ) ( ℝ1 s_hi ) ) ] - , ("s", s_mid) ) - in recur (IntervalNewtonStep (IntervalNewtonBisection whatBis)) bisGuesses + !gsGuesses = gaussSeidel a b + ( coerce ( (-) @( 𝕀 Double ) ) t i_t_mid + , coerce ( (-) @( 𝕀 Double ) ) s i_s_mid ) + in if any ( smaller . fst ) gsGuesses + then + -- If the Gauss–Seidel step was a contraction, then the box + -- contains a unique solution (by the Banach fixed point theorem). + -- + -- These boxes can thus be directly added to the solution set: + -- Newton's method is guaranteed to converge to the unique solution. + let !(done, todo) = bimap ( map ( mkGuess . fst ) ) ( map ( mkGuess . fst ) ) + $ partition snd gsGuesses + in do tell ([], done) + case todo of + [] -> return [ IntervalNewtonLeaf $ NoSolution "GaussSeidel" cand ] + _ -> recur evalStrokeDataAndGo ( IntervalNewtonStep ( IntervalNewtonContraction done ) ) + todo + else + -- Gauss–Seidel failed to shrink the boxes, so bisect instead. + -- We have to pick along which dimension to bisect: + -- - if bisecting along a particular dimension discards one of + -- the boxes, do that; + -- - otherwise, bisect along the dimension j that maximises + -- width(x_j) * || J_j ||. + let l_t = 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_mid ) + r_t = 𝕀 ( ℝ1 t_mid ) ( ℝ1 t_hi ) + d_s = 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_mid ) + u_s = 𝕀 ( ℝ1 s_mid ) ( ℝ1 s_hi ) + l = ( l_t, i, s ) + r = ( r_t, i, s ) + d = ( t, i, d_s ) + u = ( t, i, u_s ) + l_dat = ( eqs l_t `Seq.index` i ) s + r_dat = ( eqs r_t `Seq.index` i ) s + d_dat = ( eqs t `Seq.index` i ) d_s + u_dat = ( eqs t `Seq.index` i ) u_s + l_skip = + not ( ee_potential_zero l_dat ) + || not ( 𝛿E𝛿sdcdt_potential_zero l_dat ) + r_skip = + not ( ee_potential_zero r_dat ) + || not ( 𝛿E𝛿sdcdt_potential_zero r_dat ) + d_skip = + not ( ee_potential_zero d_dat ) + || not ( 𝛿E𝛿sdcdt_potential_zero d_dat ) + u_skip = + not ( ee_potential_zero u_dat ) + || not ( 𝛿E𝛿sdcdt_potential_zero u_dat ) + tJWidth = ( t_hi - t_lo ) * normVI3 ee_t f_t + sJWidth = ( s_hi - s_lo ) * normVI3 ee_s f_s + ( bisGuesses, whatBis ) + | l_skip && r_skip + = ( [], ( "lr", t_mid ) ) + | d_skip && u_skip + = ( [], ( "du", s_mid ) ) + | l_skip + = ( [ ( r, r_dat ) ], ( "r", t_mid ) ) + | r_skip + = ( [ ( l, l_dat ) ], ( "l", t_mid ) ) + | d_skip + = ( [ ( u, u_dat ) ], ( "u", s_mid ) ) + | u_skip + = ( [ ( d, d_dat ) ], ( "d", s_mid ) ) + | tJWidth >= sJWidth + -- t_hi - t_lo >= s_hi - s_lo + = ( [ ( l, l_dat ), ( r, r_dat ) ], ( "t", t_mid ) ) + | otherwise + = ( [ ( d, d_dat ), ( u, u_dat ) ], ( "s", s_mid ) ) + in recur go ( IntervalNewtonStep ( IntervalNewtonBisection whatBis ) . map (first fst) ) + bisGuesses -- Box doesn't contain a solution: discard it. | otherwise - -> return [ IntervalNewtonLeaf $ NoSolution (if ee_potential_zero then "dc/dt" else "ee") cand ] + = return [ IntervalNewtonLeaf $ NoSolution ( if ee_potential_zero sd then "dc/dt" else "ee" ) cand ] where midI :: 𝕀ℝ 2 -> 𝕀ℝ 2 midI ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = singleton $ ℝ2 ( 0.5 * ( x_lo + x_hi ) ) ( 0.5 * ( y_lo + y_hi ) ) + --width :: 𝕀ℝ 1 -> Double + --width ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) = hi - lo + --normI :: 𝕀ℝ 1 -> Double + --normI ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) = sqrt $ sup $ ( 𝕀 lo hi ) Ring.^ 2 + --normVI :: 𝕀ℝ 2 -> Double + --normVI ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = + -- sqrt $ sup $ ( 𝕀 x_lo x_hi ) Ring.^ 2 + ( 𝕀 y_lo y_hi ) Ring.^ 2 + normVI3 :: 𝕀ℝ 1 -> 𝕀ℝ 2 -> Double + normVI3 ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) + = sqrt $ max ( abs lo ) ( abs hi ) Ring.^ 2 + + max ( abs x_lo ) ( abs x_hi ) Ring.^ 2 + + max ( abs y_lo ) ( abs y_hi ) Ring.^ 2 t_mid = 0.5 * ( t_lo + t_hi ) s_mid = 0.5 * ( s_lo + s_hi ) i_t_mid = singleton ( ℝ1 t_mid ) @@ -1490,6 +1656,16 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox = !( 𝕀 y'_lo y'_hi ) = negate $ 𝕀 y_lo y_hi in 𝕀 ( ℝ2 x'_lo y'_lo ) ( ℝ2 x'_hi y'_hi ) +ee_potential_zero :: StrokeDatum 3 𝕀 -> Bool +ee_potential_zero dat = + inf ( _D22_v $ ee dat ) <= ℝ1 0 + && sup ( _D22_v $ ee dat ) >= ℝ1 0 +𝛿E𝛿sdcdt_potential_zero :: StrokeDatum 3 𝕀 -> Bool +𝛿E𝛿sdcdt_potential_zero dat = + cmpℝ2 (<=) ( inf $ unT $ _D12_v $ 𝛿E𝛿sdcdt dat ) ( ℝ2 0 0 ) + && cmpℝ2 (>=) ( sup $ unT $ _D12_v $ 𝛿E𝛿sdcdt dat ) ( ℝ2 0 0 ) + + zero, one :: Double zero = 1e-6 one = 1 - zero diff --git a/brush-strokes/src/Math/Bezier/Stroke/EnvelopeEquation.hs b/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs similarity index 89% rename from brush-strokes/src/Math/Bezier/Stroke/EnvelopeEquation.hs rename to brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs index 9c53198..d5a8adb 100644 --- a/brush-strokes/src/Math/Bezier/Stroke/EnvelopeEquation.hs +++ b/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs @@ -332,3 +332,31 @@ evaluateQuadratic bez t = maxs = fmap (Quadratic.bezier @( T Double ) sup_bez) $ inf t :| ( sup t : filter ( inside t ) ( Quadratic.extrema sup_bez ) ) in 𝕀 ( minimum mins ) ( maximum maxs ) + +{- + +evaluateCubic :: Cubic.Bezier ( 𝕀 Double ) -> 𝕀 Double -> 𝕀 Double +evaluateCubic bez t = + -- assert (inf t >= 0 && sup t <= 1) "evaluateCubic: t ⊊ [0,1]" $ -- Requires t βŠ† [0,1] + let inf_bez = Cubic.restrict @( T Double ) ( fmap inf bez ) ( inf t, sup t ) + sup_bez = Cubic.restrict @( T Double ) ( fmap sup bez ) ( inf t, sup t ) + mins = fmap (Cubic.bezier @( T Double ) inf_bez) + $ 0 :| ( 1 : Cubic.extrema inf_bez ) + maxs = fmap (Cubic.bezier @( T Double ) sup_bez) + $ 0 :| ( 1 : Cubic.extrema sup_bez ) + in 𝕀 ( minimum mins ) ( maximum maxs ) + +-- | Evaluate a quadratic BΓ©zier curve, when both the coefficients and the +-- parameter are intervals. +evaluateQuadratic :: Quadratic.Bezier ( 𝕀 Double ) -> 𝕀 Double -> 𝕀 Double +evaluateQuadratic bez t = + -- assert (inf t >= 0 && sup t <= 1) "evaluateCubic: t ⊊ [0,1]" $ -- Requires t βŠ† [0,1] + let inf_bez = Quadratic.restrict @( T Double ) ( fmap inf bez ) ( inf t, sup t ) + sup_bez = Quadratic.restrict @( T Double ) ( fmap sup bez ) ( inf t, sup t ) + mins = fmap (Quadratic.bezier @( T Double ) inf_bez) + $ 0 :| ( 1 : Quadratic.extrema inf_bez ) + maxs = fmap (Quadratic.bezier @( T Double ) sup_bez) + $ 0 :| ( 1 : Quadratic.extrema sup_bez ) + in 𝕀 ( minimum mins ) ( maximum maxs ) + +-} \ No newline at end of file diff --git a/brush-strokes/src/Math/Differentiable.hs b/brush-strokes/src/lib/Math/Differentiable.hs similarity index 100% rename from brush-strokes/src/Math/Differentiable.hs rename to brush-strokes/src/lib/Math/Differentiable.hs diff --git a/brush-strokes/src/Math/Epsilon.hs b/brush-strokes/src/lib/Math/Epsilon.hs similarity index 100% rename from brush-strokes/src/Math/Epsilon.hs rename to brush-strokes/src/lib/Math/Epsilon.hs diff --git a/brush-strokes/src/Math/Interval.hs b/brush-strokes/src/lib/Math/Interval.hs similarity index 100% rename from brush-strokes/src/Math/Interval.hs rename to brush-strokes/src/lib/Math/Interval.hs diff --git a/brush-strokes/src/Math/Interval/FMA.hs b/brush-strokes/src/lib/Math/Interval/FMA.hs similarity index 100% rename from brush-strokes/src/Math/Interval/FMA.hs rename to brush-strokes/src/lib/Math/Interval/FMA.hs diff --git a/brush-strokes/src/Math/Interval/Internal.hs b/brush-strokes/src/lib/Math/Interval/Internal.hs similarity index 94% rename from brush-strokes/src/Math/Interval/Internal.hs rename to brush-strokes/src/lib/Math/Interval/Internal.hs index b170745..d58e8e6 100644 --- a/brush-strokes/src/Math/Interval/Internal.hs +++ b/brush-strokes/src/lib/Math/Interval/Internal.hs @@ -120,19 +120,7 @@ instance Prelude.Fractional ( 𝕀 Double ) where | otherwise = error "BAD interval recip; should use extendedRecip instead" -- #endif - 𝕀 x_lo x_hi / 𝕀 y_lo y_hi --- #ifdef ASSERTS - | y_lo == 0 - = 𝕀 ( fst $ divI x_lo y_hi ) ( 1 / 0 ) - | y_hi == 0 - = 𝕀 ( -1 / 0 ) ( snd $ divI x_hi y_lo ) - | y_lo > 0 || y_hi < 0 --- #endif - = 𝕀 ( fst $ divI x_lo y_hi ) ( snd $ divI x_hi y_lo ) --- #ifdef ASSERTS - | otherwise - = error "BAD interval division; should use extendedRecip instead" --- #endif + p / q = p * recip q instance Floating ( 𝕀 Double ) where sqrt = withHW Prelude.sqrt diff --git a/brush-strokes/src/Math/Linear.hs b/brush-strokes/src/lib/Math/Linear.hs similarity index 100% rename from brush-strokes/src/Math/Linear.hs rename to brush-strokes/src/lib/Math/Linear.hs diff --git a/brush-strokes/src/Math/Linear/Internal.hs b/brush-strokes/src/lib/Math/Linear/Internal.hs similarity index 100% rename from brush-strokes/src/Math/Linear/Internal.hs rename to brush-strokes/src/lib/Math/Linear/Internal.hs diff --git a/brush-strokes/src/Math/Linear/Solve.hs b/brush-strokes/src/lib/Math/Linear/Solve.hs similarity index 100% rename from brush-strokes/src/Math/Linear/Solve.hs rename to brush-strokes/src/lib/Math/Linear/Solve.hs diff --git a/brush-strokes/src/Math/Module.hs b/brush-strokes/src/lib/Math/Module.hs similarity index 100% rename from brush-strokes/src/Math/Module.hs rename to brush-strokes/src/lib/Math/Module.hs diff --git a/brush-strokes/src/Math/Module/Internal.hs b/brush-strokes/src/lib/Math/Module/Internal.hs similarity index 100% rename from brush-strokes/src/Math/Module/Internal.hs rename to brush-strokes/src/lib/Math/Module/Internal.hs diff --git a/brush-strokes/src/Math/Monomial.hs b/brush-strokes/src/lib/Math/Monomial.hs similarity index 100% rename from brush-strokes/src/Math/Monomial.hs rename to brush-strokes/src/lib/Math/Monomial.hs diff --git a/brush-strokes/src/Math/Orientation.hs b/brush-strokes/src/lib/Math/Orientation.hs similarity index 100% rename from brush-strokes/src/Math/Orientation.hs rename to brush-strokes/src/lib/Math/Orientation.hs diff --git a/brush-strokes/src/Math/Ring.hs b/brush-strokes/src/lib/Math/Ring.hs similarity index 100% rename from brush-strokes/src/Math/Ring.hs rename to brush-strokes/src/lib/Math/Ring.hs diff --git a/brush-strokes/src/Math/Roots.hs b/brush-strokes/src/lib/Math/Roots.hs similarity index 100% rename from brush-strokes/src/Math/Roots.hs rename to brush-strokes/src/lib/Math/Roots.hs diff --git a/brush-strokes/src/TH/Utils.hs b/brush-strokes/src/lib/TH/Utils.hs similarity index 100% rename from brush-strokes/src/TH/Utils.hs rename to brush-strokes/src/lib/TH/Utils.hs diff --git a/src/convert/Main.hs b/brush-strokes/src/metafont/Main.hs similarity index 100% rename from src/convert/Main.hs rename to brush-strokes/src/metafont/Main.hs diff --git a/src/convert/MetaBrush/MetaFont/Convert.hs b/brush-strokes/src/metafont/MetaBrush/MetaFont/Convert.hs similarity index 100% rename from src/convert/MetaBrush/MetaFont/Convert.hs rename to brush-strokes/src/metafont/MetaBrush/MetaFont/Convert.hs diff --git a/src/app/MetaBrush/Application.hs b/src/app/MetaBrush/Application.hs index 5f62bed..253626f 100644 --- a/src/app/MetaBrush/Application.hs +++ b/src/app/MetaBrush/Application.hs @@ -158,7 +158,7 @@ runApplication application = do { strokeName = "Stroke 1" , strokeVisible = True , strokeUnique = strokeUnique - , strokeBrush = Just Asset.Brushes.tearDrop + , strokeBrush = Just Asset.Brushes.ellipse --tearDrop , strokeSpline = -- Spline -- { splineStart = mkPoint ( ℝ2 -20 -20 ) 5 @@ -169,7 +169,7 @@ runApplication application = do Spline { splineStart = mkPoint ( ℝ2 0 0 ) 10 25 0 , splineCurves = OpenCurves $ Seq.fromList - [ LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 100 0 ) 15 40 pi ), curveData = invalidateCache undefined } + [ LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 100 0 ) 15 40 (0.1 * pi) ), curveData = invalidateCache undefined } --, LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 -10 10 ) 8 5 ( pi / 4 ) ), curveData = invalidateCache undefined } --, LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 -10 -20 ) 10 7 ( pi / 2 ) ), curveData = invalidateCache undefined } ] @@ -179,12 +179,12 @@ runApplication application = do ) ] where - --mkPoint :: ℝ 2 -> Double -> Double -> Double -> PointData ( Record Asset.Brushes.EllipseBrushFields ) - --mkPoint pt a b phi = PointData pt Normal ( MkR $ ℝ3 a b phi ) + mkPoint :: ℝ 2 -> Double -> Double -> Double -> PointData ( Record Asset.Brushes.EllipseBrushFields ) + mkPoint pt a b phi = PointData pt Normal ( MkR $ ℝ3 a b phi ) --mkPoint :: ℝ 2 -> Double -> PointData ( Record Asset.Brushes.CircleBrushFields ) --mkPoint pt r = PointData pt Normal ( MkR $ ℝ1 r ) - mkPoint :: ℝ 2 -> Double -> Double -> Double -> PointData ( Record Asset.Brushes.TearDropBrushFields ) - mkPoint pt w h phi = PointData pt Normal ( MkR $ ℝ3 w h phi ) + --mkPoint :: ℝ 2 -> Double -> Double -> Double -> PointData ( Record Asset.Brushes.TearDropBrushFields ) + --mkPoint pt w h phi = PointData pt Normal ( MkR $ ℝ3 w h phi ) recomputeStrokesTVar <- STM.newTVarIO @Bool False documentRenderTVar <- STM.newTVarIO @( ( Int32, Int32 ) -> Cairo.Render () ) ( const $ pure () ) diff --git a/src/metabrushes/MetaBrush/Records.hs b/src/metabrushes/MetaBrush/Records.hs index 0d86207..d154ac2 100644 --- a/src/metabrushes/MetaBrush/Records.hs +++ b/src/metabrushes/MetaBrush/Records.hs @@ -139,12 +139,8 @@ instance ( Torsor ( T ( 𝕀ℝ ( Length ks ) ) ) ( 𝕀ℝ ( Length ks ) ) T ( 𝕀 ( MkR c_lo ) ( MkR c_hi ) ) type instance RepDim ( Record ks ) = Length 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 #-} +deriving newtype instance Representable r ( ℝ ( Length ks ) ) + => Representable r ( Record ks ) type instance D k ( Record ks ) = D k ( ℝ ( Length ks ) ) deriving newtype instance HasChainRule Double 2 ( ℝ ( Length ks ) )