diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs index 76fb50e..17eb870 100644 --- a/brush-strokes/src/cusps/bench/Main.hs +++ b/brush-strokes/src/cusps/bench/Main.hs @@ -69,7 +69,7 @@ benchTestCase testName ( TestCase { testDescription, testBrushStroke, testCuspOp ( dunno, sols ) = foldMap ( \ ( i, ( _trees, DoneBoxes { doneSolBoxes = defCusps, doneGiveUpBoxes = mbCusps } ) ) -> - ( map ( ( i , ) . snd ) mbCusps, map ( i, ) defCusps ) ) $ + ( map ( ( i , ) . fst ) mbCusps, map ( i, ) defCusps ) ) $ IntMap.toList $ findCuspsIn testCuspOptions testStrokeFnI $ IntMap.fromList diff --git a/brush-strokes/src/lib/Math/Linear.hs b/brush-strokes/src/lib/Math/Linear.hs index b246be4..7948193 100644 --- a/brush-strokes/src/lib/Math/Linear.hs +++ b/brush-strokes/src/lib/Math/Linear.hs @@ -128,8 +128,17 @@ pattern V4 x y z w = T ( ℝ4 x y z w ) type Vec :: Nat -> Type -> Type newtype Vec n a = Vec { vecList :: [ a ] } - deriving newtype ( Show, Eq, Ord, Functor, Foldable ) - deriving Applicative via ZipList +type role Vec nominal representational + +deriving newtype instance Show a => Show ( Vec n a ) +deriving newtype instance Eq a => Eq ( Vec n a ) +deriving newtype instance Ord a => Ord ( Vec n a ) +deriving newtype instance Functor ( Vec n ) +deriving newtype instance Foldable ( Vec n ) +deriving via ZipList + instance Applicative ( Vec n ) +instance Traversable ( Vec n ) where + traverse f ( Vec as ) = Vec <$> traverse f as universe :: forall n. KnownNat n => Vec n ( Fin n ) universe = Vec [ Fin i | i <- [ 1 .. fromIntegral ( natVal' @n proxy# ) ] ] diff --git a/brush-strokes/src/lib/Math/Root/Isolation.hs b/brush-strokes/src/lib/Math/Root/Isolation.hs index 579fbc5..b21d15b 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation.hs @@ -138,9 +138,9 @@ defaultRootIsolationAlgorithms minWidth ε_eq history box = where verySmall = and $ ( \ cd -> width cd <= minWidth ) <$> coordinates box - _bisOptions = defaultBisectionOptions minWidth ε_eq box - _gsOptions = defaultGaussSeidelOptions history - _box1Options = defaultBox1Options minWidth ε_eq + _bisOptions = defaultBisectionOptions @n @d minWidth ε_eq box + _gsOptions = defaultGaussSeidelOptions @n @d history + _box1Options = defaultBox1Options @n @d minWidth ε_eq -- Did we reduce the box width by at least ε_eq -- in at least one of the coordinates? diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Core.hs b/brush-strokes/src/lib/Math/Root/Isolation/Core.hs index 2f28a2d..f3e2294 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation/Core.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation/Core.hs @@ -89,7 +89,6 @@ type BoxCt n d = , Vars ( D 1 ( ℝ n ) ) ~ n , Module Double ( T ( ℝ n ) ) , Module ( 𝕀 Double ) ( T ( 𝕀ℝ n ) ) - , Applicative ( Vec n ) , Ord ( ℝ d ) , Module Double ( T ( ℝ d ) ) diff --git a/brush-strokes/src/lib/Math/Root/Isolation/GaussSeidel.hs b/brush-strokes/src/lib/Math/Root/Isolation/GaussSeidel.hs index 75bfeca..aca5411 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation/GaussSeidel.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation/GaussSeidel.hs @@ -12,6 +12,7 @@ module Math.Root.Isolation.GaussSeidel where -- base +import Prelude hiding ( unzip ) import Control.Arrow ( first ) import Data.Bifunctor @@ -24,6 +25,8 @@ import Data.Foldable ( toList ) import Data.List ( partition ) +import Data.List.NonEmpty + ( unzip ) import Data.Proxy ( Proxy(..) ) import Data.Type.Ord @@ -71,12 +74,21 @@ instance RootIsolationAlgorithm GaussSeidel where type GaussSeidelOptions :: Nat -> Nat -> Type data GaussSeidelOptions n d = GaussSeidelOptions - { -- | Which preconditioner to user? + { -- | Which preconditioner to use? 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 ) } + , gsPickEqs :: !( 𝕀ℝ d -> 𝕀ℝ n ) + -- | Whether to use a partial or a complete Gauss–Seidel update + , gsUpdate :: !GaussSeidelUpdateMethod + } + +-- | Whether to use a partial or a complete Gauss–Seidel update. +data GaussSeidelUpdateMethod + = GS_Partial + | GS_Complete + deriving stock Show -- | Default options for the interval Gauss–Seidel method. defaultGaussSeidelOptions @@ -90,7 +102,7 @@ defaultGaussSeidelOptions -> GaussSeidelOptions n d defaultGaussSeidelOptions history = GaussSeidelOptions - { gsPreconditioner = InverseMidJacobian + { gsPreconditioner = InverseMidpoint , gsPickEqs = case cmpNat @n @d Proxy Proxy of EQI -> id @@ -103,13 +115,14 @@ defaultGaussSeidelOptions history = choice :: Vec n ( Fin d ) choice = choices !! ( length history `mod` length choices ) in \ u -> tabulate \ i -> index u ( choice ! i ) + , gsUpdate = GS_Complete } {-# INLINEABLE defaultGaussSeidelOptions #-} -- | Preconditioner to use with the interval Gauss–Seidel method. data Preconditioner = NoPreconditioning - | InverseMidJacobian + | InverseMidpoint deriving stock Show -- | Interval Newton method with Gauss–Seidel step. @@ -123,44 +136,66 @@ intervalGaussSeidel -- ^ box -> Writer ( DoneBoxes n ) [ 𝕀ℝ n ] intervalGaussSeidel - ( GaussSeidelOptions { gsPreconditioner = precondMeth, gsPickEqs = pickEqs } ) + ( GaussSeidelOptions + { gsPreconditioner = precondMeth + , gsPickEqs = pickEqs + , gsUpdate + } ) 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 + x + | let x_mid = singleton $ boxMidpoint x + f'_x :: Vec n ( 𝕀ℝ n ) + f'_x = fmap ( \ i -> pickEqs $ eqs x `monIndex` linearMonomial i ) ( universe @n ) + f_x_mid = pickEqs $ eqs x_mid `monIndex` zeroMonomial = let -- Interval Newton method: take one Gauss–Seidel step -- for the system of equations f'(x) ( x - x_mid ) = - f(x_mid). + minus_f_x_mid = unT $ -1 *^ T ( boxMidpoint f_x_mid ) + + -- Precondition the above linear system into A ( x - x_mid ) = B. ( a, b ) = precondition precondMeth - ( fmap boxMidpoint f' ) - f' ( singleton $ unT $ -1 *^ T ( boxMidpoint f_mid ) ) + f'_x ( singleton minus_f_x_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 ) + gsGuesses = map ( first ( \ x' -> unT $ x' ^+^ T x_mid ) ) + $ gaussSeidelUpdate gsUpdate a b ( T x ^-^ T x_mid ) 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 + let ( done, todo ) = bimap ( map fst ) ( map fst ) + $ partition snd gsGuesses in do tell $ noDoneBoxes { doneSolBoxes = done } return todo where {-# INLINEABLE intervalGaussSeidel #-} +-- | A partial or complete Gauss–Seidel step for the equation \( A X = B \), +-- refining the initial guess box for \( X \) into up to \( 2^n \) (disjoint) new boxes. +gaussSeidelUpdate + :: forall n + . ( Representable Double ( ℝ n ), n ~ RepDim ( ℝ n ), Eq ( ℝ n ) ) + => GaussSeidelUpdateMethod -- ^ which step method to use + -> Vec n ( 𝕀ℝ n ) -- ^ columns of \( A \) + -> 𝕀ℝ n -- ^ \( B \) + -> T ( 𝕀ℝ n ) -- ^ initial box \( X \) + -> [ ( T ( 𝕀ℝ n ), Bool ) ] +gaussSeidelUpdate upd as b x = + case upd of + GS_Partial -> gaussSeidelStep as b x + GS_Complete -> gaussSeidelStep_Complete as b x +{-# INLINEABLE gaussSeidelUpdate #-} + -- | 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 ) ) + . ( Representable Double ( ℝ n ), n ~ RepDim ( ℝ n ), Eq ( ℝ n ) ) => Vec n ( 𝕀ℝ n ) -- ^ columns of \( A \) -> 𝕀ℝ n -- ^ \( B \) -> T ( 𝕀ℝ n ) -- ^ initial box \( X \) @@ -168,15 +203,80 @@ gaussSeidelStep 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" + let s = b `index` i - sum [ ( as ! j ) `index` i * x `index` j + | j <- toList ( universe @n ), j /= i ] + x_i = x `index` i + a_ii = ( as ! i ) `index` i + -- Take a shortcut before performing the division if possible. + if | not $ 0 `inside` ( s - a_ii * x_i ) + -- No solutions: don't bother performing a division. + -> [ ] + | 0 `inside` s && 0 `inside` a_ii + -- The division would produce [-oo,+oo]: don't do anything. + -> [ ( x, False ) ] + -- Otherwise, perform the division. + | otherwise + -> do + x_i'0 <- s `extendedDivide` a_ii + ( x_i', sub_i ) <- x_i'0 `intersect` x_i + return $ ( set i x_i' x, sub_i && contraction ) {-# INLINEABLE gaussSeidelStep #-} +-- | The complete interval-union Gauss–Seidel step. +-- +-- Algorithm 2 from: +-- "Using interval unions to solve linear systems of equations with uncertainties" +-- (Montanher, Domes, Schichl, Neumaier) (2017) +gaussSeidelStep_Complete + :: forall n + . ( Representable Double ( ℝ n ), n ~ RepDim ( ℝ n ) ) + => Vec n ( 𝕀ℝ n ) -- ^ columns of \( A \) + -> 𝕀ℝ n -- ^ \( B \) + -> T ( 𝕀ℝ n ) -- ^ initial box \( X \) + -> [ ( T ( 𝕀ℝ n ), Bool ) ] +gaussSeidelStep_Complete as b ( T x0 ) = coerce $ do + ( x', subs ) <- + forEachCoord @n ( x0, pure False ) $ \ i ( x, contractions ) -> do + let s = b `index` i - sum [ ( as ! k ) `index` i * x `index` k + | k <- toList ( universe @n ) ] + x_i = x `index` i + a_ii = ( as ! i ) `index` i + ( x', subs ) <- fromComponents \ j -> do + let x_j = x `index` j + a_ij = ( as ! j ) `index` i + s_j = s `ominus` ( a_ij * x_j ) + -- Shortcut division if possible (see gaussSeidelStep for commentary). + if | not $ 0 `inside` ( s_j - a_ii * x_i ) + -> [ ] + | 0 `inside` s_j && 0 `inside` a_ij + -> [ ( x_j, False ) ] + | otherwise + -> do + x_j'0 <- s_j `extendedDivide` a_ij + ( x_j', sub_j ) <- x_j'0 `intersect` x_j + return $ ( x_j', sub_j ) + return ( x', (||) <$> contractions <*> subs ) + return ( x', and subs ) +{-# INLINEABLE gaussSeidelStep_Complete #-} + +fromComponents + :: forall n + . ( Representable Double ( ℝ n ), n ~ RepDim ( ℝ n ) ) + => ( Fin n -> [ ( 𝕀 Double, Bool ) ] ) -> [ ( 𝕀ℝ n, Vec n Bool ) ] +fromComponents f = do + ( xs, bs ) <- unzip <$> traverse f ( universe @n ) + return $ ( tabulate $ \ i -> xs ! i, bs ) + -- TODO: this could be more efficient. +{-# INLINEABLE fromComponents #-} + +infixl 6 `ominus` +ominus :: 𝕀 Double -> 𝕀 Double -> 𝕀 Double +ominus a@( 𝕀 lo1 hi1 ) b@( 𝕀 lo2 hi2 ) + | width a >= width b + = 𝕀 ( lo1 - lo2 ) ( hi1 - hi2 ) + | otherwise + = 𝕀 ( hi1 - hi2 ) ( lo1 - lo2 ) + -- | The midpoint of a box. boxMidpoint :: Representable Double ( ℝ n ) => 𝕀ℝ n -> ℝ n boxMidpoint box = @@ -191,16 +291,15 @@ 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 = +precondition meth as b = case meth of NoPreconditioning -> ( as, b ) - InverseMidJacobian - | let mat = toEigen jac_mid + InverseMidpoint + | let mat = toEigen $ fmap boxMidpoint as det = Eigen.determinant mat , not $ nearZero det -- (TODO: a bit wasteful to compute determinant then inverse.)