WIP on monomial bases

MetaBrush.cabal

@@ -95,7 +95,6 @@ common common
         -fexcess-precision
         -fspecialise-aggressively
         -optc-O3
-        -optc-ffast-math
         -Wall
         -Wcompat
         -fwarn-missing-local-signatures

src/splines/Math/Bezier/Stroke.hs

@@ -26,7 +26,7 @@ import Control.Arrow
 import Control.Applicative
   ( Applicative(..) )
 import Control.Monad
-  ( unless )
+  ( guard, unless )
 import Control.Monad.ST
   ( RealWorld, ST )
 import Data.Bifunctor
@@ -34,7 +34,7 @@ import Data.Bifunctor
 import Data.Coerce
   ( Coercible, coerce )
 import Data.Foldable
-  ( for_ )
+  ( for_, toList )
 import Data.Functor.Identity
   ( Identity(..) )
 import Data.List.NonEmpty
@@ -82,6 +82,13 @@ import Data.Generics.Internal.VL
 import qualified Control.Parallel.Strategies as Strats
   ( rdeepseq, parTuple2, using )
 
+-- rounded-hw
+import Numeric.Rounded.Hardware
+  ( Rounded(..) )
+import Numeric.Rounded.Hardware.Interval.NonEmpty
+  ( Interval(I) )
+import qualified Numeric.Rounded.Hardware.Interval.NonEmpty as Interval
+
 -- transformers
 import Control.Monad.Trans.Class
   ( lift )
@@ -122,7 +129,7 @@ import Math.Roots
 import Math.Linear
 import Math.Linear.Dual
 
---import Debug.Utils
+import Debug.Utils
 
 --------------------------------------------------------------------------------
 
@@ -482,7 +489,6 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv =
     path :: ℝ 1 ~> ℝ 2
     ( path, usedParams ) = pathAndUsedParams @Point id
 
-{-
     curvesI :: 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 'Interval )
     curvesI = brushStrokeData @'Interval @brushParams
                pathI
@@ -492,7 +498,6 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv =
     usedParamsI :: 𝕀ℝ 1 ~> 𝕀 usedParams
     pathI :: 𝕀ℝ 1 ~> 𝕀ℝ 2
     ( pathI, usedParamsI ) = pathAndUsedParams @'Interval singleton
--}
 
     fwdBwd :: OutlineFn
     fwdBwd t
@@ -517,7 +522,14 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv =
                             $ runD ( brushFromParams @Point proxy# id )
                             $ toBrushParams params_t
 
-  in fwdBwd
+    bisSols = bisection 0.0001 curvesI
+
+  in trace
+         ( unlines $
+           ( "bisectionMethod: #(possible zeroes) = " ++ show ( length bisSols ) ) :
+           "" :
+           map show bisSols )
+         fwdBwd
 
 -----------------------------------
 -- Various utility functions
@@ -844,13 +856,14 @@ envelopeEquation :: forall i
                     , Fractional ( I i Double )
                     )
                  => D ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) )
