Slight refactor of bisection dimension choosing logic

brush-strokes/src/cusps/bench/Main.hs

@@ -182,7 +182,7 @@ benchCases :: [ TestCase ]
 benchCases =
   [ ellipseTestCase opts ("minWidth=" ++ show minWidth ++ ",ε=" ++ show narrowAbs) ( 0, 1 ) pi $ defaultStartBoxes [ 2 ]
   | minWidth <- [ 1e-5 ]
-  , narrowAbs <- [ 1e-5, 1e-4, 1e-3, 5e-3, 8e-3, 1e-2, 2e-2, 3e-2, 4e-2, 5e-2, 1e-1 ]
+  , narrowAbs <- [ 5e-2 ]
   , let opts =
           RootIsolationOptions
             { minWidth

brush-strokes/src/lib/Math/Linear.hs

@@ -130,6 +130,7 @@ newtype Vec n a = Vec { vecList :: [ a ] }
 
 universe :: forall n. KnownNat n => Vec n ( Fin n )
 universe = Vec [ Fin i | i <- [ 1 .. fromIntegral ( natVal' @n proxy# ) ] ]
+{-# INLINEABLE universe #-}
 
 coordinates :: forall r u. ( Representable r u ) => u -> Vec ( RepDim u ) r
 coordinates u = fmap ( index u ) $ universe @( RepDim u )
@@ -160,7 +161,7 @@ zipIndices ( Vec v ) = zipIndices_ 1 v
 --------------------------------------------------------------------------------
 
 -- | Rotate a vector by the given angle (counter-clockwise),
--- given the cosine and sine of the angle (in that order)
+-- given the cosine and sine of the angle (in that order).
 rotate :: ( Representable r m, RepDim m ~ 2, Ring r )
        => r -- \( \cos \theta \)
        -> r -- \( \sin \theta \)

brush-strokes/src/lib/Math/Root/Isolation.hs

@@ -44,12 +44,14 @@ import Data.Kind
   ( Type )
 import Data.Foldable
  ( toList )
+import Data.Functor
+  ( (<&>) )
 import Data.List
-  ( partition, sort )
+  ( partition )
 import Data.List.NonEmpty
   ( NonEmpty )
 import qualified Data.List.NonEmpty as NE
-  ( NonEmpty(..), last, singleton )
+  ( NonEmpty(..), cons, filter, fromList, last, singleton, sort )
 import Data.Semigroup
   ( Arg(..), Dual(..) )
 import Numeric
@@ -182,9 +184,17 @@ type BisectionOptions :: Nat -> Nat -> Type
 data BisectionOptions n d =
   BisectionOptions
   { canHaveSols :: !( ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) -> Box n -> Bool )
-  , fallbackBisectionDim :: !( [ ( RootIsolationStep, Box n ) ] -> BoxHistory n -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) ) -> ( Fin n, String ) )
+  , fallbackBisectionDim :: !( BisectionDimPicker n d )
   }
 
+type BisectionDimPicker n d
+  =  forall r
+  .  [ ( RootIsolationStep, Box n ) ]
+  -> BoxHistory n
+  -> ( 𝕀ℝ n -> D 1 ( 𝕀ℝ n ) ( 𝕀ℝ d ) )
+  -> NE.NonEmpty ( Fin n, r )
+  -> ( r, String )
+
 -- | Options for the interval Gauss–Seidel method.
 type GaussSeidelOptions :: Nat -> Nat -> Type
 data GaussSeidelOptions n d =
@@ -192,7 +202,7 @@ data GaussSeidelOptions n d =
     { -- | 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
+      -- (the identity for a well-determined problem, or a projection for
       -- an overdetermined system).
     , gsDims :: ( 𝕀ℝ d -> 𝕀ℝ n ) }
 
@@ -246,7 +256,7 @@ defaultRootIsolationAlgorithms minWidth narrowAbs box history
         -> 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
+      _ -> GaussSeidel _gaussSeidelOptions
               NE.:| [ Box1 _box1Options ]
 
   where
@@ -254,7 +264,7 @@ defaultRootIsolationAlgorithms minWidth narrowAbs box history
     _box1Options =
       Box1Options
         { box1EpsEq = narrowAbs
-        , box1CoordsToNarrow = toList $ universe @n -- [ Fin 1, Fin 2 ]
+        , box1CoordsToNarrow = toList $ universe @n
         , box1EqsToUse = toList $ universe @d
         }
 
@@ -267,23 +277,27 @@ defaultRootIsolationAlgorithms minWidth narrowAbs box history
         , box2EqsToUse = toList $ universe @d
         }
 
+    _gaussSeidelOptions :: GaussSeidelOptions n d
+    _gaussSeidelOptions =
+      GaussSeidelOptions
+        { gsPreconditioner = InverseMidJacobian
+        , gsDims =
+            \ ( 𝕀 ( ℝ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 )
+        }
+
     -- 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 )
+        ( \ cd1 cd2 -> width cd2 - width cd1 >= narrowAbs )
           <$> coordinates b1
           <*> coordinates b2
 {-# INLINEABLE defaultRootIsolationAlgorithms #-}
 
-defaultGaussSeidelOptions :: GaussSeidelOptions N 3
-defaultGaussSeidelOptions =
-  GaussSeidelOptions
-    { gsPreconditioner = InverseMidJacobian
-    , gsDims = \ ( 𝕀 ( ℝ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
   :: forall n d
   .  ( 1 <= n, BoxCt n d )
@@ -295,7 +309,7 @@ defaultBisectionOptions minWidth _narrowAbs box =
         \ eqs box' ->
           -- box(0)-consistency
           let iRange' :: Box d
-              iRange' = ( `monIndex` zeroMonomial ) $ eqs box'
+              iRange' = eqs box' `monIndex` zeroMonomial
           in unT ( origin @Double ) `inside` iRange'
 
           -- box(1)-consistency
@@ -309,52 +323,52 @@ defaultBisectionOptions minWidth _narrowAbs box =
           --    iRange'' = eqs box'' `monIndex` zeroMonomial
           --in unT ( origin @Double ) `inside` iRange''
     , fallbackBisectionDim =
-        \ _roundHist _prevRoundsHist eqs ->
+        \ _roundHist _prevRoundsHist eqs possibleCoordChoices ->
-          let df = eqs box
+          let datPerCoord =
-              datPerCoord =
+                possibleCoordChoices <&> \ ( i, r ) ->
-                [ CoordBisectionData
+                  CoordBisectionData
-                   { coordIndex = i
+                    { coordIndex = i
-                   , coordInterval = box `index` i
+                    , coordInterval = box `index` i
-                   , coordJacobianColumn = df `monIndex` ( linearMonomial i )
+                    , coordJacobianColumn = eqs box `monIndex` ( linearMonomial i )
-                   }
+                    , coordBisectionData = r
-                | i <- toList $ universe @n ]
+                    }
 
-              -- First, check if the largest dimension is over 10 times larger
+              -- First, check if the largest dimension is over 20 times larger
               -- than the smallest dimension; if so bisect along that coordinate.
-          in case sortOnArg ( width . coordInterval ) datPerCoord of
+          in case sortOnArgNE ( width . coordInterval ) datPerCoord of
-                [] -> error "impossible: dimension 0"
+                Arg _ d NE.:| [] -> ( coordBisectionData d, "" )
-                [Arg _ d] -> (coordIndex d, "")
+                Arg w0 _ NE.:| ds ->
-                Arg w0 _ : ds ->
                   let Arg w1 d1 = last ds
-                  in if w1 >= 10 * w0
+                  in if w1 >= 20 * w0
-                     then ( coordIndex d1, "tooWide" )
+                     then ( coordBisectionData d1, "tooWide" )
                      -- Otherwise, pick the dimension with the largest spread = width * Jacobian column norm
-                     else case filter ( not . isTooSmall ) $ sortOnArg ( Dual . spread ) datPerCoord of
+                     else
-                       [] -> ( coordIndex d1, "tooWide'" )
+                      let isTooSmall ( Arg ( Dual w ) _ ) = w < minWidth
-                       Arg _ d : _ -> ( coordIndex d, "spread" )
+                      in case NE.filter ( not . isTooSmall ) $ sortOnArgNE ( Dual . spread ) datPerCoord of
+                        [] -> ( coordBisectionData d1, "tooWide'" )
+                        Arg _ d : _ -> ( coordBisectionData d, "spread" )
     }
-  where
-    isTooSmall ( Arg ( Dual w ) _ ) = w < minWidth
 {-# INLINEABLE defaultBisectionOptions #-}
 
-sortOnArg :: Ord b => (a -> b) -> [a] -> [Arg b a]
+sortOnArgNE :: Ord b => ( a -> b ) -> NE.NonEmpty a -> NE.NonEmpty ( Arg b a )
-sortOnArg f = sort . map ( \ a -> Arg ( f a ) a )
+sortOnArgNE f = NE.sort . fmap ( \ a -> Arg ( f a ) a )
-{-# INLINEABLE sortOnArg #-}
+{-# INLINEABLE sortOnArgNE #-}
 
-type CoordBisectionData :: Nat -> Nat -> Type
+type CoordBisectionData :: Nat -> Nat -> Type -> Type
-data CoordBisectionData n d =
+data CoordBisectionData n d r =
   CoordBisectionData
     { coordIndex :: !( Fin n )
     , coordInterval :: !( 𝕀 Double )
     , coordJacobianColumn :: !( 𝕀ℝ d )
+    , coordBisectionData :: !r
     }
-deriving stock instance Show ( ℝ d )
+deriving stock instance ( Show ( ℝ d ), Show r )
-                     => Show ( CoordBisectionData n d )
+                     => Show ( CoordBisectionData n d r )
 
 spread :: ( BoxCt n d, Representable Double ( ℝ d ) )
-       => CoordBisectionData n d -> Double
+       => CoordBisectionData n d r -> Double
 spread ( CoordBisectionData { coordInterval = cd, coordJacobianColumn = j_cd } )
-  = width cd * normVI j_cd
+  = width cd * normVI j_cd --maxVI j_cd
 {-# INLINEABLE spread #-}
 
 -- | Use the following algorithms using interval arithmetic in order
@@ -400,13 +414,12 @@ isolateRootsIn
       = 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,
+            -- 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
         iRange :: Box d
-        iRange = ( `monIndex` zeroMonomial ) $ eqs cand
+        iRange = eqs cand `monIndex` zeroMonomial
 
     -- Run a round of cusp finding strategies, then recur.
     doStrategies
@@ -477,56 +490,57 @@ doStrategy roundHistory previousRoundsHistory eqs minWidth algo box =
 -- (The difficult part lies in determining along which dimension to bisect.)
 bisect
   :: forall n
-  .  ( KnownNat n, Representable Double ( ℝ n ) )
+  .  ( 1 <= n, KnownNat n, Representable Double ( ℝ n ) )
   => ( Box n -> Bool )
      -- ^ how to check whether a box contains solutions
-  -> ( Fin n, String )
+  -> ( forall r. NE.NonEmpty ( Fin n, r ) -> ( r, String ) )
-      -- ^ fallback choice of dimension (and "why" we chose it)
+      -- ^ heuristic bisection coordinate picker
   -> Box n
+      -- ^ the box to bisect
   -> ( [ Box n ], ( String, Double ) )
-bisect canHaveSols ( fallbackDim, why ) box =
+bisect canHaveSols pickFallBackBisCoord 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.
   case findFewestSols solsList of
-    Just ( Arg _nbSols ( i, mid, oks ) ) ->
+    -- If there is a coordinate for which bisection results in no solutions,
-      ( oks, ("cd = " ++ show i, mid ) )
+    -- or in fewer sub-boxes with solutions than any other coordinate choice,
-    Nothing ->
+    -- pick that coordinate for bisection.
-      case bisectInCoord box fallbackDim of
+    ( i, ( mid, subBoxesWithSols ) ) NE.:| [] ->
-        ( mid, ( lo, hi ) ) ->
+      ( subBoxesWithSols, ( "cd = " ++ show i, mid ) )
-          ( [ lo, hi ], ( why, mid ) )
+    -- Otherwise, fall back to the provided heuristic.
+    is ->
+      let ( ( mid, subBoxesWithSols ), why ) = pickFallBackBisCoord is
+      in ( subBoxesWithSols, ( why, mid ) )
   where
     solsList =
-      [ Arg ( fromIntegral $ length oks ) ( i, mid, oks )
+      NE.fromList -- (n >= 1 so definitely non-empty)
-      | i <- toList $ universe @n
+        [ Arg ( fromIntegral $ length subBoxesWithSols ) ( i, ( mid, subBoxesWithSols ) )
-      , let (mid, (lo, hi)) = bisectInCoord box i
+        | i <- toList $ universe @n
-            lo_ok = canHaveSols lo
+        , let ( mid, ( loBox, hiBox ) ) = bisectInCoord box i
-            hi_ok = canHaveSols hi
+              loBox_ok = canHaveSols loBox
-            oks = [ lo | lo_ok ] ++ [ hi | hi_ok ]
+              hiBox_ok = canHaveSols hiBox
-      ]
+              subBoxesWithSols = [ loBox | loBox_ok ] ++ [ hiBox | hiBox_ok ]
+        ]
 {-# INLINEABLE bisect #-}
 
--- | Find any element with the least argument.
+-- | Return the elements 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 :: forall a. NE.NonEmpty ( Arg Word a ) -> NE.NonEmpty a
-findFewestSols [] = Nothing
+findFewestSols ( ( Arg nbSols arg ) NE.:| args )
-findFewestSols ( arg@( Arg nbSols _ ) : args )
   | nbSols == 0
-  = Just arg
+  = NE.singleton arg
   | otherwise
-  = Just $ go arg args
+  = go nbSols ( NE.singleton arg ) args
   where
-    go :: Arg Word a -> [ Arg Word a ] -> Arg Word a
+    go :: Word -> NE.NonEmpty a -> [ Arg Word a ] -> NE.NonEmpty a
-    go bestSoFar [] = bestSoFar
+    go _               bestSoFar [] = bestSoFar
-    go bestSoFar@( Arg bestNbSolsSoFar _ ) ( arg'@( Arg nbSols' _ ) : args' )
+    go bestNbSolsSoFar bestSoFar ( ( Arg nbSols' arg' ) : args' )
       | nbSols' == 0
-      = arg'
+      = NE.singleton arg'
-      | nbSols' < bestNbSolsSoFar
-      = go arg' args'
       | otherwise
-      = go bestSoFar args'
+      = case compare nbSols' bestNbSolsSoFar of
+          LT -> go nbSols' ( NE.singleton arg' ) args'
+          GT -> go bestNbSolsSoFar bestSoFar args'
+          EQ -> go bestNbSolsSoFar ( arg `NE.cons` bestSoFar ) args'
 
 bisectInCoord
   :: Representable Double ( ℝ n )
@@ -755,8 +769,8 @@ matMulVec
   -> 𝕀ℝ m          -- ^ vector \( v \)
   -> 𝕀ℝ n
 matMulVec as v = tabulate $ \ r ->
-  sum [ scaleInterval ( ( as ! c ) `index` r ) ( index v c )
+  sum [ scaleInterval ( a `index` r ) ( index v c )
-      | c <- toList ( universe @m )
+      | ( c, a ) <- toList ( (,) <$> universe @m <*> as )
       ]
 {-# INLINEABLE matMulVec #-}
 
@@ -768,6 +782,14 @@ normVI ( 𝕀 los his ) =
       nm1 lo hi = max ( abs lo ) ( abs hi ) Ring.^ 2
 {-# INLINEABLE normVI #-}
 
+maxVI :: ( Applicative ( Vec d ), Representable Double ( ℝ d ) ) => 𝕀ℝ d -> Double
+maxVI ( 𝕀 los his ) =
+  maximum ( maxAbs <$> coordinates los <*> coordinates his )
+    where
+      maxAbs :: Double -> Double -> Double
+      maxAbs lo hi = max ( abs lo ) ( abs hi )
+{-# INLINEABLE maxVI  #-}
+
 -- Use the univariate interval Newton method to narrow from the left
 -- a candidate interval.
 --