From 61671dc280d8dee8aeabb977ab9b68e8fc4716b2 Mon Sep 17 00:00:00 2001 From: sheaf Date: Wed, 13 Mar 2024 18:00:37 +0100 Subject: [PATCH] Implement box(1) consistency check for cusp finding --- brush-strokes/src/cusps/bench/Main.hs | 27 +- brush-strokes/src/lib/Math/Bezier/Stroke.hs | 269 ++++++++++++++++---- brush-strokes/src/lib/Math/Interval.hs | 4 + 3 files changed, 233 insertions(+), 67 deletions(-) diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs index c58128f..1fb8cc4 100644 --- a/brush-strokes/src/cusps/bench/Main.hs +++ b/brush-strokes/src/cusps/bench/Main.hs @@ -610,32 +610,13 @@ x2 = 𝕀 x2_lo x2_hi -} {- -SMOKING GUN: -> (f, fI) = brushStrokeFunctions $ ellipseBrushStroke (0,1) pi -> (t, i, s) = mkBox (0.54, 0.55) 2 (0.5480, 0.5481) -> _D32_dxdxdx $ dbrush $ eval fI (t,i,s) -> [ℝ2 54187.61174031432 -11626.47655724034, ℝ2 55257.54384135226 -9759.09092706883] +bound consistency / box consistency -In Mathematica: +ghci> (f, fI) = brushStrokeFunctions $ ellipseBrushStroke (0,1) pi +ghci> ee $ eval fI $ mkBox (0, succFP 0) 2 (0.5,0.75) +D22 {_D22_v = [ℝ1 -4379.198802308688, ℝ1 -1791.4548164846506], _D22_dx = T[ℝ1 -3417.0023260831304, ℝ1 6215.362085329192], _D22_dy = T[ℝ1 3524.130928930924, ℝ1 8849.881134814816], _D22_dxdx = T[ℝ1 -13397.019488083239, ℝ1 61439.0409103 - FindMinimum[{D[b[vt, s, myPb0, myPb1][[1]], {vt, 3}] /. (vt -> t), - 0.54 <= t <= 0.55 && 0.548 <= s <= 0.5481}, {{t, 0.545}, {s, - 0.54805}}] -> { -503.445, {t -> 0.54, s -> 0.548} } - FindMaximum[{D[b[vt, s, myPb0, myPb1][[1]], {vt, 3}] /. (vt -> t), - 0.54 <= t <= 0.55 && 0.548 <= s <= 0.5481}, {{t, 0.545}, {s, - 0.54805}}] -> { -480.086, {t -> 0.55, s -> 0.5481} } - -so the actual range is [-503.445,-480.086] but we have computed a range of -[54187.61174031432, 55257.54384135226] which does not contain the actual range! - --} - -{- - -Math.Ring.cos ( linearD @_ @3 (\i -> 𝕀 (getR1 $ inf i) (getR1 $ sup i)) (𝕀 (ℝ1 $ 0.548 * pi) (ℝ1 $ 0.5481 * pi)) :: D3𝔸1 ( 𝕀 Double ) ) -} \ No newline at end of file diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke.hs b/brush-strokes/src/lib/Math/Bezier/Stroke.hs index 6113b74..38d9dcd 100644 --- a/brush-strokes/src/lib/Math/Bezier/Stroke.hs +++ b/brush-strokes/src/lib/Math/Bezier/Stroke.hs @@ -35,7 +35,7 @@ import Control.Arrow import Control.Concurrent.MVar ( MVar, newMVar ) import Control.Monad - ( unless ) + ( unless, when ) import Control.Monad.ST ( RealWorld, ST ) import Data.Bifunctor @@ -58,6 +58,8 @@ import Data.Maybe ( fromMaybe, isJust, listToMaybe, mapMaybe ) import Data.Semigroup ( sconcat ) +import Data.Traversable + ( for ) import GHC.Exts ( newMutVar#, runRW#, inline , Proxy#, proxy# @@ -108,8 +110,8 @@ import Control.Monad.Trans.Class ( lift ) import Control.Monad.Trans.Except ( Except, runExcept, throwE ) -import Control.Monad.Trans.State.Strict - ( StateT, runStateT, evalStateT, get, put ) +import Control.Monad.Trans.State.Strict as State + ( StateT, State, runStateT, evalState, evalStateT, get, put ) import Control.Monad.Trans.Writer.CPS ( Writer, WriterT, execWriterT, runWriter, tell ) @@ -138,6 +140,8 @@ import Math.Differentiable ( Differentiable, DiffInterp, I ) import Math.Epsilon ( epsilon, nearZero ) +import Math.Float.Utils + ( succFP, prevFP ) import Math.Interval import Math.Linear import Math.Module @@ -1493,7 +1497,7 @@ intervalNewtonGSFrom -> ( [ ( Box, IntervalNewtonTree Box ) ], ( [ Box ], [ Box ] ) ) intervalNewtonGSFrom precondMethod minWidth eqs initBox = runWriter $ - map ( initBox , ) <$> evalStrokeDataAndGo initBox + map ( initBox , ) <$> go initBox where @@ -1505,19 +1509,13 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox = rest <- traverse ( \ c -> do { trees <- r c; return [ (c, tree) | tree <- trees ] } ) cands return [ f $ concat rest ] - evalStrokeDataAndGo :: Box -> Writer ( [Box], [Box] ) [ IntervalNewtonTree Box ] - evalStrokeDataAndGo box@( t, i, s ) = go ( box, ( eqs t `Seq.index` i ) s ) - - go :: ( Box, StrokeDatum 3 𝕀 ) -- box to work on + go :: Box -- box to work on -> Writer ( [ Box ], [ Box ] ) [ IntervalNewtonTree Box ] - go ( cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) ) - , i - , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) ) - ) - , sd@( StrokeDatum { ee = D22 _ ( T ee_t ) (T ee_s ) _ _ _ - , 𝛿E𝛿sdcdt = D12 _ ( T ( T f_t ) ) ( T ( T f_s ) ) }) - ) + go cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) ) + , i + , s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) ) + ) -- Box is small: stop processing it. | t_hi - t_lo < minWidth && s_hi - s_lo < minWidth = do let res = TooSmall cand @@ -1527,6 +1525,24 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox = ee_potential_zero sd -- Check the 𝛿E𝛿sdcdt interval contains zero. , 𝛿E𝛿sdcdt_potential_zero sd + = continue cand + -- Box doesn't contain a solution: discard it. + | otherwise + = return [ IntervalNewtonLeaf $ NoSolution ( if ee_potential_zero sd then "dc/dt" else "ee" ) 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 @@ -1564,7 +1580,7 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox = in do tell ([], done) case todo of [] -> return [ IntervalNewtonLeaf $ NoSolution "GaussSeidel" cand ] - _ -> recur evalStrokeDataAndGo ( IntervalNewtonStep ( IntervalNewtonContraction done ) ) + _ -> recur go ( IntervalNewtonStep ( IntervalNewtonContraction done ) ) todo else -- Gauss–Seidel failed to shrink the boxes, so bisect instead. @@ -1605,40 +1621,23 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox = | d_skip && u_skip = ( [], ( "du", s_mid ) ) | l_skip - = ( [ ( r, r_dat ) ], ( "r", t_mid ) ) + = ( [ r ], ( "r", t_mid ) ) | r_skip - = ( [ ( l, l_dat ) ], ( "l", t_mid ) ) + = ( [ l ], ( "l", t_mid ) ) | d_skip - = ( [ ( u, u_dat ) ], ( "u", s_mid ) ) + = ( [ u ], ( "u", s_mid ) ) | u_skip - = ( [ ( d, d_dat ) ], ( "d", s_mid ) ) + = ( [ d ], ( "d", s_mid ) ) | tJWidth >= sJWidth - -- t_hi - t_lo >= s_hi - s_lo - = ( [ ( l, l_dat ), ( r, r_dat ) ], ( "t", t_mid ) ) + -- ... 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, d_dat ), ( u, u_dat ) ], ( "s", s_mid ) ) - in recur go ( IntervalNewtonStep ( IntervalNewtonBisection whatBis ) . map (first fst) ) - bisGuesses - - -- Box doesn't contain a solution: discard it. - | otherwise - = return [ IntervalNewtonLeaf $ NoSolution ( if ee_potential_zero sd then "dc/dt" else "ee" ) cand ] + = ( [ d, u ], ( "s", s_mid ) ) + in recur go ( IntervalNewtonStep ( IntervalNewtonBisection whatBis ) ) + ( makeBox1Consistent =<< bisGuesses ) where - midI :: 𝕀ℝ 2 -> 𝕀ℝ 2 - midI ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = - singleton $ ℝ2 ( 0.5 * ( x_lo + x_hi ) ) ( 0.5 * ( y_lo + y_hi ) ) - --width :: 𝕀ℝ 1 -> Double - --width ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) = hi - lo - --normI :: 𝕀ℝ 1 -> Double - --normI ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) = sqrt $ sup $ ( 𝕀 lo hi ) Ring.^ 2 - --normVI :: 𝕀ℝ 2 -> Double - --normVI ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = - -- sqrt $ sup $ ( 𝕀 x_lo x_hi ) Ring.^ 2 + ( 𝕀 y_lo y_hi ) Ring.^ 2 - normVI3 :: 𝕀ℝ 1 -> 𝕀ℝ 2 -> Double - normVI3 ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) - = sqrt $ max ( abs lo ) ( abs hi ) Ring.^ 2 - + max ( abs x_lo ) ( abs x_hi ) Ring.^ 2 - + max ( abs y_lo ) ( abs y_hi ) Ring.^ 2 t_mid = 0.5 * ( t_lo + t_hi ) s_mid = 0.5 * ( s_lo + s_hi ) i_t_mid = singleton ( ℝ1 t_mid ) @@ -1656,6 +1655,19 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox = !( 𝕀 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 + ee_potential_zero :: StrokeDatum 3 𝕀 -> Bool ee_potential_zero dat = inf ( _D22_v $ ee dat ) <= ℝ1 0 @@ -1716,3 +1728,172 @@ matMulVec ( ℝ2 a11 a21, ℝ2 a12 a22 ) ( 𝕀 ( ℝ2 u_lo v_lo ) ( ℝ2 u_hi v cmpℝ2 :: ( Double -> Double -> Bool ) -> ℝ 2 -> ℝ 2 -> Bool cmpℝ2 cmp ( ℝ2 x1 y1 ) ( ℝ2 x2 y2 ) = cmp x1 x2 && cmp y1 y2 + +midI :: 𝕀ℝ 2 -> 𝕀ℝ 2 +midI ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = + singleton $ ℝ2 ( 0.5 * ( x_lo + x_hi ) ) ( 0.5 * ( y_lo + y_hi ) ) +--width :: 𝕀ℝ 1 -> Double +--width ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) = hi - lo +--normI :: 𝕀ℝ 1 -> Double +--normI ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) = sqrt $ sup $ ( 𝕀 lo hi ) Ring.^ 2 +--normVI :: 𝕀ℝ 2 -> Double +--normVI ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = +-- sqrt $ sup $ ( 𝕀 x_lo x_hi ) Ring.^ 2 + ( 𝕀 y_lo y_hi ) Ring.^ 2 +normVI3 :: 𝕀ℝ 1 -> 𝕀ℝ 2 -> Double +normVI3 ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = + sqrt $ max ( abs lo ) ( abs hi ) Ring.^ 2 + + max ( abs x_lo ) ( abs x_hi ) Ring.^ 2 + + max ( abs y_lo ) ( abs y_hi ) Ring.^ 2 + +-- Use the univariate interval Newton method to narrow from the left +-- a candidate interval. +-- +-- See Algorithm 5 (Procedure left_narrow) in +-- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems" +-- by BartΕ‚omiej Jacek Kubica, 2017 +leftNarrow :: Double + -> Double + -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) ) + -> 𝕀 Double + -> [ 𝕀 Double ] +leftNarrow eps_equal eps_bisection ff' = left_narrow + where + left_narrow ( 𝕀 x_inf x_sup ) = + go x_sup ( 𝕀 x_inf ( if x_inf == x_sup then x_inf else succFP x_inf ) ) + go x_sup x_left = + let ( f_x_left, _f'_x_left ) = ff' x_left + in + if inf f_x_left <= 0 && sup f_x_left >= 0 + then [ 𝕀 ( inf x_left ) x_sup ] + else + let x = 𝕀 ( sup x_left ) x_sup + ( _f_x, f'_x ) = ff' x + x's = do delta <- f_x_left `extendedDivide` f'_x --f'_x_left + let x_new = x_left - delta + map fst $ x_new `intersect` x + in + if | null x's + -> [] + | ( width x - sum ( map width x's ) ) < eps_equal + -- TODO: do a check with the relative reduction in size? + -> x's + | otherwise + -> do + x' <- x's + if sup x' - inf x' < eps_bisection + then return x' + else left_narrow =<< bisect x' + +-- TODO: de-duplicate with 'leftNarrow'? +rightNarrow :: Double + -> Double + -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) ) + -> 𝕀 Double + -> [ 𝕀 Double ] +rightNarrow eps_equal eps_bisection ff' = right_narrow + where + right_narrow ( 𝕀 x_inf x_sup ) = + go x_inf ( 𝕀 ( if x_inf == x_sup then x_sup else prevFP x_sup ) x_sup ) + go x_inf x_right = + let ( f_x_right, _f'_x_right ) = ff' x_right + in + if inf f_x_right <= 0 && sup f_x_right >= 0 + then [ 𝕀 x_inf ( sup x_right ) ] + else + let x = 𝕀 x_inf ( inf x_right ) + ( _f_x, f'_x ) = ff' x + x's = do delta <- f_x_right `extendedDivide` f'_x --f'_x_right + let x_new = x_right - delta + map fst $ x_new `intersect` x + in + if | null x's + -> [] + | ( width x - sum ( map width x's ) ) < eps_equal + -- TODO: do a check with the relative reduction in size? + -> x's + | otherwise + -> do + x' <- x's + if sup x' - inf x' < eps_bisection + then return x' + else right_narrow =<< bisect x' + +bisect :: 𝕀 Double -> [ 𝕀 Double ] +bisect x@( 𝕀 x_inf x_sup ) + | x_inf == x_sup + = [ x ] + | otherwise + = [ 𝕀 x_inf x_mid, 𝕀 x_mid x_sup ] + where x_mid = 0.5 * ( x_inf + x_sup ) + +-- | Apply each function in turn, piping the output of one function into +-- the next. +-- +-- Once all the functions have been applied, check whether the Bool is True. +-- If it is, go around again with all the functions; otherwise, stop. +pipeFunctionsWhileTrue :: [ a -> State Bool [ a ] ] -> a -> State Bool [ a ] +pipeFunctionsWhileTrue fns = go fns + where + go [] x = do + doAnotherRound <- State.get + if doAnotherRound + then do { State.put False ; go fns x } + else return [ x ] + go ( f : fs ) x = do + xs <- f x + concat <$> traverse ( go fs ) xs + +allNarrowingOperators + :: Double + -> Double + -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) + -> [ Box -> State Bool [ Box ] ] +allNarrowingOperators eps_eq eps_bis eqs = + [ \ cand -> + let newCands = map ( `setter` cand ) $ narrowFn eps_eq eps_bis ( \ x0 -> ff' $ setter x0 cand ) ( getter cand ) + in do + when ( ( width ( getter cand ) - sum ( map ( width . getter ) newCands ) ) > eps_eq ) $ + -- TODO: do a check with the relative reduction in size? + State.put True + return newCands + | narrowFn <- [ leftNarrow, rightNarrow ] + , ( ( getter, setter ), ff' ) <- + [ ( ( get_t, set_t ), ff' ) | ff' <- [ ff'_t, g1g1'_t, g2g2'_t ] ] + ++ + [ ( ( get_s, set_s ), ff' ) | ff' <- [ ff'_s, g1g1'_s, g2g2'_s ] ] + ] + where + get_t ( 𝕀 ( ℝ1 t_inf ) ( ℝ1 t_sup ), _, _ ) = 𝕀 t_inf t_sup + 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 ) ) + ff'_t (t, i, s) = + 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 + in ( 𝕀 ee_inf ee_sup, 𝕀 ee'_inf ee'_sup ) + ff'_s (t, i, s) = + 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 + in ( 𝕀 ee_inf ee_sup, 𝕀 ee'_inf ee'_sup ) + g1g1'_t (t, i, s) = + 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 + in ( 𝕀 x_inf x_sup, 𝕀 x'_inf x'_sup ) + g1g1'_s (t, i, s) = + 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 + in ( 𝕀 x_inf x_sup, 𝕀 x'_inf x'_sup ) + g2g2'_t (t, i, s) = + 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 + in ( 𝕀 y_inf y_sup, 𝕀 y'_inf y'_sup ) + g2g2'_s (t, i, s) = + 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 + in ( 𝕀 y_inf y_sup, 𝕀 y'_inf y'_sup ) diff --git a/brush-strokes/src/lib/Math/Interval.hs b/brush-strokes/src/lib/Math/Interval.hs index df141a6..b003d06 100644 --- a/brush-strokes/src/lib/Math/Interval.hs +++ b/brush-strokes/src/lib/Math/Interval.hs @@ -6,6 +6,7 @@ module Math.Interval ( 𝕀(𝕀), inf, sup + , width , scaleInterval , 𝕀ℝ , singleton, nonDecreasing @@ -51,6 +52,9 @@ type instance D k ( 𝕀 v ) = D k v singleton :: a -> 𝕀 a singleton a = 𝕀 a a +width :: AbelianGroup a => 𝕀 a -> a +width ( 𝕀 lo hi ) = hi - lo + -- | Turn a non-decreasing function into a function on intervals. nonDecreasing :: ( a -> b ) -> 𝕀 a -> 𝕀 b nonDecreasing f ( 𝕀 lo hi ) = 𝕀 ( f lo ) ( f hi )