brush-strokes-0.1.0.0: src/lib/Math/Interval.hs

{-# LANGUAGE CPP #-}

{-# LANGUAGE ScopedTypeVariables  #-}
{-# LANGUAGE TemplateHaskell      #-}
{-# LANGUAGE UndecidableInstances #-}

{-# OPTIONS_GHC -Wno-orphans #-}

module Math.Interval
  ( 𝕀(𝕀), inf, sup
  , midpoint , scaleInterval, width

  , hull
  , addWithErr, subWithErr, prodWithErr, divWithErr

  , isCanonical
  , (∩), (⊘), (∈), (∈∈), (⊖)
  , extendedRecip

  , singleton, point
  , nonDecreasing
  , aabb, bisect

  , 𝕀ℝ(..)
  )
  where

-- base
import Prelude hiding ( Num(..), Fractional(..) )
import qualified Data.List.NonEmpty as NE

-- acts
import Data.Act
  ( Act((•)), Torsor((-->)) )

-- groups
import Data.Group
  ( Group(..) )

-- groups-generic
import Data.Group.Generics
  ( )

-- brush-strokes
import Math.Algebra.Dual
import Math.Float.Utils
  ( prevFP )
import Math.Interval.Internal
  ( 𝕀(𝕀), inf, sup
  , scaleInterval, width, singleton
  , hull
  , addWithErr, subWithErr, prodWithErr
  , 𝕀ℝ(..), divWithErr
  )
import Math.Linear
  ( ℝ(..), T(..)
  , RepresentableQ(..), Representable(..), coordinates
  )
import Math.Module
import Math.Monomial
import Math.Ring

--------------------------------------------------------------------------------
-- Interval arithmetic.

type instance Dim ( 𝕀ℝ n ) = n

-- | Turn a non-decreasing function into a function on intervals.
nonDecreasing :: forall n m
              .  ( Representable Double ( ℝ n )
                 , Representable 𝕀 ( 𝕀ℝ n )
                 , Representable Double ( ℝ m )
                 , Representable 𝕀 ( 𝕀ℝ m ) )
              => ( ℝ n -> ℝ m ) -> 𝕀ℝ n -> 𝕀ℝ m
nonDecreasing f u =
  let
    lo, hi :: ℝ n
    lo = tabulate \ i -> inf ( u `index` i )
    hi = tabulate \ i -> sup ( u `index` i )
    lo', hi' :: ℝ m
    lo' = f lo
    hi' = f hi
  in tabulate \ j -> 𝕀 ( lo' `index` j ) ( hi' `index` j )
{-# INLINEABLE nonDecreasing #-}

-- | Does the given value lie inside the specified interval?
(∈) :: Double -> 𝕀 -> Bool
x ∈ 𝕀 lo hi = x >= lo && x <= hi
infix 4 ∈

-- | Does the given point lie inside the specified box?
(∈∈) :: ( Representable Double ( ℝ n ), Representable 𝕀 ( 𝕀ℝ n ) ) => ℝ n -> 𝕀ℝ n -> Bool
p ∈∈ i = and $ (∈) <$> coordinates p <*> coordinates i
infix 4 ∈∈
{-# INLINEABLE (∈∈) #-}
{-# SPECIALISE INLINE (∈∈) @1 #-}
{-# SPECIALISE INLINE (∈∈) @2 #-}
{-# SPECIALISE INLINE (∈∈) @3 #-}
{-# SPECIALISE INLINE (∈∈) @4 #-}

-- | Is this interval canonical, i.e. it consists of either 1 or 2 floating point
-- values only?
isCanonical :: 𝕀 -> Bool
isCanonical ( 𝕀 x_inf x_sup ) = x_inf >= prevFP x_sup

-- | The midpoint of an interval.
--
-- NB: assumes the endpoints are neither @NaN@ nor infinity.
midpoint :: 𝕀 -> Double
midpoint ( 𝕀 x_inf x_sup ) = 0.5 * ( x_inf + x_sup )

point :: ( Representable Double ( ℝ n ), Representable 𝕀 ( 𝕀ℝ n ) )
      => ℝ n -> 𝕀ℝ n
point x = tabulate \ i -> singleton ( x `index` i )
{-# INLINEABLE point #-}

-- | Compute the intersection of two intervals.
--
-- Returns whether the first interval is a strict subset of the second interval
-- (or the intersection is a single point).
(∩) :: 𝕀 -> 𝕀 -> [ ( 𝕀, Bool ) ]
𝕀 lo1 hi1 ∩ 𝕀 lo2 hi2
  | lo > hi
  = [ ]
  | otherwise
  = [ ( 𝕀 lo hi, ( lo1 > lo2 && hi1 < hi2 ) || ( lo == hi ) ) ]
  where
    lo = max lo1 lo2
    hi = min hi1 hi2
infix 3 ∩

infixl 6 ⊖
(⊖) :: 𝕀 -> 𝕀 -> 𝕀
(⊖) a@( 𝕀 lo1 hi1 ) b@( 𝕀 lo2 hi2 )
  | width a >= width b
  = 𝕀 ( lo1 - lo2 ) ( hi1 - hi2 )
  | otherwise
  = 𝕀 ( hi1 - hi2 ) ( lo1 - lo2 )

-- | Bisect an interval.
--
-- Generally returns two sub-intervals, but returns the input interval if
-- it is canonical (and thus bisection would be pointless).
bisect :: 𝕀 -> NE.NonEmpty 𝕀
bisect x@( 𝕀 x_inf x_sup )
  | isCanonical x
  = NE.singleton x
  | otherwise
  = 𝕀 x_inf x_mid NE.:| [ 𝕀 x_mid x_sup ]
  where x_mid = midpoint x

deriving via ViaAbelianGroup ( T 𝕀 )
  instance Semigroup ( T 𝕀 )
deriving via ViaAbelianGroup ( T 𝕀 )
  instance Monoid ( T 𝕀 )
deriving via ViaAbelianGroup ( T 𝕀 )
  instance Group ( T 𝕀 )

instance Act   ( T 𝕀 ) 𝕀 where
  T g • a = g + a
instance Torsor ( T 𝕀 ) 𝕀 where
  a --> b = T $ b - a

--------------------------------------------------------------------------------
-- Extended division

-- | Extended division, returning either 1 or 2 intervals.
--
-- NB: this function can return overlapping intervals if both arguments contain 0.
-- Otherwise, the returned intervals are guaranteed to be disjoint.
(⊘) :: 𝕀 -> 𝕀 -> [ 𝕀 ]
x ⊘ y
#ifdef ASSERTS
  | 0 ∈ x && 0 ∈ y
  = error $ unlines
      [ "x ⊘ y: both arguments contain zero"
      , "x: " ++ show x, "y: " ++ show y ]
  | otherwise
#endif
  = map ( x * ) ( extendedRecip y )
infixl 7 ⊘

-- | Extended reciprocal, returning either 1 or 2 intervals.
extendedRecip :: 𝕀 -> [ 𝕀 ]
extendedRecip x@( 𝕀 lo hi )
  | lo == 0 && hi == 0
  = [ 𝕀 negInf posInf ]
  | lo >= 0 || hi <= 0
  = [ recip x ]
  | otherwise
  = [ 𝕀 negInf ( recip lo ), 𝕀 ( recip hi ) posInf ]
  where

negInf, posInf :: Double
negInf = -1 / 0
posInf =  1 / 0

-------------------------------------------------------------------------------
-- Lattices.

aabb :: ( Representable 𝕀 ( 𝕀ℝ n ), Functor f )
     => f ( 𝕀ℝ n ) -> ( f 𝕀 -> 𝕀 ) -> 𝕀ℝ n
aabb fv extrema = tabulate \ i -> extrema ( fmap ( `index` i ) fv )
{-# INLINEABLE aabb #-}

--------------------------------------------------------------------------------

instance Module 𝕀 ( T ( 𝕀ℝ 0 ) ) where
  origin  = T 𝕀ℝ0
  _ ^+^ _ = T 𝕀ℝ0
  _ ^-^ _ = T 𝕀ℝ0
  _ *^  _ = T 𝕀ℝ0

instance Module 𝕀 ( T ( 𝕀ℝ 1 ) ) where
  origin      = T $ 𝕀ℝ1 ( unT origin )
  T ( 𝕀ℝ1 x1 ) ^+^ T ( 𝕀ℝ1 x2 ) = T $ 𝕀ℝ1 ( x1 + x2 )
  T ( 𝕀ℝ1 x1 ) ^-^ T ( 𝕀ℝ1 x2 ) = T $ 𝕀ℝ1 ( x1 - x2 )
  k *^ T ( 𝕀ℝ1 x1 ) = T $ 𝕀ℝ1 ( k * x1 )

instance Module 𝕀 ( T ( 𝕀ℝ 2 ) ) where
  origin      = T $ 𝕀ℝ2 ( unT origin ) ( unT origin )
  T ( 𝕀ℝ2 x1 y1 ) ^+^ T ( 𝕀ℝ2 x2 y2 ) = T $ 𝕀ℝ2 ( x1 + x2 ) ( y1 + y2 )
  T ( 𝕀ℝ2 x1 y1 ) ^-^ T ( 𝕀ℝ2 x2 y2 ) = T $ 𝕀ℝ2 ( x1 - x2 ) ( y1 - y2 )
  k *^ T ( 𝕀ℝ2 x1 y1 ) = T $ 𝕀ℝ2 ( k * x1 ) ( k * y1 )

instance Module 𝕀 ( T ( 𝕀ℝ 3 ) ) where
  origin      = T $ 𝕀ℝ3 ( unT origin ) ( unT origin ) ( unT origin )
  T ( 𝕀ℝ3 x1 y1 z1 ) ^+^ T ( 𝕀ℝ3 x2 y2 z2 ) = T $ 𝕀ℝ3 ( x1 + x2 ) ( y1 + y2 ) ( z1 + z2 )
  T ( 𝕀ℝ3 x1 y1 z1 ) ^-^ T ( 𝕀ℝ3 x2 y2 z2 ) = T $ 𝕀ℝ3 ( x1 - x2 ) ( y1 - y2 ) ( z1 - z2 )
  k *^ T ( 𝕀ℝ3 x1 y1 z1 ) = T $ 𝕀ℝ3 ( k * x1 ) ( k * y1 ) ( k * z1 )

instance Module 𝕀 ( T ( 𝕀ℝ 4 ) ) where
  origin      = T $ 𝕀ℝ4 ( unT origin ) ( unT origin ) ( unT origin ) ( unT origin )
  T ( 𝕀ℝ4 x1 y1 z1 w1 ) ^+^ T ( 𝕀ℝ4 x2 y2 z2 w2 ) = T $ 𝕀ℝ4 ( x1 + x2 ) ( y1 + y2 ) ( z1 + z2 ) ( w1 + w2 )
  T ( 𝕀ℝ4 x1 y1 z1 w1 ) ^-^ T ( 𝕀ℝ4 x2 y2 z2 w2 ) = T $ 𝕀ℝ4 ( x1 - x2 ) ( y1 - y2 ) ( z1 - z2 ) ( w1 - w2 )
  k *^ T ( 𝕀ℝ4 x1 y1 z1 w1 ) = T $ 𝕀ℝ4 ( k * x1 ) ( k * y1 ) ( k * z1 ) ( k * w1 )

instance Inner  𝕀 ( T ( 𝕀ℝ 2 ) ) where
  T ( 𝕀ℝ2 x1 y1 ) ^.^ T ( 𝕀ℝ2 x2 y2 ) = x1 * x2 + y1 * y2

instance Cross  𝕀 ( T ( 𝕀ℝ 2 ) ) where
  T ( 𝕀ℝ2 x1 y1 ) × T ( 𝕀ℝ2 x2 y2 ) = x1 * y2 - x2 * y1

deriving via ViaModule 𝕀 ( T ( 𝕀ℝ n ) )
  instance Module 𝕀 ( T ( 𝕀ℝ n ) ) => Semigroup ( T ( 𝕀ℝ n ) )
deriving via ViaModule 𝕀 ( T ( 𝕀ℝ n ) )
  instance Module 𝕀 ( T ( 𝕀ℝ n ) ) => Monoid ( T ( 𝕀ℝ n ) )
deriving via ViaModule 𝕀 ( T ( 𝕀ℝ n ) )
  instance Module 𝕀 ( T ( 𝕀ℝ n ) ) => Group ( T ( 𝕀ℝ n ) )

deriving via ViaModule 𝕀 ( 𝕀ℝ n )
  instance Module 𝕀 ( T ( 𝕀ℝ n ) ) => Act ( T ( 𝕀ℝ n ) ) ( 𝕀ℝ n )
deriving via ( ViaModule 𝕀 ( 𝕀ℝ n ) )
  instance Module 𝕀 ( T ( 𝕀ℝ n ) ) => Torsor ( T ( 𝕀ℝ n ) ) ( 𝕀ℝ n )

--------------------------------------------------------------------------------
-- HasChainRule instances.

instance HasChainRule 𝕀 2 ( 𝕀ℝ 0 ) where
  linearD f v = D0 ( f v )
  chain _ ( D0 gfx ) = D21 gfx origin origin

instance HasChainRule 𝕀 3 ( 𝕀ℝ 0 ) where
  linearD f v = D0 ( f v )
  chain _ ( D0 gfx ) = D31 gfx origin origin origin

instance HasChainRule 𝕀 2 ( 𝕀ℝ 1 ) where

  linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 1 -> w ) -> 𝕀ℝ 1 -> D2𝔸1 w
  linearD f v =
    let !o = origin @𝕀 @( T w )
    in $$( monTabulateQ \ mon ->
              if | isZeroMonomial mon
                 -> [|| f v ||]
                 | Just i <- isLinear mon
                 -> [|| f $$( tabulateQ \ j ->
                                if | j == i
                                   -> [|| 𝕀 1 1 ||]
                                   | otherwise
                                   -> [|| 𝕀 0 0 ||]
                            ) ||]
                 | otherwise
                 -> [|| unT o ||]
         )

  chain :: forall w. Module 𝕀 ( T w ) => D2𝔸1 ( 𝕀ℝ 1 ) -> D2𝔸1 w -> D2𝔸1 w
  chain !df !dg =
    let !o = origin @𝕀 @( T w )
        !p = (^+^)  @𝕀 @( T w )
        !s = (^*)   @𝕀 @( T w )
    in $$( chainRule1NQ
         [|| o ||] [|| p ||] [|| s ||]
         [|| df ||] [|| dg ||] )

instance HasChainRule 𝕀 3 ( 𝕀ℝ 1 ) where

  linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 1 -> w ) -> 𝕀ℝ 1 -> D3𝔸1 w
  linearD f v =
    let !o = origin @𝕀 @( T w )
    in $$( monTabulateQ \ mon ->
              if | isZeroMonomial mon
                 -> [|| f v ||]
                 | Just i <- isLinear mon
                 -> [|| f $$( tabulateQ \ j ->
                                if | j == i
                                   -> [|| 𝕀 1 1 ||]
                                   | otherwise
                                   -> [|| 𝕀 0 0 ||]
                            ) ||]
                 | otherwise
                 -> [|| unT o ||]
         )

  chain :: forall w. Module 𝕀 ( T w ) => D3𝔸1 ( 𝕀ℝ 1 ) -> D3𝔸1 w -> D3𝔸1 w
  chain !df !dg =
    let !o = origin @𝕀 @( T w )
        !p = (^+^)  @𝕀 @( T w )
        !s = (^*)   @𝕀 @( T w )
    in $$( chainRule1NQ
         [|| o ||] [|| p ||] [|| s ||]
         [|| df ||] [|| dg ||] )

instance HasChainRule 𝕀 2 ( 𝕀ℝ 2 ) where

  linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 2 -> w ) -> 𝕀ℝ 2 -> D2𝔸2 w
  linearD f v =
    let !o = origin @𝕀 @( T w )
    in $$( monTabulateQ \ mon ->
              if | isZeroMonomial mon
                 -> [|| f v ||]
                 | Just i <- isLinear mon
                 -> [|| f $$( tabulateQ \ j ->
                                if | j == i
                                   -> [|| 𝕀 1 1 ||]
                                   | otherwise
                                   -> [|| 𝕀 0 0 ||]
                            ) ||]
                 | otherwise
                 -> [|| unT o ||]
         )

  chain :: forall w. Module 𝕀 ( T w ) => D2𝔸1 ( 𝕀ℝ 2 ) -> D2𝔸2 w -> D2𝔸1 w
  chain !df !dg =
    let !o = origin @𝕀 @( T w )
        !p = (^+^)  @𝕀 @( T w )
        !s = (^*)   @𝕀 @( T w )
    in $$( chainRule1NQ
         [|| o ||] [|| p ||] [|| s ||]
         [|| df ||] [|| dg ||] )

instance HasChainRule 𝕀 3 ( 𝕀ℝ 2 ) where

  linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 2 -> w ) -> 𝕀ℝ 2 -> D3𝔸2 w
  linearD f v =
    let !o = origin @𝕀 @( T w )
    in $$( monTabulateQ \ mon ->
              if | isZeroMonomial mon
                 -> [|| f v ||]
                 | Just i <- isLinear mon
                 -> [|| f $$( tabulateQ \ j ->
                                if | j == i
                                   -> [|| 𝕀 1 1 ||]
                                   | otherwise
                                   -> [|| 𝕀 0 0 ||]
                            ) ||]
                 | otherwise
                 -> [|| unT o ||]
         )

  chain :: forall w. Module 𝕀 ( T w ) => D3𝔸1 ( 𝕀ℝ 2 ) -> D3𝔸2 w -> D3𝔸1 w
  chain !df !dg =
    let !o = origin @𝕀 @( T w )
        !p = (^+^)  @𝕀 @( T w )
        !s = (^*)   @𝕀 @( T w )
    in $$( chainRule1NQ
         [|| o ||] [|| p ||] [|| s ||]
         [|| df ||] [|| dg ||] )

instance HasChainRule 𝕀 2 ( 𝕀ℝ 3 ) where

  linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 3 -> w ) -> 𝕀ℝ 3 -> D2𝔸3 w
  linearD f v =
    let !o = origin @𝕀 @( T w )
    in $$( monTabulateQ \ mon ->
              if | isZeroMonomial mon
                 -> [|| f v ||]
                 | Just i <- isLinear mon
                 -> [|| f $$( tabulateQ \ j ->
                                if | j == i
                                   -> [|| 𝕀 1 1 ||]
                                   | otherwise
                                   -> [|| 𝕀 0 0 ||]
                            ) ||]
                 | otherwise
                 -> [|| unT o ||]
         )

  chain :: forall w. Module 𝕀 ( T w ) => D2𝔸1 ( 𝕀ℝ 3 ) -> D2𝔸3 w -> D2𝔸1 w
  chain !df !dg =
    let !o = origin @𝕀 @( T w )
        !p = (^+^)  @𝕀 @( T w )
        !s = (^*)   @𝕀 @( T w )
    in $$( chainRule1NQ
         [|| o ||] [|| p ||] [|| s ||]
         [|| df ||] [|| dg ||] )

instance HasChainRule 𝕀 3 ( 𝕀ℝ 3 ) where

  linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 3 -> w ) -> 𝕀ℝ 3 -> D3𝔸3 w
  linearD f v =
    let !o = origin @𝕀 @( T w )
    in $$( monTabulateQ \ mon ->
              if | isZeroMonomial mon
                 -> [|| f v ||]
                 | Just i <- isLinear mon
                 -> [|| f $$( tabulateQ \ j ->
                                if | j == i
                                   -> [|| 𝕀 1 1 ||]
                                   | otherwise
                                   -> [|| 𝕀 0 0 ||]
                            ) ||]
                 | otherwise
                 -> [|| unT o ||]
         )

  chain :: forall w. Module 𝕀 ( T w ) => D3𝔸1 ( 𝕀ℝ 3 ) -> D3𝔸3 w -> D3𝔸1 w
  chain !df !dg =
    let !o = origin @𝕀 @( T w )
        !p = (^+^)  @𝕀 @( T w )
        !s = (^*)   @𝕀 @( T w )
    in $$( chainRule1NQ
         [|| o ||] [|| p ||] [|| s ||]
         [|| df ||] [|| dg ||] )

instance HasChainRule 𝕀 2 ( 𝕀ℝ 4 ) where

  linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 4 -> w ) -> 𝕀ℝ 4 -> D2𝔸4 w
  linearD f v =
    let !o = origin @𝕀 @( T w )
    in $$( monTabulateQ \ mon ->
              if | isZeroMonomial mon
                 -> [|| f v ||]
                 | Just i <- isLinear mon
                 -> [|| f $$( tabulateQ \ j ->
                                if | j == i
                                   -> [|| 𝕀 1 1 ||]
                                   | otherwise
                                   -> [|| 𝕀 0 0 ||]
                            ) ||]
                 | otherwise
                 -> [|| unT o ||]
         )

  chain :: forall w. Module 𝕀 ( T w ) => D2𝔸1 ( 𝕀ℝ 4 ) -> D2𝔸4 w -> D2𝔸1 w
  chain !df !dg =
    let !o = origin @𝕀 @( T w )
        !p = (^+^)  @𝕀 @( T w )
        !s = (^*)   @𝕀 @( T w )
    in $$( chainRule1NQ
         [|| o ||] [|| p ||] [|| s ||]
         [|| df ||] [|| dg ||] )

instance HasChainRule 𝕀 3 ( 𝕀ℝ 4 ) where

  linearD :: forall w. Module 𝕀 ( T w ) => ( 𝕀ℝ 4 -> w ) -> 𝕀ℝ 4 -> D3𝔸4 w
  linearD f v =
    let !o = origin @𝕀 @( T w )
    in $$( monTabulateQ \ mon ->
              if | isZeroMonomial mon
                 -> [|| f v ||]
                 | Just i <- isLinear mon
                 -> [|| f $$( tabulateQ \ j ->
                                if | j == i
                                   -> [|| 𝕀 1 1 ||]
                                   | otherwise
                                   -> [|| 𝕀 0 0 ||]
                            ) ||]
                 | otherwise
                 -> [|| unT o ||]
         )

  chain :: forall w. Module 𝕀 ( T w ) => D3𝔸1 ( 𝕀ℝ 4 ) -> D3𝔸4 w -> D3𝔸1 w
  chain !df !dg =
    let !o = origin @𝕀 @( T w )
        !p = (^+^)  @𝕀 @( T w )
        !s = (^*)   @𝕀 @( T w )
    in $$( chainRule1NQ
         [|| o ||] [|| p ||] [|| s ||]
         [|| df ||] [|| dg ||] )