make cusp-finding algorithm choice more configurable

brush-strokes/brush-strokes.cabal

@@ -250,7 +250,8 @@ benchmark cusps
     Main.hs
 
   build-depends:
-    brush-strokes
+      brush-strokes
+    , deepseq
 
   type:
     exitcode-stdio-1.0

brush-strokes/src/cusps/bench/Main.hs

@@ -25,7 +25,11 @@ import Data.Coerce
 import Data.Foldable
   ( for_ )
 import Data.List
-  ( intercalate )
+  ( intercalate, sortOn )
+import Data.Traversable
+  ( for )
+import GHC.Clock
+  ( getMonotonicTime )
 import GHC.Exts
   ( Proxy#, proxy# )
 import GHC.Generics
@@ -43,6 +47,13 @@ import qualified Data.Sequence as Seq
 import Data.Tree
   ( foldTree )
 
+-- deepseq
+import Control.DeepSeq
+  ( rnf )
+
+-- gauge
+--import qualified Gauge
+
 -- tree-view
 import Data.Tree.View
   ( showTree )
@@ -64,11 +75,66 @@ import Math.Ring
 
 --------------------------------------------------------------------------------
 
+{-
+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
-  for_ testCases $ \ testCase@( TestCase { testName, testBrushStroke, testAlgorithmParams, testStartBoxes } ) -> 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 (testCuspOptions bestTest)
+    ]
+
+benchTestCase :: TestCase -> IO Double
+benchTestCase ( TestCase { testName, testBrushStroke, testCuspOptions, testStartBoxes } ) = do
+  putStrLn $ unlines
+    [ replicate 40 '='
+    , "Test case " ++ testName ]
+  before <- getMonotonicTime
+  let ( _, testStrokeFnI ) = brushStrokeFunctions testBrushStroke
+      ( _newtTrees, ( dunno, sols ) ) = computeCusps testCuspOptions testStrokeFnI testStartBoxes
+  rnf dunno `seq` rnf sols `seq` return ()
+  after <- getMonotonicTime
+  let dt = after - before
+  putStrLn $ unlines
+    [ "  -  #sols: " ++ show (length 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 testAlgorithmParams testStrokeFnI testStartBoxes
+        ( newtTrees, ( dunno, sols ) ) = computeCusps testCuspOptions testStrokeFnI testStartBoxes
         showedTrees = map ( uncurry showIntervalNewtonTree ) newtTrees
         testHeader =
           [ "", "Test case '" ++ testName ++ "':" ]
@@ -84,6 +150,7 @@ main = do
   --  logFileMVar <- newMVar "logs/fnData.log"
   --  logToFile logFileMVar (unlines logLines)
   --    `seq` return ()
+-}
 
 testCases :: [ TestCase ]
 testCases = benchCases
@@ -103,7 +170,12 @@ testCases = benchCases
 -}
 
 benchCases :: [ TestCase ]
-benchCases = [ ellipseTestCase "full" ( 0, 1 ) pi $ defaultStartBoxes [ 0 .. 3 ] ]
+benchCases =
+  [ ellipseTestCase opts "full" ( 0, 1 ) pi $ defaultStartBoxes [ 0 .. 3 ] ]
+  where
+    opts = defaultCuspFindingOptions
+
+
 
 --------------------------------------------------------------------------------
 
@@ -117,10 +189,10 @@ data BrushStroke =
 
 data TestCase =
   TestCase
-  { testName            :: String
+  { testName        :: String
-  , testBrushStroke     :: BrushStroke
+  , testBrushStroke :: BrushStroke
-  , testAlgorithmParams :: CuspAlgorithmParams
+  , testCuspOptions :: CuspFindingOptions
-  , testStartBoxes      :: [ Box ]
+  , testStartBoxes  :: [ Box ]
   }
 
 brushStrokeFunctions
@@ -176,7 +248,7 @@ take 10 $ Data.List.sortOn ( \ ( _, ℝ1 e, v) -> abs e + norm v ) [ let { v = m
 potentialCusp $ eval fI $ mkBox (0.5798, 0.5799) 3 (0.26798, 0.26799)
 > True
 
-let nbPotentialSols b = let ( _newtTrees, ( dunno, sols ) ) = intervalNewtonGSFrom NoPreconditioning 1e-7 fI b in length dunno + length sols
+let nbPotentialSols b = let ( _newtTrees, ( dunno, sols ) ) = findCuspsFrom NoPreconditioning 1e-7 fI b in length dunno + length sols
 
 nbPotentialSols $ mkBox (0.5798, 0.5799) 3 (0.26798, 0.26799)
 1
@@ -184,7 +256,7 @@ nbPotentialSols $ mkBox (0.5798, 0.5799) 3 (0.26798, 0.26799)
 nbPotentialSols $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799)
 0
 
-let showTrees b = map ( uncurry showIntervalNewtonTree ) $ fst $ intervalNewtonGSFrom NoPreconditioning 1e-7 fI b
+let showTrees b = map ( uncurry showIntervalNewtonTree ) $ fst $ findCuspsFrom NoPreconditioning 1e-7 fI b
 
 putStrLn $ unlines $ map Data.Tree.View.showTree $ showTrees $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799)
 
@@ -390,7 +462,7 @@ maximum [ _y $ sup $ unT $ _D12_v $ dEdsdcdt $ eval fI $ mkBox (t, t) 2 (s, s) |
 
 
 
-let showTrees b = map ( uncurry showIntervalNewtonTree ) $ fst $ intervalNewtonGSFrom NoPreconditioning 1e-7 fI b
+let showTrees b = map ( uncurry showIntervalNewtonTree ) $ fst $ findCuspsFrom NoPreconditioning 1e-7 fI b
 
 putStrLn $ unlines $ map Data.Tree.View.showTree $ showTrees $ mkBox (0.5486101933245248, 0.5486102071493622) 2 (0.548095036738487, 0.5480952)
 
@@ -408,16 +480,12 @@ showD float = showFFloat (Just 6) float ""
 
 --------------------------------------------------------------------------------
 
-ellipseTestCase :: String -> ( Double, Double ) -> Double -> [ Box ] -> TestCase
+ellipseTestCase :: CuspFindingOptions -> String -> ( Double, Double ) -> Double -> [ Box ] -> TestCase
-ellipseTestCase str k0k1 rot startBoxes =
+ellipseTestCase opts str k0k1 rot startBoxes =
   TestCase
     { testName = "ellipse (" ++ str ++ ")"
     , testBrushStroke = ellipseBrushStroke k0k1 rot
-    , testAlgorithmParams =
+    , testCuspOptions = opts
-        CuspAlgorithmParams
-         { preconditioning = InverseMidJacobian
-         , maxWidth = 1e-7
-         }
     , testStartBoxes = startBoxes
     }
 
@@ -442,11 +510,7 @@ trickyCusp2TestCase =
   TestCase
     { testName = "trickyCusp2"
     , testBrushStroke = trickyCusp2BrushStroke
-    , testAlgorithmParams =
+    , testCuspOptions = defaultCuspFindingOptions
-        CuspAlgorithmParams
-         { preconditioning = InverseMidJacobian
-         , maxWidth = 1e-7
-         }
     , testStartBoxes = defaultStartBoxes [ 0 .. 3 ]
     }
 
@@ -496,13 +560,6 @@ data Point nbParams =
   deriving stock Generic
 deriving stock instance Show ( ℝ nbParams ) => Show ( Point nbParams )
 
-data CuspAlgorithmParams =
-  CuspAlgorithmParams
-    { preconditioning :: !Preconditioner
-    , maxWidth :: !Double
-    }
-  deriving stock Show
-
 type Brush nbParams
   = forall {t} k (i :: t)
   .  DiffInterp k i ( ℝ nbParams )
@@ -543,14 +600,12 @@ getStrokeFunctions brush sp0 crv =
 {-# INLINEABLE getStrokeFunctions #-}
 
 computeCusps
-  :: CuspAlgorithmParams
+  :: CuspFindingOptions
   -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
   -> [ Box ]
-  -> ( [ ( Box, IntervalNewtonTree Box ) ], ( [ Box ], [ Box ] ) )
+  -> ( [ ( Box, RootFindingTree Box ) ], ( [ Box ], [ Box ] ) )
-computeCusps params eqs startBoxes =
+computeCusps cuspOpts eqs startBoxes =
-  foldMap
+  foldMap ( findCuspsFrom cuspOpts eqs ) startBoxes
-    ( intervalNewtonGSFrom ( preconditioning params ) ( maxWidth params ) eqs )
-    startBoxes
 
 defaultStartBoxes :: [ Int ] -> [ Box ]
 defaultStartBoxes is =
@@ -561,8 +616,8 @@ getR1 (ℝ1 u) = u
 {-
 
 (f, fI) = brushStrokeFunctions $ ellipseBrushStroke (0,1) pi
-nbPotentialSols box = let ( _newtTrees, ( dunno, sols ) ) = intervalNewtonGSFrom NoPreconditioning 1e-7 fI box in length dunno + length sols
+nbPotentialSols box = let ( _newtTrees, ( dunno, sols ) ) = findCuspsFrom NoPreconditioning 1e-7 fI box in length dunno + length sols
-showTrees box = putStrLn $ unlines $ map Data.Tree.View.showTree $ map ( uncurry showIntervalNewtonTree ) $ fst $ intervalNewtonGSFrom NoPreconditioning 1e-7 fI box
+showTrees box = putStrLn $ unlines $ map Data.Tree.View.showTree $ map ( uncurry showIntervalNewtonTree ) $ fst $ findCuspsFrom NoPreconditioning 1e-7 fI box
 
 sol_t = 0.5486100729150693796677845183880669025324233347060776339185 :: Double
 sol_s = 0.5480950141859386853197594577293968665598143630958601978245 :: Double

brush-strokes/src/lib/Math/Bezier/Stroke.hs

@@ -19,13 +19,14 @@ module Math.Bezier.Stroke
   , brushStrokeData, pathAndUsedParams
 
   -- ** Cusp finding
-  , Preconditioner(..)
+  , RootFindingTree(..), showRootFindingTree
-  , IntervalNewtonTree(..), showIntervalNewtonTree
+  , RootFindingStep(..)
-  , IntervalNewtonStep(..)
+  , RootFindingLeaf(..)
-  , IntervalNewtonLeaf(..)
   , Box
-  , intervalNewtonGS, intervalNewtonGSFrom, gaussSeidel
+  , CuspFindingOptions(..), defaultCuspFindingOptions
-
+  , RootFindingAlgorithm(..), defaultCuspFindingAlgorithms
+  , Preconditioner(..)
+  , findCusps, findCuspsFrom, gaussSeidel
   )
   where
 
@@ -53,13 +54,11 @@ import Data.List
 import Data.List.NonEmpty
   ( NonEmpty )
 import qualified Data.List.NonEmpty as NE
-  ( cons, singleton, unzip )
+  ( NonEmpty(..), cons, last, singleton, unzip )
 import Data.Maybe
   ( fromMaybe, isJust, listToMaybe, mapMaybe )
 import Data.Semigroup
   ( sconcat )
-import Data.Traversable
-  ( for )
 import GHC.Exts
   ( newMutVar#, runRW#, inline
   , Proxy#, proxy#
@@ -263,6 +262,7 @@ computeStrokeOutline ::
 
      )
   => RootSolvingAlgorithm
+  -> Maybe CuspFindingOptions
   -> FitParameters
   -> ( ptData -> usedParams )
   -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing
@@ -279,7 +279,7 @@ computeStrokeOutline ::
         , Seq FitPoint
         , [ Cusp ]
         )
-computeStrokeOutline rootAlgo fitParams ptParams toBrushParams brushFn spline@( Spline { splineStart = spt0 } ) = case ssplineType @clo of
+computeStrokeOutline rootAlgo mbCuspOptions fitParams ptParams toBrushParams brushFn spline@( Spline { splineStart = spt0 } ) = case ssplineType @clo of
   -- Open brush path with at least one segment.
   -- Need to add caps at both ends of the path.
   SOpen
@@ -381,7 +381,7 @@ computeStrokeOutline rootAlgo fitParams ptParams toBrushParams brushFn spline@(
   where
 
     outlineInfo :: ptData -> Curve Open crvData ptData -> OutlineInfo
-    outlineInfo = inline ( outlineFunction rootAlgo ptParams toBrushParams brushFn )
+    outlineInfo = inline ( outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams brushFn )
 
     outlineFns :: Seq OutlineInfo
     outlineFns = go spt0 ( openCurves $ splineCurves ( adjustSplineType @Open spline ) )
@@ -541,6 +541,7 @@ outlineFunction
      , Show ptData, Show brushParams
      )
   => RootSolvingAlgorithm
+  -> Maybe CuspFindingOptions
   -> ( ptData -> usedParams )
   -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing
   -> ( forall {t} k (i :: t)
@@ -553,7 +554,7 @@ outlineFunction
   -> ptData
   -> Curve Open crvData ptData
   -> OutlineInfo
-outlineFunction rootAlgo ptParams toBrushParams brushFromParams = \ sp0 crv ->
+outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams brushFromParams = \ sp0 crv ->
   let
 
     usedParams :: C 2 ( ℝ 1 ) usedParams
@@ -604,10 +605,11 @@ outlineFunction rootAlgo ptParams toBrushParams brushFromParams = \ sp0 crv ->
                              $ toBrushParams params_t
 
     ( newtTrees, ( newtDunno, newtSols ) ) =
-      intervalNewtonGS
+      case mbCuspOptions of
-        InverseMidJacobian
+        Just cuspOptions ->
-        1e-7
+          findCusps cuspOptions curvesI
-        curvesI
+        Nothing ->
+          ( [], ( [], [] ) )
 
     showD :: Double -> String
     showD float = showFFloat (Just 6) float ""
@@ -627,7 +629,7 @@ outlineFunction rootAlgo ptParams toBrushParams brushFromParams = \ sp0 crv ->
          ]
         ) ++ "}"
       ]
-    showedTrees = map ( uncurry showIntervalNewtonTree ) newtTrees
+    showedTrees = map ( uncurry showRootFindingTree ) newtTrees
     solLines =
       [ " #sols: " ++ show (length newtSols)
       , "#dunno: " ++ show (length newtDunno)
@@ -638,7 +640,7 @@ outlineFunction rootAlgo ptParams toBrushParams brushFromParams = \ sp0 crv ->
 
     logContents = unlines $ functionDataLines ++ treeLines
 
-  in trace (unlines solLines) $
+  in --trace (unlines solLines) $
      --logToFile cuspFindingMVar logContents `seq`
      OutlineInfo
        { outlineFn = fwdBwd
@@ -1032,16 +1034,6 @@ solveEnvelopeEquations rootAlgo _t path_t path'_t ( fwdOffset, bwdOffset ) strok
 
   where
 
-    -- !_ = trace
-    --      ( unlines
-    --        [ "solveEnvelopeEquation"
-    --        , "  t: " ++ show _t
-    --        , "  pt: " ++ show path_t
-    --        , "  tgt: "  ++ show path'_t
-    --        , "  fwdOffset: " ++ show fwdOffset
-    --        , "  bwdOffset: " ++ show bwdOffset ]
-    --      ) True
-
     fwdSol            = findSolFrom "fwd" fwdOffset
     ( bwdPt, bwdTgt ) = findSolFrom "bwd" bwdOffset
 
@@ -1278,115 +1270,9 @@ gaussSeidel
 
         return ( ( 𝕀 ( ℝ1 x1'_lo ) ( ℝ1 x1'_hi ), 𝕀 ( ℝ1 x2'_lo ) ( ℝ1 x2'_hi ) )
                , sub_x1 && sub_x2 )
-
+-- TODO: try implementing the complete interval union Gauss–Seidel algorithm.
-{-
+-- See "Algorithm 2" in
-gaussSeidel2 :: Int
+-- "Using interval unions to solve linear systems of equations with uncertainties"
-             -> Double
-             -> Double
-             -> ( 𝕀ℝ 2, 𝕀ℝ 2 ) -- ^ columns of \( A \)
-             -> 𝕀ℝ 2           -- ^ \( B \)
-             -> ( 𝕀ℝ 1, 𝕀ℝ 1 ) -- ^ initial box \( X \)
-             -> [ ( ( 𝕀ℝ 1, 𝕀ℝ 1 ), Bool ) ]
-gaussSeidel2 maxIters eps_abs eps_rel
-  ( 𝕀 ( ℝ2 a11_lo a21_lo ) ( ℝ2 a11_hi a21_hi )
-  , 𝕀 ( ℝ2 a12_lo a22_lo ) ( ℝ2 a12_hi a22_hi ) )
-  ( 𝕀 ( ℝ2 b1_lo b2_lo ) ( ℝ2 b1_hi b2_hi ) )
-  x0
-    = let !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
-
-          -- See "Algorithm 2" in
-          --  "Using interval unions to solve linear systems of equations with uncertainties"
-          iter ( 𝕀 ( ℝ1 x1_lo ) ( ℝ1 x1_hi ), 𝕀 ( ℝ1 x2_lo ) ( ℝ1 x2_hi ) ) = do
-            let
-              !x1  = 𝕀 x1_lo x1_hi
-              !x2  = 𝕀 x2_lo x2_hi
-              blah1 = do
-                let s = b1 - a12 * x2
-                let s1 = s `monus` ( a11 * x1 )
-                y1 <-
-                  if not $ containsZero ( s1 - a11 * x1 )
-                  then []
-                  else
-                    if containsZero s1 && containsZero a11
-                    then [ ( x1, False ) ]
-                    else do
-                      x1'' <- s1 `extendedDivide` a11
-                      x1'' `intersect` x1
-
-                let s2 = s `monus` ( a12 * x2 )
-                y2 <-
-                  if not $ containsZero ( s2 - a11 * x1 )
-                  then []
-                  else
-                    if containsZero s2 && containsZero a11
-                    then [ ( x1, False ) ]
-                    else do
-                      x1'' <- s2 `extendedDivide` a12
-                      x1'' `intersect` x1
-                return ( y1 `cart` y2 )
-
-              blah2 = do
-                let s = b2 - a21 * x1
-                let s1 = s `monus` ( a21 * x1 )
-                y1 <-
-                  if not $ containsZero ( s1 - a22 * x2 )
-                  then []
-                  else
-                    if containsZero s1 && containsZero a22
-                    then [ ( x1, False ) ]
-                    else do
-                      x1'' <- s1 `extendedDivide` a21
-                      x1'' `intersect` x1
-
-                let s2 = s `monus` ( a12 * x2 )
-                y2 <-
-                  if not $ containsZero ( s2 - a22 * x2 )
-                  then []
-                  else
-                    if containsZero s2 && containsZero a22
-                    then [ ( x1, False ) ]
-                    else do
-                      x1'' <- s2 `extendedDivide` a22
-                      x1'' `intersect` x1
-                return ( y1 `cart` y2 )
-
-            blah1 ++ blah2
-
-          go :: Int -> ( 𝕀ℝ 1, 𝕀ℝ 1 ) -> [ ( ( 𝕀ℝ 1, 𝕀ℝ 1 ), Bool ) ]
-          go !i x
-            = do { nxt@( x', sub ) <- iter x
-                 ; let dw_abs = maxWidth x - maxWidth x'
-                       dw_rel = 1 - ( maxWidth x' / maxWidth x )
-                 ; if sub
-                   || i >= maxIters
-                   || dw_abs < eps_abs
-                   || dw_rel < eps_rel
-                   then return nxt
-                   else go ( i + 1 ) x'
-                 }
-
-      in go 1 x0
-  where
-    maxWidth :: ( 𝕀ℝ 1, 𝕀ℝ 1 ) -> Double
-    maxWidth ( 𝕀 ( ℝ1 x1_lo ) ( ℝ1 x1_hi ), 𝕀 ( ℝ1 x2_lo ) ( ℝ1 x2_hi ) )
-      = max ( x1_hi - x1_lo ) ( x2_hi - x2_lo )
-    containsZero :: 𝕀 Double -> Bool
-    containsZero ( 𝕀 lo hi ) = lo <= 0 && hi >= 0
-    monus :: 𝕀 Double -> 𝕀 Double -> 𝕀 Double
-    monus ( 𝕀 lo1 hi1 ) ( 𝕀 lo2 hi2 )
-      | hi1 - lo1 >= hi2 - lo2
-      = 𝕀 ( lo1 - lo2 ) ( hi1 - hi2 )
-      | otherwise
-      = 𝕀 ( hi1 - hi2 ) ( lo1 - lo2 )
-    cart :: ( 𝕀 Double, Bool ) -> ( 𝕀 Double, Bool ) -> ( ( 𝕀ℝ 1, 𝕀ℝ 1 ), Bool )
-    cart ( 𝕀 lo1 hi1, sub1 ) ( 𝕀 lo2 hi2, sub2 ) =
-      ( ( 𝕀 ( ℝ1 lo1 ) ( ℝ1 hi1 ), 𝕀 ( ℝ1 lo2 ) ( ℝ1 hi2 ) ), sub1 && sub2 )
--}
 
 -- | Compute the intersection of two intervals.
 --
@@ -1404,6 +1290,9 @@ intersect ( 𝕀 lo1 hi1 ) ( 𝕀 lo2 hi2 )
 
 -- | Computes the brush stroke coordinates of a cusp from
 -- the @(t,s)@ parameter values.
+--
+-- TODO: use Newton's method starting at the midpoint of the box,
+-- instead of just taking the midpoint of the box.
 cuspCoords :: ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) )
            -> Box
            -> Cusp
@@ -1421,42 +1310,45 @@ cuspCoords eqs ( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ), i, 𝕀 ( ℝ1 s_lo ) ( ℝ1
     t_mid = 0.5 * ( t_lo + t_hi )
     s_mid = 0.5 * ( s_lo + s_hi )
 
-
+-- | Preconditioner to use with the interval Gauss–Seidel method.
 data Preconditioner
   = NoPreconditioning
   | InverseMidJacobian
   deriving stock Show
 
--- | A tree recording the steps taken with the interval Newton method.
+-- | A tree recording the steps taken when doing cusp finding.
-data IntervalNewtonTree d
+data RootFindingTree d
-  = IntervalNewtonLeaf (IntervalNewtonLeaf d)
+  = RootFindingLeaf (RootFindingLeaf d)
-  | IntervalNewtonStep (IntervalNewtonStep d) [(d, IntervalNewtonTree d)]
+  | RootFindingStep RootFindingStep [(d, RootFindingTree d)]
 
-data IntervalNewtonStep d
+-- | A description of a step taken in cusp-finding.
-  = IntervalNewtonContraction [d]
+data RootFindingStep
-  | IntervalNewtonBisection (String, Double)
+  = GaussSeidelStep
+  | BisectionStep (String, Double)
+  | Box1Step
+  | Box2Step
 
-instance Show d => Show (IntervalNewtonStep d) where
+instance Show RootFindingStep where
-  showsPrec _ ( IntervalNewtonContraction v )
+  showsPrec _ GaussSeidelStep = showString "GS"
-    = showString "N "
+  showsPrec _ ( BisectionStep (coord, w) )
-    . showsPrec 0 v
+    = showString "bis "
-  showsPrec _ ( IntervalNewtonBisection (coord, w) )
-    = showString "B "
     . showParen True
       ( showString coord
       . showString " = "
       . showsPrec 0 w
       )
+  showsPrec _ Box1Step = showString "box(1)"
+  showsPrec _ Box2Step = showString "box(2)"
 
-data IntervalNewtonLeaf d
+data RootFindingLeaf d
   = TooSmall d
   | NoSolution String d
   deriving stock Show
 
-showIntervalNewtonTree :: Box -> IntervalNewtonTree Box -> Tree String
+showRootFindingTree :: Box -> RootFindingTree Box -> Tree String
-showIntervalNewtonTree cand (IntervalNewtonLeaf l) = Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ show l) []
+showRootFindingTree cand (RootFindingLeaf l) = Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ show l) []
-showIntervalNewtonTree cand (IntervalNewtonStep s ts)
+showRootFindingTree cand (RootFindingStep s ts)
-  = Node (show cand ++ " abc " ++ showArea (boxArea cand) ++ " " ++ show s) $ map (\ (c,t) -> showIntervalNewtonTree c t) ts
+  = Node (show cand ++ " abc " ++ showArea (boxArea cand) ++ " " ++ show s) $ map (\ (c,t) -> showRootFindingTree c t) ts
 
 boxArea :: Box -> Double
 boxArea ( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ), _, 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
@@ -1466,53 +1358,169 @@ showArea :: Double -> [Char]
 showArea area = "(area " ++ showFFloat (Just 6) area "" ++ ")"
 
 type Box = ( 𝕀ℝ 1, Int, 𝕀ℝ 1 )
+type BoxHistory = [ NE.NonEmpty ( RootFindingStep, Box ) ]
 
-intervalNewtonGS
+-- | Options for the cusp-finding methods in 'findCusps' and 'findCuspsFrom'.
-  :: Preconditioner
+data CuspFindingOptions
-  -> Double
+  = CuspFindingOptions
+  { -- | Minimum width of a box.
+    minWidth :: !Double
+    -- | Given a box and its history, return a round of cusp-finding strategies
+    -- to run, in sequence, on this box.
+  , cuspFindingAlgorithms :: !( Box -> BoxHistory -> NonEmpty RootFindingAlgorithm )
+  }
+
+data RootFindingAlgorithm
+  -- | Gauss–Seidel step with the given preconditioning method.
+  = GaussSeidel GaussSeidelOptions
+  -- | @box(1)@-consistency.
+  | Box1 Box1Options
+  -- | @box(2)@-consistency.
+  | Box2 Box2Options
+  -- | Bisection step.
+  | Bisection BisectionOptions
+
+-- | Options for the interval Gauss–Seidel method.
+data GaussSeidelOptions = GaussSeidelOptions { gsPreconditioner :: Preconditioner }
+
+-- | Options for the @box(1)@-consistency method.
+data Box1Options = Box1Options { box1EpsEq :: !Double }
+
+-- | Options for the @box(2)@-consistency method.
+data Box2Options = Box2Options { box2EpsEq :: !Double, box2LambdaMin :: !Double }
+
+-- | Options for the bisection method.
+data BisectionOptions =
+  BisectionOptions
+  { canHaveSols :: !( ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) -> Box -> Bool )
+  , fallbackBisectionDim :: !( [ ( RootFindingStep, Box ) ] -> BoxHistory -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) -> (Int, String) )
+  }
+
+defaultCuspFindingOptions :: CuspFindingOptions
+defaultCuspFindingOptions =
+  CuspFindingOptions
+    { minWidth
+    , cuspFindingAlgorithms = defaultCuspFindingAlgorithms minWidth 5e-3
+    }
+  where
+    minWidth = 1e-5
+
+defaultCuspFindingAlgorithms :: Double -> Double -> Box -> BoxHistory -> NonEmpty RootFindingAlgorithm
+defaultCuspFindingAlgorithms minWidth narrowAbs box history =
+  case history of
+    lastRoundBoxes : _
+      -- If, in the last round of strategies, we didn't try bisection...
+      | any ( \case { BisectionStep {} -> False; _ -> True } . fst ) lastRoundBoxes
+      , (_, lastRoundFirstBox ) <- NE.last lastRoundBoxes
+      -- ...and the last round didn't sufficiently reduce the size of the box...
+      , not $ box `sufficientlySmallerThan` lastRoundFirstBox
+      -- ...then try bisecting the box.
+      -> Bisection ( defaultBisectionOptions minWidth narrowAbs box ) NE.:| []
+
+    -- Otherwise, do a normal round: try an interval Gauss–Seidel step followed by box(1)-consistency.
+    _ -> GaussSeidel defaultGaussSeidelOptions
+          NE.:| [ Box1 ( Box1Options { box1EpsEq = narrowAbs } ) ]
+  where
+    sufficientlySmallerThan :: Box -> Box -> Bool
+    sufficientlySmallerThan ( t1, _, s1 ) ( t2, _, s2 ) =
+      ( wd t1 - wd t2 > narrowAbs )
+        ||
+      ( wd s1 - wd s2 > narrowAbs )
+    wd :: 𝕀ℝ 1 -> Double
+    wd ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) = hi - lo
+
+defaultGaussSeidelOptions :: GaussSeidelOptions
+defaultGaussSeidelOptions = GaussSeidelOptions { gsPreconditioner = InverseMidJacobian }
+
+defaultBisectionOptions :: Double -> Double -> Box -> BisectionOptions
+defaultBisectionOptions
+  _minWidth
+  _narrowAbs
+  ( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) )
+  , i
+  , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
+  )
+  = BisectionOptions
+      { canHaveSols =
+          \ eqs _box'@( _t', _i', _s' ) ->
+            -- box(0)-consistency
+            let dat' = ( eqs _t' `Seq.index` _i' ) _s'
+            in ee_potential_zero dat' && 𝛿E𝛿sdcdt_potential_zero dat'
+
+            -- box(1)-consistency
+            --not $ null $ makeBox1Consistent eqs _minWidth _narrowAbs _box'
+
+            -- box(2)-consistency
+            --let ( t'', i'', s'') = makeBox2Consistent eqs _minWidth _narrowAbs 0.2 _box'
+            --    dat'' = ( eqs t'' `Seq.index` i'' ) s''
+            --in ee_potential_zero dat'' && 𝛿E𝛿sdcdt_potential_zero dat''
+      , fallbackBisectionDim =
+          \ _roundHist _prevRoundsHist eqs ->
+            let StrokeDatum { ee = D22 _ ( T ee_t ) (T ee_s ) _ _ _
+                            , 𝛿E𝛿sdcdt = D12 _ ( T ( T f_t ) ) ( T ( T f_s ) ) }
+                  = ( eqs t `Seq.index` i ) s
+                wd_t = t_hi - t_lo
+                wd_s = s_hi - s_lo
+                tJWidth = wd_t * normVI3 ee_t f_t
+                sJWidth = wd_s * normVI3 ee_s f_s
+            in if tJWidth >= sJWidth
+                    -- bisect along dimension that maximises  width(x_j) * || J_j || ...
+                    -- ... but don't allow thin boxes
+                    || ( wd_t >= 10 * wd_s )
+                    && not ( wd_s >= 10 * wd_t )
+               then (0, "")
+               else (1, "")
+      }
+
+-- | Find cusps in the envelope for values of the parameters in
+-- \( 0 \leqslant t, s \leqslant 1 \), using interval arithmetic.
+--
+-- See 'findCuspsFrom' for details.
+findCusps
+  :: CuspFindingOptions
   -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-  -> ( [ ( Box, IntervalNewtonTree Box ) ], ( [ Box ], [ Box ] ) )
