From 43098dec0124af354e2093d4040d1a2acb23a59b Mon Sep 17 00:00:00 2001 From: sheaf Date: Mon, 9 Jan 2023 23:59:00 +0100 Subject: [PATCH 1/2] Newton-Raphson implementation --- src/splines/Math/Roots.hs | 122 +++++++++++++++++++++++++++++++++++++- 1 file changed, 121 insertions(+), 1 deletion(-) diff --git a/src/splines/Math/Roots.hs b/src/splines/Math/Roots.hs index 10f0ac5..c3fd8bb 100644 --- a/src/splines/Math/Roots.hs +++ b/src/splines/Math/Roots.hs @@ -1,4 +1,7 @@ -{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE DuplicateRecordFields #-} +{-# LANGUAGE ScopedTypeVariables #-} + +{-# OPTIONS_GHC -Wno-name-shadowing #-} module Math.Roots where @@ -274,3 +277,120 @@ derivative p = do go ( i + 1 ) go 0 pure p' + +-------------------------------------------------------------------------------- + +data NewtonRaphson + = NewtonRaphson + { f :: Double -> (# Double, Double #) + , maxIters :: Word + , min_x, max_x :: Double + , digits :: Int } + +-- Newton–Raphson implementation taken from Boost C++ library: +-- https://github.com/boostorg/math/blob/0dc6a70caa6bbec2b6ae25eede36c430f0ccae13/include/boost/math/tools/roots.hpp#L217 + +newtonRaphson :: NewtonRaphson + -> Double + -> Maybe Double +newtonRaphson ( NewtonRaphson { f, maxIters, min_x, max_x, digits } ) x0 = + go $ + NewtonRaphsonState + { f_x_prev = 0 + , x = x0 + , δ = maxDouble, δ1 = maxDouble + , iters = 0 + , min_x, max_x + , min_f_x = 0, max_f_x = 0 } + where + factor = encodeFloat 1 ( 1 - digits ) + go ( NewtonRaphsonState { f_x_prev, x, δ = δ1, δ1 = δ2, iters, min_x, max_x, min_f_x, max_f_x } ) + = case f x of + (# f_x, f'_x #) + | f_x == 0 + -> Just x + | δ <- + if | f'_x == 0 + -> handleZeroDerivative f_x_prev x f_x δ1 min_x max_x + | otherwise + -> f_x / f'_x + , (# δ, δ1 #) <- + if | abs ( δ * δ ) > abs δ2 + , let shift = if δ > 0 then 0.5 * ( x - min_x ) else 0.5 * ( x - max_x ) + , δ <- if x /= 0 && ( abs shift > abs x ) + then signum δ * abs x * 1.1 + else shift + -> (# δ, 3 * δ #) + | otherwise + -> (# δ, δ1 #) + , let new_x = x - δ + , (# δ, new_x #) <- + if | new_x <= min_x + , δ <- 0.5 * ( x - min_x ) + -> (# δ, x - δ #) + | new_x >= max_x + , δ <- 0.5 * ( x - max_x ) + -> (# δ, x - δ #) + | otherwise + -> (# δ, new_x #) + -> if + | abs δ <= abs ( new_x * factor ) + || new_x == min_x || new_x == max_x + -> Just x + | iters >= maxIters + -> Nothing + | (# min_x, max_x, min_f_x, max_f_x #) <- + if δ > 0 + then (# min_x, x, min_f_x, f_x #) + else (# x, max_x, f_x, max_f_x #) + -> if min_f_x * max_f_x > 0 + then Nothing + else + go $ NewtonRaphsonState + { f_x_prev = f_x, x = new_x + , δ, δ1 + , iters = iters + 1 + , min_x, max_x + , min_f_x, max_f_x + } + + handleZeroDerivative :: Double + -> Double -> Double + -> Double + -> Double -> Double + -> Double + handleZeroDerivative f_x_prev x f_x δ min_x max_x + -- Handle zero derivative on first iteration. + | f_x_prev == 0 + , x_prev <- if x <= 0.5 * ( min_x + max_x ) then max_x else min_x + , (# f_x_prev, _ #) <- f x_prev + , δ <- x_prev - x + = finish f_x_prev δ + | otherwise + = finish f_x_prev δ + + where + finish f_x_prev δ + | signum f_x_prev * signum f_x < 0 + = if δ < 0 then 0.5 * ( x - min_x ) else 0.5 * ( x - max_x ) + | otherwise + = if δ < 0 then 0.5 * ( x - max_x ) else 0.5 * ( x - min_x ) + +-- | Loop state for the 'newtonRaphson' function. +data NewtonRaphsonState = + NewtonRaphsonState + { f_x_prev :: !Double + , x :: !Double + , δ, δ1 :: !Double + , iters :: !Word + , min_x, max_x :: !Double + , min_f_x, max_f_x :: !Double } + +maxDouble :: Double +maxDouble = encodeFloat m n + where + b = floatRadix ( 0 :: Double ) + e = floatDigits ( 0 :: Double ) + (_, e') = floatRange ( 0 :: Double ) + m = b ^ e - 1 + n = e' - e From 8fc5e6c9b8e86e053bd77947a41172a7a27c4d49 Mon Sep 17 00:00:00 2001 From: sheaf Date: Mon, 9 Jan 2023 20:54:19 +0100 Subject: [PATCH 2/2] compute brush strokes using the envelope equation --- MetaBrush.cabal | 4 + app/Main.hs | 8 +- src/app/MetaBrush/Action.hs | 2 +- src/app/MetaBrush/Application.hs | 19 +- src/app/MetaBrush/Render/Document.hs | 13 +- src/app/MetaBrush/UI/ToolBar.hs | 4 +- src/metabrushes/MetaBrush/Asset/Brushes.hs | 41 +- src/metabrushes/MetaBrush/Brush.hs | 12 +- .../MetaBrush/Document/SubdivideStroke.hs | 4 +- src/metabrushes/MetaBrush/Records.hs | 35 +- src/splines/Debug/Utils.hs | 26 + src/splines/Math/Bezier/Cubic/Fit.hs | 14 +- src/splines/Math/Bezier/Stroke.hs | 531 +++++++++++------- src/splines/Math/Linear.hs | 10 +- src/splines/Math/Linear/Dual.hs | 162 +++--- src/splines/Math/Module.hs | 10 + src/splines/Math/Roots.hs | 248 ++++---- 17 files changed, 683 insertions(+), 460 deletions(-) create mode 100644 src/splines/Debug/Utils.hs diff --git a/MetaBrush.cabal b/MetaBrush.cabal index 7d1f5f5..a156589 100644 --- a/MetaBrush.cabal +++ b/MetaBrush.cabal @@ -36,6 +36,8 @@ common common >= 4.17 && < 4.18 , acts ^>= 0.3.1.0 + , code-page + ^>= 0.2.1 , containers >= 0.6.0.1 && < 0.7 , deepseq @@ -112,6 +114,7 @@ common common common extras build-depends: + directory >= 1.3.4.0 && < 1.4 , filepath @@ -180,6 +183,7 @@ library splines , Math.Module , Math.Orientation , Math.Roots + , Debug.Utils build-depends: bifunctors diff --git a/app/Main.hs b/app/Main.hs index ad2335f..dd18056 100644 --- a/app/Main.hs +++ b/app/Main.hs @@ -1,6 +1,4 @@ -{-# LANGUAGE BlockArguments #-} {-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE TypeApplications #-} module Main ( main ) @@ -14,6 +12,10 @@ import System.Exit import GHC.Conc ( getNumProcessors, setNumCapabilities ) +-- code-page +import System.IO.CodePage + ( withCP65001 ) + -- gi-gio import qualified GI.Gio as GIO @@ -30,7 +32,7 @@ import MetaBrush.Application -------------------------------------------------------------------------------- main :: IO () -main = do +main = withCP65001 do procs <- getNumProcessors let diff --git a/src/app/MetaBrush/Action.hs b/src/app/MetaBrush/Action.hs index 9668caf..735e931 100644 --- a/src/app/MetaBrush/Action.hs +++ b/src/app/MetaBrush/Action.hs @@ -821,7 +821,7 @@ instance HandleAction MouseClick where changeText :: Text changeText = "Subdivide " <> loc pure ( UpdateDoc . UpdateDocumentTo $ HistoryChange {..} ) - + -- Ignore double click event otherwise. _ -> pure Don'tModifyDoc diff --git a/src/app/MetaBrush/Application.hs b/src/app/MetaBrush/Application.hs index 2954617..57db0eb 100644 --- a/src/app/MetaBrush/Application.hs +++ b/src/app/MetaBrush/Application.hs @@ -84,7 +84,6 @@ import Math.Linear import MetaBrush.Action ( ActionOrigin(..) ) import qualified MetaBrush.Asset.Brushes as Asset.Brushes - ( EllipseBrushFields, ellipse ) import MetaBrush.Asset.Colours ( getColours ) import MetaBrush.Asset.Logo @@ -161,6 +160,12 @@ runApplication application = do , strokeUnique = strokeUnique , strokeBrush = Just Asset.Brushes.ellipse , strokeSpline = +-- Spline +-- { splineStart = mkPoint ( ℝ2 -20 -20 ) 5 +-- , splineCurves = OpenCurves $ Seq.fromList +-- [ LineTo { curveEnd = NextPoint ( mkPoint ( ℝ2 20 20 ) 5 ), curveData = invalidateCache undefined } +-- ] +-- } Spline { splineStart = mkPoint ( ℝ2 10 -20 ) 2 1 0 , splineCurves = OpenCurves $ Seq.fromList @@ -176,6 +181,8 @@ 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 -> PointData ( Record Asset.Brushes.CircleBrushFields ) +-- mkPoint pt r = PointData pt Normal ( MkR $ ℝ1 r ) recomputeStrokesTVar <- STM.newTVarIO @Bool False documentRenderTVar <- STM.newTVarIO @( ( Int32, Int32 ) -> Cairo.Render () ) ( const $ pure () ) @@ -193,11 +200,11 @@ runApplication application = do maxHistorySizeTVar <- STM.newTVarIO @Int 1000 fitParametersTVar <- STM.newTVarIO @FitParameters ( FitParameters - { maxSubdiv = 6 - , nbSegments = 12 - , dist_tol = 5e-3 - , t_tol = 1e-4 - , maxIters = 100 + { maxSubdiv = 2 --3 -- 6 + , nbSegments = 3 --6 -- 12 + , dist_tol = 0.1 -- 5e-3 + , t_tol = 0.1 -- 1e-4 + , maxIters = 2 -- 100 } ) diff --git a/src/app/MetaBrush/Render/Document.hs b/src/app/MetaBrush/Render/Document.hs index c7762ee..8f98239 100644 --- a/src/app/MetaBrush/Render/Document.hs +++ b/src/app/MetaBrush/Render/Document.hs @@ -292,17 +292,18 @@ strokeRenderData fitParams { inject , project = toUsedParams :: Record pointFields -> Record usedFields } -> do - let embedUsedParams r = inject r brush_defaults + let embedUsedParams = inject brush_defaults -- Compute the outline using the brush function. ( outline, fitPts ) <- - computeStrokeOutline @( T ( Record usedFields) ) @clo - fitParams ( toUsedParams . brushParams ) ( fun brushFn . embedUsedParams ) spline + computeStrokeOutline @clo fitParams + ( toUsedParams . brushParams ) embedUsedParams brushFn + spline pure $ StrokeWithOutlineRenderData { strokeDataSpline = spline , strokeOutlineData = ( outline, fitPts ) - , strokeBrushFunction = fun brushFn . embedUsedParams . toUsedParams + , strokeBrushFunction = fun brushFn . fun embedUsedParams . toUsedParams } _ -> pure $ StrokeRenderData @@ -604,7 +605,7 @@ drawFitPoint _ zoom ( FitPoint { fitPoint = ℝ2 x y } ) = do lift do Cairo.save Cairo.translate x y - Cairo.arc 0 0 ( 2 / zoom ) 0 ( 2 * pi ) + Cairo.arc 0 0 ( 6 / zoom ) 0 ( 2 * pi ) -- 2 / zoom Cairo.setSourceRGBA r g b 0.8 Cairo.fill Cairo.restore @@ -621,7 +622,7 @@ drawFitPoint _ zoom ( FitTangent { fitPoint = ℝ2 x y, fitTangent = V2 tx ty } Cairo.translate x y Cairo.moveTo 0 0 Cairo.lineTo ( 0.05 * tx ) ( 0.05 * ty ) - Cairo.setLineWidth ( 1 / zoom ) + Cairo.setLineWidth ( 4 / zoom ) -- 1 / zoom Cairo.setSourceRGBA r g b 0.8 Cairo.stroke Cairo.arc 0 0 ( 2 / zoom ) 0 ( 2 * pi ) diff --git a/src/app/MetaBrush/UI/ToolBar.hs b/src/app/MetaBrush/UI/ToolBar.hs index 0a60fbb..d1bbeba 100644 --- a/src/app/MetaBrush/UI/ToolBar.hs +++ b/src/app/MetaBrush/UI/ToolBar.hs @@ -55,12 +55,12 @@ data Mode data ToolBar = ToolBar { selectionTool, penTool, pathTool, brushTool, metaTool, debugTool - :: !GTK.ToggleButton + :: !GTK.ToggleButton } createToolBar :: Variables -> Colours -> GTK.Box -> IO ToolBar createToolBar ( Variables {..} ) colours toolBar = do - + widgetAddClass toolBar "toolBar" GTK.widgetSetValign toolBar GTK.AlignStart diff --git a/src/metabrushes/MetaBrush/Asset/Brushes.hs b/src/metabrushes/MetaBrush/Asset/Brushes.hs index a283370..130fe58 100644 --- a/src/metabrushes/MetaBrush/Asset/Brushes.hs +++ b/src/metabrushes/MetaBrush/Asset/Brushes.hs @@ -3,10 +3,6 @@ module MetaBrush.Asset.Brushes where --- base -import Data.Coerce - ( coerce ) - -- containers import qualified Data.Sequence as Seq ( fromList ) @@ -24,9 +20,9 @@ import qualified Data.HashMap.Strict as HashMap -- MetaBrush import Math.Bezier.Spline import Math.Linear - ( ℝ(..), T(..) ) + ( ℝ(..), Fin(..) ) import Math.Linear.Dual - ( D, type (~>)(..), Var(var), konst ) + ( D, type (~>)(..), var, konst ) import Math.Module ( Module((^+^), (*^)) ) import MetaBrush.Brush @@ -86,17 +82,16 @@ circleBrush :: Record CircleBrushFields ~> Spline 'Closed () ( ℝ 2 ) circleBrush = D \ params -> let r :: D ( Record CircleBrushFields ) Double - r = runD ( var @1 ) params + r = runD ( var ( Fin 1## ) ) params mkPt :: Double -> Double -> D ( Record CircleBrushFields ) ( ℝ 2 ) mkPt ( kon -> x ) ( kon -> y ) - = fmap coerce - $ ( x * r ) *^ e_x + = ( x * r ) *^ e_x ^+^ ( y * r ) *^ e_y in circleSpline @( Record CircleBrushFields ) mkPt where - e_x, e_y :: D ( Record CircleBrushFields ) ( T ( ℝ 2 ) ) - e_x = pure $ T $ ℝ2 1 0 - e_y = pure $ T $ ℝ2 0 1 + e_x, e_y :: D ( Record CircleBrushFields ) ( ℝ 2 ) + e_x = pure $ ℝ2 1 0 + e_y = pure $ ℝ2 0 1 kon = konst @( Record CircleBrushFields ) @@ -104,25 +99,17 @@ ellipseBrush :: Record EllipseBrushFields ~> Spline 'Closed () ( ℝ 2 ) ellipseBrush = D \ params -> let a, b, phi :: D ( Record EllipseBrushFields ) Double - a = runD ( var @1 ) params - b = runD ( var @2 ) params - phi = runD ( var @3 ) params + a = runD ( var ( Fin 1## ) ) params + b = runD ( var ( Fin 2## ) ) params + phi = runD ( var ( Fin 3## ) ) params mkPt :: Double -> Double -> D ( Record EllipseBrushFields ) ( ℝ 2 ) mkPt ( kon -> x ) ( kon -> y ) - = fmap coerce - $ ( x * a * cos phi - y * b * sin phi ) *^ e_x + = ( x * a * cos phi - y * b * sin phi ) *^ e_x ^+^ ( y * b * cos phi + x * a * sin phi ) *^ e_y in circleSpline @( Record EllipseBrushFields ) mkPt where - e_x, e_y :: D ( Record EllipseBrushFields ) ( T ( ℝ 2 ) ) - e_x = pure $ T $ ℝ2 1 0 - e_y = pure $ T $ ℝ2 0 1 + e_x, e_y :: D ( Record EllipseBrushFields ) ( ℝ 2 ) + e_x = pure $ ℝ2 1 0 + e_y = pure $ ℝ2 0 1 kon = konst @( Record EllipseBrushFields ) - ---ellipseArc :: ℝ 2 ~> ℝ 2 ---ellipseArc = brushStroke ( linear myPath ) ( uncurryD $ fmap bezier3 myBrush ) - ---testing :: Double -> Double -> (# Double, T ( ℝ 2 ) #) ---testing :: Double -> Double -> (# Double, T (ℝ 2) #) ---testing t s = envelopeEquation ellipseArc t s diff --git a/src/metabrushes/MetaBrush/Brush.hs b/src/metabrushes/MetaBrush/Brush.hs index a3fdec0..e32fca4 100644 --- a/src/metabrushes/MetaBrush/Brush.hs +++ b/src/metabrushes/MetaBrush/Brush.hs @@ -1,6 +1,7 @@ -{-# LANGUAGE AllowAmbiguousTypes #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE QuantifiedConstraints #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE UndecidableInstances #-} module MetaBrush.Brush ( Brush(..), SomeBrush(..), BrushFunction @@ -38,12 +39,17 @@ import qualified Data.Text as Text -- MetaBrush import Math.Linear + ( ℝ, Representable ) import Math.Linear.Dual ( Diffy ) +import Math.Module + ( Interpolatable ) import Math.Bezier.Spline ( SplineType(Closed), SplinePts) import MetaBrush.Records + ( KnownSymbols, Length, Record, WithParams ) import MetaBrush.Serialisable + ( Serialisable ) -------------------------------------------------------------------------------- diff --git a/src/metabrushes/MetaBrush/Document/SubdivideStroke.hs b/src/metabrushes/MetaBrush/Document/SubdivideStroke.hs index e17b851..fdf5016 100644 --- a/src/metabrushes/MetaBrush/Document/SubdivideStroke.hs +++ b/src/metabrushes/MetaBrush/Document/SubdivideStroke.hs @@ -49,7 +49,7 @@ import Math.Bezier.Spline import Math.Bezier.Stroke ( CachedStroke(..), invalidateCache ) import Math.Module - ( lerp, quadrance, closestPointOnSegment ) + ( Interpolatable, lerp, quadrance, closestPointOnSegment ) import Math.Linear ( Segment(..), ℝ(..), T(..) ) import MetaBrush.Document @@ -57,8 +57,6 @@ import MetaBrush.Document , PointData(..), DiffPointData(..) , coords, _strokeSpline ) -import MetaBrush.Records - ( Interpolatable ) -------------------------------------------------------------------------------- diff --git a/src/metabrushes/MetaBrush/Records.hs b/src/metabrushes/MetaBrush/Records.hs index 250bba8..c80c117 100644 --- a/src/metabrushes/MetaBrush/Records.hs +++ b/src/metabrushes/MetaBrush/Records.hs @@ -54,13 +54,6 @@ import Math.Module -------------------------------------------------------------------------------- --- | A convenient constraint synonym for types that support interpolation. -type Interpolatable :: Type -> Constraint -class ( Torsor ( T r ) r, Module Double ( T r ) ) => Interpolatable r -instance ( Torsor ( T r ) r, Module Double ( T r ) ) => Interpolatable r - --------------------------------------------------------------------------------- - -- | A function from a given record type, with provided default values -- that can be overridden. type WithParams :: [ Symbol ] -> Type -> Type @@ -123,9 +116,11 @@ instance ( Torsor ( T ( ℝ ( Length ks ) ) ) ( ℝ ( Length ks ) ) => Torsor ( T ( Record ks ) ) ( Record ks ) where MkR g --> MkR a = T $ MkR $ unT $ g --> a -type instance D ( Record ks ) = D ( ℝ ( Length ks ) ) +deriving newtype + instance Representable ( ℝ ( Length ks ) ) + => Representable ( Record ks ) -deriving newtype instance Var n ( ℝ ( Length ks ) ) => Var n ( Record ks ) +type instance D ( Record ks ) = D ( ℝ ( Length ks ) ) deriving newtype instance Diffy ( ℝ ( Length ks ) ) => Diffy ( Record ks ) -------------------------------------------------------------------------------- @@ -154,31 +149,35 @@ intersect :: forall r1 r2 l1 l2 , KnownSymbols r1, KnownSymbols r2 , l1 ~ Length r1, l2 ~ Length r2 , Representable ( ℝ l1 ), Representable ( ℝ l2 ) - , Interpolatable ( Record r1 ) + , Interpolatable ( Record r1 ), Diffy ( ℝ l2 ) ) => Intersection r1 r2 intersect -- Shortcut when the two rows are equal. | Just Refl <- eqT @r1 @r2 , Refl <- ( unsafeCoerce Refl :: r1 :~: Intersect r1 r2 ) - = Intersection { project = id, inject = const } + = Intersection { project = id, inject = \ _ -> linear id } | otherwise = doIntersection @r1 @r2 \ ( _ :: Proxy# r1r2 ) r1_idxs r2_idxs -> let project :: Record r1 -> Record r1r2 project = \ ( MkR r1 ) -> MkR $ projection ( (!) r1_idxs ) r1 - inject :: Record r1r2 -> Record r2 -> Record r2 - inject = \ ( MkR r1r2 ) ( MkR r2 ) -> MkR $ injection ( find eqFin r2_idxs ) r1r2 r2 + inject :: Record r2 -> Record r1r2 ~> Record r2 + inject = \ ( MkR r2 ) -> linear \ ( MkR r1r2 ) -> MkR $ injection ( find eqFin r2_idxs ) r1r2 r2 in Intersection { project, inject } data Intersection r1 r2 where Intersection - :: forall r1r2 r1 r2 - . ( KnownSymbols r1r2, Representable ( ℝ ( Length r1r2 ) ) + :: forall r1r2 r1 r2 l12 + . ( l12 ~ Length r1r2 + , KnownSymbols r1r2 + , Representable ( ℝ l12 ) + , Diffy ( ℝ l12 ) , Interpolatable ( Record r1r2 ) ) - => { project :: Record r1 -> Record r1r2 - , inject :: Record r1r2 -> Record r2 -> Record r2 + => { project :: Record r1 -> Record r1r2 + , inject :: Record r2 -> Record r1r2 ~> Record r2 + -- ^ overrides the components of the first record with the second } -> Intersection r1 r2 {-# INLINE doIntersection #-} @@ -190,7 +189,7 @@ doIntersection ) => ( forall r1r2 l12. ( r1r2 ~ Intersect r1 r2, l12 ~ Length r1r2 - , Representable ( ℝ l12 ), Interpolatable ( ℝ l12 ) + , Representable ( ℝ l12 ), Diffy ( ℝ l12 ), Interpolatable ( ℝ l12 ) , KnownSymbols r1r2, Representable ( ℝ ( Length r1r2 ) ) ) => Proxy# r1r2 -> Vec l12 ( Fin l1 ) -> Vec l12 ( Fin l2 ) -> kont ) diff --git a/src/splines/Debug/Utils.hs b/src/splines/Debug/Utils.hs new file mode 100644 index 0000000..eb85aec --- /dev/null +++ b/src/splines/Debug/Utils.hs @@ -0,0 +1,26 @@ +module Debug.Utils ( trace ) where + + +-- base +import System.IO + ( BufferMode(..), hSetBuffering, hFlush, hPutStrLn, stdout ) +import System.IO.Unsafe + ( unsafePerformIO ) + +-- deepseq +import Control.DeepSeq + ( force ) + +-- code-page +import System.IO.CodePage + ( withCP65001 ) + +-------------------------------------------------------------------------------- + +trace :: String -> a -> a +trace ( force -> !str ) expr = + unsafePerformIO $ withCP65001 do + hSetBuffering stdout ( BlockBuffering $ Just 50 ) + hPutStrLn stdout str + hFlush stdout + return expr diff --git a/src/splines/Math/Bezier/Cubic/Fit.hs b/src/splines/Math/Bezier/Cubic/Fit.hs index 2969289..3283647 100644 --- a/src/splines/Math/Bezier/Cubic/Fit.hs +++ b/src/splines/Math/Bezier/Cubic/Fit.hs @@ -120,7 +120,7 @@ data FitPoint -- including the meaning of \( \texttt{t_tol} \) and \( \texttt{maxIters} \). fitSpline :: FitParameters - -> ( Double -> ( ℝ 2, T ( ℝ 2 ) ) ) -- ^ curve \( t \mapsto ( C(t), C'(t) ) \) to fit + -> ( ℝ 1 -> ( ℝ 2, T ( ℝ 2 ) ) ) -- ^ curve \( t \mapsto ( C(t), C'(t) ) \) to fit -> ( SplinePts Open, Seq FitPoint ) fitSpline ( FitParameters {..} ) = go 0 where @@ -128,16 +128,16 @@ fitSpline ( FitParameters {..} ) = go 0 dt = recip ( fromIntegral nbSegments ) go :: Int - -> ( Double -> ( ℝ 2, T ( ℝ 2 ) ) ) + -> ( ℝ 1 -> ( ℝ 2, T ( ℝ 2 ) ) ) -> ( SplinePts Open, Seq FitPoint ) go subdiv curveFn = let p, r :: ℝ 2 tp, tr :: T ( ℝ 2 ) qs :: [ ℝ 2 ] - (p, tp) = curveFn 0 - (r, tr) = curveFn 1 - qs = [ fst $ curveFn ( dt * fromIntegral j ) | j <- [ 1 .. nbSegments - 1 ] ] + (p, tp) = curveFn $ ℝ1 0 + (r, tr) = curveFn $ ℝ1 1 + qs = [ fst $ curveFn ( ℝ1 $ dt * fromIntegral j ) | j <- [ 1 .. nbSegments - 1 ] ] in case fitPiece dist_tol t_tol maxIters p tp qs r tr of ( bez, Max ( Arg max_sq_error t_split ) ) @@ -150,8 +150,8 @@ fitSpline ( FitParameters {..} ) = go 0 c1, c2 :: SplinePts Open ps1, ps2 :: Seq FitPoint ( ( c1, ps1 ), ( c2, ps2 ) ) - = ( go ( subdiv + 1 ) ( \ t -> curveFn $ t * t_split_eff ) - , go ( subdiv + 1 ) ( \ t -> curveFn $ t_split_eff + t * ( 1 - t_split_eff ) ) + = ( go ( subdiv + 1 ) ( \ ( ℝ1 t ) -> curveFn $ ℝ1 $ t * t_split_eff ) + , go ( subdiv + 1 ) ( \ ( ℝ1 t ) -> curveFn $ ℝ1 $ t_split_eff + t * ( 1 - t_split_eff ) ) ) `Parallel.Strategy.using` ( Parallel.Strategy.parTuple2 Parallel.Strategy.rdeepseq Parallel.Strategy.rdeepseq ) -> ( c1 <> c2, ps1 <> ps2 ) diff --git a/src/splines/Math/Bezier/Stroke.hs b/src/splines/Math/Bezier/Stroke.hs index d0b70cc..6fd5020 100644 --- a/src/splines/Math/Bezier/Stroke.hs +++ b/src/splines/Math/Bezier/Stroke.hs @@ -1,5 +1,6 @@ -{-# LANGUAGE AllowAmbiguousTypes #-} -{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE QuantifiedConstraints #-} +{-# LANGUAGE ScopedTypeVariables #-} module Math.Bezier.Stroke ( Offset(..) @@ -9,10 +10,8 @@ module Math.Bezier.Stroke -- * Brush stroking - -- $brushes , brushStroke, envelopeEquation - , linear, bezier2, bezier3 --- , uncurryD + , line, bezier2, bezier3 ) where @@ -22,6 +21,8 @@ import Prelude hiding ( unzip ) import Control.Arrow ( first, (***) ) +import Control.Applicative + ( Applicative(..) ) import Control.Monad ( unless ) import Control.Monad.ST @@ -30,6 +31,8 @@ import Data.Bifunctor ( Bifunctor(bimap) ) import Data.Foldable ( for_ ) +import Data.Functor.Identity + ( Identity(..) ) import Data.List.NonEmpty ( unzip ) import Data.Maybe @@ -53,11 +56,15 @@ import Data.Act import Data.Sequence ( Seq(..) ) import qualified Data.Sequence as Seq + --( empty, reverse, singleton, index ) +import Data.Set + ( Set ) +import qualified Data.Set as Set ( singleton ) -- deepseq import Control.DeepSeq - ( NFData(..), NFData1, deepseq ) + ( NFData(..), NFData1, deepseq, force ) -- generic-lens import Data.Generics.Product.Typed @@ -65,10 +72,6 @@ import Data.Generics.Product.Typed import Data.Generics.Internal.VL ( set, view ) --- groups -import Data.Group - ( Group ) - -- parallel import qualified Control.Parallel.Strategies as Strats ( rdeepseq, parTuple2, using ) @@ -101,8 +104,8 @@ import qualified Math.Bezier.Quadratic as Quadratic import Math.Epsilon ( epsilon ) import Math.Module - ( Module(..), Inner((^.^)) - , lerp, squaredNorm, cross + ( Module(..), Inner((^.^)), Interpolatable + , lerp, cross , convexCombination, strictlyParallel ) import Math.Orientation @@ -110,10 +113,12 @@ import Math.Orientation , between ) import Math.Roots - ( solveQuadratic ) + ( solveQuadratic, newtonRaphson ) import Math.Linear import Math.Linear.Dual +--import Debug.Utils + -------------------------------------------------------------------------------- data Offset @@ -168,42 +173,42 @@ coords = view typed -------------------------------------------------------------------------------- +type OutlineFn = ℝ 1 -> ( ( ℝ 2, T ( ℝ 2 ) ), ( ℝ 2, T ( ℝ 2 ) ) ) + computeStrokeOutline :: - forall diffParams ( clo :: SplineType ) brushParams crvData ptData s + forall ( clo :: SplineType ) usedParams brushParams crvData ptData s . ( KnownSplineType clo - , Group diffParams - , Module Double diffParams - , Torsor diffParams brushParams + , Interpolatable usedParams + , Diffy usedParams, Diffy brushParams , HasType ( ℝ 2 ) ptData , HasType ( CachedStroke s ) crvData , NFData ptData, NFData crvData -- Debugging. - , Show ptData + , Show ptData, Show brushParams ) => FitParameters - -> ( ptData -> brushParams ) - -> ( brushParams -> SplinePts Closed ) + -> ( ptData -> usedParams ) + -> ( usedParams ~> brushParams ) + -> ( brushParams ~> SplinePts Closed ) -> Spline clo crvData ptData -> ST s ( Either ( SplinePts Closed ) ( SplinePts Closed, SplinePts Closed ) , Seq FitPoint ) -computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart = spt0 } ) = case ssplineType @clo of +computeStrokeOutline fitParams ptParams toBrushParams brushFn spline@( Spline { splineStart = spt0 } ) = case ssplineType @clo of -- Open brush path with at least one segment. -- Need to add caps at both ends of the path. SOpen | firstCurve :<| _ <- openCurves $ splineCurves spline , prevCurves :|> lastCurve <- openCurves $ splineCurves spline - , ( firstOutlineFwd, firstOutlineBwd ) :<| _ <- outlineFns - , _ :|> ( lastOutlineFwd, lastOutlineBwd ) <- outlineFns + , firstOutlineFn :<| _ <- outlineFns + , _ :|> lastOutlineFn <- outlineFns , let endPt :: ptData endPt = openCurveEnd lastCurve startTgtFwd, startTgtBwd, endTgtFwd, endTgtBwd :: T ( ℝ 2 ) - startTgtFwd = snd ( firstOutlineFwd 0 ) - startTgtBwd = -1 *^ snd ( firstOutlineBwd 1 ) - endTgtFwd = snd ( lastOutlineFwd 1 ) - endTgtBwd = -1 *^ snd ( lastOutlineBwd 0 ) + ( ( _, startTgtFwd), ( _, startTgtBwd ) ) = firstOutlineFn $ ℝ1 0 + ( ( _, endTgtFwd ), ( _, endTgtBwd ) ) = lastOutlineFn $ ℝ1 1 startBrush, endBrush :: SplinePts Closed startBrush = brushShape spt0 endBrush = brushShape endPt @@ -252,8 +257,8 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart = -- Add forward and backward caps at the start. SClosed | ClosedCurves prevCurves lastCurve <- splineCurves spline - , ( firstOutlineFwd, firstOutlineBwd ) :<| _ <- outlineFns - , _ :|> ( lastOutlineFwd, lastOutlineBwd ) <- outlineFns + , firstOutlineFn :<| _ <- outlineFns + , _ :|> lastOutlineFn <- outlineFns , let startTgt, endTgt, startTgtFwd, startTgtBwd, endTgtFwd, endTgtBwd :: T ( ℝ 2 ) startTgt = case prevCurves of @@ -262,10 +267,8 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart = endTgt = case prevCurves of Empty -> endTangent spt0 spt0 lastCurve _ :|> lastPrev -> endTangent spt0 ( openCurveEnd lastPrev ) lastCurve - startTgtFwd = snd ( firstOutlineFwd 0 ) - startTgtBwd = -1 *^ snd ( firstOutlineBwd 1 ) - endTgtFwd = snd ( lastOutlineFwd 1 ) - endTgtBwd = -1 *^ snd ( lastOutlineBwd 0 ) + ( ( _, startTgtFwd), ( _, startTgtBwd ) ) = firstOutlineFn $ ℝ1 0 + ( ( _, endTgtFwd ), ( _, endTgtBwd ) ) = lastOutlineFn $ ℝ1 1 fwdStartCap, bwdStartCap :: SplinePts Open TwoSided fwdStartCap bwdStartCap = fmap fst . snd . runWriter @@ -284,26 +287,19 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart = ) where - outlineFns - :: Seq - ( Double -> ( ℝ 2, T ( ℝ 2 ) ) - , Double -> ( ℝ 2, T ( ℝ 2 ) ) - ) + outlineFns :: Seq OutlineFn outlineFns = go spt0 ( openCurves $ splineCurves ( adjustSplineType @Open spline ) ) where go :: ptData -> Seq ( Curve Open crvData ptData ) - -> Seq - ( Double -> ( ℝ 2, T ( ℝ 2 ) ) - , Double -> ( ℝ 2, T ( ℝ 2 ) ) - ) + -> Seq OutlineFn go _ Empty = Empty go p0 ( crv :<| crvs ) = - outlineFunctions @diffParams ptParams brushFn p0 crv :<| go ( openCurveEnd crv ) crvs + outlineFunction ptParams toBrushParams brushFn p0 crv :<| go ( openCurveEnd crv ) crvs brushShape :: ptData -> SplinePts Closed - brushShape pt = brushFn ( ptParams pt ) + brushShape pt = fun brushFn $ fun toBrushParams $ ptParams pt updateSpline :: ( T ( ℝ 2 ), T ( ℝ 2 ), T ( ℝ 2 ) ) -> ST s OutlineData updateSpline ( lastTgt, lastTgtFwd, lastTgtBwd ) @@ -313,17 +309,15 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart = ( \ ptData curve -> do ( prevTgt, prev_tgtFwd, prev_tgtBwd ) <- get let - fwd, bwd :: Double -> ( ℝ 2, T ( ℝ 2 ) ) - ( fwd, bwd ) = outlineFunctions @diffParams ptParams brushFn ptData curve + fwdBwd :: OutlineFn + fwdBwd = outlineFunction ptParams toBrushParams brushFn ptData curve tgt, next_tgt, tgtFwd, next_tgtFwd, tgtBwd, next_tgtBwd :: T ( ℝ 2 ) - tgt = startTangent spt0 ptData curve - next_tgt = endTangent spt0 ptData curve - tgtFwd = snd ( fwd 0 ) - next_tgtFwd = snd ( fwd 1 ) - tgtBwd = -1 *^ snd ( bwd 1 ) - next_tgtBwd = -1 *^ snd ( bwd 0 ) + tgt = startTangent spt0 ptData curve + next_tgt = endTangent spt0 ptData curve + ( ( _, tgtFwd ), ( _, tgtBwd ) ) = fwdBwd $ ℝ1 0 + ( ( _, next_tgtFwd ), ( _, next_tgtBwd ) ) = fwdBwd $ ℝ1 1 lift $ tellBrushJoin ( prevTgt, prev_tgtFwd, tgtBwd ) ptData ( tgt, tgtFwd, prev_tgtBwd ) - lift $ updateCurveData ( curveData curve ) fwd bwd + lift $ updateCurveData ( curveData curve ) fwdBwd put ( next_tgt, next_tgtFwd, next_tgtBwd ) ) ( const ( pure () ) ) @@ -331,10 +325,9 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart = updateCurveData :: crvData - -> ( Double -> ( ℝ 2, T ( ℝ 2 ) ) ) - -> ( Double -> ( ℝ 2, T ( ℝ 2 ) ) ) + -> OutlineFn -> WriterT OutlineData ( ST s ) () - updateCurveData ( view ( typed @( CachedStroke s ) ) -> CachedStroke { cachedStrokeRef } ) fwd bwd = do + updateCurveData ( view ( typed @( CachedStroke s ) ) -> CachedStroke { cachedStrokeRef } ) fwdBwd = do mbOutline <- lift ( readSTRef cachedStrokeRef ) case mbOutline of -- Cached fit data is available: use it. @@ -346,12 +339,12 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart = let fwdData, bwdData :: ( SplinePts Open, Seq FitPoint ) ( fwdData, bwdData ) = - ( fitSpline fitParams fwd, fitSpline fitParams bwd ) + ( fitSpline fitParams ( fst . fwdBwd ), fitSpline fitParams ( snd . fwdBwd ) ) `Strats.using` ( Strats.parTuple2 Strats.rdeepseq Strats.rdeepseq ) outlineData :: OutlineData - outlineData = TwoSided fwdData bwdData - outlineData `deepseq` tell ( outlineData ) + outlineData = TwoSided fwdData ( bimap reverseSpline Seq.reverse bwdData ) + outlineData `deepseq` tell outlineData lift $ writeSTRef cachedStrokeRef ( Just outlineData ) -- Connecting paths at a point of discontinuity of the tangent vector direction (G1 discontinuity). @@ -416,107 +409,55 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart = $ joinWithBrush brush0 tgtBwd prevTgtBwd -- | Computes the forward and backward stroke outline functions for a single curve. -outlineFunctions - :: forall diffParams brushParams crvData ptData - . ( Group diffParams, Module Double diffParams - , Torsor diffParams brushParams +outlineFunction + :: forall usedParams brushParams crvData ptData + . ( Interpolatable usedParams + , Diffy usedParams, Diffy brushParams , HasType ( ℝ 2 ) ptData -- Debugging. - , Show ptData + , Show ptData, Show brushParams ) - => ( ptData -> brushParams ) - -> ( brushParams -> SplinePts Closed ) + => ( ptData -> usedParams ) + -> ( usedParams ~> brushParams ) + -> ( brushParams ~> SplinePts Closed ) -> ptData -> Curve Open crvData ptData - -> ( Double -> ( ℝ 2, T ( ℝ 2 ) ) - , Double -> ( ℝ 2, T ( ℝ 2 ) ) - ) -outlineFunctions ptParams brushFn sp0 crv = + -> OutlineFn +outlineFunction ptParams toBrushParams brushFromParams sp0 crv = let - p0 :: ℝ 2 - p0 = coords sp0 - brush :: Double -> SplinePts Closed - f :: Double -> ℝ 2 - f' :: Double -> T ( ℝ 2 ) - ( brush, f, f' ) = case crv of - LineTo { curveEnd = NextPoint sp1 } - | let - p1 :: ℝ 2 - p1 = coords sp1 - tgt :: T ( ℝ 2 ) - tgt = p0 --> p1 - brush1 :: Double -> SplinePts Closed - brush1 t = brushFn ( lerp @diffParams t ( ptParams sp0 ) ( ptParams sp1 ) ) - -> ( brush1, \ t -> lerp @( T ( ℝ 2 ) ) t p0 p1, const tgt ) - Bezier2To { controlPoint = sp1, curveEnd = NextPoint sp2 } - | let - p1, p2 :: ℝ 2 - p1 = coords sp1 - p2 = coords sp2 - bez :: Quadratic.Bezier ( ℝ 2 ) - bez = Quadratic.Bezier {..} - brush2 :: Double -> SplinePts Closed - brush2 t = - brushFn $ - Quadratic.bezier @diffParams - ( Quadratic.Bezier ( ptParams sp0 ) ( ptParams sp1 ) ( ptParams sp2 ) ) t - -> ( brush2, Quadratic.bezier @( T ( ℝ 2 ) ) bez, Quadratic.bezier' bez ) - Bezier3To { controlPoint1 = sp1, controlPoint2 = sp2, curveEnd = NextPoint sp3 } - | let - p1, p2, p3 :: ℝ 2 - p1 = coords sp1 - p2 = coords sp2 - p3 = coords sp3 - bez :: Cubic.Bezier ( ℝ 2 ) - bez = Cubic.Bezier {..} - brush3 :: Double -> SplinePts Closed - brush3 t = - brushFn $ - Cubic.bezier @diffParams - ( Cubic.Bezier ( ptParams sp0 ) ( ptParams sp1 ) ( ptParams sp2 ) ( ptParams sp3 ) ) t - -> ( brush3, Cubic.bezier @( T ( ℝ 2 ) ) bez, Cubic.bezier' bez ) - fwd, bwd :: Double -> ( ℝ 2, T ( ℝ 2 ) ) - fwd t - = ( off t --offset ( withTangent ( fwd' t ) ( brush t ) ) • f t - , fwd' t - ) + usedParams :: ℝ 1 ~> usedParams + path :: ℝ 1 ~> ℝ 2 + ( path, usedParams ) = + case crv of + LineTo { curveEnd = NextPoint sp1 } + | let seg = Segment sp0 sp1 + -> ( line ( fmap coords seg ), line ( fmap ptParams seg ) ) + Bezier2To { controlPoint = sp1, curveEnd = NextPoint sp2 } + | let bez2 = Quadratic.Bezier sp0 sp1 sp2 + -> ( bezier2 ( fmap coords bez2 ), bezier2 ( fmap ptParams bez2 ) ) + Bezier3To { controlPoint1 = sp1, controlPoint2 = sp2, curveEnd = NextPoint sp3 } + | let bez3 = Cubic.Bezier sp0 sp1 sp2 sp3 + -> ( bezier3 ( fmap coords bez3 ), bezier3 ( fmap ptParams bez3 ) ) + + fwdBwd :: OutlineFn + fwdBwd t + = solveEnvelopeEquations path_t path'_t ( fwdOffset, bwdOffset ) curves + -- = ( ( offset fwdOffset • path_t, path'_t ) + -- , ( offset bwdOffset • path_t, -1 *^ path'_t ) ) where - off :: Double -> ℝ 2 - off u = offset ( withTangent ( f' u ) ( brush u ) ) • f u - offTgt :: Double -> T ( ℝ 2 ) - offTgt u - | u < 0.5 - = 1e9 *^ ( off u --> off (u + 1e-9) ) - | otherwise - = 1e9 *^ ( off (u - 1e-9) --> off u ) - fwd' :: Double -> T ( ℝ 2 ) - fwd' u - | squaredNorm ( offTgt u ) < epsilon - = f' u - | otherwise - = offTgt u - bwd t - = ( off s --offset ( withTangent ( -1 *^ bwd' s ) ( brush s ) ) • f s - , bwd' s - ) - where - s :: Double - s = 1 - t - off :: Double -> ℝ 2 - off u = offset ( withTangent ( -1 *^ f' u ) ( brush u ) ) • f u - offTgt :: Double -> T ( ℝ 2 ) - offTgt u - | u < 0.5 - = 1e9 *^ ( off u --> off (u + 1e-9) ) - | otherwise - = 1e9 *^ ( off (u - 1e-9) --> off u ) - bwd' :: Double -> T ( ℝ 2 ) - bwd' u - | squaredNorm ( offTgt u ) < epsilon - = -1 *^ f' u - | otherwise - = offTgt u - in ( fwd, bwd ) + + curves :: Seq ( ℝ 1 -> StrokeDatum ) + curves = brushStrokeData path ( usedParams `chainRule` toBrushParams ) brushFromParams t + + fwdOffset = withTangent path'_t brush_t + bwdOffset = withTangent ( -1 *^ path'_t ) brush_t + + D1 path_t path'_t _ = runD path t + D1 params_t _ _ = runD usedParams t + brush_t = value @brushParams + $ runD brushFromParams ( fun toBrushParams params_t ) + + in fwdBwd ----------------------------------- -- Various utility functions @@ -616,7 +557,7 @@ discardCurveData = bimap ( const () ) coords dropFirstPiece :: HasType ( ℝ 2 ) ptData => Spline Open crvData ptData -> Maybe ( SplinePts Open ) dropFirstPiece ( Spline { splineCurves = OpenCurves curves } ) = case curves of Empty -> Nothing - fstPiece :<| laterPieces -> + fstPiece :<| laterPieces -> Just $ Spline { splineStart = coords ( openCurveEnd fstPiece ) , splineCurves = OpenCurves $ fmap discardCurveData laterPieces @@ -799,38 +740,6 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) -------------------------------------------------------------------------------- --- $brushes --- --- You can compute the envelope equation for a brush stroke by using --- the functions 'linear', 'bezier2' and 'bezier3' in conjunction with --- the 'brushStroke' function, e.g. --- --- > brushStroke ( bezier2 path ) ( uncurryD $ fmap bezier3 brush ) - --- | Linear interpolation, as a differentiable function. -linear :: forall b. ( Module Double ( T b ), Torsor ( T b ) b ) - => Segment b -> ℝ 1 ~> b -linear ( Segment a b ) = D \ ( ℝ1 t ) -> - D1 ( lerp @( T b ) t a b ) - ( a --> b ) - origin - --- | A quadratic Bézier curve, as a differentiable function. -bezier2 :: forall b. ( Module Double ( T b ), Torsor ( T b ) b ) - => Quadratic.Bezier b -> ℝ 1 ~> b -bezier2 bez = D \ ( ℝ1 t ) -> - D1 ( Quadratic.bezier @( T b ) bez t ) - ( Quadratic.bezier' bez t ) - ( Quadratic.bezier'' bez ) - --- | A cubic Bézier curve, as a differentiable function. -bezier3 :: forall b. ( Module Double ( T b ), Torsor ( T b ) b ) - => Cubic.Bezier b -> ℝ 1 ~> b -bezier3 bez = D \ ( ℝ1 t ) -> - D1 ( Cubic.bezier @( T b ) bez t ) - ( Cubic.bezier' bez t ) - ( Cubic.bezier'' bez t ) - -- | A brush stroke, as described by the equation -- -- \[ c(t,s) = p(t) + b(t,s) \] @@ -839,15 +748,10 @@ bezier3 bez = D \ ( ℝ1 t ) -> -- -- - \( p(t) \) is the path that the brush follows, and -- - \( b(t,s) \) is the brush shape, as it varies along the path. -brushStroke :: ℝ 1 ~> ℝ 2 -- ^ stroke path \( p(t) \) - -> ℝ 2 ~> ℝ 2 -- ^ brush \( b(t,s) \) - -> ℝ 2 ~> ℝ 2 -brushStroke ( D f_p ) ( D f_b ) = D \ ( ℝ2 t0 s0 ) -> - let !( D1 p dpdt d2pdt2 ) - = f_p ( ℝ1 t0 ) - !( D2 b dbdt dbds d2bdt2 d2bdtds d2bds2 ) - = f_b ( ℝ2 t0 s0 ) - in +brushStroke :: D ( ℝ 1 ) ( ℝ 2 ) -- ^ stroke path \( p(t) \) + -> D ( ℝ 2 ) ( ℝ 2 ) -- ^ brush \( b(t,s) \) + -> D ( ℝ 2 ) ( ℝ 2 ) +brushStroke ( D1 p dpdt d2pdt2 ) ( D2 b dbdt dbds d2bdt2 d2bdtds d2bds2 ) = D2 ( unT $ T p ^+^ T b ) -- c = p + b @@ -868,14 +772,16 @@ brushStroke ( D f_p ) ( D f_b ) = D \ ( ℝ2 t0 s0 ) -> -- -- \[ \frac{\partial E}{\partial s} \frac{\mathrm{d} c}{\mathrm{d} t} \] -- --- whose roots correspond to cusps in the envelope. -envelopeEquation :: ℝ 2 ~> ℝ 2 -> Double -> Double -> (# Double, T ( ℝ 2 ) #) -envelopeEquation ( D c ) t s = - case c ( ℝ2 t s ) of - D2 _ dcdt dcds d2cdt2 d2cdtds d2cds2 -> - let dEdt = d2cdt2 `cross2` dcds + dcdt `cross2` d2cdtds - dEds = d2cdtds `cross2` dcds + dcdt `cross2` d2cds2 - in (# dcdt `cross2` dcds, dEds *^ dcdt ^-^ dEdt *^ dcds #) +-- whose roots correspond to cusps in the envelope, and +-- +-- \[ \frac{\partial E}{\partial s}. \] +envelopeEquation :: D ( ℝ 2 ) ( ℝ 2 ) -> ( Double, T ( ℝ 2 ), Double, Double ) +envelopeEquation ( D2 _ dcdt dcds d2cdt2 d2cdtds d2cds2 ) = + let ee = dcdt `cross` dcds + dEdt = d2cdt2 `cross` dcds + dcdt `cross` d2cdtds + dEds = d2cdtds `cross` dcds + dcdt `cross` d2cds2 + tot = dcdt ^-^ ( dEdt / dEds ) *^ dcds --dEds *^ dcdt ^-^ dEdt *^ dcds + in ( ee, tot, dEdt, dEds ) -- Computation of total derivative dc/dt: -- -- dc/dt = ∂c/∂t + ∂c/∂s ∂s/∂t @@ -883,7 +789,218 @@ envelopeEquation ( D c ) t s = -- -- ( ∂E / ∂s ) dc/dt = ( ∂E / ∂s ) ∂c/∂t - ( ∂E / ∂t ) ∂c/∂s. --- | Cross-product of two 2D vectors. -cross2 :: T ( ℝ 2 ) -> T ( ℝ 2 ) -> Double -cross2 ( T ( ℝ2 x1 y1 ) ) ( T ( ℝ2 x2 y2 ) ) - = x1 * y2 - x2 * y1 + +-- | Linear interpolation, as a differentiable function. +line :: forall b. ( Module Double ( T b ), Torsor ( T b ) b ) + => Segment b -> ℝ 1 ~> b +line ( Segment a b ) = D \ ( ℝ1 t ) -> + D1 ( lerp @( T b ) t a b ) + ( a --> b ) + origin + +-- | A quadratic Bézier curve, as a differentiable function. +bezier2 :: forall b. ( Module Double ( T b ), Torsor ( T b ) b ) + => Quadratic.Bezier b -> ℝ 1 ~> b +bezier2 bez = D \ ( ℝ1 t ) -> + D1 ( Quadratic.bezier @( T b ) bez t ) + ( Quadratic.bezier' bez t ) + ( Quadratic.bezier'' bez ) + +-- | A cubic Bézier curve, as a differentiable function. +bezier3 :: forall b. ( Module Double ( T b ), Torsor ( T b ) b ) + => Cubic.Bezier b -> ℝ 1 ~> b +bezier3 bez = D \ ( ℝ1 t ) -> + D1 ( Cubic.bezier @( T b ) bez t ) + ( Cubic.bezier' bez t ) + ( Cubic.bezier'' bez t ) + +splineCurveFns :: SplinePts Closed -> Seq ( ℝ 1 ~> ℝ 2 ) +splineCurveFns spls + = runIdentity + . bifoldSpline + ( \ pt crv -> Identity . Seq.singleton $ curveFn pt crv ) + ( const $ Identity Seq.empty ) + . adjustSplineType @Open + $ spls + where + curveFn :: ℝ 2 -> Curve Open () ( ℝ 2 ) -> ( ℝ 1 ~> ℝ 2 ) + curveFn p0 = \case + LineTo { curveEnd = NextPoint p1 } + -> line $ Segment p0 p1 + Bezier2To { controlPoint = p1, curveEnd = NextPoint p2 } + -> bezier2 $ Quadratic.Bezier p0 p1 p2 + Bezier3To { controlPoint1 = p1, controlPoint2 = p2, curveEnd = NextPoint p3 } + -> bezier3 $ Cubic.Bezier p0 p1 p2 p3 + +-- | 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 :: ℝ 2 + -> T ( ℝ 2 ) + -> ( Offset, Offset ) + -> Seq ( ℝ 1 -> StrokeDatum ) + -> ( ( ℝ 2, T ( ℝ 2 ) ), ( ℝ 2, T ( ℝ 2 ) ) ) +solveEnvelopeEquations path_t path'_t ( fwdOffset, bwdOffset ) strokeData + = ( fwdSol, ( bwdPt, -1 *^ bwdTgt ) ) + + where + + -- !_ = trace + -- ( unlines + -- [ "solveEnvelopeEquation" + -- , " pt: " ++ show path_t + -- , " tgt: " ++ show path'_t + -- , " fwdOffset: " ++ show fwdOffset + -- , " bwdOffset: " ++ show bwdOffset ] + -- ) True + + fwdSol = findSolFrom False fwdOffset + ( bwdPt, bwdTgt ) = findSolFrom True bwdOffset + + n :: Int + n = length strokeData + + findSolFrom :: Bool -> Offset -> ( ℝ 2, T ( ℝ 2 ) ) + findSolFrom goingBwds ( Offset { offsetIndex = i00, offsetParameter = s00, offset = off } ) + = go ( Set.singleton i00 ) i00 ( fromMaybe 0.5 s00 ) + where + + plausibleTangent :: T ( ℝ 2 ) -> Bool + plausibleTangent tgt + | goingBwds + = path'_t ^.^ tgt < 0 + | otherwise + = path'_t ^.^ tgt > 0 + + go :: Set Int -> Int -> Double -> ( ℝ 2, T ( ℝ 2 ) ) + go seen i0 s0 = + case sol s0 ( strokeData `Seq.index` i0 ) of + ( good, s, pt, tgt ) + | True --good && plausibleTangent tgt + -> ( pt, tgt ) + -- | let ( i', s0' ) = mbNextPoint i0 ( unℝ1 s ) + -- , not ( i' `Set.member` seen ) + -- -> go ( Set.insert i' seen ) i' s0' + | otherwise + -> ( off • path_t, path'_t ) + + mbNextPoint :: Int -> Double -> ( Int, Double ) + mbNextPoint i0 s0 + | s0 <= 0.5 + = ( prev i0, 0.9 ) + | otherwise + = ( next i0, 0.1 ) + + prev, next :: Int -> Int + prev i + | i == 0 + = n - 1 + | otherwise + = i - 1 + next i + | i == n - 1 + = 0 + | otherwise + = i + 1 + + sol :: Double -> ( ℝ 1 -> StrokeDatum ) -> ( Bool, ℝ 1, ℝ 2, T ( ℝ 2 ) ) + sol initialGuess f = + let (good, s) = case newtonRaphson maxIters precision domain ( eqn f ) initialGuess of + Nothing -> ( False, initialGuess ) + Just s0 -> ( True , s0 ) + in case f ( ℝ1 s ) of -- TODO: a bit redundant to have to compute this again... + StrokeDatum { dstroke, dcdt } -> + ( good, ℝ1 s, value @( ℝ 2 ) dstroke, dcdt ) + + + eqn :: ( ℝ 1 -> StrokeDatum ) -> ( Double -> ( Double, Double ) ) + eqn f s = + case f ( ℝ1 s ) of + StrokeDatum { ee, 𝛿E𝛿s } -> + ( ee, 𝛿E𝛿s ) + + maxIters :: Word + maxIters = 5 --30 + precision :: Int + precision = 10 + domain :: ( Double, Double ) + domain = ( 0, 1 ) + +newtype ZipSeq a = ZipSeq { getZipSeq :: Seq a } + deriving stock Functor +instance Applicative ZipSeq where + pure _ = error "no pure for ZipSeq" + liftA2 f ( ZipSeq xs ) ( ZipSeq ys ) = ZipSeq ( Seq.zipWith f xs ys ) + +brushStrokeData :: forall brushParams + . ( Diffy brushParams, Show brushParams ) + => ( ℝ 1 ~> ℝ 2 ) -- ^ path + -> ( ℝ 1 ~> brushParams ) -- ^ brush parameters + -> ( brushParams ~> SplinePts Closed ) -- ^ brush from parameters + -> ( ℝ 1 -> Seq ( ℝ 1 -> StrokeDatum ) ) +brushStrokeData path params brush = + \ t -> + let + dpath_t :: D ( ℝ 1 ) ( ℝ 2 ) + !dpath_t = runD path t + dparams_t :: D ( ℝ 1 ) brushParams + !dparams_t@( D1 { v = params_t } ) = runD params t + dbrush_params :: D brushParams ( SplinePts Closed ) + !dbrush_params = runD brush params_t + splines :: Seq ( D brushParams ( ℝ 1 ~> ℝ 2 ) ) + !splines = getZipSeq $ traverse @_ @ZipSeq ( ZipSeq . splineCurveFns ) dbrush_params + dbrushes_t :: Seq ( ℝ 1 -> D ( ℝ 2 ) ( ℝ 2 ) ) + !dbrushes_t = force $ fmap ( uncurryD' . ( dparams_t `chain` ) ) splines + + in fmap ( mkStrokeDatum dpath_t ) dbrushes_t + where + + mkStrokeDatum :: D ( ℝ 1 ) ( ℝ 2 ) + -> ( ℝ 1 -> D ( ℝ 2 ) ( ℝ 2 ) ) + -> ( ℝ 1 -> StrokeDatum ) + mkStrokeDatum dpath_t dbrush_t s = + let dbrush_t_s = dbrush_t s + dstroke@( D2 _c _𝛿c𝛿t _𝛿c𝛿s _ _ _ ) = brushStroke dpath_t dbrush_t_s + ( ee, dcdt, _𝛿E𝛿t, 𝛿E𝛿s ) = envelopeEquation dstroke + in -- trace + -- ( unlines + -- [ "envelopeEquation:" + -- , " s = " ++ show s + -- , " c = " ++ show _c + -- , " ∂c/∂t = " ++ show _𝛿c𝛿t + -- , " ∂c/∂s = " ++ show _𝛿c𝛿s + -- , " E = " ++ show ee + -- , " ∂E/∂t = " ++ show _𝛿E𝛿t + -- , " ∂E/∂s = " ++ show 𝛿E𝛿s + -- , " dc/dt = " ++ show dcdt ] ) $ + StrokeDatum + { dpath = dpath_t + , dbrush = dbrush_t_s + , dstroke + , ee, dcdt, 𝛿E𝛿s } + +-- | The value and derivative of a brush stroke at a given coordinate +-- \( (t_0, s_0) \), together with the value of the envelope equation at that +-- point. +data StrokeDatum + = StrokeDatum + { -- | Path \( p(t_0) \) (with its 1st and 2nd derivatives). + dpath :: D ( ℝ 1 ) ( ℝ 2 ) + -- | Brush shape \( b(t_0, s_0) \) (with its 1st and 2nd derivatives). + , dbrush :: D ( ℝ 2 ) ( ℝ 2 ) + + -- Everything below can be computed in terms of the first two fields. + + -- | Stroke \( c(t_0,s_0) = p(t_0) + b(t_0,s_0) \) (with its 1st and 2nd derivatives). + , dstroke :: D ( ℝ 2 ) ( ℝ 2 ) + -- | Envelope + -- + -- \[ E(t_0,s_0) = \left ( \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} \right )_{(t_0,s_0)}. \] + , ee :: Double + -- | Total derivative + -- + -- \[ \left ( \frac{\mathrm{d} c}{\mathrm{d} t} \right )_{(t_0,s_0)}. \] + , dcdt :: T ( ℝ 2 ) + -- | \( \left ( \frac{\partial E}{\partial s} \right )_{(t_0,s_0)}. \) + , 𝛿E𝛿s :: Double + } + deriving stock Show diff --git a/src/splines/Math/Linear.hs b/src/splines/Math/Linear.hs index cf771d0..5f1c62b 100644 --- a/src/splines/Math/Linear.hs +++ b/src/splines/Math/Linear.hs @@ -71,7 +71,7 @@ data family ℝ n data instance ℝ 0 = ℝ0 deriving stock ( Show, Eq, Ord, Generic ) deriving anyclass NFData -newtype instance ℝ 1 = ℝ1 Double +newtype instance ℝ 1 = ℝ1 { unℝ1 :: Double } deriving stock ( Generic ) deriving newtype ( Show, Eq, Ord, NFData ) data instance ℝ 2 = ℝ2 {-# UNPACK #-} !Double {-# UNPACK #-} !Double @@ -86,8 +86,12 @@ data instance ℝ 3 = ℝ3 {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-# -- | Tangent space to Euclidean space. type T :: Type -> Type newtype T e = T { unT :: e } - deriving stock ( Eq, Functor ) - deriving newtype ( Show, NFData ) -- newtype Show instance for debugging... + deriving stock ( Eq, Functor, Foldable, Traversable ) + deriving newtype NFData -- newtype Show instance for debugging... + +instance {-# OVERLAPPING #-} Show ( ℝ n ) => Show ( T ( ℝ n ) ) where + show ( T p ) = "V" ++ drop 1 ( show p ) +deriving stock instance Show v => Show ( T v ) instance Applicative T where pure = T diff --git a/src/splines/Math/Linear/Dual.hs b/src/splines/Math/Linear/Dual.hs index 149d846..94d7c47 100644 --- a/src/splines/Math/Linear/Dual.hs +++ b/src/splines/Math/Linear/Dual.hs @@ -1,5 +1,6 @@ {-# LANGUAGE AllowAmbiguousTypes #-} {-# LANGUAGE DuplicateRecordFields #-} +{-# LANGUAGE QuantifiedConstraints #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE UndecidableInstances #-} @@ -8,12 +9,12 @@ module Math.Linear.Dual where -- base import Control.Applicative ( liftA2 ) +import Data.Coerce + ( coerce ) import Data.Kind ( Type, Constraint ) import GHC.Generics - ( Generic, Generic1, Generically1(..) ) -import GHC.TypeNats - ( Nat ) + ( Generic, Generic1(..), Generically1(..) ) -- MetaBrush import Math.Module @@ -24,8 +25,14 @@ import Math.Linear -- | Differentiable mappings between spaces. infixr 0 ~> type (~>) :: Type -> Type -> Type -newtype a ~> b = D { runD :: a -> D a b } -deriving stock instance Functor ( D a ) => Functor ( (~>) a ) +newtype u ~> v = D { runD :: u -> D u v } +deriving stock instance Functor ( D u ) => Functor ( (~>) u ) + +instance ( Applicative ( D u ), Module ( D u Double ) ( D u v ) ) => Module Double ( T ( u ~> v ) ) where + origin = T $ D \ _ -> origin + T ( D f ) ^+^ T ( D g ) = T $ D \ t -> f t ^+^ g t + T ( D f ) ^-^ T ( D g ) = T $ D \ t -> f t ^-^ g t + a *^ T ( D f ) = T $ D \ t -> pure a *^ f t -- | @D ( ℝ n ) v@ is \( \mathbb{R}[x_1, \ldots, x_n]/(x_1, \ldots, x_n)^3 \otimes_\mathbb{R} v \) type D :: Type -> Type -> Type @@ -36,22 +43,25 @@ type instance D ( ℝ 2 ) = Dℝ2 type instance D ( ℝ 3 ) = Dℝ3 newtype Dℝ0 v = D0 { v :: v } - deriving stock ( Show, Eq, Functor, Generic, Generic1 ) + deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) deriving newtype ( Num, Fractional, Floating ) deriving Applicative via Generically1 Dℝ0 -data Dℝ1 v = D1 { v :: !v, dx :: !( T v ), ddx :: !( T v ) } - deriving stock ( Show, Eq, Functor, Generic, Generic1 ) +data Dℝ1 v = D1 { v :: !v, df :: !( T v ), d²f :: !( T v ) } + deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 ) deriving Applicative via Generically1 Dℝ1 +deriving stock instance ( Show v, Show ( T v ) ) => Show ( Dℝ1 v ) data Dℝ2 v = D2 { v :: !v, dx, dy :: !( T v ), ddx, dxdy, ddy :: !( T v ) } - deriving stock ( Show, Eq, Functor, Generic, Generic1 ) + deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 ) deriving Applicative via Generically1 Dℝ2 +deriving stock instance ( Show v, Show ( T v ) ) => Show ( Dℝ2 v ) data Dℝ3 v = D3 { v :: !v, dx, dy, dz :: !( T v ), ddx, dxdy, ddy, dxdz, dydz, ddz :: !( T v ) } - deriving stock ( Show, Eq, Functor, Generic, Generic1 ) + deriving stock ( Eq, Functor, Foldable, Traversable, Generic, Generic1 ) deriving Applicative via Generically1 Dℝ3 +deriving stock instance ( Show v, Show ( T v ) ) => Show ( Dℝ3 v ) instance Num ( Dℝ1 Double ) where (+) = liftA2 (+) @@ -110,29 +120,29 @@ instance Num ( Dℝ3 Double ) where ( T $ dz1 * dz2 + v1 * ddz2 + ddz1 * v2) -instance Module Double v => Module ( Dℝ0 Double ) ( Dℝ0 v ) where - (^+^) = liftA2 (^+^) - (^-^) = liftA2 (^-^) - origin = pure origin - (*^) = liftA2 (*^) +instance Module Double ( T v ) => Module ( Dℝ0 Double ) ( Dℝ0 v ) where + (^+^) = liftA2 ( coerce $ (^+^) @Double @( T v ) ) + (^-^) = liftA2 ( coerce $ (^-^) @Double @( T v ) ) + origin = pure ( coerce $ origin @Double @( T v ) ) + (*^) = liftA2 ( coerce $ (*^) @Double @( T v ) ) -instance Module Double v => Module ( Dℝ1 Double ) ( Dℝ1 v ) where - (^+^) = liftA2 (^+^) - (^-^) = liftA2 (^-^) - origin = pure origin - (*^) = liftA2 (*^) +instance Module Double ( T v ) => Module ( Dℝ1 Double ) ( Dℝ1 v ) where + (^+^) = liftA2 ( coerce $ (^+^) @Double @( T v ) ) + (^-^) = liftA2 ( coerce $ (^-^) @Double @( T v ) ) + origin = pure ( coerce $ origin @Double @( T v ) ) + (*^) = liftA2 ( coerce $ (*^) @Double @( T v ) ) -instance Module Double v => Module ( Dℝ2 Double ) ( Dℝ2 v ) where - (^+^) = liftA2 (^+^) - (^-^) = liftA2 (^-^) - origin = pure origin - (*^) = liftA2 (*^) +instance Module Double ( T v ) => Module ( Dℝ2 Double ) ( Dℝ2 v ) where + (^+^) = liftA2 ( coerce $ (^+^) @Double @( T v ) ) + (^-^) = liftA2 ( coerce $ (^-^) @Double @( T v ) ) + origin = pure ( coerce $ origin @Double @( T v ) ) + (*^) = liftA2 ( coerce $ (*^) @Double @( T v ) ) -instance Module Double v => Module ( Dℝ3 Double ) ( Dℝ3 v ) where - (^+^) = liftA2 (^+^) - (^-^) = liftA2 (^-^) - origin = pure origin - (*^) = liftA2 (*^) +instance Module Double ( T v ) => Module ( Dℝ3 Double ) ( Dℝ3 v ) where + (^+^) = liftA2 ( coerce $ (^+^) @Double @( T v ) ) + (^-^) = liftA2 ( coerce $ (^-^) @Double @( T v ) ) + origin = pure ( coerce $ origin @Double @( T v ) ) + (*^) = liftA2 ( coerce $ (*^) @Double @( T v ) ) instance Fractional ( Dℝ1 Double ) where (/) = error "I haven't yet defined (/) for Dℝ1" @@ -204,48 +214,67 @@ instance Floating ( Dℝ3 Double ) where -------------------------------------------------------------------------------- uncurryD :: ( ℝ 1 ~> ℝ 1 ~> b ) -> ( ℝ 2 ~> b ) -uncurryD ( D b ) = D \ ( ℝ2 t0 s0 ) -> - let !(D1 ( D b_t0 ) ( T ( D dbdt_t0 ) ) ( T ( D d2bdt2_t0 ) ) ) = b ( ℝ1 t0 ) - !(D1 b_t0s0 dbds_t0s0 d2bds2_t0s0 ) = b_t0 ( ℝ1 s0 ) - !(D1 dbdt_t0s0 d2bdtds_t0s0 _ ) = dbdt_t0 ( ℝ1 s0 ) - !(D1 d2bdt2_t0s0 _ _ ) = d2bdt2_t0 ( ℝ1 s0 ) +uncurryD ( D b ) = D \ ( ℝ2 t0 s0 ) -> uncurryD' ( b ( ℝ1 t0 ) ) ( ℝ1 s0 ) + +uncurryD' :: D ( ℝ 1 ) ( ℝ 1 ~> b ) -> ℝ 1 -> D ( ℝ 2 ) b +uncurryD' ( D1 ( D b_t0 ) ( T ( D dbdt_t0 ) ) ( T ( D d2bdt2_t0 ) ) ) ( ℝ1 s0 ) = + let !( D1 b_t0s0 dbds_t0s0 d2bds2_t0s0 ) = b_t0 ( ℝ1 s0 ) + !( D1 dbdt_t0s0 d2bdtds_t0s0 _ ) = dbdt_t0 ( ℝ1 s0 ) + !( D1 d2bdt2_t0s0 _ _ ) = d2bdt2_t0 ( ℝ1 s0 ) in D2 b_t0s0 ( T dbdt_t0s0 ) dbds_t0s0 ( T d2bdt2_t0s0 ) d2bdtds_t0s0 d2bds2_t0s0 -------------------------------------------------------------------------------- +type Diffy :: Type -> Constraint +class ( Applicative ( D v ) + , Traversable ( D v ) + , Module Double ( T v ) ) + => Diffy v where + chain :: ( Module Double ( T w ) ) + => D ( ℝ 1 ) v -> D v w -> D ( ℝ 1 ) w + konst :: Module Double ( T w ) => w -> D v w + value :: D v w -> w + linear :: Module Double ( T w ) => ( v -> w ) -> ( v ~> w ) + chainRule :: ( Diffy v, Module Double ( T w ) ) - => ( ( ℝ 1 ) ~> v ) - -> ( v ~> w ) - -> ( ( ℝ 1 ) ~> w ) -chainRule ( D f ) ( D g ) = + => ( ( ℝ 1 ) ~> v ) -> ( v ~> w ) -> ( ( ℝ 1 ) ~> w ) +chainRule ( D df ) ( D dg ) = D \ x -> - case f x of - df@( D1 { v = f_x } ) -> - chain df ( g f_x ) + case df x of + df_x@( D1 { v = f_x } ) -> + chain df_x ( dg f_x ) -- | Recover the underlying function, discarding all infinitesimal information. fun :: forall v w. Diffy v => ( v ~> w ) -> ( v -> w ) -fun ( D f ) = value @v . f +fun ( D df ) = value @v . df +{-# INLINE fun #-} -type Diffy :: Type -> Constraint -class Diffy v where - chain :: ( Module Double ( T w ) ) - => Dℝ1 v -> D v w -> D ( ℝ 1 ) w - konst :: Module Double ( T w ) => w -> D v w - value :: D v w -> w +var :: ( Representable u, Diffy u ) => Fin ( Dim u ) -> ( u ~> Double ) +var i = linear ( `index` i ) +{-# INLINE var #-} instance Diffy ( ℝ 0 ) where chain _ ( D0 w ) = D1 w origin origin + {-# INLINE chain #-} konst k = D0 k + {-# INLINE konst #-} value ( D0 w ) = w + {-# INLINE value #-} + linear f = D \ _ -> D0 ( f ℝ0 ) + {-# INLINE linear #-} instance Diffy ( ℝ 1 ) where chain ( D1 _ ( T ( ℝ1 x' ) ) ( T ( ℝ1 x'' ) ) ) ( D1 v g_x g_xx ) = D1 v ( x' *^ g_x ) ( x'' *^ g_x ^+^ ( x' * x' ) *^ g_xx ) + {-# INLINE chain #-} konst k = D1 k origin origin + {-# INLINE konst #-} value ( D1 { v } ) = v + {-# INLINE value #-} + linear f = D \ u -> D1 ( f u ) ( T $ f u ) origin + {-# INLINE linear #-} instance Diffy ( ℝ 2 ) where chain ( D1 _ ( T ( ℝ2 x' y' ) ) ( T ( ℝ2 x'' y'' ) ) ) ( D2 v g_x g_y g_xx g_xy g_yy ) @@ -254,8 +283,16 @@ instance Diffy ( ℝ 2 ) where ( x'' *^ g_x ^+^ y'' *^ g_y ^+^ ( x' * x' ) *^ g_xx ^+^ ( y' * y' ) *^ g_yy ^+^ 2 *^ ( ( x' * y' ) *^ g_xy ) ) + {-# INLINE chain #-} konst k = D2 k origin origin origin origin origin + {-# INLINE konst #-} value ( D2 { v } ) = v + {-# INLINE value #-} + linear f = D \ u -> + D2 ( f u ) + ( T $ f ( ℝ2 1 0 ) ) ( T $ f ( ℝ2 0 1 ) ) + origin origin origin + {-# INLINE linear #-} instance Diffy ( ℝ 3 ) where chain ( D1 _ ( T ( ℝ3 x' y' z' ) ) ( T ( ℝ3 x'' y'' z'' ) ) ) @@ -265,24 +302,13 @@ instance Diffy ( ℝ 3 ) where ( x'' *^ g_x ^+^ y'' *^ g_y ^+^ z'' *^ g_z ^+^ ( x' * x' ) *^ g_xx ^+^ ( y' * y' ) *^ g_yy ^+^ ( z' * z' ) *^ g_zz ^+^ 2 *^ ( ( x' * y' ) *^ g_xy ) ^+^ ( x' * z' ) *^ g_xz ^+^ ( y' * z' ) *^ g_yz ) + {-# INLINE chain #-} konst k = D3 k origin origin origin origin origin origin origin origin origin + {-# INLINE konst #-} value ( D3 { v } ) = v - --------------------------------------------------------------------------------- - -type Var :: Nat -> Type -> Constraint -class Var n v where - var :: v ~> Double - -instance Var 1 ( ℝ 1 ) where - var = D \ ( ℝ1 x ) -> D1 x ( T 1 ) ( T 0 ) -instance Var 1 ( ℝ 2 ) where - var = D \ ( ℝ2 x _ ) -> D2 x ( T 1 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) -instance Var 2 ( ℝ 2 ) where - var = D \ ( ℝ2 _ y ) -> D2 y ( T 0 ) ( T 1 ) ( T 0 ) ( T 0 ) ( T 0 ) -instance Var 1 ( ℝ 3 ) where - var = D \ ( ℝ3 x _ _ ) -> D3 x ( T 1 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) -instance Var 2 ( ℝ 3 ) where - var = D \ ( ℝ3 _ y _ ) -> D3 y ( T 0 ) ( T 1 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) -instance Var 3 ( ℝ 3 ) where - var = D \ ( ℝ3 _ _ z ) -> D3 z ( T 0 ) ( T 0 ) ( T 1 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) ( T 0 ) + {-# INLINE value #-} + linear f = D \ u -> + D3 ( f u ) + ( T $ f ( ℝ3 1 0 0 ) ) ( T $ f ( ℝ3 0 1 0 ) ) ( T $ f ( ℝ3 0 0 1 ) ) + origin origin origin origin origin origin + {-# INLINE linear #-} diff --git a/src/splines/Math/Module.hs b/src/splines/Math/Module.hs index db952dd..b5dbc0d 100644 --- a/src/splines/Math/Module.hs +++ b/src/splines/Math/Module.hs @@ -5,6 +5,7 @@ module Math.Module ( Module(..), lerp , Inner(..) + , Interpolatable , norm, squaredNorm, quadrance, distance , proj, projC, closestPointOnSegment , cross @@ -17,6 +18,8 @@ import Control.Applicative ( liftA2 ) import Control.Monad ( guard ) +import Data.Kind + ( Type, Constraint ) import Data.Monoid ( Ap(..), Sum(..) ) @@ -109,6 +112,13 @@ closestPointOnSegment c ( Segment p0 p1 ) -------------------------------------------------------------------------------- +-- | A convenient constraint synonym for types that support interpolation. +type Interpolatable :: Type -> Constraint +class ( Torsor ( T r ) r, Module Double ( T r ) ) => Interpolatable r +instance ( Torsor ( T r ) r, Module Double ( T r ) ) => Interpolatable r + +-------------------------------------------------------------------------------- + instance Num a => Module a ( Sum a ) where origin = Sum 0 diff --git a/src/splines/Math/Roots.hs b/src/splines/Math/Roots.hs index c3fd8bb..2f679c7 100644 --- a/src/splines/Math/Roots.hs +++ b/src/splines/Math/Roots.hs @@ -3,7 +3,18 @@ {-# OPTIONS_GHC -Wno-name-shadowing #-} -module Math.Roots where +module Math.Roots + ( -- * Quadratic solver + solveQuadratic + + -- * Laguerre method + , roots, realRoots + , laguerre, eval, derivative + + -- * Newton–Raphson + , newtonRaphson + ) + where -- base import Control.Monad @@ -72,6 +83,9 @@ solveQuadratic a0 a1 a2 disc :: a disc = a1 * a1 - 4 * a0 * a2 +-------------------------------------------------------------------------------- +-- Root finding using Laguerre's method + -- | Find real roots of a polynomial with real coefficients. -- -- Coefficients are given in order of decreasing degree, e.g.: @@ -279,118 +293,140 @@ derivative p = do pure p' -------------------------------------------------------------------------------- - -data NewtonRaphson - = NewtonRaphson - { f :: Double -> (# Double, Double #) - , maxIters :: Word - , min_x, max_x :: Double - , digits :: Int } +-- Newton–Raphson -- Newton–Raphson implementation taken from Boost C++ library: -- https://github.com/boostorg/math/blob/0dc6a70caa6bbec2b6ae25eede36c430f0ccae13/include/boost/math/tools/roots.hpp#L217 - -newtonRaphson :: NewtonRaphson - -> Double - -> Maybe Double -newtonRaphson ( NewtonRaphson { f, maxIters, min_x, max_x, digits } ) x0 = - go $ - NewtonRaphsonState - { f_x_prev = 0 - , x = x0 - , δ = maxDouble, δ1 = maxDouble - , iters = 0 - , min_x, max_x - , min_f_x = 0, max_f_x = 0 } +{-# SPECIALISE + newtonRaphson + :: Word -> Int -> ( Double, Double ) -> ( Double -> ( Double, Double ) ) -> Double -> Maybe Double + #-} +{-# INLINEABLE newtonRaphson #-} +newtonRaphson :: ( RealFloat r, Show r ) + => Word -- ^ maximum number of iterations + -> Int -- ^ desired digits of precision + -> ( r, r ) -- ^ @(min_x, max_x)@. + -> ( r -> ( r, r ) ) -- ^ function and its derivative + -> r -- ^ initial guess + -> Maybe r +newtonRaphson maxIters digits ( min_x, max_x ) f x0 = + doNewtonRaphson f maxIters factor min_x max_x 0 0 0 x0 maxRealFloat maxRealFloat where - factor = encodeFloat 1 ( 1 - digits ) - go ( NewtonRaphsonState { f_x_prev, x, δ = δ1, δ1 = δ2, iters, min_x, max_x, min_f_x, max_f_x } ) - = case f x of - (# f_x, f'_x #) - | f_x == 0 + !factor = encodeFloat 1 ( 1 - digits ) + +doNewtonRaphson :: ( Fractional r, Ord r, Show r ) + => ( r -> (r, r) ) + -> Word + -> r + -> r -> r + -> r -> r + -> r + -> r + -> r -> r + -> Maybe r +doNewtonRaphson f maxIters factor min_x max_x min_f_x max_f_x f_x_prev x δ1 δ2 = + go min_x max_x min_f_x max_f_x f_x_prev 0 x δ1 δ2 + where + go min_x max_x min_f_x max_f_x f_x_prev !iters !x !δ1 !δ2 = + case f x of + ( f_x, f'_x ) + | f_x == 0 + -> Just x + | ( new_x, δ, δ1 ) <- newtonRaphsonStep f min_x max_x f_x_prev x f_x f'_x δ1 δ2 + -> if + | abs δ <= abs ( new_x * factor ) + || new_x == min_x || new_x == max_x -> Just x - | δ <- - if | f'_x == 0 - -> handleZeroDerivative f_x_prev x f_x δ1 min_x max_x - | otherwise - -> f_x / f'_x - , (# δ, δ1 #) <- - if | abs ( δ * δ ) > abs δ2 - , let shift = if δ > 0 then 0.5 * ( x - min_x ) else 0.5 * ( x - max_x ) - , δ <- if x /= 0 && ( abs shift > abs x ) - then signum δ * abs x * 1.1 - else shift - -> (# δ, 3 * δ #) - | otherwise - -> (# δ, δ1 #) - , let new_x = x - δ - , (# δ, new_x #) <- - if | new_x <= min_x - , δ <- 0.5 * ( x - min_x ) - -> (# δ, x - δ #) - | new_x >= max_x - , δ <- 0.5 * ( x - max_x ) - -> (# δ, x - δ #) - | otherwise - -> (# δ, new_x #) - -> if - | abs δ <= abs ( new_x * factor ) - || new_x == min_x || new_x == max_x - -> Just x - | iters >= maxIters - -> Nothing - | (# min_x, max_x, min_f_x, max_f_x #) <- - if δ > 0 - then (# min_x, x, min_f_x, f_x #) - else (# x, max_x, f_x, max_f_x #) - -> if min_f_x * max_f_x > 0 - then Nothing - else - go $ NewtonRaphsonState - { f_x_prev = f_x, x = new_x - , δ, δ1 - , iters = iters + 1 - , min_x, max_x - , min_f_x, max_f_x - } + | iters >= maxIters + -> Nothing + | ( min_x, max_x, min_f_x, max_f_x ) <- + if δ > 0 + then ( min_x, x, min_f_x, f_x ) + else ( x, max_x, f_x, max_f_x ) + -> if min_f_x * max_f_x > 0 + then Nothing + else + go min_x max_x min_f_x max_f_x f_x ( iters + 1 ) new_x δ δ1 - handleZeroDerivative :: Double - -> Double -> Double - -> Double - -> Double -> Double - -> Double - handleZeroDerivative f_x_prev x f_x δ min_x max_x - -- Handle zero derivative on first iteration. - | f_x_prev == 0 - , x_prev <- if x <= 0.5 * ( min_x + max_x ) then max_x else min_x - , (# f_x_prev, _ #) <- f x_prev - , δ <- x_prev - x - = finish f_x_prev δ - | otherwise - = finish f_x_prev δ +newtonRaphsonStep :: ( Fractional r, Ord r, Show r ) + => ( r -> ( r, r ) ) + -> r + -> r -> r -> r + -> r -> r + -> r -> r + -> ( r, r, r ) +newtonRaphsonStep f min_x max_x f_x_prev x f_x f'_x δ1 δ2 + | δ <- + if | f'_x == 0 + -> zeroDerivativeStep f min_x max_x f_x_prev x f_x δ1 + | otherwise + -> f_x / f'_x + , ( δ, δ1 ) <- + if | abs ( δ * δ ) > abs δ2 + -> bisectStep min_x max_x x δ + | otherwise + -> ( δ, δ1 ) + , let !new_x = x - δ + , ( δ, new_x ) <- + if | new_x <= min_x + , δ <- 0.5 * ( x - min_x ) + -> ( δ, x - δ ) + | new_x >= max_x + , δ <- 0.5 * ( x - max_x ) + -> ( δ, x - δ ) + | otherwise + -> ( δ, new_x ) + = ( new_x, δ, δ1 ) - where - finish f_x_prev δ - | signum f_x_prev * signum f_x < 0 - = if δ < 0 then 0.5 * ( x - min_x ) else 0.5 * ( x - max_x ) - | otherwise - = if δ < 0 then 0.5 * ( x - max_x ) else 0.5 * ( x - min_x ) - --- | Loop state for the 'newtonRaphson' function. -data NewtonRaphsonState = - NewtonRaphsonState - { f_x_prev :: !Double - , x :: !Double - , δ, δ1 :: !Double - , iters :: !Word - , min_x, max_x :: !Double - , min_f_x, max_f_x :: !Double } - -maxDouble :: Double -maxDouble = encodeFloat m n +-- Handle \( f'(x_0) = 0 \). +zeroDerivativeStep :: ( Fractional r, Ord r, Show r ) + => ( r -> ( r, r ) ) + -> r + -> r -> r + -> r + -> r -> r + -> r +zeroDerivativeStep f min_x max_x f_x_prev x f_x δ + -- Handle zero derivative on first iteration. + | f_x_prev == 0 + , x_prev <- if x <= 0.5 * ( min_x + max_x ) then max_x else min_x + , ( f_x_prev, _ ) <- f x_prev + , δ <- x_prev - x + = finish f_x_prev δ + | otherwise + = finish f_x_prev δ where - b = floatRadix ( 0 :: Double ) - e = floatDigits ( 0 :: Double ) - (_, e') = floatRange ( 0 :: Double ) + finish f_x_prev δ + | signum f_x_prev * signum f_x < 0 + = if δ < 0 then 0.5 * ( x - min_x ) else 0.5 * ( x - max_x ) + | otherwise + = if δ < 0 then 0.5 * ( x - max_x ) else 0.5 * ( x - min_x ) + +-- Handle failure of convergence in the past two steps. +bisectStep :: ( Fractional r, Ord r, Show r ) => r -> r -> r -> r -> ( r, r ) +bisectStep min_x max_x x δ + | let !shift = if δ > 0 then 0.5 * ( x - min_x ) else 0.5 * ( x - max_x ) + , δ <- if x /= 0 && ( abs shift > abs x ) + then signum δ * abs x * 1.1 + else shift + = ( δ, 3 * δ ) + +{-# SPECIALISE maxRealFloat :: Float #-} +{-# SPECIALISE maxRealFloat :: Double #-} +maxRealFloat :: forall r. RealFloat r => r +maxRealFloat = encodeFloat m n + where + b = floatRadix @r 0 + e = floatDigits @r 0 + (_, e') = floatRange @r 0 m = b ^ e - 1 n = e' - e + +{- + +for some reason GHC isn't able to optimise maxRealFloat @Double into + +maxDouble :: Double +maxDouble = D# 1.7976931348623157e308## + +-}