hashirama/metabrush
1
0
Fork 0
mirror of https://gitlab.com/sheaf/metabrush.git synced 2024-11-27 09:24:08 +00:00
Code Issues Actions Projects Releases Wiki Activity

start modularising cusp finding code

Browse source
This commit is contained in:
sheaf 2024-04-05 17:48:53 +02:00
parent 55470d1f0e
commit bd468fcf82
2 changed files with 141 additions and 237 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

29
brush-strokes/src/cusps/bench/Main.hs
View file

@ -118,7 +118,9 @@ benchTestCase ( TestCase { testName, testBrushStroke, testCuspOptions, testStart
, "Test case " ++ testName ]
before <- getMonotonicTime
let ( _, testStrokeFnI ) = brushStrokeFunctions testBrushStroke
( _newtTrees, ( dunno, sols ) ) = computeCusps testCuspOptions testStrokeFnI testStartBoxes
( dunno, sols ) =
foldMap ( \ ( i, ( _trees, ( mbCusps, defCusps ) ) ) -> ( map ( i , ) mbCusps, map ( i, ) defCusps ) ) $
computeCusps testCuspOptions testStrokeFnI testStartBoxes
rnf dunno `seq` rnf sols `seq` return ()
after <- getMonotonicTime
let dt = after - before
@ -192,7 +194,7 @@ data TestCase =
{ testName :: String
, testBrushStroke :: BrushStroke
, testCuspOptions :: CuspFindingOptions
, testStartBoxes :: [ Box ]
, testStartBoxes :: [ ( Int, Box ) ]
}
brushStrokeFunctions
@ -216,9 +218,9 @@ mkVal t i s = ( ℝ1 t, i, ℝ1 s )
mkI :: ( Double, Double ) -> 𝕀 Double
mkI ( lo, hi ) = 𝕀 lo hi
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 ) )
mkBox :: ( Double, Double ) -> ( Double, Double ) -> Box
mkBox ( t_min, t_max ) ( s_min, s_max ) =
( 𝕀 ( ℝ2 t_min s_min ) ( ℝ2 t_max s_max ) )
zero, one :: Double
zero = 1e-6
@ -480,7 +482,7 @@ showD float = showFFloat (Just 6) float ""
--------------------------------------------------------------------------------
ellipseTestCase :: CuspFindingOptions -> String -> ( Double, Double ) -> Double -> [ Box ] -> TestCase
ellipseTestCase :: CuspFindingOptions -> String -> ( Double, Double ) -> Double -> [ ( Int, Box ) ] -> TestCase
ellipseTestCase opts str k0k1 rot startBoxes =
TestCase
{ testName = "ellipse (" ++ str ++ ")"
@ -602,14 +604,19 @@ getStrokeFunctions brush sp0 crv =
computeCusps
:: CuspFindingOptions
-> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-> [ Box ]
-> ( [ ( Box, RootFindingTree Box ) ], ( [ Box ], [ Box ] ) )
-> [ ( Int, Box ) ]
-> [ ( Int, ( [ ( Box, RootFindingTree Box ) ], ( [ Box ], [ Box ] ) ) ) ]
computeCusps cuspOpts eqs startBoxes =
foldMap ( findCuspsFrom cuspOpts eqs ) startBoxes
map ( \ ( i, box ) -> ( i , ) $ findCuspsFrom cuspOpts ( mkEqn i ) box ) startBoxes
where
mkEqn 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 )
in ( eqs t `Seq.index` i ) s
defaultStartBoxes :: [ Int ] -> [ Box ]
defaultStartBoxes :: [ Int ] -> [ ( Int, Box ) ]
defaultStartBoxes is =
[ mkBox (zero, one) i (zero, one) | i <- is ]
[ ( i, mkBox (zero, one) (zero, one) ) | i <- is ]
getR1 (ℝ1 u) = u

349
brush-strokes/src/lib/Math/Bezier/Stroke.hs
View file

