From 131753da8269edd4f3e65c61dfe382e58478294a Mon Sep 17 00:00:00 2001 From: sheaf Date: Wed, 17 Apr 2024 20:41:21 +0200 Subject: [PATCH] Split off root-isolation algorithms into Math.Root.Isolation --- brush-strokes/brush-strokes.cabal | 1 + brush-strokes/src/cusps/bench/Main.hs | 43 +- brush-strokes/src/lib/Math/Bezier/Stroke.hs | 970 +++--------------- .../Math/Bezier/Stroke/EnvelopeEquation.hs | 8 +- brush-strokes/src/lib/Math/Interval.hs | 4 +- brush-strokes/src/lib/Math/Linear.hs | 2 +- brush-strokes/src/lib/Math/Linear/Internal.hs | 36 +- brush-strokes/src/lib/Math/Root/Isolation.hs | 746 ++++++++++++++ 8 files changed, 938 insertions(+), 872 deletions(-) create mode 100644 brush-strokes/src/lib/Math/Root/Isolation.hs diff --git a/brush-strokes/brush-strokes.cabal b/brush-strokes/brush-strokes.cabal index 3fe1b13..0bf66a7 100644 --- a/brush-strokes/brush-strokes.cabal +++ b/brush-strokes/brush-strokes.cabal @@ -124,6 +124,7 @@ library , Math.Orientation , Math.Ring , Math.Roots + , Math.Root.Isolation , Debug.Utils other-modules: diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs index 5459d36..3fa5fcd 100644 --- a/brush-strokes/src/cusps/bench/Main.hs +++ b/brush-strokes/src/cusps/bench/Main.hs @@ -40,6 +40,8 @@ import Numeric ( showFFloat ) -- containers +import qualified Data.IntMap.Strict as IntMap + ( fromList, toList ) import Data.Sequence ( Seq ) import qualified Data.Sequence as Seq @@ -72,6 +74,7 @@ import Math.Linear import Math.Module import Math.Ring ( Transcendental ) +import Math.Root.Isolation -------------------------------------------------------------------------------- @@ -120,7 +123,12 @@ benchTestCase ( TestCase { testName, testBrushStroke, testCuspOptions, testStart let ( _, testStrokeFnI ) = brushStrokeFunctions testBrushStroke ( dunno, sols ) = foldMap ( \ ( i, ( _trees, ( mbCusps, defCusps ) ) ) -> ( map ( i , ) mbCusps, map ( i, ) defCusps ) ) $ - computeCusps testCuspOptions testStrokeFnI testStartBoxes + IntMap.toList $ + findCuspsIn testCuspOptions testStrokeFnI $ + IntMap.fromList + [ ( i, [ box ] ) + | ( i, box ) <- testStartBoxes + ] rnf dunno `seq` rnf sols `seq` return () after <- getMonotonicTime let dt = after - before @@ -175,9 +183,7 @@ benchCases :: [ TestCase ] benchCases = [ ellipseTestCase opts "full" ( 0, 1 ) pi $ defaultStartBoxes [ 0 .. 3 ] ] where - opts = defaultCuspFindingOptions - - + opts = defaultRootIsolationOptions -------------------------------------------------------------------------------- @@ -193,7 +199,7 @@ data TestCase = TestCase { testName :: String , testBrushStroke :: BrushStroke - , testCuspOptions :: CuspFindingOptions + , testCuspOptions :: RootIsolationOptions , testStartBoxes :: [ ( Int, Box ) ] } @@ -250,7 +256,7 @@ take 10 $ Data.List.sortOn ( \ ( _, ℝ1 e, v) -> abs e + norm v ) [ let { v = m potentialCusp $ eval fI $ mkBox (0.5798, 0.5799) 3 (0.26798, 0.26799) > True -let nbPotentialSols b = let ( _newtTrees, ( dunno, sols ) ) = findCuspsFrom NoPreconditioning 1e-7 fI b in length dunno + length sols +let nbPotentialSols b = let ( _newtTrees, ( dunno, sols ) ) = isolateRootsIn NoPreconditioning 1e-7 fI b in length dunno + length sols nbPotentialSols $ mkBox (0.5798, 0.5799) 3 (0.26798, 0.26799) 1 @@ -258,7 +264,7 @@ nbPotentialSols $ mkBox (0.5798, 0.5799) 3 (0.26798, 0.26799) nbPotentialSols $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799) 0 -let showTrees b = map ( uncurry showIntervalNewtonTree ) $ fst $ findCuspsFrom NoPreconditioning 1e-7 fI b +let showTrees b = map ( uncurry showIntervalNewtonTree ) $ fst $ isolateRootsIn NoPreconditioning 1e-7 fI b putStrLn $ unlines $ map Data.Tree.View.showTree $ showTrees $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799) @@ -464,7 +470,7 @@ maximum [ _y $ sup $ unT $ _D12_v $ dEdsdcdt $ eval fI $ mkBox (t, t) 2 (s, s) | -let showTrees b = map ( uncurry showIntervalNewtonTree ) $ fst $ findCuspsFrom NoPreconditioning 1e-7 fI b +let showTrees b = map ( uncurry showIntervalNewtonTree ) $ fst $ isolateRootsIn NoPreconditioning 1e-7 fI b putStrLn $ unlines $ map Data.Tree.View.showTree $ showTrees $ mkBox (0.5486101933245248, 0.5486102071493622) 2 (0.548095036738487, 0.5480952) @@ -482,7 +488,7 @@ showD float = showFFloat (Just 6) float "" -------------------------------------------------------------------------------- -ellipseTestCase :: CuspFindingOptions -> String -> ( Double, Double ) -> Double -> [ ( Int, Box ) ] -> TestCase +ellipseTestCase :: RootIsolationOptions -> String -> ( Double, Double ) -> Double -> [ ( Int, Box ) ] -> TestCase ellipseTestCase opts str k0k1 rot startBoxes = TestCase { testName = "ellipse (" ++ str ++ ")" @@ -512,7 +518,7 @@ trickyCusp2TestCase = TestCase { testName = "trickyCusp2" , testBrushStroke = trickyCusp2BrushStroke - , testCuspOptions = defaultCuspFindingOptions + , testCuspOptions = defaultRootIsolationOptions , testStartBoxes = defaultStartBoxes [ 0 .. 3 ] } @@ -601,19 +607,6 @@ getStrokeFunctions brush sp0 crv = brush @3 @𝕀 proxy# singleton ) {-# INLINEABLE getStrokeFunctions #-} -computeCusps - :: CuspFindingOptions - -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) - -> [ ( Int, Box ) ] - -> [ ( Int, ( [ ( Box, RootFindingTree Box ) ], ( [ Box ], [ Box ] ) ) ) ] -computeCusps cuspOpts eqs startBoxes = - map ( \ ( i, box ) -> ( i , ) $ findCuspsFrom cuspOpts ( mkEqn i ) box ) startBoxes - where - mkEqn i ( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) = - let t = 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) - s = 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) - in ( eqs t `Seq.index` i ) s - defaultStartBoxes :: [ Int ] -> [ ( Int, Box ) ] defaultStartBoxes is = [ ( i, mkBox (zero, one) (zero, one) ) | i <- is ] @@ -623,8 +616,8 @@ getR1 (ℝ1 u) = u {- (f, fI) = brushStrokeFunctions $ ellipseBrushStroke (0,1) pi -nbPotentialSols box = let ( _newtTrees, ( dunno, sols ) ) = findCuspsFrom NoPreconditioning 1e-7 fI box in length dunno + length sols -showTrees box = putStrLn $ unlines $ map Data.Tree.View.showTree $ map ( uncurry showIntervalNewtonTree ) $ fst $ findCuspsFrom NoPreconditioning 1e-7 fI box +nbPotentialSols box = let ( _newtTrees, ( dunno, sols ) ) = isolateRootsIn NoPreconditioning 1e-7 fI box in length dunno + length sols +showTrees box = putStrLn $ unlines $ map Data.Tree.View.showTree $ map ( uncurry showIntervalNewtonTree ) $ fst $ isolateRootsIn NoPreconditioning 1e-7 fI box sol_t = 0.5486100729150693796677845183880669025324233347060776339185 :: Double sol_s = 0.5480950141859386853197594577293968665598143630958601978245 :: Double diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke.hs b/brush-strokes/src/lib/Math/Bezier/Stroke.hs index d483e24..93b3140 100644 --- a/brush-strokes/src/lib/Math/Bezier/Stroke.hs +++ b/brush-strokes/src/lib/Math/Bezier/Stroke.hs @@ -1,8 +1,6 @@ {-# LANGUAGE AllowAmbiguousTypes #-} {-# LANGUAGE PolyKinds #-} -{-# LANGUAGE QuantifiedConstraints #-} {-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE UndecidableInstances #-} module Math.Bezier.Stroke ( Offset(..), Cusp(..) @@ -19,14 +17,7 @@ module Math.Bezier.Stroke , brushStrokeData, pathAndUsedParams -- ** Cusp finding - , RootFindingTree(..), showRootFindingTree - , RootFindingStep(..) - , RootFindingLeaf(..) - , Box - , CuspFindingOptions(..), defaultCuspFindingOptions - , RootFindingAlgorithm(..), defaultCuspFindingAlgorithms - , Preconditioner(..) - , findCusps, findCuspsFrom, gaussSeidel + , findCusps, findCuspsIn ) where @@ -34,7 +25,7 @@ module Math.Bezier.Stroke import Control.Arrow ( first, (***) ) import Control.Monad - ( unless, when ) + ( unless ) import Control.Monad.ST ( RealWorld, ST ) import Data.Bifunctor @@ -48,11 +39,11 @@ import Data.Foldable import Data.Functor.Identity ( Identity(..) ) import Data.List - ( nub, partition, sort ) + ( sort ) import Data.List.NonEmpty ( NonEmpty ) import qualified Data.List.NonEmpty as NE - ( NonEmpty(..), cons, last, singleton, unzip ) + ( cons, singleton, unzip ) import Data.Maybe ( fromMaybe, isJust, listToMaybe, mapMaybe ) import Data.Semigroup @@ -67,8 +58,6 @@ import GHC.Generics ( Generic, Generic1, Generically(..) ) import GHC.TypeNats ( type (-) ) -import Numeric - ( showFFloat ) -- acts import Data.Act @@ -82,13 +71,11 @@ import Data.Act import Data.IntMap.Strict ( IntMap ) import qualified Data.IntMap.Strict as IntMap - ( fromList, toList ) + ( fromList, mapWithKey, toList ) import Data.Sequence ( Seq(..) ) import qualified Data.Sequence as Seq ( empty, index, length, reverse, singleton, zipWith ) -import Data.Tree - ( Tree(..) ) -- deepseq import Control.DeepSeq @@ -97,7 +84,7 @@ import Control.DeepSeq -- generic-lens import Data.Generics.Product.Typed ( HasType(typed) ) -import Data.Generics.Internal.VL +import qualified Data.Generics.Internal.VL as Lens ( set, view ) -- parallel @@ -110,9 +97,9 @@ import Control.Monad.Trans.Class import Control.Monad.Trans.Except ( Except, runExcept, throwE ) import Control.Monad.Trans.State.Strict as State - ( StateT, State, runStateT, evalState, evalStateT, get, put ) + ( StateT, runStateT, evalStateT, get, put ) import Control.Monad.Trans.Writer.CPS - ( Writer, WriterT, execWriterT, runWriter, tell ) + ( WriterT, execWriterT, runWriter, tell ) -- MetaBrush import Math.Algebra.Dual @@ -134,9 +121,7 @@ import Math.Bezier.Stroke.EnvelopeEquation import Math.Differentiable ( Differentiable, DiffInterp, I ) import Math.Epsilon - ( epsilon, nearZero ) -import Math.Float.Utils - ( succFP, prevFP ) + ( epsilon ) import Math.Interval import Math.Linear import Math.Module @@ -147,8 +132,8 @@ import Math.Orientation ( Orientation(..), splineOrientation , between ) -import qualified Math.Ring as Ring import Math.Roots +import Math.Root.Isolation --import Debug.Utils @@ -198,7 +183,7 @@ instance NFData ( CachedStroke s ) where rnf _ = () discardCache :: forall crvData s. HasType ( CachedStroke s ) crvData => crvData -> ST s () -discardCache ( view ( typed @( CachedStroke s ) ) -> CachedStroke { cachedStrokeRef } ) = +discardCache ( Lens.view ( typed @( CachedStroke s ) ) -> CachedStroke { cachedStrokeRef } ) = writeSTRef cachedStrokeRef Nothing {-# INLINE invalidateCache #-} @@ -206,11 +191,11 @@ invalidateCache :: forall crvData. HasType ( CachedStroke RealWorld ) crvData => invalidateCache = runRW# \ s -> case newMutVar# Nothing s of (# _, mutVar #) -> - set ( typed @( CachedStroke RealWorld ) ) + Lens.set ( typed @( CachedStroke RealWorld ) ) ( CachedStroke $ STRef mutVar ) coords :: forall ptData. HasType ( ℝ 2 ) ptData => ptData -> ℝ 2 -coords = view typed +coords = Lens.view typed -------------------------------------------------------------------------------- @@ -258,7 +243,7 @@ computeStrokeOutline :: ) => RootSolvingAlgorithm - -> Maybe CuspFindingOptions + -> Maybe RootIsolationOptions -> FitParameters -> ( ptData -> usedParams ) -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing @@ -419,7 +404,7 @@ computeStrokeOutline rootAlgo mbCuspOptions fitParams ptParams toBrushParams bru :: crvData -> OutlineInfo -> WriterT OutlineData ( ST s ) () - updateCurveData ( view ( typed @( CachedStroke s ) ) -> CachedStroke { cachedStrokeRef } ) fwdBwd = do + updateCurveData ( Lens.view ( typed @( CachedStroke s ) ) -> CachedStroke { cachedStrokeRef } ) fwdBwd = do mbOutline <- lift ( readSTRef cachedStrokeRef ) case mbOutline of -- Cached fit data is available: use it. @@ -537,7 +522,7 @@ outlineFunction , Show ptData, Show brushParams ) => RootSolvingAlgorithm - -> Maybe CuspFindingOptions + -> Maybe RootIsolationOptions -> ( ptData -> usedParams ) -> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing -> ( forall {t} k (i :: t) @@ -945,6 +930,8 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) , offset = T $ Cubic.bezier @( T ( ℝ 2 ) ) bez t } +-------------------------------------------------------------------------------- + splineCurveFns :: forall k i . ( HasBΓ©zier k i , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) @@ -972,6 +959,86 @@ splineCurveFns co spls Bezier3To { controlPoint1 = p1, controlPoint2 = p2, curveEnd = NextPoint p3 } -> bezier3 @k @i @( ℝ 2 ) co $ Cubic.Bezier p0 p1 p2 p3 +newtype ZipSeq a = ZipSeq { getZipSeq :: Seq a } + deriving stock Functor +instance Applicative ZipSeq where + pure _ = error "only use Applicative ZipSeq with non-empty Traversable functors" + liftA2 f ( ZipSeq xs ) ( ZipSeq ys ) = ZipSeq ( Seq.zipWith f xs ys ) + {-# INLINE liftA2 #-} + +brushStrokeData :: forall k brushParams i arr + . ( arr ~ C k + , HasBΓ©zier k i, HasEnvelopeEquation k + , Differentiable k i brushParams + , HasChainRule ( I i Double ) k ( I i ( ℝ 1 ) ) + , Applicative ( D k ( ℝ 1 ) ) + + , D ( k - 2 ) ( I i ( ℝ 2 ) ) ~ D ( k - 2 ) ( ℝ 2 ) + , D ( k - 1 ) ( I i ( ℝ 2 ) ) ~ D ( k - 1 ) ( ℝ 2 ) + , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) + , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) + , D k ( I i ( ℝ 2 ) ) ~ D k ( ℝ 2 ) + , Cross ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) + , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) ) + , Show brushParams + ) + => ( I i Double -> I i ( ℝ 1 ) ) + -> ( I i ( ℝ 1 ) -> I i Double ) + -> ( I i ( ℝ 1 ) `arr` I i ( ℝ 2 ) ) + -- ^ path + -> ( I i ( ℝ 1 ) `arr` I i brushParams ) + -- ^ brush parameters + -> ( I i brushParams `arr` Spline Closed () ( I i ( ℝ 2 ) ) ) + -- ^ brush from parameters + -> ( I i ( ℝ 1 ) -> Seq ( I i ( ℝ 1 ) -> StrokeDatum k i ) ) +brushStrokeData co1 co2 path params brush = + \ t -> + let + dpath_t :: D k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) + !dpath_t = runD path t + dparams_t :: D k ( I i ( ℝ 1 ) ) ( I i brushParams ) + !dparams_t = runD params t + dbrush_params :: D k ( I i brushParams ) ( Spline Closed () ( I i ( ℝ 2 ) ) ) + !dbrush_params = runD brush $ value @( I i Double ) @k @( I i ( ℝ 1 ) ) dparams_t + splines :: Seq ( D k ( I i brushParams ) ( I i ( ℝ 1 ) `arr` I i ( ℝ 2 ) ) ) + !splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns @k @i co2 ) dbrush_params + dbrushes_t :: Seq ( I i ( ℝ 1 ) -> D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) ) + !dbrushes_t = force $ fmap ( uncurryD @k . chain @( I i Double ) @k dparams_t ) splines + -- This is the crucial use of the chain rule. + + in fmap ( mkStrokeDatum dpath_t ) dbrushes_t + where + + mkStrokeDatum :: D k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) + -> ( I i ( ℝ 1 ) -> D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) ) + -> ( I i ( ℝ 1 ) -> StrokeDatum k i ) + mkStrokeDatum dpath_t dbrush_t s = + let dbrush_t_s = dbrush_t s + dstroke = brushStroke @k dpath_t dbrush_t_s + ( ee, 𝛿E𝛿sdcdt ) = envelopeEquation @k @i co1 dstroke + in -- trace + -- ( unlines + -- [ "envelopeEquation:" + -- , " t = " ++ show t + -- , " s = " ++ show s + -- , " c = " ++ show _c + -- , " βˆ‚c/βˆ‚t = " ++ show _𝛿c𝛿t + -- , " βˆ‚c/βˆ‚s = " ++ show _𝛿c𝛿s + -- , " E = " ++ show ee + -- , " βˆ‚E/βˆ‚t = " ++ show _𝛿E𝛿t + -- , " βˆ‚E/βˆ‚s = " ++ show ee_s + -- , " dc/dt = " ++ show dcdt ] ) $ + StrokeDatum + { dpath = dpath_t + , dbrush = dbrush_t_s + , dstroke + , ee + , 𝛿E𝛿sdcdt + } + +-------------------------------------------------------------------------------- +-- Solving the envelolpe equation: root-finding. + -- | Which method to use to solve the envelope equation? data RootSolvingAlgorithm -- | Use the Newton–Raphson method. @@ -1071,134 +1138,8 @@ solveEnvelopeEquations rootAlgo _t path_t path'_t ( fwdOffset, bwdOffset ) strok let ( i0, s ) = is `divMod'` 1 in ( i0 `mod` n, s ) -newtype ZipSeq a = ZipSeq { getZipSeq :: Seq a } - deriving stock Functor -instance Applicative ZipSeq where - pure _ = error "only use Applicative ZipSeq with non-empty Traversable functors" - liftA2 f ( ZipSeq xs ) ( ZipSeq ys ) = ZipSeq ( Seq.zipWith f xs ys ) - {-# INLINE liftA2 #-} - -brushStrokeData :: forall k brushParams i arr - . ( arr ~ C k - , HasBΓ©zier k i, HasEnvelopeEquation k - , Differentiable k i brushParams - , HasChainRule ( I i Double ) k ( I i ( ℝ 1 ) ) - , Applicative ( D k ( ℝ 1 ) ) - - , D ( k - 2 ) ( I i ( ℝ 2 ) ) ~ D ( k - 2 ) ( ℝ 2 ) - , D ( k - 1 ) ( I i ( ℝ 2 ) ) ~ D ( k - 1 ) ( ℝ 2 ) - , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) - , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) - , D k ( I i ( ℝ 2 ) ) ~ D k ( ℝ 2 ) - , Cross ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) - , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) ) - , Show brushParams - ) - => ( I i Double -> I i ( ℝ 1 ) ) - -> ( I i ( ℝ 1 ) -> I i Double ) - -> ( I i ( ℝ 1 ) `arr` I i ( ℝ 2 ) ) - -- ^ path - -> ( I i ( ℝ 1 ) `arr` I i brushParams ) - -- ^ brush parameters - -> ( I i brushParams `arr` Spline Closed () ( I i ( ℝ 2 ) ) ) - -- ^ brush from parameters - -> ( I i ( ℝ 1 ) -> Seq ( I i ( ℝ 1 ) -> StrokeDatum k i ) ) -brushStrokeData co1 co2 path params brush = - \ t -> - let - dpath_t :: D k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) - !dpath_t = runD path t - dparams_t :: D k ( I i ( ℝ 1 ) ) ( I i brushParams ) - !dparams_t = runD params t - dbrush_params :: D k ( I i brushParams ) ( Spline Closed () ( I i ( ℝ 2 ) ) ) - !dbrush_params = runD brush $ value @( I i Double ) @k @( I i ( ℝ 1 ) ) dparams_t - splines :: Seq ( D k ( I i brushParams ) ( I i ( ℝ 1 ) `arr` I i ( ℝ 2 ) ) ) - !splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns @k @i co2 ) dbrush_params - dbrushes_t :: Seq ( I i ( ℝ 1 ) -> D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) ) - !dbrushes_t = force $ fmap ( uncurryD @k . chain @( I i Double ) @k dparams_t ) splines - -- This is the crucial use of the chain rule. - - in fmap ( mkStrokeDatum dpath_t ) dbrushes_t - where - - mkStrokeDatum :: D k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) - -> ( I i ( ℝ 1 ) -> D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) ) - -> ( I i ( ℝ 1 ) -> StrokeDatum k i ) - mkStrokeDatum dpath_t dbrush_t s = - let dbrush_t_s = dbrush_t s - dstroke = brushStroke @k dpath_t dbrush_t_s - ( ee, 𝛿E𝛿sdcdt ) = envelopeEquation @k @i co1 dstroke - in -- trace - -- ( unlines - -- [ "envelopeEquation:" - -- , " t = " ++ show t - -- , " s = " ++ show s - -- , " c = " ++ show _c - -- , " βˆ‚c/βˆ‚t = " ++ show _𝛿c𝛿t - -- , " βˆ‚c/βˆ‚s = " ++ show _𝛿c𝛿s - -- , " E = " ++ show ee - -- , " βˆ‚E/βˆ‚t = " ++ show _𝛿E𝛿t - -- , " βˆ‚E/βˆ‚s = " ++ show ee_s - -- , " dc/dt = " ++ show dcdt ] ) $ - StrokeDatum - { dpath = dpath_t - , dbrush = dbrush_t_s - , dstroke - , ee - , 𝛿E𝛿sdcdt - } - -------------------------------------------------------------------------------- - --- Take one interval Gauss–Seidel step for the equation \( A X = B \), --- refining the initial guess box for \( X \) into up to four (disjoint) new boxes. --- --- The boolean indicates whether the Gauss–Seidel step was a contraction. -gaussSeidel :: ( 𝕀ℝ 2, 𝕀ℝ 2 ) -- ^ columns of \( A \) - -> 𝕀ℝ 2 -- ^ \( B \) - -> 𝕀ℝ 2 -- ^ initial box \( X \) - -> [ ( 𝕀ℝ 2, Bool ) ] -gaussSeidel - ( 𝕀 ( ℝ2 a11_lo a21_lo ) ( ℝ2 a11_hi a21_hi ) - , 𝕀 ( ℝ2 a12_lo a22_lo ) ( ℝ2 a12_hi a22_hi ) ) - ( 𝕀 ( ℝ2 b1_lo b2_lo ) ( ℝ2 b1_hi b2_hi ) ) - ( 𝕀 ( ℝ2 x1_lo x2_lo ) ( ℝ2 x1_hi x2_hi ) ) - = let !a11 = 𝕀 a11_lo a11_hi - !a12 = 𝕀 a12_lo a12_hi - !a21 = 𝕀 a21_lo a21_hi - !a22 = 𝕀 a22_lo a22_hi - !b1 = 𝕀 b1_lo b1_hi - !b2 = 𝕀 b2_lo b2_hi - !x1 = 𝕀 x1_lo x1_hi - !x2 = 𝕀 x2_lo x2_hi - in nub $ do - - -- x1' = ( ( b1 - a12 * x2 ) / a11 ) ∩ x1 - x1'_0 <- ( b1 - a12 * x2 ) `extendedDivide` a11 - ( x1'@( 𝕀 x1'_lo x1'_hi ), sub_x1 ) <- x1'_0 `intersect` x1 - -- x2' = ( ( b2 - a21 * x1' ) / a22 ) ∩ x2 - x2'_0 <- ( b2 - a21 * x1' ) `extendedDivide` a22 - ( 𝕀 x2'_lo x2'_hi, sub_x2 ) <- x2'_0 `intersect` x2 - - return ( ( 𝕀 ( ℝ2 x1'_lo x2'_lo) ( ℝ2 x1'_hi x2'_hi ) ) - , sub_x1 && sub_x2 ) --- 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" - --- | 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 +-- Finding cusps in the envelope equation using root isolation. -- | Computes the brush stroke coordinates of a cusp from -- the @(t,s)@ parameter values. @@ -1222,684 +1163,51 @@ cuspCoords eqs ( i, 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) t_mid = 0.5 * ( t_lo + t_hi ) s_mid = 0.5 * ( s_lo + s_hi ) --- | Preconditioner to use with the interval Gauss–Seidel method. -data Preconditioner - = NoPreconditioning - | InverseMidJacobian - deriving stock Show - --- | A tree recording the steps taken when doing cusp finding. -data RootFindingTree d - = RootFindingLeaf (RootFindingLeaf d) - | RootFindingStep RootFindingStep [(d, RootFindingTree d)] - --- | A description of a step taken in cusp-finding. -data RootFindingStep - = GaussSeidelStep - | BisectionStep (String, Double) - | Box1Step - | Box2Step - -instance Show RootFindingStep 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)" - -data RootFindingLeaf d - = TooSmall d - deriving stock Show - -showRootFindingTree :: Box -> RootFindingTree Box -> Tree String -showRootFindingTree cand (RootFindingLeaf l) = Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ show l) [] -showRootFindingTree cand (RootFindingStep s ts) - = Node (show cand ++ " abc " ++ showArea (boxArea cand) ++ " " ++ show s) $ map (\ (c,t) -> showRootFindingTree c t) ts - -boxArea :: Box -> Double -boxArea ( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) - = abs ( t_hi - t_lo ) * abs ( s_hi - s_lo ) - -showArea :: Double -> [Char] -showArea area = "(area " ++ showFFloat (Just 6) area "" ++ ")" - -type Box = 𝕀ℝ 2 -type BoxHistory = [ NE.NonEmpty ( RootFindingStep, Box ) ] - --- | Options for the cusp-finding methods in 'findCusps' and 'findCuspsFrom'. -data CuspFindingOptions - = CuspFindingOptions - { -- | Minimum width of a box. - minWidth :: !Double - -- | Given a box and its history, return a round of cusp-finding strategies - -- to run, in sequence, on this box. - , cuspFindingAlgorithms :: !( Box -> BoxHistory -> NonEmpty RootFindingAlgorithm ) - } - -data RootFindingAlgorithm - -- | Bisection step. - = Bisection !BisectionOptions - -- | Gauss–Seidel step with the given preconditioning method. - | GaussSeidel !GaussSeidelOptions - -- | @box(1)@-consistency. - | Box1 !Box1Options - -- | @box(2)@-consistency. - | Box2 !Box2Options - --- | Options for the bisection method. -data BisectionOptions = - BisectionOptions - { canHaveSols :: !( ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 ) -> Box -> Bool ) - , fallbackBisectionDim :: !( [ ( RootFindingStep, Box ) ] -> BoxHistory -> ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 ) -> ( Int, String ) ) - } - --- | Options for the interval Gauss–Seidel method. -data GaussSeidelOptions = GaussSeidelOptions { gsPreconditioner :: !Preconditioner } - --- | Options for the @box(1)@-consistency method. -data Box1Options = Box1Options { box1EpsEq :: !Double } - --- | Options for the @box(2)@-consistency method. -data Box2Options = Box2Options { box2EpsEq :: !Double, box2LambdaMin :: !Double } - -defaultCuspFindingOptions :: CuspFindingOptions -defaultCuspFindingOptions = - CuspFindingOptions - { minWidth - , cuspFindingAlgorithms = defaultCuspFindingAlgorithms minWidth narrowAbs - } - where - minWidth = 1e-5 - narrowAbs = 5e-3 - -defaultCuspFindingAlgorithms :: Double -> Double -> Box -> BoxHistory -> NonEmpty RootFindingAlgorithm -defaultCuspFindingAlgorithms minWidth narrowAbs box history = - case history of - lastRoundBoxes : _ - -- If, in the last round of strategies, we didn't try bisection... - | any ( \case { BisectionStep {} -> False; _ -> True } . fst ) lastRoundBoxes - , (_, lastRoundFirstBox ) <- NE.last lastRoundBoxes - -- ...and the last round didn't sufficiently reduce the size of the box... - , not $ box `sufficientlySmallerThan` lastRoundFirstBox - -- ...then try bisecting the box. - -> Bisection ( defaultBisectionOptions minWidth narrowAbs box ) NE.:| [] - - -- Otherwise, do a normal round. - -- Currently: we try an interval Gauss–Seidel step followed by box(1)-consistency. - _ -> GaussSeidel defaultGaussSeidelOptions - NE.:| [ Box1 ( Box1Options { box1EpsEq = narrowAbs } ) ] - where - sufficientlySmallerThan :: Box -> Box -> Bool - sufficientlySmallerThan - ( 𝕀 ( ℝ2 t1_lo s1_lo ) ( ℝ2 t1_hi s1_hi ) ) - ( 𝕀 ( ℝ2 t2_lo s2_lo ) ( ℝ2 t2_hi s2_hi ) ) = - ( ( t1_hi - t1_lo ) - ( t2_hi - t2_lo ) > narrowAbs ) - || - ( ( s1_hi - s1_lo ) - ( s2_hi - s2_lo ) > narrowAbs ) - -defaultGaussSeidelOptions :: GaussSeidelOptions -defaultGaussSeidelOptions = GaussSeidelOptions { gsPreconditioner = InverseMidJacobian } - -defaultBisectionOptions :: Double -> Double -> Box -> BisectionOptions -defaultBisectionOptions - _minWidth - _narrowAbs - box@( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) - = BisectionOptions - { canHaveSols = - \ eqs box' -> - -- box(0)-consistency - let dat' = eqs box' - in ee_potential_zero dat' && 𝛿E𝛿sdcdt_potential_zero dat' - - -- box(1)-consistency - --not $ null $ makeBox1Consistent eqs _minWidth _narrowAbs _box' - - -- box(2)-consistency - --let ( t'', i'', s'') = makeBox2Consistent eqs _minWidth _narrowAbs 0.2 _box' - -- dat'' = ( eqs t'' `Seq.index` i'' ) s'' - --in ee_potential_zero dat'' && 𝛿E𝛿sdcdt_potential_zero dat'' - , fallbackBisectionDim = - \ _roundHist _prevRoundsHist eqs -> - let StrokeDatum { ee = D22 _ ( T ee_t ) (T ee_s ) _ _ _ - , 𝛿E𝛿sdcdt = D12 _ ( T ( T f_t ) ) ( T ( T f_s ) ) } - = eqs box - wd_t = t_hi - t_lo - wd_s = s_hi - s_lo - tJWidth = wd_t * normVI3 ee_t f_t - sJWidth = wd_s * normVI3 ee_s f_s - in if tJWidth >= sJWidth - -- bisect along dimension that maximises width(x_j) * || J_j || ... - -- ... but don't allow thin boxes - || ( wd_t >= 10 * wd_s ) - && not ( wd_s >= 10 * wd_t ) - then (0, "") - else (1, "") - } - -- | Find cusps in the envelope for values of the parameters in -- \( 0 \leqslant t, s \leqslant 1 \), using interval arithmetic. -- --- See 'findCuspsFrom' for details. +-- See 'isolateRootsIn' for details of the algorithms used. findCusps - :: CuspFindingOptions + :: RootIsolationOptions -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) - -> IntMap ( [ ( Box, RootFindingTree Box ) ], ( [ Box ], [ Box ] ) ) -findCusps opts eqs = - IntMap.fromList - [ ( i, findCuspsFrom opts eqnPiece initBox ) - | i <- [ 0 .. length ( eqs ( 𝕀 ( ℝ1 zero ) ( ℝ1 one ) ) ) - 1 ] - , let initBox = 𝕀 ( ℝ2 zero zero ) ( ℝ2 one one ) - eqnPiece ( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) = - let t = 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) - s = 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) - in ( eqs t `Seq.index` i ) s - ] - --- | Use the following algorithms using interval arithmetic in order --- to find cusps in the envelope: --- --- - interval Newton method with Gauss–Seidel step for inversion --- of the interval Jacobian, --- - coordinate-wise Newton method (@box(1)@-consistency algorithm), --- - @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. -findCuspsFrom - :: CuspFindingOptions - -> ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 ) - -> Box - -> ( [ ( Box, RootFindingTree Box ) ], ( [ Box ], [ Box ] ) ) -findCuspsFrom - ( CuspFindingOptions - { minWidth - , cuspFindingAlgorithms - } - ) - eqs initBox = runWriter $ map ( initBox , ) <$> go [] initBox - - where - - go :: BoxHistory - -> Box -- box to work on - -> Writer ( [ Box ], [ Box ] ) - [ RootFindingTree Box ] - go history - cand@( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) - | -- Check the envelope equation interval contains zero. - not ( ee_potential_zero sd) - || - -- Check the 𝛿E𝛿sdcdt interval contains zero. - not ( 𝛿E𝛿sdcdt_potential_zero sd ) - -- Box doesn't contain a solution: discard it. - = return [] - -- Box is small: stop processing it. - | t_hi - t_lo < minWidth && s_hi - s_lo < minWidth - = do let res = TooSmall cand - tell ( [ cand ], [] ) - return [ RootFindingLeaf res ] - | otherwise - = doStrategies history ( cuspFindingAlgorithms cand history ) cand - where - sd = eqs $ 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) - - -- Run a round of cusp finding strategies, then recur. - doStrategies - :: BoxHistory - -> NonEmpty RootFindingAlgorithm - -> Box - -> Writer ( [ Box ], [ Box ] ) - [ RootFindingTree Box ] - doStrategies prevRoundsHist = do_strats [] - where - do_strats :: [ ( RootFindingStep, Box ) ] - -> NE.NonEmpty RootFindingAlgorithm - -> Box - -> Writer ( [ Box ], [ Box ] ) - [ RootFindingTree Box ] - do_strats thisRoundHist ( algo NE.:| algos ) box = do - -- Run one strategy in the round. - ( step, boxes ) <- doStrategy thisRoundHist prevRoundsHist eqs minWidth algo box - case algos of - -- If there are other algorithms to run in this round, run them next. - nextAlgo : otherAlgos -> - let thisRoundHist' = ( step, box ) : thisRoundHist - in recur step ( do_strats thisRoundHist' ( nextAlgo NE.:| otherAlgos ) ) boxes - -- Otherwise, we have done one full round of strategies. - -- Recur back to the top (calling 'go'). - [] -> - let thisRoundHist' = ( step, box ) NE.:| thisRoundHist - history = thisRoundHist' : prevRoundsHist - in recur step ( go history ) boxes - - recur :: RootFindingStep - -> ( Box -> Writer ( [ Box ], [ Box ] ) [ RootFindingTree Box ] ) - -> [ Box ] - -> Writer ( [Box], [Box] ) [ RootFindingTree Box ] - recur step r cands = do - rest <- traverse ( \ c -> do { trees <- r c; return [ (c, tree) | tree <- trees ] } ) cands - return [ RootFindingStep step $ concat rest ] - - --- | Execute a cusp-finding strategy, replacing the input box with --- (hopefully smaller) output boxes. -doStrategy :: [ ( RootFindingStep, Box ) ] - -> BoxHistory - -> ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 ) - -> Double - -> RootFindingAlgorithm - -> Box - -> Writer ( [ Box ], [ Box ] ) - ( RootFindingStep, [ Box ] ) -doStrategy roundHistory previousRoundsHistory eqs minWidth algo box = - case algo of - GaussSeidel ( GaussSeidelOptions { gsPreconditioner } ) -> do - boxes <- intervalGaussSeidel eqs gsPreconditioner box - return ( GaussSeidelStep, boxes ) - Box1 ( Box1Options { box1EpsEq } ) -> - return ( Box1Step, makeBox1Consistent eqs minWidth box1EpsEq box ) - Box2 ( Box2Options { box2LambdaMin, box2EpsEq } ) -> - return ( Box2Step, [ makeBox2Consistent eqs minWidth box2EpsEq box2LambdaMin box ] ) - Bisection ( BisectionOptions { canHaveSols, fallbackBisectionDim } ) -> do - let ( boxes, whatBis ) = bisect ( canHaveSols eqs ) ( fallbackBisectionDim roundHistory previousRoundsHistory eqs ) box - return ( BisectionStep whatBis, boxes ) - --- | Bisect the given box. --- --- (The difficult part lies in determining along which dimension to bisect.) -bisect :: ( Box -> Bool ) - -- ^ how to check whether a box contains solutions - -> ( Int, String ) - -- ^ fall-back choice of dimension (and "why" we chose it) - -> Box - -> ( [ Box ], ( String, Double ) ) -bisect canHaveSols ( dim, why ) - ( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) - -- 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. - = let t_mid = 0.5 * ( t_lo + t_hi ) - s_mid = 0.5 * ( s_lo + s_hi ) - l = 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_mid s_hi ) - r = 𝕀 ( ℝ2 t_mid s_lo ) ( ℝ2 t_hi s_hi ) - d = 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_mid ) - u = 𝕀 ( ℝ2 t_lo s_mid ) ( ℝ2 t_hi s_hi ) - l_skip = not $ canHaveSols l - r_skip = not $ canHaveSols r - d_skip = not $ canHaveSols d - u_skip = not $ canHaveSols u - - ( bisGuesses, whatBis ) - | l_skip && r_skip - = ( [], ( "lr", t_mid ) ) - | d_skip && u_skip - = ( [], ( "du", s_mid ) ) - | l_skip - = ( [ r ], ( "r", t_mid ) ) - | r_skip - = ( [ l ], ( "l", t_mid ) ) - | d_skip - = ( [ u ], ( "u", s_mid ) ) - | u_skip - = ( [ d ], ( "d", s_mid ) ) - | otherwise - = let why' = case why of - "" -> "" - _ -> " (" ++ why ++ ")" - in case dim of - 0 -> ( [ l, r ], ( "t" ++ why', t_mid ) ) - 1 -> ( [ d, u ], ( "s" ++ why', s_mid ) ) - _ -> error "bisect: fall-back bisection dimension should be either 0 or 1" - in ( bisGuesses, whatBis ) - --- | Interval Newton method with Gauss–Seidel step. -intervalGaussSeidel - :: ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 ) - -- ^ equations - -> Preconditioner - -- ^ preconditioner to use for the Gauss–Seidel step - -> 𝕀ℝ 2 - -> Writer ( [ 𝕀ℝ 2 ], [ 𝕀ℝ 2 ] ) [ 𝕀ℝ 2 ] -intervalGaussSeidel - eqs - precondMethod - ts@( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) - | StrokeDatum { 𝛿E𝛿sdcdt = D12 _ ( T ( T f_t ) ) ( T ( T f_s ) ) } - <- eqs ts - , StrokeDatum { ee = D22 _ee_mid _ _ _ _ _ - , 𝛿E𝛿sdcdt = D12 ( T f_mid ) ( T ( T _f_t_mid ) ) ( T ( T _f_s_mid ) ) } - <- eqs ts_mid - = let -- Interval Newton method: take one Gauss–Seidel step - -- for the equation f'(x) v = - f(x_mid), - -- where f = 𝛿E/𝛿s * dc/dt - ( a, b ) = precondition precondMethod - ( midI f_t, midI f_s ) - ( f_t, f_s ) ( neg f_mid ) - --(a, b) - -- | 𝕀 (ℝ1 ee_lo) (ℝ1 ee_hi) <- ee_mid - -- , 𝕀 (ℝ1 ee_t_lo) (ℝ1 ee_t_hi) <- ee_t - -- , 𝕀 (ℝ1 ee_s_lo) (ℝ1 ee_s_hi) <- ee_s - -- , 𝕀 (ℝ2 fx_lo fy_lo) (ℝ2 fx_hi fy_hi) <- f_mid - -- , 𝕀 (ℝ2 f_tx_lo f_ty_lo) (ℝ2 f_tx_hi f_ty_hi) <- f_t - -- , 𝕀 (ℝ2 f_sx_lo f_sy_lo) (ℝ2 f_sx_hi f_sy_hi) <- f_s - -- = ( ( 𝕀 (ℝ2 f_tx_lo ee_t_lo) (ℝ2 f_tx_hi ee_t_hi) - -- , 𝕀 (ℝ2 f_sx_lo ee_s_lo) (ℝ2 f_sx_hi ee_s_hi) - -- ) - -- , neg $ 𝕀 (ℝ2 fx_lo ee_lo) (ℝ2 fx_hi ee_hi) - -- ) - - gsGuesses = gaussSeidel a b ( unT $ T ts ^-^ T ts_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 ( mkGuess . fst ) ) ( map ( mkGuess . fst ) ) - $ partition snd gsGuesses - in do tell ([], done) - return todo - where - t_mid = 0.5 * ( t_lo + t_hi ) - s_mid = 0.5 * ( s_lo + s_hi ) - ts_mid = singleton ( ℝ2 t_mid s_mid ) - mkGuess ts_g = unT $ T ts_g ^+^ T ts_mid - neg ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) - = let !( 𝕀 x'_lo x'_hi ) = negate $ 𝕀 x_lo x_hi - !( 𝕀 y'_lo y'_hi ) = negate $ 𝕀 y_lo y_hi - in 𝕀 ( ℝ2 x'_lo y'_lo ) ( ℝ2 x'_hi y'_hi ) - --- | 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 - :: ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 ) - -> Double -> Double - -> Box -> [ Box ] -makeBox1Consistent eqs minWidth epsEq x = - ( `State.evalState` False ) $ - pipeFunctionsWhileTrue (allNarrowingOperators epsEq minWidth 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 - :: ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 ) - -> Double -> Double -> Double - -> Box -> Box -makeBox2Consistent eqs minWidth epsEq lambdaMin x0 = ( `State.evalState` False ) $ doLoop 0.25 x0 - where - doLoop :: Double -> Box -> State Bool Box - doLoop lambda x = do - x' <- boundConsistency get_t set_t lambda x - x'' <- boundConsistency get_s set_s lambda x' - modified <- State.get - let lambda' = if modified then lambda else 0.5 * lambda - if lambda' < lambdaMin - then return x'' - else do { State.put False ; doLoop lambda' x'' } - - boundConsistency :: ( Box -> 𝕀 Double ) - -> ( 𝕀 Double -> Box -> Box ) - -> Double -> Box -> State Bool Box - boundConsistency getter setter lambda box = do - let x@( 𝕀 x_inf x_sup ) = getter box - c1 = ( 1 - lambda ) * x_inf + lambda * x_sup - c2 = lambda * x_inf + ( 1 - lambda ) * x_sup - x'_inf = - case makeBox1Consistent eqs minWidth epsEq ( setter ( 𝕀 x_inf c1 ) box ) of - [] -> c1 - x's -> minimum $ map ( inf . getter ) x's - x'_sup = - case makeBox1Consistent eqs minWidth epsEq ( 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' >= epsEq ) $ - State.put True - return $ setter x' box - -ee_potential_zero :: StrokeDatum 3 𝕀 -> Bool -ee_potential_zero dat = - inf ( _D22_v $ ee dat ) <= ℝ1 0 - && sup ( _D22_v $ ee dat ) >= ℝ1 0 -𝛿E𝛿sdcdt_potential_zero :: StrokeDatum 3 𝕀 -> Bool -𝛿E𝛿sdcdt_potential_zero dat = - cmpℝ2 (<=) ( inf $ unT $ _D12_v $ 𝛿E𝛿sdcdt dat ) ( ℝ2 0 0 ) - && cmpℝ2 (>=) ( sup $ unT $ _D12_v $ 𝛿E𝛿sdcdt dat ) ( ℝ2 0 0 ) - -zero, one :: Double -zero = 1e-6 -one = 1 - zero -{-# INLINE zero #-} -{-# INLINE one #-} - -precondition :: Preconditioner - -> ( 𝕀ℝ 2, 𝕀ℝ 2 ) - -> ( 𝕀ℝ 2, 𝕀ℝ 2 ) - -> 𝕀ℝ 2 - -> ( ( 𝕀ℝ 2, 𝕀ℝ 2 ), 𝕀ℝ 2 ) -precondition meth jac_mid a@( a1, a2 ) b = - case meth of - NoPreconditioning - -> ( a, b ) - InverseMidJacobian - | ( 𝕀 ( ℝ2 a11_lo a21_lo ) ( ℝ2 a11_hi a21_hi ) - , 𝕀 ( ℝ2 a12_lo a22_lo ) ( ℝ2 a12_hi a22_hi ) ) <- jac_mid - , let !a11 = 0.5 * ( a11_lo + a11_hi ) - !a12 = 0.5 * ( a12_lo + a12_hi ) - !a21 = 0.5 * ( a21_lo + a21_hi ) - !a22 = 0.5 * ( a22_lo + a22_hi ) - !d = a11 * a22 - a12 * a21 - , not ( nearZero d ) - , let !precond = ( ℝ2 a22 -a12, ℝ2 -a21 a11 ) - !inv = recip d - f x = scale inv $ matMulVec precond x - -> ( ( f a1, f a2 ), f b ) - | otherwise - -> ( a, b ) - -scale :: Double -> 𝕀ℝ 2 -> 𝕀ℝ 2 -scale s ( 𝕀 ( ℝ2 a1_lo a2_lo ) ( ℝ2 a1_hi a2_hi ) ) - = 𝕀 ( ℝ2 b1_lo b2_lo ) ( ℝ2 b1_hi b2_hi ) - where - 𝕀 b1_lo b1_hi = scaleInterval s ( 𝕀 a1_lo a1_hi ) - 𝕀 b2_lo b2_hi = scaleInterval s ( 𝕀 a2_lo a2_hi ) - -matMulVec :: ( ℝ 2, ℝ 2 ) -> 𝕀ℝ 2 -> 𝕀ℝ 2 -matMulVec ( ℝ2 a11 a21, ℝ2 a12 a22 ) ( 𝕀 ( ℝ2 u_lo v_lo ) ( ℝ2 u_hi v_hi ) ) - = 𝕀 ( ℝ2 u'_lo v'_lo ) ( ℝ2 u'_hi v'_hi ) - where - u = 𝕀 u_lo u_hi - v = 𝕀 v_lo v_hi - 𝕀 u'_lo u'_hi = scaleInterval a11 u + scaleInterval a12 v - 𝕀 v'_lo v'_hi = scaleInterval a21 u + scaleInterval a22 v - -cmpℝ2 :: ( Double -> Double -> Bool ) -> ℝ 2 -> ℝ 2 -> Bool -cmpℝ2 cmp ( ℝ2 x1 y1 ) ( ℝ2 x2 y2 ) - = cmp x1 x2 && cmp y1 y2 - -midI :: 𝕀ℝ 2 -> 𝕀ℝ 2 -midI ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = - singleton $ ℝ2 ( 0.5 * ( x_lo + x_hi ) ) ( 0.5 * ( y_lo + y_hi ) ) ---width :: 𝕀ℝ 1 -> Double ---width ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) = hi - lo ---normI :: 𝕀ℝ 1 -> Double ---normI ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) = sqrt $ sup $ ( 𝕀 lo hi ) Ring.^ 2 ---normVI :: 𝕀ℝ 2 -> Double ---normVI ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = --- sqrt $ sup $ ( 𝕀 x_lo x_hi ) Ring.^ 2 + ( 𝕀 y_lo y_hi ) Ring.^ 2 -normVI3 :: 𝕀ℝ 1 -> 𝕀ℝ 2 -> Double -normVI3 ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = - sqrt $ max ( abs lo ) ( abs hi ) Ring.^ 2 - + max ( abs x_lo ) ( abs x_hi ) Ring.^ 2 - + max ( abs y_lo ) ( abs y_hi ) Ring.^ 2 - --- 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 eps_equal eps_bisection 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 delta <- f_x_left `extendedDivide` f'_x --f'_x_left - let x_new = x_left - delta - map fst $ x_new `intersect` x - in - if | null x's - -> [] - | ( width x - sum ( map width x's ) ) < eps_equal - -- TODO: do a check with the relative reduction in size? - -> x's - | otherwise - -> do - x' <- x's - if sup x' - inf x' < eps_bisection - then return x' - else left_narrow =<< bisectI x' - --- TODO: de-duplicate with 'leftNarrow'? -rightNarrow :: Double - -> Double - -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) ) - -> 𝕀 Double - -> [ 𝕀 Double ] -rightNarrow eps_equal eps_bisection 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 delta <- f_x_right `extendedDivide` f'_x --f'_x_right - let x_new = x_right - delta - map fst $ x_new `intersect` x - in - if | null x's - -> [] - | ( width x - sum ( map width x's ) ) < eps_equal - -- TODO: do a check with the relative reduction in size? - -> x's - | otherwise - -> do - x' <- x's - if sup x' - inf x' < eps_bisection - 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 - :: Double - -> Double - -> ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 ) - -> [ Box -> State Bool [ Box ] ] -allNarrowingOperators eps_eq eps_bis eqs = - [ \ cand -> - let newCands = map ( `setter` cand ) $ narrowFn eps_eq eps_bis ( \ x0 -> ff' $ setter x0 cand ) ( getter cand ) - in do - when ( ( width ( getter cand ) - sum ( map ( width . getter ) newCands ) ) > eps_eq ) $ - -- TODO: do a check with the relative reduction in size? - State.put True - return newCands - | narrowFn <- [ leftNarrow, rightNarrow ] - , ( ( getter, setter ), ff' ) <- - [ ( ( get_t, set_t ), ff' ) | ff' <- [ ff'_t, g1g1'_t, g2g2'_t ] ] - ++ - [ ( ( get_s, set_s ), ff' ) | ff' <- [ ff'_s, g1g1'_s, g2g2'_s ] ] + -> IntMap ( [ ( Box, RootIsolationTree Box ) ], ( [ Box ], [ Box ] ) ) +findCusps opts boxStrokeData = + findCuspsIn opts boxStrokeData $ + IntMap.fromList + [ ( i, [ 𝕀 ( ℝ2 zero zero ) ( ℝ2 one one ) ] ) + | i <- [ 0 .. length ( boxStrokeData ( 𝕀 ( ℝ1 zero ) ( ℝ1 one ) ) ) - 1 ] ] where - ff'_t ts = - let D22 { _D22_v = 𝕀 ( ℝ1 ee_inf ) ( ℝ1 ee_sup ) - , _D22_dx = T ( 𝕀 ( ℝ1 ee'_inf ) ( ℝ1 ee'_sup ) ) - } = ee $ eqs ts - in ( 𝕀 ee_inf ee_sup, 𝕀 ee'_inf ee'_sup ) - ff'_s ts = - let D22 { _D22_v = 𝕀 ( ℝ1 ee_inf ) ( ℝ1 ee_sup ) - , _D22_dy = T ( 𝕀 ( ℝ1 ee'_inf ) ( ℝ1 ee'_sup ) ) - } = ee $ eqs ts - in ( 𝕀 ee_inf ee_sup, 𝕀 ee'_inf ee'_sup ) - g1g1'_t ts = - let D12 { _D12_v = T ( 𝕀 ( ℝ2 x_inf _ ) ( ℝ2 x_sup _ ) ) - , _D12_dx = T ( T ( 𝕀 ( ℝ2 x'_inf _ ) ( ℝ2 x'_sup _ ) ) ) - } = 𝛿E𝛿sdcdt $ eqs ts - in ( 𝕀 x_inf x_sup, 𝕀 x'_inf x'_sup ) - g1g1'_s ts = - let D12 { _D12_v = T ( 𝕀 ( ℝ2 x_inf _ ) ( ℝ2 x_sup _ ) ) - , _D12_dy = T ( T ( 𝕀 ( ℝ2 x'_inf _ ) ( ℝ2 x'_sup _ ) ) ) - } = 𝛿E𝛿sdcdt $ eqs ts - in ( 𝕀 x_inf x_sup, 𝕀 x'_inf x'_sup ) - g2g2'_t ts = - let D12 { _D12_v = T ( 𝕀 ( ℝ2 _ y_inf ) ( ℝ2 _ y_sup ) ) - , _D12_dx = T ( T ( 𝕀 ( ℝ2 _ y'_inf ) ( ℝ2 _ y'_sup ) ) ) - } = 𝛿E𝛿sdcdt $ eqs ts - in ( 𝕀 y_inf y_sup, 𝕀 y'_inf y'_sup ) - g2g2'_s ts = - let D12 { _D12_v = T ( 𝕀 ( ℝ2 _ y_inf ) ( ℝ2 _ y_sup ) ) - , _D12_dy = T ( T ( 𝕀 ( ℝ2 _ y'_inf ) ( ℝ2 _ y'_sup ) ) ) - } = 𝛿E𝛿sdcdt $ eqs ts - in ( 𝕀 y_inf y_sup, 𝕀 y'_inf y'_sup ) + zero, one :: Double + zero = 1e-6 + one = 1 - zero + {-# INLINE zero #-} + {-# INLINE one #-} --- TODO: replace this with Fin n and representable please -get_t, get_s :: Box -> 𝕀 Double -get_t ( 𝕀 ( ℝ2 t_inf _ ) ( ℝ2 t_sup _ ) ) = 𝕀 t_inf t_sup -get_s ( 𝕀 ( ℝ2 _ s_inf ) ( ℝ2 _ s_sup ) ) = 𝕀 s_inf s_sup - -set_t, set_s :: 𝕀 Double -> Box -> Box -set_t ( 𝕀 t_inf t_sup ) ( 𝕀 ( ℝ2 _ s_inf ) ( ℝ2 _ s_sup ) ) = 𝕀 ( ℝ2 t_inf s_inf ) ( ℝ2 t_sup s_sup ) -set_s ( 𝕀 s_inf s_sup ) ( 𝕀 ( ℝ2 t_inf _ ) ( ℝ2 t_sup _ ) ) = 𝕀 ( ℝ2 t_inf s_inf ) ( ℝ2 t_sup s_sup ) +-- | Like 'findCusps', but parametrised over the initial boxes for the +-- root isolation method. +findCuspsIn + :: RootIsolationOptions + -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) + -> IntMap [ Box ] + -> IntMap ( [ ( Box, RootIsolationTree Box ) ], ( [ Box ], [ Box ] ) ) +findCuspsIn opts boxStrokeData initBoxes = + IntMap.mapWithKey ( \ i -> foldMap ( isolateRootsIn opts ( eqnPiece i ) ) ) initBoxes + where + eqnPiece i ( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) = + let t = 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) + s = 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) + -- This conversion is quite painful, but hey... + StrokeDatum + { ee = + D22 { _D22_v = 𝕀 ( ℝ1 ee_lo ) ( ℝ1 ee_hi ) + , _D22_dx = T ( 𝕀 ( ℝ1 ee_t_lo ) ( ℝ1 ee_t_hi ) ) + , _D22_dy = T ( 𝕀 ( ℝ1 ee_s_lo ) ( ℝ1 ee_s_hi ) ) } + , 𝛿E𝛿sdcdt = + D12 { _D12_v = T ( 𝕀 ( ℝ2 f_lo g_lo ) ( ℝ2 f_hi g_hi ) ) + , _D12_dx = T ( T ( 𝕀 ( ℝ2 f_t_lo g_t_lo ) ( ℝ2 f_t_hi g_t_hi ) ) ) + , _D12_dy = T ( T ( 𝕀 ( ℝ2 f_s_lo g_s_lo ) ( ℝ2 f_s_hi g_s_hi ) ) ) } + } = ( boxStrokeData t `Seq.index` i ) s + in D12 ( 𝕀 ( ℝ3 ee_lo f_lo g_lo ) ( ℝ3 ee_hi f_hi g_hi ) ) + ( T $ 𝕀 ( ℝ3 ee_t_lo f_t_lo g_t_lo ) ( ℝ3 ee_t_hi f_t_hi g_t_hi ) ) + ( T $ 𝕀 ( ℝ3 ee_s_lo f_s_lo g_s_lo ) ( ℝ3 ee_s_hi f_s_hi g_s_hi ) ) diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs b/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs index b05f051..9c69a6c 100644 --- a/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs +++ b/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs @@ -317,9 +317,9 @@ evaluateCubic bez t = let inf_bez = fmap inf bez sup_bez = fmap sup bez mins = fmap (Cubic.bezier @( T Double ) inf_bez) - $ inf t :| ( sup t : filter ( inside t ) ( Cubic.extrema inf_bez ) ) + $ inf t :| ( sup t : filter ( `inside` t ) ( Cubic.extrema inf_bez ) ) maxs = fmap (Cubic.bezier @( T Double ) sup_bez) - $ inf t :| ( sup t : filter ( inside t ) ( Cubic.extrema sup_bez ) ) + $ inf t :| ( sup t : filter ( `inside` t ) ( Cubic.extrema sup_bez ) ) in 𝕀 ( minimum mins ) ( maximum maxs ) -- | Evaluate a quadratic BΓ©zier curve, when both the coefficients and the @@ -330,9 +330,9 @@ evaluateQuadratic bez t = let inf_bez = fmap inf bez sup_bez = fmap sup bez mins = fmap (Quadratic.bezier @( T Double ) inf_bez) - $ inf t :| ( sup t : filter ( inside t ) ( Quadratic.extrema inf_bez ) ) + $ inf t :| ( sup t : filter ( `inside` t ) ( Quadratic.extrema inf_bez ) ) maxs = fmap (Quadratic.bezier @( T Double ) sup_bez) - $ inf t :| ( sup t : filter ( inside t ) ( Quadratic.extrema sup_bez ) ) + $ inf t :| ( sup t : filter ( `inside` t ) ( Quadratic.extrema sup_bez ) ) in 𝕀 ( minimum mins ) ( maximum maxs ) {- diff --git a/brush-strokes/src/lib/Math/Interval.hs b/brush-strokes/src/lib/Math/Interval.hs index b003d06..1ca0582 100644 --- a/brush-strokes/src/lib/Math/Interval.hs +++ b/brush-strokes/src/lib/Math/Interval.hs @@ -59,8 +59,8 @@ width ( 𝕀 lo hi ) = hi - lo nonDecreasing :: ( a -> b ) -> 𝕀 a -> 𝕀 b nonDecreasing f ( 𝕀 lo hi ) = 𝕀 ( f lo ) ( f hi ) -inside :: Ord a => 𝕀 a -> a -> Bool -inside ( 𝕀 lo hi ) x = x >= lo && x <= hi +inside :: Ord a => a -> 𝕀 a -> Bool +inside x ( 𝕀 lo hi ) = x >= lo && x <= hi deriving via ViaAbelianGroup ( T ( 𝕀 Double ) ) instance Semigroup ( T ( 𝕀 Double ) ) diff --git a/brush-strokes/src/lib/Math/Linear.hs b/brush-strokes/src/lib/Math/Linear.hs index e0260a2..d49c687 100644 --- a/brush-strokes/src/lib/Math/Linear.hs +++ b/brush-strokes/src/lib/Math/Linear.hs @@ -12,7 +12,7 @@ module Math.Linear , ℝ(..), T(.., V2, V3, V4) , Fin(..), MFin(..) , RepDim, RepresentableQ(..) - , Representable(..), injection, projection + , Representable(..), set, injection, projection , Vec(..), (!), find, zipIndices, unzip ) where diff --git a/brush-strokes/src/lib/Math/Linear/Internal.hs b/brush-strokes/src/lib/Math/Linear/Internal.hs index 8830247..f0b1bbf 100644 --- a/brush-strokes/src/lib/Math/Linear/Internal.hs +++ b/brush-strokes/src/lib/Math/Linear/Internal.hs @@ -6,7 +6,9 @@ module Math.Linear.Internal , Fin(..), MFin(..) , RepDim , RepresentableQ(..) - , Representable(..), projection, injection + , Representable(..) + , set + , projection, injection ) where @@ -111,6 +113,12 @@ class Representable r v | v -> r where tabulate :: ( Fin ( RepDim v ) -> r ) -> v index :: v -> Fin ( RepDim v ) -> r +{-# INLINEABLE set #-} +set :: Representable r v => Fin ( RepDim v ) -> r -> v -> v +set i r u = tabulate \ j -> + if i == j + then r + else index u j projection :: ( Representable r u, Representable r v ) => ( Fin ( RepDim v ) -> Fin ( RepDim u ) ) @@ -138,20 +146,24 @@ instance RepresentableQ Double ( ℝ 0 ) where instance RepresentableQ Double ( ℝ 1 ) where tabulateQ f = [|| ℝ1 $$( f ( Fin 1 ) ) ||] - indexQ p _ = [|| unℝ1 $$p ||] + indexQ p = \ case + Fin 1 -> [|| unℝ1 $$p ||] + Fin i -> error $ "invalid index for ℝ 1: " ++ show i instance RepresentableQ Double ( ℝ 2 ) where tabulateQ f = [|| ℝ2 $$( f ( Fin 1 ) ) $$( f ( Fin 2 ) ) ||] indexQ p = \ case Fin 1 -> [|| _ℝ2_x $$p ||] - _ -> [|| _ℝ2_y $$p ||] + Fin 2 -> [|| _ℝ2_y $$p ||] + Fin i -> error $ "invalid index for ℝ 2: " ++ show i instance RepresentableQ Double ( ℝ 3 ) where tabulateQ f = [|| ℝ3 $$( f ( Fin 1 ) ) $$( f ( Fin 2 ) ) $$( f ( Fin 3 ) ) ||] indexQ p = \ case Fin 1 -> [|| _ℝ3_x $$p ||] Fin 2 -> [|| _ℝ3_y $$p ||] - _ -> [|| _ℝ3_z $$p ||] + Fin 3 -> [|| _ℝ3_z $$p ||] + Fin i -> error $ "invalid index for ℝ 3: " ++ show i instance RepresentableQ Double ( ℝ 4 ) where tabulateQ f = [|| ℝ4 $$( f ( Fin 1 ) ) $$( f ( Fin 2 ) ) $$( f ( Fin 3 ) ) $$( f ( Fin 4 ) ) ||] @@ -159,7 +171,8 @@ instance RepresentableQ Double ( ℝ 4 ) where Fin 1 -> [|| _ℝ4_x $$p ||] Fin 2 -> [|| _ℝ4_y $$p ||] Fin 3 -> [|| _ℝ4_z $$p ||] - _ -> [|| _ℝ4_w $$p ||] + Fin 4 -> [|| _ℝ4_w $$p ||] + Fin i -> error $ "invalid index for ℝ 4: " ++ show i instance Representable Double ( ℝ 0 ) where tabulate _ = ℝ0 @@ -170,7 +183,9 @@ instance Representable Double ( ℝ 0 ) where instance Representable Double ( ℝ 1 ) where tabulate f = ℝ1 ( f ( Fin 1 ) ) {-# INLINE tabulate #-} - index p _ = unℝ1 p + index p = \ case + Fin 1 -> unℝ1 p + Fin i -> error $ "invalid index for ℝ 1: " ++ show i {-# INLINE index #-} instance Representable Double ( ℝ 2 ) where @@ -178,7 +193,8 @@ instance Representable Double ( ℝ 2 ) where {-# INLINE tabulate #-} index p = \ case Fin 1 -> _ℝ2_x p - _ -> _ℝ2_y p + Fin 2 -> _ℝ2_y p + Fin i -> error $ "invalid index for ℝ 2: " ++ show i {-# INLINE index #-} instance Representable Double ( ℝ 3 ) where @@ -187,7 +203,8 @@ instance Representable Double ( ℝ 3 ) where index p = \ case Fin 1 -> _ℝ3_x p Fin 2 -> _ℝ3_y p - _ -> _ℝ3_z p + Fin 3 -> _ℝ3_z p + Fin i -> error $ "invalid index for ℝ 3: " ++ show i {-# INLINE index #-} instance Representable Double ( ℝ 4 ) where @@ -197,5 +214,6 @@ instance Representable Double ( ℝ 4 ) where Fin 1 -> _ℝ4_x p Fin 2 -> _ℝ4_y p Fin 3 -> _ℝ4_z p - _ -> _ℝ4_w p + Fin 4 -> _ℝ4_w p + Fin i -> error $ "invalid index for ℝ 4: " ++ show i {-# INLINE index #-} diff --git a/brush-strokes/src/lib/Math/Root/Isolation.hs b/brush-strokes/src/lib/Math/Root/Isolation.hs new file mode 100644 index 0000000..c603b08 --- /dev/null +++ b/brush-strokes/src/lib/Math/Root/Isolation.hs @@ -0,0 +1,746 @@ +module Math.Root.Isolation + ( -- * Root isolation using interval arithmetic + Box, BoxHistory + , isolateRootsIn + + -- ** Configuration of the root isolation methods + , RootIsolationOptions(..), defaultRootIsolationOptions + , RootIsolationAlgorithm(..), defaultRootIsolationAlgorithms + + -- *** 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 + , RootIsolationStep(..) + , RootIsolationLeaf(..) + ) + where + +-- base +import Control.Arrow + ( first ) +import Control.Monad + ( when ) +import Data.Bifunctor + ( Bifunctor(bimap) ) +import Data.List + ( nub, partition) +import Data.List.NonEmpty + ( NonEmpty ) +import qualified Data.List.NonEmpty as NE + ( NonEmpty(..), last ) +import Numeric + ( showFFloat ) + +-- containers +import Data.Tree + ( Tree(..) ) + +-- 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, D1𝔸2(..) ) +import Math.Epsilon + ( nearZero ) +import Math.Float.Utils + ( succFP, prevFP ) +import Math.Interval +import Math.Linear +import Math.Module + ( Module(..) ) +import qualified Math.Ring as Ring + +--import Debug.Utils + +-------------------------------------------------------------------------------- + +-- | A tree recording the steps taken when doing cusp finding. +data RootIsolationTree d + = RootIsolationLeaf (RootIsolationLeaf d) + | RootIsolationStep RootIsolationStep [(d, RootIsolationTree d)] + +-- | A description of a step taken in cusp-finding. +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)" + +data RootIsolationLeaf d + = TooSmall d + deriving stock Show + +showRootIsolationTree :: Box -> RootIsolationTree Box -> Tree String +showRootIsolationTree cand (RootIsolationLeaf l) = Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ show l) [] +showRootIsolationTree cand (RootIsolationStep s ts) + = Node (show cand ++ " abc " ++ showArea (boxArea cand) ++ " " ++ show s) $ map (\ (c,t) -> showRootIsolationTree c t) ts + +boxArea :: Box -> Double +boxArea ( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) + = abs ( t_hi - t_lo ) * abs ( s_hi - s_lo ) + +showArea :: Double -> String +showArea area = "(area " ++ showFFloat (Just 6) area "" ++ ")" + +type Box = 𝕀ℝ 2 +type BoxHistory = [ NE.NonEmpty ( RootIsolationStep, Box ) ] + +-- | Options for the root isolation methods in 'findCusps' and 'isolateRootsIn'. +data RootIsolationOptions + = RootIsolationOptions + { -- | Minimum width of a box. + minWidth :: !Double + -- | Given a box and its history, return a round of cusp-finding strategies + -- to run, in sequence, on this box. + , cuspFindingAlgorithms :: !( Box -> BoxHistory -> NonEmpty RootIsolationAlgorithm ) + } + +data RootIsolationAlgorithm + -- | Bisection step. + = Bisection !BisectionOptions + -- | Gauss–Seidel step with the given preconditioning method. + | GaussSeidel !GaussSeidelOptions + -- | @box(1)@-consistency. + | Box1 !Box1Options + -- | @box(2)@-consistency. + | Box2 !Box2Options + +-- | Options for the bisection method. +data BisectionOptions = + BisectionOptions + { canHaveSols :: !( ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) ) -> Box -> Bool ) + , fallbackBisectionDim :: !( [ ( RootIsolationStep, Box ) ] -> BoxHistory -> ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) ) -> ( Int, String ) ) + } + +-- | Options for the interval Gauss–Seidel method. +data GaussSeidelOptions = + GaussSeidelOptions + { gsPreconditioner :: !Preconditioner } + +-- | Options for the @box(1)@-consistency method. +data Box1Options = + Box1Options + { box1EpsEq :: !Double } + +-- | Options for the @box(2)@-consistency method. +data Box2Options = + Box2Options + { box2EpsEq :: !Double + , box2LambdaMin :: !Double + } + +defaultRootIsolationOptions :: RootIsolationOptions +defaultRootIsolationOptions = + RootIsolationOptions + { minWidth + , cuspFindingAlgorithms = defaultRootIsolationAlgorithms minWidth narrowAbs + } + where + minWidth = 1e-5 + narrowAbs = 5e-3 + +defaultRootIsolationAlgorithms :: Double -> Double -> Box -> BoxHistory -> NonEmpty RootIsolationAlgorithm +defaultRootIsolationAlgorithms minWidth narrowAbs box history = + case history of + lastRoundBoxes : _ + -- If, in the last round of strategies, we didn't try bisection... + | any ( \case { BisectionStep {} -> False; _ -> True } . fst ) lastRoundBoxes + , (_, lastRoundFirstBox ) <- NE.last lastRoundBoxes + -- ...and the last round didn't sufficiently reduce the size of the box... + , not $ box `sufficientlySmallerThan` lastRoundFirstBox + -- ...then try bisecting the box. + -> Bisection ( defaultBisectionOptions minWidth narrowAbs box ) NE.:| [] + + -- Otherwise, do a normal round. + -- Currently: we try an interval Gauss–Seidel step followed by box(1)-consistency. + _ -> GaussSeidel defaultGaussSeidelOptions + NE.:| [ Box1 ( Box1Options { box1EpsEq = narrowAbs } ) ] + where + sufficientlySmallerThan :: Box -> Box -> Bool + sufficientlySmallerThan + ( 𝕀 ( ℝ2 t1_lo s1_lo ) ( ℝ2 t1_hi s1_hi ) ) + ( 𝕀 ( ℝ2 t2_lo s2_lo ) ( ℝ2 t2_hi s2_hi ) ) = + ( ( t1_hi - t1_lo ) - ( t2_hi - t2_lo ) > narrowAbs ) + || + ( ( s1_hi - s1_lo ) - ( s2_hi - s2_lo ) > narrowAbs ) + +defaultGaussSeidelOptions :: GaussSeidelOptions +defaultGaussSeidelOptions = GaussSeidelOptions { gsPreconditioner = InverseMidJacobian } + +defaultBisectionOptions :: Double -> Double -> Box -> BisectionOptions +defaultBisectionOptions + _minWidth + _narrowAbs + box@( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) + = BisectionOptions + { canHaveSols = + \ eqs box' -> + -- box(0)-consistency + let D12 iRange _ _ = eqs box' + in ℝ3 0 0 0 `inside` iRange + + -- box(1)-consistency + --not $ null $ makeBox1Consistent eqs _minWidth _narrowAbs box' + + -- box(2)-consistency + --let box'' = makeBox2Consistent eqs _minWidth _narrowAbs 0.2 box' + -- D12 iRange'' = eqs box'' + --in ℝ3 0 0 0 `inside` iRange'' + , fallbackBisectionDim = + \ _roundHist _prevRoundsHist eqs -> + let D12 _ ( T f_t ) ( T f_s ) = eqs box + wd_t = t_hi - t_lo + wd_s = s_hi - s_lo + tJWidth = wd_t * normVI3 f_t + sJWidth = wd_s * normVI3 f_s + in if tJWidth >= sJWidth + -- bisect along dimension that maximises width(x_j) * || J_j || ... + -- ... but don't allow thin boxes + || ( wd_t >= 10 * wd_s ) + && not ( wd_s >= 10 * wd_t ) + then (0, "") + else (1, "") + } + +-- | Use the following algorithms using interval arithmetic in order +-- to isolate roots: +-- +-- - interval Newton method with Gauss–Seidel step for inversion +-- of the interval Jacobian, +-- - coordinate-wise Newton method (@box(1)@-consistency algorithm), +-- - @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. +isolateRootsIn + :: RootIsolationOptions + -> ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) ) + -> Box + -> ( [ ( Box, RootIsolationTree Box ) ], ( [ Box ], [ Box ] ) ) +isolateRootsIn + ( RootIsolationOptions + { minWidth + , cuspFindingAlgorithms + } + ) + eqs initBox = runWriter $ map ( initBox , ) <$> go [] initBox + + where + + go :: BoxHistory + -> Box -- box to work on + -> Writer ( [ Box ], [ Box ] ) + [ RootIsolationTree Box ] + go history + cand@( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) + | -- Check the envelope equation interval contains zero. + not $ ( ℝ3 0 0 0 `inside` iRange ) + -- Box doesn't contain a solution: discard it. + = return [] + -- Box is small: stop processing it. + | t_hi - t_lo < minWidth && s_hi - s_lo < minWidth + = do let res = TooSmall cand + tell ( [ cand ], [] ) + return [ RootIsolationLeaf res ] + | otherwise + = doStrategies history ( cuspFindingAlgorithms cand history ) cand + where + D12 iRange _ _ = eqs $ 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) + + -- Run a round of cusp finding strategies, then recur. + doStrategies + :: BoxHistory + -> NonEmpty RootIsolationAlgorithm + -> Box + -> Writer ( [ Box ], [ Box ] ) + [ RootIsolationTree Box ] + doStrategies prevRoundsHist = do_strats [] + where + do_strats :: [ ( RootIsolationStep, Box ) ] + -> NE.NonEmpty RootIsolationAlgorithm + -> Box + -> Writer ( [ Box ], [ Box ] ) + [ RootIsolationTree Box ] + do_strats thisRoundHist ( algo NE.:| algos ) box = do + -- Run one strategy in the round. + ( step, boxes ) <- doStrategy thisRoundHist prevRoundsHist eqs minWidth algo box + case algos of + -- If there are other algorithms to run in this round, run them next. + nextAlgo : otherAlgos -> + let thisRoundHist' = ( step, box ) : thisRoundHist + in recur step ( do_strats thisRoundHist' ( nextAlgo NE.:| otherAlgos ) ) boxes + -- Otherwise, we have done one full round of strategies. + -- Recur back to the top (calling 'go'). + [] -> + let thisRoundHist' = ( step, box ) NE.:| thisRoundHist + history = thisRoundHist' : prevRoundsHist + in recur step ( go history ) boxes + + recur :: RootIsolationStep + -> ( Box -> Writer ( [ Box ], [ Box ] ) [ RootIsolationTree Box ] ) + -> [ Box ] + -> Writer ( [Box], [Box] ) [ RootIsolationTree Box ] + recur step r cands = do + rest <- traverse ( \ c -> do { trees <- r c; return [ (c, tree) | tree <- trees ] } ) cands + return [ RootIsolationStep step $ concat rest ] + + +-- | Execute a cusp-finding strategy, replacing the input box with +-- (hopefully smaller) output boxes. +doStrategy :: [ ( RootIsolationStep, Box ) ] + -> BoxHistory + -> ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) ) + -> Double + -> RootIsolationAlgorithm + -> Box + -> Writer ( [ Box ], [ Box ] ) + ( RootIsolationStep, [ Box ] ) +doStrategy roundHistory previousRoundsHistory eqs minWidth algo box = + case algo of + GaussSeidel ( GaussSeidelOptions { gsPreconditioner } ) -> do + boxes <- intervalGaussSeidel eqs gsPreconditioner box + return ( GaussSeidelStep, boxes ) + Box1 ( Box1Options { box1EpsEq } ) -> + return ( Box1Step, makeBox1Consistent eqs minWidth box1EpsEq box ) + Box2 ( Box2Options { box2LambdaMin, box2EpsEq } ) -> + return ( Box2Step, [ makeBox2Consistent eqs minWidth box2EpsEq box2LambdaMin box ] ) + Bisection ( BisectionOptions { canHaveSols, fallbackBisectionDim } ) -> do + let ( boxes, whatBis ) = bisect ( canHaveSols eqs ) ( fallbackBisectionDim roundHistory previousRoundsHistory eqs ) box + return ( BisectionStep whatBis, boxes ) + +-- | Bisect the given box. +-- +-- (The difficult part lies in determining along which dimension to bisect.) +bisect :: ( Box -> Bool ) + -- ^ how to check whether a box contains solutions + -> ( Int, String ) + -- ^ fall-back choice of dimension (and "why" we chose it) + -> Box + -> ( [ Box ], ( String, Double ) ) +bisect canHaveSols ( dim, why ) + ( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) + -- 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. + = let t_mid = 0.5 * ( t_lo + t_hi ) + s_mid = 0.5 * ( s_lo + s_hi ) + l = 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_mid s_hi ) + r = 𝕀 ( ℝ2 t_mid s_lo ) ( ℝ2 t_hi s_hi ) + d = 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_mid ) + u = 𝕀 ( ℝ2 t_lo s_mid ) ( ℝ2 t_hi s_hi ) + l_skip = not $ canHaveSols l + r_skip = not $ canHaveSols r + d_skip = not $ canHaveSols d + u_skip = not $ canHaveSols u + + ( bisGuesses, whatBis ) + | l_skip && r_skip + = ( [], ( "lr", t_mid ) ) + | d_skip && u_skip + = ( [], ( "du", s_mid ) ) + | l_skip + = ( [ r ], ( "r", t_mid ) ) + | r_skip + = ( [ l ], ( "l", t_mid ) ) + | d_skip + = ( [ u ], ( "u", s_mid ) ) + | u_skip + = ( [ d ], ( "d", s_mid ) ) + | otherwise + = let why' = case why of + "" -> "" + _ -> " (" ++ why ++ ")" + in case dim of + 0 -> ( [ l, r ], ( "t" ++ why', t_mid ) ) + 1 -> ( [ d, u ], ( "s" ++ why', s_mid ) ) + _ -> error "bisect: fall-back bisection dimension should be either 0 or 1" + in ( bisGuesses, whatBis ) + +-- | Interval Newton method with Gauss–Seidel step. +intervalGaussSeidel + :: ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) ) + -- ^ equations + -> Preconditioner + -- ^ preconditioner to use for the Gauss–Seidel step + -> 𝕀ℝ 2 + -> Writer ( [ 𝕀ℝ 2 ], [ 𝕀ℝ 2 ] ) [ 𝕀ℝ 2 ] +intervalGaussSeidel + eqs + precondMethod + ts@( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) ) + | D12 _f ( T ee_f_t ) ( T ee_f_s ) + <- eqs ts + , let f_t = take2 ee_f_t + f_s = take2 ee_f_s + , D12 ee_f_mid ( T _ee_f_t_mid ) ( T _ee_f_s_mid ) + <- eqs ts_mid + , let f_mid = take2 ee_f_mid + + = let -- Interval Newton method: take one Gauss–Seidel step + -- for the equation f'(x) v = - f(x_mid), + -- where f = 𝛿E/𝛿s * dc/dt + ( a, b ) = precondition precondMethod + ( midI f_t, midI f_s ) + ( f_t, f_s ) ( neg f_mid ) + + -- NB: we have to re-center around zero to take a Gauss–Seidel step. + gsGuesses = map ( first ( \ ts' -> unT $ T ts' ^+^ T ts_mid ) ) + $ gaussSeidelStep a b ( unT $ T ts ^-^ T ts_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 + in do tell ([], done) + return todo + where + t_mid = 0.5 * ( t_lo + t_hi ) + s_mid = 0.5 * ( s_lo + s_hi ) + ts_mid = singleton ( ℝ2 t_mid s_mid ) + neg ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) + = let !( 𝕀 x'_lo x'_hi ) = negate $ 𝕀 x_lo x_hi + !( 𝕀 y'_lo y'_hi ) = negate $ 𝕀 y_lo y_hi + in 𝕀 ( ℝ2 x'_lo y'_lo ) ( ℝ2 x'_hi y'_hi ) + take2 :: 𝕀ℝ 3 -> 𝕀ℝ 2 + take2 ( 𝕀 ( ℝ3 _a_lo b_lo c_lo ) ( ℝ3 _a_hi b_hi c_hi ) ) + = 𝕀 ( ℝ2 b_lo c_lo ) ( ℝ2 b_hi c_hi ) + -- TODO: always using the last 2 coordinates, instead of any + -- pair of two coordinates. + +-- Take one interval Gauss–Seidel step for the equation \( A X = B \), +-- refining the initial guess box for \( X \) into up to four (disjoint) new boxes. +-- +-- The boolean indicates whether the Gauss–Seidel step was a contraction. +gaussSeidelStep :: ( 𝕀ℝ 2, 𝕀ℝ 2 ) -- ^ columns of \( A \) + -> 𝕀ℝ 2 -- ^ \( B \) + -> 𝕀ℝ 2 -- ^ initial box \( X \) + -> [ ( 𝕀ℝ 2, Bool ) ] +gaussSeidelStep + ( 𝕀 ( ℝ2 a11_lo a21_lo ) ( ℝ2 a11_hi a21_hi ) + , 𝕀 ( ℝ2 a12_lo a22_lo ) ( ℝ2 a12_hi a22_hi ) ) + ( 𝕀 ( ℝ2 b1_lo b2_lo ) ( ℝ2 b1_hi b2_hi ) ) + ( 𝕀 ( ℝ2 x1_lo x2_lo ) ( ℝ2 x1_hi x2_hi ) ) + = let !a11 = 𝕀 a11_lo a11_hi + !a12 = 𝕀 a12_lo a12_hi + !a21 = 𝕀 a21_lo a21_hi + !a22 = 𝕀 a22_lo a22_hi + !b1 = 𝕀 b1_lo b1_hi + !b2 = 𝕀 b2_lo b2_hi + !x1 = 𝕀 x1_lo x1_hi + !x2 = 𝕀 x2_lo x2_hi + in nub $ do + + -- x1' = ( ( b1 - a12 * x2 ) / a11 ) ∩ x1 + x1'_0 <- ( b1 - a12 * x2 ) `extendedDivide` a11 + ( x1'@( 𝕀 x1'_lo x1'_hi ), sub_x1 ) <- x1'_0 `intersect` x1 + -- x2' = ( ( b2 - a21 * x1' ) / a22 ) ∩ x2 + x2'_0 <- ( b2 - a21 * x1' ) `extendedDivide` a22 + ( 𝕀 x2'_lo x2'_hi, sub_x2 ) <- x2'_0 `intersect` x2 + + return ( ( 𝕀 ( ℝ2 x1'_lo x2'_lo) ( ℝ2 x1'_hi x2'_hi ) ) + , sub_x1 && sub_x2 ) +-- 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" + +-- | 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 + :: ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) ) + -> Double -> Double + -> Box -> [ Box ] +makeBox1Consistent eqs minWidth epsEq x = + ( `State.evalState` False ) $ + pipeFunctionsWhileTrue ( allNarrowingOperators epsEq minWidth 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 + :: ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) ) + -> Double -> Double -> Double + -> Box -> Box +makeBox2Consistent eqs minWidth epsEq lambdaMin x0 = ( `State.evalState` False ) $ doLoop 0.25 x0 + where + doLoop :: Double -> Box -> State Bool Box + doLoop lambda x = do + x' <- boundConsistency get_t set_t lambda x + x'' <- boundConsistency get_s set_s lambda x' + modified <- State.get + let lambda' = if modified then lambda else 0.5 * lambda + if lambda' < lambdaMin + then return x'' + else do { State.put False ; doLoop lambda' x'' } + + boundConsistency :: ( Box -> 𝕀 Double ) + -> ( 𝕀 Double -> Box -> Box ) + -> Double -> Box -> State Bool Box + boundConsistency getter setter lambda box = do + let x@( 𝕀 x_inf x_sup ) = getter box + c1 = ( 1 - lambda ) * x_inf + lambda * x_sup + c2 = lambda * x_inf + ( 1 - lambda ) * x_sup + x'_inf = + case makeBox1Consistent eqs minWidth epsEq ( setter ( 𝕀 x_inf c1 ) box ) of + [] -> c1 + x's -> minimum $ map ( inf . getter ) x's + x'_sup = + case makeBox1Consistent eqs minWidth epsEq ( 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' >= epsEq ) $ + State.put True + return $ setter x' box + +precondition :: Preconditioner + -> ( 𝕀ℝ 2, 𝕀ℝ 2 ) + -> ( 𝕀ℝ 2, 𝕀ℝ 2 ) + -> 𝕀ℝ 2 + -> ( ( 𝕀ℝ 2, 𝕀ℝ 2 ), 𝕀ℝ 2 ) +precondition meth jac_mid a@( a1, a2 ) b = + case meth of + NoPreconditioning + -> ( a, b ) + InverseMidJacobian + | ( 𝕀 ( ℝ2 a11_lo a21_lo ) ( ℝ2 a11_hi a21_hi ) + , 𝕀 ( ℝ2 a12_lo a22_lo ) ( ℝ2 a12_hi a22_hi ) ) <- jac_mid + , let !a11 = 0.5 * ( a11_lo + a11_hi ) + !a12 = 0.5 * ( a12_lo + a12_hi ) + !a21 = 0.5 * ( a21_lo + a21_hi ) + !a22 = 0.5 * ( a22_lo + a22_hi ) + !d = a11 * a22 - a12 * a21 + , not ( nearZero d ) + , let !precond = ( ℝ2 a22 -a12, ℝ2 -a21 a11 ) + !inv = recip d + f x = scale inv $ matMulVec precond x + -> ( ( f a1, f a2 ), f b ) + | otherwise + -> ( a, b ) + +scale :: Double -> 𝕀ℝ 2 -> 𝕀ℝ 2 +scale s ( 𝕀 ( ℝ2 a1_lo a2_lo ) ( ℝ2 a1_hi a2_hi ) ) + = 𝕀 ( ℝ2 b1_lo b2_lo ) ( ℝ2 b1_hi b2_hi ) + where + 𝕀 b1_lo b1_hi = scaleInterval s ( 𝕀 a1_lo a1_hi ) + 𝕀 b2_lo b2_hi = scaleInterval s ( 𝕀 a2_lo a2_hi ) + +matMulVec :: ( ℝ 2, ℝ 2 ) -> 𝕀ℝ 2 -> 𝕀ℝ 2 +matMulVec ( ℝ2 a11 a21, ℝ2 a12 a22 ) ( 𝕀 ( ℝ2 u_lo v_lo ) ( ℝ2 u_hi v_hi ) ) + = 𝕀 ( ℝ2 u'_lo v'_lo ) ( ℝ2 u'_hi v'_hi ) + where + u = 𝕀 u_lo u_hi + v = 𝕀 v_lo v_hi + 𝕀 u'_lo u'_hi = scaleInterval a11 u + scaleInterval a12 v + 𝕀 v'_lo v'_hi = scaleInterval a21 u + scaleInterval a22 v + +midI :: 𝕀ℝ 2 -> 𝕀ℝ 2 +midI ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = + singleton $ ℝ2 ( 0.5 * ( x_lo + x_hi ) ) ( 0.5 * ( y_lo + y_hi ) ) +--width :: 𝕀ℝ 1 -> Double +--width ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) = hi - lo +--normI :: 𝕀ℝ 1 -> Double +--normI ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) = sqrt $ sup $ ( 𝕀 lo hi ) Ring.^ 2 +--normVI :: 𝕀ℝ 2 -> Double +--normVI ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = +-- sqrt $ sup $ ( 𝕀 x_lo x_hi ) Ring.^ 2 + ( 𝕀 y_lo y_hi ) Ring.^ 2 +normVI3 :: 𝕀ℝ 3 -> Double +normVI3 ( 𝕀 ( ℝ3 lo x_lo y_lo ) ( ℝ3 hi x_hi y_hi ) ) = + sqrt $ max ( abs lo ) ( abs hi ) Ring.^ 2 + + max ( abs x_lo ) ( abs x_hi ) Ring.^ 2 + + max ( abs y_lo ) ( abs y_hi ) Ring.^ 2 + +-- 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 eps_equal eps_bisection 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 delta <- f_x_left `extendedDivide` f'_x --f'_x_left + let x_new = x_left - delta + map fst $ x_new `intersect` x + in + if | null x's + -> [] + | ( width x - sum ( map width x's ) ) < eps_equal + -- TODO: do a check with the relative reduction in size? + -> x's + | otherwise + -> do + x' <- x's + if sup x' - inf x' < eps_bisection + then return x' + else left_narrow =<< bisectI x' + +-- TODO: de-duplicate with 'leftNarrow'? +rightNarrow :: Double + -> Double + -> ( 𝕀 Double -> ( 𝕀 Double, 𝕀 Double ) ) + -> 𝕀 Double + -> [ 𝕀 Double ] +rightNarrow eps_equal eps_bisection 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 delta <- f_x_right `extendedDivide` f'_x --f'_x_right + let x_new = x_right - delta + map fst $ x_new `intersect` x + in + if | null x's + -> [] + | ( width x - sum ( map width x's ) ) < eps_equal + -- TODO: do a check with the relative reduction in size? + -> x's + | otherwise + -> do + x' <- x's + if sup x' - inf x' < eps_bisection + 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 + :: Double + -> Double + -> ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) ) + -> [ Box -> State Bool [ Box ] ] +allNarrowingOperators eps_eq eps_bis eqs = + [ \ cand -> + let getter = ( `index` coordIndex ) + setter = set coordIndex + newCands = map ( `setter` cand ) $ narrowFn eps_eq eps_bis ( \ x0 -> ff' $ setter x0 cand ) ( getter cand ) + in do + when ( ( width ( getter cand ) - sum ( map ( width . getter ) newCands ) ) > eps_eq ) $ + -- TODO: do a check with the relative reduction in size? + State.put True + return newCands + | narrowFn <- [ leftNarrow, rightNarrow ] + , ( coordIndex, ff' ) <- + [ ( Fin 1, ff'_t i ) | i <- [ Fin 1, Fin 2, Fin 3 ] ] + ++ + [ ( Fin 2, ff'_s i ) | i <- [ Fin 1, Fin 2, Fin 3 ] ] + ] + where + + ff'_t i ts = + let D12 { _D12_v = f + , _D12_dx = T f_t + } = eqs ts + in ( f `index` i, f_t `index` i ) + ff'_s i ts = + let D12 { _D12_v = f + , _D12_dy = T f_t + } = eqs ts + in ( f `index` i, f_t `index` i ) + +get_t, get_s :: Box -> 𝕀 Double +get_t = ( `index` ( Fin 1 ) ) +get_s = ( `index` ( Fin 2 ) ) + +set_t, set_s :: 𝕀 Double -> Box -> Box +set_t = set ( Fin 1 ) +set_s = set ( Fin 2 )