mirror of
https://gitlab.com/sheaf/metabrush.git
synced 2024-11-27 17:34:08 +00:00
Improve specialisation + add degree computation
This commit is contained in:
parent
6450859e3c
commit
c89fba7fd2
|
@ -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
|
||||
|
|
|
@ -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
|
||||
|
|
|
@ -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
|
||||
}
|
||||
]
|
||||
)
|
||||
|
|
|
@ -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 ) )
|
||||
}
|
||||
|
||||
--------------------------------------------------------------------------------
|
||||
|
|
|
@ -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.
|
||||
|
|
|
@ -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 #-}
|
||||
|
|
|
@ -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.
|
||||
|
|
|
@ -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 )
|
||||
|
|
|
@ -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
|
||||
|
|
|
@ -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 #-}
|
||||
|
|
|
@ -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
|
||||
|
|
|
@ -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
|
||||
|
|
126
brush-strokes/src/lib/Math/Root/Isolation/Degree.hs
Normal file
126
brush-strokes/src/lib/Math/Root/Isolation/Degree.hs
Normal file
|
@ -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 #)
|
|
@ -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 =
|
||||
|
|
|
@ -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 =
|
||||
|
|
|
@ -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
|
||||
|
|
|
@ -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.
|
||||
|
|
Loading…
Reference in a new issue