Make it easier to switch between 2 and 3-dim

This commit makes it easier to switch between the 2-dim and 3-dim
formulations of the cusp-finding problem. This is still work in progress,
trying to improve the performance of the cusp-finding algorithms.
This commit is contained in:
sheaf 2024-04-18 21:33:55 +02:00
parent 0160081e80
commit 2b167f594a
4 changed files with 85 additions and 29 deletions

View file

@ -35,7 +35,7 @@ import GHC.Exts
import GHC.Generics
( Generic )
import GHC.TypeNats
( type (-) )
( Nat, type (-) )
import Numeric
( showFFloat )
@ -126,8 +126,8 @@ benchTestCase ( TestCase { testName, testBrushStroke, testCuspOptions, testStart
IntMap.toList $
findCuspsIn testCuspOptions testStrokeFnI $
IntMap.fromList
[ ( i, [ box ] )
| ( i, box ) <- testStartBoxes
[ ( i, boxes )
| ( i, boxes ) <- testStartBoxes
]
rnf dunno `seq` rnf sols `seq` return ()
after <- getMonotonicTime
@ -199,8 +199,8 @@ data TestCase =
TestCase
{ testName :: String
, testBrushStroke :: BrushStroke
, testCuspOptions :: RootIsolationOptions 2 3
, testStartBoxes :: [ ( Int, Box 2 ) ]
, testCuspOptions :: RootIsolationOptions N 3
, testStartBoxes :: [ ( Int, [ Box 2 ] ) ]
}
brushStrokeFunctions
@ -488,7 +488,7 @@ showD float = showFFloat (Just 6) float ""
--------------------------------------------------------------------------------
ellipseTestCase :: RootIsolationOptions 2 3 -> String -> ( Double, Double ) -> Double -> [ ( Int, Box 2 ) ] -> TestCase
ellipseTestCase :: RootIsolationOptions N 3 -> String -> ( Double, Double ) -> Double -> [ ( Int, [ Box 2 ] ) ] -> TestCase
ellipseTestCase opts str k0k1 rot startBoxes =
TestCase
{ testName = "ellipse (" ++ str ++ ")"
@ -607,11 +607,15 @@ getStrokeFunctions brush sp0 crv =
brush @3 @𝕀 proxy# singleton )
{-# INLINEABLE getStrokeFunctions #-}
defaultStartBoxes :: [ Int ] -> [ ( Int, Box 2 ) ]
defaultStartBoxes :: [ Int ] -> [ ( Int, [ Box 2 ] ) ]
defaultStartBoxes is =
[ ( i, 𝕀 ( 2 zero zero ) ( 2 one one ) )
| i <- is
]
[ ( i, [ 𝕀 ( 2 zero zero ) ( 2 one one ) ] ) | i <- is ]
-- [ ( i, [ 𝕀 ( 2 t s ) ( 2 ( t + 0.099999 ) ( s + 0.099999 ) )
-- | t <- [ 0.00001, 0.1 .. 0.9 ]
-- , s <- [ 0.00001, 0.1 .. 0.9 ]
-- ] )
-- | i <- is
-- ]
getR1 (1 u) = u

View file

@ -18,6 +18,9 @@ module Math.Bezier.Stroke
-- ** Cusp finding
, findCusps, findCuspsIn
-- TODO: hack for switching between 2 and 3 dim formulations of cusp-finding
, N
)
where
@ -243,7 +246,7 @@ computeStrokeOutline ::
)
=> RootSolvingAlgorithm
-> Maybe ( RootIsolationOptions 2 3 )
-> Maybe ( RootIsolationOptions N 3 )
-> FitParameters
-> ( ptData -> usedParams )
-> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing
@ -522,7 +525,7 @@ outlineFunction
, Show ptData, Show brushParams
)
=> RootSolvingAlgorithm
-> Maybe ( RootIsolationOptions 2 3 )
-> Maybe ( RootIsolationOptions N 3 )
-> ( ptData -> usedParams )
-> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing
-> ( forall {t} k (i :: t)
@ -1147,9 +1150,9 @@ solveEnvelopeEquations rootAlgo _t path_t path'_t ( fwdOffset, bwdOffset ) strok
-- TODO: use Newton's method starting at the midpoint of the box,
-- instead of just taking the midpoint of the box.
cuspCoords :: ( 1 -> Seq ( 1 -> StrokeDatum 2 () ) )
-> ( Int, Box 2 )
-> ( Int, Box N )
-> Cusp
cuspCoords eqs ( i, 𝕀 ( 2 t_lo s_lo ) ( 2 t_hi s_hi ) )
cuspCoords eqs ( i, box )
| StrokeDatum
{ dpath
, dstroke = D22 { _D22_v = stroke } }
@ -1160,17 +1163,21 @@ cuspCoords eqs ( i, 𝕀 ( 2 t_lo s_lo ) ( 2 t_hi s_hi ) )
, cuspStrokeCoords = coerce stroke
}
where
𝕀 t_lo t_hi = box `index` Fin 1
𝕀 s_lo s_hi = box `index` Fin 2
t_mid = 0.5 * ( t_lo + t_hi )
s_mid = 0.5 * ( s_lo + s_hi )
type N = 2
-- | Find cusps in the envelope for values of the parameters in
-- \( 0 \leqslant t, s \leqslant 1 \), using interval arithmetic.
--
-- See 'isolateRootsIn' for details of the algorithms used.
findCusps
:: RootIsolationOptions 2 3
:: RootIsolationOptions N 3
-> ( 𝕀 1 -> Seq ( 𝕀 1 -> StrokeDatum 3 𝕀 ) )
-> IntMap ( [ ( Box 2, RootIsolationTree ( Box 2 ) ) ], ( [ Box 2 ], [ Box 2 ] ) )
-> IntMap ( [ ( Box N, RootIsolationTree ( Box N ) ) ], ( [ Box N ], [ Box N ] ) )
findCusps opts boxStrokeData =
findCuspsIn opts boxStrokeData $
IntMap.fromList
@ -1187,10 +1194,10 @@ findCusps opts boxStrokeData =
-- | Like 'findCusps', but parametrised over the initial boxes for the
-- root isolation method.
findCuspsIn
:: RootIsolationOptions 2 3
:: RootIsolationOptions N 3
-> ( 𝕀 1 -> Seq ( 𝕀 1 -> StrokeDatum 3 𝕀 ) )
-> IntMap [ Box 2 ]
-> IntMap ( [ ( Box 2, RootIsolationTree ( Box 2 ) ) ], ( [ Box 2 ], [ Box 2 ] ) )
-> IntMap ( [ ( Box N, RootIsolationTree ( Box N ) ) ], ( [ Box N ], [ Box N ] ) )
findCuspsIn opts boxStrokeData initBoxes =
IntMap.mapWithKey ( \ i -> foldMap ( isolateRootsIn opts ( eqnPiece i ) ) ) initBoxes
where
@ -1214,13 +1221,42 @@ findCuspsIn opts boxStrokeData initBoxes =
{-
findCuspsIn
:: RootIsolationOptions 3 3
:: RootIsolationOptions N 3
-> ( 𝕀 1 -> Seq ( 𝕀 1 -> StrokeDatum 3 𝕀 ) )
-> IntMap [ Box 3 ]
-> IntMap [ Box 2 ]
-> IntMap ( [ ( Box 3, RootIsolationTree ( Box 3 ) ) ], ( [ Box 3 ], [ Box 3 ] ) )
findCuspsIn opts boxStrokeData initBoxes =
IntMap.mapWithKey ( \ i -> foldMap ( isolateRootsIn opts ( eqnPiece i ) ) ) initBoxes
IntMap.mapWithKey ( \ i boxes -> foldMap ( isolateRootsIn opts ( eqnPiece i ) ) $ concatMap ( mkInitBox i ) boxes ) initBoxes
where
mkInitBox :: Int -> Box 2 -> [ Box 3 ]
mkInitBox 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 )
StrokeDatum
{ dstroke =
D32
{ _D32_dx = T ( 𝕀 ( 2 cx_t_lo cy_t_lo ) ( 2 cx_t_hi cy_t_hi ) )
, _D32_dy = T ( 𝕀 ( 2 cx_s_lo cy_s_lo ) ( 2 cx_s_hi cy_s_hi ) )}
, ee =
D22
{ _D22_dx = T ( 𝕀 ( 1 ee_t_lo ) ( 1 ee_t_hi ) )
, _D22_dy = T ( 𝕀 ( 1 ee_s_lo ) ( 1 ee_s_hi ) ) }
} = ( boxStrokeData t `Seq.index` i ) s
-- λ = ∂E/∂t / ∂E/∂s
λ1 = 𝕀 ee_t_lo ee_t_hi `extendedDivide` 𝕀 ee_s_lo ee_s_hi
-- λ = ∂c/∂t / ∂c/∂s
λ2 = 𝕀 cx_t_lo cx_t_hi `extendedDivide` 𝕀 cx_s_lo cx_s_hi
λ3 = 𝕀 cy_t_lo cy_t_hi `extendedDivide` 𝕀 cy_s_lo cy_s_hi
λ = [ 𝕀 ( recip -0 ) ( recip 0 ) ]
`intersectMany` λ1
`intersectMany` λ2
`intersectMany` λ3
in
let boxes = [ 𝕀 ( 3 t_lo s_lo λ_lo ) ( 3 t_hi s_hi λ_hi )
| 𝕀 λ_lo' λ_hi' <- λ
, let λ_lo = max λ_lo' ( min ( λ_hi' - 10 ) -100 )
λ_hi = min λ_hi' ( max ( λ_lo' + 10 ) 100 ) ]
in boxes
eqnPiece :: Int -> 𝕀 ( 3 ) -> D1𝔸3 ( 𝕀 ( 3 ) )
eqnPiece i ( 𝕀 ( 3 t_lo s_lo λ_lo ) ( 3 t_hi s_hi λ_hi ) ) =
let t = 𝕀 ( 1 t_lo ) ( 1 t_hi )
@ -1253,5 +1289,20 @@ findCuspsIn opts boxStrokeData initBoxes =
( T $ 𝕀 ( 3 f1_s_lo f2_s_lo f3_s_lo ) ( 3 f1_s_hi f2_s_hi f3_s_hi ) )
( T $ 𝕀 ( 3 f1_λ_lo f2_λ_lo f3_λ_lo ) ( 3 f1_λ_hi f2_λ_hi f3_λ_hi ) )
intersectMany :: [𝕀 Double] -> [𝕀 Double] -> [𝕀 Double]
intersectMany _ [] = []
intersectMany is (j : js) = intersectOne is j ++ intersectMany is js
intersectOne :: [ 𝕀 Double ] -> 𝕀 Double -> [ 𝕀 Double ]
intersectOne is i = concatMap ( intersect i ) is
intersect :: 𝕀 Double -> 𝕀 Double -> [ 𝕀 Double ]
intersect ( 𝕀 lo1 hi1 ) ( 𝕀 lo2 hi2 )
| lo > hi
= [ ]
| otherwise
= [ 𝕀 lo hi ]
where
lo = max lo1 lo2
hi = min hi1 hi2
-}

