From 26cfdada8f58c5b9b31e1521baab9ebdbeb5d14f Mon Sep 17 00:00:00 2001 From: sheaf Date: Sat, 24 Feb 2024 19:37:34 +0100 Subject: [PATCH] Fix extendedRecip (negative infinity) --- brush-strokes/bench/Main.hs | 98 +++++++++++++++++++++ brush-strokes/src/Math/Bezier/Stroke.hs | 21 ++--- brush-strokes/src/Math/Interval.hs | 23 ++++- brush-strokes/src/Math/Interval/Internal.hs | 12 ++- 4 files changed, 136 insertions(+), 18 deletions(-) diff --git a/brush-strokes/bench/Main.hs b/brush-strokes/bench/Main.hs index 3df3cb3..929fb61 100644 --- a/brush-strokes/bench/Main.hs +++ b/brush-strokes/bench/Main.hs @@ -17,16 +17,22 @@ module Main where -- base +import Control.Concurrent.MVar + ( newMVar ) import Data.Coerce ( coerce ) import Data.Foldable ( for_ ) +import Data.List + ( intercalate ) import GHC.Exts ( Proxy#, proxy# ) import GHC.Generics ( Generic ) import GHC.TypeNats ( type (-) ) +import Numeric + ( showFFloat ) -- containers import Data.Sequence @@ -38,6 +44,8 @@ import Data.Tree -- brush-strokes import Calligraphy.Brushes +import Debug.Utils + ( logToFile ) import Math.Algebra.Dual import Math.Bezier.Spline import Math.Bezier.Stroke @@ -66,6 +74,9 @@ main = for_ testCases $ \ testCase@( TestCase { testName, testAlgorithmParams } , " dunno: " ++ show dunno , " sols: " ++ show sols ] + --logFileMVar <- newMVar "logs/trickyCusp.log" + --logToFile logFileMVar (unlines logLines) + -- `seq` return () testCases :: [ TestCase ] testCases = [ ellipse , trickyCusp2 ] @@ -140,9 +151,96 @@ putStrLn $ unlines $ map Data.Tree.View.showTree $ showTrees $ mkBox (0.5798, 0. ([ℝ1 0.5798, ℝ1 0.675],3,[ℝ1 0.26798, ℝ1 0.26799]) (area 0.000001) N [] └─ ([ℝ1 0.5973000285624527, ℝ1 0.6750000000000002],3,[ℝ1 0.26798, ℝ1 0.26799000000000006]) (area 0.000001) NoSolution "ee" ([ℝ1 0.5973000285624527, ℝ1 0.6750000000000002],3,[ℝ1 0.26798, ℝ1 0.26799000000000006]) +eval fI $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799) +> D12 {_D12_v = T[ℝ2 -10088.6674944889 -3281.3820867312834, ℝ2 4124.668381545453 4524.807156085763], _D12_dx = TT[ℝ2 -173746.97965005718 -33281.18494907289, ℝ2 298.2609121556852 23639.772884799597], _D12_dy = TT[ℝ2 -18454.27716258352 -28337.509817580823, ℝ2 1163.6949532017436 -13936.383137525536]}} +i.e. +> f = [ℝ2 -10088.6674944889 -3281.3820867312834, ℝ2 4124.668381545453 4524.807156085763] +> f_t = [ℝ2 -173746.97965005718 -33281.18494907289, ℝ2 298.2609121556852 23639.772884799597] +> f_s = [ℝ2 -18454.27716258352 -28337.509817580823, ℝ2 1163.6949532017436 -13936.383137525536] + +(f, fI) = testCaseStrokeFunctions trickyCusp2 +t = 𝕀 (ℝ1 0.5798) (ℝ1 0.675) +s = 𝕀 (ℝ1 0.26798) (ℝ1 0.26799) +t_mid = 0.5 * ( 0.5798 + 0.675 ) +s_mid = 0.5 * ( 0.26798 + 0.26799 ) +D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) = dEdsdcdt $ eval fI (t, 3, s) +t' = coerce ( (-) @( 𝕀 Double ) ) t ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1 +s' = coerce ( (-) @( 𝕀 Double ) ) s ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1 +a = ( f_t, f_s ) +b = negV2 $ singleton $ midV2 f +[((t2', s2'), isContr)] = gaussSeidel a b (t', s') +t2 = coerce ( (+) @( 𝕀 Double ) ) t2' ( singleton ( ℝ1 t_mid ) ) :: 𝕀ℝ 1 +s2 = coerce ( (+) @( 𝕀 Double ) ) s2' ( singleton ( ℝ1 s_mid ) ) :: 𝕀ℝ 1 + +t2 +> [ℝ1 0.6102365832093095, ℝ1 0.6750000000000002] +s2 +> [ℝ1 0.26798, ℝ1 0.26799000000000006] + + +let ( 𝕀 ( ℝ2 a11_lo a21_lo ) ( ℝ2 a11_hi a21_hi ), 𝕀 ( ℝ2 a12_lo a22_lo ) ( ℝ2 a12_hi a22_hi ) ) = a +let ( 𝕀 ( ℝ2 b1_lo b2_lo ) ( ℝ2 b1_hi b2_hi ) ) = b +let ( 𝕀 ( ℝ1 x1_lo ) ( ℝ1 x1_hi ), 𝕀 ( ℝ1 x2_lo ) ( ℝ1 x2_hi ) ) = ( t', s' ) + +a11 = 𝕀 a11_lo a11_hi +a12 = 𝕀 a12_lo a12_hi +a21 = 𝕀 a21_lo a21_hi +a22 = 𝕀 a22_lo a22_hi +b1 = 𝕀 b1_lo b1_hi +b2 = 𝕀 b2_lo b2_hi +x1 = 𝕀 x1_lo x1_hi +x2 = 𝕀 x2_lo x2_hi + +( b1 - a12 * x2 ) +> [2981.90728508591, 2982.0918278575364] + +extendedRecip a11 -} +negV2 :: 𝕀ℝ 2 -> 𝕀ℝ 2 +negV2 ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = + let !( 𝕀 x'_lo x'_hi ) = negate $ 𝕀 x_lo x_hi + !( 𝕀 y'_lo y'_hi ) = negate $ 𝕀 y_lo y_hi + in 𝕀 ( ℝ2 x'_lo y'_lo ) ( ℝ2 x'_hi y'_hi ) + +midV2 :: 𝕀ℝ 2 -> ℝ 2 +midV2 ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = + ℝ2 ( 0.5 * ( x_lo + x_hi ) ) ( 0.5 * ( y_lo + y_hi ) ) + +logLines :: [ String ] +logLines = + [ "f = dE/ds * dc/dt: f, df/dt, df/ds" + , "{" ++ + (intercalate "," + [ "{" ++ showD (midPoint t) ++ "," ++ showD (midPoint s) ++ ",{" ++ intercalate "," vals ++ "}}" + | t <- map ( around 0.5798 ) [-0.05, -0.049.. 0.05] + , let i = 3 + , s <- map ( around 0.26798 ) [-0.05, -0.049.. 0.05] + , let StrokeDatum + { 𝛿E𝛿sdcdt = D12 (T f) (T (T f_t)) (T (T f_s)) + } = (curvesI t `Seq.index` i) s + ℝ2 vx vy = midPoint2 f + ℝ2 vx_t vy_t = midPoint2 f_t + ℝ2 vx_s vy_s = midPoint2 f_s + vals = [ "{" ++ showD vx ++ "," ++ showD vy ++ "}" + , "{" ++ showD vx_t ++ "," ++ showD vy_t ++ "}" + , "{" ++ showD vx_s ++ "," ++ showD vy_s ++ "}" + ] + ] + ) ++ "}" + ] + where + around :: Double -> Double -> 𝕀ℝ 1 + around z0 z = 𝕀 ( ℝ1 ( z + z0 - 1e-6 ) ) ( ℝ1 ( z + z0 + 1e-6 ) ) + ( _, curvesI ) = testCaseStrokeFunctions trickyCusp2 + midPoint (𝕀 (ℝ1 lo) (ℝ1 hi)) = 0.5 * ( lo + hi ) + midPoint2 (𝕀 (ℝ2 lo_x lo_y) (ℝ2 hi_x hi_y)) + = ℝ2 ( 0.5 * ( lo_x + hi_x ) ) ( 0.5 * ( lo_y + hi_y ) ) + +showD :: Double -> String +showD float = showFFloat (Just 6) float "" + -------------------------------------------------------------------------------- ellipse :: TestCase diff --git a/brush-strokes/src/Math/Bezier/Stroke.hs b/brush-strokes/src/Math/Bezier/Stroke.hs index 0c138cc..0db8ed5 100644 --- a/brush-strokes/src/Math/Bezier/Stroke.hs +++ b/brush-strokes/src/Math/Bezier/Stroke.hs @@ -24,7 +24,7 @@ module Math.Bezier.Stroke , IntervalNewtonStep(..) , IntervalNewtonLeaf(..) , Box - , intervalNewtonGS, intervalNewtonGSFrom + , intervalNewtonGS, intervalNewtonGSFrom, gaussSeidel ) where @@ -1288,18 +1288,6 @@ intersect ( 𝕀 lo1 hi1 ) ( 𝕀 lo2 hi2 ) lo = max lo1 lo2 hi = min hi1 hi2 -extendedDivide :: 𝕀 Double -> 𝕀 Double -> [ 𝕀 Double ] -extendedDivide x y = map ( x * ) ( extendedRecip y ) - -extendedRecip :: 𝕀 Double -> [ 𝕀 Double ] -extendedRecip x@( 𝕀 lo hi ) - | lo == 0 && hi == 0 - = [ 𝕀 ( -1 / 0 ) ( 1 / 0 ) ] - | lo >= 0 || hi <= 0 - = [ recip x ] - | otherwise - = [ recip ( 𝕀 lo 0 ), recip ( 𝕀 0 hi ) ] - -- | Computes the brush stroke coordinates of a cusp from -- the @(t,s)@ parameter values. cuspCoords :: ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum 2 () ) ) @@ -1421,7 +1409,7 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox = <- ( eqs t `Seq.index` i ) s , StrokeDatum { ee = D22 _ee_mid _ _ _ _ _ - , 𝛿E𝛿sdcdt = D12 ( T f_mid ) ( T ( T f_t_mid ) ) ( T ( T f_s_mid ) ) } + , 𝛿E𝛿sdcdt = D12 ( T f_mid ) ( T ( T _f_t_mid ) ) ( T ( T _f_s_mid ) ) } <- ( eqs i_t_mid `Seq.index` i ) i_s_mid , let ee_potential_zero = inf ee <= ℝ1 0 && sup ee >= ℝ1 0 𝛿E𝛿sdcdt_potential_zero = cmpℝ2 (<=) ( inf f ) ( ℝ2 0 0 ) && cmpℝ2 (>=) ( sup f ) ( ℝ2 0 0 ) @@ -1433,7 +1421,7 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox = -- for the equation f'(x) v = - f(x_mid), -- where f = 𝛿E/𝛿s * dc/dt !( a, b ) = precondition precondMethod - ( f_t_mid, f_s_mid ) + ( midI f_t, midI f_s ) ( f_t, f_s ) ( neg f_mid ) --(a, b) -- | 𝕀 (ℝ1 ee_lo) (ℝ1 ee_hi) <- ee_mid @@ -1482,6 +1470,9 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox = | otherwise -> return [ IntervalNewtonLeaf $ NoSolution (if ee_potential_zero then "dc/dt" else "ee") cand ] where + midI :: 𝕀ℝ 2 -> 𝕀ℝ 2 + midI ( 𝕀 ( ℝ2 x_lo y_lo ) ( ℝ2 x_hi y_hi ) ) = + singleton $ ℝ2 ( 0.5 * ( x_lo + x_hi ) ) ( 0.5 * ( y_lo + y_hi ) ) t_mid = 0.5 * ( t_lo + t_hi ) s_mid = 0.5 * ( s_lo + s_hi ) i_t_mid = singleton ( ℝ1 t_mid ) diff --git a/brush-strokes/src/Math/Interval.hs b/brush-strokes/src/Math/Interval.hs index f3f2d4a..70e6232 100644 --- a/brush-strokes/src/Math/Interval.hs +++ b/brush-strokes/src/Math/Interval.hs @@ -11,11 +11,12 @@ module Math.Interval , singleton, nonDecreasing , inside , aabb + , extendedDivide, extendedRecip ) where -- base -import Prelude hiding ( Num(..) ) +import Prelude hiding ( Num(..), Fractional(..) ) -- acts import Data.Act @@ -69,6 +70,26 @@ instance Act ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where instance Torsor ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where a --> b = T $ b - a +------------------------------------------------------------------------------- +-- Extended division + +extendedDivide :: ( Field d, Field ( 𝕀 d ), Ord d ) => 𝕀 d -> 𝕀 d -> [ 𝕀 d ] +extendedDivide x y = map ( x * ) ( extendedRecip y ) +{-# SPECIALISE extendedDivide :: 𝕀 Double -> 𝕀 Double -> [ 𝕀 Double ] #-} + +extendedRecip :: ( Field d, Field ( 𝕀 d ), Ord d ) => 𝕀 d -> [ 𝕀 d ] +extendedRecip x@( 𝕀 lo hi ) + | lo == fromInteger 0 && hi == fromInteger 0 + = [ 𝕀 negInf posInf ] + | lo >= fromInteger 0 || hi <= fromInteger 0 + = [ recip x ] + | otherwise + = [ 𝕀 negInf ( recip lo ), 𝕀 ( recip hi ) posInf ] + where + negInf = fromInteger (-1) / fromInteger 0 + posInf = fromInteger 1 / fromInteger 0 +{-# SPECIALISE extendedRecip :: 𝕀 Double -> [ 𝕀 Double ] #-} + ------------------------------------------------------------------------------- -- Lattices. diff --git a/brush-strokes/src/Math/Interval/Internal.hs b/brush-strokes/src/Math/Interval/Internal.hs index aa649fe..b170745 100644 --- a/brush-strokes/src/Math/Interval/Internal.hs +++ b/brush-strokes/src/Math/Interval/Internal.hs @@ -109,7 +109,11 @@ instance Prelude.Fractional ( 𝕀 Double ) where in 𝕀 q q recip (𝕀 lo hi) -- #ifdef ASSERTS - | lo >= 0 || hi <= 0 + | lo == 0 + = 𝕀 ( fst $ divI 1 hi ) ( 1 / 0 ) + | hi == 0 + = 𝕀 ( -1 / 0 ) ( snd $ divI 1 lo ) + | lo > 0 || hi < 0 -- #endif = 𝕀 ( fst $ divI 1 hi ) ( snd $ divI 1 lo ) -- #ifdef ASSERTS @@ -118,7 +122,11 @@ instance Prelude.Fractional ( 𝕀 Double ) where -- #endif 𝕀 x_lo x_hi / 𝕀 y_lo y_hi -- #ifdef ASSERTS - | y_lo >= 0 || y_hi <= 0 + | y_lo == 0 + = 𝕀 ( fst $ divI x_lo y_hi ) ( 1 / 0 ) + | y_hi == 0 + = 𝕀 ( -1 / 0 ) ( snd $ divI x_hi y_lo ) + | y_lo > 0 || y_hi < 0 -- #endif = 𝕀 ( fst $ divI x_lo y_hi ) ( snd $ divI x_hi y_lo ) -- #ifdef ASSERTS