-                 -> ( I i Double, T ( I i ( ℝ 2 ) ), I i Double, I i Double )
+                 -> ( I i Double, T ( I i ( ℝ 2 ) ), T ( I i ( ℝ 2 ) ), I i Double, I 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
-  in ( ee, tot, dEdt, dEds )
+      tot = dcdt -- ^-^ ( dEdt / dEds ) *^ dcds
+      dEdsTot = dEds *^ dcdt ^-^ dEdt *^ dcds
+  in ( ee, tot, dEdsTot, dEdt, dEds )
   -- Computation of total derivative dc/dt:
   --
   -- dc/dt = ∂c/∂t + ∂c/∂s ∂s/∂t
@@ -1072,7 +1085,7 @@ 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 @i dstroke
+          ( ee, dcdt, 𝛿E𝛿sdcdt, 𝛿E𝛿t, 𝛿E𝛿s ) = envelopeEquation @i dstroke
       in  -- trace
           --   ( unlines
           --     [ "envelopeEquation:"
@@ -1089,7 +1102,7 @@ brushStrokeData path params brush =
             { dpath  = dpath_t
             , dbrush = dbrush_t_s
             , dstroke
-            , ee, dcdt, 𝛿E𝛿t, 𝛿E𝛿s }
+            , ee, dcdt, 𝛿E𝛿sdcdt, 𝛿E𝛿t, 𝛿E𝛿s }
 
 
 -- | The value and derivative of a brush stroke at a given coordinate
@@ -1119,9 +1132,83 @@ data StrokeDatum i
     -- \[ \left ( \frac{\mathrm{d} c}{\mathrm{d} t} \right )_{(t_0,s_0)}. \]
     --
     -- This is ill-defined when \( \frac{\partial E}{\partial s} = 0 \).
-  , dcdt :: T ( I i ( ℝ 2 ) )
+  , dcdt, 𝛿E𝛿sdcdt :: T ( I i ( ℝ 2 ) )
+
 
   }
 
 deriving stock instance Show ( StrokeDatum 'Point    )
 deriving stock instance Show ( StrokeDatum 'Interval )
+
+--------------------------------------------------------------------------------
+
+bisection :: Double
+          -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 'Interval ) )
+          -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1, 𝕀 Double, 𝕀ℝ 2 ) ]
+bisection minWidth eqs = bisect initialCands [] []
+  where
+
+    bisect :: [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1, 𝕀 Double, 𝕀ℝ 2 ) ] -- have solutions, need bisection to refine
+           -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] -- have been bisected, don't know if they contain solutions yet
+           -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1, 𝕀 Double, 𝕀ℝ 2 ) ] -- have solutions, don't bisect further
+           -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1, 𝕀 Double, 𝕀ℝ 2 ) ]
+    bisect [] [] sols = sols
+    bisect cands ( ( t, i, s ) : toTry ) sols
+      | Just ( ee, 𝛿E𝛿sdcdt ) <- isCand t i s
+      = bisect ( ( t, i, s, ee, 𝛿E𝛿sdcdt ) : cands ) toTry sols
+      | otherwise
+      = bisect cands toTry sols
+
+    bisect ( cand@( t@( I ( Rounded ( ℝ1 t_lo ) ) ( Rounded ( ℝ1 t_hi ) ) )
+                  , i
+                  , s@( I ( Rounded ( ℝ1 s_lo ) ) ( Rounded ( ℝ1 s_hi ) ) )
+                  , _, _
+                  ) : cands )
+           toTry
+           sols
+      -- If the box is small, don't bisect it further, and store it as a candidate solution.
+      | t_hi - t_lo < minWidth && s_hi - s_lo < minWidth
+      = trace ( "bisection sol: " ++ show cand ++ "\nnbCands = " ++ show ( length cands ) ++ "\nnbToTry = " ++ show ( length toTry ) )
+      $ bisect cands toTry ( cand : sols )
+      -- Otherwise, bisect in its longest direction and add the two resulting
+      -- boxes to the list of boxes to try.
+      | otherwise
+      = let newToTry
+             | t_hi - t_lo > s_hi - s_lo
+             , let t_mid = 0.5 * ( t_lo + t_hi )
+             =   ( I ( Rounded ( ℝ1 t_lo  ) ) ( Rounded ( ℝ1 t_mid ) ), i, s )
+               : ( I ( Rounded ( ℝ1 t_mid ) ) ( Rounded ( ℝ1 t_hi  ) ), i, s )
+               : toTry
+             | let s_mid = 0.5 * ( s_lo + s_hi )
+             =   ( t, i, I ( Rounded ( ℝ1 s_lo  ) ) ( Rounded ( ℝ1 s_mid ) ) )
+               : ( t, i, I ( Rounded ( ℝ1 s_mid ) ) ( Rounded ( ℝ1 s_hi  ) ) )
+               : toTry
+        in bisect cands newToTry sols
+
+    initialCands =
+      getCands
+        ( I ( Rounded $ ℝ1 0 ) ( Rounded $ ℝ1 1 ) )
+        ( I ( Rounded $ ℝ1 0 ) ( Rounded $ ℝ1 1 ) )
+
+    getCands t s =
+      [ (t, i, s, ee, 𝛿E𝛿sdcdt )
+      | let !eqs_t = eqs t
+      , ( eq_t, i ) <- zip ( toList eqs_t ) ( [0,1..] :: [Int] )
+      , let !( StrokeDatum { ee, 𝛿E𝛿sdcdt = T 𝛿E𝛿sdcdt } ) = eq_t s
+      , Interval.inf ee < 0 && Interval.sup ee > 0
+      ,    cmpℝ2 (<) ( getRounded ( Interval.inf 𝛿E𝛿sdcdt ) ) ( ℝ2 0 0 )
+        && cmpℝ2 (>) ( getRounded ( Interval.sup 𝛿E𝛿sdcdt ) ) ( ℝ2 0 0 )
+      ]
+
+    isCand :: 𝕀ℝ 1 -> Int -> 𝕀ℝ 1 -> Maybe ( 𝕀 Double, 𝕀ℝ 2 )
+    isCand t i s = case ( ( eqs t ) `Seq.index` i ) s of
+      StrokeDatum { ee, 𝛿E𝛿sdcdt = T 𝛿E𝛿sdcdt } ->
+        do guard $
+              Interval.inf ee < 0 && Interval.sup ee > 0
+              && cmpℝ2 (<) ( getRounded ( Interval.inf 𝛿E𝛿sdcdt ) ) ( ℝ2 0 0 )
+              && cmpℝ2 (>) ( getRounded ( Interval.sup 𝛿E𝛿sdcdt ) ) ( ℝ2 0 0 )
+           return ( ee, 𝛿E𝛿sdcdt )
+
+cmpℝ2 :: (Double -> Double -> Bool) -> ℝ 2 -> ℝ 2 -> Bool
+cmpℝ2 cmp ( ℝ2 x1 y1 ) ( ℝ2 x2 y2 )
+  = cmp x1 x2 && cmp y1 y2

