From a4e9c1cf32bcb0a2e0bf48ecb4a649c027856f4d Mon Sep 17 00:00:00 2001 From: sheaf Date: Fri, 13 Jan 2023 06:32:34 +0100 Subject: [PATCH] preparation for interval arithmetic --- MetaBrush.cabal | 2 + cabal.project | 4 +- src/app/MetaBrush/Action.hs | 2 +- src/app/MetaBrush/Application.hs | 4 +- src/app/MetaBrush/Asset/InfoBar.hs | 2 +- src/app/MetaBrush/Asset/Tools.hs | 2 +- src/app/MetaBrush/UI/InfoBar.hs | 2 +- src/metabrushes/MetaBrush/Asset/Brushes.hs | 4 +- src/metabrushes/MetaBrush/Brush.hs | 8 +- src/metabrushes/MetaBrush/Records.hs | 23 +++--- src/metabrushes/MetaBrush/Serialisable.hs | 3 +- src/splines/Math/Bezier/Stroke.hs | 91 +++++++++++++++------- src/splines/Math/Linear.hs | 33 +++++--- src/splines/Math/Linear/Dual.hs | 81 +++++++++---------- src/splines/Math/Module.hs | 74 ++++++++++++++++-- 15 files changed, 224 insertions(+), 111 deletions(-) diff --git a/MetaBrush.cabal b/MetaBrush.cabal index a156589..e5f9374 100644 --- a/MetaBrush.cabal +++ b/MetaBrush.cabal @@ -194,6 +194,8 @@ library splines ^>= 3.2.2.0 , prim-instances ^>= 0.2 + , rounded-hw + ^>= 0.3 , vector >= 0.12.1.2 && < 0.14 diff --git a/cabal.project b/cabal.project index 31b1714..73511f0 100644 --- a/cabal.project +++ b/cabal.project @@ -1,7 +1,9 @@ packages: . constraints: - acts -finitary + acts -finitary, + rounded-hw -pure-hs -c99 +avx512 +ghc-prim -x87-long-double + -- Fix a severe bug in Waargonaut (no corresponding Hackage release???) source-repository-package diff --git a/src/app/MetaBrush/Action.hs b/src/app/MetaBrush/Action.hs index 735e931..d951267 100644 --- a/src/app/MetaBrush/Action.hs +++ b/src/app/MetaBrush/Action.hs @@ -1078,7 +1078,7 @@ instance HandleAction Scroll where | dy > 0 = max 0.0078125 ( oldZoomFactor / sqrt 2 ) | otherwise - = min 256 ( oldZoomFactor * sqrt 2 ) + = min 4096 ( oldZoomFactor * sqrt 2 ) newCenter :: ℝ 2 newCenter = ( 1 - oldZoomFactor / newZoomFactor ) *^ ( oldCenter --> mousePos :: T ( ℝ 2 ) ) diff --git a/src/app/MetaBrush/Application.hs b/src/app/MetaBrush/Application.hs index 57db0eb..a59b501 100644 --- a/src/app/MetaBrush/Application.hs +++ b/src/app/MetaBrush/Application.hs @@ -200,8 +200,8 @@ runApplication application = do maxHistorySizeTVar <- STM.newTVarIO @Int 1000 fitParametersTVar <- STM.newTVarIO @FitParameters ( FitParameters - { maxSubdiv = 2 --3 -- 6 - , nbSegments = 3 --6 -- 12 + { maxSubdiv = 3 --2 --3 -- 6 + , nbSegments = 40 --3 --6 -- 12 , dist_tol = 0.1 -- 5e-3 , t_tol = 0.1 -- 1e-4 , maxIters = 2 -- 100 diff --git a/src/app/MetaBrush/Asset/InfoBar.hs b/src/app/MetaBrush/Asset/InfoBar.hs index ae271f5..8d00c1b 100644 --- a/src/app/MetaBrush/Asset/InfoBar.hs +++ b/src/app/MetaBrush/Asset/InfoBar.hs @@ -76,7 +76,7 @@ drawTopLeftCornerRect ( Colours { bg, viewport } ) = do Cairo.moveTo 12 24 Cairo.lineTo 12 14 Cairo.lineTo 24 14 - + Cairo.setLineWidth 4 withRGBA bg Cairo.setSourceRGBA Cairo.stroke diff --git a/src/app/MetaBrush/Asset/Tools.hs b/src/app/MetaBrush/Asset/Tools.hs index b6a9465..1369b95 100644 --- a/src/app/MetaBrush/Asset/Tools.hs +++ b/src/app/MetaBrush/Asset/Tools.hs @@ -129,7 +129,7 @@ drawPath ( Colours {..} ) = do Cairo.setLineWidth 2 - Cairo.newPath + Cairo.newPath Cairo.moveTo 7.179688 14.882813 Cairo.lineTo 31.820313 8.4375 Cairo.stroke diff --git a/src/app/MetaBrush/UI/InfoBar.hs b/src/app/MetaBrush/UI/InfoBar.hs index 84d54d9..23fc7ab 100644 --- a/src/app/MetaBrush/UI/InfoBar.hs +++ b/src/app/MetaBrush/UI/InfoBar.hs @@ -172,7 +172,7 @@ updateInfoBar viewportDrawingArea ( InfoBar {..} ) ( Variables { mousePosTVar } GTK.labelSetText zoomText $ Text.pack ( fixed 5 2 ( 100 * zoomFactor ) <> "%" ) case mbMousePos of Just ( ℝ2 mx my ) -> - GTK.labelSetText cursorPosText $ Text.pack ( "x: " <> fixed 6 2 mx <> "\ny: " <> fixed 6 2 my ) + GTK.labelSetText cursorPosText $ Text.pack ( "x: " <> fixed 6 4 mx <> "\ny: " <> fixed 6 4 my ) Nothing -> GTK.labelSetText cursorPosText $ Text.pack ( "x: " <> na <> "\ny: " <> na ) GTK.labelSetText topLeftPosText $ Text.pack ( "x: " <> fixed 6 2 l <> "\ny: " <> fixed 6 2 t ) diff --git a/src/metabrushes/MetaBrush/Asset/Brushes.hs b/src/metabrushes/MetaBrush/Asset/Brushes.hs index 130fe58..81ffcd3 100644 --- a/src/metabrushes/MetaBrush/Asset/Brushes.hs +++ b/src/metabrushes/MetaBrush/Asset/Brushes.hs @@ -93,7 +93,7 @@ circleBrush = e_x = pure $ ℝ2 1 0 e_y = pure $ ℝ2 0 1 - kon = konst @( Record CircleBrushFields ) + kon = konst @Double @( Record CircleBrushFields ) ellipseBrush :: Record EllipseBrushFields ~> Spline 'Closed () ( ℝ 2 ) ellipseBrush = @@ -112,4 +112,4 @@ ellipseBrush = e_x = pure $ ℝ2 1 0 e_y = pure $ ℝ2 0 1 - kon = konst @( Record EllipseBrushFields ) + kon = konst @Double @( Record EllipseBrushFields ) diff --git a/src/metabrushes/MetaBrush/Brush.hs b/src/metabrushes/MetaBrush/Brush.hs index e32fca4..16391e2 100644 --- a/src/metabrushes/MetaBrush/Brush.hs +++ b/src/metabrushes/MetaBrush/Brush.hs @@ -62,8 +62,8 @@ data Brush brushFields where BrushData :: forall brushFields . ( KnownSymbols brushFields - , Representable ( ℝ ( Length brushFields) ) - , Diffy ( ℝ ( Length brushFields) ) + , Representable Double ( ℝ ( Length brushFields) ) + , Diffy Double ( ℝ ( Length brushFields) ) , Typeable brushFields ) => { brushName :: !Text , brushFunction :: BrushFunction brushFields @@ -98,7 +98,7 @@ class ( KnownSymbols pointFields, Typeable pointFields , Show ( Record pointFields ) , NFData ( Record pointFields ) , Interpolatable ( Record pointFields ) - , Representable ( ℝ ( Length pointFields ) ) + , Representable Double ( ℝ ( Length pointFields ) ) ) => PointFields pointFields where { } instance ( KnownSymbols pointFields, Typeable pointFields @@ -106,7 +106,7 @@ instance ( KnownSymbols pointFields, Typeable pointFields , Show ( Record pointFields ) , NFData ( Record pointFields ) , Interpolatable ( Record pointFields ) - , Representable ( ℝ ( Length pointFields ) ) + , Representable Double ( ℝ ( Length pointFields ) ) ) => PointFields pointFields where { } diff --git a/src/metabrushes/MetaBrush/Records.hs b/src/metabrushes/MetaBrush/Records.hs index c80c117..1708c9e 100644 --- a/src/metabrushes/MetaBrush/Records.hs +++ b/src/metabrushes/MetaBrush/Records.hs @@ -80,7 +80,7 @@ deriving newtype => NFData ( Record ks ) -- | Show a record, using the given type-level field names. -instance ( KnownSymbols ks, Representable ( ℝ ( Length ks ) ) ) +instance ( KnownSymbols ks, Representable Double ( ℝ ( Length ks ) ) ) => Show ( Record ks ) where showsPrec p ( MkR r ) = showParen ( p >= 11 ) @@ -117,11 +117,11 @@ instance ( Torsor ( T ( ℝ ( Length ks ) ) ) ( ℝ ( Length ks ) ) MkR g --> MkR a = T $ MkR $ unT $ g --> a deriving newtype - instance Representable ( ℝ ( Length ks ) ) - => Representable ( Record ks ) + instance Representable r ( ℝ ( Length ks ) ) + => Representable r ( Record ks ) type instance D ( Record ks ) = D ( ℝ ( Length ks ) ) -deriving newtype instance Diffy ( ℝ ( Length ks ) ) => Diffy ( Record ks ) +deriving newtype instance Diffy Double ( ℝ ( Length ks ) ) => Diffy Double ( Record ks ) -------------------------------------------------------------------------------- @@ -148,8 +148,8 @@ intersect :: forall r1 r2 l1 l2 . ( Typeable r1, Typeable r2 , KnownSymbols r1, KnownSymbols r2 , l1 ~ Length r1, l2 ~ Length r2 - , Representable ( ℝ l1 ), Representable ( ℝ l2 ) - , Interpolatable ( Record r1 ), Diffy ( ℝ l2 ) + , Representable Double ( ℝ l1 ), Representable Double ( ℝ l2 ) + , Interpolatable ( Record r1 ), Diffy Double ( ℝ l2 ) ) => Intersection r1 r2 intersect @@ -172,8 +172,8 @@ data Intersection r1 r2 where :: forall r1r2 r1 r2 l12 . ( l12 ~ Length r1r2 , KnownSymbols r1r2 - , Representable ( ℝ l12 ) - , Diffy ( ℝ l12 ) + , Representable Double ( ℝ l12 ) + , Diffy Double ( ℝ l12 ) , Interpolatable ( Record r1r2 ) ) => { project :: Record r1 -> Record r1r2 , inject :: Record r2 -> Record r1r2 ~> Record r2 @@ -185,12 +185,13 @@ doIntersection :: forall r1 r2 l1 l2 kont . ( KnownSymbols r1, KnownSymbols r2 , l1 ~ Length r1, l2 ~ Length r2 - , Representable ( ℝ l1 ), Representable ( ℝ l2 ) + , Representable Double ( ℝ l1 ), Representable Double ( ℝ l2 ) ) => ( forall r1r2 l12. ( r1r2 ~ Intersect r1 r2, l12 ~ Length r1r2 - , Representable ( ℝ l12 ), Diffy ( ℝ l12 ), Interpolatable ( ℝ l12 ) - , KnownSymbols r1r2, Representable ( ℝ ( Length r1r2 ) ) + , Representable Double ( ℝ l12 ), Diffy Double ( ℝ l12 ) + , Interpolatable ( ℝ l12 ) + , KnownSymbols r1r2, Representable Double ( ℝ ( Length r1r2 ) ) ) => Proxy# r1r2 -> Vec l12 ( Fin l1 ) -> Vec l12 ( Fin l2 ) -> kont ) -> kont diff --git a/src/metabrushes/MetaBrush/Serialisable.hs b/src/metabrushes/MetaBrush/Serialisable.hs index d240600..dd953a2 100644 --- a/src/metabrushes/MetaBrush/Serialisable.hs +++ b/src/metabrushes/MetaBrush/Serialisable.hs @@ -110,7 +110,8 @@ instance Serialisable ( T ( ℝ 2 ) ) where JSON.Encoder.atKey' "x" encoder x . JSON.Encoder.atKey' "y" encoder y decoder = V2 <$> JSON.Decoder.atKey "x" decoder <*> JSON.Decoder.atKey "y" decoder -instance ( KnownSymbols ks, Representable ( ℝ ( Length ks ) ) ) => Serialisable ( Record ks ) where +instance ( KnownSymbols ks, Representable Double ( ℝ ( Length ks ) ) ) + => Serialisable ( Record ks ) where encoder = contramap encodeFields ( JSON.Encoder.keyValueTupleFoldable ( encoder @Double ) ) where encodeFields :: Record ks -> [ ( Text, Double ) ] diff --git a/src/splines/Math/Bezier/Stroke.hs b/src/splines/Math/Bezier/Stroke.hs index bcc367e..411264b 100644 --- a/src/splines/Math/Bezier/Stroke.hs +++ b/src/splines/Math/Bezier/Stroke.hs @@ -2,6 +2,8 @@ {-# LANGUAGE QuantifiedConstraints #-} {-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE DuplicateRecordFields #-} + module Math.Bezier.Stroke ( Offset(..) , CachedStroke(..), discardCache, invalidateCache @@ -106,9 +108,8 @@ import qualified Math.Bezier.Quadratic as Quadratic import Math.Epsilon ( epsilon ) import Math.Module - ( Module(..), Inner((^.^)), Interpolatable - , lerp, cross - , convexCombination, strictlyParallel + ( Module(..), Inner((^.^)), Cross(cross), Interpolatable + , lerp, convexCombination, strictlyParallel ) import Math.Orientation ( Orientation(..), splineOrientation @@ -181,7 +182,7 @@ computeStrokeOutline :: forall ( clo :: SplineType ) usedParams brushParams crvData ptData s . ( KnownSplineType clo , Interpolatable usedParams - , Diffy usedParams, Diffy brushParams + , Diffy Double usedParams, Diffy Double brushParams , HasType ( ℝ 2 ) ptData , HasType ( CachedStroke s ) crvData , NFData ptData, NFData crvData @@ -414,7 +415,7 @@ computeStrokeOutline fitParams ptParams toBrushParams brushFn spline@( Spline { outlineFunction :: forall usedParams brushParams crvData ptData . ( Interpolatable usedParams - , Diffy usedParams, Diffy brushParams + , Diffy Double usedParams, Diffy Double brushParams , HasType ( ℝ 2 ) ptData -- Debugging. , Show ptData, Show brushParams @@ -456,8 +457,9 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv = D1 path_t path'_t _ = runD path t D1 params_t _ _ = runD usedParams t - brush_t = value @brushParams - $ runD brushFromParams ( fun toBrushParams params_t ) + brush_t = value @Double @brushParams + $ runD brushFromParams + $ fun toBrushParams params_t in fwdBwd @@ -750,9 +752,10 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) -- -- - \( p(t) \) is the path that the brush follows, and -- - \( b(t,s) \) is the brush shape, as it varies along the path. -brushStroke :: D ( ℝ 1 ) ( ℝ 2 ) -- ^ stroke path \( p(t) \) - -> D ( ℝ 2 ) ( ℝ 2 ) -- ^ brush \( b(t,s) \) - -> D ( ℝ 2 ) ( ℝ 2 ) +brushStroke :: Module r ( T v ) + => D ( ℝ 1 ) v -- ^ stroke path \( p(t) \) + -> D ( ℝ 2 ) v -- ^ brush \( b(t,s) \) + -> D ( ℝ 2 ) v brushStroke ( D1 p dpdt d2pdt2 ) ( D2 b dbdt dbds d2bdt2 d2bdtds d2bds2 ) = D2 ( unT $ T p ^+^ T b ) -- c = p + b @@ -770,27 +773,33 @@ brushStroke ( D1 p dpdt d2pdt2 ) ( D2 b dbdt dbds d2bdt2 d2bdtds d2bds2 ) = -- -- \[ E = \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} = 0, ] -- --- as well as the vector +-- together with the total derivative -- --- \[ \frac{\partial E}{\partial s} \frac{\mathrm{d} c}{\mathrm{d} t} \] +-- \[ \frac{\mathrm{d} c}{\mathrm{d} t}, \] -- --- whose roots correspond to cusps in the envelope, and +-- and the partial derivatives -- --- \[ \frac{\partial E}{\partial s}. \] -envelopeEquation :: D ( ℝ 2 ) ( ℝ 2 ) -> ( Double, T ( ℝ 2 ), Double, Double ) +-- \[ \frac{\partial E}{\partial s}, \qquad \frac{\partial E}{\partial s}. \] +-- +-- NB: if \( \frac{\partial E}{\partial s} \) is zero, the total derivative is ill-defined. +envelopeEquation :: ( D ( i ( ℝ 2 ) ) ~ D ( ℝ 2 ) + , Cross ( i Double ) ( T ( i ( ℝ 2 ) ) ) + , Fractional ( i Double ) + ) + => D ( i ( ℝ 2 ) ) ( i ( ℝ 2 ) ) + -> ( i Double, T ( i ( ℝ 2 ) ), i Double, i 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 + tot = dcdt ^-^ ( dEdt / dEds ) *^ dcds in ( ee, tot, dEdt, dEds ) -- Computation of total derivative dc/dt: -- -- dc/dt = ∂c/∂t + ∂c/∂s ∂s/∂t - -- ∂s/∂t = - ( ∂E / ∂t ) / ( ∂E / ∂s ) + -- ∂s/∂t = - ∂E/∂t / ∂E/∂s -- - -- ( ∂E / ∂s ) dc/dt = ( ∂E / ∂s ) ∂c/∂t - ( ∂E / ∂t ) ∂c/∂s. - + -- ∂E/∂s dc/dt = ∂E/∂s ∂c/∂t - ∂E/∂t ∂c/∂s. -- | Linear interpolation, as a differentiable function. line :: forall b. ( Module Double ( T b ), Torsor ( T b ) b ) @@ -906,9 +915,27 @@ solveEnvelopeEquations path_t path'_t ( fwdOffset, bwdOffset ) strokeData 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 ) - + StrokeDatum { ee = _ee, dstroke, 𝛿E𝛿t = _𝛿E𝛿t, 𝛿E𝛿s, dcdt } -> + -- The total derivative dc/dt is computed by dividing by ∂E/∂s, + -- so check it isn't zero first. This corresponds to cusps in the envelope. + let totDeriv + | abs 𝛿E𝛿s < epsilon + , let s' = if s >= 0.5 then s - 1e-9 else s + 1e-9 + = case f ( ℝ1 s' ) of { StrokeDatum { dcdt = dcdt_s' } -> dcdt_s' } + | otherwise + = dcdt + in --trace + -- ( unlines + -- [ "solveEnvelopeEquations" + -- , " t = " ++ show t + -- , " s = " ++ show s + -- , " c = " ++ show dstroke + -- , " E = " ++ show _ee + -- , " ∂E/∂t = " ++ show _𝛿E𝛿t + -- , " ∂E/∂s = " ++ show 𝛿E𝛿s + -- , " dc/dt = " ++ show totDeriv + -- ] ) + ( good, ℝ1 s, value @Double @( ℝ 2 ) dstroke, totDeriv ) eqn :: ( ℝ 1 -> StrokeDatum ) -> ( Double -> ( Double, Double ) ) eqn f s = @@ -930,7 +957,7 @@ instance Applicative ZipSeq where liftA2 f ( ZipSeq xs ) ( ZipSeq ys ) = ZipSeq ( Seq.zipWith f xs ys ) brushStrokeData :: forall brushParams - . ( Diffy brushParams, Show brushParams ) + . ( Diffy Double brushParams, Show brushParams ) => ( ℝ 1 ~> ℝ 2 ) -- ^ path -> ( ℝ 1 ~> brushParams ) -- ^ brush parameters -> ( brushParams ~> SplinePts Closed ) -- ^ brush from parameters @@ -947,7 +974,7 @@ brushStrokeData path params brush = splines :: Seq ( D brushParams ( ℝ 1 ~> ℝ 2 ) ) !splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns ) dbrush_params dbrushes_t :: Seq ( ℝ 1 -> D ( ℝ 2 ) ( ℝ 2 ) ) - !dbrushes_t = force $ fmap ( uncurryD' . ( dparams_t `chain` ) ) splines + !dbrushes_t = force $ fmap ( uncurryD . ( dparams_t `chain` ) ) splines in fmap ( mkStrokeDatum dpath_t ) dbrushes_t where @@ -958,10 +985,11 @@ brushStrokeData path params brush = 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 + ( ee, dcdt, 𝛿E𝛿t, 𝛿E𝛿s ) = coerce $ envelopeEquation @Identity $ coerce dstroke in -- trace -- ( unlines -- [ "envelopeEquation:" + -- , " t = " ++ show t -- , " s = " ++ show s -- , " c = " ++ show _c -- , " ∂c/∂t = " ++ show _𝛿c𝛿t @@ -974,7 +1002,7 @@ brushStrokeData path params brush = { dpath = dpath_t , dbrush = dbrush_t_s , dstroke - , ee, dcdt, 𝛿E𝛿s } + , ee, dcdt, 𝛿E𝛿t, 𝛿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 @@ -994,11 +1022,16 @@ data StrokeDatum -- -- \[ E(t_0,s_0) = \left ( \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} \right )_{(t_0,s_0)}. \] , ee :: Double + -- | \( \left ( \frac{\partial E}{\partial s} \right )_{(t_0,s_0)}. \) + , 𝛿E𝛿s :: Double + -- | \( \left ( \frac{\partial E}{\partial t} \right )_{(t_0,s_0)}. \) + , 𝛿E𝛿t :: 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 + -- + -- This is ill-defined when \( \frac{\partial E}{\partial s} = 0 \). + , dcdt :: T ( ℝ 2 ) + } deriving stock Show diff --git a/src/splines/Math/Linear.hs b/src/splines/Math/Linear.hs index 5f1c62b..77848c0 100644 --- a/src/splines/Math/Linear.hs +++ b/src/splines/Math/Linear.hs @@ -12,6 +12,9 @@ module Math.Linear , Fin(..), eqFin, MFin(..) , Dim, Representable(..), injection, projection , Vec(..), (!), find + + -- * Intervals + , 𝕀, 𝕀ℝ ) where -- base @@ -44,6 +47,10 @@ import Data.Group import Data.Group.Generics ( ) +-- rounded-hw +import Numeric.Rounded.Hardware.Interval.NonEmpty + ( Interval ) + -------------------------------------------------------------------------------- data Mat22 = Mat22 !Double !Double !Double !Double @@ -171,24 +178,24 @@ type family Dim v type instance Dim ( ℝ n ) = n -type Representable :: Type -> Constraint -class Representable v where - tabulate :: ( Fin ( Dim v ) -> Double ) -> v - index :: v -> Fin ( Dim v ) -> Double +type Representable :: Type -> Type -> Constraint +class Representable r v | v -> r where + tabulate :: ( Fin ( Dim v ) -> r ) -> v + index :: v -> Fin ( Dim v ) -> r -instance Representable ( ℝ 0 ) where +instance Representable Double ( ℝ 0 ) where {-# INLINE tabulate #-} tabulate _ = ℝ0 {-# INLINE index #-} index _ _ = 0 -instance Representable ( ℝ 1 ) where +instance Representable Double ( ℝ 1 ) where {-# INLINE tabulate #-} tabulate f = ℝ1 ( f ( Fin 1## ) ) {-# INLINE index #-} index ( ℝ1 x ) _ = x -instance Representable ( ℝ 2 ) where +instance Representable Double ( ℝ 2 ) where {-# INLINE tabulate #-} tabulate f = ℝ2 ( f ( Fin 1## ) ) ( f ( Fin 2## ) ) {-# INLINE index #-} @@ -196,7 +203,7 @@ instance Representable ( ℝ 2 ) where Fin 1## -> x _ -> y -instance Representable ( ℝ 3 ) where +instance Representable Double ( ℝ 3 ) where {-# INLINE tabulate #-} tabulate f = ℝ3 ( f ( Fin 1## ) ) ( f ( Fin 2## ) ) ( f ( Fin 3## ) ) {-# INLINE index #-} @@ -206,14 +213,14 @@ instance Representable ( ℝ 3 ) where _ -> z {-# INLINE projection #-} -projection :: ( Representable u, Representable v ) +projection :: ( Representable r u, Representable r v ) => ( Fin ( Dim v ) -> Fin ( Dim u ) ) -> u -> v projection f = \ u -> tabulate \ i -> index u ( f i ) {-# INLINE injection #-} -injection :: ( Representable u, Representable v ) +injection :: ( Representable r u, Representable r v ) => ( Fin ( Dim v ) -> MFin ( Dim u ) ) -> u -> v -> v injection f = \ u v -> @@ -242,3 +249,9 @@ find eq v b = MFin ( go 1## v ) | otherwise = go ( j `plusWord#` 1## ) as go _ VZ = 0## + +-------------------------------------------------------------------------------- +-- Intervals. + +type 𝕀 = Interval +type 𝕀ℝ n = 𝕀 ( ℝ n ) diff --git a/src/splines/Math/Linear/Dual.hs b/src/splines/Math/Linear/Dual.hs index 94d7c47..ee66200 100644 --- a/src/splines/Math/Linear/Dual.hs +++ b/src/splines/Math/Linear/Dual.hs @@ -11,6 +11,8 @@ import Control.Applicative ( liftA2 ) import Data.Coerce ( coerce ) +import Data.Functor.Identity + ( Identity(..) ) import Data.Kind ( Type, Constraint ) import GHC.Generics @@ -28,11 +30,12 @@ type (~>) :: Type -> Type -> Type 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 +instance ( Applicative ( D u ), Module r ( T v ) ) + => Module r ( T ( u ~> v ) ) where + origin = T $ D \ _ -> pure $ coerce $ origin @r @( T v ) + T ( D f ) ^+^ T ( D g ) = T $ D \ t -> liftA2 ( coerce $ (^+^) @r @( T v ) ) ( f t ) ( g t ) + T ( D f ) ^-^ T ( D g ) = T $ D \ t -> liftA2 ( coerce $ (^-^) @r @( T v ) ) ( f t ) ( g t ) + a *^ T ( D f ) = T $ D \ t -> fmap ( coerce $ (*^) @r @( T v ) 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 @@ -42,6 +45,9 @@ type instance D ( ℝ 1 ) = Dℝ1 type instance D ( ℝ 2 ) = Dℝ2 type instance D ( ℝ 3 ) = Dℝ3 +type instance D ( Identity a ) = D a +type instance D ( 𝕀 u ) = D u + newtype Dℝ0 v = D0 { v :: v } deriving stock ( Show, Eq, Functor, Foldable, Traversable, Generic, Generic1 ) deriving newtype ( Num, Fractional, Floating ) @@ -67,7 +73,7 @@ instance Num ( Dℝ1 Double ) where (+) = liftA2 (+) (-) = liftA2 (-) negate = fmap negate - fromInteger = konst @( ℝ 1 ) . fromInteger + fromInteger = konst @Double @( ℝ 1 ) . fromInteger abs = error "no" signum = error "no" @@ -82,7 +88,7 @@ instance Num ( Dℝ2 Double ) where (+) = liftA2 (+) (-) = liftA2 (-) negate = fmap negate - fromInteger = konst @( ℝ 2 ) . fromInteger + fromInteger = konst @Double @( ℝ 2 ) . fromInteger abs = error "no" signum = error "no" @@ -101,7 +107,7 @@ instance Num ( Dℝ3 Double ) where (+) = liftA2 (+) (-) = liftA2 (-) negate = fmap negate - fromInteger = konst @( ℝ 3 ) . fromInteger + fromInteger = konst @Double @( ℝ 3 ) . fromInteger abs = error "no" signum = error "no" @@ -146,9 +152,9 @@ instance Module Double ( T v ) => Module ( Dℝ3 Double ) ( Dℝ3 v ) where instance Fractional ( Dℝ1 Double ) where (/) = error "I haven't yet defined (/) for Dℝ1" - fromRational = konst @( ℝ 1 ) . fromRational + fromRational = konst @Double @( ℝ 1 ) . fromRational instance Floating ( Dℝ1 Double ) where - pi = konst @( ℝ 1 ) pi + pi = konst @Double @( ℝ 1 ) pi sin ( D1 v ( T dx ) ( T ddx ) ) = let !s = sin v !c = cos v @@ -161,9 +167,9 @@ instance Floating ( Dℝ1 Double ) where instance Fractional ( Dℝ2 Double ) where (/) = error "I haven't yet defined (/) for Dℝ2" - fromRational = konst @( ℝ 2 ) . fromRational + fromRational = konst @Double @( ℝ 2 ) . fromRational instance Floating ( Dℝ2 Double ) where - pi = konst @( ℝ 2 ) pi + pi = konst @Double @( ℝ 2 ) pi sin ( D2 v ( T dx ) ( T dy ) ( T ddx ) ( T dxdy ) ( T ddy ) ) = let !s = sin v !c = cos v @@ -184,9 +190,9 @@ instance Floating ( Dℝ2 Double ) where instance Fractional ( Dℝ3 Double ) where (/) = error "I haven't yet defined (/) for Dℝ3" - fromRational = konst @( ℝ 3 ) . fromRational + fromRational = konst @Double @( ℝ 3 ) . fromRational instance Floating ( Dℝ3 Double ) where - pi = konst @( ℝ 3 ) pi + pi = konst @Double @( ℝ 3 ) pi sin ( D3 v ( T dx ) ( T dy ) ( T dz ) ( T ddx ) ( T dxdy ) ( T ddy ) ( T dxdz ) ( T dydz ) ( T ddz ) ) = let !s = sin v !c = cos v @@ -213,30 +219,24 @@ instance Floating ( Dℝ3 Double ) where -------------------------------------------------------------------------------- -uncurryD :: ( ℝ 1 ~> ℝ 1 ~> b ) -> ( ℝ 2 ~> b ) -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 ) +uncurryD :: D a ~ D ( ℝ 1 ) + => D ( ℝ 1 ) ( a ~> b ) -> a -> D ( ℝ 2 ) b +uncurryD ( D1 ( D b_t0 ) ( T ( D dbdt_t0 ) ) ( T ( D d2bdt2_t0 ) ) ) s0 = + let !( D1 b_t0s0 dbds_t0s0 d2bds2_t0s0 ) = b_t0 s0 + !( D1 dbdt_t0s0 d2bdtds_t0s0 _ ) = dbdt_t0 s0 + !( D1 d2bdt2_t0s0 _ _ ) = d2bdt2_t0 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 ) +type Diffy :: Type -> Type -> Constraint +class ( Traversable ( D v ), Module r ( T v ) ) => Diffy r v where + chain :: ( Module r ( T w ) ) => D ( ℝ 1 ) v -> D v w -> D ( ℝ 1 ) w + konst :: Module r ( T w ) => w -> D v w + value :: D v w -> w + linear :: Module r ( T w ) => ( v -> w ) -> ( v ~> w ) -chainRule :: ( Diffy v, Module Double ( T w ) ) +chainRule :: ( Diffy r v, Module r ( T w ) ) => ( ( ℝ 1 ) ~> v ) -> ( v ~> w ) -> ( ( ℝ 1 ) ~> w ) chainRule ( D df ) ( D dg ) = D \ x -> @@ -245,15 +245,16 @@ chainRule ( D df ) ( D dg ) = 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 df ) = value @v . df +fun :: forall v w r. Diffy r v => ( v ~> w ) -> ( v -> w ) +fun ( D df ) = value @r @v . df {-# INLINE fun #-} -var :: ( Representable u, Diffy u ) => Fin ( Dim u ) -> ( u ~> Double ) +var :: forall u r. ( Module r ( T r ), Representable r u, Diffy r u ) + => Fin ( Dim u ) -> ( u ~> r ) var i = linear ( `index` i ) {-# INLINE var #-} -instance Diffy ( ℝ 0 ) where +instance Diffy Double ( ℝ 0 ) where chain _ ( D0 w ) = D1 w origin origin {-# INLINE chain #-} konst k = D0 k @@ -263,7 +264,7 @@ instance Diffy ( ℝ 0 ) where linear f = D \ _ -> D0 ( f ℝ0 ) {-# INLINE linear #-} -instance Diffy ( ℝ 1 ) where +instance Diffy Double ( ℝ 1 ) where chain ( D1 _ ( T ( ℝ1 x' ) ) ( T ( ℝ1 x'' ) ) ) ( D1 v g_x g_xx ) = D1 v ( x' *^ g_x ) @@ -276,7 +277,7 @@ instance Diffy ( ℝ 1 ) where linear f = D \ u -> D1 ( f u ) ( T $ f u ) origin {-# INLINE linear #-} -instance Diffy ( ℝ 2 ) where +instance Diffy Double ( ℝ 2 ) where chain ( D1 _ ( T ( ℝ2 x' y' ) ) ( T ( ℝ2 x'' y'' ) ) ) ( D2 v g_x g_y g_xx g_xy g_yy ) = D1 v ( x' *^ g_x ^+^ y' *^ g_y ) @@ -294,7 +295,7 @@ instance Diffy ( ℝ 2 ) where origin origin origin {-# INLINE linear #-} -instance Diffy ( ℝ 3 ) where +instance Diffy Double ( ℝ 3 ) where chain ( D1 _ ( T ( ℝ3 x' y' z' ) ) ( T ( ℝ3 x'' y'' z'' ) ) ) ( D3 v g_x g_y g_z g_xx g_xy g_yy g_xz g_yz g_zz ) = D1 v diff --git a/src/splines/Math/Module.hs b/src/splines/Math/Module.hs index b5dbc0d..92d44f0 100644 --- a/src/splines/Math/Module.hs +++ b/src/splines/Math/Module.hs @@ -4,12 +4,12 @@ module Math.Module ( Module(..), lerp - , Inner(..) + , Inner(..), Cross(..) , Interpolatable , norm, squaredNorm, quadrance, distance , proj, projC, closestPointOnSegment - , cross , strictlyParallel, convexCombination + , 𝕀 ) where @@ -18,6 +18,10 @@ import Control.Applicative ( liftA2 ) import Control.Monad ( guard ) +import Data.Coerce + ( coerce ) +import Data.Functor.Identity + ( Identity(..) ) import Data.Kind ( Type, Constraint ) import Data.Monoid @@ -31,6 +35,14 @@ import Data.Act ( (-->) ) ) +-- rounded-hw +import Numeric.Rounded.Hardware + ( Rounded(..) ) +import Numeric.Rounded.Hardware.Interval.NonEmpty + ( Interval(..) ) +import Numeric.Rounded.Hardware.Class + ( intervalAdd, intervalSub, intervalMul ) + -- MetaBrush import Math.Epsilon ( epsilon ) @@ -69,6 +81,9 @@ infixl 8 ^.^ class Module r m => Inner r m where (^.^) :: m -> m -> r +class Module r m => Cross r m where + cross :: m -> m -> r + -- | Norm of a vector, computed using the inner product. norm :: forall m r. ( Floating r, Inner r m ) => m -> r norm = sqrt . squaredNorm @@ -114,8 +129,8 @@ 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 +class ( Torsor ( T u ) u, Module Double ( T u ) ) => Interpolatable u +instance ( Torsor ( T u ) u, Module Double ( T u ) ) => Interpolatable u -------------------------------------------------------------------------------- @@ -158,9 +173,18 @@ instance Module Double ( T ( ℝ 3 ) ) where instance Inner Double ( T ( ℝ 2 ) ) where V2 x1 y1 ^.^ V2 x2 y2 = x1 * x2 + y1 * y2 --- | Cross-product of two 2D vectors. -cross :: T ( ℝ 2 ) -> T ( ℝ 2 ) -> Double -cross ( V2 x1 y1 ) ( V2 x2 y2 ) = x1 * y2 - x2 * y1 +instance Cross Double ( T ( ℝ 2 ) ) where + cross ( V2 x1 y1 ) ( V2 x2 y2 ) = x1 * y2 - x2 * y1 + +instance Module r ( T m ) => Module ( Identity r ) ( T ( Identity m ) ) where + origin = coerce $ origin @r @( T m ) + (^+^) = coerce $ (^+^) @r @( T m ) + (^-^) = coerce $ (^-^) @r @( T m ) + (*^) = coerce $ (*^) @r @( T m ) +instance Inner r ( T m ) => Inner ( Identity r ) ( T ( Identity m ) ) where + (^.^) = coerce $ (^.^) @r @( T m ) +instance Cross r ( T m ) => Cross ( Identity r ) ( T ( Identity m ) ) where + cross = coerce $ cross @r @( T m ) -- | Compute whether two vectors point in the same direction, -- that is, whether each vector is a (strictly) positive multiple of the other. @@ -199,3 +223,39 @@ convexCombination v0 v1 u c0, c10 :: Double c0 = v0 `cross` u c10 = ( v0 ^-^ v1 ) `cross` u + +-------------------------------------------------------------------------------- +-- Interval arithmetic using rounded-hw library. + +instance Module ( 𝕀 Double ) ( T ( 𝕀ℝ 2 ) ) where + origin = T ( I ( Rounded ( ℝ2 0 0 ) ) ( Rounded ( ℝ2 0 0 ) ) ) + T ( I ( Rounded ( ℝ2 x1_lo y1_lo ) ) ( Rounded ( ℝ2 x1_hi y1_hi ) ) ) ^+^ + T ( I ( Rounded ( ℝ2 x2_lo y2_lo ) ) ( Rounded ( ℝ2 x2_hi y2_hi ) ) ) + = let !( Rounded x_lo, Rounded x_hi ) = intervalAdd ( Rounded x1_lo ) ( Rounded x1_hi ) ( Rounded x2_lo ) ( Rounded x2_hi ) + !( Rounded y_lo, Rounded y_hi ) = intervalAdd ( Rounded y1_lo ) ( Rounded y1_hi ) ( Rounded y2_lo ) ( Rounded y2_hi ) + in T ( I ( Rounded ( ℝ2 x_lo y_lo ) ) ( Rounded ( ℝ2 x_hi y_hi ) ) ) + T ( I ( Rounded ( ℝ2 x1_lo y1_lo ) ) ( Rounded ( ℝ2 x1_hi y1_hi ) ) ) ^-^ + T ( I ( Rounded ( ℝ2 x2_lo y2_lo ) ) ( Rounded ( ℝ2 x2_hi y2_hi ) ) ) + = let !( Rounded x_lo, Rounded x_hi ) = intervalSub ( Rounded x1_lo ) ( Rounded x1_hi ) ( Rounded x2_lo ) ( Rounded x2_hi ) + !( Rounded y_lo, Rounded y_hi ) = intervalSub ( Rounded y1_lo ) ( Rounded y1_hi ) ( Rounded y2_lo ) ( Rounded y2_hi ) + in T ( I ( Rounded ( ℝ2 x_lo y_lo ) ) ( Rounded ( ℝ2 x_hi y_hi ) ) ) + I ( Rounded k_lo ) ( Rounded k_hi ) *^ T ( I ( Rounded ( ℝ2 x1_lo y1_lo ) ) ( Rounded ( ℝ2 x1_hi y1_hi ) ) ) + = let !( Rounded x_lo, Rounded x_hi ) = intervalMul ( Rounded k_lo ) ( Rounded k_hi ) ( Rounded x1_lo ) ( Rounded x1_hi ) + !( Rounded y_lo, Rounded y_hi ) = intervalMul ( Rounded k_lo ) ( Rounded k_hi ) ( Rounded y1_lo ) ( Rounded y1_hi ) + in T ( I ( Rounded ( ℝ2 x_lo y_lo ) ) ( Rounded ( ℝ2 x_hi y_hi ) ) ) + +instance Inner ( 𝕀 Double ) ( T ( 𝕀ℝ 2 ) ) where + T ( I ( Rounded ( ℝ2 x1_lo y1_lo ) ) ( Rounded ( ℝ2 x1_hi y1_hi ) ) ) ^.^ + T ( I ( Rounded ( ℝ2 x2_lo y2_lo ) ) ( Rounded ( ℝ2 x2_hi y2_hi ) ) ) + = let !( x_lo, x_hi ) = intervalMul ( Rounded x1_lo ) ( Rounded x1_hi ) ( Rounded x2_lo ) ( Rounded x2_hi ) + !( y_lo, y_hi ) = intervalMul ( Rounded y1_lo ) ( Rounded y1_hi ) ( Rounded y2_lo ) ( Rounded y2_hi ) + !( z_lo, z_hi ) = intervalAdd x_lo x_hi y_lo y_hi + in I z_lo z_hi + +instance Cross ( 𝕀 Double ) ( T ( 𝕀ℝ 2 ) ) where + T ( I ( Rounded ( ℝ2 x1_lo y1_lo ) ) ( Rounded ( ℝ2 x1_hi y1_hi ) ) ) `cross` + T ( I ( Rounded ( ℝ2 x2_lo y2_lo ) ) ( Rounded ( ℝ2 x2_hi y2_hi ) ) ) + = let !( x_lo, x_hi ) = intervalMul ( Rounded x1_lo ) ( Rounded x1_hi ) ( Rounded y2_lo ) ( Rounded y2_hi ) + !( y_lo, y_hi ) = intervalMul ( Rounded x2_lo ) ( Rounded x2_hi ) ( Rounded y1_lo ) ( Rounded y1_hi ) + !( z_lo, z_hi ) = intervalSub x_lo x_hi y_lo y_hi + in I z_lo z_hi