+  -> ( [ ( Box, RootFindingTree Box ) ], ( [ Box ], [ Box ] ) )
-intervalNewtonGS precondMethod minWidth eqs =
+findCusps opts eqs = foldMap ( findCuspsFrom opts eqs ) initBoxes
-  foldMap
-    ( intervalNewtonGSFrom precondMethod minWidth eqs )
-    initBoxes
   where
     initBoxes =
       [ ( 𝕀 ( ℝ1 zero ) ( ℝ1 one ), i, 𝕀 ( ℝ1 zero ) ( ℝ1 one ) )
       | i <- [ 0 .. length ( eqs ( 𝕀 ( ℝ1 zero ) ( ℝ1 one ) ) ) - 1 ]
       ]
 
--- | Interval Newton method with Gauss–Seidel step for inversion
+-- | Use the following algorithms using interal arithmetic in order
--- of the interval Jacobian.
+-- to find cusps in the envelope:
+--
+--  - interval Newton method with Gauss–Seidel step for inversion
+--    of the interval Jacobian,
+--  - coordinate-wise Newton method (@box(1)@-consistency algorithm),
+--  - @box(2)@-consistency algorithm,
+--  - bisection.
 --
 -- Returns @(tree, (dunno, sols))@ where @tree@ is the full tree (useful for debugging),
 -- @sols@ are boxes that contain a unique solution (and to which Newton's method
 -- will converge starting from anywhere inside the box), and @dunno@ are small
 -- boxes which might or might not contain solutions.
