Documentation in root isolation module

brush-strokes/src/lib/Math/Root/Isolation.hs

@@ -97,12 +97,12 @@ import qualified Math.Ring as Ring
 
 --------------------------------------------------------------------------------
 
--- | A tree recording the steps taken when doing cusp finding.
+-- | A tree recording the steps taken when isolating roots.
 data RootIsolationTree d
   = RootIsolationLeaf String d
   | RootIsolationStep RootIsolationStep [(d, RootIsolationTree d)]
 
--- | A description of a step taken in cusp-finding.
+-- | A description of a step taken when isolating roots.
 data RootIsolationStep
   = GaussSeidelStep
   | BisectionStep String Double
@@ -155,17 +155,19 @@ type BoxCt n d = ( n ~ N, d ~ 3 )
   )
 -}
 
--- | Options for the root isolation methods in 'findCusps' and 'isolateRootsIn'.
+-- | Options for the root isolation methods in 'isolateRootsIn'.
 type RootIsolationOptions :: Nat -> Nat -> Type
-data RootIsolationOptions n d
+newtype RootIsolationOptions n d
   = RootIsolationOptions
-  { -- | Minimum width of a box.
+  { -- | Given a box and its history, return a round of root isolation strategies
-    minWidth :: !Double
-    -- | Given a box and its history, return a round of cusp-finding strategies
     -- to run, in sequence, on this box.
     --
-    -- Returning an empty list means: give up on this box.
+    -- A return value of @Left whyStop@ indicates stopping any further iterations
-  , cuspFindingAlgorithms :: !( Box n -> BoxHistory n -> Either String ( NE.NonEmpty ( RootIsolationAlgorithm n d ) ) )
+    -- on this box, e.g. because it is too small.
+   rootIsolationAlgorithms
+     :: Box n
+     -> BoxHistory n
+     -> Either String ( NE.NonEmpty ( RootIsolationAlgorithm n d ) )
   }
 
 type RootIsolationAlgorithm :: Nat -> Nat -> Type
@@ -184,16 +186,16 @@ type BisectionOptions :: Nat -> Nat -> Type
 data BisectionOptions n d =
   BisectionOptions
   { canHaveSols :: !( ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) -> Box n -> Bool )
-  , fallbackBisectionDim :: !( BisectionDimPicker n d )
+  , fallbackBisectionCoord :: !( BisectionCoordPicker n d )
   }
 
-type BisectionDimPicker n d
+-- | A function to choose along which coordination dimension we should
-  =  forall r
+-- perform a bisection.
-  .  [ ( RootIsolationStep, Box n ) ]
+type BisectionCoordPicker n d
+  =  [ ( RootIsolationStep, Box n ) ]
   -> BoxHistory n
   -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
-  -> NE.NonEmpty ( Fin n, r )
+  -> forall r. ( NE.NonEmpty ( Fin n, r ) -> ( r, String ) )
-  -> ( r, String )
 
 -- | Options for the interval Gauss–Seidel method.
 type GaussSeidelOptions :: Nat -> Nat -> Type
@@ -204,12 +206,13 @@ data GaussSeidelOptions n d =
       -- | Function that projects over the equations we will consider
       -- (the identity for a well-determined problem, or a projection for
       -- an overdetermined system).
-    , gsDims :: ( 𝕀ℝ d -> 𝕀ℝ n ) }
+    , gsPickEqs :: ( 𝕀ℝ d -> 𝕀ℝ n ) }
 
 -- | Options for the @box(1)@-consistency method.
 data Box1Options n d =
   Box1Options
     { box1EpsEq :: !Double
+    , box1EpsBis :: !Double
     , box1CoordsToNarrow :: [ Fin n ]
     , box1EqsToUse :: [ Fin d ]
     }
