From ac9deb968ac9daee34ce4c1b8b242da0f8607d33 Mon Sep 17 00:00:00 2001 From: sheaf Date: Mon, 22 Apr 2024 22:11:07 +0200 Subject: [PATCH] Clean up benchmarking component --- MetaBrush.cabal | 1 - brush-strokes/brush-strokes.cabal | 17 +- brush-strokes/src/cusps/bench/Bench/Cases.hs | 78 +++ brush-strokes/src/cusps/bench/Bench/Types.hs | 120 ++++ brush-strokes/src/cusps/bench/Main.hs | 622 ++---------------- .../src/lib/Math/Root/Isolation/Bisection.hs | 1 + .../src/lib/Math/Root/Isolation/Core.hs | 8 +- .../src/lib/Math/Root/Isolation/Narrowing.hs | 2 +- 8 files changed, 286 insertions(+), 563 deletions(-) create mode 100644 brush-strokes/src/cusps/bench/Bench/Cases.hs create mode 100644 brush-strokes/src/cusps/bench/Bench/Types.hs diff --git a/MetaBrush.cabal b/MetaBrush.cabal index adaed84..c73d9ef 100644 --- a/MetaBrush.cabal +++ b/MetaBrush.cabal @@ -81,7 +81,6 @@ common common StandaloneDeriving StandaloneKindSignatures TupleSections - TypeAbstractions TypeApplications TypeFamilyDependencies TypeOperators diff --git a/brush-strokes/brush-strokes.cabal b/brush-strokes/brush-strokes.cabal index 1baf32c..23cf924 100644 --- a/brush-strokes/brush-strokes.cabal +++ b/brush-strokes/brush-strokes.cabal @@ -19,6 +19,8 @@ common common build-depends: base >= 4.17 && < 4.20 + , code-page + ^>= 0.2.1 , containers >= 0.6.0.1 && < 0.8 , tree-view @@ -52,17 +54,20 @@ common common StandaloneDeriving StandaloneKindSignatures TupleSections - TypeAbstractions TypeApplications TypeFamilyDependencies TypeOperators UnboxedTuples ViewPatterns + if impl(ghc >= 9.8) + default-extensions: + TypeAbstractions + ghc-options: -O2 -fexpose-all-unfoldings - -- -funfolding-use-threshold=1000 + -- -funfolding-use-threshold=1000 -fspecialise-aggressively -flate-dmd-anal -fmax-worker-args=200 @@ -127,9 +132,9 @@ library , Math.Roots , Math.Root.Isolation , Math.Root.Isolation.Bisection + , Math.Root.Isolation.Core , Math.Root.Isolation.GaussSeidel , Math.Root.Isolation.Narrowing - , Math.Root.Isolation.Core , Debug.Utils other-modules: @@ -148,8 +153,6 @@ library >= 2.18 && < 2.22 , bifunctors >= 5.5.4 && < 5.7 - , code-page - ^>= 0.2.1 , deepseq >= 1.4.4.0 && < 1.6 , directory @@ -257,6 +260,10 @@ benchmark cusps main-is: Main.hs + other-modules: + Bench.Types + Bench.Cases + build-depends: brush-strokes , deepseq diff --git a/brush-strokes/src/cusps/bench/Bench/Cases.hs b/brush-strokes/src/cusps/bench/Bench/Cases.hs new file mode 100644 index 0000000..3d206f5 --- /dev/null +++ b/brush-strokes/src/cusps/bench/Bench/Cases.hs @@ -0,0 +1,78 @@ +module Bench.Cases where + +-- brush-strokes +import Calligraphy.Brushes +import Math.Bezier.Spline +import Math.Bezier.Stroke +import Math.Interval +import Math.Linear +import Math.Module +import Math.Root.Isolation + +-- brush-strokes bench:cusps +import Bench.Types + +-------------------------------------------------------------------------------- + +ellipseTestCase :: RootIsolationOptions N 3 -> String -> ( Double, Double ) -> Double -> [ ( Int, [ Box 2 ] ) ] -> TestCase +ellipseTestCase opts descr k0k1 rot startBoxes = + TestCase + { testDescription = descr + , testBrushStroke = ellipseBrushStroke k0k1 rot + , testCuspOptions = opts + , 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 :: Double -> Double -> Double -> Double + 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 + { testDescription = "" + , testBrushStroke = trickyCusp2BrushStroke + , testCuspOptions = defaultRootIsolationOptions + , 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 + + +defaultStartBoxes :: [ Int ] -> [ ( Int, [ Box 2 ] ) ] +defaultStartBoxes is = + [ ( i, [ 𝕀 ( ℝ2 zero zero ) ( ℝ2 one one ) ] ) | i <- is ] + +zero, one :: Double +zero = 1e-6 +one = 1 - zero +{-# INLINE zero #-} +{-# INLINE one #-} diff --git a/brush-strokes/src/cusps/bench/Bench/Types.hs b/brush-strokes/src/cusps/bench/Bench/Types.hs new file mode 100644 index 0000000..5c559d7 --- /dev/null +++ b/brush-strokes/src/cusps/bench/Bench/Types.hs @@ -0,0 +1,120 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE UndecidableInstances #-} + +module Bench.Types + ( BrushStroke(..) + , TestCase(..) + , Point(..) + , Params(..) + , brushStrokeFunctions + ) where + +-- base +import Data.Coerce + ( coerce ) +import GHC.Generics + ( Generic ) + +-- containers +import Data.Sequence + ( Seq ) + +-- brush-strokes +import Calligraphy.Brushes +import Math.Algebra.Dual +import Math.Bezier.Spline +import Math.Bezier.Stroke +import Math.Bezier.Stroke.EnvelopeEquation +import Math.Interval +import Math.Linear +import Math.Module +import Math.Ring + ( Transcendental ) +import Math.Root.Isolation + +-------------------------------------------------------------------------------- + +data BrushStroke = + forall nbParams. ParamsCt nbParams => + BrushStroke + { brush :: !( Brush ( ℝ nbParams ) ) + , stroke :: !( Point nbParams, Curve Open () ( Point nbParams ) ) + } + +data TestCase = + TestCase + { testDescription :: String + , testBrushStroke :: BrushStroke + , testCuspOptions :: RootIsolationOptions N 3 + , testStartBoxes :: [ ( Int, [ Box 2 ] ) ] + } + +-------------------------------------------------------------------------------- + +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 ) + +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 brushShape brushShapeI mbRot ) 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 + brushShape + mbRot + , brushStrokeData @3 @( ℝ nbParams ) coerce coerce + pathI usedParamsI + brushShapeI + ( fmap nonDecreasing mbRot ) + ) +{-# INLINEABLE getStrokeFunctions #-} + +brushStrokeFunctions + :: BrushStroke + -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) + , 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) +brushStrokeFunctions ( BrushStroke { stroke = ( sp0, crv ), brush } ) = + getStrokeFunctions brush sp0 crv diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs index f3f4f60..76fb50e 100644 --- a/brush-strokes/src/cusps/bench/Main.hs +++ b/brush-strokes/src/cusps/bench/Main.hs @@ -1,123 +1,69 @@ -{-# LANGUAGE PolyKinds #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE UndecidableInstances #-} - -module Main - ( main - - -- Testing - , TestCase(..) - , testCases - , BrushStroke(..) - , brushStrokeFunctions - , eval - , mkVal, mkBox - , potentialCusp - , dEdsdcdt - ) - where +module Main where -- base -import Control.Concurrent.MVar - ( newMVar ) -import Data.Coerce - ( coerce ) +import Control.Monad + ( when ) import Data.Foldable ( for_ ) -import Data.List - ( intercalate, sortOn ) +import qualified Data.List.NonEmpty as NE + ( NonEmpty(..) + , fromList, head, length, sort + ) +import Data.Semigroup + ( Arg(..) ) import Data.Traversable ( for ) import GHC.Clock ( getMonotonicTime ) -import GHC.Exts - ( Proxy#, proxy# ) -import GHC.Generics - ( Generic ) -import GHC.TypeNats - ( Nat, type (-) ) -import Numeric - ( showFFloat ) + +-- code-page +import System.IO.CodePage + ( withCP65001 ) -- containers import qualified Data.IntMap.Strict as IntMap ( fromList, toList ) -import Data.Sequence - ( Seq ) -import qualified Data.Sequence as Seq - ( index ) -import Data.Tree - ( foldTree ) -- deepseq import Control.DeepSeq ( rnf ) --- gauge ---import qualified Gauge - --- 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 ) import Math.Root.Isolation +-- brush-strokes bench:cusps +import Bench.Cases +import Bench.Types + -------------------------------------------------------------------------------- -{- main :: IO () -main = Gauge.defaultMainWith benchConfig - [ Gauge.bgroup "Cusp finding" - [ Gauge.bench testName $ Gauge.nf benchTest benchCase - | benchCase@( TestCase { testName } ) <- benchCases - ] - ] - -benchConfig :: Gauge.Config -benchConfig = Gauge.defaultConfig - { Gauge.timeLimit = Just 0 - , Gauge.minSamples = Nothing - , Gauge.minDuration = 0 - } - -benchTest :: TestCase -> () -benchTest ( TestCase { testBrushStroke, testCuspOptions, testStartBoxes } ) = - let ( _, testStrokeFnI ) = brushStrokeFunctions testBrushStroke - ( _, ( dunno, sols ) ) = computeCusps testCuspOptions testStrokeFnI testStartBoxes - in dunno `seq` sols `seq` () --} - -main :: IO () -main = do +main = withCP65001 $ do putStrLn "Running cusp-finding benchmarks." - testsWithTime <- - for testCases $ \ tst -> do { dt <- benchTestCase tst; return ( tst, dt ) } - let (bestTest, bestTime) = head $ sortOn snd testsWithTime - putStrLn $ unlines - [ replicate 40 '=' - , "Best time: " ++ show bestTime ++ "s" - , "Best parameters: " ++ show (testName bestTest) - ] + for_ benchGroups $ \ ( groupName, benchGroupCases ) -> do + putStrLn $ unlines + [ replicate 40 '=' + , "Benchmark group '" ++ groupName ++ "':" ] + testsWithTime <- for benchGroupCases $ \ tst -> do + dt <- benchTestCase groupName tst + return $ Arg dt tst + let Arg bestTime bestTest = NE.head $ NE.sort testsWithTime + when ( NE.length benchGroupCases >= 1 ) $ + putStrLn $ unlines $ + [ "Best time in '" ++ groupName ++ "' group: " ++ show bestTime ++ "s" ] + ++ [ " (" ++ descr ++ ")" + | let descr = testDescription bestTest + , not ( null descr ) + ] -benchTestCase :: TestCase -> IO Double -benchTestCase ( TestCase { testName, testBrushStroke, testCuspOptions, testStartBoxes } ) = do - putStrLn $ unlines - [ replicate 40 '=' - , "Test case " ++ testName ] +-- TODO: run the benchmark multiple times to reduce noise. +benchTestCase :: String -> TestCase -> IO Double +benchTestCase testName ( TestCase { testDescription, testBrushStroke, testCuspOptions, testStartBoxes } ) = do + putStr $ unlines + [ " " ++ replicate 38 '-' + , " -- Test case: " ++ testName ++ ( if null testDescription then "" else " (" ++ testDescription ++ ")" ) + , " --" ] before <- getMonotonicTime let ( _, testStrokeFnI ) = brushStrokeFunctions testBrushStroke ( dunno, sols ) = @@ -133,90 +79,40 @@ benchTestCase ( TestCase { testName, testBrushStroke, testCuspOptions, testStart rnf dunno `seq` rnf sols `seq` return () after <- getMonotonicTime let dt = after - before - putStrLn $ unlines - [ " - #sols: " ++ show sols --( length sols ) - , " - #dunno: " ++ show ( length dunno ) - , " - Time elapsed: " ++ show dt ++ "s" - , "" ] + putStrLn $ unlines $ + [ " -- sols:" ++ ( if null sols then " ∅" else "" ) ] + ++ + ( map ( \ sol -> " -- • " ++ show sol ) sols ) + ++ + [ " -- #dunno: " ++ show ( length dunno ) + , " -- Time elapsed: " ++ show dt ++ "s" ] return dt -{- -main :: IO () -main = - for_ testCases $ \ testCase@( TestCase { testName, testBrushStroke, testCuspOptions, testStartBoxes } ) -> do - let ( _, testStrokeFnI ) = brushStrokeFunctions testBrushStroke - ( newtTrees, ( dunno, sols ) ) = computeCusps testCuspOptions 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 opts ("minWidth=" ++ show minWidth ++ ",ε=" ++ show narrowAbs) ( 0, 1 ) pi $ defaultStartBoxes [ 2 ] - | minWidth <- [ 1e-5 ] - , narrowAbs <- [ 5e-2 ] - , let opts = - RootIsolationOptions - { rootIsolationAlgorithms = defaultRootIsolationAlgorithms minWidth narrowAbs - } +benchGroups :: [ ( String, NE.NonEmpty TestCase ) ] +benchGroups = + [ ( "ellipse" + , NE.fromList + [ ellipseTestCase opts ("ε=" ++ show narrowAbs) + ( 0, 1 ) pi + ( defaultStartBoxes [ 0 .. 3 ] ) + | narrowAbs <- [ 5e-2, 1e-6 ] + , let opts = + RootIsolationOptions + { rootIsolationAlgorithms = + defaultRootIsolationAlgorithms minWidth narrowAbs + } + ] + ) + , ( "trickyCusp2", NE.fromList [ trickyCusp2TestCase ] ) ] + where + minWidth = 1e-5 -------------------------------------------------------------------------------- -data BrushStroke = - forall nbParams. ParamsCt nbParams => - BrushStroke - { brush :: !( Brush ( ℝ nbParams ) ) - , stroke :: !( Point nbParams, Curve Open () ( Point nbParams ) ) - } +{- -- Old testing code. - -data TestCase = - TestCase - { testName :: String - , testBrushStroke :: BrushStroke - , testCuspOptions :: RootIsolationOptions N 3 - , testStartBoxes :: [ ( Int, [ Box 2 ] ) ] - } - -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. +getR1 (ℝ1 u) = u eval :: ( I i ( ℝ 1 ) -> Seq ( I i ( ℝ 1 ) -> StrokeDatum k i ) ) @@ -234,12 +130,6 @@ mkBox :: ( Double, Double ) -> ( Double, Double ) -> Box 2 mkBox ( t_min, t_max ) ( s_min, s_max ) = ( 𝕀 ( ℝ2 t_min s_min ) ( ℝ2 t_max s_max ) ) -zero, one :: Double -zero = 1e-6 -one = 1 - zero -{-# INLINE zero #-} -{-# INLINE one #-} - potentialCusp :: StrokeDatum 3 𝕀 -> Bool potentialCusp ( StrokeDatum @@ -253,118 +143,6 @@ potentialCusp 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 ) ) = isolateRootsIn 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 $ isolateRootsIn 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 @@ -409,276 +187,10 @@ logLines = 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 $ isolateRootsIn 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 :: RootIsolationOptions N 3 -> String -> ( Double, Double ) -> Double -> [ ( Int, [ Box 2 ] ) ] -> TestCase -ellipseTestCase opts str k0k1 rot startBoxes = - TestCase - { testName = "ellipse (" ++ str ++ ")" - , testBrushStroke = ellipseBrushStroke k0k1 rot - , testCuspOptions = opts - , 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 - , testCuspOptions = defaultRootIsolationOptions - , 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 ) - -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 brushShape brushShapeI mbRot ) 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 - brushShape - mbRot - , brushStrokeData @3 @( ℝ nbParams ) coerce coerce - pathI usedParamsI - brushShapeI - ( fmap nonDecreasing mbRot ) - ) -{-# INLINEABLE getStrokeFunctions #-} - -defaultStartBoxes :: [ Int ] -> [ ( Int, [ Box 2 ] ) ] -defaultStartBoxes is = - [ ( i, [ 𝕀 ( ℝ2 zero zero ) ( ℝ2 one one ) ] ) | i <- is ] --- [ ( i, [ 𝕀 ( ℝ2 t s ) ( ℝ2 ( t + 0.099999 ) ( s + 0.099999 ) ) --- | t <- [ 0.00001, 0.1 .. 0.9 ] --- , s <- [ 0.00001, 0.1 .. 0.9 ] --- ] ) --- | i <- is --- ] - -getR1 (ℝ1 u) = u - -{- - -(f, fI) = brushStrokeFunctions $ ellipseBrushStroke (0,1) pi -nbPotentialSols box = let ( _newtTrees, ( dunno, sols ) ) = isolateRootsIn NoPreconditioning 1e-7 fI box in length dunno + length sols -showTrees box = putStrLn $ unlines $ map Data.Tree.View.showTree $ map ( uncurry showIntervalNewtonTree ) $ fst $ isolateRootsIn NoPreconditioning 1e-7 fI box - -sol_t = 0.5486100729150693796677845183880669025324233347060776339185 :: Double -sol_s = 0.5480950141859386853197594577293968665598143630958601978245 :: Double -containsSol (t, _i, s) = getR1 (inf t) <= sol_t && getR1 (sup t) >= sol_t && getR1 (inf s) <= sol_s && getR1 (sup s) >= sol_s - -(t, i, s) = mkBox (0.54, 0.55) 2 (0.5480, 0.5481) -containsSol (t, i, s) -nbPotentialSols (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) -D0 ( T f_mid ) = dEdsdcdt $ eval f $ mkVal t_mid 2 s_mid -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 f_mid -[((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.5451193766263323, ℝ1 0.545225929860598] -s2 -> [ℝ1 0.548, ℝ1 0.5481] -containsSol (t2, i, s2) -> 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 - -( b1 - a12 * x2 ) `extendedDivide` a11 --} - -{- - -bound consistency / box consistency - -ghci> (f, fI) = brushStrokeFunctions $ ellipseBrushStroke (0,1) pi -ghci> ee $ eval fI $ mkBox (0, succFP 0) 2 (0.5,0.75) -D22 {_D22_v = [ℝ1 -4379.198802308688, ℝ1 -1791.4548164846506], _D22_dx = T[ℝ1 -3417.0023260831304, ℝ1 6215.362085329192], _D22_dy = T[ℝ1 3524.130928930924, ℝ1 8849.881134814816], _D22_dxdx = T[ℝ1 -13397.019488083239, ℝ1 61439.0409103 - - - -} \ No newline at end of file diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Bisection.hs b/brush-strokes/src/lib/Math/Root/Isolation/Bisection.hs index c133bd8..f51927d 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation/Bisection.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation/Bisection.hs @@ -64,6 +64,7 @@ instance RootIsolationAlgorithm Bisection where ( fallbackBisectionCoord thisRoundHist prevRoundsHist eqs ) box return ( whatBis, boxes ) + {-# INLINEABLE rootIsolationAlgorithm #-} -- | Options for the bisection method. type BisectionOptions :: Nat -> Nat -> Type diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Core.hs b/brush-strokes/src/lib/Math/Root/Isolation/Core.hs index 498a455..2f28a2d 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation/Core.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation/Core.hs @@ -234,7 +234,11 @@ showArea :: Double -> String showArea area = "(area " ++ showFFloat (Just 6) area "" ++ ")" -------------------------------------------------------------------------------- --- Utilities. +-- Utilities for feeding one computation into the next. + +-- NB: I tried adding RULES on these functions in order to perform +-- "loop unrolling"-like optimisations, but this does not seem to +-- improve the performance. -- | Run an effectful computation several times in sequence, piping its output -- into the next input, once for each coordinate dimension. @@ -263,3 +267,5 @@ pipeFunctionsWhileTrue fns = go fns go ( f : fs ) x = do xs <- f x concat <$> traverse ( go fs ) xs + +-------------------------------------------------------------------------------- diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Narrowing.hs b/brush-strokes/src/lib/Math/Root/Isolation/Narrowing.hs index d292718..9ff0dd5 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation/Narrowing.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation/Narrowing.hs @@ -227,7 +227,7 @@ narrowingMethods ε_eq ε_bis ( AdaptiveShaving opts ) -- TODO: haven't implemented right shaving yet narrowingMethods _ε_eq ε_bis TwoSidedShaving = [ sbc ε_bis ] -{-# INLINE narrowingMethods #-} +{-# INLINEABLE narrowingMethods #-} allNarrowingOperators :: forall n d