Fixes and restructuring

This commit is contained in:
sheaf 2024-02-28 17:20:34 +01:00
parent 26cfdada8f
commit 2289468a84
35 changed files with 1021 additions and 557 deletions

View file

@ -187,49 +187,6 @@ library metabrushes
, bytestring
>= 0.10.10.0 && < 0.12
executable cusps
import:
common
hs-source-dirs:
src/cusps
default-language:
Haskell2010
main-is:
Main.hs
other-modules:
Math.Interval.Abstract
executable convert-metafont
import:
common, extras
hs-source-dirs:
src/convert
default-language:
Haskell2010
main-is:
Main.hs
other-modules:
MetaBrush.MetaFont.Convert
build-depends:
metabrushes,
diagrams-contrib,
diagrams-lib,
linear,
parsec
executable MetaBrush
import:

View file

@ -1,370 +0,0 @@
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE UndecidableInstances #-}
module Main
( main
-- Testing
, TestCase(..)
, testCases
, testCaseStrokeFunctions
, eval
, mkVal, mkBox
, potentialCusp
, dEdsdcdt
)
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
( Seq )
import qualified Data.Sequence as Seq
( index )
import Data.Tree
( foldTree )
-- brush-strokes
import Calligraphy.Brushes
import Debug.Utils
( logToFile )
import Math.Algebra.Dual
import Math.Bezier.Spline
import Math.Bezier.Stroke
import Math.Bezier.Stroke.EnvelopeEquation
import Math.Differentiable
import Math.Interval
import Math.Linear
import Math.Module
import Math.Ring
( Transcendental )
--------------------------------------------------------------------------------
main :: IO ()
main = for_ testCases $ \ testCase@( TestCase { testName, testAlgorithmParams } ) -> do
let ( _, testStrokeFnI ) = testCaseStrokeFunctions testCase
( newtTrees, ( dunno, sols ) ) = computeCusps testAlgorithmParams testStrokeFnI
showedTrees = map ( uncurry showIntervalNewtonTree ) newtTrees
putStrLn $ unlines $
[ ""
, "Test case '" ++ testName ++ "':" ] ++
map ( " " ++ )
[ " #sols: " ++ show (length sols)
, "#dunno: " ++ show (length dunno)
, "#trees: " ++ show @Int (sum @_ @Int $ map (foldTree ( \ _ bs -> 1 + sum bs )) showedTrees)
, " dunno: " ++ show dunno
, " sols: " ++ show sols
]
--logFileMVar <- newMVar "logs/trickyCusp.log"
--logToFile logFileMVar (unlines logLines)
-- `seq` return ()
testCases :: [ TestCase ]
testCases = [ ellipse , trickyCusp2 ]
--------------------------------------------------------------------------------
data TestCase =
forall nbParams. ParamsCt nbParams =>
TestCase
{ testName :: !String
, testBrush :: !( Brush nbParams )
, testStroke :: !( Point nbParams, Curve Open () ( Point nbParams ))
, testAlgorithmParams :: !CuspAlgorithmParams
}
testCaseStrokeFunctions
:: TestCase
-> ( 1 -> Seq ( 1 -> StrokeDatum 2 () )
, 𝕀 1 -> Seq ( 𝕀 1 -> StrokeDatum 3 𝕀 ) )
testCaseStrokeFunctions ( TestCase { testStroke = ( sp0, crv ), testBrush } ) =
getStrokeFunctions testBrush sp0 crv
-- Utilities to use in GHCi to help debugging.
eval
:: ( I i ( 1 ) -> Seq ( I i ( 1 ) -> StrokeDatum k i ) )
-> ( I i ( 1 ), Int, I i ( 1 ) )
-> StrokeDatum k i
eval f ( t, i, s ) = ( f t `Seq.index` i ) s
mkVal :: Double -> Int -> Double -> ( 1, Int, 1 )
mkVal t i s = ( 1 t, i, 1 s )
mkBox :: ( Double, Double ) -> Int -> ( Double, Double ) -> Box
mkBox ( t_min, t_max ) i ( s_min, s_max ) =
( 𝕀 ( 1 t_min ) ( 1 t_max ) , i, 𝕀 ( 1 s_min ) ( 1 s_max ) )
potentialCusp :: StrokeDatum 3 𝕀 -> Bool
potentialCusp
( StrokeDatum
{ ee = D22 { _D22_v = 𝕀 ( 1 ee_min ) ( 1 ee_max ) }
, 𝛿E𝛿sdcdt = D12 { _D12_v = T ( 𝕀 ( 2 vx_min vy_min ) ( 2 vx_max vy_max ) )}
}
) = ee_min <= 0 && ee_max >= 0
&& vx_min <= 0 && vx_max >= 0
&& vy_min <= 0 && vy_max >= 0
dEdsdcdt :: StrokeDatum k i -> D ( k - 2 ) ( I i ( 2 ) ) ( T ( I i ( 2 ) ) )
dEdsdcdt ( StrokeDatum { 𝛿E𝛿sdcdt = v } ) = v
{-
let (f, fI) = testCaseStrokeFunctions trickyCusp2
take 10 $ Data.List.sortOn ( \ ( _, 1 e, v) -> abs e + norm v ) [ let { v = mkVal x 3 y; d = eval f v } in ( v, _D12_v $ ee d, _D0_v $ dEdsdcdt d ) | x <- [0.57,0.5701 .. 0.58], y <- [0.29,0.291..0.3] ]
> ((1 0.5798800000000057,3,1 0.267980000000008),1 -2.8596965543670194e-4,V2 7.79559474412963e-2 2.0389671921293484e-2)
potentialCusp $ eval fI $ mkBox (0.5798, 0.5799) 3 (0.26798, 0.26799)
> True
let nbPotentialSols b = let ( _newtTrees, ( dunno, sols ) ) = intervalNewtonGSFrom NoPreconditioning 1e-7 fI b in length dunno + length sols
nbPotentialSols $ mkBox (0.5798, 0.5799) 3 (0.26798, 0.26799)
1
nbPotentialSols $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799)
0
let showTrees b = map ( uncurry showIntervalNewtonTree ) $ fst $ intervalNewtonGSFrom NoPreconditioning 1e-7 fI b
putStrLn $ unlines $ map Data.Tree.View.showTree $ showTrees $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799)
([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
ellipse =
TestCase
{ testName = "ellipse"
, testBrush = ellipseBrush
, testStroke = ( p0, LineTo ( NextPoint p1 ) () )
, testAlgorithmParams =
CuspAlgorithmParams
{ preconditioning = NoPreconditioning
, maxWidth = 1e-7
}
}
where
mkPt x y w h phi =
Point
{ pointCoords = 2 x y
, pointParams = Params $ 3 w h phi
}
p0 = mkPt 0 0 10 25 0
p1 = mkPt 100 0 15 40 pi
trickyCusp2 :: TestCase
trickyCusp2 =
TestCase
{ testName = "trickyCusp2"
, testBrush = circleBrush
, testStroke = ( p0, Bezier3To p1 p2 ( NextPoint p3 ) () )
, testAlgorithmParams =
CuspAlgorithmParams
{ preconditioning = NoPreconditioning
, maxWidth = 1e-7
}
}
where
mkPt x y =
Point
{ pointCoords = 2 x y
, pointParams = Params $ 1 5.0
}
p0 = mkPt 5e+1 -5e+1
p1 = mkPt -7.72994362904069e+1 -3.124468786098509e+1
p2 = mkPt -5.1505430313958364e+1 -3.9826386521527986e+1
p3 = mkPt -5e+1 -5e+1
--------------------------------------------------------------------------------
type ParamsCt nbParams
= ( Show ( nbParams )
, HasChainRule Double 2 ( nbParams )
, HasChainRule ( 𝕀 Double ) 3 ( 𝕀 ( nbParams ) )
, Applicative ( D 2 ( nbParams ) )
, Applicative ( D 3 ( nbParams ) )
, Traversable ( D 2 ( nbParams ) )
, Traversable ( D 3 ( nbParams ) )
, Representable Double ( nbParams )
, Module Double ( T ( nbParams ) )
, Module ( 𝕀 Double ) ( T ( 𝕀 ( nbParams ) ) )
, Module ( D 2 ( nbParams ) Double ) ( D 2 ( nbParams ) ( 2 ) )
, Module ( D 3 ( nbParams ) ( 𝕀 Double ) ) ( D 3 ( nbParams ) ( 𝕀 ( 2 ) ) )
, Transcendental ( D 2 ( nbParams ) Double )
, Transcendental ( D 3 ( nbParams ) ( 𝕀 Double ) )
)
newtype Params nbParams = Params { getParams :: ( nbParams ) }
deriving newtype instance Show ( nbParams ) => Show ( Params nbParams )
data Point nbParams =
Point
{ pointCoords :: !( 2 )
, pointParams :: !( Params nbParams ) }
deriving stock Generic
deriving stock instance Show ( nbParams ) => Show ( Point nbParams )
data CuspAlgorithmParams =
CuspAlgorithmParams
{ preconditioning :: !Preconditioner
, maxWidth :: !Double
}
deriving stock Show
type Brush nbParams
= forall {t} k (i :: t)
. DiffInterp k i ( nbParams )
=> Proxy# i
-> ( forall a. a -> I i a )
-> C k ( I i ( nbParams ) )
( Spline Closed () ( I i ( 2 ) ) )
getStrokeFunctions
:: forall nbParams
. ParamsCt nbParams
=> Brush nbParams
-- ^ brush shape
-> Point nbParams
-- ^ start point
-> Curve Open () ( Point nbParams )
-- ^ curve points
-> ( 1 -> Seq ( 1 -> StrokeDatum 2 () )
, 𝕀 1 -> Seq ( 𝕀 1 -> StrokeDatum 3 𝕀 ) )
getStrokeFunctions brush sp0 crv =
let
usedParams :: C 2 ( 1 ) ( nbParams )
path :: C 2 ( 1 ) ( 2 )
( path, usedParams ) =
pathAndUsedParams @2 @() coerce id ( getParams . pointParams )
sp0 crv
usedParamsI :: C 3 ( 𝕀 1 ) ( 𝕀 nbParams )
pathI :: C 3 ( 𝕀 1 ) ( 𝕀 2 )
( pathI, usedParamsI ) =
pathAndUsedParams @3 @𝕀 coerce singleton ( getParams . pointParams )
sp0 crv
in ( brushStrokeData @2 @( nbParams ) coerce coerce
path usedParams $
brush @2 @() proxy# id
, brushStrokeData @3 @( nbParams ) coerce coerce
pathI usedParamsI $
brush @3 @𝕀 proxy# singleton )
{-# INLINEABLE getStrokeFunctions #-}
computeCusps
:: CuspAlgorithmParams
-> ( 𝕀 1 -> Seq ( 𝕀 1 -> StrokeDatum 3 𝕀 ) )
-> ( [ ( Box, IntervalNewtonTree Box ) ], ( [ Box ], [ Box ] ) )
computeCusps params =
intervalNewtonGS ( preconditioning params ) ( maxWidth params )

View file

@ -83,13 +83,23 @@ common common
base
>= 4.19
common extra
build-depends:
acts
^>= 0.3.1.0
, generic-lens
>= 2.2 && < 2.3
, groups
^>= 0.5.3
library
import:
common
common, extra
hs-source-dirs:
src
src/lib
default-language:
Haskell2010
@ -129,9 +139,6 @@ library
build-depends:
template-haskell
>= 2.18 && < 2.22
, acts
^>= 0.3.1.0
, bifunctors
>= 5.5.4 && < 5.7
, code-page
@ -144,10 +151,6 @@ library
^>= 3.3.7.0
, filepath
>= 1.4 && < 1.6
, generic-lens
>= 2.2 && < 2.3
, groups
^>= 0.5.3
, groups-generic
^>= 0.3.1.0
, parallel
@ -161,13 +164,58 @@ library
, transformers
>= 0.5.6.2 && < 0.7
--executable convert-metafont
--
-- import:
-- common
--
-- hs-source-dirs:
-- src/metafont
--
-- default-language:
-- Haskell2010
--
-- main-is:
-- Main.hs
--
-- other-modules:
-- Calligraphy.MetaFont.Convert
--
-- build-depends:
-- diagrams-contrib,
-- diagrams-lib,
-- linear,
-- parsec
executable inspect
import:
common, extra
hs-source-dirs:
src/cusps/inspect
default-language:
Haskell2010
main-is:
Main.hs
other-modules:
Math.Interval.Abstract
build-depends:
brush-strokes,
data-reify
^>= 0.6.3
benchmark cusps
import:
common
hs-source-dirs:
bench
src/cusps/bench
main-is:
Main.hs

View file

@ -0,0 +1,616 @@
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE UndecidableInstances #-}
module Main
( main
-- Testing
, TestCase(..)
, testCases
, BrushStroke(..)
, brushStrokeFunctions
, eval
, mkVal, mkBox
, potentialCusp
, dEdsdcdt
)
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
( Seq )
import qualified Data.Sequence as Seq
( index )
import Data.Tree
( foldTree )
-- tree-view
import Data.Tree.View
( showTree )
-- brush-strokes
import Calligraphy.Brushes
import Debug.Utils
( logToFile )
import Math.Algebra.Dual
import Math.Bezier.Spline
import Math.Bezier.Stroke
import Math.Bezier.Stroke.EnvelopeEquation
import Math.Differentiable
import Math.Interval
import Math.Linear
import Math.Module
import Math.Ring
( Transcendental )
--------------------------------------------------------------------------------
main :: IO ()
main = do
for_ testCases $ \ testCase@( TestCase { testName, testBrushStroke, testAlgorithmParams, testStartBoxes } ) -> do
let ( _, testStrokeFnI ) = brushStrokeFunctions testBrushStroke
( newtTrees, ( dunno, sols ) ) = computeCusps testAlgorithmParams testStrokeFnI testStartBoxes
showedTrees = map ( uncurry showIntervalNewtonTree ) newtTrees
testHeader =
[ "", "Test case '" ++ testName ++ "':" ]
putStrLn $ unlines $
testHeader ++
map ( " " ++ )
[ " #sols: " ++ show (length sols)
, "#dunno: " ++ show (length dunno)
, "#trees: " ++ show @Int (sum @_ @Int $ map (foldTree ( \ _ bs -> 1 + sum bs )) showedTrees)
, " dunno: " ++ show dunno
, " sols: " ++ show sols
]
-- logFileMVar <- newMVar "logs/fnData.log"
-- logToFile logFileMVar (unlines logLines)
-- `seq` return ()
testCases :: [ TestCase ]
testCases = benchCases
-- [ --trickyCusp2
-- ellipse "full" (0,1) pi $ defaultStartBoxes [ 2 ]
-- ]
-- ++
{-
[ ellipse ( "(k1, k2) = " ++ show (k1, k2) ) (k1, k2) pi $ defaultStartBoxes [ 2 ]
| (k1, k2) <-
[(0.5,0.6), (0.55, 0.56)]
] ++
[ ellipse ( "'(k1, k2) = " ++ show (k1, k2) ) (0,1) pi [ mkBox (k1 + zero, k2 + zero) 2 (zero, one) ]
| (k1, k2) <-
[(0.5,0.6), (0.55, 0.56)]
]
-}
benchCases :: [ TestCase ]
benchCases = [ ellipseTestCase "full" ( 0, 1 ) pi $ defaultStartBoxes [ 0 .. 3 ] ]
--------------------------------------------------------------------------------
data BrushStroke =
forall nbParams. ParamsCt nbParams =>
BrushStroke
{ brush :: !( Brush nbParams )
, stroke :: !( Point nbParams, Curve Open () ( Point nbParams ) )
}
data TestCase =
TestCase
{ testName :: String
, testBrushStroke :: BrushStroke
, testAlgorithmParams :: CuspAlgorithmParams
, testStartBoxes :: [ Box ]
}
brushStrokeFunctions
:: BrushStroke
-> ( 1 -> Seq ( 1 -> StrokeDatum 2 () )
, 𝕀 1 -> Seq ( 𝕀 1 -> StrokeDatum 3 𝕀 ) )
brushStrokeFunctions ( BrushStroke { stroke = ( sp0, crv ), brush } ) =
getStrokeFunctions brush sp0 crv
-- Utilities to use in GHCi to help debugging.
eval
:: ( I i ( 1 ) -> Seq ( I i ( 1 ) -> StrokeDatum k i ) )
-> ( I i ( 1 ), Int, I i ( 1 ) )
-> StrokeDatum k i
eval f ( t, i, s ) = ( f t `Seq.index` i ) s
mkVal :: Double -> Int -> Double -> ( 1, Int, 1 )
mkVal t i s = ( 1 t, i, 1 s )
mkBox :: ( Double, Double ) -> Int -> ( Double, Double ) -> Box
mkBox ( t_min, t_max ) i ( s_min, s_max ) =
( 𝕀 ( 1 t_min ) ( 1 t_max ) , i, 𝕀 ( 1 s_min ) ( 1 s_max ) )
zero, one :: Double
zero = 1e-6
one = 1 - zero
{-# INLINE zero #-}
{-# INLINE one #-}
potentialCusp :: StrokeDatum 3 𝕀 -> Bool
potentialCusp
( StrokeDatum
{ ee = D22 { _D22_v = 𝕀 ( 1 ee_min ) ( 1 ee_max ) }
, 𝛿E𝛿sdcdt = D12 { _D12_v = T ( 𝕀 ( 2 vx_min vy_min ) ( 2 vx_max vy_max ) )}
}
) = ee_min <= 0 && ee_max >= 0
&& vx_min <= 0 && vx_max >= 0
&& vy_min <= 0 && vy_max >= 0
dEdsdcdt :: StrokeDatum k i -> D ( k - 2 ) ( I i ( 2 ) ) ( T ( I i ( 2 ) ) )
dEdsdcdt ( StrokeDatum { 𝛿E𝛿sdcdt = v } ) = v
{-
let (f, fI) = testCaseStrokeFunctions trickyCusp2
take 10 $ Data.List.sortOn ( \ ( _, 1 e, v) -> abs e + norm v ) [ let { v = mkVal x 3 y; d = eval f v } in ( v, _D12_v $ ee d, _D0_v $ dEdsdcdt d ) | x <- [0.57,0.5701 .. 0.58], y <- [0.29,0.291..0.3] ]
> ((1 0.5798800000000057,3,1 0.267980000000008),1 -2.8596965543670194e-4,V2 7.79559474412963e-2 2.0389671921293484e-2)
potentialCusp $ eval fI $ mkBox (0.5798, 0.5799) 3 (0.26798, 0.26799)
> True
let nbPotentialSols b = let ( _newtTrees, ( dunno, sols ) ) = intervalNewtonGSFrom NoPreconditioning 1e-7 fI b in length dunno + length sols
nbPotentialSols $ mkBox (0.5798, 0.5799) 3 (0.26798, 0.26799)
1
nbPotentialSols $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799)
0
let showTrees b = map ( uncurry showIntervalNewtonTree ) $ fst $ intervalNewtonGSFrom NoPreconditioning 1e-7 fI b
putStrLn $ unlines $ map Data.Tree.View.showTree $ showTrees $ mkBox (0.5798, 0.675) 3 (0.26798, 0.26799)
([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
-----------------------
fI k rot = snd $ testCaseStrokeFunctions $ ellipse k rot
d k rot = width $ _D22_v $ ee $ eval (fI k rot) $ mkBox (1e-6, 1-1e-6) 3 (1e-6, 1-1e-6)
-----------------------
(f, fI) = testCaseStrokeFunctions $ ellipse 1 pi
t = 𝕀 (1 0.001) (1 0.099)
s = 𝕀 (1 0.001) (1 0.999)
t_mid = 0.5 * ( 0.001 + 0.099 )
s_mid = 0.5 * ( 0.001 + 0.999 )
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
> t = [0.001, 0.099]
> t' = [-0.049, 0.049]
(f, fI) = testCaseStrokeFunctions $ ellipse 0.1 pi
t = 𝕀 (1 0.001) (1 0.999)
s = 𝕀 (1 0.001) (1 0.999)
t_mid = 0.5 * ( 0.001 + 0.999 )
s_mid = 0.5 * ( 0.001 + 0.999 )
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
> t = [0.001, 0.999]
> t' = [-0.499, 0.499]
-}
width :: 𝕀 1 -> Double
width (𝕀 (1 lo) (1 hi)) = hi - lo
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 =
[ "E, dE/ds * dc/dt"
, "{" ++
(intercalate ","
[ "{" ++ showD t ++ "," ++ showD s ++ ",{" ++ intercalate "," vals ++ "}}"
| t <- [ 0.5484, 0.5484 + 0.00001 .. 0.5488 ]
, s <- [ 0.5479, 0.5479 + 0.00001 .. 0.5483 ]
, let StrokeDatum
{ ee = D22 ee _ _ _ _ _
, 𝛿E𝛿sdcdt = D12 (T f) (T (T f_t)) (T (T f_s))
} = (curvesI (singleton (1 t)) `Seq.index` i) (singleton (1 s))
2 vx vy = midPoint2 f
--2 vx_t vy_t = midPoint2 f_t
--2 vx_s vy_s = midPoint2 f_s
vals = [ showD ( midPoint ee )
, "{" ++ showD vx ++ "," ++ showD vy ++ "}"
-- , "{" ++ showD vx_t ++ "," ++ showD vy_t ++ "}"
-- , "{" ++ showD vx_s ++ "," ++ showD vy_s ++ "}"
]
]
) ++ "}"
]
where
i = 2
( curves, curvesI ) = brushStrokeFunctions $ ellipseBrushStroke ( 0, 1 ) pi
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 ) )
-- t = 0.5486102, s = 0.5480951
bloo =
[ ( e * e + vx * vx + vy * vy, ( t, s ) )
| t <- [ 0.548609, 0.548609 + 0.0000001 .. 0.54862 ]
, s <- [ 0.548094, 0.548094 + 0.0000001 .. 0.548096 ]
, let StrokeDatum
{ ee = D22 ee _ _ _ _ _
, 𝛿E𝛿sdcdt = D12 (T f) _ _
} = (curvesI (singleton (1 t)) `Seq.index` i) (singleton (1 s))
e = midPoint ee
2 vx vy = midPoint2 f
vals = [ showD ( midPoint ee )
, "{" ++ showD vx ++ "," ++ showD vy ++ "}"
]
]
where
i = 2
( curves, curvesI ) = brushStrokeFunctions $ ellipseBrushStroke ( 0, 1 ) pi
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 ) )
{-
(f, fI) = brushStrokeFunctions $ ellipseBrushStroke ( 0, 1 ) pi
_D22_v $ ee $ eval fI $ mkBox (0.5, 0.6) 2 (0,1)
> [1 -9531.427315889887, 1 10135.074304695485]
minimum $ map inf $ [ _D22_v $ ee $ eval fI $ mkBox (t, t + 0.01) 2 (0,1) | t <- [ 0.5, 0.51 .. 0.59 ] ]
> 1 -5718.905635365308
maximum $ map sup $ [ _D22_v $ ee $ eval fI $ mkBox (t, t + 0.01) 2 (0,1) | t <- [ 0.5, 0.51 .. 0.59 ] ]
> 1 5099.008191092755
minimum $ map inf $ [ _D22_v $ ee $ eval fI $ mkBox (t, t) 2 (s, s) | t <- [ 0.5, 0.501 .. 0.6 ], s <- [ zero, zero + 0.001 .. one ] ]
> 1 -675.9595496147458
maximum $ map sup $ [ _D22_v $ ee $ eval fI $ mkBox (t, t) 2 (s, s) | t <- [ 0.5, 0.501 .. 0.6 ], s <- [ zero, zero + 0.001 .. one ] ]
> 1 2401.9644509525997
_D12_v $ dEdsdcdt $ eval fI $ mkBox (0.5, 0.6) 2 (0,1)
> T[2 -1.7300637136531524e7 -1.262824151868635e7, 2 1.632868898735965e7 1.1869759856947478e7]
minimum [ _x $ inf $ unT $ _D12_v $ dEdsdcdt $ eval fI $ mkBox (t, t + 0.01) 2 (0,1) | t <- [ 0.5, 0.51 .. 0.59 ] ]
-5606615.948203902
maximum [ _x $ sup $ unT $ _D12_v $ dEdsdcdt $ eval fI $ mkBox (t, t + 0.01) 2 (0,1) | t <- [ 0.5, 0.51 .. 0.59 ] ]
4340858.832347277
minimum [ _x $ inf $ unT $ _D12_v $ dEdsdcdt $ eval fI $ mkBox (t, t) 2 (s, s) | t <- [ 0.5, 0.501 .. 0.6 ], s <- [ zero, zero + 0.001 .. one ] ]
-1785730.2396688666
maximum [ _x $ sup $ unT $ _D12_v $ dEdsdcdt $ eval fI $ mkBox (t, t) 2 (s, s) | t <- [ 0.5, 0.501 .. 0.6 ], s <- [ zero, zero + 0.001 .. one ] ]
974842.6547409865
maximum [ _y $ sup $ unT $ _D12_v $ dEdsdcdt $ eval fI $ mkBox (t, t) 2 (s, s) | t <- [ 0.5, 0.501 .. 0.6 ], s <- [ zero, zero + 0.001 .. one ] ]
845211.4833711373
let showTrees b = map ( uncurry showIntervalNewtonTree ) $ fst $ intervalNewtonGSFrom NoPreconditioning 1e-7 fI b
putStrLn $ unlines $ map Data.Tree.View.showTree $ showTrees $ mkBox (0.5486101933245248, 0.5486102071493622) 2 (0.548095036738487, 0.5480952)
([1 0.5486101960904595, 1 0.5486102071493623],2,[1 0.5480950771755867, 1 0.5480952000000001])
-}
_x ( 2 x _ ) = x
_y ( 2 _ y ) = y
showD :: Double -> String
showD float = showFFloat (Just 6) float ""
--------------------------------------------------------------------------------
ellipseTestCase :: String -> ( Double, Double ) -> Double -> [ Box ] -> TestCase
ellipseTestCase str k0k1 rot startBoxes =
TestCase
{ testName = "ellipse (" ++ str ++ ")"
, testBrushStroke = ellipseBrushStroke k0k1 rot
, testAlgorithmParams =
CuspAlgorithmParams
{ preconditioning = InverseMidJacobian
, maxWidth = 1e-7
}
, testStartBoxes = startBoxes
}
ellipseBrushStroke :: ( Double, Double ) -> Double -> BrushStroke
ellipseBrushStroke ( k0, k1 ) rot =
BrushStroke
{ brush = ellipseBrush
, stroke = ( p0, LineTo ( NextPoint p1 ) () ) }
where
mkPt x y w h phi =
Point
{ pointCoords = 2 x y
, pointParams = Params $ 3 w h phi
}
l k = lerp @( T Double ) k
p k = mkPt ( l k 0 100 ) 0 ( l k 10 15 ) ( l k 25 40 ) ( l k 0 rot )
p0 = p k0
p1 = p k1
trickyCusp2TestCase :: TestCase
trickyCusp2TestCase =
TestCase
{ testName = "trickyCusp2"
, testBrushStroke = trickyCusp2BrushStroke
, testAlgorithmParams =
CuspAlgorithmParams
{ preconditioning = InverseMidJacobian
, maxWidth = 1e-7
}
, testStartBoxes = defaultStartBoxes [ 0 .. 3 ]
}
trickyCusp2BrushStroke :: BrushStroke
trickyCusp2BrushStroke =
BrushStroke
{ brush = circleBrush
, stroke = ( p0, Bezier3To p1 p2 ( NextPoint p3 ) () )
}
where
mkPt x y =
Point
{ pointCoords = 2 x y
, pointParams = Params $ 1 5.0
}
p0 = mkPt 5e+1 -5e+1
p1 = mkPt -7.72994362904069e+1 -3.124468786098509e+1
p2 = mkPt -5.1505430313958364e+1 -3.9826386521527986e+1
p3 = mkPt -5e+1 -5e+1
--------------------------------------------------------------------------------
type ParamsCt nbParams
= ( Show ( nbParams )
, HasChainRule Double 2 ( nbParams )
, HasChainRule ( 𝕀 Double ) 3 ( 𝕀 ( nbParams ) )
, Applicative ( D 2 ( nbParams ) )
, Applicative ( D 3 ( nbParams ) )
, Traversable ( D 2 ( nbParams ) )
, Traversable ( D 3 ( nbParams ) )
, Representable Double ( nbParams )
, Module Double ( T ( nbParams ) )
, Module ( 𝕀 Double ) ( T ( 𝕀 ( nbParams ) ) )
, Module ( D 2 ( nbParams ) Double ) ( D 2 ( nbParams ) ( 2 ) )
, Module ( D 3 ( nbParams ) ( 𝕀 Double ) ) ( D 3 ( nbParams ) ( 𝕀 ( 2 ) ) )
, Transcendental ( D 2 ( nbParams ) Double )
, Transcendental ( D 3 ( nbParams ) ( 𝕀 Double ) )
)
newtype Params nbParams = Params { getParams :: ( nbParams ) }
deriving newtype instance Show ( nbParams ) => Show ( Params nbParams )
data Point nbParams =
Point
{ pointCoords :: !( 2 )
, pointParams :: !( Params nbParams ) }
deriving stock Generic
deriving stock instance Show ( nbParams ) => Show ( Point nbParams )
data CuspAlgorithmParams =
CuspAlgorithmParams
{ preconditioning :: !Preconditioner
, maxWidth :: !Double
}
deriving stock Show
type Brush nbParams
= forall {t} k (i :: t)
. DiffInterp k i ( nbParams )
=> Proxy# i
-> ( forall a. a -> I i a )
-> C k ( I i ( nbParams ) )
( Spline Closed () ( I i ( 2 ) ) )
getStrokeFunctions
:: forall nbParams
. ParamsCt nbParams
=> Brush nbParams
-- ^ brush shape
-> Point nbParams
-- ^ start point
-> Curve Open () ( Point nbParams )
-- ^ curve points
-> ( 1 -> Seq ( 1 -> StrokeDatum 2 () )
, 𝕀 1 -> Seq ( 𝕀 1 -> StrokeDatum 3 𝕀 ) )
getStrokeFunctions brush sp0 crv =
let
usedParams :: C 2 ( 1 ) ( nbParams )
path :: C 2 ( 1 ) ( 2 )
( path, usedParams ) =
pathAndUsedParams @2 @() coerce id ( getParams . pointParams )
sp0 crv
usedParamsI :: C 3 ( 𝕀 1 ) ( 𝕀 nbParams )
pathI :: C 3 ( 𝕀 1 ) ( 𝕀 2 )
( pathI, usedParamsI ) =
pathAndUsedParams @3 @𝕀 coerce singleton ( getParams . pointParams )
sp0 crv
in ( brushStrokeData @2 @( nbParams ) coerce coerce
path usedParams $
brush @2 @() proxy# id
, brushStrokeData @3 @( nbParams ) coerce coerce
pathI usedParamsI $
brush @3 @𝕀 proxy# singleton )
{-# INLINEABLE getStrokeFunctions #-}
computeCusps
:: CuspAlgorithmParams
-> ( 𝕀 1 -> Seq ( 𝕀 1 -> StrokeDatum 3 𝕀 ) )
-> [ Box ]
-> ( [ ( Box, IntervalNewtonTree Box ) ], ( [ Box ], [ Box ] ) )
computeCusps params eqs startBoxes =
foldMap
( intervalNewtonGSFrom ( preconditioning params ) ( maxWidth params ) eqs )
startBoxes
defaultStartBoxes :: [ Int ] -> [ Box ]
defaultStartBoxes is =
[ mkBox (zero, one) i (zero, one) | i <- is ]
getR1 (1 u) = u
{-
(f, fI) = brushStrokeFunctions $ ellipseBrushStroke (0,1) pi
nbPotentialSols box = let ( _newtTrees, ( dunno, sols ) ) = intervalNewtonGSFrom NoPreconditioning 1e-7 fI box in length dunno + length sols
showTrees box = map ( uncurry showIntervalNewtonTree ) $ fst $ intervalNewtonGSFrom NoPreconditioning 1e-7 fI box
(t, i, s) = mkBox (0.548610200176363, 0.5486102071493623) 2 (0.5480950215354709, 0.5480952)
putStrLn $ unlines $ map Data.Tree.View.showTree $ showTrees (t,i,s)
t_mid = 0.5 * ( getR1 ( inf t ) + getR1 ( sup t ) )
s_mid = 0.5 * ( getR1 ( inf s ) + getR1 ( sup s ) )
D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) = dEdsdcdt $ eval fI (t, i, 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.548610200176363, 1 0.5486102071493624]
s2
> [1 0.5480950911334656, 1 0.5480952000000001]
mkBox (0.548610200176363, 0.5486102071493624) i (0.5480950911334656, 0.5480952000000001)
t inf (no change)
t sup (no change)
s_inf:
0.5480950215354709
0.5480950911334656
s_sup (no change)
ghci> potentialCusp $ eval fI $ mkBox (0.548610200176363, 0.5486102071493623) 2 (0.54809502, 0.5480952)
True
ghci> potentialCusp $ eval fI $ mkBox (0.548610200176363, 0.5486102071493623) 2 (0.54809503, 0.5480952)
False
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
-}

View file

@ -32,7 +32,7 @@ import qualified Data.Sequence as Seq
import Data.Generics.Product.Typed
( HasType )
-- MetaBrush
-- brush-strokes
import Math.Algebra.Dual
import qualified Math.Bezier.Quadratic as Quadratic
import qualified Math.Bezier.Cubic as Cubic
@ -58,6 +58,7 @@ import Math.Ring
, ifThenElse
)
-- brush-strokes:inspect-cusps
import Math.Interval.Abstract
--------------------------------------------------------------------------------

View file

@ -43,7 +43,7 @@ import qualified Data.Map.Strict as Map
import Data.Group
( Group(..) )
-- MetaBrush
-- brush-strokes
import Math.Algebra.Dual
import qualified Math.Bezier.Quadratic as Quadratic
import qualified Math.Bezier.Cubic as Cubic

View file

@ -66,3 +66,5 @@ logToFile logFileMVar logContents =
FilePath.takeDirectory logFile
appendFile logFile logContentsWithHeader
return logFile
{-# NOINLINE logToFile #-}

View file

@ -60,8 +60,8 @@ newtype D𝔸0 v = D0 { _D0_v :: v }
-- | \( \mathbb{Z}[x]/(x)^2 \).
data D1𝔸1 v =
D11 { _D11_v :: !v
, _D11_dx :: !( T v )
D11 { _D11_v :: v
, _D11_dx :: ( T v )
}
deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
deriving anyclass NFData
@ -70,9 +70,9 @@ data D1𝔸1 v =
-- | \( \mathbb{Z}[x]/(x)^3 \).
data D2𝔸1 v =
D21 { _D21_v :: !v
, _D21_dx :: !( T v )
, _D21_dxdx :: !( T v )
D21 { _D21_v :: v
, _D21_dx :: ( T v )
, _D21_dxdx :: ( T v )
}
deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
deriving anyclass NFData
@ -81,10 +81,10 @@ data D2𝔸1 v =
-- | \( \mathbb{Z}[x]/(x)^4 \).
data D3𝔸1 v =
D31 { _D31_v :: !v
, _D31_dx :: !( T v )
, _D31_dxdx :: !( T v )
, _D31_dxdxdx :: !( T v )
D31 { _D31_v :: v
, _D31_dx :: ( T v )
, _D31_dxdx :: ( T v )
, _D31_dxdxdx :: ( T v )
}
deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
deriving anyclass NFData
@ -93,8 +93,8 @@ data D3𝔸1 v =
-- | \( \mathbb{Z}[x, y]/(x, y)^2 \).
data D1𝔸2 v =
D12 { _D12_v :: !v
, _D12_dx, _D12_dy :: !( T v )
D12 { _D12_v :: v
, _D12_dx, _D12_dy :: ( T v )
}
deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
deriving anyclass NFData
@ -103,9 +103,9 @@ data D1𝔸2 v =
-- | \( \mathbb{Z}[x, y]/(x, y)^3 \).
data D2𝔸2 v =
D22 { _D22_v :: !v
, _D22_dx, _D22_dy :: !( T v )
, _D22_dxdx, _D22_dxdy, _D22_dydy :: !( T v )
D22 { _D22_v :: v
, _D22_dx, _D22_dy :: ( T v )
, _D22_dxdx, _D22_dxdy, _D22_dydy :: ( T v )
}
deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
deriving anyclass NFData
@ -114,10 +114,10 @@ data D2𝔸2 v =
-- | \( \mathbb{Z}[x, y]/(x, y)^4 \).
data D3𝔸2 v =
D32 { _D32_v :: !v
, _D32_dx, _D32_dy :: !( T v )
, _D32_dxdx, _D32_dxdy, _D32_dydy :: !( T v )
, _D32_dxdxdx, _D32_dxdxdy, _D32_dxdydy, _D32_dydydy :: !( T v )
D32 { _D32_v :: v
, _D32_dx, _D32_dy :: ( T v )
, _D32_dxdx, _D32_dxdy, _D32_dydy :: ( T v )
, _D32_dxdxdx, _D32_dxdxdy, _D32_dxdydy, _D32_dydydy :: ( T v )
}
deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
deriving anyclass NFData
@ -126,8 +126,8 @@ data D3𝔸2 v =
-- | \( \mathbb{Z}[x, y, z]/(x, y, z)^2 \).
data D1𝔸3 v =
D13 { _D13_v :: !v
, _D13_dx, _D13_dy, _D13_dz :: !( T v )
D13 { _D13_v :: v
, _D13_dx, _D13_dy, _D13_dz :: ( T v )
}
deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
deriving anyclass NFData
@ -136,8 +136,8 @@ data D1𝔸3 v =
-- | \( \mathbb{Z}[x, y, z]/(x, y, z)^3 \).
data D2𝔸3 v =
D23 { _D23_v :: !v
, _D23_dx, _D23_dy, _D23_dz :: !( T v )
D23 { _D23_v :: v
, _D23_dx, _D23_dy, _D23_dz :: ( T v )
, _D23_dxdx, _D23_dxdy, _D23_dydy, _D23_dxdz, _D23_dydz, _D23_dzdz :: !( T v )
}
deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
@ -147,8 +147,8 @@ data D2𝔸3 v =
-- | \( \mathbb{Z}[x, y, z]/(x, y, z)^4 \).
data D3𝔸3 v =
D33 { _D33_v :: !v
, _D33_dx, _D33_dy, _D33_dz :: !( T v )
D33 { _D33_v :: v
, _D33_dx, _D33_dy, _D33_dz :: ( T v )
, _D33_dxdx, _D33_dxdy, _D33_dydy, _D33_dxdz, _D33_dydz, _D33_dzdz :: !( T v )
, _D33_dxdxdx, _D33_dxdxdy, _D33_dxdydy, _D33_dydydy
, _D33_dxdxdz, _D33_dxdydz, _D33_dxdzdz, _D33_dydydz, _D33_dydzdz, _D33_dzdzdz :: !( T v )
@ -160,8 +160,8 @@ data D3𝔸3 v =
-- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^2 \).
data D1𝔸4 v =
D14 { _D14_v :: !v
, _D14_dx, _D14_dy, _D14_dz, _D14_dw :: !( T v )
D14 { _D14_v :: v
, _D14_dx, _D14_dy, _D14_dz, _D14_dw :: ( T v )
}
deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
deriving anyclass NFData
@ -170,10 +170,10 @@ data D1𝔸4 v =
-- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^3 \).
data D2𝔸4 v =
D24 { _D24_v :: !v
, _D24_dx, _D24_dy, _D24_dz, _D24_dw :: !( T v )
D24 { _D24_v :: v
, _D24_dx, _D24_dy, _D24_dz, _D24_dw :: ( T v )
, _D24_dxdx, _D24_dxdy, _D24_dydy, _D24_dxdz
, _D24_dydz, _D24_dzdz, _D24_dxdw, _D24_dydw, _D24_dzdw, _D24_dwdw :: !( T v )
, _D24_dydz, _D24_dzdz, _D24_dxdw, _D24_dydw, _D24_dzdw, _D24_dwdw :: ( T v )
}
deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
deriving anyclass NFData
@ -182,14 +182,14 @@ data D2𝔸4 v =
-- | \( \mathbb{Z}[x, y, z, w]/(x, y, z, w)^4 \).
data D3𝔸4 v =
D34 { _D34_v :: !v
, _D34_dx, _D34_dy, _D34_dz, _D34_dw :: !( T v )
D34 { _D34_v :: v
, _D34_dx, _D34_dy, _D34_dz, _D34_dw :: ( T v )
, _D34_dxdx, _D34_dxdy, _D34_dydy, _D34_dxdz, _D34_dydz, _D34_dzdz
, _D34_dxdw, _D34_dydw, _D34_dzdw, _D34_dwdw :: !( T v )
, _D34_dxdw, _D34_dydw, _D34_dzdw, _D34_dwdw :: ( T v )
, _D34_dxdxdx, _D34_dxdxdy, _D34_dxdydy, _D34_dydydy,
_D34_dxdxdz, _D34_dxdydz, _D34_dxdzdz, _D34_dydzdz, _D34_dydydz, _D34_dzdzdz,
_D34_dxdxdw, _D34_dxdydw, _D34_dydydw, _D34_dxdzdw, _D34_dydzdw, _D34_dzdzdw,
_D34_dxdwdw, _D34_dydwdw, _D34_dzdwdw, _D34_dwdwdw :: !( T v )
_D34_dxdwdw, _D34_dydwdw, _D34_dzdwdw, _D34_dwdwdw :: ( T v )
}
deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 )
deriving anyclass NFData

View file

@ -8,7 +8,7 @@ module Math.Bezier.Cubic
, bezier, bezier', bezier'', bezier'''
, derivative
, curvature, squaredCurvature, signedCurvature
, subdivide
, subdivide, restrict
, ddist, closestPoint
, drag, selfIntersectionParameters
, extrema
@ -180,6 +180,17 @@ subdivide ( Bezier {..} ) t = ( Bezier p0 q1 q2 pt, Bezier pt r1 r2 p3 )
pt = lerp @v t q2 r1
{-# INLINEABLE subdivide #-}
-- | Restrict a cubic Bézier curve to a sub-interval, re-parametrising
-- to \( [0,1] \).
restrict :: forall v r p. ( Torsor v p, Ring.Field r, Module r v ) => Bezier p -> ( r , r ) -> Bezier p
restrict bez ( a, b ) = fst $ ( flip ( subdivide @v ) b' ) $ snd $ subdivide @v bez a
where
b' = ( b Ring.- a ) Ring./ ( Ring.fromInteger 1 Ring.- a )
-- TODO: this could be made more efficient.
-- See e.g. "https://math.stackexchange.com/questions/4172835/cubic-b%C3%A9zier-spline-multiple-split"
-- or the paper "On the numerical condition of Bernstein-Bézier subdivision process".
{-# INLINEABLE restrict #-}
-- | Polynomial coefficients of the derivative of the distance to a cubic Bézier curve.
ddist :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> p -> [ r ]
ddist ( Bezier {..} ) c = [ a5, a4, a3, a2, a1, a0 ]

View file

@ -6,7 +6,7 @@ module Math.Bezier.Quadratic
( Bezier(..)
, bezier, bezier', bezier''
, curvature, squaredCurvature, signedCurvature
, subdivide
, subdivide, restrict
, ddist, closestPoint
, interpolate
, extrema
@ -141,6 +141,17 @@ subdivide ( Bezier {..} ) t = ( Bezier p0 q1 pt, Bezier pt r1 p2 )
pt = lerp @v t q1 r1
{-# INLINEABLE subdivide #-}
-- | Restrict a quadratic Bézier curve to a sub-interval, re-parametrising
-- to \( [0,1] \).
restrict :: forall v r p. ( Torsor v p, Ring.Field r, Module r v ) => Bezier p -> ( r , r ) -> Bezier p
restrict bez ( a, b ) = fst $ ( flip ( subdivide @v ) b' ) $ snd $ subdivide @v bez a
where
b' = ( b Ring.- a ) Ring./ ( Ring.fromInteger 1 Ring.- a )
-- TODO: this could be made more efficient.
-- See e.g. "https://math.stackexchange.com/questions/4172835/cubic-b%C3%A9zier-spline-multiple-split"
-- or the paper "On the numerical condition of Bernstein-Bézier subdivision process".
{-# INLINEABLE restrict #-}
-- | Polynomial coefficients of the derivative of the distance to a quadratic Bézier curve.
ddist :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> p -> [ r ]
ddist ( Bezier {..} ) c = [ a3, a2, a1, a0 ]

View file

@ -148,6 +148,7 @@ import Math.Orientation
( Orientation(..), splineOrientation
, between
)
import qualified Math.Ring as Ring
import Math.Roots
import Debug.Utils
@ -600,7 +601,7 @@ outlineFunction rootAlgo ptParams toBrushParams brushFromParams = \ sp0 crv ->
( newtTrees, ( newtDunno, newtSols ) ) =
intervalNewtonGS
NoPreconditioning --InverseMidJacobian
InverseMidJacobian
1e-7
curvesI
@ -634,7 +635,7 @@ outlineFunction rootAlgo ptParams toBrushParams brushFromParams = \ sp0 crv ->
logContents = unlines $ functionDataLines ++ treeLines
in trace (unlines solLines) $
logToFile cuspFindingMVar logContents `seq`
--logToFile cuspFindingMVar logContents `seq`
OutlineInfo
{ outlineFn = fwdBwd
, outlineDefiniteCusps = map ( cuspCoords curves ) newtSols
@ -1112,7 +1113,7 @@ solveEnvelopeEquations rootAlgo _t path_t path'_t ( fwdOffset, bwdOffset ) strok
++ "}"
]
in logToFile rootSolvingMVar logContents `seq`
in --logToFile rootSolvingMVar logContents `seq`
( good, ds, dcdt )
(runSolveMethod, methodName) = case rootAlgo of
@ -1206,7 +1207,7 @@ brushStrokeData co1 co2 path params brush =
splines :: Seq ( D k ( I i brushParams ) ( I i ( 1 ) `arr` I i ( 2 ) ) )
!splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns @k @i co2 ) dbrush_params
dbrushes_t :: Seq ( I i ( 1 ) -> D k ( I i ( 2 ) ) ( I i ( 2 ) ) )
!dbrushes_t = force $ fmap ( uncurryD @k . chain @(I i Double) @k dparams_t ) splines
!dbrushes_t = force $ fmap ( uncurryD @k . chain @( I i Double ) @k dparams_t ) splines
-- This is the crucial use of the chain rule.
in fmap ( mkStrokeDatum dpath_t ) dbrushes_t
@ -1274,6 +1275,115 @@ gaussSeidel
return ( ( 𝕀 ( 1 x1'_lo ) ( 1 x1'_hi ), 𝕀 ( 1 x2'_lo ) ( 1 x2'_hi ) )
, sub_x1 && sub_x2 )
{-
gaussSeidel2 :: Int
-> Double
-> Double
-> ( 𝕀 2, 𝕀 2 ) -- ^ columns of \( A \)
-> 𝕀 2 -- ^ \( B \)
-> ( 𝕀 1, 𝕀 1 ) -- ^ initial box \( X \)
-> [ ( ( 𝕀 1, 𝕀 1 ), Bool ) ]
gaussSeidel2 maxIters eps_abs eps_rel
( 𝕀 ( 2 a11_lo a21_lo ) ( 2 a11_hi a21_hi )
, 𝕀 ( 2 a12_lo a22_lo ) ( 2 a12_hi a22_hi ) )
( 𝕀 ( 2 b1_lo b2_lo ) ( 2 b1_hi b2_hi ) )
x0
= let !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
-- See "Algorithm 2" in
-- "Using interval unions to solve linear systems of equations with uncertainties"
iter ( 𝕀 ( 1 x1_lo ) ( 1 x1_hi ), 𝕀 ( 1 x2_lo ) ( 1 x2_hi ) ) = do
let
!x1 = 𝕀 x1_lo x1_hi
!x2 = 𝕀 x2_lo x2_hi
blah1 = do
let s = b1 - a12 * x2
let s1 = s `monus` ( a11 * x1 )
y1 <-
if not $ containsZero ( s1 - a11 * x1 )
then []
else
if containsZero s1 && containsZero a11
then [ ( x1, False ) ]
else do
x1'' <- s1 `extendedDivide` a11
x1'' `intersect` x1
let s2 = s `monus` ( a12 * x2 )
y2 <-
if not $ containsZero ( s2 - a11 * x1 )
then []
else
if containsZero s2 && containsZero a11
then [ ( x1, False ) ]
else do
x1'' <- s2 `extendedDivide` a12
x1'' `intersect` x1
return ( y1 `cart` y2 )
blah2 = do
let s = b2 - a21 * x1
let s1 = s `monus` ( a21 * x1 )
y1 <-
if not $ containsZero ( s1 - a22 * x2 )
then []
else
if containsZero s1 && containsZero a22
then [ ( x1, False ) ]
else do
x1'' <- s1 `extendedDivide` a21
x1'' `intersect` x1
let s2 = s `monus` ( a12 * x2 )
y2 <-
if not $ containsZero ( s2 - a22 * x2 )
then []
else
if containsZero s2 && containsZero a22
then [ ( x1, False ) ]
else do
x1'' <- s2 `extendedDivide` a22
x1'' `intersect` x1
return ( y1 `cart` y2 )
blah1 ++ blah2
go :: Int -> ( 𝕀 1, 𝕀 1 ) -> [ ( ( 𝕀 1, 𝕀 1 ), Bool ) ]
go !i x
= do { nxt@( x', sub ) <- iter x
; let dw_abs = maxWidth x - maxWidth x'
dw_rel = 1 - ( maxWidth x' / maxWidth x )
; if sub
|| i >= maxIters
|| dw_abs < eps_abs
|| dw_rel < eps_rel
then return nxt
else go ( i + 1 ) x'
}
in go 1 x0
where
maxWidth :: ( 𝕀 1, 𝕀 1 ) -> Double
maxWidth ( 𝕀 ( 1 x1_lo ) ( 1 x1_hi ), 𝕀 ( 1 x2_lo ) ( 1 x2_hi ) )
= max ( x1_hi - x1_lo ) ( x2_hi - x2_lo )
containsZero :: 𝕀 Double -> Bool
containsZero ( 𝕀 lo hi ) = lo <= 0 && hi >= 0
monus :: 𝕀 Double -> 𝕀 Double -> 𝕀 Double
monus ( 𝕀 lo1 hi1 ) ( 𝕀 lo2 hi2 )
| hi1 - lo1 >= hi2 - lo2
= 𝕀 ( lo1 - lo2 ) ( hi1 - hi2 )
| otherwise
= 𝕀 ( hi1 - hi2 ) ( lo1 - lo2 )
cart :: ( 𝕀 Double, Bool ) -> ( 𝕀 Double, Bool ) -> ( ( 𝕀 1, 𝕀 1 ), Bool )
cart ( 𝕀 lo1 hi1, sub1 ) ( 𝕀 lo2 hi2, sub2 ) =
( ( 𝕀 ( 1 lo1 ) ( 1 hi1 ), 𝕀 ( 1 lo2 ) ( 1 hi2 ) ), sub1 && sub2 )
-}
-- | Compute the intersection of two intervals.
--
-- Returns whether the first interval is a strict subset of the second interval,
@ -1342,7 +1452,7 @@ data IntervalNewtonLeaf d
showIntervalNewtonTree :: Box -> IntervalNewtonTree Box -> Tree String
showIntervalNewtonTree cand (IntervalNewtonLeaf l) = Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ show l) []
showIntervalNewtonTree cand (IntervalNewtonStep s ts)
= Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ show s) $ map (\ (c,t) -> showIntervalNewtonTree c t) ts
= Node (show cand ++ " abc " ++ showArea (boxArea cand) ++ " " ++ show s) $ map (\ (c,t) -> showIntervalNewtonTree c t) ts
boxArea :: Box -> Double
boxArea ( 𝕀 ( 1 t_lo ) ( 1 t_hi ), _, 𝕀 ( 1 s_lo ) ( 1 s_hi ) )
@ -1383,96 +1493,152 @@ intervalNewtonGSFrom
-> ( [ ( Box, IntervalNewtonTree Box ) ], ( [ Box ], [ Box ] ) )
intervalNewtonGSFrom precondMethod minWidth eqs initBox =
runWriter $
map ( initBox , ) <$> go initBox
map ( initBox , ) <$> evalStrokeDataAndGo initBox
where
recur f cands = do
rest <- traverse ( \ c -> do { trees <- go c; return [ (c, tree) | tree <- trees ] } ) cands
recur :: ( cand -> Writer ( [ Box ], [ Box ] ) [ IntervalNewtonTree Box ] )
-> ( [ ( cand, IntervalNewtonTree Box ) ] -> IntervalNewtonTree Box )
-> [ cand ]
-> Writer ( [Box], [Box] ) [ IntervalNewtonTree Box ]
recur r f cands = do
rest <- traverse ( \ c -> do { trees <- r c; return [ (c, tree) | tree <- trees ] } ) cands
return [ f $ concat rest ]
go :: Box -- box to work on
evalStrokeDataAndGo :: Box -> Writer ( [Box], [Box] ) [ IntervalNewtonTree Box ]
evalStrokeDataAndGo box@( t, i, s ) = go ( box, ( eqs t `Seq.index` i ) s )
go :: ( Box, StrokeDatum 3 𝕀 ) -- box to work on
-> Writer ( [ Box ], [ Box ] )
[ IntervalNewtonTree Box ]
go cand@( t@( 𝕀 ( 1 t_lo ) ( 1 t_hi ) )
, i
, s@( 𝕀 ( 1 s_lo ) ( 1 s_hi ) )
)
go ( cand@( t@( 𝕀 ( 1 t_lo ) ( 1 t_hi ) )
, i
, s@( 𝕀 ( 1 s_lo ) ( 1 s_hi ) )
)
, sd@( StrokeDatum { ee = D22 _ ( T ee_t ) (T ee_s ) _ _ _
, 𝛿E𝛿sdcdt = D12 _ ( T ( T f_t ) ) ( T ( T f_s ) ) })
)
-- Box is small: stop processing it.
| t_hi - t_lo < minWidth && s_hi - s_lo < minWidth
= do let res = TooSmall cand
tell ( [ cand ], [] )
return [ IntervalNewtonLeaf res ]
| StrokeDatum { ee = D22 ee ( T _ee_t ) ( T _ee_s ) _ _ _
, 𝛿E𝛿sdcdt = D12 ( T f ) ( T ( T f_t ) ) ( T ( T f_s ) ) }
<- ( eqs t `Seq.index` i ) s
| -- Check the envelope equation interval contains zero.
ee_potential_zero sd
-- Check the 𝛿E𝛿sdcdt interval contains zero.
, 𝛿E𝛿sdcdt_potential_zero sd
, StrokeDatum { ee = D22 _ee_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 = cmp2 (<=) ( inf f ) ( 2 0 0 ) && cmp2 (>=) ( sup f ) ( 2 0 0 )
= if | -- Check the envelope equation interval contains zero.
ee_potential_zero
-- Check the 𝛿E𝛿sdcdt interval contains zero.
, 𝛿E𝛿sdcdt_potential_zero
-> let -- Interval Newton method: take one GaussSeidel step
-- for the equation f'(x) v = - f(x_mid),
-- where f = 𝛿E/𝛿s * dc/dt
!( a, b ) = precondition precondMethod
( midI f_t, midI f_s )
( f_t, f_s ) ( neg f_mid )
--(a, b)
-- | 𝕀 (1 ee_lo) (1 ee_hi) <- ee_mid
-- , 𝕀 (1 ee_t_lo) (1 ee_t_hi) <- ee_t
-- , 𝕀 (1 ee_s_lo) (1 ee_s_hi) <- ee_s
-- , 𝕀 (2 fx_lo fy_lo) (2 fx_hi fy_hi) <- f_mid
-- , 𝕀 (2 f_tx_lo f_ty_lo) (2 f_tx_hi f_ty_hi) <- f_t
-- , 𝕀 (2 f_sx_lo f_sy_lo) (2 f_sx_hi f_sy_hi) <- f_s
-- = ( ( 𝕀 (2 f_tx_lo ee_t_lo) (2 f_tx_hi ee_t_hi)
-- , 𝕀 (2 f_sx_lo ee_s_lo) (2 f_sx_hi ee_s_hi)
-- )
-- , neg $ 𝕀 (2 fx_lo ee_lo) (2 fx_hi ee_hi)
-- )
= let -- Interval Newton method: take one GaussSeidel step
-- for the equation f'(x) v = - f(x_mid),
-- where f = 𝛿E/𝛿s * dc/dt
!( a, b ) = precondition precondMethod
( midI f_t, midI f_s )
( f_t, f_s ) ( neg f_mid )
--(a, b)
-- | 𝕀 (1 ee_lo) (1 ee_hi) <- ee_mid
-- , 𝕀 (1 ee_t_lo) (1 ee_t_hi) <- ee_t
-- , 𝕀 (1 ee_s_lo) (1 ee_s_hi) <- ee_s
-- , 𝕀 (2 fx_lo fy_lo) (2 fx_hi fy_hi) <- f_mid
-- , 𝕀 (2 f_tx_lo f_ty_lo) (2 f_tx_hi f_ty_hi) <- f_t
-- , 𝕀 (2 f_sx_lo f_sy_lo) (2 f_sx_hi f_sy_hi) <- f_s
-- = ( ( 𝕀 (2 f_tx_lo ee_t_lo) (2 f_tx_hi ee_t_hi)
-- , 𝕀 (2 f_sx_lo ee_s_lo) (2 f_sx_hi ee_s_hi)
-- )
-- , neg $ 𝕀 (2 fx_lo ee_lo) (2 fx_hi ee_hi)
-- )
!gsGuesses = gaussSeidel a b
( coerce ( (-) @( 𝕀 Double ) ) t i_t_mid
, coerce ( (-) @( 𝕀 Double ) ) s i_s_mid )
in if all ( smaller . fst ) gsGuesses
then
-- If the GaussSeidel 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 ( mkGuess . fst ) ) ( map ( mkGuess . fst ) )
$ partition snd gsGuesses
in do tell ([], done)
case todo of
[] -> return [ IntervalNewtonLeaf $ NoSolution "GaussSeidel" cand ]
_ -> recur (IntervalNewtonStep (IntervalNewtonContraction done)) todo
else
-- GaussSeidel failed to shrink the boxes.
-- Bisect along the widest dimension instead.
let (bisGuesses, whatBis)
| t_hi - t_lo > s_hi - s_lo
= ( [ ( 𝕀 ( 1 t_lo ) ( 1 t_mid ), i, s )
, ( 𝕀 ( 1 t_mid ) ( 1 t_hi ), i, s ) ]
, ("t", t_mid) )
| otherwise
= ( [ ( t, i, 𝕀 ( 1 s_lo ) ( 1 s_mid ) )
, ( t, i, 𝕀 ( 1 s_mid ) ( 1 s_hi ) ) ]
, ("s", s_mid) )
in recur (IntervalNewtonStep (IntervalNewtonBisection whatBis)) bisGuesses
!gsGuesses = gaussSeidel a b
( coerce ( (-) @( 𝕀 Double ) ) t i_t_mid
, coerce ( (-) @( 𝕀 Double ) ) s i_s_mid )
in if any ( smaller . fst ) gsGuesses
then
-- If the GaussSeidel 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 ( mkGuess . fst ) ) ( map ( mkGuess . fst ) )
$ partition snd gsGuesses
in do tell ([], done)
case todo of
[] -> return [ IntervalNewtonLeaf $ NoSolution "GaussSeidel" cand ]
_ -> recur evalStrokeDataAndGo ( IntervalNewtonStep ( IntervalNewtonContraction done ) )
todo
else
-- GaussSeidel failed to shrink the boxes, so bisect instead.
-- We have to pick along which dimension to bisect:
-- - if bisecting along a particular dimension discards one of
-- the boxes, do that;
-- - otherwise, bisect along the dimension j that maximises
-- width(x_j) * || J_j ||.
let l_t = 𝕀 ( 1 t_lo ) ( 1 t_mid )
r_t = 𝕀 ( 1 t_mid ) ( 1 t_hi )
d_s = 𝕀 ( 1 s_lo ) ( 1 s_mid )
u_s = 𝕀 ( 1 s_mid ) ( 1 s_hi )
l = ( l_t, i, s )
r = ( r_t, i, s )
d = ( t, i, d_s )
u = ( t, i, u_s )
l_dat = ( eqs l_t `Seq.index` i ) s
r_dat = ( eqs r_t `Seq.index` i ) s
d_dat = ( eqs t `Seq.index` i ) d_s
u_dat = ( eqs t `Seq.index` i ) u_s
l_skip =
not ( ee_potential_zero l_dat )
|| not ( 𝛿E𝛿sdcdt_potential_zero l_dat )
r_skip =
not ( ee_potential_zero r_dat )
|| not ( 𝛿E𝛿sdcdt_potential_zero r_dat )
d_skip =
not ( ee_potential_zero d_dat )
|| not ( 𝛿E𝛿sdcdt_potential_zero d_dat )
u_skip =
not ( ee_potential_zero u_dat )
|| not ( 𝛿E𝛿sdcdt_potential_zero u_dat )
tJWidth = ( t_hi - t_lo ) * normVI3 ee_t f_t
sJWidth = ( s_hi - s_lo ) * normVI3 ee_s f_s
( bisGuesses, whatBis )
| l_skip && r_skip
= ( [], ( "lr", t_mid ) )
| d_skip && u_skip
= ( [], ( "du", s_mid ) )
| l_skip
= ( [ ( r, r_dat ) ], ( "r", t_mid ) )
| r_skip
= ( [ ( l, l_dat ) ], ( "l", t_mid ) )
| d_skip
= ( [ ( u, u_dat ) ], ( "u", s_mid ) )
| u_skip
= ( [ ( d, d_dat ) ], ( "d", s_mid ) )
| tJWidth >= sJWidth
-- t_hi - t_lo >= s_hi - s_lo
= ( [ ( l, l_dat ), ( r, r_dat ) ], ( "t", t_mid ) )
| otherwise
= ( [ ( d, d_dat ), ( u, u_dat ) ], ( "s", s_mid ) )
in recur go ( IntervalNewtonStep ( IntervalNewtonBisection whatBis ) . map (first fst) )
bisGuesses
-- Box doesn't contain a solution: discard it.
| otherwise
-> return [ IntervalNewtonLeaf $ NoSolution (if ee_potential_zero then "dc/dt" else "ee") cand ]
= return [ IntervalNewtonLeaf $ NoSolution ( if ee_potential_zero sd 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 ) )
--width :: 𝕀 1 -> Double
--width ( 𝕀 ( 1 lo ) ( 1 hi ) ) = hi - lo
--normI :: 𝕀 1 -> Double
--normI ( 𝕀 ( 1 lo ) ( 1 hi ) ) = sqrt $ sup $ ( 𝕀 lo hi ) Ring.^ 2
--normVI :: 𝕀 2 -> Double
--normVI ( 𝕀 ( 2 x_lo y_lo ) ( 2 x_hi y_hi ) ) =
-- sqrt $ sup $ ( 𝕀 x_lo x_hi ) Ring.^ 2 + ( 𝕀 y_lo y_hi ) Ring.^ 2
normVI3 :: 𝕀 1 -> 𝕀 2 -> Double
normVI3 ( 𝕀 ( 1 lo ) ( 1 hi ) ) ( 𝕀 ( 2 x_lo y_lo ) ( 2 x_hi y_hi ) )
= sqrt $ max ( abs lo ) ( abs hi ) Ring.^ 2
+ max ( abs x_lo ) ( abs x_hi ) Ring.^ 2
+ max ( abs y_lo ) ( abs y_hi ) Ring.^ 2
t_mid = 0.5 * ( t_lo + t_hi )
s_mid = 0.5 * ( s_lo + s_hi )
i_t_mid = singleton ( 1 t_mid )
@ -1490,6 +1656,16 @@ intervalNewtonGSFrom precondMethod minWidth eqs initBox =
!( 𝕀 y'_lo y'_hi ) = negate $ 𝕀 y_lo y_hi
in 𝕀 ( 2 x'_lo y'_lo ) ( 2 x'_hi y'_hi )
ee_potential_zero :: StrokeDatum 3 𝕀 -> Bool
ee_potential_zero dat =
inf ( _D22_v $ ee dat ) <= 1 0
&& sup ( _D22_v $ ee dat ) >= 1 0
𝛿E𝛿sdcdt_potential_zero :: StrokeDatum 3 𝕀 -> Bool
𝛿E𝛿sdcdt_potential_zero dat =
cmp2 (<=) ( inf $ unT $ _D12_v $ 𝛿E𝛿sdcdt dat ) ( 2 0 0 )
&& cmp2 (>=) ( sup $ unT $ _D12_v $ 𝛿E𝛿sdcdt dat ) ( 2 0 0 )
zero, one :: Double
zero = 1e-6
one = 1 - zero

View file

@ -332,3 +332,31 @@ evaluateQuadratic bez t =
maxs = fmap (Quadratic.bezier @( T Double ) sup_bez)
$ inf t :| ( sup t : filter ( inside t ) ( Quadratic.extrema sup_bez ) )
in 𝕀 ( minimum mins ) ( maximum maxs )
{-
evaluateCubic :: Cubic.Bezier ( 𝕀 Double ) -> 𝕀 Double -> 𝕀 Double
evaluateCubic bez t =
-- assert (inf t >= 0 && sup t <= 1) "evaluateCubic: t ⊊ [0,1]" $ -- Requires t ⊆ [0,1]
let inf_bez = Cubic.restrict @( T Double ) ( fmap inf bez ) ( inf t, sup t )
sup_bez = Cubic.restrict @( T Double ) ( fmap sup bez ) ( inf t, sup t )
mins = fmap (Cubic.bezier @( T Double ) inf_bez)
$ 0 :| ( 1 : Cubic.extrema inf_bez )
maxs = fmap (Cubic.bezier @( T Double ) sup_bez)
$ 0 :| ( 1 : Cubic.extrema sup_bez )
in 𝕀 ( minimum mins ) ( maximum maxs )
-- | Evaluate a quadratic Bézier curve, when both the coefficients and the
-- parameter are intervals.
evaluateQuadratic :: Quadratic.Bezier ( 𝕀 Double ) -> 𝕀 Double -> 𝕀 Double
evaluateQuadratic bez t =
-- assert (inf t >= 0 && sup t <= 1) "evaluateCubic: t ⊊ [0,1]" $ -- Requires t ⊆ [0,1]
let inf_bez = Quadratic.restrict @( T Double ) ( fmap inf bez ) ( inf t, sup t )
sup_bez = Quadratic.restrict @( T Double ) ( fmap sup bez ) ( inf t, sup t )
mins = fmap (Quadratic.bezier @( T Double ) inf_bez)
$ 0 :| ( 1 : Quadratic.extrema inf_bez )
maxs = fmap (Quadratic.bezier @( T Double ) sup_bez)
$ 0 :| ( 1 : Quadratic.extrema sup_bez )
in 𝕀 ( minimum mins ) ( maximum maxs )
-}

View file

@ -120,19 +120,7 @@ instance Prelude.Fractional ( 𝕀 Double ) where
| otherwise
= error "BAD interval recip; should use extendedRecip instead"
-- #endif
𝕀 x_lo x_hi / 𝕀 y_lo y_hi
-- #ifdef ASSERTS
| 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
| otherwise
= error "BAD interval division; should use extendedRecip instead"
-- #endif
p / q = p * recip q
instance Floating ( 𝕀 Double ) where
sqrt = withHW Prelude.sqrt

View file

@ -158,7 +158,7 @@ runApplication application = do
{ strokeName = "Stroke 1"
, strokeVisible = True
, strokeUnique = strokeUnique
, strokeBrush = Just Asset.Brushes.tearDrop
, strokeBrush = Just Asset.Brushes.ellipse --tearDrop
, strokeSpline =
-- Spline
-- { splineStart = mkPoint ( 2 -20 -20 ) 5
@ -169,7 +169,7 @@ runApplication application = do
Spline
{ splineStart = mkPoint ( 2 0 0 ) 10 25 0
, splineCurves = OpenCurves $ Seq.fromList
[ LineTo { curveEnd = NextPoint ( mkPoint ( 2 100 0 ) 15 40 pi ), curveData = invalidateCache undefined }
[ LineTo { curveEnd = NextPoint ( mkPoint ( 2 100 0 ) 15 40 (0.1 * pi) ), curveData = invalidateCache undefined }
--, LineTo { curveEnd = NextPoint ( mkPoint ( 2 -10 10 ) 8 5 ( pi / 4 ) ), curveData = invalidateCache undefined }
--, LineTo { curveEnd = NextPoint ( mkPoint ( 2 -10 -20 ) 10 7 ( pi / 2 ) ), curveData = invalidateCache undefined }
]
@ -179,12 +179,12 @@ runApplication application = do
)
]
where
--mkPoint :: 2 -> Double -> Double -> Double -> PointData ( Record Asset.Brushes.EllipseBrushFields )
--mkPoint pt a b phi = PointData pt Normal ( MkR $ 3 a b phi )
mkPoint :: 2 -> Double -> Double -> Double -> PointData ( Record Asset.Brushes.EllipseBrushFields )
mkPoint pt a b phi = PointData pt Normal ( MkR $ 3 a b phi )
--mkPoint :: 2 -> Double -> PointData ( Record Asset.Brushes.CircleBrushFields )
--mkPoint pt r = PointData pt Normal ( MkR $ 1 r )
mkPoint :: 2 -> Double -> Double -> Double -> PointData ( Record Asset.Brushes.TearDropBrushFields )
mkPoint pt w h phi = PointData pt Normal ( MkR $ 3 w h phi )
--mkPoint :: 2 -> Double -> Double -> Double -> PointData ( Record Asset.Brushes.TearDropBrushFields )
--mkPoint pt w h phi = PointData pt Normal ( MkR $ 3 w h phi )
recomputeStrokesTVar <- STM.newTVarIO @Bool False
documentRenderTVar <- STM.newTVarIO @( ( Int32, Int32 ) -> Cairo.Render () ) ( const $ pure () )

View file

@ -139,12 +139,8 @@ instance ( Torsor ( T ( 𝕀 ( Length ks ) ) ) ( 𝕀 ( Length ks ) )
T ( 𝕀 ( MkR c_lo ) ( MkR c_hi ) )
type instance RepDim ( Record ks ) = Length ks
instance Representable r ( ( Length ks ) )
=> Representable r ( Record ks ) where
tabulate f = Record $ tabulate f
{-# INLINE tabulate #-}
index f (Record r) = index f r
{-# INLINE index #-}
deriving newtype instance Representable r ( ( Length ks ) )
=> Representable r ( Record ks )
type instance D k ( Record ks ) = D k ( ( Length ks ) )
deriving newtype instance HasChainRule Double 2 ( ( Length ks ) )