src/splines/Math/Linear.hs

@@ -10,10 +10,10 @@ module Math.Linear
 
   -- * Points and vectors (second version)
   , ℝ(..), T(.., V2, V3)
-  , Fin(..), eqFin, MFin(..)
+  , Fin(..), MFin(..)
   , Dim, Representable(..), ApRep(..)
   , injection, projection
-  , Vec(..), (!), find
+  , Vec(..), (!), find, zipIndices
 
   -- * Intervals
   , 𝕀, 𝕀ℝ, singleton, nonDecreasing
@@ -27,10 +27,6 @@ import Data.Kind
   ( Type, Constraint )
 import Data.Monoid
   ( Sum(..) )
-import GHC.Exts
-  ( TYPE, RuntimeRep(..)
-  , Word#, plusWord#, minusWord#, isTrue#, eqWord#
-  )
 import GHC.Generics
   ( Generic, Generic1, Generically(..), Generically1(..) )
 import GHC.TypeNats
@@ -99,6 +95,11 @@ data    instance ℝ 3 = ℝ3 {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-#
   deriving anyclass NFData
   deriving stock    ( Show, Eq, Ord )
 
+data    instance ℝ 4 = ℝ4 {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-# UNPACK #-} !Double
+  deriving stock    Generic
+  deriving anyclass NFData
+  deriving stock    ( Show, Eq, Ord )
+
 deriving via ApRep ( Sum Double ) ( ℝ n )
   instance Representable Double ( ℝ n ) => Semigroup ( T ( ℝ n ) )
 deriving via ApRep ( Sum Double ) ( ℝ n )
@@ -151,17 +152,13 @@ pattern V3 x y z = T ( ℝ3 x y z )
 --------------------------------------------------------------------------------
 
 -- | 1, ..., n
-type Fin :: Nat -> TYPE WordRep
-newtype Fin n = Fin Word#
-
-{-# INLINE eqFin #-}
-eqFin :: Fin n -> Fin n -> Bool
-eqFin ( Fin i ) ( Fin j ) = isTrue# ( i `eqWord#` j )
+type Fin :: Nat -> Type
+newtype Fin n = Fin Word
+  deriving stock Eq
 
 -- | 0, ..., n
-type MFin :: Nat -> TYPE WordRep
-newtype MFin n = MFin Word#
-
+type MFin :: Nat -> Type
+newtype MFin n = MFin Word
 
 type Dim :: k -> Nat
 type family Dim v
@@ -181,26 +178,36 @@ instance Representable Double ( ℝ 0 ) where
 
 instance Representable Double ( ℝ 1 ) where
   {-# INLINE tabulate #-}
-  tabulate f = ℝ1 ( f ( Fin 1## ) )
+  tabulate f = ℝ1 ( f ( Fin 1 ) )
   {-# INLINE index #-}
   index ( ℝ1 x ) _ = x
 
 instance Representable Double ( ℝ 2 ) where
   {-# INLINE tabulate #-}
-  tabulate f = ℝ2 ( f ( Fin 1## ) ) ( f ( Fin 2## ) )
+  tabulate f = ℝ2 ( f ( Fin 1 ) ) ( f ( Fin 2 ) )
   {-# INLINE index #-}
   index ( ℝ2 x y ) = \ case
-    Fin 1## -> x
-    _       -> y
+    Fin 1 -> x
+    _     -> y
 
 instance Representable Double ( ℝ 3 ) where
   {-# INLINE tabulate #-}
-  tabulate f = ℝ3 ( f ( Fin 1## ) ) ( f ( Fin 2## ) ) ( f ( Fin 3## ) )
+  tabulate f = ℝ3 ( f ( Fin 1 ) ) ( f ( Fin 2 ) ) ( f ( Fin 3 ) )
   {-# INLINE index #-}
   index ( ℝ3 x y z ) = \ case
-    Fin 1## -> x
-    Fin 2## -> y
-    _       -> z
+    Fin 1 -> x
+    Fin 2 -> y
+    _     -> z
+
+instance Representable Double ( ℝ 4 ) where
+  {-# INLINE tabulate #-}
+  tabulate f = ℝ4 ( f ( Fin 1 ) ) ( f ( Fin 2 ) ) ( f ( Fin 3 ) ) ( f ( Fin 4 ) )
+  {-# INLINE index #-}
+  index ( ℝ4 x y z w ) = \ case
+    Fin 1 -> x
+    Fin 2 -> y
+    Fin 3 -> z
+    _     -> w
 
 {-# INLINE projection #-}
 projection :: ( Representable r u, Representable r v )
@@ -215,30 +222,42 @@ injection :: ( Representable r u, Representable r v )
           -> u -> v -> v
 injection f = \ u v ->
   tabulate \ i -> case f i of
-    MFin 0## -> index v i
-    MFin j   -> index u ( Fin j )
+    MFin 0 -> index v i
+    MFin j -> index u ( Fin j )
 
-type Vec :: Nat -> TYPE WordRep -> Type
+infixr 5 `VS`
+type Vec :: Nat -> Type -> Type
 data Vec n a where
   VZ :: Vec 0 a
   VS :: a -> Vec n a -> Vec ( 1 + n ) a
 
+deriving stock instance Functor     ( Vec n )
+deriving stock instance Foldable    ( Vec n )
+deriving stock instance Traversable ( Vec n )
+
 infixl 9 !
 (!) :: forall l a. Vec l a -> Fin l -> a
-VS a _ ! Fin 1##  = a
-VS _ a ! Fin i    = a ! Fin ( i `minusWord#` 1## )
-_      ! _        = error "impossible: Fin 0 is uninhabited"
+VS a _ ! Fin 1 = a
+VS _ a ! Fin i = a ! Fin ( i - 1 )
+_      ! _     = error "impossible: Fin 0 is uninhabited"
 
 find :: forall l a. ( a -> a -> Bool ) -> Vec l a -> a -> MFin l
-find eq v b = MFin ( go 1## v )
+find eq v b = MFin ( go 1 v )
   where
-    go :: Word# -> Vec n a -> Word#
+    go :: Word -> Vec n a -> Word
     go j ( VS a as )
       | a `eq` b
       = j
       | otherwise
-      = go ( j `plusWord#` 1## ) as
-    go _ VZ = 0##
+      = go ( j + 1 ) as
+    go _ VZ = 0
+
+zipIndices :: forall n a. Vec n a -> [ ( Fin n, a ) ]
+zipIndices = go 0
+  where
+    go :: forall i. Word -> Vec i a -> [ ( Fin n, a ) ]
+    go _ VZ = []
+    go i (a `VS` as) = ( Fin i, a ) : go ( i + 1 ) as
 
 --------------------------------------------------------------------------------
 -- Instances in terms of representable.