diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs index 67b8853..9e4366d 100644 --- a/brush-strokes/src/cusps/bench/Main.hs +++ b/brush-strokes/src/cusps/bench/Main.hs @@ -182,7 +182,7 @@ benchCases :: [ TestCase ] benchCases = [ ellipseTestCase opts ("minWidth=" ++ show minWidth ++ ",Ξ΅=" ++ show narrowAbs) ( 0, 1 ) pi $ defaultStartBoxes [ 2 ] | minWidth <- [ 1e-5 ] - , narrowAbs <- [ 1e-5, 1e-4, 1e-3, 5e-3, 8e-3, 1e-2, 2e-2, 3e-2, 4e-2, 5e-2, 1e-1 ] + , narrowAbs <- [ 5e-2 ] , let opts = RootIsolationOptions { minWidth diff --git a/brush-strokes/src/lib/Math/Linear.hs b/brush-strokes/src/lib/Math/Linear.hs index 7225b35..3d500f3 100644 --- a/brush-strokes/src/lib/Math/Linear.hs +++ b/brush-strokes/src/lib/Math/Linear.hs @@ -130,6 +130,7 @@ newtype Vec n a = Vec { vecList :: [ a ] } universe :: forall n. KnownNat n => Vec n ( Fin n ) universe = Vec [ Fin i | i <- [ 1 .. fromIntegral ( natVal' @n proxy# ) ] ] +{-# INLINEABLE universe #-} coordinates :: forall r u. ( Representable r u ) => u -> Vec ( RepDim u ) r coordinates u = fmap ( index u ) $ universe @( RepDim u ) @@ -160,7 +161,7 @@ zipIndices ( Vec v ) = zipIndices_ 1 v -------------------------------------------------------------------------------- -- | Rotate a vector by the given angle (counter-clockwise), --- given the cosine and sine of the angle (in that order) +-- given the cosine and sine of the angle (in that order). rotate :: ( Representable r m, RepDim m ~ 2, Ring r ) => r -- \( \cos \theta \) -> r -- \( \sin \theta \) diff --git a/brush-strokes/src/lib/Math/Root/Isolation.hs b/brush-strokes/src/lib/Math/Root/Isolation.hs index 5bd3d3a..89248c2 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation.hs @@ -44,12 +44,14 @@ import Data.Kind ( Type ) import Data.Foldable ( toList ) +import Data.Functor + ( (<&>) ) import Data.List - ( partition, sort ) + ( partition ) import Data.List.NonEmpty ( NonEmpty ) import qualified Data.List.NonEmpty as NE - ( NonEmpty(..), last, singleton ) + ( NonEmpty(..), cons, filter, fromList, last, singleton, sort ) import Data.Semigroup ( Arg(..), Dual(..) ) import Numeric @@ -182,9 +184,17 @@ type BisectionOptions :: Nat -> Nat -> Type data BisectionOptions n d = BisectionOptions { canHaveSols :: !( ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) -> Box n -> Bool ) - , fallbackBisectionDim :: !( [ ( RootIsolationStep, Box n ) ] -> BoxHistory n -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) -> ( Fin n, String ) ) + , fallbackBisectionDim :: !( BisectionDimPicker n d ) } +type BisectionDimPicker n d + = forall r + . [ ( RootIsolationStep, Box n ) ] + -> BoxHistory n + -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) + -> NE.NonEmpty ( Fin n, r ) + -> ( r, String ) + -- | Options for the interval Gauss–Seidel method. type GaussSeidelOptions :: Nat -> Nat -> Type data GaussSeidelOptions n d = @@ -192,7 +202,7 @@ data GaussSeidelOptions n d = { -- | Which preconditioner to user? gsPreconditioner :: !Preconditioner -- | Function that projects over the equations we will consider - -- (the identity for a well-determined problem, or a project for + -- (the identity for a well-determined problem, or a projection for -- an overdetermined system). , gsDims :: ( 𝕀ℝ d -> 𝕀ℝ n ) } @@ -246,7 +256,7 @@ defaultRootIsolationAlgorithms minWidth narrowAbs box history -> NE.singleton $ Bisection ( defaultBisectionOptions minWidth narrowAbs box ) -- Otherwise, do a normal round. -- Currently: we try an interval Gauss–Seidel step followed by box(1)-consistency. - _ -> GaussSeidel defaultGaussSeidelOptions + _ -> GaussSeidel _gaussSeidelOptions NE.:| [ Box1 _box1Options ] where @@ -254,7 +264,7 @@ defaultRootIsolationAlgorithms minWidth narrowAbs box history _box1Options = Box1Options { box1EpsEq = narrowAbs - , box1CoordsToNarrow = toList $ universe @n -- [ Fin 1, Fin 2 ] + , box1CoordsToNarrow = toList $ universe @n , box1EqsToUse = toList $ universe @d } @@ -267,23 +277,27 @@ defaultRootIsolationAlgorithms minWidth narrowAbs box history , box2EqsToUse = toList $ universe @d } + _gaussSeidelOptions :: GaussSeidelOptions n d + _gaussSeidelOptions = + GaussSeidelOptions + { gsPreconditioner = InverseMidJacobian + , gsDims = + \ ( 𝕀 ( ℝ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 ) + 1 -> 𝕀 ( ℝ2 b_lo c_lo ) ( ℝ2 b_hi c_hi ) + _ -> 𝕀 ( ℝ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? sufficientlySmallerThan :: Box n -> Box n -> Bool b1 `sufficientlySmallerThan` b2 = or $ - ( \ cd1 cd2 -> width cd2 - width cd1 > narrowAbs ) + ( \ cd1 cd2 -> width cd2 - width cd1 >= narrowAbs ) <$> coordinates b1 <*> coordinates b2 {-# INLINEABLE defaultRootIsolationAlgorithms #-} -defaultGaussSeidelOptions :: GaussSeidelOptions N 3 -defaultGaussSeidelOptions = - GaussSeidelOptions - { gsPreconditioner = InverseMidJacobian - , gsDims = \ ( 𝕀 ( ℝ3 _a_lo b_lo c_lo ) ( ℝ3 _a_hi b_hi c_hi ) ) - -> 𝕀 ( ℝ2 b_lo c_lo ) ( ℝ2 b_hi c_hi ) -} - defaultBisectionOptions :: forall n d . ( 1 <= n, BoxCt n d ) @@ -295,7 +309,7 @@ defaultBisectionOptions minWidth _narrowAbs box = \ eqs box' -> -- box(0)-consistency let iRange' :: Box d - iRange' = ( `monIndex` zeroMonomial ) $ eqs box' + iRange' = eqs box' `monIndex` zeroMonomial in unT ( origin @Double ) `inside` iRange' -- box(1)-consistency @@ -309,52 +323,52 @@ defaultBisectionOptions minWidth _narrowAbs box = -- iRange'' = eqs box'' `monIndex` zeroMonomial --in unT ( origin @Double ) `inside` iRange'' , fallbackBisectionDim = - \ _roundHist _prevRoundsHist eqs -> - let df = eqs box - datPerCoord = - [ CoordBisectionData - { coordIndex = i - , coordInterval = box `index` i - , coordJacobianColumn = df `monIndex` ( linearMonomial i ) - } - | i <- toList $ universe @n ] + \ _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 10 times larger + -- First, check if the largest dimension is over 20 times larger -- than the smallest dimension; if so bisect along that coordinate. - in case sortOnArg ( width . coordInterval ) datPerCoord of - [] -> error "impossible: dimension 0" - [Arg _ d] -> (coordIndex d, "") - Arg w0 _ : ds -> + 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 >= 10 * w0 - then ( coordIndex d1, "tooWide" ) + in if w1 >= 20 * w0 + then ( coordBisectionData d1, "tooWide" ) -- Otherwise, pick the dimension with the largest spread = width * Jacobian column norm - else case filter ( not . isTooSmall ) $ sortOnArg ( Dual . spread ) datPerCoord of - [] -> ( coordIndex d1, "tooWide'" ) - Arg _ d : _ -> ( coordIndex d, "spread" ) + 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" ) } - where - isTooSmall ( Arg ( Dual w ) _ ) = w < minWidth {-# INLINEABLE defaultBisectionOptions #-} -sortOnArg :: Ord b => (a -> b) -> [a] -> [Arg b a] -sortOnArg f = sort . map ( \ a -> Arg ( f a ) a ) -{-# INLINEABLE sortOnArg #-} +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 -data CoordBisectionData n d = +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 ( CoordBisectionData n d ) +deriving stock instance ( Show ( ℝ d ), Show r ) + => Show ( CoordBisectionData n d r ) spread :: ( BoxCt n d, Representable Double ( ℝ d ) ) - => CoordBisectionData n d -> Double + => CoordBisectionData n d r -> Double spread ( CoordBisectionData { coordInterval = cd, coordJacobianColumn = j_cd } ) - = width cd * normVI j_cd + = width cd * normVI j_cd --maxVI j_cd {-# INLINEABLE spread #-} -- | Use the following algorithms using interval arithmetic in order @@ -400,13 +414,12 @@ isolateRootsIn = case cuspFindingAlgorithms 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, - -- or we have reached an iteration depth). + -- We are giving up on this box (e.g. because it is too small). tell ( [ cand ], [] ) return [ RootIsolationLeaf whyStop cand ] where iRange :: Box d - iRange = ( `monIndex` zeroMonomial ) $ eqs cand + iRange = eqs cand `monIndex` zeroMonomial -- Run a round of cusp finding strategies, then recur. doStrategies @@ -477,56 +490,57 @@ doStrategy roundHistory previousRoundsHistory eqs minWidth algo box = -- (The difficult part lies in determining along which dimension to bisect.) bisect :: forall n - . ( KnownNat n, Representable Double ( ℝ n ) ) + . ( 1 <= n, KnownNat n, Representable Double ( ℝ n ) ) => ( Box n -> Bool ) -- ^ how to check whether a box contains solutions - -> ( Fin n, String ) - -- ^ fallback choice of dimension (and "why" we chose it) + -> ( 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 ( fallbackDim, why ) box = - -- We try to bisect along a dimension which eliminates zeros from one of the - -- sub-regions. - -- - -- If this fails, we fall back to the provided dimension choice. +bisect canHaveSols pickFallBackBisCoord box = case findFewestSols solsList of - Just ( Arg _nbSols ( i, mid, oks ) ) -> - ( oks, ("cd = " ++ show i, mid ) ) - Nothing -> - case bisectInCoord box fallbackDim of - ( mid, ( lo, hi ) ) -> - ( [ lo, hi ], ( why, mid ) ) + -- 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 = - [ Arg ( fromIntegral $ length oks ) ( i, mid, oks ) - | i <- toList $ universe @n - , let (mid, (lo, hi)) = bisectInCoord box i - lo_ok = canHaveSols lo - hi_ok = canHaveSols hi - oks = [ lo | lo_ok ] ++ [ hi | hi_ok ] - ] + 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 #-} --- | Find any element with the least argument. +-- | Return the elements with the least argument. -- -- NB: this function shortcuts as soon as it finds an element with argument 0. -findFewestSols :: forall a. [ Arg Word a ] -> Maybe ( Arg Word a ) -findFewestSols [] = Nothing -findFewestSols ( arg@( Arg nbSols _ ) : args ) +findFewestSols :: forall a. NE.NonEmpty ( Arg Word a ) -> NE.NonEmpty a +findFewestSols ( ( Arg nbSols arg ) NE.:| args ) | nbSols == 0 - = Just arg + = NE.singleton arg | otherwise - = Just $ go arg args + = go nbSols ( NE.singleton arg ) args where - go :: Arg Word a -> [ Arg Word a ] -> Arg Word a - go bestSoFar [] = bestSoFar - go bestSoFar@( Arg bestNbSolsSoFar _ ) ( arg'@( Arg nbSols' _ ) : args' ) + go :: Word -> NE.NonEmpty a -> [ Arg Word a ] -> NE.NonEmpty a + go _ bestSoFar [] = bestSoFar + go bestNbSolsSoFar bestSoFar ( ( Arg nbSols' arg' ) : args' ) | nbSols' == 0 - = arg' - | nbSols' < bestNbSolsSoFar - = go arg' args' + = NE.singleton arg' | otherwise - = go bestSoFar args' + = 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 ) @@ -755,8 +769,8 @@ matMulVec -> 𝕀ℝ m -- ^ vector \( v \) -> 𝕀ℝ n matMulVec as v = tabulate $ \ r -> - sum [ scaleInterval ( ( as ! c ) `index` r ) ( index v c ) - | c <- toList ( universe @m ) + sum [ scaleInterval ( a `index` r ) ( index v c ) + | ( c, a ) <- toList ( (,) <$> universe @m <*> as ) ] {-# INLINEABLE matMulVec #-} @@ -768,6 +782,14 @@ normVI ( 𝕀 los his ) = 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. --