diff --git a/brush-strokes/brush-strokes.cabal b/brush-strokes/brush-strokes.cabal index 6362515..3fe1b13 100644 --- a/brush-strokes/brush-strokes.cabal +++ b/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 diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs index 1fb8cc4..760fcfb 100644 --- a/brush-strokes/src/cusps/bench/Main.hs +++ b/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 - , testBrushStroke :: BrushStroke - , testAlgorithmParams :: CuspAlgorithmParams - , testStartBoxes :: [ Box ] + { testName :: String + , testBrushStroke :: BrushStroke + , testCuspOptions :: CuspFindingOptions + , 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 str k0k1 rot startBoxes = +ellipseTestCase :: CuspFindingOptions -> String -> ( Double, Double ) -> Double -> [ Box ] -> TestCase +ellipseTestCase opts str k0k1 rot startBoxes = TestCase { testName = "ellipse (" ++ str ++ ")" , testBrushStroke = ellipseBrushStroke k0k1 rot - , testAlgorithmParams = - CuspAlgorithmParams - { preconditioning = InverseMidJacobian - , maxWidth = 1e-7 - } + , testCuspOptions = opts , testStartBoxes = startBoxes } @@ -442,11 +510,7 @@ trickyCusp2TestCase = TestCase { testName = "trickyCusp2" , testBrushStroke = trickyCusp2BrushStroke - , testAlgorithmParams = - CuspAlgorithmParams - { preconditioning = InverseMidJacobian - , maxWidth = 1e-7 - } + , testCuspOptions = defaultCuspFindingOptions , 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 ] ) ) -computeCusps params eqs startBoxes = - foldMap - ( intervalNewtonGSFrom ( preconditioning params ) ( maxWidth params ) eqs ) - startBoxes + -> ( [ ( Box, RootFindingTree Box ) ], ( [ Box ], [ Box ] ) ) +computeCusps cuspOpts eqs startBoxes = + foldMap ( findCuspsFrom cuspOpts 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 -showTrees box = putStrLn $ unlines $ map Data.Tree.View.showTree $ map ( uncurry showIntervalNewtonTree ) $ fst $ intervalNewtonGSFrom NoPreconditioning 1e-7 fI box +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 $ findCuspsFrom NoPreconditioning 1e-7 fI box sol_t = 0.5486100729150693796677845183880669025324233347060776339185 :: Double sol_s = 0.5480950141859386853197594577293968665598143630958601978245 :: Double diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke.hs b/brush-strokes/src/lib/Math/Bezier/Stroke.hs index 579c655..d6b4373 100644 --- a/brush-strokes/src/lib/Math/Bezier/Stroke.hs +++ b/brush-strokes/src/lib/Math/Bezier/Stroke.hs @@ -19,13 +19,14 @@ module Math.Bezier.Stroke , brushStrokeData, pathAndUsedParams -- ** Cusp finding - , Preconditioner(..) - , IntervalNewtonTree(..), showIntervalNewtonTree - , IntervalNewtonStep(..) - , IntervalNewtonLeaf(..) + , RootFindingTree(..), showRootFindingTree + , RootFindingStep(..) + , RootFindingLeaf(..) , 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 - InverseMidJacobian - 1e-7 - curvesI + case mbCuspOptions of + Just cuspOptions -> + findCusps cuspOptions 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 ) - -{- -gaussSeidel2 :: Int - -> 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 ) --} +-- TODO: try implementing the complete interval union Gauss–Seidel algorithm. +-- See "Algorithm 2" in +-- "Using interval unions to solve linear systems of equations with uncertainties" -- | 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. -data IntervalNewtonTree d - = IntervalNewtonLeaf (IntervalNewtonLeaf d) - | IntervalNewtonStep (IntervalNewtonStep d) [(d, IntervalNewtonTree d)] +-- | A tree recording the steps taken when doing cusp finding. +data RootFindingTree d + = RootFindingLeaf (RootFindingLeaf d) + | RootFindingStep RootFindingStep [(d, RootFindingTree d)] -data IntervalNewtonStep d - = IntervalNewtonContraction [d] - | IntervalNewtonBisection (String, Double) +-- | A description of a step taken in cusp-finding. +data RootFindingStep + = GaussSeidelStep + | BisectionStep (String, Double) + | Box1Step + | Box2Step -instance Show d => Show (IntervalNewtonStep d) where - showsPrec _ ( IntervalNewtonContraction v ) - = showString "N " - . showsPrec 0 v - showsPrec _ ( IntervalNewtonBisection (coord, w) ) - = showString "B " +instance Show RootFindingStep where + showsPrec _ GaussSeidelStep = showString "GS" + showsPrec _ ( BisectionStep (coord, w) ) + = showString "bis " . 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 -showIntervalNewtonTree cand (IntervalNewtonLeaf l) = Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ show l) [] -showIntervalNewtonTree cand (IntervalNewtonStep s ts) - = Node (show cand ++ " abc " ++ showArea (boxArea cand) ++ " " ++ show s) $ map (\ (c,t) -> showIntervalNewtonTree c t) ts +showRootFindingTree :: Box -> RootFindingTree Box -> Tree String +showRootFindingTree cand (RootFindingLeaf l) = Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ show l) [] +showRootFindingTree cand (RootFindingStep s 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 - :: Preconditioner - -> Double +-- | Options for the cusp-finding methods in 'findCusps' and 'findCuspsFrom'. +data CuspFindingOptions + = 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 ] ) ) -intervalNewtonGS precondMethod minWidth eqs = - foldMap - ( intervalNewtonGSFrom precondMethod minWidth eqs ) - initBoxes + -> ( [ ( Box, RootFindingTree Box ) ], ( [ Box ], [ Box ] ) ) +findCusps opts eqs = foldMap ( findCuspsFrom opts 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 --- of the interval Jacobian. +-- | Use the following algorithms using interal arithmetic in order +-- 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 - :: Preconditioner - -> Double +findCuspsFrom + :: CuspFindingOptions -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) -> Box - -> ( [ ( Box, IntervalNewtonTree Box ) ], ( [ Box ], [ Box ] ) ) -intervalNewtonGSFrom precondMethod minWidth eqs initBox = - runWriter $ - map ( initBox , ) <$> go initBox + -> ( [ ( Box, RootFindingTree Box ) ], ( [ Box ], [ Box ] ) ) +findCuspsFrom + ( CuspFindingOptions + { minWidth + , cuspFindingAlgorithms + } + ) + eqs initBox = runWriter $ map ( initBox , ) <$> go [] initBox where - recur :: ( cand -> Writer ( [ Box ], [ Box ] ) [ IntervalNewtonTree Box ] ) - -> ( [ ( cand, IntervalNewtonTree Box ) ] -> IntervalNewtonTree Box ) - -> [ 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 + go :: BoxHistory + -> Box -- box to work on -> Writer ( [ Box ], [ Box ] ) - [ IntervalNewtonTree Box ] - go cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) ) + [ RootFindingTree Box ] + 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 - -> Writer ( [ Box ], [ Box ] ) - [ IntervalNewtonTree Box ] - continue - cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) ) - , i - , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) ) - ) - | 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 ) + -- Run a round of cusp finding strategies, then recur. + doStrategies + :: BoxHistory + -> NonEmpty RootFindingAlgorithm + -> Box + -> Writer ( [ Box ], [ Box ] ) + [ RootFindingTree Box ] + doStrategies hist0 = do_strats [] 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 ) - i_s_mid = singleton ( ℝ1 s_mid ) - mkGuess ( t0, s0 ) = ( coerce ( (+) @( 𝕀 Double ) ) t0 i_t_mid - , i - , coerce ( (+) @( 𝕀 Double ) ) s0 i_s_mid ) - smaller ( 𝕀 ( ℝ1 t0_lo ) ( ℝ1 t0_hi ), 𝕀 ( ℝ1 s0_lo ) ( ℝ1 s0_hi ) ) - = ( t0_lo + t_mid ) > t_lo + 0.25 * minWidth - || ( t0_hi + t_mid ) < t_hi - 0.25 * minWidth - || ( s0_lo + s_mid ) > s_lo + 0.25 * minWidth - || ( s0_hi + s_mid ) < s_hi - 0.25 * minWidth - 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 ) + 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_skip = not $ canHaveSols l + r_skip = not $ canHaveSols r + d_skip = not $ canHaveSols d + u_skip = not $ canHaveSols u - -- Attempting to implement Algorithm 6 "Heuristic to apply bound-consistency" - -- from the paper - -- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems" - doStrategy box@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ), _, 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) ) - | eps_box1 > eps_box2 - = if | w_t > ebs_box1 || w_s > ebs_box1 - -> makeBox1Consistent box - | w_t > ebs_box2 || w_s > ebs_box2 - -> [ makeBox2Consistent box ] - | otherwise - -> [ box ] - | otherwise - = if | w_t > ebs_box2 || w_s > ebs_box2 - -> [ makeBox2Consistent box ] - | w_t > ebs_box1 || w_s > ebs_box1 - -> makeBox1Consistent box - | otherwise - -> [ box ] - where - eps_box1 = 0.4 - eps_box2 = 0.1 - w_t = t_hi - t_lo - w_s = s_hi - s_lo + ( 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 ) ) + | otherwise + = let why' = case why of + "" -> "" + _ -> " (" ++ why ++ ")" + in case dim of + 0 -> ( [ l, r ], ( "t" ++ why', t_mid ) ) + 1 -> ( [ d, u ], ( "s" ++ why', s_mid ) ) + _ -> error "bisect: fall-back bisection dimension should be either 0 or 1" + in ( bisGuesses, whatBis ) - -- An implementation of "bound-consistency" from the paper - -- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems" - makeBox2Consistent :: Box -> Box - makeBox2Consistent x = ( `State.evalState` False ) $ doLoop 0.25 x - 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' < 0.001 - then return x'' - else do { State.put False ; doLoop lambda' x'' } +-- | Interval Newton method with Gauss–Seidel step. +intervalGaussSeidel + :: ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) + -- ^ equations + -> Preconditioner + -- ^ preconditioner to use for the Gauss–Seidel step + -> Box + -> Writer ( [ Box ], [ Box ] ) [ Box ] +intervalGaussSeidel + eqs + precondMethod + ( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) ) + , i + , 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 ) - -> ( 𝕀 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 ( setter ( 𝕀 x_inf c1 ) box ) of - [] -> c1 - x's -> minimum $ map ( inf . getter ) x's - x'_sup = - case makeBox1Consistent ( 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' >= eps_eq ) $ - State.put True - return $ setter x' box + gsGuesses = gaussSeidel a b + ( coerce ( (-) @( 𝕀 Double ) ) t i_t_mid + , coerce ( (-) @( 𝕀 Double ) ) s i_s_mid ) + in + -- 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) + return todo + where + t_mid = 0.5 * ( t_lo + t_hi ) + s_mid = 0.5 * ( s_lo + s_hi ) + i_t_mid = singleton ( ℝ1 t_mid ) + i_s_mid = singleton ( ℝ1 s_mid ) + 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 - -- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems" - -- - -- See also - -- "Presentation of a highly tuned multithreaded interval solver for underdetermined and well-determined nonlinear systems" - makeBox1Consistent :: Box -> [ Box ] - makeBox1Consistent x = - ( `State.evalState` False ) $ - pipeFunctionsWhileTrue (allNarrowingOperators eps_eq eps_bis eqs) x - eps_eq = 1e-5 -- TODO: this is an additive constant, - -- but I think we also want a multiplicative constant? - eps_bis = minWidth +-- | An implementation of "bc_enforce" from the paper +-- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems" +-- +-- See also +-- "Presentation of a highly tuned multithreaded interval solver for underdetermined and well-determined nonlinear systems" +makeBox1Consistent + :: ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) + -> Double -> Double + -> Box -> [ Box ] +makeBox1Consistent eqs minWidth epsEq x = + ( `State.evalState` False ) $ + 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 ] -bisect x@( 𝕀 x_inf x_sup ) +bisectI :: 𝕀 Double -> [ 𝕀 Double ] +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 ) ) \ No newline at end of file +set_s ( 𝕀 s_inf s_sup ) ( t, i, _ ) = ( t, i, 𝕀 ( ℝ1 s_inf ) ( ℝ1 s_sup ) ) diff --git a/src/app/MetaBrush/Application.hs b/src/app/MetaBrush/Application.hs index 253626f..08d9617 100644 --- a/src/app/MetaBrush/Application.hs +++ b/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 ) diff --git a/src/app/MetaBrush/Render/Document.hs b/src/app/MetaBrush/Render/Document.hs index 8c78b44..7fc2b1f 100644 --- a/src/app/MetaBrush/Render/Document.hs +++ b/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 rootAlgo fitParams +strokeRenderData + :: 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 $