From bd468fcf82e9275680d812e80ef921104eea6afd Mon Sep 17 00:00:00 2001 From: sheaf Date: Fri, 5 Apr 2024 17:48:53 +0200 Subject: [PATCH] start modularising cusp finding code --- brush-strokes/src/cusps/bench/Main.hs | 29 +- brush-strokes/src/lib/Math/Bezier/Stroke.hs | 349 +++++++------------- 2 files changed, 141 insertions(+), 237 deletions(-) diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs index 760fcfb..5459d36 100644 --- a/brush-strokes/src/cusps/bench/Main.hs +++ b/brush-strokes/src/cusps/bench/Main.hs @@ -118,7 +118,9 @@ benchTestCase ( TestCase { testName, testBrushStroke, testCuspOptions, testStart , "Test case " ++ testName ] before <- getMonotonicTime let ( _, testStrokeFnI ) = brushStrokeFunctions testBrushStroke - ( _newtTrees, ( dunno, sols ) ) = computeCusps testCuspOptions testStrokeFnI testStartBoxes + ( dunno, sols ) = + foldMap ( \ ( i, ( _trees, ( mbCusps, defCusps ) ) ) -> ( map ( i , ) mbCusps, map ( i, ) defCusps ) ) $ + computeCusps testCuspOptions testStrokeFnI testStartBoxes rnf dunno `seq` rnf sols `seq` return () after <- getMonotonicTime let dt = after - before @@ -192,7 +194,7 @@ data TestCase = { testName :: String , testBrushStroke :: BrushStroke , testCuspOptions :: CuspFindingOptions - , testStartBoxes :: [ Box ] + , testStartBoxes :: [ ( Int, Box ) ] } brushStrokeFunctions @@ -216,9 +218,9 @@ mkVal t i s = ( ℝ1 t, i, ℝ1 s ) mkI :: ( Double, Double ) -> 𝕀 Double mkI ( lo, hi ) = 𝕀 lo hi -mkBox :: ( Double, Double ) -> Int -> ( Double, Double ) -> Box -mkBox ( t_min, t_max ) i ( s_min, s_max ) = - ( 𝕀 ( ℝ1 t_min ) ( ℝ1 t_max ) , i, 𝕀 ( ℝ1 s_min ) ( ℝ1 s_max ) ) +mkBox :: ( Double, Double ) -> ( Double, Double ) -> Box +mkBox ( t_min, t_max ) ( s_min, s_max ) = + ( 𝕀 ( ℝ2 t_min s_min ) ( ℝ2 t_max s_max ) ) zero, one :: Double zero = 1e-6 @@ -480,7 +482,7 @@ showD float = showFFloat (Just 6) float "" -------------------------------------------------------------------------------- -ellipseTestCase :: CuspFindingOptions -> String -> ( Double, Double ) -> Double -> [ Box ] -> TestCase +ellipseTestCase :: CuspFindingOptions -> String -> ( Double, Double ) -> Double -> [ ( Int, Box ) ] -> TestCase ellipseTestCase opts str k0k1 rot startBoxes = TestCase { testName = "ellipse (" ++ str ++ ")" @@ -602,14 +604,19 @@ getStrokeFunctions brush sp0 crv = computeCusps :: CuspFindingOptions -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) - -> [ Box ] - -> ( [ ( Box, RootFindingTree Box ) ], ( [ Box ], [ Box ] ) ) + -> [ ( Int, Box ) ] + -> [ ( Int, ( [ ( Box, RootFindingTree Box ) ], ( [ Box ], [ Box ] ) ) ) ] computeCusps cuspOpts eqs startBoxes = - foldMap ( findCuspsFrom cuspOpts eqs ) startBoxes + map ( \ ( i, box ) -> ( i , ) $ findCuspsFrom cuspOpts ( mkEqn i ) box ) startBoxes + where + mkEqn i ( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) = + let t = 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) + s = 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) + in ( eqs t `Seq.index` i ) s -defaultStartBoxes :: [ Int ] -> [ Box ] +defaultStartBoxes :: [ Int ] -> [ ( Int, Box ) ] defaultStartBoxes is = - [ mkBox (zero, one) i (zero, one) | i <- is ] + [ ( i, mkBox (zero, one) (zero, one) ) | i <- is ] getR1 (ℝ1 u) = u diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke.hs b/brush-strokes/src/lib/Math/Bezier/Stroke.hs index d6b4373..d483e24 100644 --- a/brush-strokes/src/lib/Math/Bezier/Stroke.hs +++ b/brush-strokes/src/lib/Math/Bezier/Stroke.hs @@ -33,8 +33,6 @@ module Math.Bezier.Stroke -- base import Control.Arrow ( first, (***) ) -import Control.Concurrent.MVar - ( MVar, newMVar ) import Control.Monad ( unless, when ) import Control.Monad.ST @@ -50,7 +48,7 @@ import Data.Foldable import Data.Functor.Identity ( Identity(..) ) import Data.List - ( intercalate, nub, partition, sort ) + ( nub, partition, sort ) import Data.List.NonEmpty ( NonEmpty ) import qualified Data.List.NonEmpty as NE @@ -71,8 +69,6 @@ import GHC.TypeNats ( type (-) ) import Numeric ( showFFloat ) -import System.IO.Unsafe - ( unsafePerformIO ) -- acts import Data.Act @@ -83,12 +79,16 @@ import Data.Act ) -- containers +import Data.IntMap.Strict + ( IntMap ) +import qualified Data.IntMap.Strict as IntMap + ( fromList, toList ) import Data.Sequence ( Seq(..) ) import qualified Data.Sequence as Seq ( empty, index, length, reverse, singleton, zipWith ) import Data.Tree - ( Tree(..), foldTree ) + ( Tree(..) ) -- deepseq import Control.DeepSeq @@ -114,10 +114,6 @@ import Control.Monad.Trans.State.Strict as State import Control.Monad.Trans.Writer.CPS ( Writer, WriterT, execWriterT, runWriter, tell ) --- tree-view -import Data.Tree.View - ( showTree ) - -- MetaBrush import Math.Algebra.Dual import qualified Math.Bezier.Cubic as Cubic @@ -154,7 +150,7 @@ import Math.Orientation import qualified Math.Ring as Ring import Math.Roots -import Debug.Utils +--import Debug.Utils -------------------------------------------------------------------------------- @@ -604,56 +600,24 @@ outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams brushFromParams = $ runD ( brushFromParams @2 @() proxy# id ) $ toBrushParams params_t - ( newtTrees, ( newtDunno, newtSols ) ) = + ( potentialCusps, definiteCusps ) = case mbCuspOptions of Just cuspOptions -> - findCusps cuspOptions curvesI + foldMap + ( \ ( i, ( _trees, ( potCusps, defCusps ) ) ) -> + ( map ( i , ) potCusps, map ( i , ) defCusps ) + ) + ( IntMap.toList $ findCusps cuspOptions curvesI ) Nothing -> - ( [], ( [], [] ) ) + ( [], [] ) - showD :: Double -> String - showD float = showFFloat (Just 6) float "" - functionDataLines = - [ "E, dE/ds, dE/ds * dc/dt:" - , "{" ++ - (intercalate "," - [ "{" ++ showD t ++ "," ++ showD (s + fromIntegral i) ++ ",{" ++ intercalate "," vals ++ "}}" - | t <- [0.001, 0.01.. 0.999] - , i <- [0 .. (fromIntegral $ length $ curves (ℝ1 0.5)) - 1] - , s <- [0.001, 0.01.. 0.999] - , let StrokeDatum - { ee = D12 (ℝ1 e) _dEdt (T (ℝ1 dEds)) - , 𝛿E𝛿sdcdt = D0 (T (ℝ2 vx vy)) - } = (curves (ℝ1 t) `Seq.index` i) (ℝ1 s) - vals = [showD e, showD dEds, "{" ++ showD vx ++ "," ++ showD vy ++ "}"] - ] - ) ++ "}" - ] - showedTrees = map ( uncurry showRootFindingTree ) newtTrees - solLines = - [ " #sols: " ++ show (length newtSols) - , "#dunno: " ++ show (length newtDunno) - , "Tree size: " ++ show (sum @_ @Int $ map (foldTree ( \ _ bs -> 1 + sum bs )) showedTrees) - ] - treeLines = - solLines ++ map showTree showedTrees - - logContents = unlines $ functionDataLines ++ treeLines - - in --trace (unlines solLines) $ - --logToFile cuspFindingMVar logContents `seq` - OutlineInfo + in OutlineInfo { outlineFn = fwdBwd - , outlineDefiniteCusps = map ( cuspCoords curves ) newtSols - , outlinePotentialCusps = map ( cuspCoords curves ) newtDunno + , outlineDefiniteCusps = map ( cuspCoords curves ) definiteCusps + , outlinePotentialCusps = map ( cuspCoords curves ) potentialCusps } {-# INLINEABLE outlineFunction #-} -cuspFindingMVar :: MVar FilePath -cuspFindingMVar = unsafePerformIO - $ newMVar "logs/cusps.txt" -{-# NOINLINE cuspFindingMVar #-} - pathAndUsedParams :: forall k i arr crvData ptData usedParams . ( HasType ( ℝ 2 ) ptData , HasBΓ©zier k i @@ -1015,11 +979,6 @@ data RootSolvingAlgorithm -- | Use the modified Halley M2 method. | HalleyM2 { maxIters :: Word, precision :: Int } -rootSolvingMVar :: MVar FilePath -rootSolvingMVar = unsafePerformIO - $ newMVar "logs/envelope.txt" -{-# NOINLINE rootSolvingMVar #-} - -- | Solve the envelope equations at a given point \( t = t_0 \), to find -- \( s_0 \) such that \( c(t_0, s_0) \) is on the envelope of the brush stroke. solveEnvelopeEquations :: RootSolvingAlgorithm @@ -1055,68 +1014,21 @@ solveEnvelopeEquations rootAlgo _t path_t path'_t ( fwdOffset, bwdOffset ) strok plausibleTangent tgt = path'_t ^.^ tgt > 0 sol :: String -> Seq ( ℝ 1 -> StrokeDatum 2 () ) -> Double -> ( Bool, ℝ 2, T ( ℝ 2 ) ) - sol desc f is0 = - let (solRes, solSteps) = runSolveMethod ( eqn f ) is0 + sol _desc f is0 = + let (solRes, _solSteps) = runSolveMethod ( eqn f ) is0 ( good, is ) = case solRes of Nothing -> ( False, is0 ) Just is1 -> ( True , is1 ) ( ds, dcdt ) = finish f is - showD :: Double -> String - showD float = showFFloat (Just 6) float "" - (x_min, x_max) = domain - logContents = unlines - [ "method = " <> methodName - , desc - , "t = " ++ show _t - , "start (n,s) = " <> show (fromDomain is0) - , "steps = {" - , unlines $ map (\ (x, y, y') -> show (x, (y, y'))) solSteps - , "}," - , "E:" - , "{" ++ - (intercalate "," - [ "{" ++ showD x ++ "," ++ showD y ++ "}" - | x <- [x_min, x_min + 0.01 .. x_max] - , let (# y, _y' #) = eqn f x - ] - ) ++ "}" - , "𝛿E𝛿s:" - , "{" ++ - (intercalate "," - [ "{" ++ showD x ++ "," ++ showD y' ++ "}" - | x <- [x_min, x_min + 0.01 .. x_max] - , let (# _y, y' #) = eqn f x - ] - ) ++ "}" - , "𝛿E𝛿sdcdt (x):" - , "{" ++ - (intercalate "," - [ "{" ++ showD x ++ "," ++ showD 𝛿E𝛿sdcdt_x ++ "}" - | x <- [x_min, x_min + 0.01 .. x_max] - , let StrokeDatum { 𝛿E𝛿sdcdt = D0 (T (ℝ2 𝛿E𝛿sdcdt_x _)) } = evalStrokeDatum f x - ] - ) - ++ "}" - , "𝛿E𝛿sdcdt (y):" - , "{" ++ - (intercalate "," - [ "{" ++ showD x ++ "," ++ showD 𝛿E𝛿sdcdt_y ++ "}" - | x <- [x_min, x_min + 0.01 .. x_max] - , let StrokeDatum { 𝛿E𝛿sdcdt = D0 (T (ℝ2 _ 𝛿E𝛿sdcdt_y)) } = evalStrokeDatum f x - ] - ) - ++ "}" - ] - in --logToFile rootSolvingMVar logContents `seq` - ( good, ds, dcdt ) + in ( good, ds, dcdt ) - (runSolveMethod, methodName) = case rootAlgo of + runSolveMethod = case rootAlgo of HalleyM2 { maxIters, precision } -> - (halleyM2 maxIters precision, "HalleyM2") + halleyM2 maxIters precision NewtonRaphson { maxIters, precision } -> - (newtonRaphson maxIters precision domain, "NewtonRaphson") + newtonRaphson maxIters precision domain finish :: Seq ( ℝ 1 -> StrokeDatum 2 () ) -> Double -> ( ℝ 2, T ( ℝ 2 ) ) finish fs is = @@ -1244,13 +1156,13 @@ brushStrokeData co1 co2 path params brush = -- The boolean indicates whether the Gauss–Seidel step was a contraction. gaussSeidel :: ( 𝕀ℝ 2, 𝕀ℝ 2 ) -- ^ columns of \( A \) -> 𝕀ℝ 2 -- ^ \( B \) - -> ( 𝕀ℝ 1, 𝕀ℝ 1 ) -- ^ initial box \( X \) - -> [ ( ( 𝕀ℝ 1, 𝕀ℝ 1 ), Bool ) ] + -> 𝕀ℝ 2 -- ^ initial box \( X \) + -> [ ( 𝕀ℝ 2, Bool ) ] gaussSeidel ( 𝕀 ( ℝ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 ) ) - ( 𝕀 ( ℝ1 x1_lo ) ( ℝ1 x1_hi ), 𝕀 ( ℝ1 x2_lo ) ( ℝ1 x2_hi ) ) + ( 𝕀 ( ℝ2 x1_lo x2_lo ) ( ℝ2 x1_hi x2_hi ) ) = let !a11 = 𝕀 a11_lo a11_hi !a12 = 𝕀 a12_lo a12_hi !a21 = 𝕀 a21_lo a21_hi @@ -1268,7 +1180,7 @@ gaussSeidel x2'_0 <- ( b2 - a21 * x1' ) `extendedDivide` a22 ( 𝕀 x2'_lo x2'_hi, sub_x2 ) <- x2'_0 `intersect` x2 - return ( ( 𝕀 ( ℝ1 x1'_lo ) ( ℝ1 x1'_hi ), 𝕀 ( ℝ1 x2'_lo ) ( ℝ1 x2'_hi ) ) + return ( ( 𝕀 ( ℝ2 x1'_lo x2'_lo) ( ℝ2 x1'_hi x2'_hi ) ) , sub_x1 && sub_x2 ) -- TODO: try implementing the complete interval union Gauss–Seidel algorithm. -- See "Algorithm 2" in @@ -1294,9 +1206,9 @@ intersect ( 𝕀 lo1 hi1 ) ( 𝕀 lo2 hi2 ) -- 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 + -> ( Int, Box ) -> Cusp -cuspCoords eqs ( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ), i, 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) ) +cuspCoords eqs ( i, 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) | StrokeDatum { dpath , dstroke = D22 { _D22_v = stroke } } @@ -1342,7 +1254,6 @@ instance Show RootFindingStep where data RootFindingLeaf d = TooSmall d - | NoSolution String d deriving stock Show showRootFindingTree :: Box -> RootFindingTree Box -> Tree String @@ -1351,13 +1262,13 @@ 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 ) ) +boxArea ( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) = abs ( t_hi - t_lo ) * abs ( s_hi - s_lo ) showArea :: Double -> [Char] showArea area = "(area " ++ showFFloat (Just 6) area "" ++ ")" -type Box = ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) +type Box = 𝕀ℝ 2 type BoxHistory = [ NE.NonEmpty ( RootFindingStep, Box ) ] -- | Options for the cusp-finding methods in 'findCusps' and 'findCuspsFrom'. @@ -1371,17 +1282,24 @@ data CuspFindingOptions } data RootFindingAlgorithm + -- | Bisection step. + = Bisection !BisectionOptions -- | Gauss–Seidel step with the given preconditioning method. - = GaussSeidel GaussSeidelOptions + | GaussSeidel !GaussSeidelOptions -- | @box(1)@-consistency. - | Box1 Box1Options + | Box1 !Box1Options -- | @box(2)@-consistency. - | Box2 Box2Options - -- | Bisection step. - | Bisection BisectionOptions + | Box2 !Box2Options + +-- | Options for the bisection method. +data BisectionOptions = + BisectionOptions + { canHaveSols :: !( ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 ) -> Box -> Bool ) + , fallbackBisectionDim :: !( [ ( RootFindingStep, Box ) ] -> BoxHistory -> ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 ) -> ( Int, String ) ) + } -- | Options for the interval Gauss–Seidel method. -data GaussSeidelOptions = GaussSeidelOptions { gsPreconditioner :: Preconditioner } +data GaussSeidelOptions = GaussSeidelOptions { gsPreconditioner :: !Preconditioner } -- | Options for the @box(1)@-consistency method. data Box1Options = Box1Options { box1EpsEq :: !Double } @@ -1389,21 +1307,15 @@ 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 + , cuspFindingAlgorithms = defaultCuspFindingAlgorithms minWidth narrowAbs } where minWidth = 1e-5 + narrowAbs = 5e-3 defaultCuspFindingAlgorithms :: Double -> Double -> Box -> BoxHistory -> NonEmpty RootFindingAlgorithm defaultCuspFindingAlgorithms minWidth narrowAbs box history = @@ -1417,17 +1329,18 @@ defaultCuspFindingAlgorithms minWidth narrowAbs box history = -- ...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. + -- Otherwise, do a normal round. + -- Currently: we 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 + sufficientlySmallerThan + ( 𝕀 ( ℝ2 t1_lo s1_lo ) ( ℝ2 t1_hi s1_hi ) ) + ( 𝕀 ( ℝ2 t2_lo s2_lo ) ( ℝ2 t2_hi s2_hi ) ) = + ( ( t1_hi - t1_lo ) - ( t2_hi - t2_lo ) > narrowAbs ) + || + ( ( s1_hi - s1_lo ) - ( s2_hi - s2_lo ) > narrowAbs ) defaultGaussSeidelOptions :: GaussSeidelOptions defaultGaussSeidelOptions = GaussSeidelOptions { gsPreconditioner = InverseMidJacobian } @@ -1436,15 +1349,12 @@ defaultBisectionOptions :: Double -> Double -> Box -> BisectionOptions defaultBisectionOptions _minWidth _narrowAbs - ( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) ) - , i - , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) ) - ) + box@( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) = BisectionOptions { canHaveSols = - \ eqs _box'@( _t', _i', _s' ) -> + \ eqs box' -> -- box(0)-consistency - let dat' = ( eqs _t' `Seq.index` _i' ) _s' + let dat' = eqs box' in ee_potential_zero dat' && 𝛿E𝛿sdcdt_potential_zero dat' -- box(1)-consistency @@ -1458,7 +1368,7 @@ defaultBisectionOptions \ _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 + = eqs box wd_t = t_hi - t_lo wd_s = s_hi - s_lo tJWidth = wd_t * normVI3 ee_t f_t @@ -1479,15 +1389,19 @@ defaultBisectionOptions findCusps :: CuspFindingOptions -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) - -> ( [ ( 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 ] - ] + -> IntMap ( [ ( Box, RootFindingTree Box ) ], ( [ Box ], [ Box ] ) ) +findCusps opts eqs = + IntMap.fromList + [ ( i, findCuspsFrom opts eqnPiece initBox ) + | i <- [ 0 .. length ( eqs ( 𝕀 ( ℝ1 zero ) ( ℝ1 one ) ) ) - 1 ] + , let initBox = 𝕀 ( ℝ2 zero zero ) ( ℝ2 one one ) + eqnPiece ( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) = + let t = 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) + s = 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) + in ( eqs t `Seq.index` i ) s + ] --- | Use the following algorithms using interal arithmetic in order +-- | Use the following algorithms using interval arithmetic in order -- to find cusps in the envelope: -- -- - interval Newton method with Gauss–Seidel step for inversion @@ -1502,7 +1416,7 @@ findCusps opts eqs = foldMap ( findCuspsFrom opts eqs ) initBoxes -- boxes which might or might not contain solutions. findCuspsFrom :: CuspFindingOptions - -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) + -> ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 ) -> Box -> ( [ ( Box, RootFindingTree Box ) ], ( [ Box ], [ Box ] ) ) findCuspsFrom @@ -1520,26 +1434,23 @@ findCuspsFrom -> Writer ( [ Box ], [ Box ] ) [ RootFindingTree Box ] go history - cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) ) - , i - , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) ) - ) + cand@( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) + | -- Check the envelope equation interval contains zero. + not ( ee_potential_zero sd) + || + -- Check the 𝛿E𝛿sdcdt interval contains zero. + not ( 𝛿E𝛿sdcdt_potential_zero sd ) + -- Box doesn't contain a solution: discard it. + = return [] -- Box is small: stop processing it. | t_hi - t_lo < minWidth && s_hi - s_lo < minWidth = do let res = TooSmall cand tell ( [ cand ], [] ) 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 - = doStrategies history ( cuspFindingAlgorithms cand history ) cand - -- Box doesn't contain a solution: discard it. | otherwise - , let whyNoSol = if ee_potential_zero sd then "dc/dt" else "ee" - = return [ RootFindingLeaf $ NoSolution whyNoSol cand ] + = doStrategies history ( cuspFindingAlgorithms cand history ) cand where - sd = ( eqs t `Seq.index` i ) s + sd = eqs $ 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) -- Run a round of cusp finding strategies, then recur. doStrategies @@ -1548,26 +1459,26 @@ findCuspsFrom -> Box -> Writer ( [ Box ], [ Box ] ) [ RootFindingTree Box ] - doStrategies hist0 = do_strats [] + doStrategies prevRoundsHist = do_strats [] where do_strats :: [ ( RootFindingStep, Box ) ] -> NE.NonEmpty RootFindingAlgorithm -> Box -> Writer ( [ Box ], [ Box ] ) [ RootFindingTree Box ] - do_strats stratsHist ( algo NE.:| algos ) box = do + do_strats thisRoundHist ( algo NE.:| algos ) box = do -- Run one strategy in the round. - ( step, boxes ) <- doStrategy stratsHist hist0 eqs minWidth algo box + ( step, boxes ) <- doStrategy thisRoundHist prevRoundsHist 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 + let thisRoundHist' = ( step, box ) : thisRoundHist + in recur step ( do_strats thisRoundHist' ( 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 + let thisRoundHist' = ( step, box ) NE.:| thisRoundHist + history = thisRoundHist' : prevRoundsHist in recur step ( go history ) boxes recur :: RootFindingStep @@ -1583,7 +1494,7 @@ findCuspsFrom -- (hopefully smaller) output boxes. doStrategy :: [ ( RootFindingStep, Box ) ] -> BoxHistory - -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) + -> ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 ) -> Double -> RootFindingAlgorithm -> Box @@ -1612,24 +1523,17 @@ bisect :: ( Box -> Bool ) -> Box -> ( [ Box ], ( String, Double ) ) bisect canHaveSols ( dim, why ) - ( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) ) - , i - , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) ) - ) + ( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi 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 ) - 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 = 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_mid s_hi ) + r = 𝕀 ( ℝ2 t_mid s_lo ) ( ℝ2 t_hi s_hi ) + d = 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_mid ) + u = 𝕀 ( ℝ2 t_lo s_mid ) ( ℝ2 t_hi s_hi ) l_skip = not $ canHaveSols l r_skip = not $ canHaveSols r d_skip = not $ canHaveSols d @@ -1660,24 +1564,21 @@ bisect canHaveSols ( dim, why ) -- | Interval Newton method with Gauss–Seidel step. intervalGaussSeidel - :: ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) + :: ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 ) -- ^ equations -> Preconditioner -- ^ preconditioner to use for the Gauss–Seidel step - -> Box - -> Writer ( [ Box ], [ Box ] ) [ Box ] + -> 𝕀ℝ 2 + -> Writer ( [ 𝕀ℝ 2 ], [ 𝕀ℝ 2 ] ) [ 𝕀ℝ 2 ] intervalGaussSeidel eqs precondMethod - ( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) ) - , i - , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) ) - ) + ts@( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) | StrokeDatum { 𝛿E𝛿sdcdt = D12 _ ( T ( T f_t ) ) ( T ( T f_s ) ) } - <- ( eqs t `Seq.index` i ) s + <- eqs ts , 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 + <- eqs ts_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 @@ -1697,9 +1598,7 @@ intervalGaussSeidel -- , 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 ) + gsGuesses = gaussSeidel a b ( unT $ T ts ^-^ T ts_mid ) in -- If the Gauss–Seidel step was a contraction, then the box -- contains a unique solution (by the Banach fixed point theorem). @@ -1713,11 +1612,8 @@ intervalGaussSeidel 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 ) + ts_mid = singleton ( ℝ2 t_mid s_mid ) + mkGuess ts_g = unT $ T ts_g ^+^ T ts_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 @@ -1729,7 +1625,7 @@ intervalGaussSeidel -- See also -- "Presentation of a highly tuned multithreaded interval solver for underdetermined and well-determined nonlinear systems" makeBox1Consistent - :: ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) + :: ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 ) -> Double -> Double -> Box -> [ Box ] makeBox1Consistent eqs minWidth epsEq x = @@ -1739,7 +1635,7 @@ makeBox1Consistent eqs minWidth epsEq 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 𝕀 ) ) + :: ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 ) -> Double -> Double -> Double -> Box -> Box makeBox2Consistent eqs minWidth epsEq lambdaMin x0 = ( `State.evalState` False ) $ doLoop 0.25 x0 @@ -1951,7 +1847,7 @@ pipeFunctionsWhileTrue fns = go fns allNarrowingOperators :: Double -> Double - -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) + -> ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 ) -> [ Box -> State Bool [ Box ] ] allNarrowingOperators eps_eq eps_bis eqs = [ \ cand -> @@ -1968,41 +1864,42 @@ allNarrowingOperators eps_eq eps_bis eqs = [ ( ( get_s, set_s ), ff' ) | ff' <- [ ff'_s, g1g1'_s, g2g2'_s ] ] ] where - ff'_t (t, i, s) = + ff'_t ts = let D22 { _D22_v = 𝕀 ( ℝ1 ee_inf ) ( ℝ1 ee_sup ) , _D22_dx = T ( 𝕀 ( ℝ1 ee'_inf ) ( ℝ1 ee'_sup ) ) - } = ee $ ( eqs t `Seq.index` i ) s + } = ee $ eqs ts in ( 𝕀 ee_inf ee_sup, 𝕀 ee'_inf ee'_sup ) - ff'_s (t, i, s) = + ff'_s ts = let D22 { _D22_v = 𝕀 ( ℝ1 ee_inf ) ( ℝ1 ee_sup ) , _D22_dy = T ( 𝕀 ( ℝ1 ee'_inf ) ( ℝ1 ee'_sup ) ) - } = ee $ ( eqs t `Seq.index` i ) s + } = ee $ eqs ts in ( 𝕀 ee_inf ee_sup, 𝕀 ee'_inf ee'_sup ) - g1g1'_t (t, i, s) = + g1g1'_t ts = let D12 { _D12_v = T ( 𝕀 ( ℝ2 x_inf _ ) ( ℝ2 x_sup _ ) ) , _D12_dx = T ( T ( 𝕀 ( ℝ2 x'_inf _ ) ( ℝ2 x'_sup _ ) ) ) - } = 𝛿E𝛿sdcdt $ ( eqs t `Seq.index` i ) s + } = 𝛿E𝛿sdcdt $ eqs ts in ( 𝕀 x_inf x_sup, 𝕀 x'_inf x'_sup ) - g1g1'_s (t, i, s) = + g1g1'_s ts = let D12 { _D12_v = T ( 𝕀 ( ℝ2 x_inf _ ) ( ℝ2 x_sup _ ) ) , _D12_dy = T ( T ( 𝕀 ( ℝ2 x'_inf _ ) ( ℝ2 x'_sup _ ) ) ) - } = 𝛿E𝛿sdcdt $ ( eqs t `Seq.index` i ) s + } = 𝛿E𝛿sdcdt $ eqs ts in ( 𝕀 x_inf x_sup, 𝕀 x'_inf x'_sup ) - g2g2'_t (t, i, s) = + g2g2'_t ts = let D12 { _D12_v = T ( 𝕀 ( ℝ2 _ y_inf ) ( ℝ2 _ y_sup ) ) , _D12_dx = T ( T ( 𝕀 ( ℝ2 _ y'_inf ) ( ℝ2 _ y'_sup ) ) ) - } = 𝛿E𝛿sdcdt $ ( eqs t `Seq.index` i ) s + } = 𝛿E𝛿sdcdt $ eqs ts in ( 𝕀 y_inf y_sup, 𝕀 y'_inf y'_sup ) - g2g2'_s (t, i, s) = + g2g2'_s ts = let D12 { _D12_v = T ( 𝕀 ( ℝ2 _ y_inf ) ( ℝ2 _ y_sup ) ) , _D12_dy = T ( T ( 𝕀 ( ℝ2 _ y'_inf ) ( ℝ2 _ y'_sup ) ) ) - } = 𝛿E𝛿sdcdt $ ( eqs t `Seq.index` i ) s + } = 𝛿E𝛿sdcdt $ eqs ts in ( 𝕀 y_inf y_sup, 𝕀 y'_inf y'_sup ) +-- TODO: replace this with Fin n and representable please get_t, get_s :: Box -> 𝕀 Double -get_t ( 𝕀 ( ℝ1 t_inf ) ( ℝ1 t_sup ), _, _ ) = 𝕀 t_inf t_sup +get_t ( 𝕀 ( ℝ2 t_inf _ ) ( ℝ2 t_sup _ ) ) = 𝕀 t_inf t_sup +get_s ( 𝕀 ( ℝ2 _ s_inf ) ( ℝ2 _ s_sup ) ) = 𝕀 s_inf s_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_t ( 𝕀 t_inf t_sup ) ( 𝕀 ( ℝ2 _ s_inf ) ( ℝ2 _ s_sup ) ) = 𝕀 ( ℝ2 t_inf s_inf ) ( ℝ2 t_sup s_sup ) +set_s ( 𝕀 s_inf s_sup ) ( 𝕀 ( ℝ2 t_inf _ ) ( ℝ2 t_sup _ ) ) = 𝕀 ( ℝ2 t_inf s_inf ) ( ℝ2 t_sup s_sup )