add more info to cusps benchmark output

brush-strokes/.gitignore

@@ -1,3 +1,5 @@
 dist-newstyle/
 logs/
 cabal.project.local
+*.log
+*.prof

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
                 }
       ]
     )

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# ) ] ]

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

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

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
+        !( !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 #-}
 {-
 

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 \)