From 2b167f594abe560e9c8d6029723d9b635fb5da63 Mon Sep 17 00:00:00 2001 From: sheaf Date: Thu, 18 Apr 2024 21:33:55 +0200 Subject: [PATCH] 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. --- brush-strokes/src/cusps/bench/Main.hs | 24 ++++--- brush-strokes/src/lib/Math/Bezier/Stroke.hs | 73 +++++++++++++++++--- brush-strokes/src/lib/Math/Interval.hs | 4 +- brush-strokes/src/lib/Math/Root/Isolation.hs | 13 ++-- 4 files changed, 85 insertions(+), 29 deletions(-) diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs index 7bd4867..82dfd88 100644 --- a/brush-strokes/src/cusps/bench/Main.hs +++ b/brush-strokes/src/cusps/bench/Main.hs @@ -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 diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke.hs b/brush-strokes/src/lib/Math/Bezier/Stroke.hs index 60e6e64..cba086e 100644 --- a/brush-strokes/src/lib/Math/Bezier/Stroke.hs +++ b/brush-strokes/src/lib/Math/Bezier/Stroke.hs @@ -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 -} \ No newline at end of file diff --git a/brush-strokes/src/lib/Math/Interval.hs b/brush-strokes/src/lib/Math/Interval.hs index 0ddbe84..39d7c0d 100644 --- a/brush-strokes/src/lib/Math/Interval.hs +++ b/brush-strokes/src/lib/Math/Interval.hs @@ -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 ) diff --git a/brush-strokes/src/lib/Math/Root/Isolation.hs b/brush-strokes/src/lib/Math/Root/Isolation.hs index 5d18c44..7613a52 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation.hs @@ -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 ] )