From 0e6b9a822ba17cbc6fe4dd94ae9c25529dacaf7b Mon Sep 17 00:00:00 2001 From: sheaf Date: Sun, 21 Apr 2024 20:35:35 +0200 Subject: [PATCH] Documentation in root isolation module --- brush-strokes/src/lib/Math/Root/Isolation.hs | 149 ++++++++++--------- 1 file changed, 75 insertions(+), 74 deletions(-) diff --git a/brush-strokes/src/lib/Math/Root/Isolation.hs b/brush-strokes/src/lib/Math/Root/Isolation.hs index 89248c2..5a8df0f 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation.hs +++ b/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. - minWidth :: !Double - -- | Given a box and its history, return a round of cusp-finding strategies + { -- | Given a box and its history, return a round of root isolation strategies -- to run, in sequence, on this box. -- - -- Returning an empty list means: give up on this box. - , cuspFindingAlgorithms :: !( Box n -> BoxHistory n -> Either String ( NE.NonEmpty ( RootIsolationAlgorithm n d ) ) ) + -- A return value of @Left whyStop@ indicates stopping any further iterations + -- 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 - = forall r - . [ ( RootIsolationStep, Box n ) ] +-- | A function to choose along which coordination dimension we should +-- perform a bisection. +type BisectionCoordPicker n d + = [ ( RootIsolationStep, Box n ) ] -> BoxHistory n -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) - -> NE.NonEmpty ( Fin n, r ) - -> ( r, String ) + -> forall r. ( NE.NonEmpty ( Fin n, r ) -> ( 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 - , cuspFindingAlgorithms = defaultRootIsolationAlgorithms minWidth narrowAbs + { rootIsolationAlgorithms = defaultRootIsolationAlgorithms minWidth Ξ΅_eq } 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 - -- Give up after depth 50 - -- | length history >= 50 - -- = Left "depth >= 50" + = Left $ "widths <= " ++ show minWidth | 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 - ( RootIsolationOptions - { minWidth - , cuspFindingAlgorithms - } - ) +isolateRootsIn ( RootIsolationOptions { rootIsolationAlgorithms } ) 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 ] ) - Bisection ( BisectionOptions { canHaveSols, fallbackBisectionDim } ) -> do - let ( boxes, ( whatBis, mid ) ) = bisect ( canHaveSols eqs ) ( fallbackBisectionDim roundHistory previousRoundsHistory eqs ) box + return ( Box2Step, [ makeBox2Consistent box2Options eqs box ] ) + Bisection ( BisectionOptions { canHaveSols, fallbackBisectionCoord } ) -> do + 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_mid = projectDims $ eqs boxMid `monIndex` zeroMonomial + f' = fmap ( \ i -> pickEqs $ eqs box `monIndex` linearMonomial i ) ( universe @n ) + 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 -forEachDim a0 f = go ( toList $ universe @n ) a0 +forEachCoord :: forall n a m. ( KnownNat n, Monad m ) => a -> ( Fin n -> a -> m a ) -> m a +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