hashirama/metabrush
1
0
Fork 0
mirror of https://gitlab.com/sheaf/metabrush.git synced 2024-11-30 19:04:06 +00:00
Code Issues Actions Projects Releases Wiki Activity

Add tests and fix MonomialBasis D3A3 instance

Browse source
This commit is contained in:
sheaf 2024-03-08 15:35:39 +01:00
parent cebfeb0b7a
commit 34c129d72a
8 changed files with 532 additions and 150 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

32
brush-strokes/brush-strokes.cabal
View file

@ -205,10 +205,38 @@ executable inspect
Math.Interval.Abstract Math.Interval.Abstract
build-depends: build-depends:
brush-strokes, brush-strokes
data-reify , data-reify
^>= 0.6.3 ^>= 0.6.3
test-suite tests
import:
common
type:
exitcode-stdio-1.0
hs-source-dirs:
src/test
default-language:
Haskell2010
main-is:
Main.hs
build-depends:
brush-strokes
, falsify
^>= 0.2
, hspray
^>= 0.1.3
, tasty
^>= 1.5
, unordered-containers
>= 0.2.15 && < 0.3
benchmark cusps benchmark cusps
import: import:

3
brush-strokes/cabal.project
View file

@ -4,10 +4,13 @@ constraints:
acts -finitary, acts -finitary,
rounded-hw -pure-hs -c99 -avx512 +ghc-prim -x87-long-double rounded-hw -pure-hs -c99 -avx512 +ghc-prim -x87-long-double
tests: True
allow-newer: allow-newer:
acts:base, acts:deepseq, acts:base, acts:deepseq,
groups-generic:base, groups-generic:base,
eigen:primitive, eigen:primitive,
falsify:base, falsify:tasty,
--package brush-strokes --package brush-strokes
profiling: True profiling: True

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

@ -608,3 +608,34 @@ x2 = 𝕀 x2_lo x2_hi
( b1 - a12 * x2 ) `extendedDivide` a11 ( b1 - a12 * x2 ) `extendedDivide` a11
-} -}
{-
SMOKING GUN:
> (f, fI) = brushStrokeFunctions $ ellipseBrushStroke (0,1) pi
> (t, i, s) = mkBox (0.54, 0.55) 2 (0.5480, 0.5481)
> _D32_dxdxdx $ dbrush $ eval fI (t,i,s)
> [ℝ2 54187.61174031432 -11626.47655724034, ℝ2 55257.54384135226 -9759.09092706883]
In Mathematica:
FindMinimum[{D[b[vt, s, myPb0, myPb1][[1]], {vt, 3}] /. (vt -> t),
0.54 <= t <= 0.55 && 0.548 <= s <= 0.5481}, {{t, 0.545}, {s,
0.54805}}]
> { -503.445, {t -> 0.54, s -> 0.548} }
FindMaximum[{D[b[vt, s, myPb0, myPb1][[1]], {vt, 3}] /. (vt -> t),
0.54 <= t <= 0.55 && 0.548 <= s <= 0.5481}, {{t, 0.545}, {s,
0.54805}}]
> { -480.086, {t -> 0.55, s -> 0.5481} }
so the actual range is [-503.445,-480.086] but we have computed a range of
[54187.61174031432, 55257.54384135226] which does not contain the actual range!
-}
{-
Math.Ring.cos ( linearD @_ @3 (\i -> 𝕀 (getR1 $ inf i) (getR1 $ sup i)) (𝕀 (ℝ1 $ 0.548 * pi) (ℝ1 $ 0.5481 * pi)) :: D3𝔸1 ( 𝕀 Double ) )
-}

14
brush-strokes/src/cusps/inspect/Main.hs
View file

@ -140,8 +140,8 @@ boo x i y z
ellipseTest :: StrokeDatum 3 AI ellipseTest :: StrokeDatum 3 AI
ellipseTest = ellipseTest =
boo ( ellipseBrush ( Pt . Val ) scaleI ) 1 boo ( ellipseBrush ( Pt . Val ) scaleI ) 1
( mkPoint ( ℝ2 0 0 ) 1.0111 1.0222 0 ) ( mkPoint ( ℝ2 0 0 ) 10 25 0 )
( LineTo { curveEnd = NextPoint $ mkPoint ( ℝ2 100 0 ) 10.0333 10.0444 ( pi Prelude./ 2 ) ( LineTo { curveEnd = NextPoint $ mkPoint ( ℝ2 100 0 ) 15 40 pi
, curveData = () } ) , curveData = () } )
where where
mkPoint :: ℝ 2 -> Double -> Double -> Double -> PointData ( ℝ 3 ) mkPoint :: ℝ 2 -> Double -> Double -> Double -> PointData ( ℝ 3 )
@ -258,7 +258,7 @@ instance HasBézier 3 AI where
( Cubic.bezier''' bez ) ( Cubic.bezier''' bez )
-------------------------------------------------------------------------------- --------------------------------------------------------------------------------
-- Brushes (TODO copied from MetaBrush.Asset.Brushes) -- Brushes (TODO copied from Calligraphy.Brushes)
κ :: Double κ :: Double
κ = 0.5519150244935105707435627227925 κ = 0.5519150244935105707435627227925
@ -304,13 +304,13 @@ ellipseBrush mkI scaleI =
mkPt x y mkPt x y
= let !x' = a `scaledBy` x = let !x' = a `scaledBy` x
!y' = b `scaledBy` y !y' = b `scaledBy` y
{- -- {-
!c = cos phi !c = cos phi
!s = sin phi !s = sin phi
in ( x' * c - y' * s ) *^ e_x in ( x' * c - y' * s ) *^ e_x
^+^ ( x' * s + y' * c ) *^ e_y ^+^ ( x' * s + y' * c ) *^ e_y
-} -- -}
-- {- {-
!r = sqrt $ x' ^ 2 + y' ^ 2 !r = sqrt $ x' ^ 2 + y' ^ 2
!arctgt = atan ( y' / x' ) !arctgt = atan ( y' / x' )
-- a and b are always strictly positive, so we can compute -- a and b are always strictly positive, so we can compute
@ -326,7 +326,7 @@ ellipseBrush mkI scaleI =
in in
( r * cos phi' ) *^ e_x ( r * cos phi' ) *^ e_x
^+^ ( r * sin phi' ) *^ e_y ^+^ ( r * sin phi' ) *^ e_y
-- -} -}
in circleSpline @AI @k @( ℝ 3 ) @( ℝ 2 ) mkPt in circleSpline @AI @k @( ℝ 3 ) @( ℝ 2 ) mkPt
where where