-intervalNewtonGSFrom
+findCuspsFrom
-  :: Preconditioner
+  :: CuspFindingOptions
-  -> Double
   -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
   -> Box
-  -> ( [ ( Box, IntervalNewtonTree Box ) ], ( [ Box ], [ Box ] ) )
+  -> ( [ ( Box, RootFindingTree Box ) ], ( [ Box ], [ Box ] ) )
-intervalNewtonGSFrom precondMethod minWidth eqs initBox =
+findCuspsFrom
-  runWriter $
+  ( CuspFindingOptions
-    map ( initBox , ) <$> go initBox
+    { minWidth
+    , cuspFindingAlgorithms
+    }
+  )
+  eqs initBox = runWriter $ map ( initBox , ) <$> go [] initBox
 
   where
 
-    recur :: ( cand -> Writer ( [ Box ], [ Box ] ) [ IntervalNewtonTree Box ] )
+    go :: BoxHistory
-          -> ( [ ( cand, IntervalNewtonTree Box ) ] -> IntervalNewtonTree Box )
+       -> Box -- box to work on
-          -> [ cand ]
-          -> Writer ( [Box], [Box] ) [ IntervalNewtonTree Box ]
-    recur r f cands = do
-      rest <- traverse ( \ c -> do { trees <- r c; return [ (c, tree) | tree <- trees ] } ) cands
-      return [ f $ concat rest ]
-
-    go :: Box -- box to work on
        -> Writer ( [ Box ], [ Box ] )
