Use 'n choose k' to choose Gauss-Seidel dimensions

This commit is contained in:
sheaf 2024-04-21 21:15:07 +02:00
parent 0e6b9a822b
commit 0e59d85143
4 changed files with 61 additions and 20 deletions

View file

@ -185,8 +185,7 @@ benchCases =
, narrowAbs <- [ 5e-2 ]
, let opts =
RootIsolationOptions
{ minWidth
, cuspFindingAlgorithms = defaultRootIsolationAlgorithms minWidth narrowAbs
{ rootIsolationAlgorithms = defaultRootIsolationAlgorithms minWidth narrowAbs
}
]

View file

@ -229,6 +229,8 @@ data OutlineInfo =
{ outlineFn :: OutlineFn
, outlineDefiniteCusps, outlinePotentialCusps :: [ Cusp ] }
type N = 2
computeStrokeOutline ::
forall ( clo :: SplineType ) usedParams brushParams crvData ptData s
. ( KnownSplineType clo

View file

@ -10,7 +10,8 @@ module Math.Linear
-- * Points and vectors (second version)
, (..), T(.., V2, V3, V4)
, Fin(..), MFin(..), universe, coordinates
, Fin(..), MFin(..)
, universe, coordinates, choose
, RepDim, RepresentableQ(..)
, Representable(..), set, injection, projection
, Vec(..), (!), find, zipIndices
@ -34,7 +35,9 @@ import GHC.Generics
import GHC.Stack
( HasCallStack )
import GHC.TypeNats
( Nat, KnownNat, natVal' )
( Nat, KnownNat, type (<=)
, natVal'
)
-- acts
import Data.Act
@ -136,6 +139,27 @@ coordinates :: forall r u. ( Representable r u ) => u -> Vec ( RepDim u ) r
coordinates u = fmap ( index u ) $ universe @( RepDim u )
{-# INLINEABLE coordinates #-}
-- | Binomial coefficient: choose all subsets of size @k@ of the given set
-- of size @n@.
choose
:: forall n k
. ( KnownNat n, KnownNat k
, 1 <= n, 1 <= k, k <= n
)
=> [ Vec k ( Fin n ) ]
choose = coerce $ go ( fromIntegral $ natVal' @n proxy# )
( fromIntegral $ natVal' @k proxy# )
where
go :: Word -> Word -> [ [ Word ] ]
go n k
| k == 1
= [ [ i ] | i <- [ 1 .. n ] ]
| n == k
= [ [ 1 .. n ] ]
go n k = go ( n - 1 ) k
++ ( map ( ++ [ n ] ) $ go ( n - 1 ) ( k - 1 ) )
{-# INLINEABLE choose #-}
infixl 9 !
(!) :: forall l a. HasCallStack => Vec l a -> Fin l -> a
( Vec l ) ! Fin i = l !! ( fromIntegral i - 1 )

View file

@ -24,10 +24,6 @@ module Math.Root.Isolation
-- ** Trees recording search space of root isolation algorithms
, RootIsolationTree(..), showRootIsolationTree
, RootIsolationStep(..)
-- * Hack for changing between 2 and 3 d formulations
-- for my personal testing
, N
)
where
@ -52,13 +48,19 @@ import Data.List.NonEmpty
( NonEmpty )
import qualified Data.List.NonEmpty as NE
( NonEmpty(..), cons, filter, fromList, last, singleton, sort )
import Data.Proxy
( Proxy(..) )
import Data.Semigroup
( Arg(..), Dual(..) )
import Data.Type.Ord
( OrderingI(..) )
import Numeric
( showFFloat )
import GHC.TypeNats
( Nat, KnownNat
, type (<=) )
, type (<=)
, cmpNat
)
-- containers
import Data.Tree
@ -88,7 +90,7 @@ import Math.Linear
import Math.Module
( Module(..) )
import Math.Monomial
( MonomialBasis(..)
( MonomialBasis(..), Deg, Vars
, linearMonomial, zeroMonomial
)
import qualified Math.Ring as Ring
@ -137,10 +139,18 @@ showArea area = "(area " ++ showFFloat (Just 6) area "" ++ ")"
type Box n = 𝕀 n
type BoxHistory n = [ NE.NonEmpty ( RootIsolationStep, Box n ) ]
type N = 2
type BoxCt n d = ( n ~ N, d ~ 3 )
{-
( Show ( 𝕀 n ), Show ( n )
-- | Dimension constraints for root isolation in a system of equations:
--
-- - @n@: number of variables
-- - @d@: number of equations
--
-- NB: we require n <= d (no support for under-constrained systems).
type BoxCt n d =
( KnownNat n, KnownNat d
, 1 <= n, 1 <= d, n <= d
, Show ( 𝕀 n ), Show ( n )
, Eq ( n )
, Representable Double ( n )
, MonomialBasis ( D 1 ( n ) )
@ -153,7 +163,6 @@ type BoxCt n d = ( n ~ N, d ~ 3 )
, Module Double ( T ( d ) )
, Representable Double ( d )
)
-}
-- | Options for the root isolation methods in 'isolateRootsIn'.
type RootIsolationOptions :: Nat -> Nat -> Type
@ -235,6 +244,7 @@ defaultRootIsolationOptions =
where
minWidth = 1e-5
ε_eq = 5e-3
{-# INLINEABLE defaultRootIsolationOptions #-}
defaultRootIsolationAlgorithms
:: forall n d
@ -284,11 +294,17 @@ defaultRootIsolationAlgorithms minWidth ε_eq box history
GaussSeidelOptions
{ gsPreconditioner = InverseMidJacobian
, gsPickEqs =
\ ( 𝕀 ( 3 a_lo b_lo c_lo ) ( 3 a_hi b_hi c_hi ) ) ->
case length history `mod` 3 of
0 -> 𝕀 ( 2 a_lo b_lo ) ( 2 a_hi b_hi )
1 -> 𝕀 ( 2 b_lo c_lo ) ( 2 b_hi c_hi )
_ -> 𝕀 ( 2 a_lo c_lo ) ( 2 a_hi c_hi )
case cmpNat @n @d Proxy Proxy of
EQI -> id
LTI ->
-- If there are more equations (d) than variables (n),
-- pick a size n subset of the variables,
-- (go through all combinations cyclically).
let choices :: [ Vec n ( Fin d ) ]
choices = choose @d @n
choice :: Vec n ( Fin d )
choice = choices !! ( length history `mod` length choices )
in \ u -> tabulate \ i -> index u ( choice ! i )
}
-- Did we reduce the box width by at least ε_eq