diff --git a/MetaBrush.cabal b/MetaBrush.cabal index 641b871..adaed84 100644 --- a/MetaBrush.cabal +++ b/MetaBrush.cabal @@ -81,8 +81,9 @@ common common StandaloneDeriving StandaloneKindSignatures TupleSections + TypeAbstractions TypeApplications - TypeFamilies + TypeFamilyDependencies TypeOperators UnboxedTuples ViewPatterns diff --git a/brush-strokes/brush-strokes.cabal b/brush-strokes/brush-strokes.cabal index cf5afde..1baf32c 100644 --- a/brush-strokes/brush-strokes.cabal +++ b/brush-strokes/brush-strokes.cabal @@ -52,8 +52,9 @@ common common StandaloneDeriving StandaloneKindSignatures TupleSections + TypeAbstractions TypeApplications - TypeFamilies + TypeFamilyDependencies TypeOperators UnboxedTuples ViewPatterns @@ -125,6 +126,10 @@ library , Math.Ring , Math.Roots , Math.Root.Isolation + , Math.Root.Isolation.Bisection + , Math.Root.Isolation.GaussSeidel + , Math.Root.Isolation.Narrowing + , Math.Root.Isolation.Core , Debug.Utils other-modules: diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs index c408119..f3f4f60 100644 --- a/brush-strokes/src/cusps/bench/Main.hs +++ b/brush-strokes/src/cusps/bench/Main.hs @@ -121,7 +121,9 @@ benchTestCase ( TestCase { testName, testBrushStroke, testCuspOptions, testStart before <- getMonotonicTime let ( _, testStrokeFnI ) = brushStrokeFunctions testBrushStroke ( dunno, sols ) = - foldMap ( \ ( i, ( _trees, ( mbCusps, defCusps ) ) ) -> ( map ( i , ) mbCusps, map ( i, ) defCusps ) ) $ + foldMap + ( \ ( i, ( _trees, DoneBoxes { doneSolBoxes = defCusps, doneGiveUpBoxes = mbCusps } ) ) -> + ( map ( ( i , ) . snd ) mbCusps, map ( i, ) defCusps ) ) $ IntMap.toList $ findCuspsIn testCuspOptions testStrokeFnI $ IntMap.fromList diff --git a/brush-strokes/src/lib/Math/Algebra/Dual.hs b/brush-strokes/src/lib/Math/Algebra/Dual.hs index fc1f342..470538f 100644 --- a/brush-strokes/src/lib/Math/Algebra/Dual.hs +++ b/brush-strokes/src/lib/Math/Algebra/Dual.hs @@ -51,7 +51,7 @@ deriving stock instance Functor ( D k u ) => Functor ( C k u ) -- | @D k u v@ is the space of @k@-th order germs of functions from @u@ to @v@, -- represented by the algebra: -- --- \[ \mathbb{Z}[x_1, \ldots, x_n]/(x_1, \ldots, x_n)^(k+1) \otimes_\mathbb{Z} v \] +-- \[ \mathbb{Z}[x_1, \ldots, x_n]/(x_1, \ldots, x_n)^{k+1} \otimes_\mathbb{Z} v \] -- -- when @u@ is of dimension @n@. type D :: Nat -> Type -> Type -> Type diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke.hs b/brush-strokes/src/lib/Math/Bezier/Stroke.hs index 737a371..e06591d 100644 --- a/brush-strokes/src/lib/Math/Bezier/Stroke.hs +++ b/brush-strokes/src/lib/Math/Bezier/Stroke.hs @@ -597,8 +597,8 @@ outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams case mbCuspOptions of Just cuspOptions -> foldMap - ( \ ( i, ( _trees, ( potCusps, defCusps ) ) ) -> - ( map ( i , ) potCusps, map ( i , ) defCusps ) + ( \ ( i, ( _trees, DoneBoxes { doneSolBoxes = defCusps, doneGiveUpBoxes = potCusps } ) ) -> + ( map ( ( i , ) . fst ) potCusps, map ( i , ) defCusps ) ) ( IntMap.toList $ findCusps cuspOptions curvesI ) Nothing -> @@ -1166,7 +1166,7 @@ cuspCoords eqs ( i, box ) findCusps :: RootIsolationOptions N 3 -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) - -> IntMap ( [ ( Box N, RootIsolationTree ( Box N ) ) ], ( [ Box N ], [ Box N ] ) ) + -> IntMap ( [ ( Box N, RootIsolationTree ( Box N ) ) ], DoneBoxes N ) findCusps opts boxStrokeData = findCuspsIn opts boxStrokeData $ IntMap.fromList @@ -1186,7 +1186,7 @@ findCuspsIn :: RootIsolationOptions N 3 -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) -> IntMap [ Box 2 ] - -> IntMap ( [ ( Box N, RootIsolationTree ( Box N ) ) ], ( [ Box N ], [ Box N ] ) ) + -> IntMap ( [ ( Box N, RootIsolationTree ( Box N ) ) ], DoneBoxes N ) findCuspsIn opts boxStrokeData initBoxes = IntMap.mapWithKey ( \ i -> foldMap ( isolateRootsIn opts ( eqnPiece i ) ) ) initBoxes where diff --git a/brush-strokes/src/lib/Math/Interval.hs b/brush-strokes/src/lib/Math/Interval.hs index 39d7c0d..cb80187 100644 --- a/brush-strokes/src/lib/Math/Interval.hs +++ b/brush-strokes/src/lib/Math/Interval.hs @@ -6,13 +6,15 @@ module Math.Interval ( 𝕀(𝕀), inf, sup - , width + , width, midpoint , scaleInterval , 𝕀ℝ + , isCanonical , singleton, nonDecreasing , inside , aabb , extendedDivide, extendedRecip + , intersect, bisect ) where @@ -20,6 +22,7 @@ module Math.Interval import Prelude hiding ( Num(..), Fractional(..) ) import Data.List ( nub ) +import qualified Data.List.NonEmpty as NE -- acts import Data.Act @@ -35,6 +38,8 @@ import Data.Group.Generics -- brush-strokes import Math.Algebra.Dual +import Math.Float.Utils + ( prevFP ) import Math.Interval.Internal ( 𝕀(𝕀), inf, sup, scaleInterval ) import Math.Linear @@ -51,18 +56,62 @@ import Math.Ring type 𝕀ℝ n = 𝕀 ( ℝ n ) type instance D k ( 𝕀 v ) = D k v +-- | An interval reduced to a single point. singleton :: a -> 𝕀 a singleton a = 𝕀 a a +-- | The width of an interval. +-- +-- NB: assumes the endpoints are neither @NaN@ nor infinity. width :: AbelianGroup a => 𝕀 a -> a width ( 𝕀 lo hi ) = hi - lo +{-# INLINEABLE width #-} -- | Turn a non-decreasing function into a function on intervals. nonDecreasing :: ( a -> b ) -> 𝕀 a -> 𝕀 b nonDecreasing f ( 𝕀 lo hi ) = 𝕀 ( f lo ) ( f hi ) +-- | Does the given value lie inside the specified interval? inside :: Ord a => a -> 𝕀 a -> Bool inside x ( 𝕀 lo hi ) = x >= lo && x <= hi +{-# INLINEABLE inside #-} + +-- | Is this interval canonical, i.e. it consists of either 1 or 2 floating point +-- values only? +isCanonical :: 𝕀 Double -> Bool +isCanonical ( 𝕀 x_inf x_sup ) = x_inf >= prevFP x_sup + +-- | The midpoint of an interval. +-- +-- NB: assumes the endpoints are neither @NaN@ nor infinity. +midpoint :: 𝕀 Double -> Double +midpoint ( 𝕀 x_inf x_sup ) = 0.5 * ( x_inf + x_sup ) + +-- | Compute the intersection of two intervals. +-- +-- Returns whether the first interval is a strict subset of the second interval +-- (or the intersection is a single point). +intersect :: 𝕀 Double -> 𝕀 Double -> [ ( 𝕀 Double, Bool ) ] +intersect ( 𝕀 lo1 hi1 ) ( 𝕀 lo2 hi2 ) + | lo > hi + = [ ] + | otherwise + = [ ( 𝕀 lo hi, ( lo1 > lo2 && hi1 < hi2 ) || ( lo == hi ) ) ] + where + lo = max lo1 lo2 + hi = min hi1 hi2 + +-- | Bisect an interval. +-- +-- Generally returns two sub-intervals, but returns the input interval if +-- it is canonical (and thus bisection would be pointless). +bisect :: 𝕀 Double -> NE.NonEmpty ( 𝕀 Double ) +bisect x@( 𝕀 x_inf x_sup ) + | isCanonical x + = NE.singleton x + | otherwise + = 𝕀 x_inf x_mid NE.:| [ 𝕀 x_mid x_sup ] + where x_mid = midpoint x deriving via ViaAbelianGroup ( T ( 𝕀 Double ) ) instance Semigroup ( T ( 𝕀 Double ) ) diff --git a/brush-strokes/src/lib/Math/Root/Isolation.hs b/brush-strokes/src/lib/Math/Root/Isolation.hs index 57a3a9a..579fbc5 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation.hs @@ -1,169 +1,89 @@ {-# LANGUAGE ScopedTypeVariables #-} module Math.Root.Isolation - ( -- * Root isolation using interval arithmetic - Box, BoxHistory - , isolateRootsIn + ( + -- * Main entry point: execute root isolation strategies + isolateRootsIn + , Box + , DoneBoxes(..) - -- ** Configuration of the root isolation methods - , RootIsolationOptions(..), defaultRootIsolationOptions - , RootIsolationAlgorithm(..), defaultRootIsolationAlgorithms + -- * General API for root isolation algorithms + , BoxCt, BoxHistory + , RootIsolationAlgorithm(..) + , RootIsolationAlgorithmWithOptions(..) - -- *** Options for the bisection method - , BisectionOptions(..) - - -- *** Options for the interval Newton method with Gauss–Seidel step - , GaussSeidelOptions(..), Preconditioner(..) - - -- *** Options for the @box(1)@-consistency algorithm - , Box1Options(..) - - -- *** Options for the @box(2)@-consistency algorithm - , Box2Options(..) - - -- ** Trees recording search space of root isolation algorithms - , RootIsolationTree(..), showRootIsolationTree + -- ** Inspecting history , RootIsolationStep(..) + + -- ** Visualising history + , RootIsolationTree(..) + , showRootIsolationTree + + -- * Configuration of the root isolation methods + , RootIsolationOptions(..), defaultRootIsolationOptions + , defaultRootIsolationAlgorithms + + -- * The bisection method + , Bisection + , BisectionOptions(..), BisectionCoordPicker + , defaultBisectionOptions + + -- * The interval Newton method with Gauss–Seidel step + , GaussSeidel + , GaussSeidelOptions(..), Preconditioner(..) + , defaultGaussSeidelOptions + + -- * Box-consistency methods + + -- ** The @box(1)@-consistency algorithm + , Box1 + , Box1Options(..) + , defaultBox1Options + + -- *** Narrowing operators for the @box(1)@-consistency algorithm + , NarrowingMethod(..) + , narrowingMethods + , AdaptiveShavingOptions(..) + , defaultAdaptiveShavingOptions + + -- ** The @box(2)@-consistency algorithm + , Box2 + , Box2Options(..) + , defaultBox2Options ) where -- base -import Control.Arrow - ( first ) -import Control.Monad - ( when ) -import Data.Bifunctor - ( Bifunctor(bimap) ) -import Data.Coerce - ( coerce ) import Data.Kind ( Type ) -import Data.Foldable - ( toList ) -import Data.Functor - ( (<&>) ) -import Data.List - ( partition ) import Data.List.NonEmpty ( NonEmpty ) import qualified Data.List.NonEmpty as NE - ( NonEmpty(..), cons, filter, fromList, last, singleton, sort ) -import Data.Proxy - ( Proxy(..) ) -import Data.Semigroup - ( Arg(..), Dual(..) ) -import Data.Type.Ord - ( OrderingI(..) ) -import Numeric - ( showFFloat ) + ( NonEmpty(..), last, singleton ) import GHC.TypeNats - ( Nat, KnownNat - , type (<=) - , cmpNat - ) - --- containers -import Data.Tree - ( Tree(..) ) - --- eigen -import qualified Eigen.Matrix as Eigen - ( Matrix - , determinant, generate, inverse, unsafeCoeff - ) + ( Nat ) -- transformers -import Control.Monad.Trans.State.Strict as State - ( State, evalState, get, put ) import Control.Monad.Trans.Writer.CPS ( Writer, runWriter, tell ) -- MetaBrush import Math.Algebra.Dual ( D ) -import Math.Epsilon - ( nearZero ) -import Math.Float.Utils - ( succFP, prevFP ) import Math.Interval import Math.Linear import Math.Module ( Module(..) ) import Math.Monomial - ( MonomialBasis(..), Deg, Vars - , linearMonomial, zeroMonomial - ) -import qualified Math.Ring as Ring + ( MonomialBasis(..), zeroMonomial ) ---import Debug.Utils +import Math.Root.Isolation.Bisection +import Math.Root.Isolation.Core +import Math.Root.Isolation.GaussSeidel +import Math.Root.Isolation.Narrowing -------------------------------------------------------------------------------- --- | 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 when isolating roots. -data RootIsolationStep - = GaussSeidelStep - | BisectionStep String Double - | Box1Step - | Box2Step - -instance Show RootIsolationStep where - showsPrec _ GaussSeidelStep = showString "GS" - showsPrec _ ( BisectionStep coord w ) - = showString "bis " - . showParen True - ( showString coord - . showString " = " - . showsPrec 0 w - ) - showsPrec _ Box1Step = showString "box(1)" - showsPrec _ Box2Step = showString "box(2)" - -showRootIsolationTree - :: ( Representable Double ( ℝ n ), Show ( Box n ) ) - => Box n -> RootIsolationTree ( Box n ) -> Tree String -showRootIsolationTree cand (RootIsolationLeaf why l) = Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ why ++ " " ++ show l) [] -showRootIsolationTree cand (RootIsolationStep s ts) - = Node (show cand ++ " abc " ++ showArea (boxArea cand) ++ " " ++ show s) $ map (\ (c,t) -> showRootIsolationTree c t) ts - -boxArea :: Representable Double ( ℝ n ) => Box n -> Double -boxArea ( 𝕀 lo hi ) = - product ( ( \ l h -> abs ( h - l ) ) <$> coordinates lo <*> coordinates hi ) - -showArea :: Double -> String -showArea area = "(area " ++ showFFloat (Just 6) area "" ++ ")" - -type Box n = 𝕀ℝ n -type BoxHistory n = [ NE.NonEmpty ( RootIsolationStep, Box n ) ] - --- | Dimension constraints for root isolation in a system of equations: --- --- - @n@: number of variables --- - @d@: number of equations --- --- NB: we require n <= d (no support for under-constrained systems). -type BoxCt n d = - ( KnownNat n, KnownNat d - , 1 <= n, 1 <= d, n <= d - - , Show ( 𝕀ℝ n ), Show ( ℝ n ) - , Eq ( ℝ n ) - , Representable Double ( ℝ n ) - , MonomialBasis ( D 1 ( ℝ n ) ) - , Deg ( D 1 ( ℝ n ) ) ~ 1 - , Vars ( D 1 ( ℝ n ) ) ~ n - , Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) - , Applicative ( Vec n ) - - , Ord ( ℝ d ) - , Module Double ( T ( ℝ d ) ) - , Representable Double ( ℝ d ) - ) - -- | Options for the root isolation methods in 'isolateRootsIn'. type RootIsolationOptions :: Nat -> Nat -> Type newtype RootIsolationOptions n d @@ -174,68 +94,11 @@ newtype RootIsolationOptions 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 ) ) + :: BoxHistory n + -> Box n + -> Either String ( NE.NonEmpty ( RootIsolationAlgorithmWithOptions n d ) ) } -type RootIsolationAlgorithm :: Nat -> Nat -> Type -data RootIsolationAlgorithm n d - -- | Bisection step. - = Bisection !( BisectionOptions n d ) - -- | Gauss–Seidel step with the given preconditioning method. - | GaussSeidel !( GaussSeidelOptions n d ) - -- | @box(1)@-consistency. - | Box1 !( Box1Options n d ) - -- | @box(2)@-consistency. - | Box2 !( Box2Options n d ) - --- | Options for the bisection method. -type BisectionOptions :: Nat -> Nat -> Type -data BisectionOptions n d = - BisectionOptions - { canHaveSols :: !( ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) -> Box n -> Bool ) - , fallbackBisectionCoord :: !( BisectionCoordPicker n d ) - } - --- | 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 ) ) - -> forall r. ( NE.NonEmpty ( Fin n, r ) -> ( r, String ) ) - --- | Options for the interval Gauss–Seidel method. -type GaussSeidelOptions :: Nat -> Nat -> Type -data GaussSeidelOptions n d = - GaussSeidelOptions - { -- | Which preconditioner to user? - gsPreconditioner :: !Preconditioner - -- | Function that projects over the equations we will consider - -- (the identity for a well-determined problem, or a projection for - -- an overdetermined system). - , 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 ] - } - --- | Options for the @box(2)@-consistency method. -data Box2Options n d = - Box2Options - { box2EpsEq :: !Double - , box2EpsBis :: !Double - , box2LambdaMin :: !Double - , box2CoordsToNarrow :: [ Fin n ] - , box2EqsToUse :: [ Fin d ] - } - defaultRootIsolationOptions :: BoxCt n d => RootIsolationOptions n d defaultRootIsolationOptions = RootIsolationOptions @@ -249,12 +112,16 @@ defaultRootIsolationOptions = defaultRootIsolationAlgorithms :: forall n d . BoxCt n d - => Double -> Double -> Box n -> BoxHistory n -> Either String ( NE.NonEmpty ( RootIsolationAlgorithm n d ) ) -defaultRootIsolationAlgorithms minWidth Ξ΅_eq box history = + => Double -- ^ minimum width of boxes (don't bisect further) + -> Double -- ^ threshold for progress + -> BoxHistory n + -> Box n + -> Either String ( NE.NonEmpty ( RootIsolationAlgorithmWithOptions n d ) ) +defaultRootIsolationAlgorithms minWidth Ξ΅_eq history box = case history of lastRoundBoxes : _ -- If, in the last round of strategies, we didn't try bisection... - | any ( \case { BisectionStep {} -> False; _ -> True } . fst ) lastRoundBoxes + | all ( \case { IsolationStep @Bisection _ -> False; _ -> True } . fst ) lastRoundBoxes , ( _step, lastRoundFirstBox ) <- NE.last lastRoundBoxes -- ...and the last round didn't sufficiently reduce the size of the box... , not $ box `sufficientlySmallerThan` lastRoundFirstBox @@ -262,51 +129,18 @@ defaultRootIsolationAlgorithms minWidth Ξ΅_eq box history = -- ...unless the box is already too small, in which case we give up. -> if verySmall then Left $ "widths <= " ++ show minWidth - else Right $ NE.singleton $ Bisection ( defaultBisectionOptions minWidth Ξ΅_eq box ) + else Right $ NE.singleton $ AlgoWithOptions @Bisection _bisOptions -- Otherwise, do a normal round. -- Currently: we try an interval Gauss–Seidel step followed by box(1)-consistency. - _ -> Right $ GaussSeidel _gaussSeidelOptions - NE.:| [ Box1 _box1Options | not verySmall ] + _ -> Right $ AlgoWithOptions @GaussSeidel _gsOptions + NE.:| [ AlgoWithOptions @Box1 _box1Options | not verySmall ] where verySmall = and $ ( \ cd -> width cd <= minWidth ) <$> coordinates box - _box1Options :: Box1Options n d - _box1Options = - Box1Options - { box1EpsEq = Ξ΅_eq - , box1EpsBis = minWidth - , box1CoordsToNarrow = toList $ universe @n - , box1EqsToUse = toList $ universe @d - } - - _box2Options :: Box2Options n d - _box2Options = - Box2Options - { box2EpsEq = Ξ΅_eq - , box2EpsBis = minWidth - , box2LambdaMin = 0.001 - , box2CoordsToNarrow = toList $ universe @n - , box2EqsToUse = toList $ universe @d - } - - _gaussSeidelOptions :: GaussSeidelOptions n d - _gaussSeidelOptions = - GaussSeidelOptions - { gsPreconditioner = InverseMidJacobian - , gsPickEqs = - case cmpNat @n @d Proxy Proxy of - EQI -> id - LTI -> - -- If there are more equations (d) than variables (n), - -- pick a size n subset of the variables, - -- (go through all combinations cyclically). - let choices :: [ Vec n ( Fin d ) ] - choices = choose @d @n - choice :: Vec n ( Fin d ) - choice = choices !! ( length history `mod` length choices ) - in \ u -> tabulate \ i -> index u ( choice ! i ) - } + _bisOptions = defaultBisectionOptions minWidth Ξ΅_eq box + _gsOptions = defaultGaussSeidelOptions history + _box1Options = defaultBox1Options minWidth Ξ΅_eq -- Did we reduce the box width by at least Ξ΅_eq -- in at least one of the coordinates? @@ -318,82 +152,13 @@ defaultRootIsolationAlgorithms minWidth Ξ΅_eq box history = <*> coordinates b2 {-# INLINEABLE defaultRootIsolationAlgorithms #-} -defaultBisectionOptions - :: forall n d - . ( 1 <= n, BoxCt n d ) - => Double -> Double - -> Box n -> BisectionOptions n d -defaultBisectionOptions minWidth _Ξ΅_eq box = - BisectionOptions - { canHaveSols = - \ eqs box' -> - -- box(0)-consistency - let iRange' :: Box d - iRange' = eqs box' `monIndex` zeroMonomial - in unT ( origin @Double ) `inside` iRange' +-------------------------------------------------------------------------------- +-- Main driver - -- box(1)-consistency - --let box1Options = Box1Options _Ξ΅_eq ( toList $ universe @n ) ( toList $ universe @d ) - --in not $ null $ makeBox1Consistent _minWidth box1Options eqs box' - - -- box(2)-consistency - --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'' - , fallbackBisectionCoord = - \ _roundHist _prevRoundsHist eqs possibleCoordChoices -> - let datPerCoord = - possibleCoordChoices <&> \ ( i, r ) -> - CoordBisectionData - { coordIndex = i - , coordInterval = box `index` i - , coordJacobianColumn = eqs box `monIndex` ( linearMonomial i ) - , coordBisectionData = r - } - - -- First, check if the largest dimension is over 20 times larger - -- 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 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 - [] -> ( coordBisectionData d1, "tooWide'" ) - Arg _ d : _ -> ( coordBisectionData d, "spread" ) - } -{-# INLINEABLE defaultBisectionOptions #-} - -sortOnArgNE :: Ord b => ( a -> b ) -> NE.NonEmpty a -> NE.NonEmpty ( Arg b a ) -sortOnArgNE f = NE.sort . fmap ( \ a -> Arg ( f a ) a ) -{-# INLINEABLE sortOnArgNE #-} - -type CoordBisectionData :: Nat -> Nat -> Type -> Type -data CoordBisectionData n d r = - CoordBisectionData - { coordIndex :: !( Fin n ) - , coordInterval :: !( 𝕀 Double ) - , coordJacobianColumn :: !( 𝕀ℝ d ) - , coordBisectionData :: !r - } -deriving stock instance ( Show ( ℝ d ), Show r ) - => Show ( CoordBisectionData n d r ) - -spread :: ( BoxCt n d, Representable Double ( ℝ d ) ) - => CoordBisectionData n d r -> Double -spread ( CoordBisectionData { coordInterval = cd, coordJacobianColumn = j_cd } ) - = width cd * normVI j_cd --maxVI j_cd -{-# INLINEABLE spread #-} - --- | Use the following algorithms using interval arithmetic in order --- to isolate roots: +-- | Use various methods using interval arithmetic in order to solve a system of +-- non-linear equations. +-- +-- Currently supported algorithms: -- -- - interval Newton method with Gauss–Seidel step for inversion -- of the interval Jacobian, @@ -401,17 +166,21 @@ spread ( CoordBisectionData { coordInterval = cd, coordJacobianColumn = j_cd } ) -- - @box(2)@-consistency algorithm, -- - bisection. -- --- Returns @(tree, (dunno, sols))@ where @tree@ is the full tree (useful for debugging), --- @sols@ are boxes that contain a unique solution (and to which Newton's method --- will converge starting from anywhere inside the box), and @dunno@ are small --- boxes which might or might not contain solutions. +-- Returns @(tree, DoneBoxes sols dunno)@ where @tree@ is the full tree +-- (useful for debugging), sols are boxes that contain a unique solution +-- (and to which Newton's method will converge starting from anywhere inside +-- the box), and @dunno@ are boxes that we have given up processing (usually +-- because they are too small); they may or may not enclose solutions. isolateRootsIn :: forall n d . BoxCt n d => RootIsolationOptions n d + -- ^ configuration (which algorithms to use, and with what parameters) -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) + -- ^ equations to solve -> Box n - -> ( [ ( Box n, RootIsolationTree ( Box n ) ) ], ( [ Box n ], [ Box n ] ) ) + -- ^ initial search domain + -> ( [ ( Box n, RootIsolationTree ( Box n ) ) ], DoneBoxes n ) isolateRootsIn ( RootIsolationOptions { rootIsolationAlgorithms } ) eqs initBox = runWriter $ map ( initBox , ) <$> go [] initBox @@ -419,7 +188,7 @@ isolateRootsIn ( RootIsolationOptions { rootIsolationAlgorithms } ) go :: BoxHistory n -> Box n -- box to work on - -> Writer ( [ Box n ], [ Box n ] ) + -> Writer ( DoneBoxes n ) [ RootIsolationTree ( Box n ) ] go history cand | -- Check the range of the equations contains zero. @@ -427,11 +196,11 @@ isolateRootsIn ( RootIsolationOptions { rootIsolationAlgorithms } ) -- Box doesn't contain a solution: discard it. = return [] | otherwise - = case rootIsolationAlgorithms cand history of + = case rootIsolationAlgorithms history cand of Right strats -> doStrategies history strats cand Left whyStop -> do -- We are giving up on this box (e.g. because it is too small). - tell ( [ cand ], [] ) + tell $ noDoneBoxes { doneGiveUpBoxes = [ ( cand, whyStop ) ] } return [ RootIsolationLeaf whyStop cand ] where iRange :: Box d @@ -440,20 +209,20 @@ isolateRootsIn ( RootIsolationOptions { rootIsolationAlgorithms } ) -- Run a round of root isolation strategies, then recur. doStrategies :: BoxHistory n - -> NonEmpty ( RootIsolationAlgorithm n d ) + -> NonEmpty ( RootIsolationAlgorithmWithOptions n d ) -> Box n - -> Writer ( [ Box n ], [ Box n ] ) + -> Writer ( DoneBoxes n ) [ RootIsolationTree ( Box n ) ] doStrategies prevRoundsHist = do_strats [] where do_strats :: [ ( RootIsolationStep, Box n ) ] - -> NE.NonEmpty ( RootIsolationAlgorithm n d ) + -> NE.NonEmpty ( RootIsolationAlgorithmWithOptions n d ) -> Box n - -> Writer ( [ Box n ], [ Box n ] ) + -> Writer ( DoneBoxes n ) [ RootIsolationTree ( Box n )] do_strats thisRoundHist ( algo NE.:| algos ) box = do -- Run one strategy in the round. - ( step, boxes ) <- doStrategy thisRoundHist prevRoundsHist eqs algo box + ( step, boxes ) <- doStrategy algo thisRoundHist prevRoundsHist eqs box case algos of -- If there are other algorithms to run in this round, run them next. nextAlgo : otherAlgos -> @@ -467,9 +236,9 @@ isolateRootsIn ( RootIsolationOptions { rootIsolationAlgorithms } ) in recur step ( go history ) boxes recur :: RootIsolationStep - -> ( Box n -> Writer ( [ Box n ], [ Box n ] ) [ RootIsolationTree ( Box n ) ] ) + -> ( Box n -> Writer ( DoneBoxes n ) [ RootIsolationTree ( Box n ) ] ) -> [ Box n ] - -> Writer ( [ Box n ], [ Box n ] ) [ RootIsolationTree ( Box n ) ] + -> Writer ( DoneBoxes n ) [ RootIsolationTree ( Box n ) ] recur step r cands = do rest <- traverse ( \ c -> do { trees <- r c; return [ (c, tree) | tree <- trees ] } ) cands return [ RootIsolationStep step $ concat rest ] @@ -479,541 +248,16 @@ isolateRootsIn ( RootIsolationOptions { rootIsolationAlgorithms } ) -- (hopefully smaller) output boxes. doStrategy :: BoxCt n d - => [ ( RootIsolationStep, Box n ) ] + => RootIsolationAlgorithmWithOptions n d + -> [ ( RootIsolationStep, Box n ) ] -> BoxHistory n -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) - -> RootIsolationAlgorithm n d -> Box n - -> Writer ( [ Box n ], [ Box n ] ) + -> Writer ( DoneBoxes n ) ( RootIsolationStep, [ Box n ] ) -doStrategy roundHistory previousRoundsHistory eqs algo box = +doStrategy algo thisRoundHist prevRoundsHist eqs cand = case algo of - GaussSeidel gsOptions -> do - boxes <- intervalGaussSeidel gsOptions eqs box - return ( GaussSeidelStep, boxes ) - Box1 box1Options -> - return ( Box1Step, makeBox1Consistent box1Options eqs box ) - Box2 box2Options -> - 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 ) + AlgoWithOptions @algo options -> do + ( step, res ) <- rootIsolationAlgorithm options thisRoundHist prevRoundsHist eqs cand + return ( SomeRootIsolationStep @algo step, res ) {-# INLINEABLE doStrategy #-} - --- | Bisect the given box. --- --- (The difficult part lies in determining along which coordinate --- dimension to bisect.) -bisect - :: forall n - . ( 1 <= n, KnownNat n, Representable Double ( ℝ n ) ) - => ( Box n -> Bool ) - -- ^ how to check whether a box contains solutions - -> ( forall r. NE.NonEmpty ( Fin n, r ) -> ( r, String ) ) - -- ^ heuristic bisection coordinate picker - -> Box n - -- ^ the box to bisect - -> ( [ Box n ], ( String, Double ) ) -bisect canHaveSols pickFallBackBisCoord box = - case findFewestSols solsList of - -- If there is a coordinate for which bisection results in no solutions, - -- or in fewer sub-boxes with solutions than any other coordinate choice, - -- pick that coordinate for bisection. - ( i, ( mid, subBoxesWithSols ) ) NE.:| [] -> - ( subBoxesWithSols, ( "cd = " ++ show i, mid ) ) - -- Otherwise, fall back to the provided heuristic. - is -> - let ( ( mid, subBoxesWithSols ), why ) = pickFallBackBisCoord is - in ( subBoxesWithSols, ( why, mid ) ) - where - solsList = - NE.fromList -- (n >= 1 so definitely non-empty) - [ Arg ( fromIntegral $ length subBoxesWithSols ) ( i, ( mid, subBoxesWithSols ) ) - | i <- toList $ universe @n - , let ( mid, ( loBox, hiBox ) ) = bisectInCoord box i - loBox_ok = canHaveSols loBox - hiBox_ok = canHaveSols hiBox - subBoxesWithSols = [ loBox | loBox_ok ] ++ [ hiBox | hiBox_ok ] - ] -{-# INLINEABLE bisect #-} - --- | Return the elements with the least argument. --- --- NB: this function shortcuts as soon as it finds an element with argument 0. -findFewestSols :: forall a. NE.NonEmpty ( Arg Word a ) -> NE.NonEmpty a -findFewestSols ( ( Arg nbSols arg ) NE.:| args ) - | nbSols == 0 - = NE.singleton arg - | otherwise - = go nbSols ( NE.singleton arg ) args - where - go :: Word -> NE.NonEmpty a -> [ Arg Word a ] -> NE.NonEmpty a - go _ bestSoFar [] = bestSoFar - go bestNbSolsSoFar bestSoFar ( ( Arg nbSols' arg' ) : args' ) - | nbSols' == 0 - = NE.singleton arg' - | otherwise - = case compare nbSols' bestNbSolsSoFar of - LT -> go nbSols' ( NE.singleton arg' ) args' - GT -> go bestNbSolsSoFar bestSoFar args' - EQ -> go bestNbSolsSoFar ( arg `NE.cons` bestSoFar ) args' - -bisectInCoord - :: Representable Double ( ℝ n ) - => Box n -> Fin n -> ( Double, ( Box n, Box n ) ) -bisectInCoord box i = - let 𝕀 z_lo z_hi = box `index` i - z_mid = 0.5 * ( z_lo + z_hi ) - in ( z_mid - , ( set i ( 𝕀 z_lo z_mid ) box - , set i ( 𝕀 z_mid z_hi ) box ) ) -{-# INLINEABLE bisectInCoord #-} - --- | Interval Newton method with Gauss–Seidel step. -intervalGaussSeidel - :: forall n d - . BoxCt n d - => GaussSeidelOptions n d - -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) - -- ^ equations - -> 𝕀ℝ n - -- ^ box - -> Writer ( [ 𝕀ℝ n ], [ 𝕀ℝ n ] ) [ 𝕀ℝ n ] -intervalGaussSeidel - ( GaussSeidelOptions { gsPreconditioner = precondMeth, gsPickEqs = pickEqs } ) - eqs - box - | let boxMid = singleton $ midpoint box - f' :: Vec n ( 𝕀ℝ n ) - 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). - ( a, b ) = precondition precondMeth - ( fmap midpoint f' ) - f' ( negateBox $ singleton $ midpoint $ f_mid ) - - -- NB: we have to change coordinates, putting the midpoint of the box - -- at the origin, in order to take a Gauss–Seidel step. - gsGuesses = map ( first ( \ box' -> unT $ box' ^+^ T boxMid ) ) - $ gaussSeidelStep a b ( T box ^-^ T boxMid ) - in - -- If the Gauss–Seidel step was a contraction, then the box - -- contains a unique solution (by the Banach fixed point theorem). - -- - -- These boxes can thus be directly added to the solution set: - -- Newton's method is guaranteed to converge to the unique solution. - let !(done, todo) = bimap ( map fst ) ( map fst ) - $ partition snd gsGuesses - in do tell ( [], done ) - return todo - where -{-# INLINEABLE intervalGaussSeidel #-} - --- | The midpoint of a box. -midpoint :: Representable Double ( ℝ n ) => 𝕀ℝ n -> ℝ n -midpoint box = - tabulate $ \ i -> - let 𝕀 z_lo z_hi = index box i - z_mid = 0.5 * ( z_lo + z_hi ) - in z_mid -{-# INLINEABLE midpoint #-} - --- | Negate a box. -negateBox :: Representable Double ( ℝ n ) => 𝕀ℝ n -> 𝕀ℝ n -negateBox v = tabulate $ \ i -> negate ( index v i ) -{-# INLINEABLE negateBox #-} - --- | Take one interval Gauss–Seidel step for the equation \( A X = B \), --- refining the initial guess box for \( X \) into up to \( 2^n \) (disjoint) new boxes. --- --- The boolean indicates whether the Gauss–Seidel step was a contraction. -gaussSeidelStep - :: forall n - . ( Representable Double ( ℝ n ), Eq ( ℝ n ) ) - => Vec n ( 𝕀ℝ n ) -- ^ columns of \( A \) - -> 𝕀ℝ n -- ^ \( B \) - -> T ( 𝕀ℝ n ) -- ^ initial box \( X \) - -> [ ( T ( 𝕀ℝ n ), Bool ) ] -gaussSeidelStep as b ( T x0 ) = coerce $ - 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 ) - ( x_i', sub_i ) <- x_i'0 `intersect` ( x `index` i ) - return $ ( set i x_i' x, sub_i && contraction ) --- TODO: try implementing the complete interval union Gauss–Seidel algorithm. --- See "Algorithm 2" in --- "Using interval unions to solve linear systems of equations with uncertainties" -{-# INLINEABLE gaussSeidelStep #-} - --- | Run an effectful computation several times in sequence, piping its output --- into the next input, once for each coordinate dimension. -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 forEachCoord #-} - --- | Compute the intersection of two intervals. --- --- Returns whether the first interval is a strict subset of the second interval, --- or the intersection is a single point. -intersect :: 𝕀 Double -> 𝕀 Double -> [ ( 𝕀 Double, Bool ) ] -intersect ( 𝕀 lo1 hi1 ) ( 𝕀 lo2 hi2 ) - | lo > hi - = [ ] - | otherwise - = [ ( 𝕀 lo hi, ( lo1 > lo2 && hi1 < hi2 ) || ( lo == hi ) ) ] - where - lo = max lo1 lo2 - hi = min hi1 hi2 - --- | Preconditioner to use with the interval Gauss–Seidel method. -data Preconditioner - = NoPreconditioning - | InverseMidJacobian - deriving stock Show - --- | 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 - :: BoxCt n d - => Box1Options n d - -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) - -> Box n -> [ Box n ] -makeBox1Consistent box1Options eqs x = - ( `State.evalState` False ) $ - 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 - => Box2Options n d - -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) - -> Box n -> Box n -makeBox2Consistent (Box2Options Ξ΅_eq Ξ΅_bis Ξ»Min coordsToNarrow eqsToUse) eqs x0 - = ( `State.evalState` False ) $ doLoop 0.25 x0 - where - box1Options :: Box1Options n d - box1Options = Box1Options Ξ΅_eq Ξ΅_bis coordsToNarrow eqsToUse - doBox1 :: Box n -> [ Box n ] - doBox1 = makeBox1Consistent box1Options eqs - doLoop :: Double -> Box n -> State Bool ( Box n ) - doLoop Ξ» x = do - x'' <- forEachCoord @n x $ boundConsistency Ξ» - modified <- State.get - let Ξ»' = if modified then Ξ» else 0.5 * Ξ» - if Ξ»' < Ξ»Min - then return x'' - else do { State.put False ; doLoop Ξ»' x'' } - - boundConsistency :: Double -> Fin n -> Box n -> State Bool ( Box n ) - boundConsistency Ξ» i box = do - let x@( 𝕀 x_inf x_sup ) = getter box - c1 = ( 1 - Ξ» ) * x_inf + Ξ» * x_sup - c2 = Ξ» * x_inf + ( 1 - Ξ» ) * x_sup - x'_inf = - case doBox1 ( setter ( 𝕀 x_inf c1 ) box ) of - [] -> c1 - x's -> minimum $ map ( inf . getter ) x's - x'_sup = - case doBox1 ( setter ( 𝕀 c2 x_sup ) box ) of - [] -> c2 - x's -> maximum $ map ( sup . getter ) x's - x' = 𝕀 x'_inf x'_sup - 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 - :: forall n - . ( KnownNat n, Representable Double ( ℝ n ) ) - => Preconditioner -- ^ pre-conditioning method to use - -> Vec n ( ℝ n ) -- ^ entry-wise midpoint matrix of the interval Jacobian matrix - -> Vec n ( 𝕀ℝ n ) -- ^ columns of \( A \) - -> 𝕀ℝ n -- ^ \( B \) - -> ( Vec n ( 𝕀ℝ n ), 𝕀ℝ n ) -precondition meth jac_mid as b = - case meth of - NoPreconditioning - -> ( as, b ) - InverseMidJacobian - | 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 ) - | otherwise - -> ( as, b ) - where - toEigen :: Vec n ( ℝ n ) -> Eigen.Matrix n n Double - toEigen cols = - Eigen.generate $ \ r c -> - ( cols ! Fin ( fromIntegral c + 1 ) ) `index` ( Fin ( fromIntegral r + 1 ) ) - - fromEigen :: Eigen.Matrix n n Double -> Vec n ( ℝ n ) - fromEigen mat = - fmap - ( \ ( Fin c ) -> - tabulate $ \ ( Fin r ) -> - Eigen.unsafeCoeff ( fromIntegral r - 1 ) ( fromIntegral c - 1 ) mat - ) - ( universe @n ) - --- | Matrix multiplication \( A v \). -matMulVec - :: forall n m - . ( Representable Double ( ℝ n ), Representable Double ( ℝ m ) ) - => Vec m ( ℝ n ) -- ^ columns of the matrix \( A ) - -> 𝕀ℝ m -- ^ vector \( v \) - -> 𝕀ℝ n -matMulVec as v = tabulate $ \ r -> - sum [ scaleInterval ( a `index` r ) ( index v c ) - | ( c, a ) <- toList ( (,) <$> universe @m <*> as ) - ] -{-# INLINEABLE matMulVec #-} - -normVI :: ( Applicative ( Vec d ), Representable Double ( ℝ d ) ) => 𝕀ℝ d -> Double -normVI ( 𝕀 los his ) = - sqrt $ sum ( nm1 <$> coordinates los <*> coordinates his ) - where - nm1 :: Double -> Double -> Double - nm1 lo hi = max ( abs lo ) ( abs hi ) Ring.^ 2 -{-# INLINEABLE normVI #-} - -maxVI :: ( Applicative ( Vec d ), Representable Double ( ℝ d ) ) => 𝕀ℝ d -> Double -maxVI ( 𝕀 los his ) = - maximum ( maxAbs <$> coordinates los <*> coordinates his ) - where - maxAbs :: Double -> Double -> Double - maxAbs lo hi = max ( abs lo ) ( abs hi ) -{-# INLINEABLE maxVI #-} - --- 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 Ξ΅_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 ) ) - 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 Ξ΄ <- 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 ) ) < Ξ΅_eq - -> x's - | otherwise - -> do - x' <- x's - if sup x' - inf x' < Ξ΅_bis - then return x' - else left_narrow =<< bisectI x' - --- TODO: de-duplicate with 'leftNarrow'? -rightNarrow :: Double - -> Double - -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) ) - -> 𝕀 Double - -> [ 𝕀 Double ] -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 ) - 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 Ξ΄ <- 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 ) ) < Ξ΅_eq - -> x's - | otherwise - -> do - x' <- x's - if sup x' - inf x' < Ξ΅_bis - then return x' - else right_narrow =<< bisectI x' - -bisectI :: 𝕀 Double -> [ 𝕀 Double ] -bisectI 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 - :: forall n d - . BoxCt n d - => Box1Options n d - -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) - -> [ Box n -> State Bool [ Box n ] ] -allNarrowingOperators ( Box1Options Ξ΅_eq Ξ΅_bis coordsToNarrow eqsToUse ) eqs = - [ \ cand -> - let getter = ( `index` coordIndex ) - setter = set coordIndex - 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 ) ) > Ξ΅_eq ) $ - -- TODO: do a check with the relative reduction in size? - State.put True - return newCands - | 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 - f, f' :: Box d - f = df `monIndex` zeroMonomial - 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' ) diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Bisection.hs b/brush-strokes/src/lib/Math/Root/Isolation/Bisection.hs new file mode 100644 index 0000000..c133bd8 --- /dev/null +++ b/brush-strokes/src/lib/Math/Root/Isolation/Bisection.hs @@ -0,0 +1,266 @@ + +{-# LANGUAGE ScopedTypeVariables #-} + +module Math.Root.Isolation.Bisection + ( -- * The bisection method for root isolation + Bisection + , bisection + + -- ** Configuration options + , BisectionCoordPicker + , BisectionOptions(..), defaultBisectionOptions + + -- *** Helper code for picking which dimension to bisect + , CoordBisectionData(..) + , spread, normVI, maxVI + , sortOnArgNE + ) + where + +-- base +import Data.Kind + ( Type ) +import Data.Foldable + ( toList ) +import Data.Functor + ( (<&>) ) +import qualified Data.List.NonEmpty as NE + ( NonEmpty(..), cons, filter + , head, nonEmpty, singleton, sort + ) +import Data.Semigroup + ( Arg(..), Dual(..) ) +import GHC.TypeNats + ( Nat, KnownNat + , type (<=) + ) + +-- MetaBrush +import Math.Algebra.Dual + ( D ) +import Math.Interval +import Math.Linear +import Math.Module + ( Module(..) ) +import Math.Monomial + ( MonomialBasis(..), linearMonomial, zeroMonomial ) +import qualified Math.Ring as Ring +import Math.Root.Isolation.Core + +-------------------------------------------------------------------------------- +-- Bisection + +-- | The bisection algorithm; see 'bisection'. +data Bisection +instance RootIsolationAlgorithm Bisection where + type instance StepDescription Bisection = ( String, Double ) + type instance RootIsolationAlgorithmOptions Bisection n d = BisectionOptions n d + rootIsolationAlgorithm + ( BisectionOptions { canHaveSols, fallbackBisectionCoord } ) + thisRoundHist prevRoundsHist eqs box = do + let ( boxes, whatBis ) = + bisection + ( canHaveSols eqs ) + ( fallbackBisectionCoord thisRoundHist prevRoundsHist eqs ) + box + return ( whatBis, boxes ) + +-- | Options for the bisection method. +type BisectionOptions :: Nat -> Nat -> Type +data BisectionOptions n d = + BisectionOptions + { -- | Custom function to check whether the given box might contain solutions to the + -- given equations. + -- + -- If you always return @True@, then we will always bisect along + -- the dimension picked by the 'fallbackBisectionCoord' function. + -- + -- NB: only return 'False' if non-existence of solutions is guaranteed + -- (otherwise, the root isolation algorithm might not be consistent). + canHaveSols :: !( ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) -> Box n -> Bool ) + -- | Heuristic to choose which coordinate dimension to bisect. + -- + -- It's only a fallback, as we prefer to bisect along coordinate dimensions + -- that minimise the number of sub-boxes created. + , fallbackBisectionCoord :: !( BisectionCoordPicker n d ) + } + +-- | A function to choose along which coordination dimension we should bisect. +type BisectionCoordPicker n d + = [ ( RootIsolationStep, Box n ) ] + -> BoxHistory n + -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) + -> forall r. ( NE.NonEmpty ( Fin n, r ) -> ( r, String ) ) + +-- | Default options for the bisection method. +defaultBisectionOptions + :: forall n d + . ( 1 <= n, BoxCt n d ) + => Double -> Double + -> Box n -> BisectionOptions n d +defaultBisectionOptions minWidth _Ξ΅_eq box = + BisectionOptions + { canHaveSols = + \ eqs box' -> + -- box(0)-consistency + let iRange' :: Box d + iRange' = eqs box' `monIndex` zeroMonomial + in unT ( origin @Double ) `inside` iRange' + + -- box(1)-consistency + --let box1Options = Box1Options _Ξ΅_eq ( toList $ universe @n ) ( toList $ universe @d ) + --in not $ null $ makeBox1Consistent _minWidth box1Options eqs box' + + -- box(2)-consistency + --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'' + , fallbackBisectionCoord = + \ _thisRoundHist _prevRoundsHist eqs possibleCoordChoices -> + let datPerCoord = + possibleCoordChoices <&> \ ( i, r ) -> + CoordBisectionData + { coordIndex = i + , coordInterval = box `index` i + , coordJacobianColumn = eqs box `monIndex` ( linearMonomial i ) + , coordBisectionData = r + } + + -- First, check if the largest dimension is over 20 times larger + -- 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 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 + [] -> ( coordBisectionData d1, "tooWide'" ) + Arg _ d : _ -> ( coordBisectionData d, "spread" ) + -- TODO: pick a dimension that previous Newton steps did not + -- manage to narrow well? + } +{-# INLINEABLE defaultBisectionOptions #-} + +sortOnArgNE :: Ord b => ( a -> b ) -> NE.NonEmpty a -> NE.NonEmpty ( Arg b a ) +sortOnArgNE f = NE.sort . fmap ( \ a -> Arg ( f a ) a ) +{-# INLINEABLE sortOnArgNE #-} + +-- | A utility datatype that is useful in implementing bisection dimension picker +-- functions ('fallbackBisectionCoord'). +type CoordBisectionData :: Nat -> Nat -> Type -> Type +data CoordBisectionData n d r = + CoordBisectionData + { coordIndex :: !( Fin n ) + , coordInterval :: !( 𝕀 Double ) + , coordJacobianColumn :: !( 𝕀ℝ d ) + , coordBisectionData :: !r + } +deriving stock instance ( Show ( ℝ d ), Show r ) + => Show ( CoordBisectionData n d r ) + +spread :: ( BoxCt n d, Representable Double ( ℝ d ) ) + => CoordBisectionData n d r -> Double +spread ( CoordBisectionData { coordInterval = cd, coordJacobianColumn = j_cd } ) + = width cd * normVI j_cd +{-# INLINEABLE spread #-} + +normVI :: ( Applicative ( Vec d ), Representable Double ( ℝ d ) ) => 𝕀ℝ d -> Double +normVI ( 𝕀 los his ) = + sqrt $ sum ( nm1 <$> coordinates los <*> coordinates his ) + where + nm1 :: Double -> Double -> Double + nm1 lo hi = max ( abs lo ) ( abs hi ) Ring.^ 2 +{-# INLINEABLE normVI #-} + +maxVI :: ( Applicative ( Vec d ), Representable Double ( ℝ d ) ) => 𝕀ℝ d -> Double +maxVI ( 𝕀 los his ) = + maximum ( maxAbs <$> coordinates los <*> coordinates his ) + where + maxAbs :: Double -> Double -> Double + maxAbs lo hi = max ( abs lo ) ( abs hi ) +{-# INLINEABLE maxVI #-} + +-------------------------------------------------------------------------------- + +-- | Bisect the given box. +-- +-- (The difficult part lies in determining along which coordinate +-- dimension to bisect.) +bisection + :: forall n + . ( 1 <= n, KnownNat n, Representable Double ( ℝ n ) ) + => ( Box n -> Bool ) + -- ^ how to check whether a box contains solutions + -> ( forall r. NE.NonEmpty ( Fin n, r ) -> ( r, String ) ) + -- ^ heuristic bisection coordinate picker + -> Box n + -- ^ the box to bisect + -> ( [ Box n ], ( String, Double ) ) +bisection canHaveSols pickFallbackBisCoord box = + case NE.nonEmpty solsList of + Nothing -> + -- We discarded dimensions along which bisection was useless + -- (because the interval was canonical in that dimension). + -- If there are no dimensions left, then don't do any bisection. + -- (TODO: we shouldn't really ever get here.) + ( [ box ], ( "noBis", 0 / 0 ) ) + Just solsNE -> + case findFewestSols solsNE of + -- If there is a coordinate for which bisection results in no solutions, + -- or in fewer sub-boxes with solutions than any other coordinate choice, + -- pick that coordinate for bisection. + Arg nbSols ( ( i, ( mid, subBoxesWithSols ) ) NE.:| [] ) -> + ( subBoxesWithSols, ( "cd = " ++ show i ++ "(#subs=" ++ show nbSols ++ ")", mid ) ) + -- Otherwise, fall back to the provided heuristic. + Arg _nbSols is -> + let ( ( mid, subBoxesWithSols ), why ) = pickFallbackBisCoord is + in ( subBoxesWithSols, ( why, mid ) ) + where + solsList = + [ Arg ( fromIntegral $ length subBoxesWithSols ) ( i, ( mid, subBoxesWithSols ) ) + | i <- toList $ universe @n + , let ( mid, subBoxes ) = bisectInCoord box i + subBoxesWithSols = NE.filter canHaveSols subBoxes + -- discard coordinate dimensions in which the box is a singleton + , length subBoxes >= 2 || null subBoxesWithSols + ] +{-# INLINEABLE bisection #-} + +-- | Bisect a box in the given coordinate dimension. +bisectInCoord + :: Representable Double ( ℝ n ) + => Box n -> Fin n -> ( Double, NE.NonEmpty ( Box n ) ) +bisectInCoord box i = + let z = box `index` i + zs' = bisect z + in ( sup ( NE.head zs' ) + , fmap ( \ z' -> set i z' box ) zs' ) +{-# INLINEABLE bisectInCoord #-} + +-- | Return the elements with the least argument. +-- +-- NB: this function shortcuts as soon as it finds an element with argument 0. +findFewestSols :: forall a. NE.NonEmpty ( Arg Word a ) -> Arg Word ( NE.NonEmpty a ) +findFewestSols ( Arg nbSols arg NE.:| args ) + | nbSols == 0 + = Arg 0 $ NE.singleton arg + | otherwise + = go nbSols ( NE.singleton arg ) args + where + go :: Word -> NE.NonEmpty a -> [ Arg Word a ] -> Arg Word ( NE.NonEmpty a ) + go bestNbSolsSoFar bestSoFar [] = Arg bestNbSolsSoFar bestSoFar + go bestNbSolsSoFar bestSoFar ( ( Arg nbSols' arg' ) : args' ) + | nbSols' == 0 + = Arg 0 $ NE.singleton arg' + | otherwise + = case compare nbSols' bestNbSolsSoFar of + LT -> go nbSols' ( NE.singleton arg' ) args' + GT -> go bestNbSolsSoFar bestSoFar args' + EQ -> go bestNbSolsSoFar ( arg `NE.cons` bestSoFar ) args' diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Core.hs b/brush-strokes/src/lib/Math/Root/Isolation/Core.hs new file mode 100644 index 0000000..498a455 --- /dev/null +++ b/brush-strokes/src/lib/Math/Root/Isolation/Core.hs @@ -0,0 +1,265 @@ + +{-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE UndecidableInstances #-} + +-- | Core definitions and utilities common to root isolation methods. +module Math.Root.Isolation.Core + ( -- * Root isolation types + Box + , DoneBoxes(..), noDoneBoxes + + -- * General typeclass for root isolation methods + , BoxCt + , RootIsolationAlgorithm(..) + , RootIsolationAlgorithmWithOptions(..) + + -- ** Inspecting history + , RootIsolationStep(IsolationStep, ..) + , BoxHistory + + -- ** Visualising history + , RootIsolationTree(..) + , showRootIsolationTree + + -- * Utility functions + , pipeFunctionsWhileTrue + , forEachCoord + ) where + +-- base +import Data.Foldable + ( toList ) +import Data.Kind + ( Type, Constraint ) +import qualified Data.List.NonEmpty as NE + ( NonEmpty ) +import Data.Type.Equality + ( (:~~:)(HRefl) ) +import Data.Typeable + ( Typeable, heqT ) +import GHC.TypeNats + ( Nat, KnownNat, type (<=) ) +import Numeric + ( showFFloat ) + +-- containers +import Data.Tree + ( Tree(..) ) + +-- transformers +import Control.Monad.Trans.State.Strict as State + ( State, get, put ) +import Control.Monad.Trans.Writer.CPS + ( Writer ) + +-- MetaBrush +import Math.Algebra.Dual + ( D ) +import Math.Interval +import Math.Linear +import Math.Module + ( Module(..) ) +import Math.Monomial + ( MonomialBasis(..), Deg, Vars ) + +-------------------------------------------------------------------------------- + +-- | An axis-aligned box in @n@-dimensions. +type Box n = 𝕀ℝ n + +-- | Dimension constraints for root isolation in a system of equations: +-- +-- - @n@: number of variables +-- - @d@: number of equations +-- +-- NB: we require n <= d (no support for under-constrained systems). +-- +-- NB: in practice, this constraint should specialise away. +type BoxCt n d = + ( KnownNat n, KnownNat d + , 1 <= n, 1 <= d, n <= d + + , Show ( 𝕀ℝ n ), Show ( ℝ n ) + , Eq ( ℝ n ) + , Representable Double ( ℝ n ) + , MonomialBasis ( D 1 ( ℝ n ) ) + , Deg ( D 1 ( ℝ n ) ) ~ 1 + , Vars ( D 1 ( ℝ n ) ) ~ n + , Module Double ( T ( ℝ n ) ) + , Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) + , Applicative ( Vec n ) + + , Ord ( ℝ d ) + , Module Double ( T ( ℝ d ) ) + , Representable Double ( ℝ d ) + ) + +-- | Boxes we are done with and will not continue processing. +data DoneBoxes n = + DoneBoxes + { -- | Boxes which definitely contain a unique solution. + doneSolBoxes :: ![ Box n ] + -- | Boxes which may or may not contain solutions, + -- and that we have stopped processing for some reason. + , doneGiveUpBoxes :: ![ ( Box n, String ) ] + } +deriving stock instance Show ( Box n ) => Show ( DoneBoxes n ) + +instance Semigroup ( DoneBoxes n ) where + DoneBoxes a1 b1 <> DoneBoxes a2 b2 = DoneBoxes ( a1 <> a2 ) ( b1 <> b2 ) +instance Monoid ( DoneBoxes n ) where + mempty = noDoneBoxes +noDoneBoxes :: DoneBoxes n +noDoneBoxes = DoneBoxes [] [] + +-------------------------------------------------------------------------------- +-- Class for all root isolation algorithms. +-- +-- This keeps the implementation open-ended, and allows inspection of +-- other root isolation methods, so that heuristics can look at +-- what happened in previous steps to decide what to do. + +-- | Existential wrapper over any root isolation algorithm, +-- with the options necessary to run it. +data RootIsolationAlgorithmWithOptions n d where + AlgoWithOptions + :: RootIsolationAlgorithm ty + => RootIsolationAlgorithmOptions ty n d + -> RootIsolationAlgorithmWithOptions n d + +-- | Type-class for root isolation algorithms. +-- +-- This design keeps the set of root isolation algorithms open-ended, +-- while retaining the ability to inspect previous steps (using the +-- 'IsolationStep' pattern). +type RootIsolationAlgorithm :: forall {k}. k -> Constraint +class ( Typeable ty, Show ( StepDescription ty ) ) + => RootIsolationAlgorithm ty where + -- | The type of additional information about an algorithm step. + -- + -- Only really useful for debugging; gets stored in 'RootIsolationTree's. + type StepDescription ty + -- | Configuration options expected by this root isolation method. + type RootIsolationAlgorithmOptions ty (n :: Nat) (d :: Nat) = r | r -> ty n d + -- | Run one step of the root isolation method. + -- + -- This gets given the equations and a box, and should attempt to + -- shrink the box in some way, returning smaller boxes. + -- + -- Should return: + -- + -- - a description of the step taken (see 'StepDescription'), + -- - new boxes to process (the return value of type @['Box' n]@), + -- which can be empty if the algorithm can prove that the input + -- bix does not contain any solutions; + -- - (as a writer side-effect) boxes to definitely stop processing; see 'DoneBoxes'. + rootIsolationAlgorithm + :: forall (n :: Nat) (d :: Nat) + . BoxCt n d + => RootIsolationAlgorithmOptions ty n d + -- ^ options for this root isolation algorithm + -> [ ( RootIsolationStep, Box n ) ] + -- ^ history of the current round + -> BoxHistory n + -- ^ previous rounds history + -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) + -- ^ equations + -> Box n + -- ^ box + -> Writer ( DoneBoxes n ) ( StepDescription ty, [ Box n ] ) + +-- | Match on an unknown root isolation algorithm step with a known algorithm. +pattern IsolationStep + :: forall (ty :: Type) + . RootIsolationAlgorithm ty + => StepDescription ty -> RootIsolationStep +pattern IsolationStep stepDescr + <- ( rootIsolationAlgorithmStep_maybe @ty -> Just stepDescr ) + where + IsolationStep stepDescr = SomeRootIsolationStep @ty stepDescr + +-- | Helper function used to define the 'IsolationStep' pattern. +-- +-- Inspects whether an existential 'RootIsolationStep' packs a step for +-- the given algorithm. +rootIsolationAlgorithmStep_maybe + :: forall ty. RootIsolationAlgorithm ty + => RootIsolationStep -> Maybe ( StepDescription ty ) +rootIsolationAlgorithmStep_maybe ( SomeRootIsolationStep @existential descr ) + | Just HRefl <- heqT @existential @ty + = Just descr + | otherwise + = Nothing +{-# INLINEABLE rootIsolationAlgorithmStep_maybe #-} + +-- | History for a given box: what was the outcome of previous root isolation +-- methods? +type BoxHistory n = [ NE.NonEmpty ( RootIsolationStep, Box n ) ] + +-- | A description of a step taken when isolating roots. +data RootIsolationStep where + SomeRootIsolationStep + :: forall step + . ( Typeable step + , Show ( StepDescription step ) + ) + => StepDescription step + -> RootIsolationStep + +instance Show RootIsolationStep where + showsPrec p ( SomeRootIsolationStep stepDescr ) = showsPrec p stepDescr + +-------------------------------------------------------------------------------- +-- Trees recording steps taken by the algorithm, for visualisation & debugging. + +-- | A tree recording the steps taken when isolating roots. +data RootIsolationTree d + = RootIsolationLeaf String d + | RootIsolationStep RootIsolationStep [ ( d, RootIsolationTree d ) ] + +showRootIsolationTree + :: ( Representable Double ( ℝ n ), Show ( Box n ) ) + => Box n -> RootIsolationTree ( Box n ) -> Tree String +showRootIsolationTree cand (RootIsolationLeaf why l) = Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ why ++ " " ++ show l) [] +showRootIsolationTree cand (RootIsolationStep s ts) + = Node (show cand ++ " abc " ++ showArea (boxArea cand) ++ " " ++ show s) $ map (\ (c,t) -> showRootIsolationTree c t) ts + +boxArea :: Representable Double ( ℝ n ) => Box n -> Double +boxArea ( 𝕀 lo hi ) = + product ( ( \ l h -> abs ( h - l ) ) <$> coordinates lo <*> coordinates hi ) + +showArea :: Double -> String +showArea area = "(area " ++ showFFloat (Just 6) area "" ++ ")" + +-------------------------------------------------------------------------------- +-- Utilities. + +-- | Run an effectful computation several times in sequence, piping its output +-- into the next input, once for each coordinate dimension. +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 forEachCoord #-} + +-- | 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 diff --git a/brush-strokes/src/lib/Math/Root/Isolation/GaussSeidel.hs b/brush-strokes/src/lib/Math/Root/Isolation/GaussSeidel.hs new file mode 100644 index 0000000..75bfeca --- /dev/null +++ b/brush-strokes/src/lib/Math/Root/Isolation/GaussSeidel.hs @@ -0,0 +1,239 @@ +{-# LANGUAGE ScopedTypeVariables #-} + +module Math.Root.Isolation.GaussSeidel + ( -- * The interval Newton method with Gauss–Seidel step + GaussSeidel + , intervalGaussSeidel + + -- ** Configuration options + , GaussSeidelOptions(..), Preconditioner(..) + , defaultGaussSeidelOptions + ) + where + +-- base +import Control.Arrow + ( first ) +import Data.Bifunctor + ( Bifunctor(bimap) ) +import Data.Coerce + ( coerce ) +import Data.Kind + ( Type ) +import Data.Foldable + ( toList ) +import Data.List + ( partition ) +import Data.Proxy + ( Proxy(..) ) +import Data.Type.Ord + ( OrderingI(..) ) +import GHC.TypeNats + ( Nat, KnownNat, type (<=), cmpNat ) + +-- eigen +import qualified Eigen.Matrix as Eigen + ( Matrix + , determinant, generate, inverse, unsafeCoeff + ) + +-- transformers +import Control.Monad.Trans.Writer.CPS + ( Writer, tell ) + +-- MetaBrush +import Math.Algebra.Dual + ( D ) +import Math.Epsilon + ( nearZero ) +import Math.Interval +import Math.Linear +import Math.Module + ( Module(..) ) +import Math.Monomial + ( MonomialBasis(..), linearMonomial, zeroMonomial ) +import Math.Root.Isolation.Core + +-------------------------------------------------------------------------------- +-- Gauss–Seidel + +-- | The interval Newton method with a Gauss–Seidel step; see 'intervalGaussSeidel'. +data GaussSeidel +instance RootIsolationAlgorithm GaussSeidel where + type instance StepDescription GaussSeidel = () + type instance RootIsolationAlgorithmOptions GaussSeidel n d = GaussSeidelOptions n d + rootIsolationAlgorithm opts _thisRoundHist _prevRoundsHist eqs box = do + res <- intervalGaussSeidel opts eqs box + return ( (), res ) + {-# INLINEABLE rootIsolationAlgorithm #-} + +-- | Options for the interval Gauss–Seidel method. +type GaussSeidelOptions :: Nat -> Nat -> Type +data GaussSeidelOptions n d = + GaussSeidelOptions + { -- | Which preconditioner to user? + gsPreconditioner :: !Preconditioner + -- | Function that projects over the equations we will consider + -- (the identity for a well-determined problem, or a projection for + -- an overdetermined system). + , gsPickEqs :: ( 𝕀ℝ d -> 𝕀ℝ n ) } + +-- | Default options for the interval Gauss–Seidel method. +defaultGaussSeidelOptions + :: forall n d + . ( KnownNat n, KnownNat d + , 1 <= n, 1 <= d, n <= d + , Representable Double ( ℝ n ) + , Representable Double ( ℝ d ) + ) + => BoxHistory n + -> GaussSeidelOptions n d +defaultGaussSeidelOptions history = + GaussSeidelOptions + { gsPreconditioner = InverseMidJacobian + , gsPickEqs = + case cmpNat @n @d Proxy Proxy of + EQI -> id + LTI -> + -- If there are more equations (d) than variables (n), + -- pick a size n subset of the variables + -- (go through all combinations cyclically). + let choices :: [ Vec n ( Fin d ) ] + choices = choose @d @n + choice :: Vec n ( Fin d ) + choice = choices !! ( length history `mod` length choices ) + in \ u -> tabulate \ i -> index u ( choice ! i ) + } +{-# INLINEABLE defaultGaussSeidelOptions #-} + +-- | Preconditioner to use with the interval Gauss–Seidel method. +data Preconditioner + = NoPreconditioning + | InverseMidJacobian + deriving stock Show + +-- | Interval Newton method with Gauss–Seidel step. +intervalGaussSeidel + :: forall n d + . BoxCt n d + => GaussSeidelOptions n d + -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) + -- ^ equations + -> 𝕀ℝ n + -- ^ box + -> Writer ( DoneBoxes n ) [ 𝕀ℝ n ] +intervalGaussSeidel + ( GaussSeidelOptions { gsPreconditioner = precondMeth, gsPickEqs = pickEqs } ) + eqs + box + | let boxMid = singleton $ boxMidpoint box + f' :: Vec n ( 𝕀ℝ n ) + 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). + ( a, b ) = precondition precondMeth + ( fmap boxMidpoint f' ) + f' ( singleton $ unT $ -1 *^ T ( boxMidpoint f_mid ) ) + + -- NB: we have to change coordinates, putting the midpoint of the box + -- at the origin, in order to take a Gauss–Seidel step. + gsGuesses = map ( first ( \ box' -> unT $ box' ^+^ T boxMid ) ) + $ gaussSeidelStep a b ( T box ^-^ T boxMid ) + in + -- If the Gauss–Seidel step was a contraction, then the box + -- contains a unique solution (by the Banach fixed point theorem). + -- + -- These boxes can thus be directly added to the solution set: + -- Newton's method is guaranteed to converge to the unique solution. + let !(done, todo) = bimap ( map fst ) ( map fst ) + $ partition snd gsGuesses + in do tell $ noDoneBoxes { doneSolBoxes = done } + return todo + where +{-# INLINEABLE intervalGaussSeidel #-} + +-- | Take one interval Gauss–Seidel step for the equation \( A X = B \), +-- refining the initial guess box for \( X \) into up to \( 2^n \) (disjoint) new boxes. +-- +-- The boolean indicates whether the Gauss–Seidel step was a contraction. +gaussSeidelStep + :: forall n + . ( Representable Double ( ℝ n ), Eq ( ℝ n ) ) + => Vec n ( 𝕀ℝ n ) -- ^ columns of \( A \) + -> 𝕀ℝ n -- ^ \( B \) + -> T ( 𝕀ℝ n ) -- ^ initial box \( X \) + -> [ ( T ( 𝕀ℝ n ), Bool ) ] +gaussSeidelStep as b ( T x0 ) = coerce $ + 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 ) + ( x_i', sub_i ) <- x_i'0 `intersect` ( x `index` i ) + return $ ( set i x_i' x, sub_i && contraction ) +-- TODO: try implementing the complete interval union Gauss–Seidel algorithm. +-- See "Algorithm 2" in +-- "Using interval unions to solve linear systems of equations with uncertainties" +{-# INLINEABLE gaussSeidelStep #-} + +-- | The midpoint of a box. +boxMidpoint :: Representable Double ( ℝ n ) => 𝕀ℝ n -> ℝ n +boxMidpoint box = + tabulate $ \ i -> + let 𝕀 z_lo z_hi = index box i + z_mid = 0.5 * ( z_lo + z_hi ) + in z_mid +{-# INLINEABLE boxMidpoint #-} + +-- | Pre-condition the system \( AX = B \). +precondition + :: forall n + . ( KnownNat n, Representable Double ( ℝ n ) ) + => Preconditioner -- ^ pre-conditioning method to use + -> Vec n ( ℝ n ) -- ^ entry-wise midpoint matrix of the interval Jacobian matrix + -> Vec n ( 𝕀ℝ n ) -- ^ columns of \( A \) + -> 𝕀ℝ n -- ^ \( B \) + -> ( Vec n ( 𝕀ℝ n ), 𝕀ℝ n ) +precondition meth jac_mid as b = + case meth of + NoPreconditioning + -> ( as, b ) + InverseMidJacobian + | 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 ) + | otherwise + -> ( as, b ) + where + toEigen :: Vec n ( ℝ n ) -> Eigen.Matrix n n Double + toEigen cols = + Eigen.generate $ \ r c -> + ( cols ! Fin ( fromIntegral c + 1 ) ) `index` ( Fin ( fromIntegral r + 1 ) ) + + fromEigen :: Eigen.Matrix n n Double -> Vec n ( ℝ n ) + fromEigen mat = + fmap + ( \ ( Fin c ) -> + tabulate $ \ ( Fin r ) -> + Eigen.unsafeCoeff ( fromIntegral r - 1 ) ( fromIntegral c - 1 ) mat + ) + ( universe @n ) +{-# INLINEABLE precondition #-} + +-- | Matrix multiplication \( A v \). +matMulVec + :: forall n m + . ( Representable Double ( ℝ n ), Representable Double ( ℝ m ) ) + => Vec m ( ℝ n ) -- ^ columns of the matrix \( A ) + -> 𝕀ℝ m -- ^ vector \( v \) + -> 𝕀ℝ n +matMulVec as v = tabulate $ \ r -> + sum [ scaleInterval ( a `index` r ) ( index v c ) + | ( c, a ) <- toList ( (,) <$> universe @m <*> as ) + ] +{-# INLINEABLE matMulVec #-} diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Narrowing.hs b/brush-strokes/src/lib/Math/Root/Isolation/Narrowing.hs new file mode 100644 index 0000000..d292718 --- /dev/null +++ b/brush-strokes/src/lib/Math/Root/Isolation/Narrowing.hs @@ -0,0 +1,490 @@ +{-# LANGUAGE ScopedTypeVariables #-} + +module Math.Root.Isolation.Narrowing + ( -- * @box(1)@-consistency + Box1 + , makeBox1Consistent + -- ** Configuration options + , Box1Options(..) + , defaultBox1Options + + , -- *** Narrowing methods for @box(1)@-consistency + NarrowingMethod(..) + , narrowingMethods + + -- **** Options for the adaptive shaving method + , AdaptiveShavingOptions(..) + , defaultAdaptiveShavingOptions + + -- * @box(2)@-consistency + , Box2 + , makeBox2Consistent + + -- ** Configuration options + , Box2Options(..) + , defaultBox2Options + ) + where + +-- base +import Control.Monad + ( when ) +import Data.Foldable + ( toList ) +import qualified Data.List.NonEmpty as NE + ( toList ) +import GHC.TypeNats + ( KnownNat ) + +-- transformers +import Control.Monad.Trans.State.Strict as State + ( State, evalState, get, put ) + +-- brush-strokes +import Math.Algebra.Dual + ( D ) +import Math.Float.Utils + ( succFP, prevFP ) +import Math.Interval +import Math.Linear + ( ℝ + , Representable + , Fin, set, index, universe + ) +import Math.Monomial + ( MonomialBasis(..), Deg, Vars + , zeroMonomial, linearMonomial + ) +import Math.Root.Isolation.Core + +-------------------------------------------------------------------------------- +-- Box-consistency driver code + +-- | A @box(1)@-consistency enforcing algorithm; see 'makeBox1Consistent'. +data Box1 +instance RootIsolationAlgorithm Box1 where + type instance StepDescription Box1 = () + type instance RootIsolationAlgorithmOptions Box1 n d = Box1Options n d + rootIsolationAlgorithm opts _thisRoundHist _prevRoundsHist eqs box = + return $ ( (), makeBox1Consistent opts eqs box ) + {-# INLINEABLE rootIsolationAlgorithm #-} + +-- | A @box(2)@-consistency enforcing algorithm; see 'makeBox1Consistent'. +data Box2 +instance RootIsolationAlgorithm Box2 where + type instance StepDescription Box2 = () + type instance RootIsolationAlgorithmOptions Box2 n d = Box2Options n d + rootIsolationAlgorithm opts _thisRoundHist _prevRoundsHist eqs box = + return ( () , [ makeBox2Consistent opts eqs box ] ) + {-# INLINEABLE rootIsolationAlgorithm #-} + +-- | Options for the @box(1)@-consistency method. +data Box1Options n d = + Box1Options + { box1EpsEq :: !Double + , box1EpsBis :: !Double + , box1CoordsToNarrow :: !( [ Fin n ] ) + , box1EqsToUse :: !( [ Fin d ] ) + , box1NarrowingMethod :: !NarrowingMethod + } + deriving stock Show + +-- | Options for the @box(2)@-consistency method. +data Box2Options n d = + Box2Options + { box2Box1Options :: !( Box1Options n d ) + , box2EpsEq :: !Double + , box2LambdaMin :: !Double + } + deriving stock Show + +-- | Default options for the @box(1)@-consistency method. +defaultBox1Options + :: forall n d + . ( KnownNat n, KnownNat d ) + => Double -- ^ minimum width of boxes (don't bisect further) + -> Double -- ^ threshold for progress + -> Box1Options n d +defaultBox1Options minWidth Ξ΅_eq = + Box1Options + { box1EpsEq = Ξ΅_eq + , box1EpsBis = minWidth + , box1CoordsToNarrow = toList $ universe @n + , box1EqsToUse = toList $ universe @d + , box1NarrowingMethod = Kubica + } +{-# INLINEABLE defaultBox1Options #-} + +-- | Default options for the @box(2)@-consistency method. +defaultBox2Options + :: forall n d + . ( KnownNat n, KnownNat d ) + => Double -- ^ minimum width of boxes (don't bisect further) + -> Double -- ^ threshold for progress + -> Box2Options n d +defaultBox2Options minWidth Ξ΅_eq = + Box2Options + { box2Box1Options = defaultBox1Options minWidth Ξ΅_eq + , box2EpsEq = Ξ΅_eq + , box2LambdaMin = 0.001 + } + +-- | 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 + :: ( KnownNat n, Representable Double ( ℝ n ) + , MonomialBasis ( D 1 ( ℝ n ) ) + , Deg ( D 1 ( ℝ n ) ) ~ 1 + , Vars ( D 1 ( ℝ n ) ) ~ n + , Representable Double ( ℝ d ) + ) + => Box1Options n d + -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) + -> 𝕀ℝ n -> [ 𝕀ℝ n ] +makeBox1Consistent box1Options eqs x = + ( `State.evalState` False ) $ + pipeFunctionsWhileTrue ( allNarrowingOperators box1Options eqs ) x +{-# INLINEABLE defaultBox2Options #-} + +-- | 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 + . ( KnownNat n, Representable Double ( ℝ n ) + , MonomialBasis ( D 1 ( ℝ n ) ) + , Deg ( D 1 ( ℝ n ) ) ~ 1 + , Vars ( D 1 ( ℝ n ) ) ~ n + , Representable Double ( ℝ d ) + ) + => Box2Options n d + -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) + -> 𝕀ℝ n -> 𝕀ℝ n +makeBox2Consistent (Box2Options box1Options Ξ΅_eq Ξ»Min) eqs x0 + = ( `State.evalState` False ) $ doLoop 0.25 x0 + where + doBox1 :: 𝕀ℝ n -> [ 𝕀ℝ n ] + doBox1 = makeBox1Consistent box1Options eqs + doLoop :: Double -> 𝕀ℝ n -> State Bool ( 𝕀ℝ n ) + doLoop Ξ» x = do + x'' <- forEachCoord @n x $ boundConsistency Ξ» + modified <- State.get + let Ξ»' = if modified then Ξ» else 0.5 * Ξ» + if Ξ»' < Ξ»Min + then return x'' + else do { State.put False ; doLoop Ξ»' x'' } + + boundConsistency :: Double -> Fin n -> 𝕀ℝ n -> State Bool ( 𝕀ℝ n ) + boundConsistency Ξ» i box = do + let x@( 𝕀 x_inf x_sup ) = getter box + c1 = ( 1 - Ξ» ) * x_inf + Ξ» * x_sup + c2 = Ξ» * x_inf + ( 1 - Ξ» ) * x_sup + x'_inf = + case doBox1 ( setter ( 𝕀 x_inf c1 ) box ) of + [] -> c1 + x's -> minimum $ map ( inf . getter ) x's + x'_sup = + case doBox1 ( setter ( 𝕀 c2 x_sup ) box ) of + [] -> c2 + x's -> maximum $ map ( sup . getter ) x's + x' = 𝕀 x'_inf x'_sup + when ( width x - width x' >= Ξ΅_eq ) $ + State.put True + return $ setter x' box + where + getter = ( `index` i ) + setter = set i + +-------------------------------------------------------------------------------- +-- Narrowing methods + +-- | The narrowing method to use to enforce @box(1)@-consistency. +data NarrowingMethod + -- | Algorithm 5 from the paper + -- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems" + -- (BartΕ‚omiej Jacek Kubica, 2017) + = Kubica + -- | Two-sided shaving @sbc@ from the paper + -- "A Data-Parallel Algorithm to Reliably Solve Systems of Nonlinear Equations", + -- (Goldsztejn & Goualard, 2008) + | TwoSidedShaving + -- | @sbc3ag@ (adaptive guessing) shaving, from the paper + -- "Box Consistency through Adaptive Shaving" + -- (Goldsztejn & Goualard, 2010). + | AdaptiveShaving AdaptiveShavingOptions + deriving stock Show + +narrowingMethods + :: Double -> Double + -> NarrowingMethod + -> [ ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) ) -> 𝕀 Double -> [ 𝕀 Double ] ] +narrowingMethods Ξ΅_eq Ξ΅_bis Kubica + = [ leftNarrow Ξ΅_eq Ξ΅_bis, rightNarrow Ξ΅_eq Ξ΅_bis ] +narrowingMethods Ξ΅_eq Ξ΅_bis ( AdaptiveShaving opts ) + = [ leftShave Ξ΅_eq Ξ΅_bis opts, rightNarrow Ξ΅_eq Ξ΅_bis ] + -- TODO: haven't implemented right shaving yet +narrowingMethods _Ξ΅_eq Ξ΅_bis TwoSidedShaving + = [ sbc Ξ΅_bis ] +{-# INLINE narrowingMethods #-} + +allNarrowingOperators + :: forall n d + . ( KnownNat n, Representable Double ( ℝ n ) + , MonomialBasis ( D 1 ( ℝ n ) ) + , Deg ( D 1 ( ℝ n ) ) ~ 1 + , Vars ( D 1 ( ℝ n ) ) ~ n + , Representable Double ( ℝ d ) + ) + => Box1Options n d + -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) + -> [ 𝕀ℝ n -> State Bool [ 𝕀ℝ n ] ] +allNarrowingOperators ( Box1Options Ξ΅_eq Ξ΅_bis coordsToNarrow eqsToUse narrowingMethod ) eqs = + [ \ cand -> + let getter = ( `index` coordIndex ) + setter = set coordIndex + newCands = map ( `setter` cand ) + $ narrowFn ( \ box -> ff' coordIndex eqnIndex $ setter box cand ) + ( getter cand ) + in do + -- Record when we achieved a meaningful reduction, + -- so that we continue trying further narrowings. + when ( ( width ( getter cand ) - sum ( map ( width . getter ) newCands ) ) >= Ξ΅_eq ) $ + -- NB: making this return 'True' less often seems slightly beneficial? + -- Further investigation needed. + State.put True + return newCands + | narrowFn <- narrowingMethods Ξ΅_eq Ξ΅_bis narrowingMethod + , coordIndex <- coordsToNarrow + , eqnIndex <- eqsToUse + ] + where + ff' :: Fin n -> Fin d -> 𝕀ℝ n -> ( 𝕀 Double, 𝕀 Double ) + ff' i d ts = + let df = eqs ts + f, f' :: 𝕀ℝ d + f = df `monIndex` zeroMonomial + f' = df `monIndex` linearMonomial i + in ( f `index` d, f' `index` d ) +{-# INLINEABLE allNarrowingOperators #-} + +-------------------------------------------------------------------------------- +-- Kubica's algorithm. + +-- 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 Ξ΅_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 ) ) + 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 Ξ΄ <- 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 ) ) < Ξ΅_eq + -> x's + | otherwise + -> do + x' <- x's + if sup x' - inf x' < Ξ΅_bis + then return x' + else left_narrow =<< NE.toList ( bisect x' ) + +-- TODO: de-duplicate with 'leftNarrow'? +rightNarrow :: Double + -> Double + -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) ) + -> 𝕀 Double + -> [ 𝕀 Double ] +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 ) + 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 Ξ΄ <- 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 ) ) < Ξ΅_eq + -> x's + | otherwise + -> do + x' <- x's + if sup x' - inf x' < Ξ΅_bis + then return x' + else right_narrow =<< NE.toList ( bisect x' ) + +-------------------------------------------------------------------------------- +-- Adaptive shaving. + +data AdaptiveShavingOptions + = AdaptiveShavingOptions + { Ξ³_init, Οƒ_good, Οƒ_bad, Ξ²_good, Ξ²_bad :: !Double + } + deriving stock Show + +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 + -> Double + -> AdaptiveShavingOptions + -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) ) + -> 𝕀 Double + -> [ 𝕀 Double ] +leftShave Ξ΅_eq Ξ΅_bis + ( 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 < Ξ΅_bis + = [ 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. + = Ξ³ + xs' = x_minus_guess : guesses' + in if ( width x' - sum ( map width xs' ) ) < Ξ΅_eq + then xs' + else left_narrow Ξ³' =<< xs' +{-# INLINEABLE leftShave #-} + +-------------------------------------------------------------------------------- +-- Two-sided shaving. + +-- | @sbc@ algorithm from the paper +-- +-- "A Data-Parallel Algorithm to Reliably Solve Systems of Nonlinear Equations", +-- (FrΓ©dΓ©ric Goualard, Alexandre Goldsztejn, 2008). +sbc :: Double + -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) ) + -> 𝕀 Double -> [ 𝕀 Double ] +sbc Ξ΅ ff' = go + where + go :: 𝕀 Double -> [ 𝕀 Double ] + go x@( 𝕀 x_l x_r ) + | width x <= Ξ΅ + = [ x ] + | otherwise + = let x_mid = 0.5 * ( x_l + x_r ) + ( left_done, left_todo ) + | 0 `inside` ( fst $ ff' ( 𝕀 x_l x_lp ) ) + = ( True, [ ] ) + | not $ 0 `inside` ( fst $ ff' i_l ) + = ( False, [ ] ) + | otherwise + = let + xls = do + let l = 𝕀 x_lp x_lp + Ξ΄ <- fst ( ff' l ) `extendedDivide` snd ( ff' i_l ) + map fst $ ( l - Ξ΄ ) `intersect` i_l + in ( False, xls ) + where x_lp = min ( succFP x_l ) x_mid + i_l = 𝕀 x_lp x_mid + ( right_done, right_todo ) + | 0 `inside` ( fst $ ff' ( 𝕀 x_rm x_r ) ) + = ( True, [ ] ) + | not $ 0 `inside` ( fst $ ff' i_r ) + = ( False, [ ] ) + | otherwise + = let + xrs = do + let r = 𝕀 x_rm x_rm + Ξ΄ <- fst ( ff' r ) `extendedDivide` snd ( ff' i_r ) + map fst $ ( r - Ξ΄ ) `intersect` i_r + in ( False, xrs ) + where x_rm = max ( prevFP x_r ) x_mid + i_r = 𝕀 x_mid x_rm + in do let lefts' = if left_done + then [ 𝕀 x_l x_mid ] + else go =<< left_todo + rights' = if right_done + then [ 𝕀 x_mid x_r ] + else go =<< right_todo + lefts' ++ rights' +{-# INLINEABLE sbc #-}