-                 [ IntervalNewtonTree Box ]
+                 [ RootFindingTree Box ]
-    go cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) )
+    go history
+       cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) )
             , i
             , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
             )
@@ -1520,212 +1528,251 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox =
       | t_hi - t_lo < minWidth && s_hi - s_lo < minWidth
       = do let res = TooSmall cand
            tell ( [ cand ], [] )
-           return [ IntervalNewtonLeaf res ]
+           return [ RootFindingLeaf res ]
       | -- Check the envelope equation interval contains zero.
         ee_potential_zero sd
         -- Check the 𝛿E𝛿sdcdt interval contains zero.
       , 𝛿E𝛿sdcdt_potential_zero sd
-      = continue cand
+      = doStrategies history ( cuspFindingAlgorithms cand history ) cand
       -- Box doesn't contain a solution: discard it.
       | otherwise
-      = return [ IntervalNewtonLeaf $ NoSolution ( if ee_potential_zero sd then "dc/dt" else "ee" ) cand ]
+      , let whyNoSol = if ee_potential_zero sd then "dc/dt" else "ee"
+      = return [ RootFindingLeaf $ NoSolution whyNoSol cand ]
       where
         sd = ( eqs t `Seq.index` i ) s
 
-    continue :: Box
+    -- Run a round of cusp finding strategies, then recur.
-             -> Writer ( [ Box ], [ Box ] )
+    doStrategies
-                       [ IntervalNewtonTree Box ]
+      :: BoxHistory
-    continue
+      -> NonEmpty RootFindingAlgorithm
-      cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) )
+      -> Box
-           , i
+      -> Writer ( [ Box ], [ Box ] )
-           , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
+                [ RootFindingTree Box ]
-           )
+    doStrategies hist0 = do_strats []
-      | StrokeDatum { ee = D22 _ ( T ee_t ) (T ee_s ) _ _ _
-                    , 𝛿E𝛿sdcdt = D12 _ ( T ( T f_t ) ) ( T ( T f_s ) ) }
-        <- ( eqs t `Seq.index` i ) s
-      , StrokeDatum { ee = D22 _ee_mid _ _ _ _ _
-                    , 𝛿E𝛿sdcdt = D12 ( T f_mid ) ( T ( T _f_t_mid ) ) ( T ( T _f_s_mid ) ) }
-        <- ( eqs i_t_mid `Seq.index` i ) i_s_mid
-      = let -- Interval Newton method: take one Gauss–Seidel step
-            -- for the equation f'(x) v = - f(x_mid),
-            -- where f = 𝛿E/𝛿s * dc/dt
-            !( a, b ) = precondition precondMethod
-                          ( midI f_t, midI f_s )
-                          ( f_t, f_s ) ( neg f_mid )
-            --(a, b)
-            --  | 𝕀 (ℝ1 ee_lo) (ℝ1 ee_hi) <- ee_mid
-            --  , 𝕀 (ℝ1 ee_t_lo) (ℝ1 ee_t_hi) <- ee_t
-            --  , 𝕀 (ℝ1 ee_s_lo) (ℝ1 ee_s_hi) <- ee_s
-            --  , 𝕀 (ℝ2 fx_lo fy_lo) (ℝ2 fx_hi fy_hi) <- f_mid
-            --  , 𝕀 (ℝ2 f_tx_lo f_ty_lo) (ℝ2 f_tx_hi f_ty_hi) <- f_t
-            --  , 𝕀 (ℝ2 f_sx_lo f_sy_lo) (ℝ2 f_sx_hi f_sy_hi) <- f_s
-            --  = ( ( 𝕀 (ℝ2 f_tx_lo ee_t_lo) (ℝ2 f_tx_hi ee_t_hi)
-            --      , 𝕀 (ℝ2 f_sx_lo ee_s_lo) (ℝ2 f_sx_hi ee_s_hi)
-            --      )
-            --    , neg $ 𝕀 (ℝ2 fx_lo ee_lo) (ℝ2 fx_hi ee_hi)
-            --    )
-
-            !gsGuesses = gaussSeidel a b
-                           ( coerce ( (-) @( 𝕀 Double ) ) t i_t_mid
-                           , coerce ( (-) @( 𝕀 Double ) ) s i_s_mid )
-        in if any ( smaller . fst ) gsGuesses
-           then
-             -- If the Gauss–Seidel step was a contraction, then the box
-             -- contains a unique solution (by the Banach fixed point theorem).
-             --
-             -- These boxes can thus be directly added to the solution set:
-             -- Newton's method is guaranteed to converge to the unique solution.
-             let !(done, todo) = bimap ( map ( mkGuess . fst ) ) ( map ( mkGuess . fst ) )
-                               $ partition snd gsGuesses
-             in do tell ([], done)
-                   case todo of
-                     [] -> return [ IntervalNewtonLeaf $ NoSolution "GaussSeidel" cand ]
-                     _ -> recur go ( IntervalNewtonStep ( IntervalNewtonContraction done ) )
-                            todo
-           else
-              -- Gauss–Seidel failed to shrink the boxes, so bisect instead.
-              -- We have to pick along which dimension to bisect:
-              --  - if bisecting along a particular dimension discards one of
-              --    the boxes, do that;
-              --  - otherwise, bisect along the dimension j that maximises
-              --    width(x_j) * || J_j ||.
-             let l_t = 𝕀 ( ℝ1 t_lo  ) ( ℝ1 t_mid )
-                 r_t = 𝕀 ( ℝ1 t_mid ) ( ℝ1 t_hi )
-                 d_s = 𝕀 ( ℝ1 s_lo  ) ( ℝ1 s_mid )
-                 u_s = 𝕀 ( ℝ1 s_mid ) ( ℝ1 s_hi  )
-                 l = ( l_t, i, s )
-                 r = ( r_t, i, s )
-                 d = ( t, i, d_s )
-                 u = ( t, i, u_s )
-                 l_dat = ( eqs l_t `Seq.index` i ) s
-                 r_dat = ( eqs r_t `Seq.index` i ) s
-                 d_dat = ( eqs t   `Seq.index` i ) d_s
-                 u_dat = ( eqs t   `Seq.index` i ) u_s
-                 l_skip =
-                     not ( ee_potential_zero       l_dat )
-                  || not ( 𝛿E𝛿sdcdt_potential_zero l_dat )
-                 r_skip =
-                      not ( ee_potential_zero       r_dat )
-                   || not ( 𝛿E𝛿sdcdt_potential_zero r_dat )
-                 d_skip =
-                      not ( ee_potential_zero       d_dat )
-                   || not ( 𝛿E𝛿sdcdt_potential_zero d_dat )
-                 u_skip =
-                      not ( ee_potential_zero       u_dat )
-                   || not ( 𝛿E𝛿sdcdt_potential_zero u_dat )
-                 tJWidth = ( t_hi - t_lo ) * normVI3 ee_t f_t
-                 sJWidth = ( s_hi - s_lo ) * normVI3 ee_s f_s
-                 ( bisGuesses, whatBis )
-                   | l_skip && r_skip
-                   = ( [], ( "lr", t_mid ) )
-                   | d_skip && u_skip
-                   = ( [], ( "du", s_mid ) )
-                   | l_skip
-                   = ( [ r ], ( "r", t_mid ) )
-                   | r_skip
-                   = ( [ l ], ( "l", t_mid ) )
-                   | d_skip
-                   = ( [ u ], ( "u", s_mid ) )
-                   | u_skip
-                   = ( [ d ], ( "d", s_mid ) )
-                   | tJWidth >= sJWidth
-                   -- ... but don't allow thin boxes.
-                   || ( t_hi - t_lo >= 10 * ( s_hi - s_lo ) )
-                   && not ( s_hi - s_lo >= 10 * ( t_hi - t_lo ) )
-                   = ( [ l, r ], ( "t", t_mid ) )
-                   | otherwise
-                   = ( [ d, u ], ( "s", s_mid ) )
-             in recur go ( IntervalNewtonStep ( IntervalNewtonBisection whatBis ) )
-                  ( doStrategy =<< bisGuesses )
       where