@@ -218,6 +221,7 @@ data Box1Options n d =
 data Box2Options n d =
   Box2Options
     { box2EpsEq :: !Double
+    , box2EpsBis :: !Double
     , box2LambdaMin :: !Double
     , box2CoordsToNarrow :: [ Fin n ]
     , box2EqsToUse :: [ Fin d ]
@@ -226,24 +230,20 @@ data Box2Options n d =
 defaultRootIsolationOptions :: BoxCt n d => RootIsolationOptions n d
 defaultRootIsolationOptions =
   RootIsolationOptions
-    { minWidth
+    { rootIsolationAlgorithms = defaultRootIsolationAlgorithms minWidth ε_eq
-    , cuspFindingAlgorithms = defaultRootIsolationAlgorithms minWidth narrowAbs
     }
   where
     minWidth = 1e-5
-    narrowAbs = 5e-3
+    ε_eq = 5e-3
 
 defaultRootIsolationAlgorithms
   :: forall n d
   .  BoxCt n d
   => Double -> Double -> Box n -> BoxHistory n -> Either String ( NE.NonEmpty ( RootIsolationAlgorithm n d ) )
-defaultRootIsolationAlgorithms minWidth narrowAbs box history
+defaultRootIsolationAlgorithms minWidth ε_eq box history
   -- All the box widths are lower than the minimum width: give up.
   | and $ ( \ cd -> width cd <= minWidth ) <$> coordinates box
-  = Left $ "sizes <= " ++ show minWidth
+  = Left $ "widths <= " ++ show minWidth
-  -- Give up after depth 50
-  -- | length history >= 50
-  -- = Left "depth >= 50"
   | otherwise
   = Right $ case history of
       lastRoundBoxes : _
@@ -253,7 +253,7 @@ defaultRootIsolationAlgorithms minWidth narrowAbs box history
         -- ...and the last round didn't sufficiently reduce the size of the box...
         , not $ box `sufficientlySmallerThan` lastRoundFirstBox
         -- ...then try bisecting the box.
-        -> NE.singleton $ Bisection ( defaultBisectionOptions minWidth narrowAbs box )
+        -> NE.singleton $ Bisection ( defaultBisectionOptions minWidth ε_eq box )
       -- Otherwise, do a normal round.
       -- Currently: we try an interval Gauss–Seidel step followed by box(1)-consistency.
       _ -> GaussSeidel _gaussSeidelOptions
@@ -263,7 +263,8 @@ defaultRootIsolationAlgorithms minWidth narrowAbs box history
     _box1Options :: Box1Options n d
     _box1Options =
       Box1Options
-        { box1EpsEq = narrowAbs
+        { box1EpsEq = ε_eq
+        , box1EpsBis = minWidth
         , box1CoordsToNarrow = toList $ universe @n
         , box1EqsToUse = toList $ universe @d
         }
@@ -271,7 +272,8 @@ defaultRootIsolationAlgorithms minWidth narrowAbs box history
     _box2Options :: Box2Options n d
     _box2Options =
       Box2Options
-        { box2EpsEq = narrowAbs
+        { box2EpsEq = ε_eq
+        , box2EpsBis = minWidth
         , box2LambdaMin = 0.001
         , box2CoordsToNarrow = toList $ universe @n
         , box2EqsToUse = toList $ universe @d
@@ -281,7 +283,7 @@ defaultRootIsolationAlgorithms minWidth narrowAbs box history
     _gaussSeidelOptions =
       GaussSeidelOptions
         { gsPreconditioner = InverseMidJacobian
-        , gsDims =
+        , gsPickEqs =
             \ ( 𝕀 ( ℝ3 a_lo b_lo c_lo ) ( ℝ3 a_hi b_hi c_hi ) ) ->
               case length history `mod` 3 of
                 0 -> 𝕀 ( ℝ2 a_lo b_lo ) ( ℝ2 a_hi b_hi )
@@ -289,11 +291,12 @@ defaultRootIsolationAlgorithms minWidth narrowAbs box history
                 _ -> 𝕀 ( ℝ2 a_lo c_lo ) ( ℝ2 a_hi c_hi )
         }
 
-    -- Did we reduce the box width by at least "narrowAbs" in at least one of the dimensions?
+    -- Did we reduce the box width by at least ε_eq
+    -- in at least one of the coordinates?
     sufficientlySmallerThan :: Box n -> Box n -> Bool
     b1 `sufficientlySmallerThan` b2 =
       or $
-        ( \ cd1 cd2 -> width cd2 - width cd1 >= narrowAbs )
+        ( \ cd1 cd2 -> width cd2 - width cd1 >= ε_eq )
           <$> coordinates b1
           <*> coordinates b2
 {-# INLINEABLE defaultRootIsolationAlgorithms #-}
@@ -303,7 +306,7 @@ defaultBisectionOptions
   .  ( 1 <= n, BoxCt n d )
   => Double -> Double
   -> Box n -> BisectionOptions n d
-defaultBisectionOptions minWidth _narrowAbs box =
+defaultBisectionOptions minWidth _ε_eq box =
   BisectionOptions
     { canHaveSols =
         \ eqs box' ->
@@ -313,16 +316,16 @@ defaultBisectionOptions minWidth _narrowAbs box =
           in unT ( origin @Double ) `inside` iRange'
 
           -- box(1)-consistency
-          --let box1Options = Box1Options _narrowAbs ( toList $ universe @n ) ( toList $ universe @d )
+          --let box1Options = Box1Options _ε_eq ( toList $ universe @n ) ( toList $ universe @d )
           --in not $ null $ makeBox1Consistent _minWidth box1Options eqs box'
 
           -- box(2)-consistency
-          --let box2Options = Box2Options _narrowAbs 0.001 ( toList $ universe @n ) ( toList $ universe @d )
+          --let box2Options = Box2Options _ε_eq 0.001 ( toList $ universe @n ) ( toList $ universe @d )
           --    box'' = makeBox2Consistent _minWidth box2Options eqs box'
           --    iRange'' :: Box d
           --    iRange'' = eqs box'' `monIndex` zeroMonomial
           --in unT ( origin @Double ) `inside` iRange''
-    , fallbackBisectionDim =
+    , fallbackBisectionCoord =
         \ _roundHist _prevRoundsHist eqs possibleCoordChoices ->
           let datPerCoord =
                 possibleCoordChoices <&> \ ( i, r ) ->
@@ -334,14 +337,15 @@ defaultBisectionOptions minWidth _narrowAbs box =
                     }
 
               -- First, check if the largest dimension is over 20 times larger
-              -- than the smallest dimension; if so bisect along that coordinate.
+              -- than the smallest dimension; if so bisect that dimension.
           in case sortOnArgNE ( width . coordInterval ) datPerCoord of
                 Arg _ d NE.:| [] -> ( coordBisectionData d, "" )
                 Arg w0 _ NE.:| ds ->
                   let Arg w1 d1 = last ds
                   in if w1 >= 20 * w0
                      then ( coordBisectionData d1, "tooWide" )
-                     -- Otherwise, pick the dimension with the largest spread = width * Jacobian column norm
+                     -- Otherwise, pick the coordination dimension with maximum spread
+                     -- (spread = width * Jacobian column norm).
                      else
                       let isTooSmall ( Arg ( Dual w ) _ ) = w < minWidth
                       in case NE.filter ( not . isTooSmall ) $ sortOnArgNE ( Dual . spread ) datPerCoord of
@@ -391,12 +395,7 @@ isolateRootsIn
   -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
   -> Box n
   -> ( [ ( Box n, RootIsolationTree ( Box n ) ) ], ( [ Box n ], [ Box n ] ) )
-isolateRootsIn
+isolateRootsIn ( RootIsolationOptions { rootIsolationAlgorithms } )
-  ( RootIsolationOptions
-    { minWidth
-    , cuspFindingAlgorithms
-    }
-  )
   eqs initBox = runWriter $ map ( initBox , ) <$> go [] initBox
 
   where
@@ -411,7 +410,7 @@ isolateRootsIn
         -- Box doesn't contain a solution: discard it.
       = return []
       | otherwise
-      = case cuspFindingAlgorithms cand history of
+      = case rootIsolationAlgorithms cand history of
           Right strats -> doStrategies history strats cand
           Left whyStop -> do
             -- We are giving up on this box (e.g. because it is too small).
@@ -421,7 +420,7 @@ isolateRootsIn
         iRange :: Box d
         iRange = eqs cand `monIndex` zeroMonomial
 
-    -- Run a round of cusp finding strategies, then recur.
+    -- Run a round of root isolation strategies, then recur.
     doStrategies
       :: BoxHistory n
       -> NonEmpty ( RootIsolationAlgorithm n d )
@@ -437,7 +436,7 @@ isolateRootsIn
                             [ RootIsolationTree ( Box n )]
         do_strats thisRoundHist ( algo NE.:| algos ) box = do
           -- Run one strategy in the round.
-          ( step, boxes ) <- doStrategy thisRoundHist prevRoundsHist eqs minWidth algo box
+          ( step, boxes ) <- doStrategy thisRoundHist prevRoundsHist eqs algo box
           case algos of
             -- If there are other algorithms to run in this round, run them next.
             nextAlgo : otherAlgos ->
@@ -459,35 +458,39 @@ isolateRootsIn
       return [ RootIsolationStep step $ concat rest ]
 {-# INLINEABLE isolateRootsIn #-}
 
--- | Execute a cusp-finding strategy, replacing the input box with
+-- | Execute a root isolation strategy, replacing the input box with
 -- (hopefully smaller) output boxes.
 doStrategy
   :: BoxCt n d
   => [ ( RootIsolationStep, Box n ) ]
   -> BoxHistory n
   -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
-  -> Double
   -> RootIsolationAlgorithm n d
   -> Box n
   -> Writer ( [ Box n ], [ Box n ] )
             ( RootIsolationStep, [ Box n ] )
-doStrategy roundHistory previousRoundsHistory eqs minWidth algo box =
+doStrategy roundHistory previousRoundsHistory eqs algo box =
   case algo of
     GaussSeidel gsOptions -> do
       boxes <- intervalGaussSeidel gsOptions eqs box
       return ( GaussSeidelStep, boxes )
     Box1 box1Options ->
-      return ( Box1Step, makeBox1Consistent minWidth box1Options eqs box )
+      return ( Box1Step, makeBox1Consistent box1Options eqs box )
     Box2 box2Options ->
-      return ( Box2Step, [ makeBox2Consistent minWidth box2Options eqs box ] )
+      return ( Box2Step, [ makeBox2Consistent box2Options eqs box ] )
-    Bisection ( BisectionOptions { canHaveSols, fallbackBisectionDim } ) -> do
+    Bisection ( BisectionOptions { canHaveSols, fallbackBisectionCoord } ) -> do
-      let ( boxes, ( whatBis, mid ) ) = bisect ( canHaveSols eqs ) ( fallbackBisectionDim roundHistory previousRoundsHistory eqs ) box
+      let ( boxes, ( whatBis, mid ) ) =
+            bisect
+              ( canHaveSols eqs )
+              ( fallbackBisectionCoord roundHistory previousRoundsHistory eqs )
+              box
       return ( BisectionStep whatBis mid, boxes )
 {-# INLINEABLE doStrategy #-}
 
 -- | Bisect the given box.
 --
--- (The difficult part lies in determining along which dimension to bisect.)
+-- (The difficult part lies in determining along which coordinate
+-- dimension to bisect.)
 bisect
   :: forall n
   .  ( 1 <= n, KnownNat n, Representable Double ( ℝ n ) )
@@ -564,14 +567,13 @@ intervalGaussSeidel
      -- ^ box
   -> Writer ( [ 𝕀ℝ n ], [ 𝕀ℝ n ] ) [ 𝕀ℝ n ]
 intervalGaussSeidel
-  ( GaussSeidelOptions { gsPreconditioner = precondMeth, gsDims = projectDims } )
+  ( GaussSeidelOptions { gsPreconditioner = precondMeth, gsPickEqs = pickEqs } )
   eqs
   box
   | let boxMid = singleton $ midpoint box
-        df = eqs box
         f' :: Vec n ( 𝕀ℝ n )
-        f' = fmap ( \ i -> projectDims $ df `monIndex` linearMonomial i ) ( universe @n )
+        f' = fmap ( \ i -> pickEqs $ eqs box `monIndex` linearMonomial i ) ( universe @n )
-        f_mid = projectDims $ eqs boxMid `monIndex` zeroMonomial
+        f_mid = pickEqs $ eqs boxMid `monIndex` zeroMonomial
 
   = let -- Interval Newton method: take one Gauss–Seidel step
         -- for the system of equations f'(x) ( x - x_mid ) = - f(x_mid).
@@ -622,7 +624,7 @@ gaussSeidelStep
   -> T ( 𝕀ℝ n )     -- ^ initial box \( X \)
   -> [ ( T ( 𝕀ℝ n ), Bool ) ]
 gaussSeidelStep as b ( T x0 ) = coerce $
-  forEachDim @n ( x0, True ) $ \ i ( x, contraction ) -> do
+  forEachCoord @n ( x0, True ) $ \ i ( x, contraction ) -> do
     -- x_i' = ( b_i - sum { j /= i } a_ij * x_j ) / a_ii
     x_i'0 <- ( b `index` i - sum [ ( as ! j ) `index` i * x `index` j | j <- toList ( universe @n ), j /= i ] )
           `extendedDivide` ( ( as ! i ) `index` i )
@@ -635,14 +637,14 @@ gaussSeidelStep as b ( T x0 ) = coerce $
 
 -- | Run an effectful computation several times in sequence, piping its output
 -- into the next input, once for each coordinate dimension.
-forEachDim :: forall n a m. ( KnownNat n, Monad m ) => a -> ( Fin n -> a -> m a ) -> m a
+forEachCoord :: forall n a m. ( KnownNat n, Monad m ) => a -> ( Fin n -> a -> m a ) -> m a
-forEachDim a0 f = go ( toList $ universe @n ) a0
+forEachCoord a0 f = go ( toList $ universe @n ) a0
   where
     go [] a = return a
     go ( i : is ) a = do
       a' <- f i a
       go is a'
-{-# INLINEABLE forEachDim #-}
+{-# INLINEABLE forEachCoord #-}
 
 -- | Compute the intersection of two intervals.
 --
@@ -671,32 +673,31 @@ data Preconditioner
 -- "Presentation of a highly tuned multithreaded interval solver for underdetermined and well-determined nonlinear systems"
 makeBox1Consistent
   :: BoxCt n d
-  => Double -> Box1Options n d
+  => Box1Options n d
   -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
   -> Box n -> [ Box n ]
-makeBox1Consistent minWidth box1Options eqs x =
+makeBox1Consistent box1Options eqs x =
   ( `State.evalState` False ) $
-    pipeFunctionsWhileTrue ( allNarrowingOperators minWidth box1Options eqs ) x
+    pipeFunctionsWhileTrue ( allNarrowingOperators box1Options 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
   :: forall n d
   .  BoxCt n d
-  => Double
+  => Box2Options n d
-  -> Box2Options n d
   -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
   -> Box n -> Box n
-makeBox2Consistent minWidth (Box2Options ε_eq λMin coordsToNarrow eqsToUse) eqs x0
+makeBox2Consistent (Box2Options ε_eq ε_bis λMin coordsToNarrow eqsToUse) eqs x0
   = ( `State.evalState` False ) $ doLoop 0.25 x0
   where
     box1Options :: Box1Options n d
-    box1Options = Box1Options ε_eq coordsToNarrow eqsToUse
+    box1Options = Box1Options ε_eq ε_bis coordsToNarrow eqsToUse
     doBox1 :: Box n -> [ Box n ]
-    doBox1 = makeBox1Consistent minWidth box1Options eqs
+    doBox1 = makeBox1Consistent box1Options eqs
     doLoop :: Double -> Box n -> State Bool ( Box n )
     doLoop λ x = do
-      x'' <- forEachDim @n x $ boundConsistency λ
+      x'' <- forEachCoord @n x $ boundConsistency λ
       modified <- State.get
       let λ' = if modified then λ else 0.5 * λ
       if λ' < λMin
@@ -741,6 +742,7 @@ precondition meth jac_mid as b =
       | let mat = toEigen jac_mid
             det = Eigen.determinant mat
       , not $ nearZero det
+      -- (TODO: a bit wasteful to compute determinant then inverse.)
       , let precond = Eigen.inverse mat
             doPrecond = matMulVec ( fromEigen precond )
       -> ( fmap doPrecond as, doPrecond b )
@@ -788,7 +790,7 @@ maxVI ( 𝕀 los his ) =
     where
       maxAbs :: Double -> Double -> Double
       maxAbs lo hi = max ( abs lo ) ( abs hi )
-{-# INLINEABLE maxVI  #-}
+{-# INLINEABLE maxVI #-}
 
 -- Use the univariate interval Newton method to narrow from the left
 -- a candidate interval.
@@ -889,11 +891,10 @@ pipeFunctionsWhileTrue fns = go fns
 allNarrowingOperators
   :: forall n d
   .  BoxCt n d
-  => Double
+  => Box1Options n d
-  -> Box1Options n d
   -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
   -> [ Box n -> State Bool [ Box n ] ]
-allNarrowingOperators ε_bis ( Box1Options ε_eq coordsToNarrow eqsToUse ) eqs =
+allNarrowingOperators ( Box1Options ε_eq ε_bis coordsToNarrow eqsToUse ) eqs =
   [ \ cand ->
     let getter = ( `index` coordIndex )
         setter = set coordIndex