brush-strokes-0.1.0.0: src/lib/Math/Root/Isolation/Degree.hs

{-# LANGUAGE ScopedTypeVariables #-}

module Math.Root.Isolation.Degree
  ( -- * Computing the topological degree
    topologicalDegree
  )
  where

-- base
import Data.Kind
  ( Type )
import GHC.TypeNats
  ( Nat )

-- primitive
import Data.Primitive.SmallArray

-- brush-strokes
import Math.Interval
import Math.Linear

--------------------------------------------------------------------------------
-- Topological degree.

-- | Compute the topological degree of the given interval function (2D only).
topologicalDegree
  :: Double
     -- ^ minimum width when bisecting facets
  -> ( 𝕀ℝ 2 -> 𝕀ℝ 2 )
     -- ^ function
  -> 𝕀ℝ 2
     -- ^ box
  -> Maybe Int
topologicalDegree ε_bis f ( 𝕀ℝ2 ( 𝕀 x_lo x_hi ) ( 𝕀 y_lo y_hi ) ) =
  topologicalDegreeFromContour $
    smallArrayFromList $
      concatMap ( tagFacet ε_bis f )
        [ ( Tag ( Fin 2 ) P, 𝕀ℝ2 ( 𝕀 x_hi x_hi ) ( 𝕀 y_lo y_hi ) )
        , ( Tag ( Fin 1 ) N, 𝕀ℝ2 ( 𝕀 x_lo x_hi ) ( 𝕀 y_hi y_hi ) )
        , ( Tag ( Fin 2 ) N, 𝕀ℝ2 ( 𝕀 x_lo x_lo ) ( 𝕀 y_lo y_hi ) )
        , ( Tag ( Fin 1 ) P, 𝕀ℝ2 ( 𝕀 x_lo x_hi ) ( 𝕀 y_lo y_lo ) ) ]

-- | Tag a facet with the sign of one of the equations that is non-zero
-- on that facet.
--
-- If all equations take the value 0 on that facet, bisect and recur.
tagFacet
  :: Double
  -> ( 𝕀ℝ 2 -> 𝕀ℝ 2 )
  -> ( Tag 2, 𝕀ℝ 2 )
  -> [ Tag 2 ]
tagFacet ε_bis f ( tag@( Tag i _ ), facet0 ) = go facet0
  where
    go facet
      | f1_b_lo > 0
      = [ Tag ( Fin 1 ) P ]
      | f1_b_hi < 0
      = [ Tag ( Fin 1 ) N ]
      | f2_b_lo > 0
      = [ Tag ( Fin 2 ) P ]
      | f2_b_hi < 0
      = [ Tag ( Fin 2 ) N ]
      | width ( facet `index` i ) < ε_bis
      = [ NoTag ] -- TODO: shortcut the whole computation;
                  -- if we fail to tag one segment, the final answer could
                  -- be off by { -2, -1, 0, 1, 2 } (I believe)
      | otherwise
      = concatMap go $ bisectFacet tag facet
      where
        𝕀 f1_b_lo f1_b_hi = f facet `index` Fin 1
        𝕀 f2_b_lo f2_b_hi = f facet `index` Fin 2
tagFacet _ _ ( NoTag, _ ) = error "impossible: input always tagged"

bisectFacet :: Tag 2 -> 𝕀ℝ 2 -> [ 𝕀ℝ 2 ]
bisectFacet ( Tag i s ) x =
  ( if s == P then id else reverse ) $ bisectInCoord i x
bisectFacet NoTag _ = error "impossible: input always tagged"

bisectInCoord :: Fin n -> 𝕀ℝ 2 -> [ 𝕀ℝ 2 ]
bisectInCoord ( Fin 1 ) ( 𝕀ℝ2 ( 𝕀 x_lo x_hi ) ( 𝕀 y_lo y_hi ) ) =
  [ 𝕀ℝ2 ( 𝕀 x_lo  x_mid ) ( 𝕀 y_lo  y_hi )
  , 𝕀ℝ2 ( 𝕀 x_mid x_hi  ) ( 𝕀 y_lo y_hi ) ]
    where x_mid = 0.5 * ( x_lo + x_hi )
bisectInCoord _ ( 𝕀ℝ2 ( 𝕀 x_lo x_hi ) ( 𝕀 y_lo y_hi ) ) =
  [ 𝕀ℝ2 ( 𝕀 x_lo x_hi  ) ( 𝕀 y_lo  y_mid )
  , 𝕀ℝ2 ( 𝕀 x_lo  x_hi ) ( 𝕀 y_mid  y_hi ) ]
    where y_mid = 0.5 * ( y_lo + y_hi )

-- | Compute the topological degree of \( f \colon \mathbb{IR}^n \to \mathbb{IR}^n \)
-- on a box \( D \subset \mathbb{IR}^n \) by using the signs of the components
-- of \( f \) on a counter-clockwise polygonal contour describing the boundary of \( D \).
topologicalDegreeFromContour :: SmallArray ( Tag 2 ) -> Maybe Int
topologicalDegreeFromContour facetTags = go 0 0
  where
    n = sizeofSmallArray facetTags
    go :: Int -> Int -> Maybe Int
    go !d !i
      | i >= n
      = Just d
      | otherwise
      = case indexSmallArray facetTags i of
          NoTag -> Nothing
          Tag ( Fin 1 ) P -> go ( d + δ ) ( i + 1 )
            where
              δ = if indexSmallArray facetTags ( if i == n - 1 then 0 else i + 1 ) == Tag ( Fin 2 ) P
                  then 1 else 0
                - if indexSmallArray facetTags ( if i == 0 then n - 1 else i - 1 ) == Tag ( Fin 2 ) P
                  then 1 else 0
          _ -> go d ( i + 1 )

--------------------------------------------------------------------------------
-- Tagging edges with signs of certain components of the function.

data Sign = P | N
  deriving stock ( Eq, Ord )
instance Show Sign where
  showsPrec _ = \case { P -> showString "+"; N -> showString "-" }
type Tag :: Nat -> Type
newtype Tag d = MkTag Int
  deriving newtype ( Eq, Ord )
pattern NoTag :: Tag d
pattern NoTag = MkTag 0
pattern Tag :: Fin d -> Sign -> Tag d
pattern Tag c s <- ( getTag -> (# c, s #) )
  where Tag ( Fin c ) s = MkTag $ case s of { P -> fromIntegral c; N -> -( fromIntegral c ) }
{-# COMPLETE Tag, NoTag #-}
getTag :: Tag d -> (# Fin d, Sign #)
getTag ( MkTag t ) = (# Fin ( fromIntegral $ abs t ), if t >= 0 then P else N #)