Further modularisation of root isolation code

The code is now generic over the dimension. There is a slight performance
loss that I need to investigate; perhaps some things are not getting
specialised? Maybe it is better to be more explicit about staging and
splice in the functions with fixed dimensions.
This commit is contained in:
sheaf 2024-04-18 20:14:19 +02:00
parent 131753da82
commit 0160081e80
10 changed files with 954 additions and 759 deletions

View file

@ -199,8 +199,8 @@ data TestCase =
TestCase
{ testName :: String
, testBrushStroke :: BrushStroke
, testCuspOptions :: RootIsolationOptions
, testStartBoxes :: [ ( Int, Box ) ]
, testCuspOptions :: RootIsolationOptions 2 3
, testStartBoxes :: [ ( Int, Box 2 ) ]
}
brushStrokeFunctions
@ -224,7 +224,7 @@ mkVal t i s = ( 1 t, i, 1 s )
mkI :: ( Double, Double ) -> 𝕀 Double
mkI ( lo, hi ) = 𝕀 lo hi
mkBox :: ( Double, Double ) -> ( Double, Double ) -> Box
mkBox :: ( Double, Double ) -> ( Double, Double ) -> Box 2
mkBox ( t_min, t_max ) ( s_min, s_max ) =
( 𝕀 ( 2 t_min s_min ) ( 2 t_max s_max ) )
@ -488,7 +488,7 @@ showD float = showFFloat (Just 6) float ""
--------------------------------------------------------------------------------
ellipseTestCase :: RootIsolationOptions -> String -> ( Double, Double ) -> Double -> [ ( Int, Box ) ] -> TestCase
ellipseTestCase :: RootIsolationOptions 2 3 -> String -> ( Double, Double ) -> Double -> [ ( Int, Box 2 ) ] -> TestCase
ellipseTestCase opts str k0k1 rot startBoxes =
TestCase
{ testName = "ellipse (" ++ str ++ ")"
@ -607,9 +607,11 @@ getStrokeFunctions brush sp0 crv =
brush @3 @𝕀 proxy# singleton )
{-# INLINEABLE getStrokeFunctions #-}
defaultStartBoxes :: [ Int ] -> [ ( Int, Box ) ]
defaultStartBoxes :: [ Int ] -> [ ( Int, Box 2 ) ]
defaultStartBoxes is =
[ ( i, mkBox (zero, one) (zero, one) ) | i <- is ]
[ ( i, 𝕀 ( 2 zero zero ) ( 2 one one ) )
| i <- is
]
getR1 (1 u) = u

View file

@ -1061,14 +1061,14 @@ instance HasChainRule Double 2 ( 1 ) where
konst :: forall w. AbelianGroup w => w -> D2𝔸1 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndex [|| df ||] zeroMonomial )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module Double ( T w ) => ( 1 -> w ) -> 1 -> D2𝔸1 w
linearD f v =
let !o = origin @Double @( T w )
in $$( monTabulate \ mon ->
in $$( monTabulateQ \ mon ->
if | isZeroMonomial mon
-> [|| f v ||]
| Just i <- isLinear mon
@ -1100,14 +1100,14 @@ instance HasChainRule Double 3 ( 1 ) where
konst :: forall w. AbelianGroup w => w -> D3𝔸1 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndex [|| df ||] zeroMonomial )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module Double ( T w ) => ( 1 -> w ) -> 1 -> D3𝔸1 w
linearD f v =
let !o = origin @Double @( T w )
in $$( monTabulate \ mon ->
in $$( monTabulateQ \ mon ->
if | isZeroMonomial mon
-> [|| f v ||]
| Just i <- isLinear mon
@ -1139,14 +1139,14 @@ instance HasChainRule Double 2 ( 2 ) where
konst :: forall w. AbelianGroup w => w -> D2𝔸2 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndex [|| df ||] zeroMonomial )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module Double ( T w ) => ( 2 -> w ) -> 2 -> D2𝔸2 w
linearD f v =
let !o = origin @Double @( T w )
in $$( monTabulate \ mon ->
in $$( monTabulateQ \ mon ->
if | isZeroMonomial mon
-> [|| f v ||]
| Just i <- isLinear mon
@ -1178,14 +1178,14 @@ instance HasChainRule Double 3 ( 2 ) where
konst :: forall w. AbelianGroup w => w -> D3𝔸2 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndex [|| df ||] zeroMonomial )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module Double ( T w ) => ( 2 -> w ) -> 2 -> D3𝔸2 w
linearD f v =
let !o = origin @Double @( T w )
in $$( monTabulate \ mon ->
in $$( monTabulateQ \ mon ->
if | isZeroMonomial mon
-> [|| f v ||]
| Just i <- isLinear mon
@ -1217,14 +1217,14 @@ instance HasChainRule Double 2 ( 3 ) where
konst :: forall w. AbelianGroup w => w -> D2𝔸3 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndex [|| df ||] zeroMonomial )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module Double ( T w ) => ( 3 -> w ) -> 3 -> D2𝔸3 w
linearD f v =
let !o = origin @Double @( T w )
in $$( monTabulate \ mon ->
in $$( monTabulateQ \ mon ->
if | isZeroMonomial mon
-> [|| f v ||]
| Just i <- isLinear mon
@ -1256,14 +1256,14 @@ instance HasChainRule Double 3 ( 3 ) where
konst :: forall w. AbelianGroup w => w -> D3𝔸3 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndex [|| df ||] zeroMonomial )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module Double ( T w ) => ( 3 -> w ) -> 3 -> D3𝔸3 w
linearD f v =
let !o = origin @Double @( T w )
in $$( monTabulate \ mon ->
in $$( monTabulateQ \ mon ->
if | isZeroMonomial mon
-> [|| f v ||]
| Just i <- isLinear mon
@ -1295,14 +1295,14 @@ instance HasChainRule Double 2 ( 4 ) where
konst :: forall w. AbelianGroup w => w -> D2𝔸4 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndex [|| df ||] zeroMonomial )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module Double ( T w ) => ( 4 -> w ) -> 4 -> D2𝔸4 w
linearD f v =
let !o = origin @Double @( T w )
in $$( monTabulate \ mon ->
in $$( monTabulateQ \ mon ->
if | isZeroMonomial mon
-> [|| f v ||]
| Just i <- isLinear mon
@ -1334,14 +1334,14 @@ instance HasChainRule Double 3 ( 4 ) where
konst :: forall w. AbelianGroup w => w -> D3𝔸4 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndex [|| df ||] zeroMonomial )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module Double ( T w ) => ( 4 -> w ) -> 4 -> D3𝔸4 w
linearD f v =
let !o = origin @Double @( T w )
in $$( monTabulate \ mon ->
in $$( monTabulateQ \ mon ->
if | isZeroMonomial mon
-> [|| f v ||]
| Just i <- isLinear mon

View file

@ -38,7 +38,8 @@ import Math.Linear
, zipIndices
)
import Math.Monomial
( Mon(..), MonomialBasis(..), Vars, Deg
( Mon(..), MonomialBasisQ(..), MonomialBasis(..)
, Vars, Deg
, mons, zeroMonomial, isZeroMonomial, totalDegree
, multiSubsetSumFaà, multiSubsetsSum
, partitionFaà, partitions
@ -203,8 +204,8 @@ data D3𝔸4 v =
-- (To be spliced in using Template Haskell.)
chainRule1NQ :: forall dr1 dv v r w d
. ( Ring r, RepresentableQ r v
, MonomialBasis dr1, Vars dr1 ~ 1
, MonomialBasis dv , Vars dv ~ RepDim v
, MonomialBasisQ dr1, Vars dr1 ~ 1
, MonomialBasisQ dv , Vars dv ~ RepDim v
, Deg dr1 ~ Deg dv
, d ~ Vars dv, KnownNat d
)
@ -215,19 +216,19 @@ chainRule1NQ :: forall dr1 dv v r w d
-> CodeQ ( dv w )
-> CodeQ ( dr1 w )
chainRule1NQ zero_w sum_w scale_w df dg =
monTabulate @dr1 \ ( Mon ( k `VS` _ ) ) ->
case k of
monTabulateQ @dr1 \ mon ->
case totalDegree mon of
-- Set the value of the composition separately,
-- as that isn't handled by the Faà di Bruno formula.
0 -> monIndex @dv dg zeroMonomial
_ ->
0 -> monIndexQ @dv dg zeroMonomial
k ->
[|| unT $ $$( foldQ sum_w zero_w
[ [|| $$scale_w ( T $$( monIndex @dv dg m_g ) )
[ [|| $$scale_w ( T $$( monIndexQ @dv dg m_g ) )
$$( foldQ [|| (Ring.+) ||] [|| Ring.fromInteger ( 0 :: Integer ) ||]
[ [|| Ring.fromInteger $$( liftTyped $ fromIntegral $ multiSubsetSumFaà k is ) Ring.*
$$( foldQ [|| (Ring.*) ||] [|| Ring.fromInteger ( 1 :: Integer ) ||]
[ foldQ [|| (Ring.*) ||] [|| Ring.fromInteger ( 1 :: Integer ) ||]
[ ( indexQ @r @v ( monIndex @dr1 df ( Mon $ f_deg `VS` VZ ) ) v_index )
[ ( indexQ @r @v ( monIndexQ @dr1 df ( Mon $ Vec [ f_deg ] ) ) v_index )
`powQ` pow
| ( f_deg, pow ) <- pows
]
@ -248,8 +249,8 @@ chainRule1NQ zero_w sum_w scale_w df dg =
--
-- (To be spliced in using Template Haskell.)
chainRuleN1Q :: forall du dr1 r
. ( MonomialBasis du
, MonomialBasis dr1, Vars dr1 ~ 1
. ( MonomialBasisQ du
, MonomialBasisQ dr1, Vars dr1 ~ 1
)
=> CodeQ ( Integer -> r ) -- ^ fromInteger (circumvent TH constraint issue)
-> CodeQ ( r -> r -> r ) -- ^ (+) (circumvent TH constraint issue)
@ -259,19 +260,19 @@ chainRuleN1Q :: forall du dr1 r
-> CodeQ ( dr1 r )
-> CodeQ ( du r )
chainRuleN1Q fromInt add times pow df dg =
monTabulate @du \ mon ->
monTabulateQ @du \ mon ->
if
| isZeroMonomial mon
-- Set the value of the composition separately,
-- as that isn't handled by the Faà di Bruno formula.
-> monIndex @dr1 dg zeroMonomial
-> monIndexQ @dr1 dg zeroMonomial
| otherwise
-> foldQ add [|| $$fromInt ( 0 :: Integer ) ||]
[ [|| $$times $$( monIndex @dr1 dg ( Mon ( k `VS` VZ ) ) )
[ [|| $$times $$( monIndexQ @dr1 dg ( Mon ( Vec [ k ] ) ) )
$$( foldQ add [|| $$fromInt ( 0 :: Integer ) ||]
[ [|| $$times ( $$fromInt $$( liftTyped $ fromIntegral $ partitionFaà mon is ) )
$$( foldQ times [|| $$fromInt ( 1 :: Integer ) ||]
[ [|| $$pow $$( monIndex @du df f_mon ) p ||]
[ [|| $$pow $$( monIndexQ @du df f_mon ) p ||]
| ( f_mon, p ) <- is
]
) ||]
@ -284,366 +285,395 @@ chainRuleN1Q fromInt add times pow df dg =
]
--------------------------------------------------------------------------------
-- MonomialBasis instances follow (nothing else).
-- MonomialBasisQ instances follow (nothing else).
type instance Deg D𝔸0 = 0
type instance Vars D𝔸0 = 0
instance MonomialBasis D𝔸0 where
monTabulate f =
instance MonomialBasisQ D𝔸0 where
monTabulateQ f =
[|| let
_D0_v = $$( f $ Mon VZ )
_D0_v = $$( f $ Mon ( Vec [ ] ) )
in D0 { .. }
||]
monIndex d _ = [|| _D0_v $$d ||]
monIndexQ d _ = [|| _D0_v $$d ||]
type instance Deg D1𝔸1 = 1
type instance Vars D1𝔸1 = 1
instance MonomialBasis D1𝔸1 where
monTabulate f =
instance MonomialBasisQ D1𝔸1 where
monTabulateQ f =
[|| let
_D11_v = $$( f $ Mon ( 0 `VS` VZ ) )
_D11_dx = T $$( f $ Mon ( 1 `VS` VZ ) )
_D11_v = $$( f $ Mon ( Vec [ 0 ] ) )
_D11_dx = T $$( f $ Mon ( Vec [ 1 ] ) )
in D11 { .. }
||]
monIndex d = \ case
Mon ( 1 `VS` VZ ) -> [|| unT $ _D11_dx $$d ||]
monIndexQ d = \ case
Mon ( Vec [ 1 ] ) -> [|| unT $ _D11_dx $$d ||]
_ -> [|| _D11_v $$d ||]
type instance Deg D2𝔸1 = 2
type instance Vars D2𝔸1 = 1
instance MonomialBasis D2𝔸1 where
monTabulate f =
instance MonomialBasisQ D2𝔸1 where
monTabulateQ f =
[|| let
_D21_v = $$( f $ Mon ( 0 `VS` VZ ) )
_D21_dx = T $$( f $ Mon ( 1 `VS` VZ ) )
_D21_dxdx = T $$( f $ Mon ( 2 `VS` VZ ) )
_D21_v = $$( f $ Mon ( Vec [ 0 ] ) )
_D21_dx = T $$( f $ Mon ( Vec [ 1 ] ) )
_D21_dxdx = T $$( f $ Mon ( Vec [ 2 ] ) )
in D21 { .. }
||]
monIndex d = \ case
Mon ( 1 `VS` VZ ) -> [|| unT $ _D21_dx $$d ||]
Mon ( 2 `VS` VZ ) -> [|| unT $ _D21_dxdx $$d ||]
monIndexQ d = \ case
Mon ( Vec [ 1 ] ) -> [|| unT $ _D21_dx $$d ||]
Mon ( Vec [ 2 ] ) -> [|| unT $ _D21_dxdx $$d ||]
_ -> [|| _D21_v $$d ||]
type instance Deg D3𝔸1 = 3
type instance Vars D3𝔸1 = 1
instance MonomialBasis D3𝔸1 where
monTabulate f =
instance MonomialBasisQ D3𝔸1 where
monTabulateQ f =
[|| let
_D31_v = $$( f $ Mon ( 0 `VS` VZ ) )
_D31_dx = T $$( f $ Mon ( 1 `VS` VZ ) )
_D31_dxdx = T $$( f $ Mon ( 2 `VS` VZ ) )
_D31_dxdxdx = T $$( f $ Mon ( 3 `VS` VZ ) )
_D31_v = $$( f $ Mon ( Vec [ 0 ] ) )
_D31_dx = T $$( f $ Mon ( Vec [ 1 ] ) )
_D31_dxdx = T $$( f $ Mon ( Vec [ 2 ] ) )
_D31_dxdxdx = T $$( f $ Mon ( Vec [ 3 ] ) )
in D31 { .. }
||]
monIndex d = \ case
Mon ( 1 `VS` VZ ) -> [|| unT $ _D31_dx $$d ||]
Mon ( 2 `VS` VZ ) -> [|| unT $ _D31_dxdx $$d ||]
Mon ( 3 `VS` VZ ) -> [|| unT $ _D31_dxdxdx $$d ||]
monIndexQ d = \ case
Mon ( Vec [ 1 ] ) -> [|| unT $ _D31_dx $$d ||]
Mon ( Vec [ 2 ] ) -> [|| unT $ _D31_dxdx $$d ||]
Mon ( Vec [ 3 ] ) -> [|| unT $ _D31_dxdxdx $$d ||]
_ -> [|| _D31_v $$d ||]
type instance Deg D1𝔸2 = 1
type instance Vars D1𝔸2 = 2
instance MonomialBasis D1𝔸2 where
monTabulate f =
instance MonomialBasisQ D1𝔸2 where
monTabulateQ f =
[|| let
_D12_v = $$( f $ Mon ( 0 `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_v = $$( f $ Mon ( Vec [ 0, 0 ] ) )
_D12_dx = T $$( f $ Mon ( Vec [ 1, 0 ] ) )
_D12_dy = T $$( f $ Mon ( Vec [ 0, 1 ] ) )
in D12 { .. }
||]
monIndexQ d = \ case
Mon ( Vec [ 1, 0 ] ) -> [|| unT $ _D12_dx $$d ||]
Mon ( Vec [ 0, 1 ] ) -> [|| unT $ _D12_dy $$d ||]
_ -> [|| _D12_v $$d ||]
instance MonomialBasis D1𝔸2 where
monTabulate f =
let
_D12_v = f $ Mon ( Vec [ 0, 0 ] )
_D12_dx = T $ f $ Mon ( Vec [ 1, 0 ] )
_D12_dy = T $ f $ Mon ( Vec [ 0, 1 ] )
in D12 { .. }
monIndex d = \ case
Mon ( 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D12_dx $$d ||]
Mon ( 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D12_dy $$d ||]
_ -> [|| _D12_v $$d ||]
Mon ( Vec [ 1, 0 ] ) -> unT $ _D12_dx d
Mon ( Vec [ 0, 1 ] ) -> unT $ _D12_dy d
_ -> _D12_v d
type instance Deg D2𝔸2 = 2
type instance Vars D2𝔸2 = 2
instance MonomialBasis D2𝔸2 where
monTabulate f =
instance MonomialBasisQ D2𝔸2 where
monTabulateQ f =
[|| let
_D22_v = $$( f $ Mon ( 0 `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_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` VZ ) )
_D22_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` VZ ) )
_D22_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` VZ ) )
_D22_v = $$( f $ Mon ( Vec [ 0, 0 ] ) )
_D22_dx = T $$( f $ Mon ( Vec [ 1, 0 ] ) )
_D22_dy = T $$( f $ Mon ( Vec [ 0, 1 ] ) )
_D22_dxdx = T $$( f $ Mon ( Vec [ 2, 0 ] ) )
_D22_dxdy = T $$( f $ Mon ( Vec [ 1, 1 ] ) )
_D22_dydy = T $$( f $ Mon ( Vec [ 0, 2 ] ) )
in D22 { .. }
||]
monIndex d = \ case
Mon ( 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D22_dx $$d ||]
Mon ( 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D22_dy $$d ||]
Mon ( 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D22_dxdx $$d ||]
Mon ( 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D22_dxdy $$d ||]
Mon ( 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D22_dydy $$d ||]
_ -> [|| _D22_v $$d ||]
monIndexQ d = \ case
Mon ( Vec [ 1, 0 ] ) -> [|| unT $ _D22_dx $$d ||]
Mon ( Vec [ 0, 1 ] ) -> [|| unT $ _D22_dy $$d ||]
Mon ( Vec [ 2, 0 ] ) -> [|| unT $ _D22_dxdx $$d ||]
Mon ( Vec [ 1, 1 ] ) -> [|| unT $ _D22_dxdy $$d ||]
Mon ( Vec [ 0, 2 ] ) -> [|| unT $ _D22_dydy $$d ||]
_ -> [|| _D22_v $$d ||]
type instance Deg D3𝔸2 = 3
type instance Vars D3𝔸2 = 2
instance MonomialBasis D3𝔸2 where
monTabulate f =
instance MonomialBasisQ D3𝔸2 where
monTabulateQ f =
[|| let
_D32_v = $$( f $ Mon ( 0 `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_dxdx = T $$( f $ Mon ( 2 `VS` 0 `VS` VZ ) )
_D32_dxdy = T $$( f $ Mon ( 1 `VS` 1 `VS` VZ ) )
_D32_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` VZ ) )
_D32_dxdxdx = T $$( f $ Mon ( 3 `VS` 0 `VS` VZ ) )
_D32_dxdxdy = T $$( f $ Mon ( 2 `VS` 1 `VS` VZ ) )
_D32_dxdydy = T $$( f $ Mon ( 1 `VS` 2 `VS` VZ ) )
_D32_dydydy = T $$( f $ Mon ( 0 `VS` 3 `VS` VZ ) )
_D32_v = $$( f $ Mon ( Vec [ 0, 0 ] ) )
_D32_dx = T $$( f $ Mon ( Vec [ 1, 0 ] ) )
_D32_dy = T $$( f $ Mon ( Vec [ 0, 1 ] ) )
_D32_dxdx = T $$( f $ Mon ( Vec [ 2, 0 ] ) )
_D32_dxdy = T $$( f $ Mon ( Vec [ 1, 1 ] ) )
_D32_dydy = T $$( f $ Mon ( Vec [ 0, 2 ] ) )
_D32_dxdxdx = T $$( f $ Mon ( Vec [ 3, 0 ] ) )
_D32_dxdxdy = T $$( f $ Mon ( Vec [ 2, 1 ] ) )
_D32_dxdydy = T $$( f $ Mon ( Vec [ 1, 2 ] ) )
_D32_dydydy = T $$( f $ Mon ( Vec [ 0, 3 ] ) )
in D32 { .. } ||]
monIndex d = \ case
Mon ( 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D32_dx $$d ||]
Mon ( 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D32_dy $$d ||]
Mon ( 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D32_dxdx $$d ||]
Mon ( 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D32_dxdy $$d ||]
Mon ( 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D32_dydy $$d ||]
Mon ( 3 `VS` 0 `VS` VZ ) -> [|| unT $ _D32_dxdxdx $$d ||]
Mon ( 2 `VS` 1 `VS` VZ ) -> [|| unT $ _D32_dxdxdy $$d ||]
Mon ( 1 `VS` 2 `VS` VZ ) -> [|| unT $ _D32_dxdydy $$d ||]
Mon ( 0 `VS` 3 `VS` VZ ) -> [|| unT $ _D32_dydydy $$d ||]
_ -> [|| _D32_v $$d ||]
monIndexQ d = \ case
Mon ( Vec [ 1, 0 ] ) -> [|| unT $ _D32_dx $$d ||]
Mon ( Vec [ 0, 1 ] ) -> [|| unT $ _D32_dy $$d ||]
Mon ( Vec [ 2, 0 ] ) -> [|| unT $ _D32_dxdx $$d ||]
Mon ( Vec [ 1, 1 ] ) -> [|| unT $ _D32_dxdy $$d ||]
Mon ( Vec [ 0, 2 ] ) -> [|| unT $ _D32_dydy $$d ||]
Mon ( Vec [ 3, 0 ] ) -> [|| unT $ _D32_dxdxdx $$d ||]
Mon ( Vec [ 2, 1 ] ) -> [|| unT $ _D32_dxdxdy $$d ||]
Mon ( Vec [ 1, 2 ] ) -> [|| unT $ _D32_dxdydy $$d ||]
Mon ( Vec [ 0, 3 ] ) -> [|| unT $ _D32_dydydy $$d ||]
_ -> [|| _D32_v $$d ||]
type instance Deg D1𝔸3 = 1
type instance Vars D1𝔸3 = 3
instance MonomialBasis D1𝔸3 where
monTabulate f =
instance MonomialBasisQ D1𝔸3 where
monTabulateQ f =
[|| let
!_D13_v = $$( f ( Mon ( 0 `VS` 0 `VS` 0 `VS` VZ ) ) )
!_D13_dx = T $$( f ( Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) ) )
!_D13_dy = T $$( f ( Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) ) )
!_D13_dz = T $$( f ( Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) ) )
!_D13_v = $$( f ( Mon ( Vec [ 0, 0, 0 ] ) ) )
!_D13_dx = T $$( f ( Mon ( Vec [ 1, 0, 0 ] ) ) )
!_D13_dy = T $$( f ( Mon ( Vec [ 0, 1, 0 ] ) ) )
!_D13_dz = T $$( f ( Mon ( Vec [ 0, 0, 1 ] ) ) )
in D13 { .. }
||]
monIndexQ d = \ case
Mon ( Vec [ 1, 0, 0 ] ) -> [|| unT $ _D13_dx $$d ||]
Mon ( Vec [ 0, 1, 0 ] ) -> [|| unT $ _D13_dy $$d ||]
Mon ( Vec [ 0, 0, 1 ] ) -> [|| unT $ _D13_dz $$d ||]
_ -> [|| _D13_v $$d ||]
instance MonomialBasis D1𝔸3 where
monTabulate f =
let
!_D13_v = f ( Mon ( Vec [ 0, 0, 0 ] ) )
!_D13_dx = T $ f ( Mon ( Vec [ 1, 0, 0 ] ) )
!_D13_dy = T $ f ( Mon ( Vec [ 0, 1, 0 ] ) )
!_D13_dz = T $ f ( Mon ( Vec [ 0, 0, 1 ] ) )
in D13 { .. }
monIndex d = \ case
Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D13_dx $$d ||]
Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D13_dy $$d ||]
Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D13_dz $$d ||]
_ -> [|| _D13_v $$d ||]
Mon ( Vec [ 1, 0, 0 ] ) -> unT $ _D13_dx d
Mon ( Vec [ 0, 1, 0 ] ) -> unT $ _D13_dy d
Mon ( Vec [ 0, 0, 1 ] ) -> unT $ _D13_dz d
_ -> _D13_v d
type instance Deg D2𝔸3 = 2
type instance Vars D2𝔸3 = 3
instance MonomialBasis D2𝔸3 where
monTabulate f =
instance MonomialBasisQ D2𝔸3 where
monTabulateQ f =
[|| let
_D23_v = $$( f $ Mon ( 0 `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_dz = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `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_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `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_dzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) )
_D23_v = $$( f $ Mon ( Vec [ 0, 0, 0 ] ) )
_D23_dx = T $$( f $ Mon ( Vec [ 1, 0, 0 ] ) )
_D23_dy = T $$( f $ Mon ( Vec [ 0, 1, 0 ] ) )
_D23_dz = T $$( f $ Mon ( Vec [ 0, 0, 1 ] ) )
_D23_dxdx = T $$( f $ Mon ( Vec [ 2, 0, 0 ] ) )
_D23_dxdy = T $$( f $ Mon ( Vec [ 1, 1, 0 ] ) )
_D23_dydy = T $$( f $ Mon ( Vec [ 0, 2, 0 ] ) )
_D23_dxdz = T $$( f $ Mon ( Vec [ 1, 0, 1 ] ) )
_D23_dydz = T $$( f $ Mon ( Vec [ 0, 1, 1 ] ) )
_D23_dzdz = T $$( f $ Mon ( Vec [ 0, 0, 2 ] ) )
in D23 { .. }
||]
monIndex d = \ case
Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D23_dx $$d ||]
Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D23_dy $$d ||]
Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D23_dz $$d ||]
Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D23_dxdx $$d ||]
Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D23_dxdy $$d ||]
Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D23_dydy $$d ||]
Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D23_dxdz $$d ||]
Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D23_dydz $$d ||]
Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D23_dzdz $$d ||]
_ -> [|| _D23_v $$d ||]
monIndexQ d = \ case
Mon ( Vec [ 1, 0, 0 ] ) -> [|| unT $ _D23_dx $$d ||]
Mon ( Vec [ 0, 1, 0 ] ) -> [|| unT $ _D23_dy $$d ||]
Mon ( Vec [ 0, 0, 1 ] ) -> [|| unT $ _D23_dz $$d ||]
Mon ( Vec [ 2, 0, 0 ] ) -> [|| unT $ _D23_dxdx $$d ||]
Mon ( Vec [ 1, 1, 0 ] ) -> [|| unT $ _D23_dxdy $$d ||]
Mon ( Vec [ 0, 2, 0 ] ) -> [|| unT $ _D23_dydy $$d ||]
Mon ( Vec [ 1, 0, 1 ] ) -> [|| unT $ _D23_dxdz $$d ||]
Mon ( Vec [ 0, 1, 1 ] ) -> [|| unT $ _D23_dydz $$d ||]
Mon ( Vec [ 0, 0, 2 ] ) -> [|| unT $ _D23_dzdz $$d ||]
_ -> [|| _D23_v $$d ||]
type instance Deg D3𝔸3 = 3
type instance Vars D3𝔸3 = 3
instance MonomialBasis D3𝔸3 where
monTabulate f =
instance MonomialBasisQ D3𝔸3 where
monTabulateQ f =
[|| let
_D33_v = $$( f $ Mon ( 0 `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_dz = T $$( f $ Mon ( 0 `VS` 0 `VS` 1 `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_dydy = T $$( f $ Mon ( 0 `VS` 2 `VS` 0 `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_dzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 2 `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_dxdydy = T $$( f $ Mon ( 1 `VS` 2 `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_dxdydz = T $$( f $ Mon ( 1 `VS` 1 `VS` 1 `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_dydzdz = T $$( f $ Mon ( 0 `VS` 1 `VS` 2 `VS` VZ ) )
_D33_dzdzdz = T $$( f $ Mon ( 0 `VS` 0 `VS` 3 `VS` VZ ) )
_D33_v = $$( f $ Mon ( Vec [ 0, 0, 0 ] ) )
_D33_dx = T $$( f $ Mon ( Vec [ 1, 0, 0 ] ) )
_D33_dy = T $$( f $ Mon ( Vec [ 0, 1, 0 ] ) )
_D33_dz = T $$( f $ Mon ( Vec [ 0, 0, 1 ] ) )
_D33_dxdx = T $$( f $ Mon ( Vec [ 2, 0, 0 ] ) )
_D33_dxdy = T $$( f $ Mon ( Vec [ 1, 1, 0 ] ) )
_D33_dydy = T $$( f $ Mon ( Vec [ 0, 2, 0 ] ) )
_D33_dxdz = T $$( f $ Mon ( Vec [ 1, 0, 1 ] ) )
_D33_dydz = T $$( f $ Mon ( Vec [ 0, 1, 1 ] ) )
_D33_dzdz = T $$( f $ Mon ( Vec [ 0, 0, 2 ] ) )
_D33_dxdxdx = T $$( f $ Mon ( Vec [ 3, 0, 0 ] ) )
_D33_dxdxdy = T $$( f $ Mon ( Vec [ 2, 1, 0 ] ) )
_D33_dxdydy = T $$( f $ Mon ( Vec [ 1, 2, 0 ] ) )
_D33_dydydy = T $$( f $ Mon ( Vec [ 0, 3, 0 ] ) )
_D33_dxdxdz = T $$( f $ Mon ( Vec [ 2, 0, 1 ] ) )
_D33_dxdydz = T $$( f $ Mon ( Vec [ 1, 1, 1 ] ) )
_D33_dxdzdz = T $$( f $ Mon ( Vec [ 1, 0, 2 ] ) )
_D33_dydydz = T $$( f $ Mon ( Vec [ 0, 2, 1 ] ) )
_D33_dydzdz = T $$( f $ Mon ( Vec [ 0, 1, 2 ] ) )
_D33_dzdzdz = T $$( f $ Mon ( Vec [ 0, 0, 3 ] ) )
in D33 { .. } ||]
monIndex d = \ case
Mon ( 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dx $$d ||]
Mon ( 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dy $$d ||]
Mon ( 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dz $$d ||]
Mon ( 2 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dxdx $$d ||]
Mon ( 1 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dxdy $$d ||]
Mon ( 0 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dydy $$d ||]
Mon ( 1 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dxdz $$d ||]
Mon ( 0 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D33_dydz $$d ||]
Mon ( 0 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D33_dzdz $$d ||]
Mon ( 3 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dxdxdx $$d ||]
Mon ( 2 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dxdxdy $$d ||]
Mon ( 1 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dxdydy $$d ||]
Mon ( 0 `VS` 3 `VS` 0 `VS` VZ ) -> [|| unT $ _D33_dydydy $$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` 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` 0 `VS` 3 `VS` VZ ) -> [|| unT $ _D33_dzdzdz $$d ||]
_ -> [|| _D33_v $$d ||]
monIndexQ d = \ case
Mon ( Vec [ 1, 0, 0 ] ) -> [|| unT $ _D33_dx $$d ||]
Mon ( Vec [ 0, 1, 0 ] ) -> [|| unT $ _D33_dy $$d ||]
Mon ( Vec [ 0, 0, 1 ] ) -> [|| unT $ _D33_dz $$d ||]
Mon ( Vec [ 2, 0, 0 ] ) -> [|| unT $ _D33_dxdx $$d ||]
Mon ( Vec [ 1, 1, 0 ] ) -> [|| unT $ _D33_dxdy $$d ||]
Mon ( Vec [ 0, 2, 0 ] ) -> [|| unT $ _D33_dydy $$d ||]
Mon ( Vec [ 1, 0, 1 ] ) -> [|| unT $ _D33_dxdz $$d ||]
Mon ( Vec [ 0, 1, 1 ] ) -> [|| unT $ _D33_dydz $$d ||]
Mon ( Vec [ 0, 0, 2 ] ) -> [|| unT $ _D33_dzdz $$d ||]
Mon ( Vec [ 3, 0, 0 ] ) -> [|| unT $ _D33_dxdxdx $$d ||]
Mon ( Vec [ 2, 1, 0 ] ) -> [|| unT $ _D33_dxdxdy $$d ||]
Mon ( Vec [ 1, 2, 0 ] ) -> [|| unT $ _D33_dxdydy $$d ||]
Mon ( Vec [ 0, 3, 0 ] ) -> [|| unT $ _D33_dydydy $$d ||]
Mon ( Vec [ 2, 0, 1 ] ) -> [|| unT $ _D33_dxdxdz $$d ||]
Mon ( Vec [ 1, 1, 1 ] ) -> [|| unT $ _D33_dxdydz $$d ||]
Mon ( Vec [ 1, 0, 2 ] ) -> [|| unT $ _D33_dxdzdz $$d ||]
Mon ( Vec [ 0, 2, 1 ] ) -> [|| unT $ _D33_dydydz $$d ||]
Mon ( Vec [ 0, 1, 2 ] ) -> [|| unT $ _D33_dydzdz $$d ||]
Mon ( Vec [ 0, 0, 3 ] ) -> [|| unT $ _D33_dzdzdz $$d ||]
_ -> [|| _D33_v $$d ||]
type instance Deg D1𝔸4 = 1
type instance Vars D1𝔸4 = 4
instance MonomialBasis D1𝔸4 where
monTabulate f =
instance MonomialBasisQ D1𝔸4 where
monTabulateQ f =
[|| let
_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_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_dw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) )
_D14_v = $$( f $ Mon ( Vec [ 0, 0, 0, 0 ] ) )
_D14_dx = T $$( f $ Mon ( Vec [ 1, 0, 0, 0 ] ) )
_D14_dy = T $$( f $ Mon ( Vec [ 0, 1, 0, 0 ] ) )
_D14_dz = T $$( f $ Mon ( Vec [ 0, 0, 1, 0 ] ) )
_D14_dw = T $$( f $ Mon ( Vec [ 0, 0, 0, 1 ] ) )
in D14 { .. } ||]
monIndex d = \ case
Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D14_dx $$d ||]
Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D14_dy $$d ||]
Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D14_dz $$d ||]
Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D14_dw $$d ||]
_ -> [|| _D14_v $$d ||]
monIndexQ d = \ case
Mon ( Vec [ 1, 0, 0, 0 ] ) -> [|| unT $ _D14_dx $$d ||]
Mon ( Vec [ 0, 1, 0, 0 ] ) -> [|| unT $ _D14_dy $$d ||]
Mon ( Vec [ 0, 0, 1, 0 ] ) -> [|| unT $ _D14_dz $$d ||]
Mon ( Vec [ 0, 0, 0, 1 ] ) -> [|| unT $ _D14_dw $$d ||]
_ -> [|| _D14_v $$d ||]
type instance Deg D2𝔸4 = 2
type instance Vars D2𝔸4 = 4
instance MonomialBasis D2𝔸4 where
monTabulate f =
instance MonomialBasisQ D2𝔸4 where
monTabulateQ f =
[|| let
_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_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_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_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_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_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_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_dwdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) )
_D24_v = $$( f $ Mon ( Vec [ 0, 0, 0, 0 ] ) )
_D24_dx = T $$( f $ Mon ( Vec [ 1, 0, 0, 0 ] ) )
_D24_dy = T $$( f $ Mon ( Vec [ 0, 1, 0, 0 ] ) )
_D24_dz = T $$( f $ Mon ( Vec [ 0, 0, 1, 0 ] ) )
_D24_dw = T $$( f $ Mon ( Vec [ 0, 0, 0, 1 ] ) )
_D24_dxdx = T $$( f $ Mon ( Vec [ 2, 0, 0, 0 ] ) )
_D24_dxdy = T $$( f $ Mon ( Vec [ 1, 1, 0, 0 ] ) )
_D24_dydy = T $$( f $ Mon ( Vec [ 0, 2, 0, 0 ] ) )
_D24_dxdz = T $$( f $ Mon ( Vec [ 1, 0, 1, 0 ] ) )
_D24_dydz = T $$( f $ Mon ( Vec [ 0, 1, 1, 0 ] ) )
_D24_dzdz = T $$( f $ Mon ( Vec [ 0, 0, 2, 0 ] ) )
_D24_dxdw = T $$( f $ Mon ( Vec [ 1, 0, 0, 1 ] ) )
_D24_dydw = T $$( f $ Mon ( Vec [ 0, 1, 0, 1 ] ) )
_D24_dzdw = T $$( f $ Mon ( Vec [ 0, 0, 1, 1 ] ) )
_D24_dwdw = T $$( f $ Mon ( Vec [ 0, 0, 0, 2 ] ) )
in D24 { .. } ||]
monIndex d = \ case
Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dx $$d ||]
Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dy $$d ||]
Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dz $$d ||]
Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D24_dw $$d ||]
Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dxdx $$d ||]
Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dxdy $$d ||]
Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dydy $$d ||]
Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dxdz $$d ||]
Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dydz $$d ||]
Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D24_dzdz $$d ||]
Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D24_dxdw $$d ||]
Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D24_dydw $$d ||]
Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D24_dzdw $$d ||]
Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D24_dwdw $$d ||]
_ -> [|| _D24_v $$d ||]
monIndexQ d = \ case
Mon ( Vec [ 1, 0, 0, 0 ] ) -> [|| unT $ _D24_dx $$d ||]
Mon ( Vec [ 0, 1, 0, 0 ] ) -> [|| unT $ _D24_dy $$d ||]
Mon ( Vec [ 0, 0, 1, 0 ] ) -> [|| unT $ _D24_dz $$d ||]
Mon ( Vec [ 0, 0, 0, 1 ] ) -> [|| unT $ _D24_dw $$d ||]
Mon ( Vec [ 2, 0, 0, 0 ] ) -> [|| unT $ _D24_dxdx $$d ||]
Mon ( Vec [ 1, 1, 0, 0 ] ) -> [|| unT $ _D24_dxdy $$d ||]
Mon ( Vec [ 0, 2, 0, 0 ] ) -> [|| unT $ _D24_dydy $$d ||]
Mon ( Vec [ 1, 0, 1, 0 ] ) -> [|| unT $ _D24_dxdz $$d ||]
Mon ( Vec [ 0, 1, 1, 0 ] ) -> [|| unT $ _D24_dydz $$d ||]
Mon ( Vec [ 0, 0, 2, 0 ] ) -> [|| unT $ _D24_dzdz $$d ||]
Mon ( Vec [ 1, 0, 0, 1 ] ) -> [|| unT $ _D24_dxdw $$d ||]
Mon ( Vec [ 0, 1, 0, 1 ] ) -> [|| unT $ _D24_dydw $$d ||]
Mon ( Vec [ 0, 0, 1, 1 ] ) -> [|| unT $ _D24_dzdw $$d ||]
Mon ( Vec [ 0, 0, 0, 2 ] ) -> [|| unT $ _D24_dwdw $$d ||]
_ -> [|| _D24_v $$d ||]
type instance Deg D3𝔸4 = 3
type instance Vars D3𝔸4 = 4
instance MonomialBasis D3𝔸4 where
monTabulate f =
instance MonomialBasisQ D3𝔸4 where
monTabulateQ f =
[|| let
_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_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_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_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_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_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_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_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_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_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_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_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_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_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_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_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_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_dwdwdw = T $$( f $ Mon ( 0 `VS` 0 `VS` 0 `VS` 3 `VS` VZ ) )
_D34_v = $$( f $ Mon ( Vec [ 0, 0, 0, 0 ] ) )
_D34_dx = T $$( f $ Mon ( Vec [ 1, 0, 0, 0 ] ) )
_D34_dy = T $$( f $ Mon ( Vec [ 0, 1, 0, 0 ] ) )
_D34_dz = T $$( f $ Mon ( Vec [ 0, 0, 1, 0 ] ) )
_D34_dw = T $$( f $ Mon ( Vec [ 0, 0, 0, 1 ] ) )
_D34_dxdx = T $$( f $ Mon ( Vec [ 2, 0, 0, 0 ] ) )
_D34_dxdy = T $$( f $ Mon ( Vec [ 1, 1, 0, 0 ] ) )
_D34_dydy = T $$( f $ Mon ( Vec [ 0, 2, 0, 0 ] ) )
_D34_dxdz = T $$( f $ Mon ( Vec [ 1, 0, 1, 0 ] ) )
_D34_dydz = T $$( f $ Mon ( Vec [ 0, 1, 1, 0 ] ) )
_D34_dzdz = T $$( f $ Mon ( Vec [ 0, 0, 2, 0 ] ) )
_D34_dxdw = T $$( f $ Mon ( Vec [ 1, 0, 0, 1 ] ) )
_D34_dydw = T $$( f $ Mon ( Vec [ 0, 1, 0, 1 ] ) )
_D34_dzdw = T $$( f $ Mon ( Vec [ 0, 0, 1, 1 ] ) )
_D34_dwdw = T $$( f $ Mon ( Vec [ 0, 0, 0, 2 ] ) )
_D34_dxdxdx = T $$( f $ Mon ( Vec [ 3, 0, 0, 0 ] ) )
_D34_dxdxdy = T $$( f $ Mon ( Vec [ 2, 1, 0, 0 ] ) )
_D34_dxdydy = T $$( f $ Mon ( Vec [ 1, 2, 0, 0 ] ) )
_D34_dydydy = T $$( f $ Mon ( Vec [ 0, 3, 0, 0 ] ) )
_D34_dxdxdz = T $$( f $ Mon ( Vec [ 2, 0, 1, 0 ] ) )
_D34_dxdydz = T $$( f $ Mon ( Vec [ 1, 1, 1, 0 ] ) )
_D34_dxdzdz = T $$( f $ Mon ( Vec [ 1, 0, 2, 0 ] ) )
_D34_dydydz = T $$( f $ Mon ( Vec [ 0, 2, 1, 0 ] ) )
_D34_dydzdz = T $$( f $ Mon ( Vec [ 0, 1, 2, 0 ] ) )
_D34_dzdzdz = T $$( f $ Mon ( Vec [ 0, 0, 3, 0 ] ) )
_D34_dxdxdw = T $$( f $ Mon ( Vec [ 2, 0, 0, 1 ] ) )
_D34_dxdydw = T $$( f $ Mon ( Vec [ 1, 1, 0, 1 ] ) )
_D34_dydydw = T $$( f $ Mon ( Vec [ 0, 2, 0, 1 ] ) )
_D34_dxdzdw = T $$( f $ Mon ( Vec [ 1, 0, 1, 1 ] ) )
_D34_dydzdw = T $$( f $ Mon ( Vec [ 0, 1, 1, 1 ] ) )
_D34_dzdzdw = T $$( f $ Mon ( Vec [ 0, 0, 2, 1 ] ) )
_D34_dxdwdw = T $$( f $ Mon ( Vec [ 1, 0, 0, 2 ] ) )
_D34_dydwdw = T $$( f $ Mon ( Vec [ 0, 1, 0, 2 ] ) )
_D34_dzdwdw = T $$( f $ Mon ( Vec [ 0, 0, 1, 2 ] ) )
_D34_dwdwdw = T $$( f $ Mon ( Vec [ 0, 0, 0, 3 ] ) )
in D34 { .. } ||]
monIndex d = \ case
Mon ( 1 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dx $$d ||]
Mon ( 0 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dy $$d ||]
Mon ( 0 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dz $$d ||]
Mon ( 0 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dw $$d ||]
Mon ( 2 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdx $$d ||]
Mon ( 1 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdy $$d ||]
Mon ( 0 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dydy $$d ||]
Mon ( 1 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdz $$d ||]
Mon ( 0 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dydz $$d ||]
Mon ( 0 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dzdz $$d ||]
Mon ( 1 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dxdw $$d ||]
Mon ( 0 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dydw $$d ||]
Mon ( 0 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dzdw $$d ||]
Mon ( 0 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D34_dwdw $$d ||]
Mon ( 3 `VS` 0 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdxdx $$d ||]
Mon ( 2 `VS` 1 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdxdy $$d ||]
Mon ( 1 `VS` 2 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdydy $$d ||]
Mon ( 0 `VS` 3 `VS` 0 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dydydy $$d ||]
Mon ( 2 `VS` 0 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdxdz $$d ||]
Mon ( 1 `VS` 1 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdydz $$d ||]
Mon ( 1 `VS` 0 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dxdzdz $$d ||]
Mon ( 0 `VS` 2 `VS` 1 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dydydz $$d ||]
Mon ( 0 `VS` 1 `VS` 2 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dydzdz $$d ||]
Mon ( 0 `VS` 0 `VS` 3 `VS` 0 `VS` VZ ) -> [|| unT $ _D34_dzdzdz $$d ||]
Mon ( 2 `VS` 0 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dxdxdw $$d ||]
Mon ( 1 `VS` 1 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dxdydw $$d ||]
Mon ( 0 `VS` 2 `VS` 0 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dydydw $$d ||]
Mon ( 1 `VS` 0 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dxdzdw $$d ||]
Mon ( 0 `VS` 1 `VS` 1 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dydzdw $$d ||]
Mon ( 0 `VS` 0 `VS` 2 `VS` 1 `VS` VZ ) -> [|| unT $ _D34_dzdzdw $$d ||]
Mon ( 1 `VS` 0 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D34_dxdwdw $$d ||]
Mon ( 0 `VS` 1 `VS` 0 `VS` 2 `VS` VZ ) -> [|| unT $ _D34_dydwdw $$d ||]
Mon ( 0 `VS` 0 `VS` 1 `VS` 2 `VS` VZ ) -> [|| unT $ _D34_dzdwdw $$d ||]
Mon ( 0 `VS` 0 `VS` 0 `VS` 3 `VS` VZ ) -> [|| unT $ _D34_dwdwdw $$d ||]
_ -> [|| _D34_v $$d ||]
monIndexQ d = \ case
Mon ( Vec [ 1, 0, 0, 0 ] ) -> [|| unT $ _D34_dx $$d ||]
Mon ( Vec [ 0, 1, 0, 0 ] ) -> [|| unT $ _D34_dy $$d ||]
Mon ( Vec [ 0, 0, 1, 0 ] ) -> [|| unT $ _D34_dz $$d ||]
Mon ( Vec [ 0, 0, 0, 1 ] ) -> [|| unT $ _D34_dw $$d ||]
Mon ( Vec [ 2, 0, 0, 0 ] ) -> [|| unT $ _D34_dxdx $$d ||]
Mon ( Vec [ 1, 1, 0, 0 ] ) -> [|| unT $ _D34_dxdy $$d ||]
Mon ( Vec [ 0, 2, 0, 0 ] ) -> [|| unT $ _D34_dydy $$d ||]
Mon ( Vec [ 1, 0, 1, 0 ] ) -> [|| unT $ _D34_dxdz $$d ||]
Mon ( Vec [ 0, 1, 1, 0 ] ) -> [|| unT $ _D34_dydz $$d ||]
Mon ( Vec [ 0, 0, 2, 0 ] ) -> [|| unT $ _D34_dzdz $$d ||]
Mon ( Vec [ 1, 0, 0, 1 ] ) -> [|| unT $ _D34_dxdw $$d ||]
Mon ( Vec [ 0, 1, 0, 1 ] ) -> [|| unT $ _D34_dydw $$d ||]
Mon ( Vec [ 0, 0, 1, 1 ] ) -> [|| unT $ _D34_dzdw $$d ||]
Mon ( Vec [ 0, 0, 0, 2 ] ) -> [|| unT $ _D34_dwdw $$d ||]
Mon ( Vec [ 3, 0, 0, 0 ] ) -> [|| unT $ _D34_dxdxdx $$d ||]
Mon ( Vec [ 2, 1, 0, 0 ] ) -> [|| unT $ _D34_dxdxdy $$d ||]
Mon ( Vec [ 1, 2, 0, 0 ] ) -> [|| unT $ _D34_dxdydy $$d ||]
Mon ( Vec [ 0, 3, 0, 0 ] ) -> [|| unT $ _D34_dydydy $$d ||]
Mon ( Vec [ 2, 0, 1, 0 ] ) -> [|| unT $ _D34_dxdxdz $$d ||]
Mon ( Vec [ 1, 1, 1, 0 ] ) -> [|| unT $ _D34_dxdydz $$d ||]
Mon ( Vec [ 1, 0, 2, 0 ] ) -> [|| unT $ _D34_dxdzdz $$d ||]
Mon ( Vec [ 0, 2, 1, 0 ] ) -> [|| unT $ _D34_dydydz $$d ||]
Mon ( Vec [ 0, 1, 2, 0 ] ) -> [|| unT $ _D34_dydzdz $$d ||]
Mon ( Vec [ 0, 0, 3, 0 ] ) -> [|| unT $ _D34_dzdzdz $$d ||]
Mon ( Vec [ 2, 0, 0, 1 ] ) -> [|| unT $ _D34_dxdxdw $$d ||]
Mon ( Vec [ 1, 1, 0, 1 ] ) -> [|| unT $ _D34_dxdydw $$d ||]
Mon ( Vec [ 0, 2, 0, 1 ] ) -> [|| unT $ _D34_dydydw $$d ||]
Mon ( Vec [ 1, 0, 1, 1 ] ) -> [|| unT $ _D34_dxdzdw $$d ||]
Mon ( Vec [ 0, 1, 1, 1 ] ) -> [|| unT $ _D34_dydzdw $$d ||]
Mon ( Vec [ 0, 0, 2, 1 ] ) -> [|| unT $ _D34_dzdzdw $$d ||]
Mon ( Vec [ 1, 0, 0, 2 ] ) -> [|| unT $ _D34_dxdwdw $$d ||]
Mon ( Vec [ 0, 1, 0, 2 ] ) -> [|| unT $ _D34_dydwdw $$d ||]
Mon ( Vec [ 0, 0, 1, 2 ] ) -> [|| unT $ _D34_dzdwdw $$d ||]
Mon ( Vec [ 0, 0, 0, 3 ] ) -> [|| unT $ _D34_dwdwdw $$d ||]
_ -> [|| _D34_v $$d ||]

View file

@ -243,7 +243,7 @@ computeStrokeOutline ::
)
=> RootSolvingAlgorithm
-> Maybe RootIsolationOptions
-> Maybe ( RootIsolationOptions 2 3 )
-> FitParameters
-> ( ptData -> usedParams )
-> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing
@ -522,7 +522,7 @@ outlineFunction
, Show ptData, Show brushParams
)
=> RootSolvingAlgorithm
-> Maybe RootIsolationOptions
-> Maybe ( RootIsolationOptions 2 3 )
-> ( ptData -> usedParams )
-> ( usedParams -> brushParams ) -- ^ assumed to be linear and non-decreasing
-> ( forall {t} k (i :: t)
@ -1147,7 +1147,7 @@ solveEnvelopeEquations rootAlgo _t path_t path'_t ( fwdOffset, bwdOffset ) strok
-- TODO: use Newton's method starting at the midpoint of the box,
-- instead of just taking the midpoint of the box.
cuspCoords :: ( 1 -> Seq ( 1 -> StrokeDatum 2 () ) )
-> ( Int, Box )
-> ( Int, Box 2 )
-> Cusp
cuspCoords eqs ( i, 𝕀 ( 2 t_lo s_lo ) ( 2 t_hi s_hi ) )
| StrokeDatum
@ -1168,13 +1168,13 @@ cuspCoords eqs ( i, 𝕀 ( 2 t_lo s_lo ) ( 2 t_hi s_hi ) )
--
-- See 'isolateRootsIn' for details of the algorithms used.
findCusps
:: RootIsolationOptions
:: RootIsolationOptions 2 3
-> ( 𝕀 1 -> Seq ( 𝕀 1 -> StrokeDatum 3 𝕀 ) )
-> IntMap ( [ ( Box, RootIsolationTree Box ) ], ( [ Box ], [ Box ] ) )
-> IntMap ( [ ( Box 2, RootIsolationTree ( Box 2 ) ) ], ( [ Box 2 ], [ Box 2 ] ) )
findCusps opts boxStrokeData =
findCuspsIn opts boxStrokeData $
IntMap.fromList
[ ( i, [ 𝕀 ( 2 zero zero ) ( 2 one one ) ] )
[ ( i, [ 𝕀 ( 2 zero zero ) ( 2 one one ) ] ) --[ 𝕀 ( 3 zero zero -1e6 ) ( 3 one one 1e6 ) ] )
| i <- [ 0 .. length ( boxStrokeData ( 𝕀 ( 1 zero ) ( 1 one ) ) ) - 1 ]
]
where
@ -1187,10 +1187,10 @@ findCusps opts boxStrokeData =
-- | Like 'findCusps', but parametrised over the initial boxes for the
-- root isolation method.
findCuspsIn
:: RootIsolationOptions
:: RootIsolationOptions 2 3
-> ( 𝕀 1 -> Seq ( 𝕀 1 -> StrokeDatum 3 𝕀 ) )
-> IntMap [ Box ]
-> IntMap ( [ ( Box, RootIsolationTree Box ) ], ( [ Box ], [ Box ] ) )
-> IntMap [ Box 2 ]
-> IntMap ( [ ( Box 2, RootIsolationTree ( Box 2 ) ) ], ( [ Box 2 ], [ Box 2 ] ) )
findCuspsIn opts boxStrokeData initBoxes =
IntMap.mapWithKey ( \ i -> foldMap ( isolateRootsIn opts ( eqnPiece i ) ) ) initBoxes
where
@ -1211,3 +1211,47 @@ findCuspsIn opts boxStrokeData initBoxes =
in D12 ( 𝕀 ( 3 ee_lo f_lo g_lo ) ( 3 ee_hi f_hi g_hi ) )
( T $ 𝕀 ( 3 ee_t_lo f_t_lo g_t_lo ) ( 3 ee_t_hi f_t_hi g_t_hi ) )
( T $ 𝕀 ( 3 ee_s_lo f_s_lo g_s_lo ) ( 3 ee_s_hi f_s_hi g_s_hi ) )
{-
findCuspsIn
:: RootIsolationOptions 3 3
-> ( 𝕀 1 -> Seq ( 𝕀 1 -> StrokeDatum 3 𝕀 ) )
-> IntMap [ Box 3 ]
-> IntMap ( [ ( Box 3, RootIsolationTree ( Box 3 ) ) ], ( [ Box 3 ], [ Box 3 ] ) )
findCuspsIn opts boxStrokeData initBoxes =
IntMap.mapWithKey ( \ i -> foldMap ( isolateRootsIn opts ( eqnPiece i ) ) ) initBoxes
where
eqnPiece :: Int -> 𝕀 ( 3 ) -> D1𝔸3 ( 𝕀 ( 3 ) )
eqnPiece i ( 𝕀 ( 3 t_lo s_lo λ_lo ) ( 3 t_hi s_hi λ_hi ) ) =
let t = 𝕀 ( 1 t_lo ) ( 1 t_hi )
s = 𝕀 ( 1 s_lo ) ( 1 s_hi )
λ = 𝕀 λ_lo λ_hi
StrokeDatum
{ dstroke =
D32 { _D32_dx = T c_t, _D32_dy = T c_s
, _D32_dxdx = T c_tt, _D32_dxdy = T c_ts, _D32_dydy = T c_ss }
, ee =
D22 { _D22_dx = T ee_t, _D22_dy = T ee_s
, _D22_dxdx = T ee_tt, _D22_dxdy = T ee_ts, _D22_dydy = T ee_ss }
} = ( boxStrokeData t `Seq.index` i ) s
-- λ ∂E/∂s - ∂E/∂t = 0
𝕀 ( 1 f1_lo ) ( 1 f1_hi ) = unT $ λ *^ T ee_s ^-^ T ee_t
𝕀 ( 1 f1_t_lo ) ( 1 f1_t_hi ) = unT $ λ *^ T ee_ts ^-^ T ee_tt
𝕀 ( 1 f1_s_lo ) ( 1 f1_s_hi ) = unT $ λ *^ T ee_ss ^-^ T ee_ts
𝕀 ( 1 f1_λ_lo ) ( 1 f1_λ_hi ) = ee_s
-- λ ∂c/∂s - ∂c/∂t = 0
𝕀 ( 2 f2_lo f3_lo ) ( 2 f2_hi f3_hi ) = unT $ λ *^ T c_s ^-^ T c_t
𝕀 ( 2 f2_t_lo f3_t_lo ) ( 2 f2_t_hi f3_t_hi ) = unT $ λ *^ T c_ts ^-^ T c_tt
𝕀 ( 2 f2_s_lo f3_s_lo ) ( 2 f2_s_hi f3_s_hi ) = unT $ λ *^ T c_ss ^-^ T c_ts
𝕀 ( 2 f2_λ_lo f3_λ_lo ) ( 2 f2_λ_hi f3_λ_hi ) = c_s
in D13 ( 𝕀 ( 3 f1_lo f2_lo f3_lo ) ( 3 f1_hi f2_hi f3_hi ) )
( T $ 𝕀 ( 3 f1_t_lo f2_t_lo f3_t_lo ) ( 3 f1_t_hi f2_t_hi f3_t_hi ) )
( T $ 𝕀 ( 3 f1_s_lo f2_s_lo f3_s_lo ) ( 3 f1_s_hi f2_s_hi f3_s_hi ) )
( T $ 𝕀 ( 3 f1_λ_lo f2_λ_lo f3_λ_lo ) ( 3 f1_λ_hi f2_λ_hi f3_λ_hi ) )
-}

View file

@ -179,14 +179,14 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀 1 ) where
konst :: forall w. AbelianGroup w => w -> D2𝔸1 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndex [|| df ||] zeroMonomial )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀 1 -> w ) -> 𝕀 1 -> D2𝔸1 w
linearD f v =
let !o = origin @( 𝕀 Double ) @( T w )
in $$( monTabulate \ mon ->
in $$( monTabulateQ \ mon ->
if | isZeroMonomial mon
-> [|| f v ||]
| Just i <- isLinear mon
@ -214,14 +214,14 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀 1 ) where
konst :: forall w. AbelianGroup w => w -> D3𝔸1 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndex [|| df ||] zeroMonomial )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀 1 -> w ) -> 𝕀 1 -> D3𝔸1 w
linearD f v =
let !o = origin @( 𝕀 Double ) @( T w )
in $$( monTabulate \ mon ->
in $$( monTabulateQ \ mon ->
if | isZeroMonomial mon
-> [|| f v ||]
| Just i <- isLinear mon
@ -249,14 +249,14 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀 2 ) where
konst :: forall w. AbelianGroup w => w -> D2𝔸2 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndex [|| df ||] zeroMonomial )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀 2 -> w ) -> 𝕀 2 -> D2𝔸2 w
linearD f v =
let !o = origin @( 𝕀 Double ) @( T w )
in $$( monTabulate \ mon ->
in $$( monTabulateQ \ mon ->
if | isZeroMonomial mon
-> [|| f v ||]
| Just i <- isLinear mon
@ -284,14 +284,14 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀 2 ) where
konst :: forall w. AbelianGroup w => w -> D3𝔸2 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndex [|| df ||] zeroMonomial )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀 2 -> w ) -> 𝕀 2 -> D3𝔸2 w
linearD f v =
let !o = origin @( 𝕀 Double ) @( T w )
in $$( monTabulate \ mon ->
in $$( monTabulateQ \ mon ->
if | isZeroMonomial mon
-> [|| f v ||]
| Just i <- isLinear mon
@ -319,14 +319,14 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀 3 ) where
konst :: forall w. AbelianGroup w => w -> D2𝔸3 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndex [|| df ||] zeroMonomial )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀 3 -> w ) -> 𝕀 3 -> D2𝔸3 w
linearD f v =
let !o = origin @( 𝕀 Double ) @( T w )
in $$( monTabulate \ mon ->
in $$( monTabulateQ \ mon ->
if | isZeroMonomial mon
-> [|| f v ||]
| Just i <- isLinear mon
@ -354,14 +354,14 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀 3 ) where
konst :: forall w. AbelianGroup w => w -> D3𝔸3 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndex [|| df ||] zeroMonomial )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀 3 -> w ) -> 𝕀 3 -> D3𝔸3 w
linearD f v =
let !o = origin @( 𝕀 Double ) @( T w )
in $$( monTabulate \ mon ->
in $$( monTabulateQ \ mon ->
if | isZeroMonomial mon
-> [|| f v ||]
| Just i <- isLinear mon
@ -389,14 +389,14 @@ instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀 4 ) where
konst :: forall w. AbelianGroup w => w -> D2𝔸4 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndex [|| df ||] zeroMonomial )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀 4 -> w ) -> 𝕀 4 -> D2𝔸4 w
linearD f v =
let !o = origin @( 𝕀 Double ) @( T w )
in $$( monTabulate \ mon ->
in $$( monTabulateQ \ mon ->
if | isZeroMonomial mon
-> [|| f v ||]
| Just i <- isLinear mon
@ -424,14 +424,14 @@ instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀 4 ) where
konst :: forall w. AbelianGroup w => w -> D3𝔸4 w
konst w =
let !o = fromInteger @w 0
in $$( monTabulate \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
in $$( monTabulateQ \ mon -> if isZeroMonomial mon then [|| w ||] else [|| o ||] )
value df = $$( monIndex [|| df ||] zeroMonomial )
value df = $$( monIndexQ [|| df ||] zeroMonomial )
linearD :: forall w. Module ( 𝕀 Double ) ( T w ) => ( 𝕀 4 -> w ) -> 𝕀 4 -> D3𝔸4 w
linearD f v =
let !o = origin @( 𝕀 Double ) @( T w )
in $$( monTabulate \ mon ->
in $$( monTabulateQ \ mon ->
if | isZeroMonomial mon
-> [|| f v ||]
| Just i <- isLinear mon

View file

@ -10,25 +10,29 @@ module Math.Linear
-- * Points and vectors (second version)
, (..), T(.., V2, V3, V4)
, Fin(..), MFin(..)
, Fin(..), MFin(..), universe, coordinates
, RepDim, RepresentableQ(..)
, Representable(..), set, injection, projection
, Vec(..), (!), find, zipIndices, unzip
, Vec(..), (!), find, zipIndices
) where
-- base
import Prelude hiding ( unzip )
import Control.Applicative
( ZipList(..) )
import Data.Coerce
( coerce )
import Data.Kind
( Type )
import GHC.Exts
( proxy# )
import GHC.Generics
( Generic, Generic1, Generically(..), Generically1(..) )
import GHC.Stack
( HasCallStack )
import GHC.TypeNats
( Nat, type (+) )
import Unsafe.Coerce
( unsafeCoerce )
( Nat, KnownNat, natVal' )
-- acts
import Data.Act
@ -114,53 +118,36 @@ pattern V4 x y z w = T ( 4 x y z w )
--------------------------------------------------------------------------------
infixr 5 `VS`
type Vec :: Nat -> Type -> Type
data Vec n a where
VZ :: Vec 0 a
VS :: forall n a. a -> Vec n a -> Vec ( 1 + n ) a
-- can't be strict, otherwise we can't conveniently
-- unsafeCoerce from lists
newtype Vec n a = Vec { vecList :: [ a ] }
deriving newtype ( Show, Eq, Ord, Functor, Foldable )
deriving Applicative via ZipList
--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 ] )
universe :: forall n. KnownNat n => Vec n ( Fin n )
universe = Vec [ Fin i | i <- [ 1 .. fromIntegral ( natVal' @n proxy# ) ] ]
deriving stock instance Functor ( Vec n )
deriving stock instance Foldable ( Vec n )
deriving stock instance Traversable ( Vec n )
instance Eq a => Eq ( Vec n a ) where
(==) = unsafeCoerce $ (==) @[a]
instance Ord a => Ord ( Vec n a ) where
compare = unsafeCoerce $ compare @[a]
(<=) = unsafeCoerce $ (<=) @[a]
coordinates :: forall r u. ( Representable r u ) => u -> Vec ( RepDim u ) r
coordinates u = fmap ( index u ) $ universe @( RepDim u )
{-# INLINEABLE coordinates #-}
infixl 9 !
(!) :: forall l a. Vec l a -> Fin l -> a
v ! Fin i = ( unsafeCoerce v :: [ a ] ) !! fromIntegral i
(!) :: forall l a. HasCallStack => Vec l a -> Fin l -> a
( Vec l ) ! Fin i = l !! ( fromIntegral i - 1 )
find :: forall l a. ( a -> Bool ) -> Vec l a -> MFin l
find f v = MFin ( find_ 1 v )
find f ( Vec v ) = MFin ( find_ 1 v )
where
find_ :: Word -> Vec n a -> Word
find_ j ( VS a as )
find_ :: Word -> [ a ] -> Word
find_ j ( a : as )
| f a
= j
| otherwise
= find_ ( j + 1 ) as
find_ _ VZ = 0
find_ _ [] = 0
zipIndices :: forall n a. Vec n a -> [ ( Fin n, a ) ]
zipIndices = zipIndices_ 1
zipIndices ( Vec v ) = zipIndices_ 1 v
where
zipIndices_ :: forall i. Word -> Vec i a -> [ ( Fin n, a ) ]
zipIndices_ _ VZ = []
zipIndices_ i (a `VS` as) = ( Fin i, a ) : zipIndices_ ( i + 1 ) as
unzip :: Vec n ( a, b ) -> ( Vec n a, Vec n b )
unzip VZ = ( VZ, VZ )
unzip ( ( a, b ) `VS` ab ) =
case unzip ab of
( as, bs ) -> ( a `VS` as, b `VS` bs )
zipIndices_ :: Word -> [ a ] -> [ ( Fin n, a ) ]
zipIndices_ _ [] = []
zipIndices_ i (a : as) = ( Fin i, a ) : zipIndices_ ( i + 1 ) as

View file

@ -1,5 +1,6 @@
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell #-}
module Math.Linear.Internal
( (..)
@ -22,7 +23,7 @@ import GHC.Generics
import GHC.Show
( showSpace )
import GHC.TypeNats
( Nat )
( Nat, KnownNat )
-- template-haskell
import Language.Haskell.TH
@ -94,7 +95,7 @@ instance Show ( 4 ) where
-- | 1, ..., n
type Fin :: Nat -> Type
newtype Fin n = Fin Word
deriving stock Eq
deriving stock ( Eq, Ord, Show )
-- | 0, ..., n
type MFin :: Nat -> Type
@ -109,16 +110,16 @@ class RepresentableQ r v | v -> r where
indexQ :: CodeQ v -> Fin ( RepDim v ) -> CodeQ r
type Representable :: Type -> Type -> Constraint
class Representable r v | v -> r where
class KnownNat ( RepDim v ) => Representable r v | v -> r where
tabulate :: ( Fin ( RepDim v ) -> r ) -> v
index :: v -> Fin ( RepDim v ) -> r
{-# INLINEABLE set #-}
set :: Representable r v => Fin ( RepDim v ) -> r -> v -> v
set i r u = tabulate \ j ->
if i == j
then r
else index u j
{-# INLINEABLE set #-}
projection :: ( Representable r u, Representable r v )
=> ( Fin ( RepDim v ) -> Fin ( RepDim u ) )

View file

@ -7,9 +7,9 @@
module Math.Monomial
( Mon(..)
, MonomialBasis(..), Deg, Vars
, MonomialBasis(..), MonomialBasisQ(..), Deg, Vars
, zeroMonomial, isZeroMonomial
, totalDegree, isLinear
, totalDegree, isLinear, linearMonomial
, split, mons
@ -57,40 +57,44 @@ type family Vars f
zeroMonomial :: forall k n. KnownNat n => Mon k n
zeroMonomial = Mon $ unsafeCoerce @[ Word ] @( Vec n Word )
$ replicate ( fromIntegral $ word @n ) 0
{-# INLINEABLE zeroMonomial #-}
isZeroMonomial :: Mon k n -> Bool
isZeroMonomial ( Mon VZ ) = True
isZeroMonomial ( Mon ( i `VS` is ) )
| i == 0
= isZeroMonomial ( Mon is )
| otherwise
= False
isZeroMonomial ( Mon ( Vec pows ) ) = all ( == 0 ) pows
totalDegree :: Mon k n -> Word
totalDegree ( Mon VZ ) = 0
totalDegree ( Mon ( i `VS` is ) ) = i + totalDegree ( Mon is )
totalDegree ( Mon ( Vec pows ) ) = sum pows
linearMonomial :: forall k n. ( KnownNat n, 1 <= k ) => Fin n -> Mon k n
linearMonomial ( Fin i' ) =
Mon $ unsafeCoerce @[ Word ] @( Vec n Word )
$ replicate ( i - 1 ) 0 ++ [ 1 ] ++ replicate ( n - i ) 0
where
n, i :: Int
n = fromIntegral $ word @n
i = fromIntegral i'
{-# INLINEABLE linearMonomial #-}
isLinear :: Mon k n -> Maybe ( Fin n )
isLinear = fmap Fin . go 1
isLinear ( Mon ( Vec pows ) ) = fmap Fin $ go 1 pows
where
go :: forall k' n'. Word -> Mon k' n' -> Maybe Word
go _ ( Mon VZ )
= Nothing
go j ( Mon ( i `VS` is ) )
go :: Word -> [ Word ] -> Maybe Word
go _ [] = Nothing
go j ( i : is )
| i == 0
= go ( j + 1 ) ( Mon is )
| i == 1 && isZeroMonomial ( Mon is )
= go ( j + 1 ) is
| i == 1 && all ( == 0 ) is
= Just j
| otherwise
= Nothing
-- | Comultiplication of monomials.
-- | Co-multiplication of monomials.
split :: Mon k n -> [ ( Mon k n, Mon k n ) ]
split ( Mon VZ ) = [ ( Mon VZ, Mon VZ ) ]
split ( Mon ( d `VS` ds ) ) =
[ ( Mon ( i `VS` as ), Mon ( ( d - i ) `VS` bs ) )
split ( Mon ( Vec [] ) ) = [ ( Mon ( Vec [] ), Mon ( Vec [] ) ) ]
split ( Mon ( Vec ( d : ds ) ) ) =
[ ( Mon ( Vec ( i : as ) ), Mon ( Vec ( ( ( d - i ) : bs ) ) ) )
| i <- [ 0 .. d ]
, ( Mon as, Mon bs ) <- ( split ( Mon ds ) )
, ( Mon ( Vec as ), Mon ( Vec bs ) ) <- ( split ( Mon ( Vec ds ) ) )
]
-- | All monomials of degree less than or equal to @k@ in @n@ variables,
@ -105,12 +109,12 @@ mons' 0 n = [ replicate ( fromIntegral n ) 0 ]
mons' k n = [ i : is | i <- reverse [ 0 .. k ], is <- mons' ( k - i ) ( n - 1 ) ]
subs :: Mon k n -> [ ( Mon k n, Word ) ]
subs ( Mon VZ ) = [ ( Mon VZ, maxBound ) ]
subs ( Mon ( i `VS` is ) )
= [ ( Mon ( j `VS` js )
subs ( Mon ( Vec []) ) = [ ( Mon ( Vec [] ), maxBound ) ]
subs ( Mon ( Vec ( i : is ) ) )
= [ ( Mon ( Vec ( j : js ) )
, if j == 0 then mult else min ( i `quot` j ) mult )
| j <- [ 0 .. i ]
, ( Mon js, mult ) <- subs ( Mon is )
, ( Mon ( Vec js ), mult ) <- subs ( Mon ( Vec is ) )
]
word :: forall n. KnownNat n => Word
@ -121,8 +125,8 @@ factorial :: Word -> Word
factorial i = product [ 1 .. i ]
vecFactorial :: Vec n Word -> Word
vecFactorial VZ = 1
vecFactorial ( i `VS` is ) = factorial i * vecFactorial is
vecFactorial ( Vec [] ) = 1
vecFactorial ( Vec ( i : is ) ) = factorial i * vecFactorial ( Vec is )
--------------------------------------------------------------------------------
-- Computations for the chain rule R^1 -> R^n -> R^m
@ -160,13 +164,13 @@ multiSubsetsSum :: forall n
multiSubsetsSum is = goMSS
where
goMSS :: forall i. Word -> Vec i Word -> [ Vec i [ ( Word, Word ) ] ]
goMSS 0 VZ = [ VZ ]
goMSS _ VZ = [ ]
goMSS s (n `VS` ns) =
[ multi `VS` rest
goMSS 0 ( Vec [] ) = [ Vec [] ]
goMSS _ ( Vec [] ) = [ ]
goMSS s ( Vec ( n : ns ) ) =
[ Vec ( multi : rest )
| s_i <- [ n * i_min .. s ]
, multi <- multiSubsetSum n s_i is
, rest <- goMSS ( s - s_i ) ns ]
, Vec rest <- goMSS ( s - s_i ) ( Vec ns ) ]
i_min = case is of
[] -> 0
_ -> max 0 $ minimum is
@ -224,19 +228,27 @@ subPart' m = zipWith ( \ i j -> i - m * j )
-- monomials in @Vars u@ variables, of degree up to (and including) @Deg u@.
type MonomialBasis :: ( Type -> Type ) -> Constraint
class MonomialBasis f where
monTabulate :: ( ( Mon ( Deg f ) ( Vars f ) ) -> CodeQ u ) -> CodeQ ( f u )
monIndex :: CodeQ ( f u ) -> Mon ( Deg f ) ( Vars f ) -> CodeQ u
monTabulate :: ( ( Mon ( Deg f ) ( Vars f ) ) -> u ) -> f u
monIndex :: f u -> Mon ( Deg f ) ( Vars f ) -> u
-- | @'MonomialBasisQ' f@ exhibits @f u@ as a free @r@-module with basis the
-- monomials in @Vars u@ variables, of degree up to (and including) @Deg u@,
-- at compile-time (to be spliced in using TemplateHaskell).
type MonomialBasisQ :: ( Type -> Type ) -> Constraint
class MonomialBasisQ f where
monTabulateQ :: ( ( Mon ( Deg f ) ( Vars f ) ) -> CodeQ u ) -> CodeQ ( f u )
monIndexQ :: CodeQ ( f u ) -> Mon ( Deg f ) ( Vars f ) -> CodeQ u
--------------------------------------------------------------------------------
prodRuleQ :: forall f r. MonomialBasis f
prodRuleQ :: forall f r. MonomialBasisQ f
=> CodeQ r -> CodeQ ( r -> r -> r ) -> CodeQ ( r -> r -> r )
-- Ring r constraint (circumvent TH constraint problem)
-> CodeQ ( f r ) -> CodeQ ( f r ) -> CodeQ ( f r )
prodRuleQ zero plus times df1 df2 =
monTabulate @f \ mon ->
monTabulateQ @f \ mon ->
[|| $$( foldQ plus zero
[ [|| $$times $$( monIndex @f df1 m1 ) $$( monIndex @f df2 m2 ) ||]
[ [|| $$times $$( monIndexQ @f df1 m1 ) $$( monIndexQ @f df2 m2 ) ||]
| (m1, m2) <- split mon
] )
||]

View file

@ -1,3 +1,5 @@
{-# LANGUAGE ScopedTypeVariables #-}
module Math.Root.Isolation
( -- * Root isolation using interval arithmetic
Box, BoxHistory
@ -22,7 +24,6 @@ module Math.Root.Isolation
-- ** Trees recording search space of root isolation algorithms
, RootIsolationTree(..), showRootIsolationTree
, RootIsolationStep(..)
, RootIsolationLeaf(..)
)
where
@ -33,19 +34,35 @@ import Control.Monad
( when )
import Data.Bifunctor
( Bifunctor(bimap) )
import Data.Coerce
( coerce )
import Data.Kind
( Type )
import Data.Foldable
( toList )
import Data.List
( nub, partition)
( nub, partition, sort )
import Data.List.NonEmpty
( NonEmpty )
import qualified Data.List.NonEmpty as NE
( NonEmpty(..), last )
( NonEmpty(..), last, singleton )
import Data.Semigroup
( Arg(..), Dual(..) )
import Numeric
( showFFloat )
import GHC.TypeNats
( Nat, KnownNat )
-- containers
import Data.Tree
( Tree(..) )
-- eigen
import qualified Eigen.Matrix as Eigen
( Matrix
, determinant, generate, inverse, unsafeCoeff
)
-- transformers
import Control.Monad.Trans.State.Strict as State
( State, evalState, get, put )
@ -54,7 +71,7 @@ import Control.Monad.Trans.Writer.CPS
-- MetaBrush
import Math.Algebra.Dual
( D, D1𝔸2(..) )
( D )
import Math.Epsilon
( nearZero )
import Math.Float.Utils
@ -63,6 +80,10 @@ import Math.Interval
import Math.Linear
import Math.Module
( Module(..) )
import Math.Monomial
( MonomialBasis(..), Deg, Vars
, linearMonomial, zeroMonomial
)
import qualified Math.Ring as Ring
--import Debug.Utils
@ -71,19 +92,19 @@ import qualified Math.Ring as Ring
-- | A tree recording the steps taken when doing cusp finding.
data RootIsolationTree d
= RootIsolationLeaf (RootIsolationLeaf d)
= RootIsolationLeaf String d
| RootIsolationStep RootIsolationStep [(d, RootIsolationTree d)]
-- | A description of a step taken in cusp-finding.
data RootIsolationStep
= GaussSeidelStep
| BisectionStep (String, Double)
| BisectionStep String Double
| Box1Step
| Box2Step
instance Show RootIsolationStep where
showsPrec _ GaussSeidelStep = showString "GS"
showsPrec _ ( BisectionStep (coord, w) )
showsPrec _ ( BisectionStep coord w )
= showString "bis "
. showParen True
( showString coord
@ -93,56 +114,81 @@ instance Show RootIsolationStep where
showsPrec _ Box1Step = showString "box(1)"
showsPrec _ Box2Step = showString "box(2)"
data RootIsolationLeaf d
= TooSmall d
deriving stock Show
showRootIsolationTree :: Box -> RootIsolationTree Box -> Tree String
showRootIsolationTree cand (RootIsolationLeaf l) = Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ show l) []
showRootIsolationTree
:: ( Representable Double ( n ), Show ( Box n ) )
=> Box n -> RootIsolationTree ( Box n ) -> Tree String
showRootIsolationTree cand (RootIsolationLeaf why l) = Node (show cand ++ " " ++ showArea (boxArea cand) ++ " " ++ why ++ " " ++ show l) []
showRootIsolationTree cand (RootIsolationStep s ts)
= Node (show cand ++ " abc " ++ showArea (boxArea cand) ++ " " ++ show s) $ map (\ (c,t) -> showRootIsolationTree c t) ts
boxArea :: Box -> Double
boxArea ( 𝕀 ( 2 t_lo s_lo ) ( 2 t_hi s_hi ) )
= abs ( t_hi - t_lo ) * abs ( s_hi - s_lo )
boxArea :: Representable Double ( n ) => Box n -> Double
boxArea ( 𝕀 lo hi ) =
product ( ( \ l h -> abs ( h - l ) ) <$> coordinates lo <*> coordinates hi )
showArea :: Double -> String
showArea area = "(area " ++ showFFloat (Just 6) area "" ++ ")"
type Box = 𝕀 2
type BoxHistory = [ NE.NonEmpty ( RootIsolationStep, Box ) ]
type Box n = 𝕀 n
type BoxHistory n = [ NE.NonEmpty ( RootIsolationStep, Box n ) ]
type BoxCt n d = ( n ~ 2, d ~ 3 )
{-
( Show ( 𝕀 n ), Show ( n )
, Eq ( n )
, Representable Double ( n )
, MonomialBasis ( D 1 ( n ) )
, Deg ( D 1 ( n ) ) ~ 1
, Vars ( D 1 ( n ) ) ~ n
, Module ( 𝕀 Double ) ( T ( 𝕀 n ) )
, Applicative ( Vec n )
, Ord ( d )
, Module Double ( T ( d ) )
, Representable Double ( d )
)
-}
-- | Options for the root isolation methods in 'findCusps' and 'isolateRootsIn'.
data RootIsolationOptions
type RootIsolationOptions :: Nat -> Nat -> Type
data RootIsolationOptions n d
= RootIsolationOptions
{ -- | Minimum width of a box.
minWidth :: !Double
-- | Given a box and its history, return a round of cusp-finding strategies
-- to run, in sequence, on this box.
, cuspFindingAlgorithms :: !( Box -> BoxHistory -> NonEmpty RootIsolationAlgorithm )
--
-- Returning an empty list means: give up on this box.
, cuspFindingAlgorithms :: !( Box n -> BoxHistory n -> Either String ( NE.NonEmpty ( RootIsolationAlgorithm n d ) ) )
}
data RootIsolationAlgorithm
type RootIsolationAlgorithm :: Nat -> Nat -> Type
data RootIsolationAlgorithm n d
-- | Bisection step.
= Bisection !BisectionOptions
= Bisection !( BisectionOptions n d )
-- | GaussSeidel step with the given preconditioning method.
| GaussSeidel !GaussSeidelOptions
| GaussSeidel !( GaussSeidelOptions n d )
-- | @box(1)@-consistency.
| Box1 !Box1Options
-- | @box(2)@-consistency.
| Box2 !Box2Options
-- | Options for the bisection method.
data BisectionOptions =
type BisectionOptions :: Nat -> Nat -> Type
data BisectionOptions n d =
BisectionOptions
{ canHaveSols :: !( ( 𝕀 2 -> D 1 ( 𝕀 2 ) ( 𝕀 3 ) ) -> Box -> Bool )
, fallbackBisectionDim :: !( [ ( RootIsolationStep, Box ) ] -> BoxHistory -> ( 𝕀 2 -> D 1 ( 𝕀 2 ) ( 𝕀 3 ) ) -> ( Int, String ) )
{ canHaveSols :: !( ( 𝕀 n -> D 1 ( 𝕀 n ) ( 𝕀 d ) ) -> Box n -> Bool )
, fallbackBisectionDim :: !( [ ( RootIsolationStep, Box n ) ] -> BoxHistory n -> ( 𝕀 n -> D 1 ( 𝕀 n ) ( 𝕀 d ) ) -> ( Fin n, String ) )
}
-- | Options for the interval GaussSeidel method.
data GaussSeidelOptions =
type GaussSeidelOptions :: Nat -> Nat -> Type
data GaussSeidelOptions n d =
GaussSeidelOptions
{ gsPreconditioner :: !Preconditioner }
{ -- | Which preconditioner to user?
gsPreconditioner :: !Preconditioner
-- | Function that projects over the equations we will consider
-- (the identity for a well-determined problem, or a project for
-- an overdetermined system).
, gsDims :: ( 𝕀 d -> 𝕀 n ) }
-- | Options for the @box(1)@-consistency method.
data Box1Options =
@ -156,7 +202,7 @@ data Box2Options =
, box2LambdaMin :: !Double
}
defaultRootIsolationOptions :: RootIsolationOptions
defaultRootIsolationOptions :: RootIsolationOptions 2 3
defaultRootIsolationOptions =
RootIsolationOptions
{ minWidth
@ -166,68 +212,119 @@ defaultRootIsolationOptions =
minWidth = 1e-5
narrowAbs = 5e-3
defaultRootIsolationAlgorithms :: Double -> Double -> Box -> BoxHistory -> NonEmpty RootIsolationAlgorithm
defaultRootIsolationAlgorithms minWidth narrowAbs box history =
case history of
lastRoundBoxes : _
-- If, in the last round of strategies, we didn't try bisection...
| any ( \case { BisectionStep {} -> False; _ -> True } . fst ) lastRoundBoxes
, (_, lastRoundFirstBox ) <- NE.last lastRoundBoxes
-- ...and the last round didn't sufficiently reduce the size of the box...
, not $ box `sufficientlySmallerThan` lastRoundFirstBox
-- ...then try bisecting the box.
-> Bisection ( defaultBisectionOptions minWidth narrowAbs box ) NE.:| []
-- Otherwise, do a normal round.
-- Currently: we try an interval GaussSeidel step followed by box(1)-consistency.
_ -> GaussSeidel defaultGaussSeidelOptions
NE.:| [ Box1 ( Box1Options { box1EpsEq = narrowAbs } ) ]
defaultRootIsolationAlgorithms
:: forall n d
. BoxCt n d
=> Double -> Double -> Box n -> BoxHistory n -> Either String ( NE.NonEmpty ( RootIsolationAlgorithm n d ) )
defaultRootIsolationAlgorithms minWidth narrowAbs box history
-- All the box widths are lower than the minimum width: give up.
| and $ ( \ cd -> width cd <= minWidth ) <$> coordinates box
= Left $ "sizes <= " ++ show minWidth
-- Give up after depth 50
-- | length history >= 50
-- = Left "depth >= 50"
| otherwise
= Right $ case history of
lastRoundBoxes : _
-- If, in the last round of strategies, we didn't try bisection...
| any ( \case { BisectionStep {} -> False; _ -> True } . fst ) lastRoundBoxes
, ( _step, lastRoundFirstBox ) <- NE.last lastRoundBoxes
-- ...and the last round didn't sufficiently reduce the size of the box...
, not $ box `sufficientlySmallerThan` lastRoundFirstBox
-- ...then try bisecting the box.
-> NE.singleton $ Bisection ( defaultBisectionOptions minWidth narrowAbs box )
-- Otherwise, do a normal round.
-- Currently: we try an interval GaussSeidel step followed by box(1)-consistency.
_ -> GaussSeidel defaultGaussSeidelOptions
NE.:| [ Box1 ( Box1Options { box1EpsEq = narrowAbs } ) ]
where
sufficientlySmallerThan :: Box -> Box -> Bool
sufficientlySmallerThan
( 𝕀 ( 2 t1_lo s1_lo ) ( 2 t1_hi s1_hi ) )
( 𝕀 ( 2 t2_lo s2_lo ) ( 2 t2_hi s2_hi ) ) =
( ( t1_hi - t1_lo ) - ( t2_hi - t2_lo ) > narrowAbs )
||
( ( s1_hi - s1_lo ) - ( s2_hi - s2_lo ) > narrowAbs )
-- Did we reduce the box width by at least "narrowAbs" in at least one of the dimensions?
sufficientlySmallerThan :: Box n -> Box n -> Bool
b1 `sufficientlySmallerThan` b2 =
or $
( \ cd1 cd2 -> width cd2 - width cd1 > narrowAbs )
<$> coordinates b1
<*> coordinates b2
{-# INLINEABLE defaultRootIsolationAlgorithms #-}
defaultGaussSeidelOptions :: GaussSeidelOptions
defaultGaussSeidelOptions = GaussSeidelOptions { gsPreconditioner = InverseMidJacobian }
defaultGaussSeidelOptions :: GaussSeidelOptions 2 3
defaultGaussSeidelOptions =
GaussSeidelOptions
{ gsPreconditioner = InverseMidJacobian
, gsDims = --id
\ ( 𝕀 ( 3 _a_lo b_lo c_lo ) ( 3 _a_hi b_hi c_hi ) )
-> 𝕀 ( 2 b_lo c_lo ) ( 2 b_hi c_hi )
}
defaultBisectionOptions :: Double -> Double -> Box -> BisectionOptions
defaultBisectionOptions
_minWidth
_narrowAbs
box@( 𝕀 ( 2 t_lo s_lo ) ( 2 t_hi s_hi ) )
= BisectionOptions
{ canHaveSols =
\ eqs box' ->
-- box(0)-consistency
let D12 iRange _ _ = eqs box'
in 3 0 0 0 `inside` iRange
:: forall n d
. BoxCt n d
=> Double -> Double
-> Box n -> BisectionOptions n d
defaultBisectionOptions _minWidth _narrowAbs box =
BisectionOptions
{ canHaveSols =
\ eqs box' ->
-- box(0)-consistency
let iRange' :: Box d
iRange' = ( `monIndex` zeroMonomial ) $ eqs box'
in unT ( origin @Double ) `inside` iRange'
-- box(1)-consistency
--not $ null $ makeBox1Consistent eqs _minWidth _narrowAbs box'
-- box(1)-consistency
--not $ null $ makeBox1Consistent eqs _minWidth _narrowAbs box'
-- box(2)-consistency
--let box'' = makeBox2Consistent eqs _minWidth _narrowAbs 0.2 box'
-- D12 iRange'' = eqs box''
--in 3 0 0 0 `inside` iRange''
, fallbackBisectionDim =
\ _roundHist _prevRoundsHist eqs ->
let D12 _ ( T f_t ) ( T f_s ) = eqs box
wd_t = t_hi - t_lo
wd_s = s_hi - s_lo
tJWidth = wd_t * normVI3 f_t
sJWidth = wd_s * normVI3 f_s
in if tJWidth >= sJWidth
-- bisect along dimension that maximises width(x_j) * || J_j || ...
-- ... but don't allow thin boxes
|| ( wd_t >= 10 * wd_s )
&& not ( wd_s >= 10 * wd_t )
then (0, "")
else (1, "")
}
-- box(2)-consistency
--let box'' = makeBox2Consistent eqs _minWidth _narrowAbs 0.2 box'
-- iRange'' :: Box d
-- iRange'' = value @Double @1 @( Box n ) $ eqs box''
--in origin @Double `inside` iRange''
, fallbackBisectionDim =
\ _roundHist _prevRoundsHist eqs ->
let df = eqs box
datPerCoord =
[ CoordBisectionData
{ coordIndex = i
, coordInterval = box `index` i
, coordJacobianColumn = df `monIndex` ( linearMonomial i )
}
| i <- toList $ universe @n ]
-- First, check if the largest dimension is over 10 times larger
-- than the smallest dimension; if so bisect along that coordinate.
in case sortOnArg ( width . coordInterval ) datPerCoord of
[] -> error "dimension 0"
[Arg _ d] -> (coordIndex d, "")
Arg w0 _ : ds ->
let Arg w1 d1 = last ds
in if w1 >= 10 * w0
then ( coordIndex d1, "tooWide")
-- Otherwise, pick the dimension with the largest spread = width * Jacobian column norm
else case sortOnArg ( Dual . spread ) datPerCoord of
[] -> error "dimension 0"
Arg _ d : _ ->
(coordIndex d, "spread")
}
{-# INLINEABLE defaultBisectionOptions #-}
sortOnArg :: Ord b => (a -> b) -> [a] -> [Arg b a]
sortOnArg f = sort . map ( \ a -> Arg ( f a ) a )
{-# INLINEABLE sortOnArg #-}
type CoordBisectionData :: Nat -> Nat -> Type
data CoordBisectionData n d =
CoordBisectionData
{ coordIndex :: !( Fin n )
, coordInterval :: !( 𝕀 Double )
, coordJacobianColumn :: !( 𝕀 d )
}
deriving stock instance Show ( d )
=> Show ( CoordBisectionData n d )
spread :: ( BoxCt n d, Representable Double ( d ) )
=> CoordBisectionData n d -> Double
spread ( CoordBisectionData { coordInterval = cd, coordJacobianColumn = j_cd } )
= width cd * normVI j_cd
{-# INLINEABLE spread #-}
-- | Use the following algorithms using interval arithmetic in order
-- to isolate roots:
@ -243,10 +340,12 @@ defaultBisectionOptions
-- will converge starting from anywhere inside the box), and @dunno@ are small
-- boxes which might or might not contain solutions.
isolateRootsIn
:: RootIsolationOptions
-> ( 𝕀 2 -> D 1 ( 𝕀 2 ) ( 𝕀 3 ) )
-> Box
-> ( [ ( Box, RootIsolationTree Box ) ], ( [ Box ], [ Box ] ) )
:: forall n d
. BoxCt n d
=> RootIsolationOptions n d
-> ( 𝕀 n -> D 1 ( 𝕀 n ) ( 𝕀 d ) )
-> Box n
-> ( [ ( Box n, RootIsolationTree ( Box n ) ) ], ( [ Box n ], [ Box n ] ) )
isolateRootsIn
( RootIsolationOptions
{ minWidth
@ -257,40 +356,41 @@ isolateRootsIn
where
go :: BoxHistory
-> Box -- box to work on
-> Writer ( [ Box ], [ Box ] )
[ RootIsolationTree Box ]
go history
cand@( 𝕀 ( 2 t_lo s_lo ) ( 2 t_hi s_hi ) )
| -- Check the envelope equation interval contains zero.
not $ ( 3 0 0 0 `inside` iRange )
-- Box doesn't contain a solution: discard it.
go :: BoxHistory n
-> Box n -- box to work on
-> Writer ( [ Box n ], [ Box n ] )
[ RootIsolationTree ( Box n ) ]
go history cand
| -- Check the range of the equations contains zero.
not $ ( unT ( origin @Double ) `inside` iRange )
-- Box doesn't contain a solution: discard it.
= return []
-- Box is small: stop processing it.
| t_hi - t_lo < minWidth && s_hi - s_lo < minWidth
= do let res = TooSmall cand
tell ( [ cand ], [] )
return [ RootIsolationLeaf res ]
| otherwise
= doStrategies history ( cuspFindingAlgorithms cand history ) cand
= case cuspFindingAlgorithms cand history of
Right strats -> doStrategies history strats cand
Left whyStop -> do
-- We are giving up on this box (e.g. because it is too small,
-- or we have reached an iteration depth).
tell ( [ cand ], [] )
return [ RootIsolationLeaf whyStop cand ]
where
D12 iRange _ _ = eqs $ 𝕀 ( 2 t_lo s_lo ) ( 2 t_hi s_hi )
iRange :: Box d
iRange = ( `monIndex` zeroMonomial ) $ eqs cand
-- Run a round of cusp finding strategies, then recur.
doStrategies
:: BoxHistory
-> NonEmpty RootIsolationAlgorithm
-> Box
-> Writer ( [ Box ], [ Box ] )
[ RootIsolationTree Box ]
:: BoxHistory n
-> NonEmpty ( RootIsolationAlgorithm n d )
-> Box n
-> Writer ( [ Box n ], [ Box n ] )
[ RootIsolationTree ( Box n ) ]
doStrategies prevRoundsHist = do_strats []
where
do_strats :: [ ( RootIsolationStep, Box ) ]
-> NE.NonEmpty RootIsolationAlgorithm
-> Box
-> Writer ( [ Box ], [ Box ] )
[ RootIsolationTree Box ]
do_strats :: [ ( RootIsolationStep, Box n ) ]
-> NE.NonEmpty ( RootIsolationAlgorithm n d )
-> Box n
-> Writer ( [ Box n ], [ Box n ] )
[ RootIsolationTree ( Box n )]
do_strats thisRoundHist ( algo NE.:| algos ) box = do
-- Run one strategy in the round.
( step, boxes ) <- doStrategy thisRoundHist prevRoundsHist eqs minWidth algo box
@ -307,116 +407,137 @@ isolateRootsIn
in recur step ( go history ) boxes
recur :: RootIsolationStep
-> ( Box -> Writer ( [ Box ], [ Box ] ) [ RootIsolationTree Box ] )
-> [ Box ]
-> Writer ( [Box], [Box] ) [ RootIsolationTree Box ]
-> ( Box n -> Writer ( [ Box n ], [ Box n ] ) [ RootIsolationTree ( Box n ) ] )
-> [ Box n ]
-> Writer ( [ Box n ], [ Box n ] ) [ RootIsolationTree ( Box n ) ]
recur step r cands = do
rest <- traverse ( \ c -> do { trees <- r c; return [ (c, tree) | tree <- trees ] } ) cands
return [ RootIsolationStep step $ concat rest ]
{-# INLINEABLE isolateRootsIn #-}
-- | Execute a cusp-finding strategy, replacing the input box with
-- (hopefully smaller) output boxes.
doStrategy :: [ ( RootIsolationStep, Box ) ]
-> BoxHistory
-> ( 𝕀 2 -> D 1 ( 𝕀 2 ) ( 𝕀 3 ) )
-> Double
-> RootIsolationAlgorithm
-> Box
-> Writer ( [ Box ], [ Box ] )
( RootIsolationStep, [ Box ] )
doStrategy
:: BoxCt n d
=> [ ( RootIsolationStep, Box n ) ]
-> BoxHistory n
-> ( 𝕀 n -> D 1 ( 𝕀 n ) ( 𝕀 d ) )
-> Double
-> RootIsolationAlgorithm n d
-> Box n
-> Writer ( [ Box n ], [ Box n ] )
( RootIsolationStep, [ Box n ] )
doStrategy roundHistory previousRoundsHistory eqs minWidth algo box =
case algo of
GaussSeidel ( GaussSeidelOptions { gsPreconditioner } ) -> do
boxes <- intervalGaussSeidel eqs gsPreconditioner box
GaussSeidel gsOptions -> do
boxes <- intervalGaussSeidel gsOptions eqs box
return ( GaussSeidelStep, boxes )
Box1 ( Box1Options { box1EpsEq } ) ->
return ( Box1Step, makeBox1Consistent eqs minWidth box1EpsEq box )
Box2 ( Box2Options { box2LambdaMin, box2EpsEq } ) ->
return ( Box2Step, [ makeBox2Consistent eqs minWidth box2EpsEq box2LambdaMin box ] )
Bisection ( BisectionOptions { canHaveSols, fallbackBisectionDim } ) -> do
let ( boxes, whatBis ) = bisect ( canHaveSols eqs ) ( fallbackBisectionDim roundHistory previousRoundsHistory eqs ) box
return ( BisectionStep whatBis, boxes )
let ( boxes, ( whatBis, mid ) ) = bisect ( canHaveSols eqs ) ( fallbackBisectionDim roundHistory previousRoundsHistory eqs ) box
return ( BisectionStep whatBis mid, boxes )
{-# INLINEABLE doStrategy #-}
-- | Bisect the given box.
--
-- (The difficult part lies in determining along which dimension to bisect.)
bisect :: ( Box -> Bool )
-- ^ how to check whether a box contains solutions
-> ( Int, String )
-- ^ fall-back choice of dimension (and "why" we chose it)
-> Box
-> ( [ Box ], ( String, Double ) )
bisect canHaveSols ( dim, why )
( 𝕀 ( 2 t_lo s_lo ) ( 2 t_hi s_hi ) )
bisect
:: forall n
. ( KnownNat n, Representable Double ( n ) )
=> ( Box n -> Bool )
-- ^ how to check whether a box contains solutions
-> ( Fin n, String )
-- ^ fallback choice of dimension (and "why" we chose it)
-> Box n
-> ( [ Box n ], ( String, Double ) )
bisect canHaveSols ( fallbackDim, why ) box =
-- We try to bisect along a dimension which eliminates zeros from one of the
-- sub-regions.
--
-- If this fails, we fall-back to the provided dimension choice.
= let t_mid = 0.5 * ( t_lo + t_hi )
s_mid = 0.5 * ( s_lo + s_hi )
l = 𝕀 ( 2 t_lo s_lo ) ( 2 t_mid s_hi )
r = 𝕀 ( 2 t_mid s_lo ) ( 2 t_hi s_hi )
d = 𝕀 ( 2 t_lo s_lo ) ( 2 t_hi s_mid )
u = 𝕀 ( 2 t_lo s_mid ) ( 2 t_hi s_hi )
l_skip = not $ canHaveSols l
r_skip = not $ canHaveSols r
d_skip = not $ canHaveSols d
u_skip = not $ canHaveSols u
-- If this fails, we fall back to the provided dimension choice.
case findFewestSols solsList of
Just ( Arg _nbSols ( i, mid, oks ) ) ->
( oks, ("cd = " ++ show i, mid ) )
Nothing ->
case bisectInCoord box fallbackDim of
( mid, (lo, hi) ) ->
( [ lo, hi ], ( why, mid ) )
where
solsList =
[ Arg ( fromIntegral $ length oks ) ( i, mid, oks )
| i <- toList $ universe @n
, let (mid, (lo, hi)) = bisectInCoord box i
lo_ok = canHaveSols lo
hi_ok = canHaveSols hi
oks = [ lo | lo_ok ] ++ [ hi | hi_ok ]
]
{-# INLINEABLE bisect #-}
( bisGuesses, whatBis )
| l_skip && r_skip
= ( [], ( "lr", t_mid ) )
| d_skip && u_skip
= ( [], ( "du", s_mid ) )
| l_skip
= ( [ r ], ( "r", t_mid ) )
| r_skip
= ( [ l ], ( "l", t_mid ) )
| d_skip
= ( [ u ], ( "u", s_mid ) )
| u_skip
= ( [ d ], ( "d", s_mid ) )
| otherwise
= let why' = case why of
"" -> ""
_ -> " (" ++ why ++ ")"
in case dim of
0 -> ( [ l, r ], ( "t" ++ why', t_mid ) )
1 -> ( [ d, u ], ( "s" ++ why', s_mid ) )
_ -> error "bisect: fall-back bisection dimension should be either 0 or 1"
in ( bisGuesses, whatBis )
-- | Find any element with the least argument.
--
-- NB: this function shortcuts as soon as it finds an element with argument 0.
findFewestSols :: forall a. [ Arg Word a ] -> Maybe ( Arg Word a )
findFewestSols [] = Nothing
findFewestSols ( arg@( Arg nbSols _ ) : args )
| nbSols == 0
= Just arg
| otherwise
= Just $ go arg args
where
go :: Arg Word a -> [ Arg Word a ] -> Arg Word a
go bestSoFar [] = bestSoFar
go bestSoFar@( Arg bestNbSolsSoFar _ ) ( arg'@( Arg nbSols' _ ) : args' )
| nbSols' == 0
= arg'
| nbSols' < bestNbSolsSoFar
= go arg' args'
| otherwise
= go bestSoFar args'
bisectInCoord
:: Representable Double ( n )
=> Box n -> Fin n -> ( Double, ( Box n, Box n ) )
bisectInCoord box i =
let 𝕀 z_lo z_hi = box `index` i
z_mid = 0.5 * ( z_lo + z_hi )
in ( z_mid
, ( set i ( 𝕀 z_lo z_mid ) box
, set i ( 𝕀 z_mid z_hi ) box ) )
{-# INLINEABLE bisectInCoord #-}
-- | Interval Newton method with GaussSeidel step.
intervalGaussSeidel
:: ( 𝕀 2 -> D 1 ( 𝕀 2 ) ( 𝕀 3 ) )
-- ^ equations
-> Preconditioner
-- ^ preconditioner to use for the GaussSeidel step
-> 𝕀 2
-> Writer ( [ 𝕀 2 ], [ 𝕀 2 ] ) [ 𝕀 2 ]
:: forall n d
. BoxCt n d
=> GaussSeidelOptions n d
-> ( 𝕀 n -> D 1 ( 𝕀 n ) ( 𝕀 d ) )
-- ^ equations
-> 𝕀 n
-- ^ box
-> Writer ( [ 𝕀 n ], [ 𝕀 n ] ) [ 𝕀 n ]
intervalGaussSeidel
( GaussSeidelOptions { gsPreconditioner = precondMeth, gsDims = projectDims } )
eqs
precondMethod
ts@( 𝕀 ( 2 t_lo s_lo ) ( 2 t_hi s_hi ) )
| D12 _f ( T ee_f_t ) ( T ee_f_s )
<- eqs ts
, let f_t = take2 ee_f_t
f_s = take2 ee_f_s
, D12 ee_f_mid ( T _ee_f_t_mid ) ( T _ee_f_s_mid )
<- eqs ts_mid
, let f_mid = take2 ee_f_mid
box
| let boxMid = singleton $ midpoint box
df = eqs box
f' :: Vec n ( 𝕀 n )
f' = fmap ( \ i -> projectDims $ df `monIndex` linearMonomial i ) ( universe @n )
f_mid = projectDims $ eqs boxMid `monIndex` zeroMonomial
= let -- Interval Newton method: take one GaussSeidel step
-- for the equation f'(x) v = - f(x_mid),
-- where f = 𝛿E/𝛿s * dc/dt
( a, b ) = precondition precondMethod
( midI f_t, midI f_s )
( f_t, f_s ) ( neg f_mid )
-- for the system of equations f'(x) ( x - x_mid ) = - f(x_mid).
( a, b ) = precondition precondMeth
( fmap midpoint f' )
f' ( negateBox $ singleton $ midpoint $ f_mid )
-- NB: we have to re-center around zero to take a GaussSeidel step.
gsGuesses = map ( first ( \ ts' -> unT $ T ts' ^+^ T ts_mid ) )
$ gaussSeidelStep a b ( unT $ T ts ^-^ T ts_mid )
-- NB: we have to change coordinates, putting the midpoint of the box
-- at the origin, in order to take a GaussSeidel step.
gsGuesses = map ( first ( \ box' -> unT $ box' ^+^ T boxMid ) )
$ gaussSeidelStep a b ( T box ^-^ T boxMid )
in
-- If the GaussSeidel step was a contraction, then the box
-- contains a unique solution (by the Banach fixed point theorem).
@ -425,57 +546,58 @@ intervalGaussSeidel
-- Newton's method is guaranteed to converge to the unique solution.
let !(done, todo) = bimap ( map fst ) ( map fst )
$ partition snd gsGuesses
in do tell ([], done)
in do tell ( [], done )
return todo
where
t_mid = 0.5 * ( t_lo + t_hi )
s_mid = 0.5 * ( s_lo + s_hi )
ts_mid = singleton ( 2 t_mid s_mid )
neg ( 𝕀 ( 2 x_lo y_lo ) ( 2 x_hi y_hi ) )
= let !( 𝕀 x'_lo x'_hi ) = negate $ 𝕀 x_lo x_hi
!( 𝕀 y'_lo y'_hi ) = negate $ 𝕀 y_lo y_hi
in 𝕀 ( 2 x'_lo y'_lo ) ( 2 x'_hi y'_hi )
take2 :: 𝕀 3 -> 𝕀 2
take2 ( 𝕀 ( 3 _a_lo b_lo c_lo ) ( 3 _a_hi b_hi c_hi ) )
= 𝕀 ( 2 b_lo c_lo ) ( 2 b_hi c_hi )
-- TODO: always using the last 2 coordinates, instead of any
-- pair of two coordinates.
{-# INLINEABLE intervalGaussSeidel #-}
-- Take one interval GaussSeidel step for the equation \( A X = B \),
-- refining the initial guess box for \( X \) into up to four (disjoint) new boxes.
-- | The midpoint of a box.
midpoint :: Representable Double ( n ) => 𝕀 n -> n
midpoint box =
tabulate $ \ i ->
let 𝕀 z_lo z_hi = index box i
z_mid = 0.5 * ( z_lo + z_hi )
in z_mid
{-# INLINEABLE midpoint #-}
-- | Negate a box.
negateBox :: Representable Double ( n ) => 𝕀 n -> 𝕀 n
negateBox v = tabulate $ \ i -> negate ( index v i )
{-# INLINEABLE negateBox #-}
-- | Take one interval GaussSeidel step for the equation \( A X = B \),
-- refining the initial guess box for \( X \) into up to \( 2^n \) (disjoint) new boxes.
--
-- The boolean indicates whether the GaussSeidel step was a contraction.
gaussSeidelStep :: ( 𝕀 2, 𝕀 2 ) -- ^ columns of \( A \)
-> 𝕀 2 -- ^ \( B \)
-> 𝕀 2 -- ^ initial box \( X \)
-> [ ( 𝕀 2, Bool ) ]
gaussSeidelStep
( 𝕀 ( 2 a11_lo a21_lo ) ( 2 a11_hi a21_hi )
, 𝕀 ( 2 a12_lo a22_lo ) ( 2 a12_hi a22_hi ) )
( 𝕀 ( 2 b1_lo b2_lo ) ( 2 b1_hi b2_hi ) )
( 𝕀 ( 2 x1_lo x2_lo ) ( 2 x1_hi x2_hi ) )
= let !a11 = 𝕀 a11_lo a11_hi
!a12 = 𝕀 a12_lo a12_hi
!a21 = 𝕀 a21_lo a21_hi
!a22 = 𝕀 a22_lo a22_hi
!b1 = 𝕀 b1_lo b1_hi
!b2 = 𝕀 b2_lo b2_hi
!x1 = 𝕀 x1_lo x1_hi
!x2 = 𝕀 x2_lo x2_hi
in nub $ do
-- x1' = ( ( b1 - a12 * x2 ) / a11 ) ∩ x1
x1'_0 <- ( b1 - a12 * x2 ) `extendedDivide` a11
( x1'@( 𝕀 x1'_lo x1'_hi ), sub_x1 ) <- x1'_0 `intersect` x1
-- x2' = ( ( b2 - a21 * x1' ) / a22 ) ∩ x2
x2'_0 <- ( b2 - a21 * x1' ) `extendedDivide` a22
( 𝕀 x2'_lo x2'_hi, sub_x2 ) <- x2'_0 `intersect` x2
return ( ( 𝕀 ( 2 x1'_lo x2'_lo) ( 2 x1'_hi x2'_hi ) )
, sub_x1 && sub_x2 )
:: forall n
. ( Representable Double ( n ), Eq ( n ) )
=> Vec n ( 𝕀 n ) -- ^ columns of \( A \)
-> 𝕀 n -- ^ \( B \)
-> T ( 𝕀 n ) -- ^ initial box \( X \)
-> [ ( T ( 𝕀 n ), Bool ) ]
gaussSeidelStep as b ( T x0 ) = coerce $ nub $
forEachDim @n ( x0, True ) $ \ i ( x, contraction ) -> do
-- x_i' = ( b_i - sum { j /= i } a_ij * x_j ) / a_ii
x_i'0 <- ( b `index` i - sum [ ( as ! j ) `index` i * x `index` j | j <- toList ( universe @n ), j /= i ] )
`extendedDivide` ( ( as ! i ) `index` i )
( x_i', sub_i ) <- x_i'0 `intersect` ( x `index` i )
return $ ( set i x_i' x, sub_i && contraction )
-- TODO: try implementing the complete interval union GaussSeidel algorithm.
-- See "Algorithm 2" in
-- "Using interval unions to solve linear systems of equations with uncertainties"
{-# INLINEABLE gaussSeidelStep #-}
-- | Run an effectful computation several times in sequence, piping its output
-- into the next input, once for each coordinate dimension.
forEachDim :: forall n a m. ( KnownNat n, Monad m ) => a -> ( Fin n -> a -> m a ) -> m a
forEachDim a0 f = go ( toList $ universe @n ) a0
where
go [] a = return a
go ( i : is ) a = do
a' <- f i a
go is a'
{-# INLINEABLE forEachDim #-}
-- | Compute the intersection of two intervals.
--
@ -503,9 +625,10 @@ data Preconditioner
-- See also
-- "Presentation of a highly tuned multithreaded interval solver for underdetermined and well-determined nonlinear systems"
makeBox1Consistent
:: ( 𝕀 2 -> D 1 ( 𝕀 2 ) ( 𝕀 3 ) )
:: BoxCt n d
=> ( 𝕀 n -> D 1 ( 𝕀 n ) ( 𝕀 d ) )
-> Double -> Double
-> Box -> [ Box ]
-> Box n -> [ Box n ]
makeBox1Consistent eqs minWidth epsEq x =
( `State.evalState` False ) $
pipeFunctionsWhileTrue ( allNarrowingOperators epsEq minWidth eqs ) x
@ -513,24 +636,26 @@ makeBox1Consistent eqs minWidth epsEq x =
-- | An implementation of "bound-consistency" from the paper
-- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
makeBox2Consistent
:: ( 𝕀 2 -> D 1 ( 𝕀 2 ) ( 𝕀 3 ) )
:: forall n d
. BoxCt n d
=> ( 𝕀 n -> D 1 ( 𝕀 n ) ( 𝕀 d ) )
-> Double -> Double -> Double
-> Box -> Box
-> Box n -> Box n
makeBox2Consistent eqs minWidth epsEq lambdaMin x0 = ( `State.evalState` False ) $ doLoop 0.25 x0
where
doLoop :: Double -> Box -> State Bool Box
doLoop :: Double -> Box n -> State Bool ( Box n )
doLoop lambda x = do
x' <- boundConsistency get_t set_t lambda x
x'' <- boundConsistency get_s set_s lambda x'
x'' <- forEachDim @n x $ \ i ->
boundConsistency ( `index` i ) ( set i ) lambda
modified <- State.get
let lambda' = if modified then lambda else 0.5 * lambda
if lambda' < lambdaMin
then return x''
else do { State.put False ; doLoop lambda' x'' }
boundConsistency :: ( Box -> 𝕀 Double )
-> ( 𝕀 Double -> Box -> Box )
-> Double -> Box -> State Bool Box
boundConsistency :: ( Box n -> 𝕀 Double )
-> ( 𝕀 Double -> Box n -> Box n )
-> Double -> Box n -> State Bool ( Box n )
boundConsistency getter setter lambda box = do
let x@( 𝕀 x_inf x_sup ) = getter box
c1 = ( 1 - lambda ) * x_inf + lambda * x_sup
@ -548,62 +673,63 @@ makeBox2Consistent eqs minWidth epsEq lambdaMin x0 = ( `State.evalState` False )
State.put True
return $ setter x' box
precondition :: Preconditioner
-> ( 𝕀 2, 𝕀 2 )
-> ( 𝕀 2, 𝕀 2 )
-> 𝕀 2
-> ( ( 𝕀 2, 𝕀 2 ), 𝕀 2 )
precondition meth jac_mid a@( a1, a2 ) b =
-- | Pre-condition the system \( AX = B \).
precondition
:: forall n
. ( KnownNat n, Representable Double ( n ) )
=> Preconditioner -- ^ pre-conditioning method to use
-> Vec n ( n ) -- ^ entry-wise midpoint matrix of the interval Jacobian matrix
-> Vec n ( 𝕀 n ) -- ^ columns of \( A \)
-> 𝕀 n -- ^ \( B \)
-> ( Vec n ( 𝕀 n ), 𝕀 n )
precondition meth jac_mid as b =
case meth of
NoPreconditioning
-> ( a, b )
-> ( as, b )
InverseMidJacobian
| ( 𝕀 ( 2 a11_lo a21_lo ) ( 2 a11_hi a21_hi )
, 𝕀 ( 2 a12_lo a22_lo ) ( 2 a12_hi a22_hi ) ) <- jac_mid
, let !a11 = 0.5 * ( a11_lo + a11_hi )
!a12 = 0.5 * ( a12_lo + a12_hi )
!a21 = 0.5 * ( a21_lo + a21_hi )
!a22 = 0.5 * ( a22_lo + a22_hi )
!d = a11 * a22 - a12 * a21
, not ( nearZero d )
, let !precond = ( 2 a22 -a12, 2 -a21 a11 )
!inv = recip d
f x = scale inv $ matMulVec precond x
-> ( ( f a1, f a2 ), f b )
| let mat = toEigen jac_mid
det = Eigen.determinant mat
, not $ nearZero det
, let precond = Eigen.inverse mat
doPrecond = matMulVec ( fromEigen precond )
-> ( fmap doPrecond as, doPrecond b )
| otherwise
-> ( a, b )
scale :: Double -> 𝕀 2 -> 𝕀 2
scale s ( 𝕀 ( 2 a1_lo a2_lo ) ( 2 a1_hi a2_hi ) )
= 𝕀 ( 2 b1_lo b2_lo ) ( 2 b1_hi b2_hi )
-> ( as, b )
where
𝕀 b1_lo b1_hi = scaleInterval s ( 𝕀 a1_lo a1_hi )
𝕀 b2_lo b2_hi = scaleInterval s ( 𝕀 a2_lo a2_hi )
toEigen :: Vec n ( n ) -> Eigen.Matrix n n Double
toEigen cols =
Eigen.generate $ \ r c ->
( cols ! Fin ( fromIntegral c + 1 ) ) `index` ( Fin ( fromIntegral r + 1 ) )
matMulVec :: ( 2, 2 ) -> 𝕀 2 -> 𝕀 2
matMulVec ( 2 a11 a21, 2 a12 a22 ) ( 𝕀 ( 2 u_lo v_lo ) ( 2 u_hi v_hi ) )
= 𝕀 ( 2 u'_lo v'_lo ) ( 2 u'_hi v'_hi )
where
u = 𝕀 u_lo u_hi
v = 𝕀 v_lo v_hi
𝕀 u'_lo u'_hi = scaleInterval a11 u + scaleInterval a12 v
𝕀 v'_lo v'_hi = scaleInterval a21 u + scaleInterval a22 v
fromEigen :: Eigen.Matrix n n Double -> Vec n ( n )
fromEigen mat =
fmap
( \ ( Fin c ) ->
tabulate $ \ ( Fin r ) ->
Eigen.unsafeCoeff ( fromIntegral r - 1 ) ( fromIntegral c - 1 ) mat
)
( universe @n )
midI :: 𝕀 2 -> 𝕀 2
midI ( 𝕀 ( 2 x_lo y_lo ) ( 2 x_hi y_hi ) ) =
singleton $ 2 ( 0.5 * ( x_lo + x_hi ) ) ( 0.5 * ( y_lo + y_hi ) )
--width :: 𝕀 1 -> Double
--width ( 𝕀 ( 1 lo ) ( 1 hi ) ) = hi - lo
--normI :: 𝕀 1 -> Double
--normI ( 𝕀 ( 1 lo ) ( 1 hi ) ) = sqrt $ sup $ ( 𝕀 lo hi ) Ring.^ 2
--normVI :: 𝕀 2 -> Double
--normVI ( 𝕀 ( 2 x_lo y_lo ) ( 2 x_hi y_hi ) ) =
-- sqrt $ sup $ ( 𝕀 x_lo x_hi ) Ring.^ 2 + ( 𝕀 y_lo y_hi ) Ring.^ 2
normVI3 :: 𝕀 3 -> Double
normVI3 ( 𝕀 ( 3 lo x_lo y_lo ) ( 3 hi x_hi y_hi ) ) =
sqrt $ max ( abs lo ) ( abs hi ) Ring.^ 2
+ max ( abs x_lo ) ( abs x_hi ) Ring.^ 2
+ max ( abs y_lo ) ( abs y_hi ) Ring.^ 2
-- | Matrix multiplication \( A v \).
matMulVec
:: forall n m
. ( Representable Double ( n ), Representable Double ( m ) )
=> Vec m ( n ) -- ^ columns of the matrix \( A )
-> 𝕀 m -- ^ vector \( v \)
-> 𝕀 n
matMulVec as v = tabulate $ \ r ->
sum [ scaleInterval ( ( as ! c ) `index` r ) ( index v c )
| c <- toList ( universe @m )
]
{-# INLINEABLE matMulVec #-}
normVI :: ( Applicative ( Vec d ), Representable Double ( d ) ) => 𝕀 d -> Double
normVI ( 𝕀 los his ) =
sqrt $ sum ( nm1 <$> coordinates los <*> coordinates his )
where
nm1 :: Double -> Double -> Double
nm1 lo hi = max ( abs lo ) ( abs hi ) Ring.^ 2
{-# INLINEABLE normVI #-}
-- Use the univariate interval Newton method to narrow from the left
-- a candidate interval.
@ -704,43 +830,36 @@ pipeFunctionsWhileTrue fns = go fns
concat <$> traverse ( go fs ) xs
allNarrowingOperators
:: Double
:: forall n d
. BoxCt n d
=> Double
-> Double
-> ( 𝕀 2 -> D 1 ( 𝕀 2 ) ( 𝕀 3 ) )
-> [ Box -> State Bool [ Box ] ]
-> ( 𝕀 n -> D 1 ( 𝕀 n ) ( 𝕀 d ) )
-> [ Box n -> State Bool [ Box n ] ]
allNarrowingOperators eps_eq eps_bis eqs =
[ \ cand ->
let getter = ( `index` coordIndex )
setter = set coordIndex
newCands = map ( `setter` cand ) $ narrowFn eps_eq eps_bis ( \ x0 -> ff' $ setter x0 cand ) ( getter cand )
newCands = map ( `setter` cand ) $ narrowFn eps_eq eps_bis ( \ x0 -> fn $ setter x0 cand ) ( getter cand )
in do
when ( ( width ( getter cand ) - sum ( map ( width . getter ) newCands ) ) > eps_eq ) $
-- TODO: do a check with the relative reduction in size?
State.put True
return newCands
| narrowFn <- [ leftNarrow, rightNarrow ]
, ( coordIndex, ff' ) <-
[ ( Fin 1, ff'_t i ) | i <- [ Fin 1, Fin 2, Fin 3 ] ]
++
[ ( Fin 2, ff'_s i ) | i <- [ Fin 1, Fin 2, Fin 3 ] ]
, ( coordIndex, fn ) <-
[ ( i, ff' i d )
| i <- toList $ universe @n
, d <- toList $ universe @d
]
]
where
ff'_t i ts =
let D12 { _D12_v = f
, _D12_dx = T f_t
} = eqs ts
in ( f `index` i, f_t `index` i )
ff'_s i ts =
let D12 { _D12_v = f
, _D12_dy = T f_t
} = eqs ts
in ( f `index` i, f_t `index` i )
get_t, get_s :: Box -> 𝕀 Double
get_t = ( `index` ( Fin 1 ) )
get_s = ( `index` ( Fin 2 ) )
set_t, set_s :: 𝕀 Double -> Box -> Box
set_t = set ( Fin 1 )
set_s = set ( Fin 2 )
ff' :: Fin n -> Fin d -> 𝕀 n -> ( 𝕀 Double, 𝕀 Double )
ff' i d ts =
let df = eqs ts
f, f' :: Box d
f = df `monIndex` zeroMonomial
f' = df `monIndex` linearMonomial i
in ( f `index` d, f' `index` d )
{-# INLINEABLE allNarrowingOperators #-}

View file

@ -25,7 +25,7 @@ import Math.Linear
import Math.Module
import Math.Monomial
( multiSubsetSum, multiSubsetsSum
, MonomialBasis ( monTabulate, monIndex )
, MonomialBasisQ ( monTabulateQ, monIndexQ )
)
import Math.Ring
@ -65,7 +65,7 @@ main =
Tasty.testGroup "brush-strokes property tests"
[ Tasty.testGroup "Automatic differentiation"
[ Tasty.testGroup "Monomial basis"
[ testProperty "Round trip D33" testMonomialBasisD33
[ testProperty "Round trip D33" testMonomialBasisQD33
]
, Tasty.testGroup "Monomials"
[ Tasty.testGroup "multiSubsetSum"
@ -172,9 +172,9 @@ testRoundTrip g roundTrip = do
.$ ("value", d )
.$ ("round tripped", roundTrip d )
testMonomialBasisD33 :: Falsify.Property ()
testMonomialBasisD33 =
testRoundTrip genD33 \ d -> $$( monTabulate \ mon -> monIndex [|| d ||] mon )
testMonomialBasisQD33 :: Falsify.Property ()
testMonomialBasisQD33 =
testRoundTrip genD33 \ d -> $$( monTabulateQ \ mon -> monIndexQ [|| d ||] mon )
where
genD33 :: Falsify.Gen ( D3𝔸3 Double )
genD33 =