Clean up benchmarking component

This commit is contained in:
sheaf 2024-04-22 22:11:07 +02:00
parent 0c54de8b1c
commit ac9deb968a
8 changed files with 286 additions and 563 deletions

View file

@ -81,7 +81,6 @@ common common
StandaloneDeriving
StandaloneKindSignatures
TupleSections
TypeAbstractions
TypeApplications
TypeFamilyDependencies
TypeOperators

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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