brush-strokes-0.1.0.0: src/arc-length/bench/BenchArcLength.hs
{-# OPTIONS_GHC -Wno-x-partial #-}
module Main where
-- base
import Control.Exception
( evaluate )
import Data.Traversable
( for )
import GHC.Clock
( getMonotonicTime )
import GHC.Exts
( SPEC(..) )
import Numeric
( showEFloat, showFFloat )
-- code-page
import System.IO.CodePage
( withCP65001 )
-- base
import Data.List
( intercalate )
-- brush-strokes
import Math.Bezier.ArcLength
( ArcLengthOptions(..), defaultArcLengthOptions
, Integrator(..)
, Regulariser(..), defaultRegulariser
, ArcLengthParametrisation(..), curveArcLengthParametrisation
)
import qualified Math.Bezier.Cubic as Cubic
import qualified Math.Bezier.Quadratic as Quadratic
import Math.Bezier.Spline
( Curve(..), NextPoint(..) )
import Math.Linear
( ℝ(..), T )
import Math.Module
( distance )
--------------------------------------------------------------------------------
-- Polyline arc-length
-- | Approximate the arc length of a parametric curve as a polyline with @n@
-- chords.
polylineArcLength :: Int -> (Double -> ℝ 2) -> Double
polylineArcLength n f =
sum ( zipWith ( distance @( T (ℝ 2) ) ) pts ( tail pts ) )
where
pts = [ f ( fromIntegral i / fromIntegral n ) | i <- [0 .. n] ]
--------------------------------------------------------------------------------
-- Arc-length test case
-- | A benchmark case for arc-length estimation
data ArcLengthCase = ArcLengthCase
{ caseName :: !String
, caseCurve :: !( Double -> ℝ 2 )
, caseIntegrate :: !( ArcLengthOptions -> Double )
}
highAccuracyRef :: ( Double -> ℝ 2 ) -> Double
highAccuracyRef = romberg 50000 1
romberg :: Int -> Int -> (Double -> ℝ 2) -> Double
romberg nBase depth f =
let
baseLengths = [ polylineArcLength (nBase * (2 ^ k)) f | k <- [0 .. depth] ]
-- Tableau generation on a tiny list of Doubles (no performance hit here)
rombergTableau :: [Double] -> [[Double]]
rombergTableau col0 = col0 : go 1 col0
where
go :: Int -> [Double] -> [[Double]]
go k col =
let factor = 4 ** fromIntegral k
nextCol = zipWith (\prev curr -> (factor * curr - prev) / (factor - 1)) col (tail col)
in nextCol : go (k + 1) nextCol
-- Grab the fully extrapolated value (the single element in the final column)
in head (rombergTableau baseLengths !! depth)
-- | Construct a quadratic Bézier test case.
quadCase :: String -> ℝ 2 -> ℝ 2 -> ℝ 2 -> ArcLengthCase
quadCase name p0 p1 p2 =
ArcLengthCase
{ caseName = name
, caseCurve = Quadratic.bezier @( T (ℝ 2) )
$ Quadratic.Bezier { p0, p1, p2 }
, caseIntegrate = \ opts ->
totalArcLength $
curveArcLengthParametrisation @( T (ℝ 2) ) opts p0
Bezier2To
{ controlPoint = p1
, curveEnd = NextPoint p2
, curveData = ()
}
}
-- | Construct a cubic Bézier test case.
cubicCase :: String -> ℝ 2 -> ℝ 2 -> ℝ 2 -> ℝ 2 -> ArcLengthCase
cubicCase name p0 p1 p2 p3 =
ArcLengthCase
{ caseName = name
, caseCurve = Cubic.bezier @( T (ℝ 2) )
$ Cubic.Bezier { p0, p1, p2, p3 }
, caseIntegrate = \ opts ->
totalArcLength $
curveArcLengthParametrisation @( T (ℝ 2) ) opts p0
Bezier3To
{ controlPoint1 = p1
, controlPoint2 = p2
, curveEnd = NextPoint p3
, curveData = ()
}
}
--------------------------------------------------------------------------------
-- Integration methods
-- | The integration methods to compare.
data Method
= Polyline !Int -- ^ polyline with @n@ chords
| Romberg !Int !Int -- ^ Romberg method
| CC !Bool !( Maybe Regulariser ) !Int -- ^ Clenshaw–Curtis: split near crits, regulariser, degree @d@
| TS !Bool !Int !Double -- ^ tanh–sinh: split near crits, @N@ nodes per side, step size @h@
| GL !Bool !( Maybe Regulariser ) !Int -- ^ Gauss–Legendre: split near crits, regulariser, @n@ nodes
| GV !Double -- ^ Gravesen adaptive polygon/chord, stopping tolerance
methodLabel :: Method -> String
methodLabel ( Polyline n ) =
"Polyline, n=" ++ show n
methodLabel ( Romberg n d ) =
"Romberg, n=" ++ show n ++ " d=" ++ show d
methodLabel ( CC cs reg d ) =
"Clenshaw-Curtis, d=" ++ show d ++ pprSplitting cs ++ pprRegulariser reg
methodLabel ( TS cs n h ) =
"TanhSinh, N=" ++ show n ++ " h=" ++ showFFloat ( Just 2 ) h "" ++ pprSplitting cs
methodLabel ( GL cs reg n ) =
"Gauss-Legendre, n=" ++ show n ++ pprSplitting cs ++ pprRegulariser reg
methodLabel ( GV eps ) =
"Gravesen, tol=" ++ showEFloat ( Just 0 ) eps ""
pprSplitting :: Bool -> String
pprSplitting False = ""
pprSplitting True = " [S]"
pprRegulariser :: Maybe Regulariser -> String
pprRegulariser Nothing = ""
pprRegulariser ( Just _ ) = " [R:√]"
-- | Compute the arc length for a given test case using the specified method.
computeArcLength :: ArcLengthCase -> Method -> Double
computeArcLength ( ArcLengthCase { caseCurve = crv, caseIntegrate } ) = \case
Polyline n ->
polylineArcLength n crv
Romberg n d ->
romberg n d crv
CC cs reg d ->
caseIntegrate
( defaultArcLengthOptions
{ arcLengthIntegrator = ClenshawCurtis
{ clenshawCurtisDegree = d
, regulariser = reg
}
, splitAtCriticalPoints = cs
} )
TS cs n h ->
caseIntegrate
( defaultArcLengthOptions
{ arcLengthIntegrator = TanhSinh { tanhSinhNodes = n, tanhSinhStepSize = h }
, splitAtCriticalPoints = cs
} )
GL cs reg n ->
caseIntegrate
( defaultArcLengthOptions
{ arcLengthIntegrator = GaussLegendre
{ gaussLegendreDegree = n
, regulariser = reg
}
, splitAtCriticalPoints = cs
} )
GV eps ->
caseIntegrate
( defaultArcLengthOptions
{ arcLengthIntegrator = Gravesen { gravesenTol = eps }
} )
--------------------------------------------------------------------------------
-- Timing
-- | Core timing loop, following the tasty-bench technique.
--
-- @f@ and @x@ are explicit parameters so that GHC's full-laziness transformation
-- cannot hoist @f x@ outside the loop.
-- The @SPEC@ argument defeats the static-argument transformation (SAT) that
-- would otherwise share the computation across recursive calls.
{-# NOINLINE benchLoop #-}
benchLoop :: SPEC -> (a -> Double) -> a -> Int -> IO Double
benchLoop !_ f x n
| n <= 1 = evaluate (f x)
| otherwise = evaluate (f x) >> benchLoop SPEC f x (n - 1)
-- | Time @f x@ for @n@ repetitions, preceded by one warm-up call.
-- Returns (result of last run, seconds per run).
timeN :: Int -> (a -> Double) -> a -> IO ( Double, Double )
timeN n f x = do
_ <- benchLoop SPEC f x 1 -- warm up
before <- getMonotonicTime
result <- benchLoop SPEC f x n
after <- getMonotonicTime
return ( result, ( after - before ) / fromIntegral n )
-- | Number of timed repetitions per (case, method) pair.
numReps :: Int
numReps = 500
-- | Time a single (case, method) pair.
timeMethod :: ArcLengthCase -> Method -> IO ( Double, Double )
timeMethod tc method = timeN numReps (computeArcLength tc) method
--------------------------------------------------------------------------------
-- Benchmark groups
-- | A named collection of test cases.
data BenchGroup = BenchGroup
{ groupName :: !String
, groupCases :: ![ArcLengthCase]
}
-- | The methods benchmarked for every (group, case) pair.
allMethods :: [Method]
allMethods =
concat
[ [ Polyline n | n <- [ 40, 200, 1000 ] ]
, [ GV e | e <- [ 1e-1, 1e-3, 1e-5 ] ]
, [ Romberg n d | (n,d) <- [ (40, 1), (20, 2), (20, 3), (100, 2) ] ]
, [ GL cs reg n | n <- [ 8, 20 ], cs <- [ False, True ], reg <- [ Nothing, Just defaultRegulariser ], maybe True ( const cs ) reg ]
, [ CC cs reg n | n <- [ 8, 20 ], cs <- [ False, True ], reg <- [ Nothing, Just defaultRegulariser ], maybe True ( const cs ) reg ]
, [ TS cs n h | (n,h) <- [ (8, 0.4), (32, 0.1) ], cs <- [ False, True ] ]
]
benchGroups :: [BenchGroup]
benchGroups =
[ BenchGroup
{ groupName = "Quadratic Béziers"
, groupCases =
[ quadCase "gentle arc" ( ℝ2 0 0 ) ( ℝ2 1 0 ) ( ℝ2 1 1 )
, quadCase "sharp arc" ( ℝ2 0 0 ) ( ℝ2 0.5 3.0 ) ( ℝ2 0.5 0 )
, quadCase "line" ( ℝ2 0 0 ) ( ℝ2 0.5 0.5 ) ( ℝ2 1 1 )
, quadCase "nearly flat" ( ℝ2 0 0 ) ( ℝ2 4.8 0 ) ( ℝ2 4 -0.01 )
, quadCase "(P1 = P0)" ( ℝ2 0 0 ) ( ℝ2 0 0 ) ( ℝ2 1 1 )
, quadCase "degenerate" ( ℝ2 0 0 ) ( ℝ2 1 1 ) ( ℝ2 -1 -1 )
]
}
, BenchGroup
{ groupName = "Cubic Béziers"
, groupCases =
[ cubicCase "simple arc" ( ℝ2 0 0 ) ( ℝ2 0 1 ) ( ℝ2 1 1 ) ( ℝ2 1 0 )
, cubicCase "circular arc" ( ℝ2 1 0 ) ( ℝ2 1 κ ) ( ℝ2 κ 1 ) ( ℝ2 0 1 )
, cubicCase "loose S" ( ℝ2 0 0 ) ( ℝ2 0.33 0 ) ( ℝ2 0.67 1 ) ( ℝ2 1 1 )
, cubicCase "tight S" ( ℝ2 0 0 ) ( ℝ2 2 0 ) ( ℝ2 -1 1 ) ( ℝ2 1 1 )
, cubicCase "inflection S" ( ℝ2 0 0 ) ( ℝ2 1 0 ) ( ℝ2 0 1 ) ( ℝ2 1 1 )
, cubicCase "loop" ( ℝ2 0 0 ) ( ℝ2 1 1 ) ( ℝ2 -1 1 ) ( ℝ2 0 0 )
, cubicCase "gamma shape" ( ℝ2 1 2 ) ( ℝ2 0.16 -1 ) ( ℝ2 3.5 1.4 ) ( ℝ2 0.25 1.46 )
, cubicCase "simple cusp" ( ℝ2 0 0 ) ( ℝ2 1 1 ) ( ℝ2 0 1 ) ( ℝ2 1 0 )
, cubicCase "tricky cusp" ( ℝ2 0 0 ) ( ℝ2 1 0 ) ( ℝ2 -0.82 -0.18 ) ( ℝ2 1 1 )
, cubicCase "tricky cusp 2" ( ℝ2 0 0 ) ( ℝ2 1 0 ) ( ℝ2 1.22 2.52 ) ( ℝ2 1 1 )
, cubicCase "near cusp" ( ℝ2 0 0 ) ( ℝ2 1 0 ) ( ℝ2 -0.7 -0.6 ) ( ℝ2 1 1 )
, cubicCase "near cusp 2" ( ℝ2 0 0 ) ( ℝ2 1 0 ) ( ℝ2 -0.25 -1.88 ) ( ℝ2 1 1 )
, cubicCase "line" ( ℝ2 0 0 ) ( ℝ2 0.3 0.3 ) ( ℝ2 0.5 0.5 ) ( ℝ2 1 1 )
, cubicCase "degenerate (collinear)" ( ℝ2 0 0 ) ( ℝ2 1 0 ) ( ℝ2 -1 0 ) ( ℝ2 0 0 )
, cubicCase "very degenerate" ( ℝ2 0 0 ) ( ℝ2 1 0 ) ( ℝ2 1 0 ) ( ℝ2 0 0 )
, cubicCase "(P0 = P1)" ( ℝ2 0 0 ) ( ℝ2 0 0 ) ( ℝ2 1 0 ) ( ℝ2 1 1 )
, cubicCase "(P1 = P2)" ( ℝ2 0 0 ) ( ℝ2 1 0 ) ( ℝ2 1 0 ) ( ℝ2 1 1 )
]
}
]
κ :: Double
κ = 0.5519150244935105707435627227925
--------------------------------------------------------------------------------
-- Benchmark result record
data BenchResult = BenchResult
{ resGroup :: !String
, resCase :: !String
, resMethod :: !String
, resTimeUs :: !Double
, resRelErr :: !Double
}
--------------------------------------------------------------------------------
-- Running and reporting
-- | Print the header row and separator for the result table.
printTableHeader :: IO ()
printTableHeader = do
let row = " "
++ lAlign 42 "Method"
++ rAlign 11 "Time"
++ rAlign 13 "Rel. error"
putStrLn row
putStrLn $ " " ++ replicate ( 42 + 11 + 13 ) '─'
-- | Print one result row.
printResultRow :: String -> Double -> Double -> Double -> IO ()
printResultRow label usPerCall _arcLen relErr =
putStrLn $ " "
++ lAlign 42 label
++ rAlign 11 ( showFFloat ( Just 2 ) usPerCall "μs" )
++ rAlign 13 ( showRelErr relErr )
-- | Run and print all benchmarks for a group, returning all results.
runGroup :: BenchGroup -> IO [BenchResult]
runGroup BenchGroup { groupName, groupCases } = do
putStrLn ""
putStrLn $ replicate 80 '='
putStrLn $ " " ++ groupName
fmap concat $ for groupCases $ \tc -> do
-- Compute the reference value without benchmarking its time.
let refLen = highAccuracyRef ( caseCurve tc )
putStrLn ""
putStrLn $ " ── " ++ caseName tc
-- putStrLn $ " reference: " ++ showFFloat ( Just 15 ) refLen ""
putStrLn ""
printTableHeader
for allMethods $ \method -> do
( result, dt ) <- timeMethod tc method
let relErr = abs ( result - refLen ) / max refLen 1e-15
usPerCall = dt * 1e6
printResultRow ( methodLabel method ) usPerCall result relErr
return BenchResult
{ resGroup = groupName
, resCase = caseName tc
, resMethod = methodLabel method
, resTimeUs = usPerCall
, resRelErr = relErr
}
main :: IO ()
main = withCP65001 $ do
putStrLn "Arc-length quadrature benchmarks"
putStrLn $ unlines
[ ""
, "Legend:"
, " R: regularisation method (change of variables)"
, " S: splitting (subdividing at near-critical points)"
]
results <- fmap concat $ mapM runGroup benchGroups
putStrLn ""
putStrLn $ replicate 80 '='
let csvPath = "bench_results.csv"
writeCSV csvPath results
putStrLn $ "CSV results written to: " ++ csvPath
-- | Write benchmark results to a CSV file.
writeCSV :: FilePath -> [BenchResult] -> IO ()
writeCSV path results = writeFile path $ unlines ( header : map row results )
where
header = "group,case,method,time_us,rel_error"
row BenchResult { resGroup, resCase, resMethod, resTimeUs, resRelErr } =
intercalate ","
[ csvQuote resGroup
, csvQuote resCase
, csvQuote resMethod
, show resTimeUs
, show resRelErr
]
csvQuote s
| any ( `elem` ",\"\n" ) s = "\"" ++ concatMap esc s ++ "\""
| otherwise = s
where esc '"' = "\"\""
esc c = [c]
--------------------------------------------------------------------------------
-- Formatting helpers
-- | Left-align a string in a field of width @n@ (pad on the right).
lAlign :: Int -> String -> String
lAlign n s = take n ( s ++ repeat ' ' )
-- | Right-align a string in a field of width @n@ (pad on the left).
rAlign :: Int -> String -> String
rAlign n s = replicate ( max 0 ( n - length s ) ) ' ' ++ s
-- | Format a relative error in scientific notation, 2 significant figures.
showRelErr :: Double -> String
showRelErr x = showEFloat ( Just 2 ) x ""