-        t_mid = 0.5 * ( t_lo + t_hi )
+        do_strats :: [ ( RootFindingStep, Box ) ]
+                  -> NE.NonEmpty RootFindingAlgorithm
+                  -> Box
+                  -> Writer ( [ Box ], [ Box ] )
+                            [ RootFindingTree Box ]
+        do_strats stratsHist ( algo NE.:| algos ) box = do
+          -- Run one strategy in the round.
+          ( step, boxes ) <- doStrategy stratsHist hist0 eqs minWidth algo box
+          case algos of
+            -- If there are other algorithms to run in this round, run them next.
+            nextAlgo : otherAlgos ->
+              let stratsHist' = ( step, box ) : stratsHist
+              in recur step ( do_strats stratsHist' ( nextAlgo NE.:| otherAlgos ) ) boxes
+            -- Otherwise, we have done one full round of strategies.
+            -- Recur back to the top (calling 'go').
+            [] ->
+              let stratsHist' = ( step, box ) NE.:| stratsHist
+                  history = stratsHist' : hist0
+              in recur step ( go history ) boxes
+
+    recur :: RootFindingStep
+          -> ( Box -> Writer ( [ Box ], [ Box ] ) [ RootFindingTree Box ] )
+          -> [ Box ]
+          -> Writer ( [Box], [Box] ) [ RootFindingTree Box ]
+    recur step r cands = do
+      rest <- traverse ( \ c -> do { trees <- r c; return [ (c, tree) | tree <- trees ] } ) cands
+      return [ RootFindingStep step $ concat rest ]
+
+
+-- | Execute a cusp-finding strategy, replacing the input box with
+-- (hopefully smaller) output boxes.
+doStrategy :: [ ( RootFindingStep, Box ) ]
+           -> BoxHistory
+           -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
+           -> Double
+           -> RootFindingAlgorithm
+           -> Box
+           -> Writer ( [ Box ], [ Box ] )
+                     ( RootFindingStep, [ Box ] )
+doStrategy roundHistory previousRoundsHistory eqs minWidth algo box =
+  case algo of
+    GaussSeidel ( GaussSeidelOptions { gsPreconditioner } ) -> do
+      boxes <- intervalGaussSeidel eqs gsPreconditioner box
+      return ( GaussSeidelStep, boxes )
+    Box1 ( Box1Options { box1EpsEq } ) ->
+      return ( Box1Step, makeBox1Consistent eqs minWidth box1EpsEq box )
+    Box2 ( Box2Options { box2LambdaMin, box2EpsEq } ) ->
+      return ( Box2Step, [ makeBox2Consistent eqs minWidth box2EpsEq box2LambdaMin box ] )
+    Bisection ( BisectionOptions { canHaveSols, fallbackBisectionDim } ) -> do
+      let ( boxes, whatBis ) = bisect ( canHaveSols eqs ) ( fallbackBisectionDim roundHistory previousRoundsHistory eqs ) box
+      return ( BisectionStep whatBis, boxes )
+
+-- | Bisect the given box.
+--
+-- (The difficult part lies in determining along which dimension to bisect.)
+bisect :: ( Box -> Bool )
+          -- ^ how to check whether a box contains solutions
+       -> ( Int, String )
+           -- ^ fall-back choice of dimension (and "why" we chose it)
+       -> Box
+       -> ( [ Box ], ( String, Double ) )
+bisect canHaveSols ( dim, why )
+  ( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) )
+  , i
+  , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
+  )
+  -- We try to bisect along a dimension which eliminates zeros from one of the
+  -- sub-regions.
+  --
+  -- If this fails, we fall-back to the provided dimension choice.
+  = let t_mid = 0.5 * ( t_lo + t_hi )
         s_mid = 0.5 * ( s_lo + s_hi )
