Implement adaptive shaving algorithm

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

@@ -110,8 +110,7 @@ main = do
   putStrLn $ unlines
     [ replicate 40 '='
     , "Best time: " ++ show bestTime ++ "s"
-    , "Best parameters:"
+    , "Best parameters: " ++ show (testName bestTest)
-    --, show (testCuspOptions bestTest)
     ]
 
 benchTestCase :: TestCase -> IO Double
@@ -181,9 +180,15 @@ testCases = benchCases
 
 benchCases :: [ TestCase ]
 benchCases =
-  [ ellipseTestCase opts "full" ( 0, 1 ) pi $ defaultStartBoxes [ 2 ] ] -- [ 0 .. 3 ] ]
+  [ ellipseTestCase opts ("minWidth=" ++ show minWidth ++ ",ε=" ++ show narrowAbs) ( 0, 1 ) pi $ defaultStartBoxes [ 2 ]
-  where
+  | minWidth <- [ 1e-5 ]
-    opts = defaultRootIsolationOptions
+  , 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
+            }
+  ]
 
 --------------------------------------------------------------------------------
 

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
+                     else case filter ( not . isTooSmall ) $ sortOnArg ( Dual . spread ) datPerCoord of
-                       [] -> error "dimension 0"
+                       [] -> ( coordIndex d1, "tooWide'" )
-                       Arg _ d : _ ->
+                       Arg _ d : _ -> ( coordIndex d, "spread" )
-                        (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
+    doLoop λ x = do
-      x'' <- forEachDim @n x $ \ i ->
+      x'' <- forEachDim @n x $ boundConsistency λ
-        boundConsistency ( `index` i ) ( set i ) lambda
       modified <- State.get
-      let lambda' = if modified then lambda else 0.5 * lambda
+      let λ' = if modified then λ else 0.5 * λ
-      if lambda' < lambdaMin
+      if λ' < λMin
       then return x''
-      else do { State.put False ; doLoop lambda' x'' }
+      else do { State.put False ; doLoop λ' x'' }
 
-    boundConsistency :: ( Box n -> 𝕀 Double )
+    boundConsistency :: Double -> Fin n -> Box n -> State Bool ( Box n )
-                     -> ( 𝕀 Double -> Box n -> Box n )
+    boundConsistency λ i box = do
-                     -> Double -> Box n -> State Bool ( Box n )
-    boundConsistency getter setter lambda box = do
       let x@( 𝕀 x_inf x_sup ) = getter box
-          c1 = ( 1 - lambda ) * x_inf + lambda * x_sup
+          c1 = ( 1 - λ ) * x_inf + λ * x_sup
-          c2 = lambda * x_inf + ( 1 - lambda ) * 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
+              x's = do δ <- f_x_left `extendedDivide` f'_x
-                       let x_new = x_left - delta
+                       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
+               | ( width x - sum ( map width x's ) ) < ε_eq
-               -- TODO: do a check with the relative reduction in size?
                -> 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
+              x's = do δ <- f_x_right `extendedDivide` f'_x
-                       let x_new = x_right - delta
+                       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
+               | ( width x - sum ( map width x's ) ) < ε_eq
-               -- TODO: do a check with the relative reduction in size?
                -> 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 ]
+  | narrowFn <- [ leftNarrow ε_eq ε_bis, rightNarrow ε_eq ε_bis ]
-  , ( coordIndex, fn ) <-
+  , coordIndex <- coordsToNarrow
-      [ ( i, ff' i d )
+  , eqnIndex <- eqsToUse
-      | i <- coordsToNarrow
-      , d <- 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' )