View file

@ -18,6 +18,8 @@ module Math.Interval
-- base
import Prelude hiding ( Num(..), Fractional(..) )
import Data.List
( nub )
-- acts
import Data.Act
@ -78,7 +80,7 @@ instance Torsor ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where
-- Extended division
extendedDivide :: 𝕀 Double -> 𝕀 Double -> [ 𝕀 Double ]
extendedDivide x y = map ( x * ) ( extendedRecip y )
extendedDivide x y = nub $ map ( x * ) ( extendedRecip y )
extendedRecip :: 𝕀 Double -> [ 𝕀 Double ]
extendedRecip x@( 𝕀 lo hi )

View file

@ -41,7 +41,7 @@ import Data.Kind
import Data.Foldable
( toList )
import Data.List
( nub, partition, sort )
( partition, sort )
import Data.List.NonEmpty
( NonEmpty )
import qualified Data.List.NonEmpty as NE
@ -81,7 +81,7 @@ import Math.Linear
import Math.Module
( Module(..) )
import Math.Monomial
( MonomialBasis(..), Deg, Vars
( MonomialBasis(..)
, linearMonomial, zeroMonomial
)
import qualified Math.Ring as Ring
@ -202,7 +202,7 @@ data Box2Options =
, box2LambdaMin :: !Double
}
defaultRootIsolationOptions :: RootIsolationOptions 2 3
defaultRootIsolationOptions :: BoxCt n d => RootIsolationOptions n d
defaultRootIsolationOptions =
RootIsolationOptions
{ minWidth
@ -251,9 +251,8 @@ defaultGaussSeidelOptions :: GaussSeidelOptions 2 3
defaultGaussSeidelOptions =
GaussSeidelOptions
{ gsPreconditioner = InverseMidJacobian
, gsDims = --id
\ ( 𝕀 ( 3 _a_lo b_lo c_lo ) ( 3 _a_hi b_hi c_hi ) )
-> 𝕀 ( 2 b_lo c_lo ) ( 2 b_hi c_hi )
, gsDims = \ ( 𝕀 ( 3 _a_lo b_lo c_lo ) ( 3 _a_hi b_hi c_hi ) )
-> 𝕀 ( 2 b_lo c_lo ) ( 2 b_hi c_hi )
}
defaultBisectionOptions
@ -576,7 +575,7 @@ gaussSeidelStep
-> 𝕀 n -- ^ \( B \)
-> T ( 𝕀 n ) -- ^ initial box \( X \)
-> [ ( T ( 𝕀 n ), Bool ) ]
gaussSeidelStep as b ( T x0 ) = coerce $ nub $
gaussSeidelStep as b ( T x0 ) = coerce $
forEachDim @n ( x0, True ) $ \ i ( x, contraction ) -> do
-- x_i' = ( b_i - sum { j /= i } a_ij * x_j ) / a_ii
x_i'0 <- ( b `index` i - sum [ ( as ! j ) `index` i * x `index` j | j <- toList ( universe @n ), j /= i ] )