-        i_t_mid = singleton ( ℝ1 t_mid )
+        l_t = 𝕀 ( ℝ1 t_lo  ) ( ℝ1 t_mid )
-        i_s_mid = singleton ( ℝ1 s_mid )
+        r_t = 𝕀 ( ℝ1 t_mid ) ( ℝ1 t_hi  )
-        mkGuess ( t0, s0 ) = ( coerce ( (+) @( 𝕀 Double ) ) t0 i_t_mid
+        d_s = 𝕀 ( ℝ1 s_lo  ) ( ℝ1 s_mid )
-                             , i
+        u_s = 𝕀 ( ℝ1 s_mid ) ( ℝ1 s_hi  )
-                             , coerce ( (+) @( 𝕀 Double ) ) s0 i_s_mid )
+        l = ( l_t, i, s   )
-        smaller ( 𝕀 ( ℝ1 t0_lo ) ( ℝ1 t0_hi ), 𝕀 ( ℝ1 s0_lo ) ( ℝ1 s0_hi ) )
+        r = ( r_t, i, s   )
-          =  ( t0_lo + t_mid ) > t_lo + 0.25 * minWidth
+        d = ( t  , i, d_s )
-          || ( t0_hi + t_mid ) < t_hi - 0.25 * minWidth
+        u = ( t  , i, u_s )
-          || ( s0_lo + s_mid ) > s_lo + 0.25 * minWidth
+        l_skip = not $ canHaveSols l
-          || ( s0_hi + s_mid ) < s_hi - 0.25 * minWidth
+        r_skip = not $ canHaveSols r
-        neg ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) )
+        d_skip = not $ canHaveSols d
-          = let !( 𝕀 x'_lo x'_hi ) = negate $ 𝕀 x_lo x_hi
+        u_skip = not $ canHaveSols u
-                !( 𝕀 y'_lo y'_hi ) = negate $ 𝕀 y_lo y_hi
-            in 𝕀 ( ℝ2 x'_lo y'_lo ) ( ℝ2 x'_hi y'_hi )
 
