From 0160081e80bd957ad0f09cd87979c58c61c0a491 Mon Sep 17 00:00:00 2001 From: sheaf Date: Thu, 18 Apr 2024 20:14:19 +0200 Subject: [PATCH] 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. --- brush-strokes/src/cusps/bench/Main.hs | 14 +- brush-strokes/src/lib/Math/Algebra/Dual.hs | 48 +- .../src/lib/Math/Algebra/Dual/Internal.hs | 600 +++++++------- brush-strokes/src/lib/Math/Bezier/Stroke.hs | 62 +- brush-strokes/src/lib/Math/Interval.hs | 48 +- brush-strokes/src/lib/Math/Linear.hs | 67 +- brush-strokes/src/lib/Math/Linear/Internal.hs | 13 +- brush-strokes/src/lib/Math/Monomial.hs | 88 +- brush-strokes/src/lib/Math/Root/Isolation.hs | 763 ++++++++++-------- brush-strokes/src/test/Main.hs | 10 +- 10 files changed, 954 insertions(+), 759 deletions(-) diff --git a/brush-strokes/src/cusps/bench/Main.hs b/brush-strokes/src/cusps/bench/Main.hs index 3fa5fcd..7bd4867 100644 --- a/brush-strokes/src/cusps/bench/Main.hs +++ b/brush-strokes/src/cusps/bench/Main.hs @@ -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 diff --git a/brush-strokes/src/lib/Math/Algebra/Dual.hs b/brush-strokes/src/lib/Math/Algebra/Dual.hs index 4f73d74..d32f2cd 100644 --- a/brush-strokes/src/lib/Math/Algebra/Dual.hs +++ b/brush-strokes/src/lib/Math/Algebra/Dual.hs @@ -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 diff --git a/brush-strokes/src/lib/Math/Algebra/Dual/Internal.hs b/brush-strokes/src/lib/Math/Algebra/Dual/Internal.hs index 2d9c446..787c03f 100644 --- a/brush-strokes/src/lib/Math/Algebra/Dual/Internal.hs +++ b/brush-strokes/src/lib/Math/Algebra/Dual/Internal.hs @@ -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 ||] diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke.hs b/brush-strokes/src/lib/Math/Bezier/Stroke.hs index 93b3140..60e6e64 100644 --- a/brush-strokes/src/lib/Math/Bezier/Stroke.hs +++ b/brush-strokes/src/lib/Math/Bezier/Stroke.hs @@ -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 ) ) + + +-} \ No newline at end of file diff --git a/brush-strokes/src/lib/Math/Interval.hs b/brush-strokes/src/lib/Math/Interval.hs index 1ca0582..0ddbe84 100644 --- a/brush-strokes/src/lib/Math/Interval.hs +++ b/brush-strokes/src/lib/Math/Interval.hs @@ -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 diff --git a/brush-strokes/src/lib/Math/Linear.hs b/brush-strokes/src/lib/Math/Linear.hs index d49c687..c5f0554 100644 --- a/brush-strokes/src/lib/Math/Linear.hs +++ b/brush-strokes/src/lib/Math/Linear.hs @@ -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 diff --git a/brush-strokes/src/lib/Math/Linear/Internal.hs b/brush-strokes/src/lib/Math/Linear/Internal.hs index f0b1bbf..d01d724 100644 --- a/brush-strokes/src/lib/Math/Linear/Internal.hs +++ b/brush-strokes/src/lib/Math/Linear/Internal.hs @@ -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 ) ) diff --git a/brush-strokes/src/lib/Math/Monomial.hs b/brush-strokes/src/lib/Math/Monomial.hs index ee2be74..2356202 100644 --- a/brush-strokes/src/lib/Math/Monomial.hs +++ b/brush-strokes/src/lib/Math/Monomial.hs @@ -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 ] ) ||] diff --git a/brush-strokes/src/lib/Math/Root/Isolation.hs b/brush-strokes/src/lib/Math/Root/Isolation.hs index c603b08..5d18c44 100644 --- a/brush-strokes/src/lib/Math/Root/Isolation.hs +++ b/brush-strokes/src/lib/Math/Root/Isolation.hs @@ -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 ) -- | Gauss–Seidel 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 Gauss–Seidel 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 Gauss–Seidel 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 Gauss–Seidel 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 Gauss–Seidel step. intervalGaussSeidel - :: ( 𝕀ℝ 2 -> D 1 ( 𝕀ℝ 2 ) ( 𝕀ℝ 3 ) ) - -- ^ equations - -> Preconditioner - -- ^ preconditioner to use for the Gauss–Seidel 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 Gauss–Seidel 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 Gauss–Seidel 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 Gauss–Seidel step. + gsGuesses = map ( first ( \ box' -> unT $ box' ^+^ T boxMid ) ) + $ gaussSeidelStep a b ( T box ^-^ T boxMid ) in -- If the Gauss–Seidel 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 Gauss–Seidel 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 Gauss–Seidel 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 Gauss–Seidel 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 Gauss–Seidel 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 #-} diff --git a/brush-strokes/src/test/Main.hs b/brush-strokes/src/test/Main.hs index b85fc22..734ef35 100644 --- a/brush-strokes/src/test/Main.hs +++ b/brush-strokes/src/test/Main.hs @@ -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 =