@ -33,8 +33,6 @@ module Math.Bezier.Stroke
-- base
import Control.Arrow
( first, (***) )
import Control.Concurrent.MVar
( MVar, newMVar )
import Control.Monad
( unless, when )
import Control.Monad.ST
@ -50,7 +48,7 @@ import Data.Foldable
import Data.Functor.Identity
( Identity(..) )
import Data.List
( intercalate, nub, partition, sort )
( nub, partition, sort )
import Data.List.NonEmpty
( NonEmpty )
import qualified Data.List.NonEmpty as NE
@ -71,8 +69,6 @@ import GHC.TypeNats
( type (-) )
import Numeric
( showFFloat )
import System.IO.Unsafe
( unsafePerformIO )
-- acts
import Data.Act
@ -83,12 +79,16 @@ import Data.Act
)
-- containers
import Data.IntMap.Strict
( IntMap )
import qualified Data.IntMap.Strict as IntMap
( fromList, toList )
import Data.Sequence
( Seq(..) )
import qualified Data.Sequence as Seq
( empty, index, length, reverse, singleton, zipWith )
import Data.Tree
( Tree(..), foldTree )
( Tree(..) )
-- deepseq
import Control.DeepSeq
@ -114,10 +114,6 @@ import Control.Monad.Trans.State.Strict as State
import Control.Monad.Trans.Writer.CPS
( Writer, WriterT, execWriterT, runWriter, tell )
-- tree-view
import Data.Tree.View
( showTree )
-- MetaBrush
import Math.Algebra.Dual
import qualified Math.Bezier.Cubic as Cubic
@ -154,7 +150,7 @@ import Math.Orientation
import qualified Math.Ring as Ring
import Math.Roots
import Debug.Utils
--import Debug.Utils
--------------------------------------------------------------------------------
@ -604,56 +600,24 @@ outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams brushFromParams =
$ runD ( brushFromParams @2 @() proxy# id )
$ toBrushParams params_t
( newtTrees, ( newtDunno, newtSols ) ) =
( potentialCusps, definiteCusps ) =
case mbCuspOptions of
Just cuspOptions ->
findCusps cuspOptions curvesI
foldMap
( \ ( i, ( _trees, ( potCusps, defCusps ) ) ) ->
( map ( i , ) potCusps, map ( i , ) defCusps )
)
( IntMap.toList $ findCusps cuspOptions curvesI )
Nothing ->
( [], ( [], [] ) )
( [], [] )
showD :: Double -> String
showD float = showFFloat (Just 6) float ""
functionDataLines =
[ "E, dE/ds, dE/ds * dc/dt:"
, "{" ++
(intercalate ","
[ "{" ++ showD t ++ "," ++ showD (s + fromIntegral i) ++ ",{" ++ intercalate "," vals ++ "}}"
| t <- [0.001, 0.01.. 0.999]
, i <- [0 .. (fromIntegral $ length $ curves (ℝ1 0.5)) - 1]
, s <- [0.001, 0.01.. 0.999]
, let StrokeDatum
{ ee = D12 (ℝ1 e) _dEdt (T (ℝ1 dEds))
, 𝛿E𝛿sdcdt = D0 (T (ℝ2 vx vy))
} = (curves (ℝ1 t) `Seq.index` i) (ℝ1 s)
vals = [showD e, showD dEds, "{" ++ showD vx ++ "," ++ showD vy ++ "}"]
]
) ++ "}"
]
showedTrees = map ( uncurry showRootFindingTree ) newtTrees
solLines =
[ " #sols: " ++ show (length newtSols)
, "#dunno: " ++ show (length newtDunno)
, "Tree size: " ++ show (sum @_ @Int $ map (foldTree ( \ _ bs -> 1 + sum bs )) showedTrees)
]
treeLines =
solLines ++ map showTree showedTrees
logContents = unlines $ functionDataLines ++ treeLines
in --trace (unlines solLines) $
--logToFile cuspFindingMVar logContents `seq`
OutlineInfo
in OutlineInfo
{ outlineFn = fwdBwd
, outlineDefiniteCusps = map ( cuspCoords curves ) newtSols
, outlinePotentialCusps = map ( cuspCoords curves ) newtDunno
, outlineDefiniteCusps = map ( cuspCoords curves ) definiteCusps
, outlinePotentialCusps = map ( cuspCoords curves ) potentialCusps
}
{-# INLINEABLE outlineFunction #-}
cuspFindingMVar :: MVar FilePath
cuspFindingMVar = unsafePerformIO
$ newMVar "logs/cusps.txt"
{-# NOINLINE cuspFindingMVar #-}
pathAndUsedParams :: forall k i arr crvData ptData usedParams
. ( HasType ( ℝ 2 ) ptData
, HasBézier k i
@ -1015,11 +979,6 @@ data RootSolvingAlgorithm
-- | Use the modified Halley M2 method.
| HalleyM2 { maxIters :: Word, precision :: Int }
rootSolvingMVar :: MVar FilePath
rootSolvingMVar = unsafePerformIO
$ newMVar "logs/envelope.txt"
{-# NOINLINE rootSolvingMVar #-}
-- | Solve the envelope equations at a given point \( t = t_0 \), to find
-- \( s_0 \) such that \( c(t_0, s_0) \) is on the envelope of the brush stroke.
solveEnvelopeEquations :: RootSolvingAlgorithm
@ -1055,68 +1014,21 @@ solveEnvelopeEquations rootAlgo _t path_t path'_t ( fwdOffset, bwdOffset ) strok
plausibleTangent tgt = path'_t ^.^ tgt > 0
sol :: String -> Seq ( ℝ 1 -> StrokeDatum 2 () ) -> Double -> ( Bool, ℝ 2, T ( ℝ 2 ) )
sol desc f is0 =
let (solRes, solSteps) = runSolveMethod ( eqn f ) is0
sol _desc f is0 =
let (solRes, _solSteps) = runSolveMethod ( eqn f ) is0
( good, is ) =
case solRes of
Nothing -> ( False, is0 )
Just is1 -> ( True , is1 )
( ds, dcdt ) = finish f is
showD :: Double -> String
showD float = showFFloat (Just 6) float ""
(x_min, x_max) = domain
logContents = unlines
[ "method = " <> methodName
, desc
, "t = " ++ show _t
, "start (n,s) = " <> show (fromDomain is0)
, "steps = {"
, unlines $ map (\ (x, y, y') -> show (x, (y, y'))) solSteps
, "},"
, "E:"
, "{" ++
(intercalate ","
[ "{" ++ showD x ++ "," ++ showD y ++ "}"
| x <- [x_min, x_min + 0.01 .. x_max]
, let (# y, _y' #) = eqn f x
]
) ++ "}"
, "𝛿E𝛿s:"
, "{" ++
(intercalate ","
[ "{" ++ showD x ++ "," ++ showD y' ++ "}"
| x <- [x_min, x_min + 0.01 .. x_max]
, let (# _y, y' #) = eqn f x
]
) ++ "}"
, "𝛿E𝛿sdcdt (x):"
, "{" ++
(intercalate ","
[ "{" ++ showD x ++ "," ++ showD 𝛿E𝛿sdcdt_x ++ "}"
| x <- [x_min, x_min + 0.01 .. x_max]
, let StrokeDatum { 𝛿E𝛿sdcdt = D0 (T (ℝ2 𝛿E𝛿sdcdt_x _)) } = evalStrokeDatum f x
]
)
++ "}"
, "𝛿E𝛿sdcdt (y):"
, "{" ++
(intercalate ","
[ "{" ++ showD x ++ "," ++ showD 𝛿E𝛿sdcdt_y ++ "}"
| x <- [x_min, x_min + 0.01 .. x_max]
, let StrokeDatum { 𝛿E𝛿sdcdt = D0 (T (ℝ2 _ 𝛿E𝛿sdcdt_y)) } = evalStrokeDatum f x
]
)
++ "}"
]
in --logToFile rootSolvingMVar logContents `seq`
( good, ds, dcdt )
in ( good, ds, dcdt )
(runSolveMethod, methodName) = case rootAlgo of
runSolveMethod = case rootAlgo of
HalleyM2 { maxIters, precision } ->
(halleyM2 maxIters precision, "HalleyM2")
halleyM2 maxIters precision
NewtonRaphson { maxIters, precision } ->
(newtonRaphson maxIters precision domain, "NewtonRaphson")
newtonRaphson maxIters precision domain
finish :: Seq ( ℝ 1 -> StrokeDatum 2 () ) -> Double -> ( ℝ 2, T ( ℝ 2 ) )
finish fs is =
@ -1244,13 +1156,13 @@ brushStrokeData co1 co2 path params brush =
-- The boolean indicates whether the Gauss–Seidel step was a contraction.
gaussSeidel :: ( 𝕀ℝ 2, 𝕀ℝ 2 ) -- ^ columns of \( A \)
-> 𝕀ℝ 2 -- ^ \( B \)
-> ( 𝕀ℝ 1, 𝕀ℝ 1 ) -- ^ initial box \( X \)
-> [ ( ( 𝕀ℝ 1, 𝕀ℝ 1 ), Bool ) ]
-> 𝕀ℝ 2 -- ^ initial box \( X \)
-> [ ( 𝕀ℝ 2, Bool ) ]
gaussSeidel
( 𝕀 ( ℝ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 ) )
( 𝕀 ( ℝ1 x1_lo ) ( ℝ1 x1_hi ), 𝕀 ( ℝ1 x2_lo ) ( ℝ1 x2_hi ) )
( 𝕀 ( ℝ2 x1_lo x2_lo ) ( ℝ2 x1_hi x2_hi ) )
= let !a11 = 𝕀 a11_lo a11_hi
!a12 = 𝕀 a12_lo a12_hi
!a21 = 𝕀 a21_lo a21_hi
@ -1268,7 +1180,7 @@ gaussSeidel
x2'_0 <- ( b2 - a21 * x1' ) `extendedDivide` a22
( 𝕀 x2'_lo x2'_hi, sub_x2 ) <- x2'_0 `intersect` x2
return ( ( 𝕀 ( ℝ1 x1'_lo ) ( ℝ1 x1'_hi ), 𝕀 ( ℝ1 x2'_lo ) ( ℝ1 x2'_hi ) )
return ( ( 𝕀 ( ℝ2 x1'_lo x2'_lo) ( ℝ2 x1'_hi x2'_hi ) )
, sub_x1 && sub_x2 )
-- TODO: try implementing the complete interval union Gauss–Seidel algorithm.
-- See "Algorithm 2" in
@ -1294,9 +1206,9 @@ intersect ( 𝕀 lo1 hi1 ) ( 𝕀 lo2 hi2 )
-- 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 () ) )
-> Box
-> ( Int, Box )
-> Cusp
cuspCoords eqs ( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ), i, 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
cuspCoords eqs ( i, 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) )
| StrokeDatum
{ dpath
, dstroke = D22 { _D22_v = stroke } }
@ -1342,7 +1254,6 @@ instance Show RootFindingStep where
data RootFindingLeaf d
= TooSmall d
| NoSolution String d
deriving stock Show
showRootFindingTree :: Box -> RootFindingTree Box -> Tree String
@ -1351,13 +1262,13 @@ showRootFindingTree cand (RootFindingStep s ts)
= Node (show cand ++ " abc " ++ showArea (boxArea cand) ++ " " ++ show s) $ map (\ (c,t) -> showRootFindingTree c t) ts
boxArea :: Box -> Double
boxArea ( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ), _, 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
boxArea ( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) )
= abs ( t_hi - t_lo ) * abs ( s_hi - s_lo )
showArea :: Double -> [Char]
showArea area = "(area " ++ showFFloat (Just 6) area "" ++ ")"
type Box = ( 𝕀ℝ 1, Int, 𝕀ℝ 1 )
type Box = 𝕀ℝ 2
type BoxHistory = [ NE.NonEmpty ( RootFindingStep, Box ) ]
-- | Options for the cusp-finding methods in 'findCusps' and 'findCuspsFrom'.
@ -1371,17 +1282,24 @@ data CuspFindingOptions
}
data RootFindingAlgorithm
-- | Bisection step.
= Bisection !BisectionOptions
-- | Gauss–Seidel step with the given preconditioning method.
= GaussSeidel GaussSeidelOptions
| GaussSeidel !GaussSeidelOptions
-- | @box(1)@-consistency.
| Box1 Box1Options
| Box1 !Box1Options
-- | @box(2)@-consistency.
| Box2 Box2Options
-- | Bisection step.
| Bisection BisectionOptions
| Box2 !Box2Options
-- | Options for the bisection method.
data BisectionOptions =
BisectionOptions
{ canHaveSols :: !( ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 ) -> Box -> Bool )
, fallbackBisectionDim :: !( [ ( RootFindingStep, Box ) ] -> BoxHistory -> ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 ) -> ( Int, String ) )
}
-- | Options for the interval Gauss–Seidel method.
data GaussSeidelOptions = GaussSeidelOptions { gsPreconditioner :: Preconditioner }
data GaussSeidelOptions = GaussSeidelOptions { gsPreconditioner :: !Preconditioner }
-- | Options for the @box(1)@-consistency method.
data Box1Options = Box1Options { box1EpsEq :: !Double }
@ -1389,21 +1307,15 @@ data Box1Options = Box1Options { box1EpsEq :: !Double }
-- | Options for the @box(2)@-consistency method.
data Box2Options = Box2Options { box2EpsEq :: !Double, box2LambdaMin :: !Double }
-- | Options for the bisection method.
data BisectionOptions =
BisectionOptions
{ canHaveSols :: !( ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) -> Box -> Bool )
, fallbackBisectionDim :: !( [ ( RootFindingStep, Box ) ] -> BoxHistory -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) ) -> (Int, String) )
}
defaultCuspFindingOptions :: CuspFindingOptions
defaultCuspFindingOptions =
CuspFindingOptions
{ minWidth
, cuspFindingAlgorithms = defaultCuspFindingAlgorithms minWidth 5e-3
, cuspFindingAlgorithms = defaultCuspFindingAlgorithms minWidth narrowAbs
}
where
minWidth = 1e-5
narrowAbs = 5e-3
defaultCuspFindingAlgorithms :: Double -> Double -> Box -> BoxHistory -> NonEmpty RootFindingAlgorithm
defaultCuspFindingAlgorithms minWidth narrowAbs box history =
@ -1417,17 +1329,18 @@ defaultCuspFindingAlgorithms minWidth narrowAbs box history =
-- ...then try bisecting the box.
-> Bisection ( defaultBisectionOptions minWidth narrowAbs box ) NE.:| []
-- Otherwise, do a normal round: try an interval Gauss–Seidel step followed by box(1)-consistency.
-- Otherwise, do a normal round.
-- Currently: we try an interval Gauss–Seidel step followed by box(1)-consistency.
_ -> GaussSeidel defaultGaussSeidelOptions
NE.:| [ Box1 ( Box1Options { box1EpsEq = narrowAbs } ) ]
where
sufficientlySmallerThan :: Box -> Box -> Bool
sufficientlySmallerThan ( t1, _, s1 ) ( t2, _, s2 ) =
( wd t1 - wd t2 > narrowAbs )
||
( wd s1 - wd s2 > narrowAbs )
wd :: 𝕀ℝ 1 -> Double
wd ( 𝕀 ( ℝ1 lo ) ( ℝ1 hi ) ) = hi - lo
sufficientlySmallerThan
( 𝕀 ( ℝ2 t1_lo s1_lo ) ( ℝ2 t1_hi s1_hi ) )
( 𝕀 ( ℝ2 t2_lo s2_lo ) ( ℝ2 t2_hi s2_hi ) ) =
( ( t1_hi - t1_lo ) - ( t2_hi - t2_lo ) > narrowAbs )
||
( ( s1_hi - s1_lo ) - ( s2_hi - s2_lo ) > narrowAbs )
defaultGaussSeidelOptions :: GaussSeidelOptions
defaultGaussSeidelOptions = GaussSeidelOptions { gsPreconditioner = InverseMidJacobian }
@ -1436,15 +1349,12 @@ defaultBisectionOptions :: Double -> Double -> Box -> BisectionOptions
defaultBisectionOptions
_minWidth
_narrowAbs
( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) )
, i
, s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
)
box@( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) )
= BisectionOptions
{ canHaveSols =
\ eqs _box'@( _t', _i', _s' ) ->
\ eqs box' ->
-- box(0)-consistency
let dat' = ( eqs _t' `Seq.index` _i' ) _s'
let dat' = eqs box'
in ee_potential_zero dat' && 𝛿E𝛿sdcdt_potential_zero dat'
-- box(1)-consistency
@ -1458,7 +1368,7 @@ defaultBisectionOptions
\ _roundHist _prevRoundsHist eqs ->
let StrokeDatum { ee = D22 _ ( T ee_t ) (T ee_s ) _ _ _
, 𝛿E𝛿sdcdt = D12 _ ( T ( T f_t ) ) ( T ( T f_s ) ) }
= ( eqs t `Seq.index` i ) s
= eqs box
wd_t = t_hi - t_lo
wd_s = s_hi - s_lo
tJWidth = wd_t * normVI3 ee_t f_t
@ -1479,15 +1389,19 @@ defaultBisectionOptions
findCusps
:: CuspFindingOptions
-> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-> ( [ ( Box, RootFindingTree Box ) ], ( [ Box ], [ Box ] ) )
findCusps opts eqs = foldMap ( findCuspsFrom opts eqs ) initBoxes
where
initBoxes =
[ ( 𝕀 ( ℝ1 zero ) ( ℝ1 one ), i, 𝕀 ( ℝ1 zero ) ( ℝ1 one ) )
| i <- [ 0 .. length ( eqs ( 𝕀 ( ℝ1 zero ) ( ℝ1 one ) ) ) - 1 ]
]
-> IntMap ( [ ( Box, RootFindingTree Box ) ], ( [ Box ], [ Box ] ) )
findCusps opts eqs =
IntMap.fromList
[ ( i, findCuspsFrom opts eqnPiece initBox )
| i <- [ 0 .. length ( eqs ( 𝕀 ( ℝ1 zero ) ( ℝ1 one ) ) ) - 1 ]
, let initBox = 𝕀 ( ℝ2 zero zero ) ( ℝ2 one one )
eqnPiece ( 𝕀 ( ℝ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 )
in ( eqs t `Seq.index` i ) s
]
-- | Use the following algorithms using interal arithmetic in order
-- | Use the following algorithms using interval arithmetic in order
-- to find cusps in the envelope:
--
-- - interval Newton method with Gauss–Seidel step for inversion
@ -1502,7 +1416,7 @@ findCusps opts eqs = foldMap ( findCuspsFrom opts eqs ) initBoxes
-- boxes which might or might not contain solutions.
findCuspsFrom
:: CuspFindingOptions
-> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-> ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 )
-> Box
-> ( [ ( Box, RootFindingTree Box ) ], ( [ Box ], [ Box ] ) )
findCuspsFrom
@ -1520,26 +1434,23 @@ findCuspsFrom
-> Writer ( [ Box ], [ Box ] )
[ RootFindingTree Box ]
go history
cand@( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) )
, i
, s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
)
cand@( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) )
| -- Check the envelope equation interval contains zero.
not ( ee_potential_zero sd)
||
-- Check the 𝛿E𝛿sdcdt interval contains zero.
not ( 𝛿E𝛿sdcdt_potential_zero sd )
-- Box doesn't contain a solution: discard it.
= return []
-- Box is small: stop processing it.
| t_hi - t_lo < minWidth && s_hi - s_lo < minWidth
= do let res = TooSmall cand
tell ( [ cand ], [] )
return [ RootFindingLeaf res ]
| -- Check the envelope equation interval contains zero.
ee_potential_zero sd
-- Check the 𝛿E𝛿sdcdt interval contains zero.
, 𝛿E𝛿sdcdt_potential_zero sd
= doStrategies history ( cuspFindingAlgorithms cand history ) cand
-- Box doesn't contain a solution: discard it.
| otherwise
, let whyNoSol = if ee_potential_zero sd then "dc/dt" else "ee"
= return [ RootFindingLeaf $ NoSolution whyNoSol cand ]
= doStrategies history ( cuspFindingAlgorithms cand history ) cand
where
sd = ( eqs t `Seq.index` i ) s
sd = eqs $ 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi )
-- Run a round of cusp finding strategies, then recur.
doStrategies
@ -1548,26 +1459,26 @@ findCuspsFrom
-> Box
-> Writer ( [ Box ], [ Box ] )
[ RootFindingTree Box ]
doStrategies hist0 = do_strats []
doStrategies prevRoundsHist = do_strats []
where
do_strats :: [ ( RootFindingStep, Box ) ]
-> NE.NonEmpty RootFindingAlgorithm
-> Box
-> Writer ( [ Box ], [ Box ] )
[ RootFindingTree Box ]
do_strats stratsHist ( algo NE.:| algos ) box = do
do_strats thisRoundHist ( algo NE.:| algos ) box = do
-- Run one strategy in the round.
( step, boxes ) <- doStrategy stratsHist hist0 eqs minWidth algo box
( step, boxes ) <- doStrategy thisRoundHist prevRoundsHist eqs minWidth algo box
case algos of
-- If there are other algorithms to run in this round, run them next.
nextAlgo : otherAlgos ->
let stratsHist' = ( step, box ) : stratsHist
in recur step ( do_strats stratsHist' ( nextAlgo NE.:| otherAlgos ) ) boxes
let thisRoundHist' = ( step, box ) : thisRoundHist
in recur step ( do_strats thisRoundHist' ( nextAlgo NE.:| otherAlgos ) ) boxes
-- Otherwise, we have done one full round of strategies.
-- Recur back to the top (calling 'go').
[] ->
let stratsHist' = ( step, box ) NE.:| stratsHist
history = stratsHist' : hist0
let thisRoundHist' = ( step, box ) NE.:| thisRoundHist
history = thisRoundHist' : prevRoundsHist
in recur step ( go history ) boxes
recur :: RootFindingStep
@ -1583,7 +1494,7 @@ findCuspsFrom
-- (hopefully smaller) output boxes.
doStrategy :: [ ( RootFindingStep, Box ) ]
-> BoxHistory
-> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-> ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 )
-> Double
-> RootFindingAlgorithm
-> Box
@ -1612,24 +1523,17 @@ bisect :: ( Box -> Bool )
-> Box
-> ( [ Box ], ( String, Double ) )
bisect canHaveSols ( dim, why )
( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) )
, i
, s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
)
( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) )
-- We try to bisect along a dimension which eliminates zeros from one of the
-- sub-regions.
--
-- If this fails, we fall-back to the provided dimension choice.
= let t_mid = 0.5 * ( t_lo + t_hi )
s_mid = 0.5 * ( s_lo + s_hi )
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 = 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_mid s_hi )
r = 𝕀 ( ℝ2 t_mid s_lo ) ( ℝ2 t_hi s_hi )
d = 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_mid )
u = 𝕀 ( ℝ2 t_lo s_mid ) ( ℝ2 t_hi s_hi )
l_skip = not $ canHaveSols l
r_skip = not $ canHaveSols r
d_skip = not $ canHaveSols d
@ -1660,24 +1564,21 @@ bisect canHaveSols ( dim, why )
-- | Interval Newton method with Gauss–Seidel step.
intervalGaussSeidel
:: ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
:: ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 )
-- ^ equations
-> Preconditioner
-- ^ preconditioner to use for the Gauss–Seidel step
-> Box
-> Writer ( [ Box ], [ Box ] ) [ Box ]
-> 𝕀ℝ 2
-> Writer ( [ 𝕀ℝ 2 ], [ 𝕀ℝ 2 ] ) [ 𝕀ℝ 2 ]
intervalGaussSeidel
eqs
precondMethod
( t@( 𝕀 ( ℝ1 t_lo ) ( ℝ1 t_hi ) )
, i
, s@( 𝕀 ( ℝ1 s_lo ) ( ℝ1 s_hi ) )
)
ts@( 𝕀 ( ℝ2 t_lo s_lo ) ( ℝ2 t_hi s_hi ) )
| StrokeDatum { 𝛿E𝛿sdcdt = D12 _ ( T ( T f_t ) ) ( T ( T f_s ) ) }
<- ( eqs t `Seq.index` i ) s
<- eqs ts
, 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
<- eqs ts_mid
= let -- Interval Newton method: take one Gauss–Seidel step
-- for the equation f'(x) v = - f(x_mid),
-- where f = 𝛿E/𝛿s * dc/dt
@ -1697,9 +1598,7 @@ intervalGaussSeidel
-- , 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 )
gsGuesses = gaussSeidel a b ( unT $ T ts ^-^ T ts_mid )
in
-- If the Gauss–Seidel step was a contraction, then the box
-- contains a unique solution (by the Banach fixed point theorem).
@ -1713,11 +1612,8 @@ intervalGaussSeidel
where
t_mid = 0.5 * ( t_lo + t_hi )
s_mid = 0.5 * ( s_lo + s_hi )
i_t_mid = singleton ( ℝ1 t_mid )
i_s_mid = singleton ( ℝ1 s_mid )
mkGuess ( t0, s0 ) = ( coerce ( (+) @( 𝕀 Double ) ) t0 i_t_mid
, i
, coerce ( (+) @( 𝕀 Double ) ) s0 i_s_mid )
ts_mid = singleton ( ℝ2 t_mid s_mid )
mkGuess ts_g = unT $ T ts_g ^+^ T ts_mid
neg ( 𝕀 ( ℝ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
@ -1729,7 +1625,7 @@ intervalGaussSeidel
-- See also
-- "Presentation of a highly tuned multithreaded interval solver for underdetermined and well-determined nonlinear systems"
makeBox1Consistent
:: ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
:: ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 )
-> Double -> Double
-> Box -> [ Box ]
makeBox1Consistent eqs minWidth epsEq x =
@ -1739,7 +1635,7 @@ makeBox1Consistent eqs minWidth epsEq x =
-- | An implementation of "bound-consistency" from the paper
-- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
makeBox2Consistent
:: ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
:: ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 )
-> Double -> Double -> Double
-> Box -> Box
makeBox2Consistent eqs minWidth epsEq lambdaMin x0 = ( `State.evalState` False ) $ doLoop 0.25 x0
@ -1951,7 +1847,7 @@ pipeFunctionsWhileTrue fns = go fns
allNarrowingOperators
:: Double
-> Double
-> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 3 𝕀 ) )
-> ( 𝕀ℝ 2 -> StrokeDatum 3 𝕀 )
-> [ Box -> State Bool [ Box ] ]
allNarrowingOperators eps_eq eps_bis eqs =
[ \ cand ->
@ -1968,41 +1864,42 @@ allNarrowingOperators eps_eq eps_bis eqs =
[ ( ( get_s, set_s ), ff' ) | ff' <- [ ff'_s, g1g1'_s, g2g2'_s ] ]
]
where
ff'_t (t, i, s) =
ff'_t ts =
let D22 { _D22_v = 𝕀 ( ℝ1 ee_inf ) ( ℝ1 ee_sup )
, _D22_dx = T ( 𝕀 ( ℝ1 ee'_inf ) ( ℝ1 ee'_sup ) )
} = ee $ ( eqs t `Seq.index` i ) s
} = ee $ eqs ts
in ( 𝕀 ee_inf ee_sup, 𝕀 ee'_inf ee'_sup )
ff'_s (t, i, s) =
ff'_s ts =
let D22 { _D22_v = 𝕀 ( ℝ1 ee_inf ) ( ℝ1 ee_sup )
, _D22_dy = T ( 𝕀 ( ℝ1 ee'_inf ) ( ℝ1 ee'_sup ) )
} = ee $ ( eqs t `Seq.index` i ) s
} = ee $ eqs ts
in ( 𝕀 ee_inf ee_sup, 𝕀 ee'_inf ee'_sup )
g1g1'_t (t, i, s) =
g1g1'_t ts =
let D12 { _D12_v = T ( 𝕀 ( ℝ2 x_inf _ ) ( ℝ2 x_sup _ ) )
, _D12_dx = T ( T ( 𝕀 ( ℝ2 x'_inf _ ) ( ℝ2 x'_sup _ ) ) )
} = 𝛿E𝛿sdcdt $ ( eqs t `Seq.index` i ) s
} = 𝛿E𝛿sdcdt $ eqs ts
in ( 𝕀 x_inf x_sup, 𝕀 x'_inf x'_sup )
g1g1'_s (t, i, s) =
g1g1'_s ts =
let D12 { _D12_v = T ( 𝕀 ( ℝ2 x_inf _ ) ( ℝ2 x_sup _ ) )
, _D12_dy = T ( T ( 𝕀 ( ℝ2 x'_inf _ ) ( ℝ2 x'_sup _ ) ) )
} = 𝛿E𝛿sdcdt $ ( eqs t `Seq.index` i ) s
} = 𝛿E𝛿sdcdt $ eqs ts
in ( 𝕀 x_inf x_sup, 𝕀 x'_inf x'_sup )
g2g2'_t (t, i, s) =
g2g2'_t ts =
let D12 { _D12_v = T ( 𝕀 ( ℝ2 _ y_inf ) ( ℝ2 _ y_sup ) )
, _D12_dx = T ( T ( 𝕀 ( ℝ2 _ y'_inf ) ( ℝ2 _ y'_sup ) ) )
} = 𝛿E𝛿sdcdt $ ( eqs t `Seq.index` i ) s
} = 𝛿E𝛿sdcdt $ eqs ts
in ( 𝕀 y_inf y_sup, 𝕀 y'_inf y'_sup )
g2g2'_s (t, i, s) =
g2g2'_s ts =
let D12 { _D12_v = T ( 𝕀 ( ℝ2 _ y_inf ) ( ℝ2 _ y_sup ) )
, _D12_dy = T ( T ( 𝕀 ( ℝ2 _ y'_inf ) ( ℝ2 _ y'_sup ) ) )
} = 𝛿E𝛿sdcdt $ ( eqs t `Seq.index` i ) s
} = 𝛿E𝛿sdcdt $ eqs ts
in ( 𝕀 y_inf y_sup, 𝕀 y'_inf y'_sup )
-- TODO: replace this with Fin n and representable please
get_t, get_s :: Box -> 𝕀 Double
get_t ( 𝕀 ( ℝ1 t_inf ) ( ℝ1 t_sup ), _, _ ) = 𝕀 t_inf t_sup
get_t ( 𝕀 ( ℝ2 t_inf _ ) ( ℝ2 t_sup _ ) ) = 𝕀 t_inf t_sup
get_s ( 𝕀 ( ℝ2 _ s_inf ) ( ℝ2 _ s_sup ) ) = 𝕀 s_inf s_sup
set_t, set_s :: 𝕀 Double -> Box -> Box
set_t ( 𝕀 t_inf t_sup ) ( _, i, s ) = ( 𝕀 ( ℝ1 t_inf ) ( ℝ1 t_sup ), i, s )
get_s ( _ , _, 𝕀 ( ℝ1 s_inf ) ( ℝ1 s_sup ) ) = 𝕀 s_inf s_sup
set_s ( 𝕀 s_inf s_sup ) ( t, i, _ ) = ( t, i, 𝕀 ( ℝ1 s_inf ) ( ℝ1 s_sup ) )
set_t ( 𝕀 t_inf t_sup ) ( 𝕀 ( ℝ2 _ s_inf ) ( ℝ2 _ s_sup ) ) = 𝕀 ( ℝ2 t_inf s_inf ) ( ℝ2 t_sup s_sup )
set_s ( 𝕀 s_inf s_sup ) ( 𝕀 ( ℝ2 t_inf _ ) ( ℝ2 t_sup _ ) ) = 𝕀 ( ℝ2 t_inf s_inf ) ( ℝ2 t_sup s_sup )