-    -- Attempting to implement Algorithm 6 "Heuristic to apply bound-consistency"
+        ( bisGuesses, whatBis )
-    -- from the paper
+          | l_skip && r_skip
-    -- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
+          = ( [], ( "lr", t_mid ) )
-    doStrategy box@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ), _, 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
+          | d_skip && u_skip
-      | eps_box1 > eps_box2
+          = ( [], ( "du", s_mid ) )
-      = if | w_t > ebs_box1 || w_s > ebs_box1
+          | l_skip
-           -> makeBox1Consistent box
+          = ( [ r ], ( "r", t_mid ) )
-           | w_t > ebs_box2 || w_s > ebs_box2
+          | r_skip
-           -> [ makeBox2Consistent box ]
+          = ( [ l ], ( "l", t_mid ) )
-           | otherwise
+          | d_skip
-           -> [ box ]
+          = ( [ u ], ( "u", s_mid ) )
-      | otherwise
+          | u_skip
-      = if | w_t > ebs_box2 || w_s > ebs_box2
+          = ( [ d ], ( "d", s_mid ) )
-           -> [ makeBox2Consistent box ]
+          | otherwise
-           | w_t > ebs_box1 || w_s > ebs_box1
+          = let why' = case why of
-           -> makeBox1Consistent box
+                  "" -> ""
-           | otherwise
+                  _ -> " (" ++ why ++ ")"
-           -> [ box ]
+            in case dim of
-      where
+                 0 -> ( [ l, r ], ( "t" ++ why', t_mid ) )
-        eps_box1 = 0.4
+                 1 -> ( [ d, u ], ( "s" ++ why', s_mid ) )
-        eps_box2 = 0.1
+                 _ -> error "bisect: fall-back bisection dimension should be either 0 or 1"
-        w_t = t_hi - t_lo
+         in ( bisGuesses, whatBis )
-        w_s = s_hi - s_lo
 
-    -- An implementation of "bound-consistency" from the paper
+-- | Interval Newton method with Gauss–Seidel step.
-    -- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
+intervalGaussSeidel
-    makeBox2Consistent :: Box -> Box
+  :: ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-    makeBox2Consistent x = ( `State.evalState` False ) $ doLoop 0.25 x
+     -- ^ equations
-      where
+  -> Preconditioner
-        doLoop :: Double -> Box -> State Bool Box
+     -- ^ preconditioner to use for the Gauss–Seidel step
-        doLoop lambda x = do
+  -> Box
-          x'  <- boundConsistency get_t set_t lambda x
+  -> Writer ( [ Box ], [ Box ] ) [ Box ]
-          x'' <- boundConsistency get_s set_s lambda x'
+intervalGaussSeidel
-          modified <- State.get
+  eqs
-          let lambda' = if modified then lambda else 0.5 * lambda
+  precondMethod
-          if lambda' < 0.001
+  ( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) )
-          then return x''
+  , i
-          else do { State.put False ; doLoop lambda' x'' }
+  , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
+  )
+  | StrokeDatum { 𝛿E𝛿sdcdt = D12 _ ( T ( T f_t ) ) ( T ( T f_s ) ) }
+    <- ( eqs t `Seq.index` i ) s
+  , StrokeDatum { ee = D22 _ee_mid _ _ _ _ _
+                , 𝛿E𝛿sdcdt = D12 ( T f_mid ) ( T ( T _f_t_mid ) ) ( T ( T _f_s_mid ) ) }
+    <- ( eqs i_t_mid `Seq.index` i ) i_s_mid
+  = let -- Interval Newton method: take one Gauss–Seidel step
+        -- for the equation f'(x) v = - f(x_mid),
+        -- where f = 𝛿E/𝛿s * dc/dt
+        ( a, b ) = precondition precondMethod
+                      ( midI f_t, midI f_s )
+                      ( f_t, f_s ) ( neg f_mid )
+        --(a, b)
+        --  | 𝕀 (ℝ1 ee_lo) (ℝ1 ee_hi) <- ee_mid
+        --  , 𝕀 (ℝ1 ee_t_lo) (ℝ1 ee_t_hi) <- ee_t
+        --  , 𝕀 (ℝ1 ee_s_lo) (ℝ1 ee_s_hi) <- ee_s
+        --  , 𝕀 (ℝ2 fx_lo fy_lo) (ℝ2 fx_hi fy_hi) <- f_mid
+        --  , 𝕀 (ℝ2 f_tx_lo f_ty_lo) (ℝ2 f_tx_hi f_ty_hi) <- f_t
+        --  , 𝕀 (ℝ2 f_sx_lo f_sy_lo) (ℝ2 f_sx_hi f_sy_hi) <- f_s
+        --  = ( ( 𝕀 (ℝ2 f_tx_lo ee_t_lo) (ℝ2 f_tx_hi ee_t_hi)
+        --      , 𝕀 (ℝ2 f_sx_lo ee_s_lo) (ℝ2 f_sx_hi ee_s_hi)
+        --      )
+        --    , neg $ 𝕀 (ℝ2 fx_lo ee_lo) (ℝ2 fx_hi ee_hi)
+        --    )
 
-        boundConsistency :: ( Box -> 𝕀 Double )
+        gsGuesses = gaussSeidel a b
-                         -> ( 𝕀 Double -> Box -> Box )
+                      ( coerce ( (-) @( 𝕀 Double ) ) t i_t_mid
-                         -> Double -> Box -> State Bool Box
+                      , coerce ( (-) @( 𝕀 Double ) ) s i_s_mid )
-        boundConsistency getter setter lambda box = do
+    in
-          let x@( 𝕀 x_inf x_sup ) = getter box
+         -- If the Gauss–Seidel step was a contraction, then the box
-              c1 = ( 1 - lambda ) * x_inf + lambda * x_sup
+         -- contains a unique solution (by the Banach fixed point theorem).
-              c2 = lambda * x_inf + ( 1 - lambda ) * x_sup
+         --
-              x'_inf =
+         -- These boxes can thus be directly added to the solution set:
-                case makeBox1Consistent ( setter ( 𝕀 x_inf c1 ) box ) of
+         -- Newton's method is guaranteed to converge to the unique solution.
-                  []  -> c1
+         let !(done, todo) = bimap ( map ( mkGuess . fst ) ) ( map ( mkGuess . fst ) )
-                  x's -> minimum $ map ( inf . getter ) x's
+                           $ partition snd gsGuesses
-              x'_sup =
+         in do tell ([], done)
-                case makeBox1Consistent ( setter ( 𝕀 c2 x_sup ) box ) of
+               return todo
-                  []  -> c2
+  where
-                  x's -> maximum $ map ( sup . getter ) x's
+    t_mid = 0.5 * ( t_lo + t_hi )
-              x' = 𝕀 x'_inf x'_sup
+    s_mid = 0.5 * ( s_lo + s_hi )
-          when ( width x - width x' >= eps_eq ) $
+    i_t_mid = singleton ( ℝ1 t_mid )
-            State.put True
+    i_s_mid = singleton ( ℝ1 s_mid )
-          return $ setter x' box
+    mkGuess ( t0, s0 ) = ( coerce ( (+) @( 𝕀 Double ) ) t0 i_t_mid
+                         , i
+                         , coerce ( (+) @( 𝕀 Double ) ) s0 i_s_mid )
+    neg ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) )
+      = let !( 𝕀 x'_lo x'_hi ) = negate $ 𝕀 x_lo x_hi
+            !( 𝕀 y'_lo y'_hi ) = negate $ 𝕀 y_lo y_hi
+        in 𝕀 ( ℝ2 x'_lo y'_lo ) ( ℝ2 x'_hi y'_hi )
 
-    -- An implementation of "bc_enforce" from the paper
+-- | An implementation of "bc_enforce" from the paper
-    -- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
+-- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
-    --
+--
-    -- See also
+-- See also
-    -- "Presentation of a highly tuned multithreaded interval solver for underdetermined and well-determined nonlinear systems"
+-- "Presentation of a highly tuned multithreaded interval solver for underdetermined and well-determined nonlinear systems"
-    makeBox1Consistent :: Box -> [ Box ]
+makeBox1Consistent
-    makeBox1Consistent x =
+  :: ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-      ( `State.evalState` False ) $
+  -> Double -> Double
-        pipeFunctionsWhileTrue (allNarrowingOperators eps_eq eps_bis eqs) x
+  -> Box -> [ Box ]
-    eps_eq = 1e-5 -- TODO: this is an additive constant,
+makeBox1Consistent eqs minWidth epsEq x =
-                  -- but I think we also want a multiplicative constant?
+  ( `State.evalState` False ) $
-    eps_bis = minWidth
+    pipeFunctionsWhileTrue (allNarrowingOperators epsEq minWidth eqs) x
+
+-- | An implementation of "bound-consistency" from the paper
+-- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
+makeBox2Consistent
+  :: ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
+  -> Double -> Double -> Double
+  -> Box -> Box
+makeBox2Consistent eqs minWidth epsEq lambdaMin x0 = ( `State.evalState` False ) $ doLoop 0.25 x0
+  where
+    doLoop :: Double -> Box -> State Bool Box
+    doLoop lambda x = do
+      x'  <- boundConsistency get_t set_t lambda x
+      x'' <- boundConsistency get_s set_s lambda x'
+      modified <- State.get
+      let lambda' = if modified then lambda else 0.5 * lambda
+      if lambda' < lambdaMin
+      then return x''
+      else do { State.put False ; doLoop lambda' x'' }
+
+    boundConsistency :: ( Box -> 𝕀 Double )
+                     -> ( 𝕀 Double -> Box -> Box )
+                     -> Double -> Box -> State Bool Box
+    boundConsistency getter setter lambda box = do
+      let x@( 𝕀 x_inf x_sup ) = getter box
+          c1 = ( 1 - lambda ) * x_inf + lambda * x_sup
+          c2 = lambda * x_inf + ( 1 - lambda ) * x_sup
+          x'_inf =
+            case makeBox1Consistent eqs minWidth epsEq ( setter ( 𝕀 x_inf c1 ) box ) of
+              []  -> c1
+              x's -> minimum $ map ( inf . getter ) x's
+          x'_sup =
+            case makeBox1Consistent eqs minWidth epsEq ( setter ( 𝕀 c2 x_sup ) box ) of
+              []  -> c2
+              x's -> maximum $ map ( sup . getter ) x's
+          x' = 𝕀 x'_inf x'_sup
+      when ( width x - width x' >= epsEq ) $
+        State.put True
+      return $ setter x' box
 
 ee_potential_zero :: StrokeDatum 3 𝕀 -> Bool
 ee_potential_zero dat =
@@ -1736,7 +1783,6 @@ ee_potential_zero dat =
      cmpℝ2 (<=) ( inf $ unT $ _D12_v $ 𝛿E𝛿sdcdt dat ) ( ℝ2 0 0 )
   && cmpℝ2 (>=) ( sup $ unT $ _D12_v $ 𝛿E𝛿sdcdt dat ) ( ℝ2 0 0 )
 
-
 zero, one :: Double
 zero = 1e-6
 one = 1 - zero
@@ -1841,7 +1887,7 @@ leftNarrow eps_equal eps_bisection ff' = left_narrow
                 x' <- x's
                 if sup x' - inf x' < eps_bisection
                 then return x'
-                else left_narrow =<< bisect x'
+                else left_narrow =<< bisectI x'
 
 -- TODO: de-duplicate with 'leftNarrow'?
 rightNarrow :: Double
@@ -1875,10 +1921,10 @@ rightNarrow eps_equal eps_bisection ff' = right_narrow
                 x' <- x's
                 if sup x' - inf x' < eps_bisection
                 then return x'
-                else right_narrow =<< bisect x'
+                else right_narrow =<< bisectI x'
 
-bisect :: 𝕀 Double -> [ 𝕀 Double ]
+bisectI :: 𝕀 Double -> [ 𝕀 Double ]
-bisect x@( 𝕀 x_inf x_sup )
+bisectI x@( 𝕀 x_inf x_sup )
   | x_inf == x_sup
   = [ x ]
   | otherwise
@@ -1953,7 +1999,10 @@ allNarrowingOperators eps_eq eps_bis eqs =
               } = 𝛿E𝛿sdcdt $ ( eqs t `Seq.index` i ) s
       in ( 𝕀 y_inf y_sup, 𝕀 y'_inf y'_sup )
 
+get_t, get_s :: Box -> 𝕀 Double
 get_t ( 𝕀 ( ℝ1 t_inf ) ( ℝ1 t_sup ), _, _ ) = 𝕀 t_inf t_sup
+
+set_t, set_s :: 𝕀 Double -> Box -> Box
 set_t ( 𝕀 t_inf t_sup ) ( _, i, s ) = ( 𝕀 ( ℝ1 t_inf ) ( ℝ1 t_sup ), i, s )
 get_s ( _ , _, 𝕀 ( ℝ1 s_inf ) ( ℝ1 s_sup ) ) = 𝕀 s_inf s_sup
-set_s ( 𝕀 s_inf s_sup ) ( t, i, _ ) = ( t, i, 𝕀 ( ℝ1 s_inf ) ( ℝ1 s_sup ) )
+set_s ( 𝕀 s_inf s_sup ) ( t, i, _ ) = ( t, i, 𝕀 ( ℝ1 s_inf ) ( ℝ1 s_sup ) )

src/app/MetaBrush/Application.hs

@@ -78,7 +78,10 @@ import Math.Bezier.Cubic.Fit
 import Math.Bezier.Spline
   ( Spline(..), Curves(..), Curve(..), NextPoint(..) )
 import Math.Bezier.Stroke
-  ( RootSolvingAlgorithm(..), invalidateCache )
+  ( RootSolvingAlgorithm(..)
+  , CuspFindingOptions(..), Preconditioner(..)
+  , invalidateCache
+  )
 import Math.Linear
   ( ℝ(..) )
 import MetaBrush.Action
@@ -212,6 +215,8 @@ runApplication application = do
                             --HalleyM2
                             NewtonRaphson
                               { maxIters = 20, precision = 8 }
+  cuspFindingOptionsTVar <- STM.newTVarIO @( Maybe CuspFindingOptions ) $
+                              Just defaultCuspFindingOptions
 
   -- Put all these stateful variables in a record for conciseness.
   let
@@ -337,6 +342,7 @@ runApplication application = do
                 debug          <- STM.readTVar debugTVar
                 fitParameters  <- STM.readTVar fitParametersTVar
                 rootsAlgo      <- STM.readTVar rootsAlgoTVar
+                mbCuspOptions  <- STM.readTVar cuspFindingOptionsTVar
                 STM.writeTVar recomputeStrokesTVar False
                 let
                   addRulers :: ( ( Int32, Int32 ) -> Cairo.Render () ) -> ( ( Int32, Int32 ) -> Cairo.Render () ) 
@@ -348,7 +354,9 @@ runApplication application = do
                       doc
                 pure
                   ( addRulers <$> getDocumentRender
-                      colours rootsAlgo fitParameters mode debug
+                      colours
+                      rootsAlgo mbCuspOptions fitParameters
+                      mode debug
                       modifiers mbMousePos mbHoldAction mbPartialPath
                       doc
                   )

src/app/MetaBrush/Render/Document.hs

@@ -74,7 +74,7 @@ import Math.Bezier.Spline
 import Math.Bezier.Stroke
   ( Cusp(..), CachedStroke(..), invalidateCache
   , computeStrokeOutline
-  , RootSolvingAlgorithm
+  , RootSolvingAlgorithm, CuspFindingOptions
   )
 import Math.Linear
   ( ℝ(..), T(..) )
@@ -148,12 +148,16 @@ blankRender :: Colours -> Cairo.Render ()
 blankRender _ = pure ()
 
 getDocumentRender
-  :: Colours -> RootSolvingAlgorithm -> FitParameters -> Mode -> Bool
+  :: Colours
+  -> RootSolvingAlgorithm -> Maybe CuspFindingOptions -> FitParameters
+  -> Mode -> Bool
   -> Set Modifier -> Maybe ( ℝ 2 ) -> Maybe HoldAction -> Maybe PartialPath
   -> Document
   -> ST RealWorld ( ( Int32, Int32 ) -> Cairo.Render () )
 getDocumentRender
-  cols rootAlgo fitParams mode debug
+  cols
+  rootAlgo mbCuspOptions fitParams
+  mode debug
   modifiers mbMousePos mbHoldEvent mbPartialPath
   doc@( Document { viewportCenter = ℝ2 cx cy, zoomFactor, documentContent = content } )
     = do
@@ -215,7 +219,10 @@ getDocumentRender
               -> previewStroke :<| foldMap visibleStrokes ( strokes content )
           _ -> foldMap visibleStrokes ( strokes content )
 
-      strokesRenderData <- traverseMaybe ( sequenceA . strokeRenderData rootAlgo fitParams ) modifiedStrokes
+      strokesRenderData <-
+        traverseMaybe
+         ( sequenceA . strokeRenderData rootAlgo mbCuspOptions fitParams )
+         modifiedStrokes
 
       let
         renderSelectionRect :: Cairo.Render ()
@@ -276,8 +283,13 @@ instance NFData StrokeRenderData where
 --    - the computed outline (using fitting algorithm),
 --    - the brush shape function.
 -- - Otherwise, this consists of the underlying spline path only.
-strokeRenderData :: RootSolvingAlgorithm -> FitParameters -> Stroke -> Maybe ( ST RealWorld StrokeRenderData )
+strokeRenderData
-strokeRenderData rootAlgo fitParams
+  :: RootSolvingAlgorithm
+  -> Maybe CuspFindingOptions
+  -> FitParameters
+  -> Stroke
+  -> Maybe ( ST RealWorld StrokeRenderData )
+strokeRenderData rootAlgo mbCuspOptions fitParams
   ( Stroke
     { strokeSpline = spline :: StrokeSpline clo ( Record pointFields )
     , strokeBrush  = ( strokeBrush :: Maybe ( Brush brushFields ) )
@@ -302,7 +314,7 @@ strokeRenderData rootAlgo fitParams
 
                   -- Compute the outline using the brush function.
                     ( outline, fitPts, cusps ) <-
-                      computeStrokeOutline @clo rootAlgo fitParams
+                      computeStrokeOutline @clo rootAlgo mbCuspOptions fitParams
                         ( toUsedParams . brushParams ) embedUsedParams brushFn
                         spline
                     pure $