231
brush-strokes/src/lib/Math/Algebra/Dual/Internal.hs
View file

@ -244,7 +244,7 @@ chainRule1NQ zero_w sum_w scale_w df dg =
] ]
) ||] ) ||]
-- | The chain rule for a composition \( \mathbb{R}^n \to \mathbb{R} \to \mathbb{R} \) -- | The chain rule for a composition \( \mathbb{R}^n \to \mathbb{R}^1 \to \mathbb{R}^1 \)
-- --
-- (To be spliced in using Template Haskell.) -- (To be spliced in using Template Haskell.)
chainRuleN1Q :: forall du dr1 r chainRuleN1Q :: forall du dr1 r
@ -291,7 +291,7 @@ type instance Vars D𝔸0 = 0
instance MonomialBasis D𝔸0 where instance MonomialBasis D𝔸0 where
monTabulate f = monTabulate f =
[|| let [|| let
!_D0_v = $$( f $ Mon VZ ) _D0_v = $$( f $ Mon VZ )
in D0 { .. } in D0 { .. }
||] ||]
@ -302,8 +302,8 @@ type instance Vars D1𝔸1 = 1
instance MonomialBasis D1𝔸1 where instance MonomialBasis D1𝔸1 where
monTabulate f = monTabulate f =
[|| let [|| let
!_D11_v = $$( f $ Mon ( 0 `VS` VZ ) ) _D11_v = $$( f $ Mon ( 0 `VS` VZ ) )
!_D11_dx = T $$( f $ Mon ( 1 `VS` VZ ) ) _D11_dx = T $$( f $ Mon ( 1 `VS` VZ ) )
in D11 { .. } in D11 { .. }
||] ||]
@ -316,9 +316,9 @@ type instance Vars D2𝔸1 = 1
instance MonomialBasis D2𝔸1 where instance MonomialBasis D2𝔸1 where
monTabulate f = monTabulate f =
[|| let [|| let
!_D21_v = $$( f $ Mon ( 0 `VS` VZ ) ) _D21_v = $$( f $ Mon ( 0 `VS` VZ ) )
!_D21_dx = T $$( f $ Mon ( 1 `VS` VZ ) ) _D21_dx = T $$( f $ Mon ( 1 `VS` VZ ) )
!_D21_dxdx = T $$( f $ Mon ( 2 `VS` VZ ) ) _D21_dxdx = T $$( f $ Mon ( 2 `VS` VZ ) )
in D21 { .. } in D21 { .. }
||] ||]
@ -332,10 +332,10 @@ type instance Vars D3𝔸1 = 1
instance MonomialBasis D3𝔸1 where instance MonomialBasis D3𝔸1 where
monTabulate f = monTabulate f =
[|| let [|| let
!_D31_v = $$( f $ Mon ( 0 `VS` VZ ) ) _D31_v = $$( f $ Mon ( 0 `VS` VZ ) )
!_D31_dx = T $$( f $ Mon ( 1 `VS` VZ ) ) _D31_dx = T $$( f $ Mon ( 1 `VS` VZ ) )
!_D31_dxdx = T $$( f $ Mon ( 2 `VS` VZ ) ) _D31_dxdx = T $$( f $ Mon ( 2 `VS` VZ ) )
!_D31_dxdxdx = T $$( f $ Mon ( 3 `VS` VZ ) ) _D31_dxdxdx = T $$( f $ Mon ( 3 `VS` VZ ) )
in D31 { .. } in D31 { .. }
||] ||]
@ -350,9 +350,9 @@ type instance Vars D1𝔸2 = 2
instance MonomialBasis D1𝔸2 where instance MonomialBasis D1𝔸2 where
monTabulate f = monTabulate f =
[|| let [|| let
!_D12_v = $$( f $ Mon ( 0 `VS` 0 `VS` VZ ) ) _D12_v = $$( f $ Mon ( 0 `VS` 0 `VS` VZ ) )
!_D12_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` VZ ) ) _D12_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` VZ ) )
!_D12_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` VZ ) ) _D12_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` VZ ) )
in D12 { .. } in D12 { .. }
||] ||]
@ -366,12 +366,12 @@ type instance Vars D2𝔸2 = 2
instance MonomialBasis D2𝔸2 where instance MonomialBasis D2𝔸2 where
monTabulate f = monTabulate f =
[|| let [|| let
!_D22_v = $$( f $ Mon ( 0 `VS` 0 `VS` VZ ) ) _D22_v = $$( f $ Mon ( 0 `VS` 0 `VS` VZ ) )
!_D22_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` VZ ) ) _D22_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` VZ ) )
!_D22_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` VZ ) ) _D22_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` VZ ) )
!_D22_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` VZ ) ) _D22_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` VZ ) )
!_D22_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` VZ ) ) _D22_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` VZ ) )
!_D22_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` VZ ) ) _D22_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` VZ ) )
in D22 { .. } in D22 { .. }
||] ||]
@ -388,16 +388,16 @@ type instance Vars D3𝔸2 = 2
instance MonomialBasis D3𝔸2 where instance MonomialBasis D3𝔸2 where
monTabulate f = monTabulate f =
[|| let [|| let
!_D32_v = $$( f $ Mon ( 0 `VS` 0 `VS` VZ ) ) _D32_v = $$( f $ Mon ( 0 `VS` 0 `VS` VZ ) )
!_D32_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` VZ ) ) _D32_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` VZ ) )
!_D32_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` VZ ) ) _D32_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` VZ ) )
!_D32_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` VZ ) ) _D32_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` VZ ) )
!_D32_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` VZ ) ) _D32_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` VZ ) )
!_D32_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` VZ ) ) _D32_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` VZ ) )
!_D32_dxdxdx = T $$( f $ Mon ( 3 `VS` 0 `VS` VZ ) ) _D32_dxdxdx = T $$( f $ Mon ( 3 `VS` 0 `VS` VZ ) )
!_D32_dxdxdy = T $$( f $ Mon ( 2 `VS` 1 `VS` VZ ) ) _D32_dxdxdy = T $$( f $ Mon ( 2 `VS` 1 `VS` VZ ) )
!_D32_dxdydy = T $$( f $ Mon ( 1 `VS` 2 `VS` VZ ) ) _D32_dxdydy = T $$( f $ Mon ( 1 `VS` 2 `VS` VZ ) )
!_D32_dydydy = T $$( f $ Mon ( 0 `VS` 3 `VS` VZ ) ) _D32_dydydy = T $$( f $ Mon ( 0 `VS` 3 `VS` VZ ) )
in D32 { .. } ||] in D32 { .. } ||]
monIndex d = \ case monIndex d = \ case
@ -435,16 +435,16 @@ type instance Vars D2𝔸3 = 3
instance MonomialBasis D2𝔸3 where instance MonomialBasis D2𝔸3 where
monTabulate f = monTabulate f =
[|| let [|| let
!_D23_v = $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` VZ ) ) _D23_v = $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` VZ ) )
!_D23_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) ) _D23_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) )
!_D23_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) ) _D23_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) )
!_D23_dz = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) ) _D23_dz = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) )
!_D23_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) ) _D23_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) )
!_D23_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) ) _D23_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) )
!_D23_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) ) _D23_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) )
!_D23_dxdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) ) _D23_dxdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) )
!_D23_dydz = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) ) _D23_dydz = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) )
!_D23_dzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) ) _D23_dzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) )
in D23 { .. } in D23 { .. }
||] ||]
@ -466,26 +466,26 @@ type instance Vars D3𝔸3 = 3
instance MonomialBasis D3𝔸3 where instance MonomialBasis D3𝔸3 where
monTabulate f = monTabulate f =
[|| let [|| let
!_D33_v = $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` VZ ) ) _D33_v = $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` VZ ) )
!_D33_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) ) _D33_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) )
!_D33_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) ) _D33_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) )
!_D33_dz = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) ) _D33_dz = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) )
!_D33_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) ) _D33_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) )
!_D33_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) ) _D33_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) )
!_D33_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) ) _D33_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) )
!_D33_dxdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) ) _D33_dxdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) )
!_D33_dydz = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) ) _D33_dydz = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) )
!_D33_dzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) ) _D33_dzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) )
!_D33_dxdxdx = T $$( f $ Mon ( 3 `VS` 0 `VS` 0 `VS` VZ ) ) _D33_dxdxdx = T $$( f $ Mon ( 3 `VS` 0 `VS` 0 `VS` VZ ) )
!_D33_dxdxdy = T $$( f $ Mon ( 2 `VS` 1 `VS` 0 `VS` VZ ) ) _D33_dxdxdy = T $$( f $ Mon ( 2 `VS` 1 `VS` 0 `VS` VZ ) )
!_D33_dxdydy = T $$( f $ Mon ( 1 `VS` 2 `VS` 0 `VS` VZ ) ) _D33_dxdydy = T $$( f $ Mon ( 1 `VS` 2 `VS` 0 `VS` VZ ) )
!_D33_dydydy = T $$( f $ Mon ( 0 `VS` 3 `VS` 0 `VS` VZ ) ) _D33_dydydy = T $$( f $ Mon ( 0 `VS` 3 `VS` 0 `VS` VZ ) )
!_D33_dxdxdz = T $$( f $ Mon ( 2 `VS` 0 `VS` 1 `VS` VZ ) ) _D33_dxdxdz = T $$( f $ Mon ( 2 `VS` 0 `VS` 1 `VS` VZ ) )
!_D33_dxdydz = T $$( f $ Mon ( 1 `VS` 1 `VS` 1 `VS` VZ ) ) _D33_dxdydz = T $$( f $ Mon ( 1 `VS` 1 `VS` 1 `VS` VZ ) )
!_D33_dxdzdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 2 `VS` VZ ) ) _D33_dxdzdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 2 `VS` VZ ) )
!_D33_dydydz = T $$( f $ Mon ( 0 `VS` 2 `VS` 1 `VS` VZ ) ) _D33_dydydz = T $$( f $ Mon ( 0 `VS` 2 `VS` 1 `VS` VZ ) )
!_D33_dydzdz = T $$( f $ Mon ( 0 `VS` 1 `VS` 2 `VS` VZ ) ) _D33_dydzdz = T $$( f $ Mon ( 0 `VS` 1 `VS` 2 `VS` VZ ) )
!_D33_dzdzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 3 `VS` VZ ) ) _D33_dzdzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 3 `VS` VZ ) )
in D33 { .. } ||] in D33 { .. } ||]
monIndex d = \ case monIndex d = \ case
@ -505,6 +505,7 @@ instance MonomialBasis D3𝔸3 where
Mon ( 2 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dxdxdz $$d ||] Mon ( 2 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dxdxdz $$d ||]
Mon ( 1 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dxdydz $$d ||] Mon ( 1 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dxdydz $$d ||]
Mon ( 1 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D33_dxdzdz $$d ||] Mon ( 1 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D33_dxdzdz $$d ||]
Mon ( 0 `VS` 2 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dydydz $$d ||]
Mon ( 0 `VS` 1 `VS` 2 `VS` VZ ) -> [|| unT $ _D33_dydzdz $$d ||] Mon ( 0 `VS` 1 `VS` 2 `VS` VZ ) -> [|| unT $ _D33_dydzdz $$d ||]
Mon ( 0 `VS` 0 `VS` 3 `VS` VZ ) -> [|| unT $ _D33_dzdzdz $$d ||] Mon ( 0 `VS` 0 `VS` 3 `VS` VZ ) -> [|| unT $ _D33_dzdzdz $$d ||]
_ -> [|| _D33_v $$d ||] _ -> [|| _D33_v $$d ||]
@ -514,11 +515,11 @@ type instance Vars D1𝔸4 = 4
instance MonomialBasis D1𝔸4 where instance MonomialBasis D1𝔸4 where
monTabulate f = monTabulate f =
[|| let [|| let
!_D14_v = $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) ) _D14_v = $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
!_D14_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) ) _D14_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
!_D14_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) ) _D14_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) )
!_D14_dz = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) ) _D14_dz = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) )
!_D14_dw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) ) _D14_dw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) )
in D14 { .. } ||] in D14 { .. } ||]
monIndex d = \ case monIndex d = \ case
@ -533,21 +534,21 @@ type instance Vars D2𝔸4 = 4
instance MonomialBasis D2𝔸4 where instance MonomialBasis D2𝔸4 where
monTabulate f = monTabulate f =
[|| let [|| let
!_D24_v = $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) ) _D24_v = $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
!_D24_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) ) _D24_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
!_D24_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) ) _D24_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) )
!_D24_dz = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) ) _D24_dz = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) )
!_D24_dw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) ) _D24_dw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) )
!_D24_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) ) _D24_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
!_D24_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) ) _D24_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) )
!_D24_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) ) _D24_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) )
!_D24_dxdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) ) _D24_dxdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) )
!_D24_dydz = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) ) _D24_dydz = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) )
!_D24_dzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) ) _D24_dzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) )
!_D24_dxdw = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) ) _D24_dxdw = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) )
!_D24_dydw = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) ) _D24_dydw = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) )
!_D24_dzdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) ) _D24_dzdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) )
!_D24_dwdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) ) _D24_dwdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) )
in D24 { .. } ||] in D24 { .. } ||]
monIndex d = \ case monIndex d = \ case
@ -573,41 +574,41 @@ type instance Vars D3𝔸4 = 4
instance MonomialBasis D3𝔸4 where instance MonomialBasis D3𝔸4 where
monTabulate f = monTabulate f =
[|| let [|| let
!_D34_v = $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) ) _D34_v = $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
!_D34_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) ) _D34_dx = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
!_D34_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) ) _D34_dy = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) )
!_D34_dz = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) ) _D34_dz = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) )
!_D34_dw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) ) _D34_dw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) )
!_D34_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) ) _D34_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
!_D34_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) ) _D34_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) )
!_D34_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) ) _D34_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) )
!_D34_dxdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) ) _D34_dxdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) )
!_D34_dydz = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) ) _D34_dydz = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) )
!_D34_dzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) ) _D34_dzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) )
!_D34_dxdw = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) ) _D34_dxdw = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) )
!_D34_dydw = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) ) _D34_dydw = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) )
!_D34_dzdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) ) _D34_dzdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) )
!_D34_dwdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) ) _D34_dwdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) )
!_D34_dxdxdx = T $$( f $ Mon ( 3 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) ) _D34_dxdxdx = T $$( f $ Mon ( 3 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) )
!_D34_dxdxdy = T $$( f $ Mon ( 2 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) ) _D34_dxdxdy = T $$( f $ Mon ( 2 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) )
!_D34_dxdydy = T $$( f $ Mon ( 1 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) ) _D34_dxdydy = T $$( f $ Mon ( 1 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) )
!_D34_dydydy = T $$( f $ Mon ( 0 `VS` 3 `VS` 0 `VS` 0 `VS` VZ ) ) _D34_dydydy = T $$( f $ Mon ( 0 `VS` 3 `VS` 0 `VS` 0 `VS` VZ ) )
!_D34_dxdxdz = T $$( f $ Mon ( 2 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) ) _D34_dxdxdz = T $$( f $ Mon ( 2 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) )
!_D34_dxdydz = T $$( f $ Mon ( 1 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) ) _D34_dxdydz = T $$( f $ Mon ( 1 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) )
!_D34_dxdzdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) ) _D34_dxdzdz = T $$( f $ Mon ( 1 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) )
!_D34_dydydz = T $$( f $ Mon ( 0 `VS` 2 `VS` 1 `VS` 0 `VS` VZ ) ) _D34_dydydz = T $$( f $ Mon ( 0 `VS` 2 `VS` 1 `VS` 0 `VS` VZ ) )
!_D34_dydzdz = T $$( f $ Mon ( 0 `VS` 1 `VS` 2 `VS` 0 `VS` VZ ) ) _D34_dydzdz = T $$( f $ Mon ( 0 `VS` 1 `VS` 2 `VS` 0 `VS` VZ ) )
!_D34_dzdzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 3 `VS` 0 `VS` VZ ) ) _D34_dzdzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 3 `VS` 0 `VS` VZ ) )
!_D34_dxdxdw = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) ) _D34_dxdxdw = T $$( f $ Mon ( 2 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) )
!_D34_dxdydw = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) ) _D34_dxdydw = T $$( f $ Mon ( 1 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) )
!_D34_dydydw = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 1 `VS` VZ ) ) _D34_dydydw = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `VS` 1 `VS` VZ ) )
!_D34_dxdzdw = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) ) _D34_dxdzdw = T $$( f $ Mon ( 1 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) )
!_D34_dydzdw = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 1 `VS` VZ ) ) _D34_dydzdw = T $$( f $ Mon ( 0 `VS` 1 `VS` 1 `VS` 1 `VS` VZ ) )
!_D34_dzdzdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 1 `VS` VZ ) ) _D34_dzdzdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` 1 `VS` VZ ) )
!_D34_dxdwdw = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) ) _D34_dxdwdw = T $$( f $ Mon ( 1 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) )
!_D34_dydwdw = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 2 `VS` VZ ) ) _D34_dydwdw = T $$( f $ Mon ( 0 `VS` 1 `VS` 0 `VS` 2 `VS` VZ ) )
!_D34_dzdwdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 2 `VS` VZ ) ) _D34_dzdwdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `VS` 2 `VS` VZ ) )
!_D34_dwdwdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 3 `VS` VZ ) ) _D34_dwdwdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 3 `VS` VZ ) )
in D34 { .. } ||] in D34 { .. } ||]
monIndex d = \ case monIndex d = \ case

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

@ -180,10 +180,10 @@ instance HasEnvelopeEquation 2 where
-- ∂²c/∂s² = ∂²b/∂s² -- ∂²c/∂s² = ∂²b/∂s²
envelopeEquation co ( D22 _ c_t c_s c_tt c_ts c_ss ) = envelopeEquation co ( D22 _ c_t c_s c_tt c_ts c_ss ) =
let !ee = c_t × c_s let ee = c_t × c_s
!ee_t = c_tt × c_s + c_t × c_ts ee_t = c_tt × c_s + c_t × c_ts
!ee_s = c_ts × c_s + c_t × c_ss ee_s = c_ts × c_s + c_t × c_ss
!𝛿E𝛿sdcdt = ee_s *^ c_t ^-^ ee_t *^ c_s 𝛿E𝛿sdcdt = ee_s *^ c_t ^-^ ee_t *^ c_s
-- TODO: we get c_t * c_t and c_s * c_s terms... -- TODO: we get c_t * c_t and c_s * c_s terms...
-- These could be squares (better with interval arithmetic)? -- These could be squares (better with interval arithmetic)?
in ( D12 ( co ee ) ( T $ co ee_t ) ( T $ co ee_s ) in ( D12 ( co ee ) ( T $ co ee_t ) ( T $ co ee_s )
@ -235,23 +235,23 @@ instance HasEnvelopeEquation 3 where
( D32 _ c_t c_s ( D32 _ c_t c_s
c_tt c_ts c_ss c_tt c_ts c_ss
c_ttt c_tts c_tss c_sss ) c_ttt c_tts c_tss c_sss )
= let !ee = c_t × c_s = let ee = c_t × c_s
!ee_t = c_tt × c_s + c_t × c_ts ee_t = c_tt × c_s + c_t × c_ts
!ee_s = c_ts × c_s + c_t × c_ss ee_s = c_ts × c_s + c_t × c_ss
!ee_tt = c_ttt × c_s ee_tt = c_ttt × c_s
+ c_tt × c_ts * 2 + c_tt × c_ts * 2
+ c_t × c_tts + c_t × c_tts
!ee_ts = c_tts × c_s ee_ts = c_tts × c_s
+ c_tt × c_ss + c_tt × c_ss
-- + c_ts × c_ts -- cancels out -- + c_ts × c_ts -- cancels out
+ c_t × c_tss + c_t × c_tss
!ee_ss = c_tss × c_s ee_ss = c_tss × c_s
+ c_ts × c_ss * 2 + c_ts × c_ss * 2
+ c_t × c_sss + c_t × c_sss
!𝛿E𝛿sdcdt = ee_s *^ c_t ^-^ ee_t *^ c_s 𝛿E𝛿sdcdt = ee_s *^ c_t ^-^ ee_t *^ c_s
!𝛿E𝛿sdcdt_t = ee_ts *^ c_t ^+^ ee_s *^ c_tt 𝛿E𝛿sdcdt_t = ee_ts *^ c_t ^+^ ee_s *^ c_tt
^-^ ( ee_tt *^ c_s ^+^ ee_t *^ c_ts ) ^-^ ( ee_tt *^ c_s ^+^ ee_t *^ c_ts )
!𝛿E𝛿sdcdt_s = ee_ss *^ c_t ^+^ ee_s *^ c_ts 𝛿E𝛿sdcdt_s = ee_ss *^ c_t ^+^ ee_s *^ c_ts
^-^ ( ee_ts *^ c_s ^+^ ee_t *^ c_ss ) ^-^ ( ee_ts *^ c_s ^+^ ee_t *^ c_ss )
in ( D22 in ( D22
( co ee ) ( co ee )

8
brush-strokes/src/lib/Math/Linear.hs
View file

@ -122,7 +122,9 @@ data Vec n a where
-- can't be strict, otherwise we can't conveniently -- can't be strict, otherwise we can't conveniently
-- unsafeCoerce from lists -- unsafeCoerce from lists
deriving stock instance Show a => Show ( Vec n a ) --deriving stock instance Show a => Show ( Vec n a )
instance Show a => Show ( Vec n a ) where
showsPrec p v = showsPrec p ( unsafeCoerce v :: [ a ] )
deriving stock instance Functor ( Vec n ) deriving stock instance Functor ( Vec n )
deriving stock instance Foldable ( Vec n ) deriving stock instance Foldable ( Vec n )
@ -137,9 +139,7 @@ instance Ord a => Ord ( Vec n a ) where
infixl 9 ! infixl 9 !
(!) :: forall l a. Vec l a -> Fin l -> a (!) :: forall l a. Vec l a -> Fin l -> a
VS a _ ! Fin 1 = a v ! Fin i = ( unsafeCoerce v :: [ a ] ) !! fromIntegral i
VS _ a ! Fin i = a ! Fin ( i - 1 )
_ ! _ = error "impossible: Fin 0 is uninhabited"
find :: forall l a. ( a -> Bool ) -> Vec l a -> MFin l find :: forall l a. ( a -> Bool ) -> Vec l a -> MFin l
find f v = MFin ( find_ 1 v ) find f v = MFin ( find_ 1 v )

319
brush-strokes/src/test/Main.hs Normal file
View file

@ -0,0 +1,319 @@
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE NumericUnderscores #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell #-}
module Main (main) where
-- base
import Prelude hiding
( Num(..), (^) )
import Data.Foldable
( toList )
import Data.List.NonEmpty
( NonEmpty(..) )
import Data.Maybe
( catMaybes )
import Data.Traversable
( for )
import Unsafe.Coerce
( unsafeCoerce )
-- brush-strokes
import Math.Algebra.Dual
import Math.Linear
import Math.Module
import Math.Monomial
( multiSubsetSum, multiSubsetsSum
, MonomialBasis ( monTabulate, monIndex )
)
import Math.Ring
-- hspray
import Math.Algebra.Hspray
( Spray )
import qualified Math.Algebra.Hspray as Spray
-- falsify
import Test.Tasty.Falsify
import qualified Test.Falsify.Generator as Falsify
( Gen )
import qualified Test.Falsify.Generator as Falsify.Gen
import Test.Falsify.Predicate
( (.$) )
import qualified Test.Falsify.Predicate as Falsify.Prop
import qualified Test.Falsify.Property as Falsify
( Property
, assert
, discard
, gen, genWith
)
import qualified Test.Falsify.Range as Falsify
-- tasty
import qualified Test.Tasty as Tasty
-- unordered-containers
import qualified Data.HashMap.Lazy as HashMap
--------------------------------------------------------------------------------
main :: IO ()
main =
Tasty.defaultMain $
Tasty.testGroup "brush-strokes property tests"
[ Tasty.testGroup "Automatic differentiation"
[ Tasty.testGroup "Monomial basis"
[ testProperty "Round trip D33" testMonomialBasisD33
]
, Tasty.testGroup "Monomials"
[ Tasty.testGroup "multiSubsetSum"
[ testProperty "multiSubsetSum valid" testMultiSubsetSumValid
, testProperty "multiSubsetSum exhaustive" testMultiSubsetSumExhaustive
]
-- , Tasty.testGroup "multiSubsetsSum"
-- [ testProperty "multiSubsetsSum exhaustive" testMultiSubsetsSumExhaustive
-- ]
]
, Tasty.testGroup "chainRule1NQ"
[ testProperty "chainRule1NQ_1" testChainRule1NQ_1
, testProperty "chainRule1NQ_2" testChainRule1NQ_2
, testProperty "chainRule1NQ_3" testChainRule1NQ_3
]
]
]
-- | Check that the 'multiSubsetSum' function returns valid answers, i.e.
-- all returned multisubsets have the desired size and sum.
testMultiSubsetSumValid :: Falsify.Property ()
testMultiSubsetSumValid = do
rg <- Falsify.genWith (\ rg -> Just $ "range = " ++ show rg ) $ Falsify.Gen.inRange $ Falsify.between ( 1, 6 )
sz <- Falsify.genWith (\ sz -> Just $ "size = " ++ show sz ) $ Falsify.Gen.inRange $ Falsify.between ( 0, 20 )
tot <- Falsify.genWith (\ tot -> Just $ "tot = " ++ show tot) $ Falsify.Gen.inRange $ Falsify.between ( sz, sz * rg )
let range = [ 1 .. rg ]
mss = multiSubsetSum sz tot range
case mss of
[] -> Falsify.discard
r:rs -> do
ms <- Falsify.gen $ Falsify.Gen.elem ( r :| rs )
Falsify.assert
$ Falsify.Prop.eq
.$ ("(sz, tot)", (sz, tot) )
.$ ("computed (sz, tot)", (size ms, total ms))
where
size, total :: [ ( Word, Word ) ] -> Word
size [] = 0
size ((_,n):ins) = n + size ins
total [] = 0
total ((i,n):ins) = i * n + total ins
-- | Check that the 'multiSubsetSum' function returns all multisubsets of
-- the given set, by generating a random multisubset, computing its size, and
-- checking it belongs to the output of the 'multiSubsetSum' function.
testMultiSubsetSumExhaustive :: Falsify.Property ()
testMultiSubsetSumExhaustive = do
rg <- Falsify.genWith (\ rg -> Just $ "range = " ++ show rg) $ Falsify.Gen.inRange $ Falsify.between ( 1, 6 )
sz <- Falsify.genWith (\ sz -> Just $ "size = " ++ show sz) $ Falsify.Gen.inRange $ Falsify.between ( 0, 10 )
let range = [ 1 .. rg ]
(multiSubset, tot) <- Falsify.genWith (\ ms -> Just $ "multisubset = " ++ show ms) $ genMultiSubset range sz
Falsify.assert
$ Falsify.Prop.elem
.$ ("all multisubsets", multiSubsetSum sz tot range )
.$ ("random multisubset", multiSubset)
genMultiSubset :: [ Word ] -> Word -> Falsify.Gen ( [ ( Word, Word ) ] , Word )
genMultiSubset [i] sz =
return $
if sz == 0
then ( [], 0 )
else ( [ ( i, sz ) ], i * sz )
genMultiSubset (i:is) sz = do
nb <- Falsify.Gen.inRange $ Falsify.between ( 0, sz )
(rest, tot) <- genMultiSubset is ( sz - nb )
return $ ( if nb == 0 then rest else ( i, nb ) : rest, tot + nb * i )
genMultiSubset [] _ = error "impossible"
coerceVec1 :: [ a ] -> Vec n a
coerceVec1 = unsafeCoerce
coerceVec2 :: Vec n a -> [ a ]
coerceVec2 = toList
-- | Check that the 'multiSubsetSums' function returns all collections of
-- multisubsets of the given set (see 'testMultiSubsetSumExhaustive').
testMultiSubsetsSumExhaustive :: Falsify.Property ()
testMultiSubsetsSumExhaustive = do
rg <- Falsify.genWith (\ rg -> Just $ "range = " ++ show rg) $ Falsify.Gen.inRange $ Falsify.between ( 1, 5 )
let range = [ 1 .. rg ]
n <- Falsify.genWith (\ n -> Just $ "n = " ++ show n ) $ Falsify.Gen.inRange $ Falsify.between ( 1, 10 )
multiSubsets <- for ( [ 0 .. n - 1 ] :: [ Word ] ) \ i -> do
sz <- Falsify.gen $ Falsify.Gen.inRange $ Falsify.between ( 0, 5 )
( ms, tot ) <- Falsify.genWith ( \ ms -> Just $ "ms_" ++ show i ++ " = " ++ show ms ) $ genMultiSubset range sz
return ( ms, sz, tot )
let mss = map ( \ (ms, _,_) -> ms ) multiSubsets
szs = map ( \ (_,sz,_) -> sz) multiSubsets
tot = sum $ map ( \(_,_,t) -> t) multiSubsets
Falsify.assert
$ Falsify.Prop.elem
.$ ("all multisubsets", map coerceVec2 $ multiSubsetsSum range tot $ coerceVec1 szs )
.$ ("random multisubset", mss)
testRoundTrip
:: ( Show a, Eq a )
=> Falsify.Gen a
-> ( a -> a )
-> Falsify.Property ()
testRoundTrip g roundTrip = do
d <- Falsify.gen g
Falsify.assert
$ Falsify.Prop.eq
.$ ("value", d )
.$ ("round tripped", roundTrip d )
testMonomialBasisD33 :: Falsify.Property ()
testMonomialBasisD33 =
testRoundTrip genD33 \ d -> $$( monTabulate \ mon -> monIndex [|| d ||] mon )
where
genD33 :: Falsify.Gen ( D3𝔸3 Double )
genD33 =
D33 <$> (unT <$> g)
<*> g <*> g <*> g
<*> g <*> g <*> g <*> g <*> g <*> g
<*> g <*> g <*> g <*> g <*> g <*> g <*> g <*> g <*> g <*> g
g :: Falsify.Gen ( T Double )
g = T . fromIntegral <$> Falsify.Gen.inRange ( Falsify.withOrigin ( -100, 100 ) ( 0 :: Int ) )
-- | Test the Faà di Bruno formula on polynomials, with a composition
-- \( g(f_1(x), f_2(x), .., f_n(x)) \).
testChainRule1NQ_1 :: Falsify.Property ()
testChainRule1NQ_1 = do
f <- genSpray "f" 1
g <- genSpray "g" 1
let gof_spray = Spray.composeSpray g [f]
gof_chain =
chain @_ @3 @( ℝ 1 ) ( ℝ1 <$> fromSpray @3 @( ℝ 1 ) f ) ( fromSpray @3 @( ℝ 1 ) g )
Falsify.assert
$ Falsify.Prop.eq
.$ ("direct", fromSpray @3 @( ℝ 1 ) gof_spray )
.$ ("chain rule", gof_chain )
-- | Test the Faà di Bruno formula on polynomials, with a composition
-- \( g(f_1(x), f_2(x), .., f_n(x)) \).
testChainRule1NQ_2 :: Falsify.Property ()
testChainRule1NQ_2 = do
f1 <- genSpray "f1" 1
f2 <- genSpray "f2" 1
g <- genSpray "g" 2
let gof_spray = Spray.composeSpray g [f1, f2]
f = ℝ2 <$> fromSpray @3 @( ℝ 1 ) f1
<*> fromSpray @3 @( ℝ 1 ) f2
gof_chain =
chain @_ @3 @( ℝ 2 ) f ( fromSpray @3 @( ℝ 2 ) g )
Falsify.assert
$ Falsify.Prop.eq
.$ ("direct", fromSpray @3 @( ℝ 1 ) gof_spray )
.$ ("chain rule", gof_chain )
-- | Test the Faà di Bruno formula on polynomials, with a composition
-- \( g(f_1(x), f_2(x), .., f_n(x)) \).
testChainRule1NQ_3 :: Falsify.Property ()
testChainRule1NQ_3 = do
f1 <- genSpray "f1" 1
f2 <- genSpray "f2" 1
f3 <- genSpray "f3" 1
g <- genSpray "g" 3
let gof_spray = Spray.composeSpray g [f1, f2, f3]
f = ℝ3 <$> fromSpray @3 @( ℝ 1 ) f1
<*> fromSpray @3 @( ℝ 1 ) f2
<*> fromSpray @3 @( ℝ 1 ) f3
gof_chain =
chain @_ @3 @( ℝ 3 ) f ( fromSpray @3 @( ℝ 3 ) g )
Falsify.assert
$ Falsify.Prop.eq
.$ ("direct", fromSpray @3 @( ℝ 1 ) gof_spray )
.$ ("chain rule", gof_chain )
class FromSpray v where
varFn :: Int -> v
linFn :: v -> Int -> Double
instance FromSpray ( ℝ 1 ) where
varFn = \case
0 -> ℝ1 1
i -> error $ "fromSpray in 1d but variable " ++ show i
linFn ( ℝ1 x ) = \case
0 -> x
i -> error $ "fromSpray in 1d but variable " ++ show i
instance FromSpray ( ℝ 2 ) where
varFn = \case
0 -> ℝ2 1 0
1 -> ℝ2 0 1
i -> error $ "fromSpray in 2d but variable " ++ show i
linFn ( ℝ2 x y ) = \case
0 -> x
1 -> y
i -> error $ "fromSpray in 2d but variable " ++ show i
instance FromSpray ( ℝ 3 ) where
varFn = \case
0 -> ℝ3 1 0 0
1 -> ℝ3 0 1 0
2 -> ℝ3 0 0 1
i -> error $ "fromSpray in 3d but variable " ++ show i
linFn ( ℝ3 x y z ) = \case
0 -> x
1 -> y
2 -> z
i -> error $ "fromSpray in 3d but variable " ++ show i
genSpray :: String -> Word -> Falsify.Property ( Spray Double )
genSpray lbl nbVars = Falsify.genWith (\ p -> Just $ lbl ++ " = " ++ Spray.prettySpray show "x" p) $ do
deg <- Falsify.Gen.inRange $ Falsify.between ( 0, 10 )
let mons = allMonomials deg nbVars
coeffs <-
for mons $ \ mon -> do
if all (== 0) mon
then return Nothing
else do
nonZero <- Falsify.Gen.bool False
if nonZero
then return Nothing
else do
-- Just use (small) integral values in tests for now,
-- to avoid errors arising from rounding.
c <- Falsify.Gen.inRange $ Falsify.withOrigin ( -100, 100 ) ( 0 :: Int )
return $ Just ( map fromIntegral mon, fromIntegral c )
return $ Spray.fromList $ catMaybes coeffs
allMonomials :: Word -> Word -> [ [ Word ] ]
allMonomials k _ | k < 0 = []
allMonomials _ 0 = [ [] ]
allMonomials 0 n = [ replicate ( fromIntegral n ) 0 ]
allMonomials k n = [ i : is | i <- reverse [ 0 .. k ], is <- allMonomials ( k - i ) ( n - 1 ) ]
-- | Convert a multivariate polynomial from the @hspray@ library to the dual algebra.
fromSpray
:: forall k v
. ( HasChainRule Double k v
, Module Double (T v)
, Applicative ( D k v )
, Ring ( D k v Double )
, FromSpray v
)
=> Spray Double
-> D k v Double
fromSpray coeffs = HashMap.foldlWithKey' addMonomial ( konst @Double @k @v $ HashMap.lookupDefault 0 (Spray.Powers mempty 0) coeffs ) coeffs
where
addMonomial :: D k v Double -> Spray.Powers -> Double -> D k v Double
addMonomial a xs c = a + monomial c ( toList $ Spray.exponents xs )
monomial :: Double -> [ Int ] -> D k v Double
monomial _ [] = konst @Double @k @v 0
monomial c is = fmap ( c * ) $ go 0 is
go :: Int -> [ Int ] -> D k v Double
go _ [] = konst @Double @k @v 1
go d (i : is) = pow d i * go ( d + 1 ) is
pow :: Int -> Int -> D k v Double
pow _ 0 = konst @Double @k @v 1
pow d i = linearD @Double @k @v ( \ x -> linFn @v x d ) ( unT origin :: v ) ^ ( fromIntegral i )