From c89fba7fd2700e7d4c9f1af221ff62cd499bf5cd Mon Sep 17 00:00:00 2001 From: sheaf Date: Thu, 25 Apr 2024 21:53:27 +0200 Subject: [PATCH] Improve specialisation + add degree computation --- brush-strokes/brush-strokes.cabal | 10 ++ brush-strokes/src/cusps/bench/Bench/Types.hs | 48 ++++++- brush-strokes/src/cusps/bench/Main.hs | 5 +- brush-strokes/src/lib/Calligraphy/Brushes.hs | 28 ++-- brush-strokes/src/lib/Math/Bezier/Stroke.hs | 65 ++++++++- .../Math/Bezier/Stroke/EnvelopeEquation.hs | 15 +++ brush-strokes/src/lib/Math/Interval.hs | 29 ++-- brush-strokes/src/lib/Math/Linear.hs | 2 +- brush-strokes/src/lib/Math/Monomial.hs | 11 +- brush-strokes/src/lib/Math/Root/Isolation.hs | 22 +-- .../src/lib/Math/Root/Isolation/Bisection.hs | 18 ++- .../src/lib/Math/Root/Isolation/Core.hs | 39 +++--- .../src/lib/Math/Root/Isolation/Degree.hs | 126 ++++++++++++++++++ .../lib/Math/Root/Isolation/GaussSeidel.hs | 75 +++++++---- .../src/lib/Math/Root/Isolation/Narrowing.hs | 27 +++- src/app/MetaBrush/Application.hs | 18 ++- src/app/MetaBrush/Render/Document.hs | 10 +- 17 files changed, 432 insertions(+), 116 deletions(-) create mode 100644 brush-strokes/src/lib/Math/Root/Isolation/Degree.hs diff --git a/brush-strokes/brush-strokes.cabal b/brush-strokes/brush-strokes.cabal index 23cf924..e63d12b 100644 --- a/brush-strokes/brush-strokes.cabal +++ b/brush-strokes/brush-strokes.cabal @@ -14,6 +14,11 @@ flag use-fma default: True manual: False +flag asserts + description: Enable debugging assertions. + default: False + manual: True + common common build-depends: @@ -89,6 +94,10 @@ common common base >= 4.19 + if flag(asserts) + cpp-options: + -DASSERTS + common extra build-depends: @@ -133,6 +142,7 @@ library , Math.Root.Isolation , Math.Root.Isolation.Bisection , Math.Root.Isolation.Core + , Math.Root.Isolation.Degree , Math.Root.Isolation.GaussSeidel , Math.Root.Isolation.Narrowing , Debug.Utils diff --git a/brush-strokes/src/cusps/bench/Bench/Types.hs b/brush-strokes/src/cusps/bench/Bench/Types.hs index 5c559d7..bc67c07 100644 --- a/brush-strokes/src/cusps/bench/Bench/Types.hs +++ b/brush-strokes/src/cusps/bench/Bench/Types.hs @@ -12,8 +12,14 @@ module Bench.Types -- base import Data.Coerce ( coerce ) +import Data.Proxy + ( Proxy(..) ) +import Data.Type.Equality + ( (:~:)(Refl) ) import GHC.Generics ( Generic ) +import GHC.TypeNats + ( sameNat ) -- containers import Data.Sequence @@ -111,10 +117,48 @@ getStrokeFunctions ( Brush brushShape brushShapeI mbRot ) sp0 crv = ( fmap nonDecreasing mbRot ) ) {-# INLINEABLE getStrokeFunctions #-} +{-# SPECIALISE getStrokeFunctions + :: Brush ( ℝ 1 ) + -> Point 1 + -> Curve Open () ( Point 1 ) + -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) + , 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) + #-} +{-# SPECIALISE getStrokeFunctions + :: Brush ( ℝ 2 ) + -> Point 2 + -> Curve Open () ( Point 2 ) + -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) + , 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) + #-} +{-# SPECIALISE getStrokeFunctions + :: Brush ( ℝ 3 ) + -> Point 3 + -> Curve Open () ( Point 3 ) + -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) + , 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) + #-} +{-# SPECIALISE getStrokeFunctions + :: Brush ( ℝ 4 ) + -> Point 4 + -> Curve Open () ( Point 4 ) + -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) + , 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) + #-} brushStrokeFunctions :: BrushStroke -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) , 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) -brushStrokeFunctions ( BrushStroke { stroke = ( sp0, crv ), brush } ) = - getStrokeFunctions brush sp0 crv +brushStrokeFunctions ( BrushStroke { stroke = ( sp0, crv ), brush = brush :: Brush ( ℝ nbParams ) } ) + -- Trick for triggering specialisation... TODO improve this. + | Just Refl <- sameNat @nbParams @1 Proxy Proxy + = getStrokeFunctions @1 brush sp0 crv + | Just Refl <- sameNat @nbParams @2 Proxy Proxy + = getStrokeFunctions @2 brush sp0 crv + | Just Refl <- sameNat @nbParams @3 Proxy Proxy + = getStrokeFunctions @3 brush sp0 crv + | Just Refl <- sameNat @nbParams @4 Proxy Proxy + = getStrokeFunctions @4 brush sp0 crv + | otherwise + = getStrokeFunctions @nbParams brush sp0 crv diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs index 06a3a34..b723f46 100644 --- a/brush-strokes/src/cusps/bench/Main.hs +++ b/brush-strokes/src/cusps/bench/Main.hs @@ -92,15 +92,14 @@ benchGroups :: [ ( String, NE.NonEmpty TestCase ) ] benchGroups = [ ( "ellipse" , NE.fromList - [ ellipseTestCase opts ("Ξ΅_bis=" ++ show Ξ΅_bis ++ if doBox1 then "box(1)" else "") + [ ellipseTestCase opts ("Ξ΅_bis=" ++ show Ξ΅_bis) ( 0, 1 ) pi ( defaultStartBoxes [ 0 .. 3 ] ) | Ξ΅_bis <- [ 1e-6, 5e-6, 1e-5, 5e-5, 1e-4, 5e-4, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 0.1, 0.2, 0.3 ] - , doBox1 <- [ False, True ] , let opts = RootIsolationOptions { rootIsolationAlgorithms = - defaultRootIsolationAlgorithms minWidth Ξ΅_bis doBox1 + defaultRootIsolationAlgorithms minWidth Ξ΅_bis } ] ) diff --git a/brush-strokes/src/lib/Calligraphy/Brushes.hs b/brush-strokes/src/lib/Calligraphy/Brushes.hs index a1d0b53..402eaae 100644 --- a/brush-strokes/src/lib/Calligraphy/Brushes.hs +++ b/brush-strokes/src/lib/Calligraphy/Brushes.hs @@ -50,9 +50,13 @@ type BrushFn i k brushParams -- | A brush, described as a base shape + an optional rotation. data Brush brushParams = Brush - { brushShape :: BrushFn () 2 brushParams - , brushShapeI :: BrushFn 𝕀 3 brushParams - , mbRotation :: Maybe ( brushParams -> Double ) + { -- | Base brush shape, before applying any rotation (if any). + brushBaseShape :: BrushFn () 2 brushParams + -- | Base brush shape, before applying any rotation (if any). + , brushBaseShapeI :: BrushFn 𝕀 3 brushParams + -- | Optional rotation angle function + -- (assumed to be a linear function). + , mbRotation :: Maybe ( brushParams -> Double ) } -------------------------------------------------------------------------------- @@ -74,27 +78,27 @@ type ParamsICt k i rec = circleBrush :: ( 1 <= RepDim params, ParamsCt params ) => Brush params circleBrush = Brush - { brushShape = circleBrushFn @() @2 proxy# id - , brushShapeI = circleBrushFn @𝕀 @3 proxy# singleton - , mbRotation = Nothing + { brushBaseShape = circleBrushFn @() @2 proxy# id + , brushBaseShapeI = circleBrushFn @𝕀 @3 proxy# singleton + , mbRotation = Nothing } {-# INLINEABLE ellipseBrush #-} ellipseBrush :: ( 3 <= RepDim params, ParamsCt params ) => Brush params ellipseBrush = Brush - { brushShape = ellipseBrushFn @() @2 proxy# id - , brushShapeI = ellipseBrushFn @𝕀 @3 proxy# singleton - , mbRotation = Just ( `index` ( Fin 3 ) ) + { brushBaseShape = ellipseBrushFn @() @2 proxy# id + , brushBaseShapeI = ellipseBrushFn @𝕀 @3 proxy# singleton + , mbRotation = Just ( `index` ( Fin 3 ) ) } {-# INLINEABLE tearDropBrush #-} tearDropBrush :: ( 3 <= RepDim params, ParamsCt params ) => Brush params tearDropBrush = Brush - { brushShape = tearDropBrushFn @() @2 proxy# id - , brushShapeI = tearDropBrushFn @𝕀 @3 proxy# singleton - , mbRotation = Just ( `index` ( Fin 3 ) ) + { brushBaseShape = tearDropBrushFn @() @2 proxy# id + , brushBaseShapeI = tearDropBrushFn @𝕀 @3 proxy# singleton + , mbRotation = Just ( `index` ( Fin 3 ) ) } -------------------------------------------------------------------------------- diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke.hs b/brush-strokes/src/lib/Math/Bezier/Stroke.hs index dc9b90a..dbccb1a 100644 --- a/brush-strokes/src/lib/Math/Bezier/Stroke.hs +++ b/brush-strokes/src/lib/Math/Bezier/Stroke.hs @@ -383,7 +383,7 @@ computeStrokeOutline rootAlgo mbCuspOptions fitParams ptParams toBrushParams bru brushShape :: ptData -> SplinePts Closed brushShape pt = - let Brush { brushShape = shapeFn, mbRotation = mbRot } = brush + let Brush { brushBaseShape = shapeFn, mbRotation = mbRot } = brush brushParams = toBrushParams $ ptParams pt shape = fun @Double shapeFn brushParams in case mbRot of @@ -546,7 +546,7 @@ outlineFunction -> Curve Open crvData ptData -> OutlineInfo outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams - ( Brush { brushShape, brushShapeI, mbRotation } ) = \ sp0 crv -> + ( Brush { brushBaseShape, brushBaseShapeI, mbRotation } ) = \ sp0 crv -> let usedParams :: C 2 ( ℝ 1 ) usedParams @@ -561,7 +561,7 @@ outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams coerce coerce path ( fmap toBrushParams usedParams ) - brushShape + brushBaseShape mbRotation curvesI :: 𝕀ℝ 1 -- t @@ -571,7 +571,7 @@ outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams coerce coerce pathI ( fmap ( nonDecreasing toBrushParams ) usedParamsI ) - brushShapeI + brushBaseShapeI ( fmap nonDecreasing mbRotation ) usedParamsI :: C 3 ( 𝕀ℝ 1 ) ( 𝕀 usedParams ) @@ -589,9 +589,19 @@ outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams D21 path_t path'_t _ = runD path t D21 params_t _ _ = runD usedParams t - brush_t = value @Double @2 @brushParams - $ runD brushShape - $ toBrushParams params_t + brush_t = + let brushParams = toBrushParams params_t + applyRot = case mbRotation of + Nothing -> id + Just getΞΈ -> + let ΞΈ = getΞΈ brushParams + cosΞΈ = cos ΞΈ + sinΞΈ = sin ΞΈ + in fmap ( unT . rotate cosΞΈ sinΞΈ . T ) + in + applyRot $ + value @Double @2 @brushParams $ + runD brushBaseShape brushParams ( potentialCusps, definiteCusps ) = case mbCuspOptions of @@ -644,6 +654,7 @@ pathAndUsedParams co toI ptParams sp0 crv = | let bez3 = Cubic.Bezier sp0 sp1 sp2 sp3 -> ( bezier3 @k @i @( ℝ 2 ) co ( fmap ( toI . coords ) bez3 ) , bezier3 @k @i @usedParams co ( fmap ( toI . ptParams ) bez3 ) ) +{-# INLINEABLE pathAndUsedParams #-} ----------------------------------- -- Various utility functions @@ -1030,6 +1041,46 @@ brushStrokeData co1 co2 path params brush mbBrushRotation = let dbrush_t_s = dbrush_t s mbRotation = mbBrushRotation <&> \ getTheta -> fmap getTheta dparams_t in envelopeEquation @k @i Proxy co1 dpath_t dbrush_t_s mbRotation +{-# INLINEABLE brushStrokeData #-} +{-# SPECIALISE brushStrokeData + :: ( I 𝕀 Double -> I 𝕀 ( ℝ 1 ) ) + -> ( I 𝕀 ( ℝ 1 ) -> I 𝕀 Double ) + -> ( C 3 ( I 𝕀 ( ℝ 1 ) ) ( I 𝕀 ( ℝ 2 ) ) ) + -> ( C 3 ( I 𝕀 ( ℝ 1 ) ) ( I 𝕀 ( ℝ 1 ) ) ) + -> ( C 3 ( I 𝕀 ( ℝ 1 ) ) ( Spline Closed () ( I 𝕀 ( ℝ 2 ) ) ) ) + -> ( Maybe ( I 𝕀 ( ℝ 1 ) -> I 𝕀 Double ) ) + -> ( I 𝕀 ( ℝ 1 ) -> Seq ( I 𝕀 ( ℝ 1 ) -> StrokeDatum 3 𝕀 ) ) + #-} +{-# SPECIALISE brushStrokeData + :: ( I 𝕀 Double -> I 𝕀 ( ℝ 1 ) ) + -> ( I 𝕀 ( ℝ 1 ) -> I 𝕀 Double ) + -> ( C 3 ( I 𝕀 ( ℝ 1 ) ) ( I 𝕀 ( ℝ 2 ) ) ) + -> ( C 3 ( I 𝕀 ( ℝ 1 ) ) ( I 𝕀 ( ℝ 2 ) ) ) + -> ( C 3 ( I 𝕀 ( ℝ 2 ) ) ( Spline Closed () ( I 𝕀 ( ℝ 2 ) ) ) ) + -> ( Maybe ( I 𝕀 ( ℝ 2 ) -> I 𝕀 Double ) ) + -> ( I 𝕀 ( ℝ 1 ) -> Seq ( I 𝕀 ( ℝ 1 ) -> StrokeDatum 3 𝕀 ) ) + #-} +{-# SPECIALISE brushStrokeData + :: ( I 𝕀 Double -> I 𝕀 ( ℝ 1 ) ) + -> ( I 𝕀 ( ℝ 1 ) -> I 𝕀 Double ) + -> ( C 3 ( I 𝕀 ( ℝ 1 ) ) ( I 𝕀 ( ℝ 2 ) ) ) + -> ( C 3 ( I 𝕀 ( ℝ 1 ) ) ( I 𝕀 ( ℝ 3 ) ) ) + -> ( C 3 ( I 𝕀 ( ℝ 3 ) ) ( Spline Closed () ( I 𝕀 ( ℝ 2 ) ) ) ) + -> ( Maybe ( I 𝕀 ( ℝ 3 ) -> I 𝕀 Double ) ) + -> ( I 𝕀 ( ℝ 1 ) -> Seq ( I 𝕀 ( ℝ 1 ) -> StrokeDatum 3 𝕀 ) ) + #-} +{-# SPECIALISE brushStrokeData + :: ( I 𝕀 Double -> I 𝕀 ( ℝ 1 ) ) + -> ( I 𝕀 ( ℝ 1 ) -> I 𝕀 Double ) + -> ( C 3 ( I 𝕀 ( ℝ 1 ) ) ( I 𝕀 ( ℝ 2 ) ) ) + -> ( C 3 ( I 𝕀 ( ℝ 1 ) ) ( I 𝕀 ( ℝ 4 ) ) ) + -> ( C 3 ( I 𝕀 ( ℝ 4 ) ) ( Spline Closed () ( I 𝕀 ( ℝ 2 ) ) ) ) + -> ( Maybe ( I 𝕀 ( ℝ 4 ) -> I 𝕀 Double ) ) + -> ( I 𝕀 ( ℝ 1 ) -> Seq ( I 𝕀 ( ℝ 1 ) -> StrokeDatum 3 𝕀 ) ) + #-} +-- TODO: these specialisations fire in the benchmarking code because +-- we instantiate brushParams with ( ℝ nbBrushParams ), but they won't fire +-- in the main app code because we are using types such as "Record EllipseBrushFields". -------------------------------------------------------------------------------- -- Solving the envelolpe equation: root-finding. diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs b/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs index 933f828..60f2710 100644 --- a/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs +++ b/brush-strokes/src/lib/Math/Bezier/Stroke/EnvelopeEquation.hs @@ -151,22 +151,26 @@ instance HasBΓ©zier 2 () where D21 ( lerp @( T b ) t a b ) ( a --> b ) origin + {-# INLINEABLE line #-} bezier2 co ( bez :: Quadratic.Bezier b ) = D \ ( co -> t ) -> D21 ( Quadratic.bezier @( T b ) bez t ) ( Quadratic.bezier' bez t ) ( Quadratic.bezier'' bez ) + {-# INLINEABLE bezier2 #-} bezier3 co ( bez :: Cubic.Bezier b ) = D \ ( co -> t ) -> D21 ( Cubic.bezier @( T b ) bez t ) ( Cubic.bezier' bez t ) ( Cubic.bezier'' bez t ) + {-# INLINEABLE bezier3 #-} instance HasEnvelopeEquation 2 where uncurryD = uncurryD2 + {-# INLINEABLE uncurryD #-} envelopeEquation ( _ :: Proxy i ) co dp@( D21 ( T -> p ) p_t p_tt ) @@ -247,6 +251,7 @@ instance HasEnvelopeEquation 2 where , D12 ( unT u ) u_t u_s , D12 ( unT v ) v_t v_s ) + {-# INLINEABLE envelopeEquation #-} instance HasBΓ©zier 3 () where @@ -256,6 +261,7 @@ instance HasBΓ©zier 3 () where ( a --> b ) origin origin + {-# INLINEABLE line #-} bezier2 co ( bez :: Quadratic.Bezier b ) = D \ ( co -> t ) -> @@ -263,6 +269,7 @@ instance HasBΓ©zier 3 () where ( Quadratic.bezier' bez t ) ( Quadratic.bezier'' bez ) origin + {-# INLINEABLE bezier2 #-} bezier3 co ( bez :: Cubic.Bezier b ) = D \ ( co -> t ) -> @@ -270,10 +277,12 @@ instance HasBΓ©zier 3 () where ( Cubic.bezier' bez t ) ( Cubic.bezier'' bez t ) ( Cubic.bezier''' bez ) + {-# INLINEABLE bezier3 #-} instance HasEnvelopeEquation 3 where uncurryD = uncurryD3 + {-# INLINEABLE uncurryD #-} envelopeEquation ( _ :: Proxy i ) co dp@( D31 ( T -> p ) p_t p_tt p_ttt ) @@ -392,6 +401,7 @@ instance HasEnvelopeEquation 3 where , D22 ( unT u ) u_t u_s u_tt u_ts u_ss , D22 ( unT v ) v_t v_s v_tt v_ts v_ss ) + {-# INLINEABLE envelopeEquation #-} instance HasBΓ©zier 3 𝕀 where @@ -401,6 +411,7 @@ instance HasBΓ©zier 3 𝕀 where ( a --> b ) origin origin + {-# INLINEABLE line #-} bezier2 co ( bez :: Quadratic.Bezier b ) = D \ ( co -> t ) -> @@ -408,6 +419,7 @@ instance HasBΓ©zier 3 𝕀 where ( Quadratic.bezier' bez t ) ( Quadratic.bezier'' bez ) origin + {-# INLINEABLE bezier2 #-} bezier3 co ( bez :: Cubic.Bezier b ) = D \ ( co -> t ) -> @@ -415,6 +427,7 @@ instance HasBΓ©zier 3 𝕀 where ( T $ aabb ( fmap unT $ Cubic.derivative ( fmap T bez ) ) ( `evaluateQuadratic` t ) ) ( Cubic.bezier'' bez t ) ( Cubic.bezier''' bez ) + {-# INLINEABLE bezier3 #-} {- Note [Computing BΓ©ziers over intervals] ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -455,6 +468,7 @@ evaluateCubic bez t = maxs = fmap (Cubic.bezier @( T Double ) sup_bez) $ inf t :| ( sup t : filter ( ∈ t ) ( Cubic.extrema sup_bez ) ) in 𝕀 ( minimum mins ) ( maximum maxs ) +{-# INLINEABLE evaluateCubic #-} -- | Evaluate a quadratic BΓ©zier curve, when both the coefficients and the -- parameter are intervals. @@ -468,3 +482,4 @@ evaluateQuadratic bez t = maxs = fmap (Quadratic.bezier @( T Double ) sup_bez) $ inf t :| ( sup t : filter ( ∈ t ) ( Quadratic.extrema sup_bez ) ) in 𝕀 ( minimum mins ) ( maximum maxs ) +{-# INLINEABLE evaluateQuadratic #-} diff --git a/brush-strokes/src/lib/Math/Interval.hs b/brush-strokes/src/lib/Math/Interval.hs index 0439ed9..a97408d 100644 --- a/brush-strokes/src/lib/Math/Interval.hs +++ b/brush-strokes/src/lib/Math/Interval.hs @@ -1,3 +1,5 @@ +{-# LANGUAGE CPP #-} + {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE UndecidableInstances #-} @@ -19,8 +21,6 @@ module Math.Interval -- base import Prelude hiding ( Num(..), Fractional(..) ) -import Data.List - ( nub ) import qualified Data.List.NonEmpty as NE -- acts @@ -138,11 +138,23 @@ instance Torsor ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where ------------------------------------------------------------------------------- -- Extended division --- | Extended division. +-- | Extended division, returning either 1 or 2 intervals. +-- +-- NB: this function can return overlapping intervals if both arguments contain 0. +-- Otherwise, the returned intervals are guaranteed to be disjoint. (⊘) :: 𝕀 Double -> 𝕀 Double -> [ 𝕀 Double ] -x ⊘ y = nub $ map ( x * ) ( extendedRecip y ) +x ⊘ y +#ifdef ASSERTS + | 0 ∈ x && 0 ∈ y + = error $ unlines + [ "x ⊘ y: both arguments contain zero" + , "x: " ++ show x, "y: " ++ show y ] + | otherwise +#endif + = map ( x * ) ( extendedRecip y ) infixl 7 ⊘ +-- | Extended reciprocal, returning either 1 or 2 intervals. extendedRecip :: 𝕀 Double -> [ 𝕀 Double ] extendedRecip x@( 𝕀 lo hi ) | lo == 0 && hi == 0 @@ -150,11 +162,12 @@ extendedRecip x@( 𝕀 lo hi ) | lo >= 0 || hi <= 0 = [ recip x ] | otherwise - = [ recip $ 𝕀 lo -0, recip $ 𝕀 0 hi ] + = [ 𝕀 negInf ( recip lo ), 𝕀 ( recip hi ) posInf ] where - negInf, posInf :: Double - negInf = -1 / 0 - posInf = 1 / 0 + +negInf, posInf :: Double +negInf = -1 / 0 +posInf = 1 / 0 ------------------------------------------------------------------------------- -- Lattices. diff --git a/brush-strokes/src/lib/Math/Linear.hs b/brush-strokes/src/lib/Math/Linear.hs index 7948193..bd4aa1a 100644 --- a/brush-strokes/src/lib/Math/Linear.hs +++ b/brush-strokes/src/lib/Math/Linear.hs @@ -133,7 +133,7 @@ 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 Functor ( Vec n ) deriving newtype instance Foldable ( Vec n ) deriving via ZipList instance Applicative ( Vec n ) diff --git a/brush-strokes/src/lib/Math/Monomial.hs b/brush-strokes/src/lib/Math/Monomial.hs index 2356202..f027a39 100644 --- a/brush-strokes/src/lib/Math/Monomial.hs +++ b/brush-strokes/src/lib/Math/Monomial.hs @@ -55,25 +55,24 @@ type family Deg f type family Vars f zeroMonomial :: forall k n. KnownNat n => Mon k n -zeroMonomial = Mon $ unsafeCoerce @[ Word ] @( Vec n Word ) - $ replicate ( fromIntegral $ word @n ) 0 -{-# INLINEABLE zeroMonomial #-} +zeroMonomial = Mon $ Vec $ replicate ( fromIntegral $ word @n ) 0 +{-# INLINE zeroMonomial #-} isZeroMonomial :: Mon k n -> Bool isZeroMonomial ( Mon ( Vec pows ) ) = all ( == 0 ) pows +{-# INLINE isZeroMonomial #-} totalDegree :: Mon k n -> Word totalDegree ( Mon ( Vec pows ) ) = sum pows linearMonomial :: forall k n. ( KnownNat n, 1 <= k ) => Fin n -> Mon k n linearMonomial ( Fin i' ) = - Mon $ unsafeCoerce @[ Word ] @( Vec n Word ) - $ replicate ( i - 1 ) 0 ++ [ 1 ] ++ replicate ( n - i ) 0 + Mon $ Vec $ replicate ( i - 1 ) 0 ++ [ 1 ] ++ replicate ( n - i ) 0 where n, i :: Int n = fromIntegral $ word @n i = fromIntegral i' -{-# INLINEABLE linearMonomial #-} +{-# INLINE linearMonomial #-} isLinear :: Mon k n -> Maybe ( Fin n ) isLinear ( Mon ( Vec pows ) ) = fmap Fin $ go 1 pows diff --git a/brush-strokes/src/lib/Math/Root/Isolation.hs b/brush-strokes/src/lib/Math/Root/Isolation.hs index e568a39..55f8e07 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation.hs @@ -102,7 +102,7 @@ newtype RootIsolationOptions n d defaultRootIsolationOptions :: BoxCt n d => RootIsolationOptions n d defaultRootIsolationOptions = RootIsolationOptions - { rootIsolationAlgorithms = defaultRootIsolationAlgorithms minWidth Ξ΅_eq True + { rootIsolationAlgorithms = defaultRootIsolationAlgorithms minWidth Ξ΅_eq } where minWidth = 1e-5 @@ -114,11 +114,10 @@ defaultRootIsolationAlgorithms . BoxCt n d => Double -- ^ minimum width of boxes (don't bisect further) -> Double -- ^ threshold for progress - -> Bool -- ^ do box1 -> BoxHistory n -> Box n -> Either String ( NE.NonEmpty ( RootIsolationAlgorithmWithOptions n d ) ) -defaultRootIsolationAlgorithms minWidth Ξ΅_eq doBox1 history box = +defaultRootIsolationAlgorithms minWidth Ξ΅_eq history box = case history of lastRoundBoxes : _ -- If, in the last round of strategies, we didn't try bisection... @@ -132,17 +131,19 @@ defaultRootIsolationAlgorithms minWidth Ξ΅_eq doBox1 history box = then Left $ "widths <= " ++ show minWidth else Right $ NE.singleton $ AlgoWithOptions @Bisection _bisOptions -- Otherwise, do a normal round. - -- Currently: we try an interval Gauss–Seidel step followed by box(1)-consistency. + -- Currently: we try an interval Gauss–Seidel. + -- (box(1)- and box(2)-consistency don't seem to help when using + -- the complete interval union Gauss–Seidel step) _ -> Right $ AlgoWithOptions @GaussSeidel _gsOptions - NE.:| [ AlgoWithOptions @Box1 _box1Options - | doBox1 && not verySmall ] + NE.:| [] where verySmall = and $ ( \ cd -> width cd <= minWidth ) <$> coordinates box _bisOptions = defaultBisectionOptions @n @d minWidth Ξ΅_eq box _gsOptions = defaultGaussSeidelOptions @n @d history - _box1Options = defaultBox1Options @n @d minWidth Ξ΅_eq + _box1Options = defaultBox1Options @n @d ( minWidth * 100 ) Ξ΅_eq + _box2Options = ( defaultBox2Options @n @d minWidth Ξ΅_eq ) { box2LambdaMin = 0.001 } -- Did we reduce the box width by at least Ξ΅_eq -- in at least one of the coordinates? @@ -249,8 +250,8 @@ isolateRootsIn ( RootIsolationOptions { rootIsolationAlgorithms } ) -- | Execute a root isolation strategy, replacing the input box with -- (hopefully smaller) output boxes. doStrategy - :: BoxCt n d - => RootIsolationAlgorithmWithOptions n d + :: forall n d + . RootIsolationAlgorithmWithOptions n d -> [ ( RootIsolationStep, Box n ) ] -> BoxHistory n -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) @@ -260,6 +261,5 @@ doStrategy doStrategy algo thisRoundHist prevRoundsHist eqs cand = case algo of AlgoWithOptions @algo options -> do - ( step, res ) <- rootIsolationAlgorithm options thisRoundHist prevRoundsHist eqs cand + ( step, res ) <- rootIsolationAlgorithm @algo options thisRoundHist prevRoundsHist eqs cand return ( SomeRootIsolationStep @algo step, res ) -{-# INLINEABLE doStrategy #-} diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Bisection.hs b/brush-strokes/src/lib/Math/Root/Isolation/Bisection.hs index 52f95d0..fc2db0a 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation/Bisection.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation/Bisection.hs @@ -1,5 +1,6 @@ -{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE UndecidableInstances #-} module Math.Root.Isolation.Bisection ( -- * The bisection method for root isolation @@ -35,6 +36,10 @@ import GHC.TypeNats , type (<=) ) +-- transformers +import Control.Monad.Trans.Writer.CPS + ( Writer ) + -- MetaBrush import Math.Algebra.Dual ( D ) @@ -52,7 +57,7 @@ import Math.Root.Isolation.Core -- | The bisection algorithm; see 'bisection'. data Bisection -instance RootIsolationAlgorithm Bisection where +instance BoxCt n d => RootIsolationAlgorithm Bisection n d where type instance StepDescription Bisection = ( String, Double ) type instance RootIsolationAlgorithmOptions Bisection n d = BisectionOptions n d rootIsolationAlgorithm @@ -65,6 +70,15 @@ instance RootIsolationAlgorithm Bisection where box return ( whatBis, boxes ) {-# INLINEABLE rootIsolationAlgorithm #-} + {-# SPECIALISE rootIsolationAlgorithm + :: RootIsolationAlgorithmOptions Bisection 2 3 + -> [ ( RootIsolationStep, Box 2 ) ] + -> BoxHistory 2 + -> ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) ) + -> Box 2 + -> Writer ( DoneBoxes 2 ) ( StepDescription Bisection, [ Box 2 ] ) #-} + -- NB: including this to be safe. The specialiser seems to sometimes + -- be able to generate this specialisation on its own, and sometimes not. -- | Options for the bisection method. type BisectionOptions :: Nat -> Nat -> Type diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Core.hs b/brush-strokes/src/lib/Math/Root/Isolation/Core.hs index f3e2294..e6bcb29 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation/Core.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation/Core.hs @@ -1,8 +1,9 @@ -{-# LANGUAGE AllowAmbiguousTypes #-} -{-# LANGUAGE PolyKinds #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE UndecidableSuperClasses #-} -- | Core definitions and utilities common to root isolation methods. module Math.Root.Isolation.Core @@ -122,9 +123,11 @@ noDoneBoxes = DoneBoxes [] [] -- | Existential wrapper over any root isolation algorithm, -- with the options necessary to run it. +type RootIsolationAlgorithmWithOptions :: Nat -> Nat -> Type data RootIsolationAlgorithmWithOptions n d where AlgoWithOptions - :: RootIsolationAlgorithm ty + :: forall {k :: Type} {n :: Nat} {d :: Nat} ( ty :: k ) + . RootIsolationAlgorithm ty n d => RootIsolationAlgorithmOptions ty n d -> RootIsolationAlgorithmWithOptions n d @@ -133,15 +136,15 @@ data RootIsolationAlgorithmWithOptions n d where -- This design keeps the set of root isolation algorithms open-ended, -- while retaining the ability to inspect previous steps (using the -- 'IsolationStep' pattern). -type RootIsolationAlgorithm :: forall {k}. k -> Constraint -class ( Typeable ty, Show ( StepDescription ty ) ) - => RootIsolationAlgorithm ty where +type RootIsolationAlgorithm :: forall {k :: Type}. k -> Nat -> Nat -> Constraint +class ( Typeable ty, Show ( StepDescription ty ), BoxCt n d ) + => RootIsolationAlgorithm ty n d where -- | The type of additional information about an algorithm step. -- -- Only really useful for debugging; gets stored in 'RootIsolationTree's. type StepDescription ty -- | Configuration options expected by this root isolation method. - type RootIsolationAlgorithmOptions ty (n :: Nat) (d :: Nat) = r | r -> ty n d + type RootIsolationAlgorithmOptions ty n d = r | r -> ty n d -- | Run one step of the root isolation method. -- -- This gets given the equations and a box, and should attempt to @@ -155,9 +158,7 @@ class ( Typeable ty, Show ( StepDescription ty ) ) -- bix does not contain any solutions; -- - (as a writer side-effect) boxes to definitely stop processing; see 'DoneBoxes'. rootIsolationAlgorithm - :: forall (n :: Nat) (d :: Nat) - . BoxCt n d - => RootIsolationAlgorithmOptions ty n d + :: RootIsolationAlgorithmOptions ty n d -- ^ options for this root isolation algorithm -> [ ( RootIsolationStep, Box n ) ] -- ^ history of the current round @@ -171,20 +172,22 @@ class ( Typeable ty, Show ( StepDescription ty ) ) -- | Match on an unknown root isolation algorithm step with a known algorithm. pattern IsolationStep - :: forall (ty :: Type) - . RootIsolationAlgorithm ty - => StepDescription ty -> RootIsolationStep + :: forall ( ty :: Type ) + . Typeable ty + => StepDescription ty + -> RootIsolationStep pattern IsolationStep stepDescr <- ( rootIsolationAlgorithmStep_maybe @ty -> Just stepDescr ) - where - IsolationStep stepDescr = SomeRootIsolationStep @ty stepDescr + -- NB: this pattern could also return @RootIsolationAlgorithm n d ty@ evidence, + -- but it's simpler to not do so for now. -- | Helper function used to define the 'IsolationStep' pattern. -- -- Inspects whether an existential 'RootIsolationStep' packs a step for -- the given algorithm. rootIsolationAlgorithmStep_maybe - :: forall ty. RootIsolationAlgorithm ty + :: forall ty + . Typeable ty => RootIsolationStep -> Maybe ( StepDescription ty ) rootIsolationAlgorithmStep_maybe ( SomeRootIsolationStep @existential descr ) | Just HRefl <- heqT @existential @ty diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Degree.hs b/brush-strokes/src/lib/Math/Root/Isolation/Degree.hs new file mode 100644 index 0000000..8546cc5 --- /dev/null +++ b/brush-strokes/src/lib/Math/Root/Isolation/Degree.hs @@ -0,0 +1,126 @@ +{-# LANGUAGE ScopedTypeVariables #-} + +module Math.Root.Isolation.Degree + ( -- * Computing the topological degree + topologicalDegree + ) + where + +-- base +import Data.Kind + ( Type ) +import GHC.TypeNats + ( Nat ) + +-- primitive +import Data.Primitive.SmallArray + +-- MetaBrush +import Math.Interval +import Math.Linear + +-------------------------------------------------------------------------------- +-- Topological degree. + +-- | Compute the topological degree of the given interval function (2D only). +topologicalDegree + :: Double + -- ^ minimum width when bisecting facets + -> ( 𝕀ℝ 2 -> 𝕀ℝ 2 ) + -- ^ function + -> 𝕀ℝ 2 + -- ^ box + -> Maybe Int +topologicalDegree Ξ΅_bis f ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = + topologicalDegreeFromContour $ + smallArrayFromList $ + concatMap ( tagFacet Ξ΅_bis f ) + [ ( Tag ( Fin 2 ) P, 𝕀 ( ℝ2 x_hi y_lo ) ( ℝ2 x_hi y_hi ) ) + , ( Tag ( Fin 1 ) N, 𝕀 ( ℝ2 x_lo y_hi ) ( ℝ2 x_hi y_hi ) ) + , ( Tag ( Fin 2 ) N, 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_lo y_hi ) ) + , ( Tag ( Fin 1 ) P, 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_lo ) ) ] + +-- | Tag a facet with the sign of one of the equations that is non-zero +-- on that facet. +-- +-- If all equations take the value 0 on that facet, bisect and recur. +tagFacet + :: Double + -> ( 𝕀ℝ 2 -> 𝕀ℝ 2 ) + -> ( Tag 2, 𝕀ℝ 2 ) + -> [ Tag 2 ] +tagFacet Ξ΅_bis f ( tag@( Tag i _ ), facet0 ) = go facet0 + where + go facet + | f1_b_lo > 0 + = [ Tag ( Fin 1 ) P ] + | f1_b_hi < 0 + = [ Tag ( Fin 1 ) N ] + | f2_b_lo > 0 + = [ Tag ( Fin 2 ) P ] + | f2_b_hi < 0 + = [ Tag ( Fin 2 ) N ] + | width ( facet `index` i ) < Ξ΅_bis + = [ NoTag ] -- TODO: shortcut the whole computation; + -- if we fail to tag one segment, the final answer could + -- be off by { -2, -1, 0, 1, 2 } (I believe) + | otherwise + = concatMap go $ bisectFacet tag facet + where + 𝕀 f1_b_lo f1_b_hi = f facet `index` Fin 1 + 𝕀 f2_b_lo f2_b_hi = f facet `index` Fin 2 +tagFacet _ _ ( NoTag, _ ) = error "impossible: input always tagged" + +bisectFacet :: Tag 2 -> 𝕀ℝ 2 -> [ 𝕀ℝ 2 ] +bisectFacet ( Tag i s ) x = + ( if s == P then id else reverse ) $ bisectInCoord i x +bisectFacet NoTag _ = error "impossible: input always tagged" + +bisectInCoord :: Fin n -> 𝕀ℝ 2 -> [ 𝕀ℝ 2 ] +bisectInCoord ( Fin 1 ) ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = + [ 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_mid y_hi ), 𝕀 ( ℝ2 x_mid y_lo ) ( ℝ2 x_hi y_hi ) ] + where x_mid = 0.5 * ( x_lo + x_hi ) +bisectInCoord _ ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = + [ 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_mid ), 𝕀 ( ℝ2 x_lo y_mid ) ( ℝ2 x_hi y_hi ) ] + where y_mid = 0.5 * ( y_lo + y_hi ) + +-- | Compute the topological degree of \( f \colon \mathbb{IR}^n \to \mathbb{IR}^n \) +-- on a box \( D \subset \mathbb{IR}^n \) by using the signs of the components +-- of \( f \) on a counter-clockwise polygonal contour describing the boundary of \( D \). +topologicalDegreeFromContour :: SmallArray ( Tag 2 ) -> Maybe Int +topologicalDegreeFromContour facetTags = go 0 0 + where + n = sizeofSmallArray facetTags + go :: Int -> Int -> Maybe Int + go !d !i + | i >= n + = Just d + | otherwise + = case indexSmallArray facetTags i of + NoTag -> Nothing + Tag ( Fin 1 ) P -> go ( d + Ξ΄ ) ( i + 1 ) + where + Ξ΄ = if indexSmallArray facetTags ( if i == n - 1 then 0 else i + 1 ) == Tag ( Fin 2 ) P + then 1 else 0 + - if indexSmallArray facetTags ( if i == 0 then n - 1 else i - 1 ) == Tag ( Fin 2 ) P + then 1 else 0 + _ -> go d ( i + 1 ) + +-------------------------------------------------------------------------------- +-- Tagging edges with signs of certain components of the function. + +data Sign = P | N + deriving stock ( Eq, Ord ) +instance Show Sign where + showsPrec _ = \case { P -> showString "+"; N -> showString "-" } +type Tag :: Nat -> Type +newtype Tag d = MkTag Int + deriving newtype ( Eq, Ord ) +pattern NoTag :: Tag d +pattern NoTag = MkTag 0 +pattern Tag :: Fin d -> Sign -> Tag d +pattern Tag c s <- ( getTag -> (# c, s #) ) + where Tag ( Fin c ) s = MkTag $ case s of { P -> fromIntegral c; N -> -( fromIntegral c ) } +{-# COMPLETE Tag, NoTag #-} +getTag :: Tag d -> (# Fin d, Sign #) +getTag ( MkTag t ) = (# Fin ( fromIntegral $ abs t ), if t >= 0 then P else N #) diff --git a/brush-strokes/src/lib/Math/Root/Isolation/GaussSeidel.hs b/brush-strokes/src/lib/Math/Root/Isolation/GaussSeidel.hs index 9fad426..b1b5ff1 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation/GaussSeidel.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation/GaussSeidel.hs @@ -1,4 +1,5 @@ -{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE UndecidableInstances #-} module Math.Root.Isolation.GaussSeidel ( -- * The interval Newton method with Gauss–Seidel step @@ -37,7 +38,7 @@ import GHC.TypeNats -- eigen import qualified Eigen.Matrix as Eigen ( Matrix - , determinant, generate, inverse, unsafeCoeff + , generate, inverse, unsafeCoeff ) -- transformers @@ -47,8 +48,6 @@ import Control.Monad.Trans.Writer.CPS -- MetaBrush import Math.Algebra.Dual ( D ) -import Math.Epsilon - ( nearZero ) import Math.Interval import Math.Linear import Math.Module @@ -62,13 +61,22 @@ import Math.Root.Isolation.Core -- | The interval Newton method with a Gauss–Seidel step; see 'intervalGaussSeidel'. data GaussSeidel -instance RootIsolationAlgorithm GaussSeidel where +instance BoxCt n d => RootIsolationAlgorithm GaussSeidel n d where type instance StepDescription GaussSeidel = () type instance RootIsolationAlgorithmOptions GaussSeidel n d = GaussSeidelOptions n d rootIsolationAlgorithm opts _thisRoundHist _prevRoundsHist eqs box = do - res <- intervalGaussSeidel opts eqs box + res <- intervalGaussSeidel @n @d opts eqs box return ( (), res ) {-# INLINEABLE rootIsolationAlgorithm #-} + {-# SPECIALISE rootIsolationAlgorithm + :: RootIsolationAlgorithmOptions GaussSeidel 2 3 + -> [ ( RootIsolationStep, Box 2 ) ] + -> BoxHistory 2 + -> ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) ) + -> Box 2 + -> Writer ( DoneBoxes 2 ) ( StepDescription GaussSeidel, [ Box 2 ] ) #-} + -- NB: including this to be safe. The specialiser seems to sometimes + -- be able to generate this specialisation on its own, and sometimes not. -- | Options for the interval Gauss–Seidel method. type GaussSeidelOptions :: Nat -> Nat -> Type @@ -144,13 +152,15 @@ intervalGaussSeidel eqs x | let x_mid = singleton $ boxMidpoint x + f :: 𝕀ℝ n -> 𝕀ℝ n + f = \ x_0 -> pickEqs $ eqs x_0 `monIndex` zeroMonomial + 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 ) + 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 @@ -160,18 +170,34 @@ intervalGaussSeidel -- at the origin, in order to take a Gauss–Seidel step. 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 - in do tell $ noDoneBoxes { doneSolBoxes = done } - return todo - where + ( done, todo ) = bimap ( map fst ) ( map fst ) + $ partition snd gsGuesses + 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. + do tell $ noDoneBoxes { doneSolBoxes = done } + return todo {-# INLINEABLE intervalGaussSeidel #-} +{- + +mbDeg = topologicalDegree 0.005 f x +det = case f'_x of + Vec [ c1, c2 ] -> + let a_11 = c1 `index` Fin 1 + a_12 = c2 `index` Fin 1 + a_21 = c1 `index` Fin 2 + a_22 = c2 `index` Fin 2 + in a_11 * a_22 - a_12 * a_21 + _ -> error "TODO: just testing n=2 here" + +if | not $ 0 ∈ det + , mbDeg == Just 0 + -> return [] + -- If the Jacobian is invertible over the box, then the topological + -- degree tells us exactly how many solutions there are in the box. +-} -- | 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. @@ -278,10 +304,7 @@ fromComponents f = do -- | The midpoint of a box. boxMidpoint :: Representable Double ( ℝ n ) => 𝕀ℝ n -> ℝ n boxMidpoint box = - tabulate $ \ i -> - let 𝕀 z_lo z_hi = index box i - z_mid = 0.5 * ( z_lo + z_hi ) - in z_mid + tabulate $ \ i -> midpoint ( box `index` i ) {-# INLINEABLE boxMidpoint #-} -- | Pre-condition the system \( AX = B \). @@ -298,14 +321,10 @@ precondition meth as b = -> ( as, b ) InverseMidpoint | let mat = toEigen $ fmap boxMidpoint as - det = Eigen.determinant mat - , not $ nearZero det - -- (TODO: a bit wasteful to compute determinant then inverse.) , let precond = Eigen.inverse mat doPrecond = matMulVec ( fromEigen precond ) + -- TODO: avoid this when condition number is small? -> ( fmap doPrecond as, doPrecond b ) - | otherwise - -> ( as, b ) where toEigen :: Vec n ( ℝ n ) -> Eigen.Matrix n n Double toEigen cols = diff --git a/brush-strokes/src/lib/Math/Root/Isolation/Narrowing.hs b/brush-strokes/src/lib/Math/Root/Isolation/Narrowing.hs index a41ace9..05c04fb 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation/Narrowing.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation/Narrowing.hs @@ -1,4 +1,5 @@ -{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE UndecidableInstances #-} module Math.Root.Isolation.Narrowing ( -- * @box(1)@-consistency @@ -39,6 +40,8 @@ import GHC.TypeNats -- transformers import Control.Monad.Trans.State.Strict as State ( State, evalState, get, put ) +import Control.Monad.Trans.Writer.CPS + ( Writer ) -- brush-strokes import Math.Algebra.Dual @@ -62,21 +65,39 @@ import Math.Root.Isolation.Core -- | A @box(1)@-consistency enforcing algorithm; see 'makeBox1Consistent'. data Box1 -instance RootIsolationAlgorithm Box1 where +instance BoxCt n d => RootIsolationAlgorithm Box1 n d where type instance StepDescription Box1 = () type instance RootIsolationAlgorithmOptions Box1 n d = Box1Options n d rootIsolationAlgorithm opts _thisRoundHist _prevRoundsHist eqs box = return $ ( (), makeBox1Consistent opts eqs box ) {-# INLINEABLE rootIsolationAlgorithm #-} + {-# SPECIALISE rootIsolationAlgorithm + :: RootIsolationAlgorithmOptions Box1 2 3 + -> [ ( RootIsolationStep, Box 2 ) ] + -> BoxHistory 2 + -> ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) ) + -> Box 2 + -> Writer ( DoneBoxes 2 ) ( StepDescription Box1, [ Box 2 ] ) #-} + -- NB: including this to be safe. The specialiser seems to sometimes + -- be able to generate this specialisation on its own, and sometimes not. -- | A @box(2)@-consistency enforcing algorithm; see 'makeBox1Consistent'. data Box2 -instance RootIsolationAlgorithm Box2 where +instance BoxCt n d => RootIsolationAlgorithm Box2 n d where type instance StepDescription Box2 = () type instance RootIsolationAlgorithmOptions Box2 n d = Box2Options n d rootIsolationAlgorithm opts _thisRoundHist _prevRoundsHist eqs box = return ( () , [ makeBox2Consistent opts eqs box ] ) {-# INLINEABLE rootIsolationAlgorithm #-} + {-# SPECIALISE rootIsolationAlgorithm + :: RootIsolationAlgorithmOptions Box2 2 3 + -> [ ( RootIsolationStep, Box 2 ) ] + -> BoxHistory 2 + -> ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) ) + -> Box 2 + -> Writer ( DoneBoxes 2 ) ( StepDescription Box2, [ Box 2 ] ) #-} + -- NB: including this to be safe. The specialiser seems to sometimes + -- be able to generate this specialisation on its own, and sometimes not. -- | Options for the @box(1)@-consistency method. data Box1Options n d = diff --git a/src/app/MetaBrush/Application.hs b/src/app/MetaBrush/Application.hs index 946b0b1..58d5e77 100644 --- a/src/app/MetaBrush/Application.hs +++ b/src/app/MetaBrush/Application.hs @@ -74,9 +74,7 @@ import Control.Monad.Trans.Reader -- MetaBrush import Math.Root.Isolation - ( RootIsolationOptions(..), defaultRootIsolationOptions - , N - ) + ( RootIsolationOptions(..), defaultRootIsolationOptions ) import Math.Bezier.Cubic.Fit ( FitParameters(..) ) import Math.Bezier.Spline @@ -208,18 +206,18 @@ runApplication application = do maxHistorySizeTVar <- STM.newTVarIO @Int 1000 fitParametersTVar <- STM.newTVarIO @FitParameters $ FitParameters - { maxSubdiv = 5 --2 --3 -- 6 - , nbSegments = 3 --3 --6 -- 12 - , dist_tol = 0.1 -- 5e-3 - , t_tol = 0.1 -- 1e-4 - , maxIters = 5 -- 100 + { maxSubdiv = 0 --5 --2 --3 -- 6 + , nbSegments = 1 --5 + , dist_tol = 5e-3 + , t_tol = 1e-4 + , maxIters = 20 } rootsAlgoTVar <- STM.newTVarIO @RootSolvingAlgorithm $ --HalleyM2 NewtonRaphson { maxIters = 20, precision = 8 } - cuspFindingOptionsTVar <- STM.newTVarIO @( Maybe ( RootIsolationOptions N 3 ) ) $ - Just defaultRootIsolationOptions + cuspFindingOptionsTVar <- STM.newTVarIO @( Maybe ( RootIsolationOptions 2 3 ) ) $ + Nothing --Just defaultRootIsolationOptions -- Put all these stateful variables in a record for conciseness. let diff --git a/src/app/MetaBrush/Render/Document.hs b/src/app/MetaBrush/Render/Document.hs index 4cb9f7f..c73a632 100644 --- a/src/app/MetaBrush/Render/Document.hs +++ b/src/app/MetaBrush/Render/Document.hs @@ -83,7 +83,7 @@ import Math.Linear import Math.Module ( Module((*^)), normalise ) import Math.Root.Isolation - ( RootIsolationOptions, N ) + ( RootIsolationOptions ) import MetaBrush.Asset.Colours ( Colours, ColourRecord(..) ) import MetaBrush.Brush @@ -153,7 +153,7 @@ blankRender _ = pure () getDocumentRender :: Colours - -> RootSolvingAlgorithm -> Maybe ( RootIsolationOptions N 3 ) -> FitParameters + -> RootSolvingAlgorithm -> Maybe ( RootIsolationOptions 2 3 ) -> FitParameters -> Mode -> Bool -> Set Modifier -> Maybe ( ℝ 2 ) -> Maybe HoldAction -> Maybe PartialPath -> Document @@ -289,7 +289,7 @@ instance NFData StrokeRenderData where -- - Otherwise, this consists of the underlying spline path only. strokeRenderData :: RootSolvingAlgorithm - -> Maybe ( RootIsolationOptions N 3 ) + -> Maybe ( RootIsolationOptions 2 3 ) -> FitParameters -> Stroke -> Maybe ( ST RealWorld StrokeRenderData ) @@ -304,7 +304,7 @@ strokeRenderData rootAlgo mbCuspOptions fitParams Just ( NamedBrush { brushFunction = fn } ) | WithParams { defaultParams = brush_defaults - , withParams = brush@( Brush { brushShape, mbRotation = mbRot } ) + , withParams = brush@( Brush { brushBaseShape, mbRotation = mbRot } ) } <- fn -> -- This is the key place where we need to perform impedance matching -- between the collection of parameters supplied along a stroke and @@ -328,7 +328,7 @@ strokeRenderData rootAlgo mbCuspOptions fitParams , strokeBrushFunction = \ params -> let brushParams = embedUsedParams $ toUsedParams params - shape = fun @Double brushShape brushParams + shape = fun @Double brushBaseShape brushParams -- TODO: remove this logic which is duplicated -- from elsewhere. The type should make it -- impossible to forget to apply the rotation.