diff --git a/brush-strokes/.gitignore b/brush-strokes/.gitignore index c8c6a76..101c5d6 100644 --- a/brush-strokes/.gitignore +++ b/brush-strokes/.gitignore @@ -1,3 +1,5 @@ dist-newstyle/ logs/ cabal.project.local +*.log +*.prof diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs index b723f46..4742526 100644 --- a/brush-strokes/src/cusps/bench/Main.hs +++ b/brush-strokes/src/cusps/bench/Main.hs @@ -1,6 +1,10 @@ +{-# LANGUAGE ScopedTypeVariables #-} + module Main where -- base +import Control.Arrow + ( second ) import Control.Monad ( when ) import Data.Foldable @@ -10,11 +14,13 @@ import qualified Data.List.NonEmpty as NE , fromList, head, length, sort ) import Data.Semigroup - ( Arg(..) ) + ( Arg(..), Sum(..), Product(..) ) import Data.Traversable ( for ) import GHC.Clock ( getMonotonicTime ) +import Numeric + ( showFFloat ) -- code-page import System.IO.CodePage @@ -23,6 +29,8 @@ import System.IO.CodePage -- containers import qualified Data.IntMap.Strict as IntMap ( fromList, toList ) +import Data.Tree + ( foldTree ) -- deepseq import Control.DeepSeq @@ -30,7 +38,12 @@ import Control.DeepSeq -- brush-strokes import Math.Bezier.Stroke +import Math.Interval +import Math.Linear import Math.Root.Isolation +import Math.Root.Isolation.Core +import Math.Root.Isolation.Newton +import Math.Root.Isolation.Newton.GaussSeidel -- brush-strokes bench:cusps import Bench.Cases @@ -66,10 +79,10 @@ benchTestCase testName ( TestCase { testDescription, testBrushStroke, testCuspOp , " --" ] before <- getMonotonicTime let ( _, testStrokeFnI ) = brushStrokeFunctions testBrushStroke - ( dunno, sols ) = + ( trees, dunno, sols ) = foldMap - ( \ ( i, ( _trees, DoneBoxes { doneSolBoxes = defCusps, doneGiveUpBoxes = mbCusps } ) ) -> - ( map ( ( i , ) . fst ) mbCusps, map ( i, ) defCusps ) ) $ + ( \ ( i, ( trees, DoneBoxes { doneSolBoxes = defCusps, doneGiveUpBoxes = mbCusps } ) ) -> + ( map ( i, ) trees, map ( ( i , ) . fst ) mbCusps, map ( i, ) defCusps ) ) $ IntMap.toList $ findCuspsIn testCuspOptions testStrokeFnI $ IntMap.fromList @@ -84,22 +97,124 @@ benchTestCase testName ( TestCase { testDescription, testBrushStroke, testCuspOp ++ ( map ( \ sol -> " -- • " ++ show sol ) sols ) ++ - [ " -- #dunno: " ++ show ( length dunno ) - , " -- Time elapsed: " ++ show dt ++ "s" ] + [ " -- #dunno: " ++ show dunno --( length dunno ) + , " -- Time elapsed: " ++ show dt ++ "s" + , " -- Tree size: " ++ show (getSum $ foldMap (foldMap sizeRootIsolationTree) trees) + , " -- Newton stats: " + ] ++ map ( " -- " ++ ) ( showNewtonStats (foldMap (foldMap newtonStats) trees) ) return dt +sizeRootIsolationTree :: ( Box 2, RootIsolationTree ( Box 2 ) ) -> Sum Int +sizeRootIsolationTree ( box, tree ) = foldMap ( const $ Sum 1 ) ( showRootIsolationTree box tree ) + + +data NewtonStats + = NewtonStats + { newtonImprovement :: !( Sum Double ) + , newtonElimination :: !( Sum Int, BigBoxes ) + , newtonTime :: !( Sum Double ) + , newtonTotal :: !( Sum Int ) + } +instance Semigroup NewtonStats where + NewtonStats imp1 el1 t1 tot1 <> NewtonStats imp2 el2 t2 tot2 = + NewtonStats ( imp1 <> imp2 ) ( el1 <> el2 ) ( t1 <> t2 ) ( tot1 <> tot2 ) +instance Monoid NewtonStats where + mempty = NewtonStats mempty mempty mempty mempty + +showNewtonStats :: NewtonStats -> [ String ] +showNewtonStats ( NewtonStats ( Sum improv ) ( Sum elims, big ) ( Sum time ) ( Sum tot ) ) + | tot == 0 + = [ ] + | otherwise + = [ "average improvement: " ++ showPercent ( improv / fromIntegral tot ) + , "average eliminations: " ++ showPercent ( fromIntegral elims / fromIntegral tot ) + , "time elapsed: " ++ show ( time / fromIntegral tot ) ++ "s" + ] + +showPercent :: Double -> String +showPercent x = fixed 3 2 ( 100 * x ) <> "%" + +fixed :: Int -> Int -> Double -> String +fixed digitsBefore digitsAfter x = + case second ( drop 1 ) . break ( == '.' ) $ showFFloat ( Just digitsAfter ) x "" of + ( as, bs ) -> + let + l, r :: Int + l = length as + r = length bs + in + replicate ( digitsBefore - l ) ' ' <> as <> "." <> bs <> replicate ( digitsAfter - r ) '0' + +data BigBoxes = + BigBoxes + { biggestArea :: Maybe ( Box 2 ) + , biggestX :: Maybe ( Box 2 ) + , biggestY :: Maybe ( Box 2 ) + } + deriving stock ( Show, Eq ) +instance Semigroup BigBoxes where + BigBoxes area1 x1 y1 <> BigBoxes area2 x2 y2 = + BigBoxes + ( pickBig boxArea area1 area2 ) + ( pickBig widthX x1 x2 ) + ( pickBig widthY y1 y2 ) + where + pickBig _ Nothing Nothing = Nothing + pickBig _ (Just x1) Nothing = Just x1 + pickBig _ Nothing (Just x2) = Just x2 + pickBig f j1@(Just x1) j2@(Just x2) + | f x1 >= f x2 + = j1 + | otherwise + = j2 + widthX ( 𝕀 ( ℝ2 x_lo _y_lo ) ( ℝ2 x_hi _y_hi ) ) + = x_hi - x_lo + widthY ( 𝕀 ( ℝ2 _x_lo y_lo ) ( ℝ2 _x_hi y_hi ) ) + = y_hi - y_lo +instance Monoid BigBoxes where + mempty = BigBoxes Nothing Nothing Nothing + + +newtonStats :: ( Box 2, RootIsolationTree ( Box 2 ) ) -> NewtonStats +newtonStats ( _box, RootIsolationLeaf {} ) = mempty +newtonStats ( box, RootIsolationStep step boxes ) + | IsolationStep @Newton ( TimeInterval dt ) <- step + = thisImprovement dt <> foldMap newtonStats boxes + | otherwise + = foldMap newtonStats boxes + where + thisImprovement dt = + NewtonStats + { newtonImprovement = Sum $ ( old - new ) / old + , newtonElimination = + if new == 0 + then ( Sum 1, BigBoxes ( Just box ) ( Just box ) ( Just box ) ) + else mempty + , newtonTime = Sum dt + , newtonTotal = Sum 1 + } + where + old = boxArea box + new = sum ( map ( boxArea . fst ) boxes ) + benchGroups :: [ ( String, NE.NonEmpty TestCase ) ] benchGroups = [ ( "ellipse" , NE.fromList - [ ellipseTestCase opts ("ε_bis=" ++ show ε_bis) + [ ellipseTestCase opts newtMeth ( 0, 1 ) pi ( defaultStartBoxes [ 0 .. 3 ] ) - | ε_bis <- [ 1e-6, 5e-6, 1e-5, 5e-5, 1e-4, 5e-4, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 0.1, 0.2, 0.3 ] + | let ε_bis = 5e-3 -- <- [ 1e-6, 5e-6, 1e-5, 5e-5, 1e-4, 5e-4, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 0.1, 0.2, 0.3 ] + , (newtMeth, newtOpts) + <- [ ( "LP", \ hist -> NewtonLP ) + , ( "GS_Complete", defaultNewtonOptions @2 @3 ) + , ( "GS_Partial", \hist -> NewtonGaussSeidel $ ( defaultGaussSeidelOptions @2 @3 hist ) { gsUpdate = GS_Partial } ) + ] , let opts = RootIsolationOptions - { rootIsolationAlgorithms = + { rootIsolationAlgorithms = \ hist -> defaultRootIsolationAlgorithms minWidth ε_bis + ( newtOpts hist ) hist } ] ) diff --git a/brush-strokes/src/lib/Math/Linear.hs b/brush-strokes/src/lib/Math/Linear.hs index bd4aa1a..95a1d6c 100644 --- a/brush-strokes/src/lib/Math/Linear.hs +++ b/brush-strokes/src/lib/Math/Linear.hs @@ -45,7 +45,7 @@ import Data.Act -- deepseq import Control.DeepSeq - ( NFData, NFData1 ) + ( NFData(..), NFData1 ) -- groups import Data.Group @@ -139,6 +139,8 @@ deriving via ZipList instance Applicative ( Vec n ) instance Traversable ( Vec n ) where traverse f ( Vec as ) = Vec <$> traverse f as +instance NFData a => NFData ( Vec n a ) where + rnf ( Vec as ) = rnf as universe :: forall n. KnownNat n => Vec n ( Fin n ) universe = Vec [ Fin i | i <- [ 1 .. fromIntegral ( natVal' @n proxy# ) ] ] diff --git a/brush-strokes/src/lib/Math/Root/Isolation.hs b/brush-strokes/src/lib/Math/Root/Isolation.hs index 9fc4d9f..70f93a6 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation.hs @@ -102,10 +102,12 @@ newtype RootIsolationOptions n d -> Either String ( NE.NonEmpty ( RootIsolationAlgorithmWithOptions n d ) ) } -defaultRootIsolationOptions :: BoxCt n d => RootIsolationOptions n d +defaultRootIsolationOptions :: forall n d. BoxCt n d => RootIsolationOptions n d defaultRootIsolationOptions = RootIsolationOptions - { rootIsolationAlgorithms = defaultRootIsolationAlgorithms minWidth ε_eq + { rootIsolationAlgorithms = \ hist -> + defaultRootIsolationAlgorithms minWidth ε_eq + ( defaultNewtonOptions @n @d hist ) hist } where minWidth = 1e-5 @@ -117,10 +119,11 @@ defaultRootIsolationAlgorithms . BoxCt n d => Double -- ^ minimum width of boxes (don't bisect further) -> Double -- ^ threshold for progress + -> NewtonOptions n d -- ^ options for Newton's method -> BoxHistory n -> Box n -> Either String ( NE.NonEmpty ( RootIsolationAlgorithmWithOptions n d ) ) -defaultRootIsolationAlgorithms minWidth ε_eq history box = +defaultRootIsolationAlgorithms minWidth ε_eq newtonOptions history box = case history of lastRoundBoxes : _ -- If, in the last round of strategies, we didn't try bisection... @@ -137,14 +140,13 @@ defaultRootIsolationAlgorithms minWidth ε_eq history box = -- Currently: we try an interval Gauss–Seidel. -- (box(1)- and box(2)-consistency don't seem to help when using -- the complete interval union Gauss–Seidel step) - _ -> Right $ AlgoWithOptions @Newton _newtonOptions + _ -> Right $ AlgoWithOptions @Newton newtonOptions NE.:| [] where verySmall = and $ ( \ cd -> width cd <= minWidth ) <$> coordinates box _bisOptions = defaultBisectionOptions @n @d minWidth ε_eq box - _newtonOptions = NewtonLP -- defaultNewtonOptions @n @d history _box1Options = defaultBox1Options @n @d ( minWidth * 100 ) ε_eq _box2Options = ( defaultBox2Options @n @d minWidth ε_eq ) { box2LambdaMin = 0.001 } @@ -200,7 +202,7 @@ isolateRootsIn ( RootIsolationOptions { rootIsolationAlgorithms } ) | -- Check the range of the equations contains zero. not $ ( unT ( origin @Double ) ∈ iRange ) -- Box doesn't contain a solution: discard it. - = return [] + = return [ RootIsolationLeaf "rangeTest" cand ] | otherwise = case rootIsolationAlgorithms history cand of Right strats -> doStrategies history strats cand diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Core.hs b/brush-strokes/src/lib/Math/Root/Isolation/Core.hs index 901443b..d436724 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation/Core.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation/Core.hs @@ -22,11 +22,15 @@ module Math.Root.Isolation.Core -- ** Visualising history , RootIsolationTree(..) + , boxArea , showRootIsolationTree -- * Utility functions , pipeFunctionsWhileTrue , forEachCoord + + -- * Timing + , TimeInterval(..), timeInterval ) where -- base @@ -40,15 +44,23 @@ import Data.Type.Equality ( (:~~:)(HRefl) ) import Data.Typeable ( Typeable, heqT ) -import GHC.TypeNats - ( Nat, KnownNat, type (<=) ) import Numeric ( showFFloat ) +import GHC.Clock + ( getMonotonicTime ) +import GHC.TypeNats + ( Nat, KnownNat, type (<=) ) +import System.IO.Unsafe + ( unsafePerformIO ) -- containers import Data.Tree ( Tree(..) ) +-- deepseq +import Control.DeepSeq + ( NFData(..), deepseq ) + -- transformers import Control.Monad.Trans.State.Strict as State ( State, get, put ) @@ -78,8 +90,7 @@ type Box n = 𝕀ℝ n -- NB: we require n <= d (no support for under-constrained systems). -- -- NB: in practice, this constraint should specialise away. -type BoxCt n d = ( n ~ 2, d ~ 3 ) -{- +type BoxCt n d = ( KnownNat n, KnownNat d , 1 <= n, 1 <= d, n <= d @@ -91,12 +102,13 @@ type BoxCt n d = ( n ~ 2, d ~ 3 ) , Vars ( D 1 ( ℝ n ) ) ~ n , Module Double ( T ( ℝ n ) ) , Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) + , NFData ( 𝕀ℝ n ) , Ord ( ℝ d ) , Module Double ( T ( ℝ d ) ) , Representable Double ( ℝ d ) + , NFData ( ℝ d ) ) --} -- | Boxes we are done with and will not continue processing. data DoneBoxes n = DoneBoxes @@ -272,3 +284,14 @@ pipeFunctionsWhileTrue fns = go fns concat <$> traverse ( go fs ) xs -------------------------------------------------------------------------------- + +newtype TimeInterval = TimeInterval Double + deriving newtype Show + +timeInterval :: NFData b => ( a -> b ) -> a -> ( b, TimeInterval ) +timeInterval f a = unsafePerformIO $ do + bef <- getMonotonicTime + let !b = f a + b `deepseq` return () + aft <- getMonotonicTime + return $ ( b, TimeInterval ( aft - bef ) ) diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Newton.hs b/brush-strokes/src/lib/Math/Root/Isolation/Newton.hs index a6a4072..acd6760 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation/Newton.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation/Newton.hs @@ -26,6 +26,10 @@ import Data.List import GHC.TypeNats ( Nat, KnownNat, type (<=) ) +-- deepseq +import Control.DeepSeq + ( force ) + -- transformers import Control.Monad.Trans.Writer.CPS ( Writer, tell ) @@ -50,11 +54,10 @@ import Math.Root.Isolation.Utils -- | The interval Newton method; see 'intervalNewton'. data Newton instance BoxCt n d => RootIsolationAlgorithm Newton n d where - type instance StepDescription Newton = () + type instance StepDescription Newton = TimeInterval type instance RootIsolationAlgorithmOptions Newton n d = NewtonOptions n d - rootIsolationAlgorithm opts _thisRoundHist _prevRoundsHist eqs box = do - res <- intervalNewton @n @d opts eqs box - return ( (), res ) + rootIsolationAlgorithm opts _thisRoundHist _prevRoundsHist eqs box = + intervalNewton @n @d opts eqs box {-# INLINEABLE rootIsolationAlgorithm #-} {-# SPECIALISE rootIsolationAlgorithm :: RootIsolationAlgorithmOptions Newton 2 3 @@ -99,7 +102,7 @@ intervalNewton -- ^ equations -> 𝕀ℝ n -- ^ box - -> Writer ( DoneBoxes n ) [ 𝕀ℝ n ] + -> Writer ( DoneBoxes n ) ( TimeInterval, [ 𝕀ℝ n ] ) intervalNewton opts eqs x = case opts of NewtonGaussSeidel ( GaussSeidelOptions @@ -118,13 +121,18 @@ intervalNewton opts eqs x = case opts of minus_f_x_mid = unT $ -1 *^ T ( boxMidpoint $ f x_mid ) -- Precondition the above linear system into A ( x - x_mid ) = B. - ( a, b ) = precondition precondMeth - f'_x ( singleton minus_f_x_mid ) + !( !a, !b ) = force $ + precondition precondMeth + f'_x ( singleton minus_f_x_mid ) + !x'_0 = force ( T x ^-^ T x_mid ) -- NB: we have to change coordinates, putting the midpoint of the box -- at the origin, in order to take a Gauss–Seidel step. - gsGuesses = map ( first ( \ x' -> unT $ x' ^+^ T x_mid ) ) - $ gaussSeidelUpdate gsUpdate a b ( T x ^-^ T x_mid ) + ( x's, dt ) = + timeInterval + ( gaussSeidelUpdate gsUpdate a b ) + x'_0 + gsGuesses = map ( first ( \ x' -> unT $ x' ^+^ T x_mid ) ) x's ( done, todo ) = bimap ( map fst ) ( map fst ) $ partition snd gsGuesses in -- If the Gauss–Seidel step was a contraction, then the box @@ -133,7 +141,7 @@ intervalNewton opts eqs x = case opts of -- These boxes can thus be directly added to the solution set: -- Newton's method is guaranteed to converge to the unique solution. do tell $ noDoneBoxes { doneSolBoxes = done } - return todo + return ( dt, todo ) NewtonLP -> -- TODO: reduce duplication with the above. let x_mid = singleton $ boxMidpoint x @@ -143,13 +151,16 @@ intervalNewton opts eqs x = case opts of f'_x = fmap ( \ i -> eqs x `monIndex` linearMonomial i ) ( universe @2 ) minus_f_x_mid = unT $ -1 *^ T ( boxMidpoint $ f x_mid ) - ( a, b ) = ( f'_x, singleton minus_f_x_mid ) - lpGuesses = map ( first ( \ x' -> unT $ x' ^+^ T x_mid ) ) - $ solveIntervalLinearEquations a b ( T x ^-^ T x_mid ) + !( !a, !b ) = force ( f'_x, singleton minus_f_x_mid ) + !x'_0 = force ( T x ^-^ T x_mid ) + ( x's, dt ) = + timeInterval + ( solveIntervalLinearEquations a b ) x'_0 + lpGuesses = map ( first ( \ x' -> unT $ x' ^+^ T x_mid ) ) x's ( done, todo ) = bimap ( map fst ) ( map fst ) $ partition snd lpGuesses in do tell $ noDoneBoxes { doneSolBoxes = done } - return todo + return ( dt, todo ) {-# INLINEABLE intervalNewton #-} {- diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Newton/LP.hs b/brush-strokes/src/lib/Math/Root/Isolation/Newton/LP.hs index fe2802d..61419a1 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation/Newton/LP.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation/Newton/LP.hs @@ -42,7 +42,7 @@ import Math.Linear -- | Solve the system of linear equations \( A X = B \) -- using linear programming. solveIntervalLinearEquations - :: ( KnownNat d, Representable Double ( ℝ d ), Show ( ℝ d ) ) + :: ( KnownNat d, Representable Double ( ℝ d ) ) => Vec 2 ( 𝕀ℝ d ) -- ^ columns of \( A \) -> 𝕀ℝ d -- ^ \( B \) -> T ( 𝕀ℝ 2 ) -- ^ initial box \( X \)