diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs index 6534688..67b8853 100644 --- a/brush-strokes/src/cusps/bench/Main.hs +++ b/brush-strokes/src/cusps/bench/Main.hs @@ -110,8 +110,7 @@ main = do putStrLn $ unlines [ replicate 40 '=' , "Best time: " ++ show bestTime ++ "s" - , "Best parameters:" - --, show (testCuspOptions bestTest) + , "Best parameters: " ++ show (testName bestTest) ] benchTestCase :: TestCase -> IO Double @@ -181,9 +180,15 @@ testCases = benchCases benchCases :: [ TestCase ] benchCases = - [ ellipseTestCase opts "full" ( 0, 1 ) pi $ defaultStartBoxes [ 2 ] ] -- [ 0 .. 3 ] ] - where - opts = defaultRootIsolationOptions + [ ellipseTestCase opts ("minWidth=" ++ show minWidth ++ ",Ξ΅=" ++ show narrowAbs) ( 0, 1 ) pi $ defaultStartBoxes [ 2 ] + | minWidth <- [ 1e-5 ] + , narrowAbs <- [ 1e-5, 1e-4, 1e-3, 5e-3, 8e-3, 1e-2, 2e-2, 3e-2, 4e-2, 5e-2, 1e-1 ] + , let opts = + RootIsolationOptions + { minWidth + , cuspFindingAlgorithms = defaultRootIsolationAlgorithms minWidth narrowAbs + } + ] -------------------------------------------------------------------------------- diff --git a/brush-strokes/src/lib/Math/Root/Isolation.hs b/brush-strokes/src/lib/Math/Root/Isolation.hs index 45dca47..5bd3d3a 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation.hs @@ -55,7 +55,8 @@ import Data.Semigroup import Numeric ( showFFloat ) import GHC.TypeNats - ( Nat, KnownNat ) + ( Nat, KnownNat + , type (<=) ) -- containers import Data.Tree @@ -285,10 +286,10 @@ defaultGaussSeidelOptions = defaultBisectionOptions :: forall n d - . BoxCt n d + . ( 1 <= n, BoxCt n d ) => Double -> Double -> Box n -> BisectionOptions n d -defaultBisectionOptions _minWidth _narrowAbs box = +defaultBisectionOptions minWidth _narrowAbs box = BisectionOptions { canHaveSols = \ eqs box' -> @@ -320,21 +321,20 @@ defaultBisectionOptions _minWidth _narrowAbs box = -- First, check if the largest dimension is over 10 times larger -- than the smallest dimension; if so bisect along that coordinate. - -- - -- TODO: filter out dimensions smaller than minimum width. in case sortOnArg ( width . coordInterval ) datPerCoord of - [] -> error "dimension 0" + [] -> error "impossible: dimension 0" [Arg _ d] -> (coordIndex d, "") Arg w0 _ : ds -> let Arg w1 d1 = last ds in if w1 >= 10 * w0 - then ( coordIndex d1, "tooWide") + then ( coordIndex d1, "tooWide" ) -- Otherwise, pick the dimension with the largest spread = width * Jacobian column norm - else case sortOnArg ( Dual . spread ) datPerCoord of - [] -> error "dimension 0" - Arg _ d : _ -> - (coordIndex d, "spread") + else case filter ( not . isTooSmall ) $ sortOnArg ( Dual . spread ) datPerCoord of + [] -> ( coordIndex d1, "tooWide'" ) + Arg _ d : _ -> ( coordIndex d, "spread" ) } + where + isTooSmall ( Arg ( Dual w ) _ ) = w < minWidth {-# INLINEABLE defaultBisectionOptions #-} sortOnArg :: Ord b => (a -> b) -> [a] -> [Arg b a] @@ -494,7 +494,7 @@ bisect canHaveSols ( fallbackDim, why ) box = ( oks, ("cd = " ++ show i, mid ) ) Nothing -> case bisectInCoord box fallbackDim of - ( mid, (lo, hi) ) -> + ( mid, ( lo, hi ) ) -> ( [ lo, hi ], ( why, mid ) ) where solsList = @@ -673,30 +673,27 @@ makeBox2Consistent -> Box2Options n d -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) -> Box n -> Box n -makeBox2Consistent minWidth (Box2Options epsEq lambdaMin coordsToNarrow eqsToUse) eqs x0 +makeBox2Consistent minWidth (Box2Options Ξ΅_eq Ξ»Min coordsToNarrow eqsToUse) eqs x0 = ( `State.evalState` False ) $ doLoop 0.25 x0 where box1Options :: Box1Options n d - box1Options = Box1Options epsEq coordsToNarrow eqsToUse + box1Options = Box1Options Ξ΅_eq coordsToNarrow eqsToUse doBox1 :: Box n -> [ Box n ] doBox1 = makeBox1Consistent minWidth box1Options eqs doLoop :: Double -> Box n -> State Bool ( Box n ) - doLoop lambda x = do - x'' <- forEachDim @n x $ \ i -> - boundConsistency ( `index` i ) ( set i ) lambda + doLoop Ξ» x = do + x'' <- forEachDim @n x $ boundConsistency Ξ» modified <- State.get - let lambda' = if modified then lambda else 0.5 * lambda - if lambda' < lambdaMin + let Ξ»' = if modified then Ξ» else 0.5 * Ξ» + if Ξ»' < Ξ»Min then return x'' - else do { State.put False ; doLoop lambda' x'' } + else do { State.put False ; doLoop Ξ»' x'' } - boundConsistency :: ( Box n -> 𝕀 Double ) - -> ( 𝕀 Double -> Box n -> Box n ) - -> Double -> Box n -> State Bool ( Box n ) - boundConsistency getter setter lambda box = do + boundConsistency :: Double -> Fin n -> Box n -> State Bool ( Box n ) + boundConsistency Ξ» i 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 + c1 = ( 1 - Ξ» ) * x_inf + Ξ» * x_sup + c2 = Ξ» * x_inf + ( 1 - Ξ» ) * x_sup x'_inf = case doBox1 ( setter ( 𝕀 x_inf c1 ) box ) of [] -> c1 @@ -706,9 +703,12 @@ makeBox2Consistent minWidth (Box2Options epsEq lambdaMin coordsToNarrow eqsToUse [] -> c2 x's -> maximum $ map ( sup . getter ) x's x' = 𝕀 x'_inf x'_sup - when ( width x - width x' >= epsEq ) $ + when ( width x - width x' >= Ξ΅_eq ) $ State.put True return $ setter x' box + where + getter = ( `index` i ) + setter = set i -- | Pre-condition the system \( AX = B \). precondition @@ -779,7 +779,7 @@ leftNarrow :: Double -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) ) -> 𝕀 Double -> [ 𝕀 Double ] -leftNarrow eps_equal eps_bisection ff' = left_narrow +leftNarrow Ξ΅_eq Ξ΅_bis 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 ) ) @@ -791,19 +791,18 @@ leftNarrow eps_equal eps_bisection ff' = left_narrow 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 + x's = do Ξ΄ <- f_x_left `extendedDivide` f'_x + let x_new = x_left - Ξ΄ 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? + | ( width x - sum ( map width x's ) ) < Ξ΅_eq -> x's | otherwise -> do x' <- x's - if sup x' - inf x' < eps_bisection + if sup x' - inf x' < Ξ΅_bis then return x' else left_narrow =<< bisectI x' @@ -813,7 +812,7 @@ rightNarrow :: Double -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) ) -> 𝕀 Double -> [ 𝕀 Double ] -rightNarrow eps_equal eps_bisection ff' = right_narrow +rightNarrow Ξ΅_eq Ξ΅_bis 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 ) @@ -825,19 +824,18 @@ rightNarrow eps_equal eps_bisection ff' = right_narrow 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 + x's = do Ξ΄ <- f_x_right `extendedDivide` f'_x + let x_new = x_right - Ξ΄ 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? + | ( width x - sum ( map width x's ) ) < Ξ΅_eq -> x's | otherwise -> do x' <- x's - if sup x' - inf x' < eps_bisection + if sup x' - inf x' < Ξ΅_bis then return x' else right_narrow =<< bisectI x' @@ -873,25 +871,21 @@ allNarrowingOperators -> Box1Options n d -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) -> [ Box n -> State Bool [ Box n ] ] -allNarrowingOperators eps_bis ( Box1Options eps_eq coordsToNarrow eqsToUse ) eqs = +allNarrowingOperators Ξ΅_bis ( Box1Options Ξ΅_eq coordsToNarrow eqsToUse ) eqs = [ \ cand -> let getter = ( `index` coordIndex ) setter = set coordIndex - newCands = map ( `setter` cand ) $ narrowFn eps_eq eps_bis ( \ x0 -> fn $ setter x0 cand ) ( getter cand ) + newCands = map ( `setter` cand ) $ narrowFn ( \ box -> ff' coordIndex eqnIndex $ setter box cand ) ( getter cand ) in do - when ( ( width ( getter cand ) - sum ( map ( width . getter ) newCands ) ) > eps_eq ) $ + when ( ( width ( getter cand ) - sum ( map ( width . getter ) newCands ) ) > Ξ΅_eq ) $ -- TODO: do a check with the relative reduction in size? State.put True return newCands - | narrowFn <- [ leftNarrow, rightNarrow ] - , ( coordIndex, fn ) <- - [ ( i, ff' i d ) - | i <- coordsToNarrow - , d <- eqsToUse - ] + | narrowFn <- [ leftNarrow Ξ΅_eq Ξ΅_bis, rightNarrow Ξ΅_eq Ξ΅_bis ] + , coordIndex <- coordsToNarrow + , eqnIndex <- eqsToUse ] where - ff' :: Fin n -> Fin d -> 𝕀ℝ n -> ( 𝕀 Double, 𝕀 Double ) ff' i d ts = let df = eqs ts @@ -900,3 +894,86 @@ allNarrowingOperators eps_bis ( Box1Options eps_eq coordsToNarrow eqsToUse ) eqs f' = df `monIndex` linearMonomial i in ( f `index` d, f' `index` d ) {-# INLINEABLE allNarrowingOperators #-} + + +data AdaptiveShavingOptions + = AdaptiveShavingOptions + { Ξ³_init, Οƒ_good, Οƒ_bad, Ξ²_good, Ξ²_bad :: !Double + } + +defaultAdaptiveShavingOptions :: AdaptiveShavingOptions +defaultAdaptiveShavingOptions = + AdaptiveShavingOptions + { Ξ³_init = 0.25 + , Οƒ_good = 0.25 + , Οƒ_bad = 0.75 + , Ξ²_good = 1.5 + , Ξ²_bad = 0.7 + } + +-- | Algorithm @lnar_sbc3ag@ (adaptive guessing) from the paper +-- "Box Consistency through Adaptive Shaving" (Goldsztejn & Goualard, 2010). +leftShave :: Double + -> AdaptiveShavingOptions + -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) ) + -> 𝕀 Double + -> [ 𝕀 Double ] +leftShave Ξ΅ -- minimum width + ( AdaptiveShavingOptions { Ξ³_init, Οƒ_good, Οƒ_bad, Ξ²_good, Ξ²_bad } ) + ff' i0 = + left_narrow Ξ³_init i0 + where + w0 = width i0 + left_narrow :: Double -> 𝕀 Double -> [ 𝕀 Double ] + left_narrow Ξ³ i@( 𝕀 x_inf x_sup ) + -- Stop if the box is too small. + | width i < Ξ΅ + = [ i ] + | otherwise + = go Ξ³ x_sup ( 𝕀 x_inf ( if x_inf == x_sup then x_inf else succFP x_inf ) ) + go :: Double -> Double -> 𝕀 Double -> [ 𝕀 Double ] + go Ξ³ x_sup x_left = + let ( f_x_left, _f'_x_left ) = ff' x_left + in + if 0 `inside` f_x_left + -- Box-consistency achieved; finish. + then [ 𝕀 ( inf x_left ) x_sup ] + else + -- Otherwise, try to shave off a chunk on the left of the interval. + let --x' = 𝕀 ( sup x_left ) x_sup + inf_guess = sup x_left + sup_guess = min x_sup ( inf_guess + Ξ³ * w0 ) -- * width x' ) + guess = 𝕀 inf_guess sup_guess + -- NB: this always uses the initial width (to "avoid asymptotic behaviour" according to the paper) + ( f_guess, f'_guess ) = ff' guess + x_minus_guess = 𝕀 ( min x_sup ( succFP $ sup guess ) ) x_sup + in if not ( 0 `inside` f_guess ) + then + -- We successfully shaved "guess" off; go round again after removing it. + -- TODO: here we could go back to the top with a new "w" maybe? + left_narrow ( Ξ²_good * Ξ³ ) x_minus_guess + else + -- Do a Newton step to try to reduce the guess interval. + -- Starting from "guess", we get a collection of sub-intervals + -- "guesses'" refining where the function can be zero. + let guesses' :: [ ( 𝕀 Double ) ] + guesses' = do + Ξ΄ <- f_x_left `extendedDivide` f'_guess + let guess' = singleton ( inf guess ) - Ξ΄ + map fst $ guess' `intersect` guess + w_guess = width guess + w_guesses' + | null guesses' + = 0 + | otherwise + = sup_guess - minimum ( map inf guesses' ) + Ξ³' | w_guesses' < Οƒ_good * w_guess + -- Good improvement, try larger guesses in the future. + = Ξ²_good * Ξ³ + | w_guesses' > Οƒ_bad * w_guess + -- Poor improvement, try smaller guesses in the future. + = Ξ²_bad * Ξ³ + | otherwise + -- Otherwise, keep the Ξ³ factor the same. + = Ξ³ + in left_narrow Ξ³' =<< ( x_minus_guess : guesses' )