brush-strokes-0.1.0.0: src/lib/Math/Root/Isolation/Narrowing.hs
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE UndecidableInstances #-}
module Math.Root.Isolation.Narrowing
( -- * @box(1)@-consistency
Box1
, makeBox1Consistent
-- ** Configuration options
, Box1Options(..)
, defaultBox1Options
, -- *** Narrowing methods for @box(1)@-consistency
NarrowingMethod(..)
, narrowingMethods
-- **** Options for the adaptive shaving method
, AdaptiveShavingOptions(..)
, defaultAdaptiveShavingOptions
-- * @box(2)@-consistency
, Box2
, makeBox2Consistent
-- ** Configuration options
, Box2Options(..)
, defaultBox2Options
)
where
-- base
import Control.Monad
( when )
import Data.Foldable
( toList )
import qualified Data.List.NonEmpty as NE
( toList )
import GHC.TypeNats
( KnownNat )
-- transformers
import Control.Monad.Trans.State.Strict as State
( State, evalState, get, put )
-- brush-strokes
import Math.Algebra.Dual
( D )
import Math.Float.Utils
( succFP, prevFP )
import Math.Interval
import Math.Linear
( ℝ
, Representable
, Fin, set, index, universe
)
import Math.Monomial
( MonomialBasis(..), Deg, Vars
, zeroMonomial, linearMonomial
)
import Math.Root.Isolation.Core
--------------------------------------------------------------------------------
-- Box-consistency driver code
-- | A @box(1)@-consistency enforcing algorithm; see 'makeBox1Consistent'.
data Box1
instance BoxCt n d => RootIsolationAlgorithm Box1 n d where
type instance StepDescription Box1 = String
type instance RootIsolationAlgorithmOptions Box1 n d = Box1Options n d
rootIsolationAlgorithm opts _thisRoundHist _prevRoundsHist eqs box =
return ( "box(1) consistency", makeBox1Consistent opts eqs box )
{-# INLINEABLE rootIsolationAlgorithm #-}
{-# SPECIALISE rootIsolationAlgorithm @2 @3 #-}
-- NB: including this to be safe. The specialiser seems to sometimes
-- be able to generate this specialisation on its own, and sometimes not.
-- | A @box(2)@-consistency enforcing algorithm; see 'makeBox1Consistent'.
data Box2
instance BoxCt n d => RootIsolationAlgorithm Box2 n d where
type instance StepDescription Box2 = String
type instance RootIsolationAlgorithmOptions Box2 n d = Box2Options n d
rootIsolationAlgorithm opts _thisRoundHist _prevRoundsHist eqs box =
return ( "box(2) consistency", [ makeBox2Consistent opts eqs box ] )
{-# INLINEABLE rootIsolationAlgorithm #-}
{-# SPECIALISE rootIsolationAlgorithm @2 @3 #-}
-- NB: including this to be safe. The specialiser seems to sometimes
-- be able to generate this specialisation on its own, and sometimes not.
-- | Options for the @box(1)@-consistency method.
data Box1Options n d =
Box1Options
{ box1EpsEq :: !Double
, box1EpsBis :: !Double
, box1CoordsToNarrow :: !( [ Fin n ] )
, box1EqsToUse :: !( [ Fin d ] )
, box1NarrowingMethod :: !NarrowingMethod
}
deriving stock Show
-- | Options for the @box(2)@-consistency method.
data Box2Options n d =
Box2Options
{ box2Box1Options :: !( Box1Options n d )
, box2EpsEq :: !Double
, box2LambdaMin :: !Double
}
deriving stock Show
-- | Default options for the @box(1)@-consistency method.
defaultBox1Options
:: forall n d
. ( KnownNat n, KnownNat d )
=> Double -- ^ minimum width of boxes (don't bisect further)
-> Double -- ^ threshold for progress
-> Box1Options n d
defaultBox1Options minWidth ε_eq =
Box1Options
{ box1EpsEq = ε_eq
, box1EpsBis = minWidth
, box1CoordsToNarrow = toList $ universe @n
, box1EqsToUse = toList $ universe @d
, box1NarrowingMethod = Kubica
}
{-# INLINEABLE defaultBox1Options #-}
-- | Default options for the @box(2)@-consistency method.
defaultBox2Options
:: forall n d
. ( KnownNat n, KnownNat d )
=> Double -- ^ minimum width of boxes (don't bisect further)
-> Double -- ^ threshold for progress
-> Box2Options n d
defaultBox2Options minWidth ε_eq =
Box2Options
{ box2Box1Options = defaultBox1Options minWidth ε_eq
, box2EpsEq = ε_eq
, box2LambdaMin = 0.001
}
{-# INLINEABLE defaultBox2Options #-}
-- | An implementation of "bc_enforce" from the paper
-- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
--
-- See also
-- "Presentation of a highly tuned multithreaded interval solver for underdetermined and well-determined nonlinear systems"
makeBox1Consistent
:: ( KnownNat n
, Representable Double ( ℝ n )
, Representable 𝕀 ( 𝕀ℝ n )
, MonomialBasis ( D 1 n )
, Deg ( D 1 n ) ~ 1
, Vars ( D 1 n ) ~ n
, Representable Double ( ℝ d )
, Representable 𝕀 ( 𝕀ℝ d )
)
=> Box1Options n d
-> ( 𝕀ℝ n -> D 1 n ( 𝕀ℝ d ) )
-> 𝕀ℝ n -> [ 𝕀ℝ n ]
makeBox1Consistent box1Options eqs x =
( `State.evalState` False ) $
pipeFunctionsWhileTrue ( allNarrowingOperators box1Options eqs ) x
{-# INLINEABLE makeBox1Consistent #-}
{-# SPECIALISE makeBox1Consistent @2 @3 #-}
-- | An implementation of "bound-consistency" from the paper
-- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
makeBox2Consistent
:: forall n d
. ( KnownNat n
, Representable Double ( ℝ n )
, Representable 𝕀 ( 𝕀ℝ n )
, MonomialBasis ( D 1 n )
, Deg ( D 1 n ) ~ 1
, Vars ( D 1 n ) ~ n
, Representable Double ( ℝ d )
, Representable 𝕀 ( 𝕀ℝ d )
)
=> Box2Options n d
-> ( 𝕀ℝ n -> D 1 n ( 𝕀ℝ d ) )
-> 𝕀ℝ n -> 𝕀ℝ n
makeBox2Consistent (Box2Options box1Options ε_eq λMin) eqs x0
= ( `State.evalState` False ) $ doLoop 0.25 x0
where
doBox1 :: 𝕀ℝ n -> [ 𝕀ℝ n ]
doBox1 = makeBox1Consistent box1Options eqs
doLoop :: Double -> 𝕀ℝ n -> State Bool ( 𝕀ℝ n )
doLoop λ x = do
x'' <- forEachCoord @n x $ boundConsistency λ
modified <- State.get
let λ' = if modified then λ else 0.5 * λ
if λ' < λMin
then return x''
else do { State.put False ; doLoop λ' x'' }
boundConsistency :: Double -> Fin n -> 𝕀ℝ n -> State Bool ( 𝕀ℝ n )
boundConsistency λ i box = do
let x@( 𝕀 x_inf x_sup ) = getter box
c1 = ( 1 - λ ) * x_inf + λ * x_sup
c2 = λ * x_inf + ( 1 - λ ) * x_sup
x'_inf =
case doBox1 ( setter ( 𝕀 x_inf c1 ) box ) of
[] -> c1
x's -> minimum $ map ( inf . getter ) x's
x'_sup =
case doBox1 ( setter ( 𝕀 c2 x_sup ) box ) of
[] -> c2
x's -> maximum $ map ( sup . getter ) x's
x' = 𝕀 x'_inf x'_sup
when ( width x - width x' >= ε_eq ) $
State.put True
return $ setter x' box
where
getter = ( `index` i )
setter = set i
{-# INLINEABLE makeBox2Consistent #-}
{-# SPECIALISE makeBox2Consistent @2 @3 #-}
--------------------------------------------------------------------------------
-- Narrowing methods
-- | The narrowing method to use to enforce @box(1)@-consistency.
data NarrowingMethod
-- | Algorithm 5 from the paper
-- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
-- (Bartłomiej Jacek Kubica, 2017)
= Kubica
-- | Two-sided shaving @sbc@ from the paper
-- "A Data-Parallel Algorithm to Reliably Solve Systems of Nonlinear Equations",
-- (Goldsztejn & Goualard, 2008)
| TwoSidedShaving
-- | @sbc3ag@ (adaptive guessing) shaving, from the paper
-- "Box Consistency through Adaptive Shaving"
-- (Goldsztejn & Goualard, 2010).
| AdaptiveShaving AdaptiveShavingOptions
deriving stock Show
narrowingMethods
:: Double -> Double
-> NarrowingMethod
-> [ ( 𝕀 -> ( 𝕀, 𝕀 ) ) -> 𝕀 -> [ 𝕀 ] ]
narrowingMethods ε_eq ε_bis Kubica
= [ leftNarrow ε_eq ε_bis, rightNarrow ε_eq ε_bis ]
narrowingMethods ε_eq ε_bis ( AdaptiveShaving opts )
= [ leftShave ε_eq ε_bis opts, rightNarrow ε_eq ε_bis ]
-- TODO: haven't implemented right shaving yet
narrowingMethods _ε_eq ε_bis TwoSidedShaving
= [ sbc ε_bis ]
{-# INLINEABLE narrowingMethods #-}
allNarrowingOperators
:: forall n d
. ( KnownNat n
, Representable Double ( ℝ n )
, Representable 𝕀 ( 𝕀ℝ n )
, MonomialBasis ( D 1 n )
, Deg ( D 1 n ) ~ 1
, Vars ( D 1 n ) ~ n
, Representable Double ( ℝ d )
, Representable 𝕀 ( 𝕀ℝ d )
)
=> Box1Options n d
-> ( 𝕀ℝ n -> D 1 n ( 𝕀ℝ d ) )
-> [ 𝕀ℝ n -> State Bool [ 𝕀ℝ n ] ]
allNarrowingOperators ( Box1Options ε_eq ε_bis coordsToNarrow eqsToUse narrowingMethod ) eqs =
[ \ cand ->
let getter = ( `index` coordIndex )
setter = set coordIndex
newCands = map ( `setter` cand )
$ narrowFn ( \ box -> ff' coordIndex eqnIndex $ setter box cand )
( getter cand )
in do
-- Record when we achieved a meaningful reduction,
-- so that we continue trying further narrowings.
when ( ( width ( getter cand ) - sum ( map ( width . getter ) newCands ) ) >= ε_eq ) $
-- NB: making this return 'True' less often seems slightly beneficial?
-- Further investigation needed.
State.put True
return newCands
| narrowFn <- narrowingMethods ε_eq ε_bis narrowingMethod
, coordIndex <- coordsToNarrow
, eqnIndex <- eqsToUse
]
where
ff' :: Fin n -> Fin d -> 𝕀ℝ n -> ( 𝕀, 𝕀 )
ff' i d ts =
let df = eqs ts
f, f' :: 𝕀ℝ d
f = df `monIndex` zeroMonomial
f' = df `monIndex` linearMonomial i
in ( f `index` d, f' `index` d )
{-# INLINEABLE allNarrowingOperators #-}
{-# SPECIALISE allNarrowingOperators @2 @3 #-}
--------------------------------------------------------------------------------
-- Kubica's algorithm.
-- Use the univariate interval Newton method to narrow from the left
-- a candidate interval.
--
-- See Algorithm 5 (Procedure left_narrow) in
-- "Parallelization of a bound-consistency enforcing procedure and its application in solving nonlinear systems"
-- by Bartłomiej Jacek Kubica, 2017
leftNarrow :: Double
-> Double
-> ( 𝕀 -> ( 𝕀, 𝕀 ) )
-> 𝕀
-> [ 𝕀 ]
leftNarrow ε_eq ε_bis ff' = left_narrow
where
left_narrow ( 𝕀 x_inf x_sup ) =
go x_sup ( 𝕀 x_inf ( if x_inf == x_sup then x_inf else succFP x_inf ) )
go x_sup x_left =
let ( f_x_left, _f'_x_left ) = ff' x_left
in
if inf f_x_left <= 0 && sup f_x_left >= 0
then [ 𝕀 ( inf x_left ) x_sup ]
else
let x = 𝕀 ( sup x_left ) x_sup
( _f_x, f'_x ) = ff' x
x's = do δ <- f_x_left ⊘ f'_x
let x_new = x_left - δ
map fst $ x_new ∩ x
in
if | null x's
-> []
| ( width x - sum ( map width x's ) ) < ε_eq
-> x's
| otherwise
-> do
x' <- x's
if sup x' - inf x' < ε_bis
then return x'
else left_narrow =<< NE.toList ( bisect x' )
-- TODO: de-duplicate with 'leftNarrow'?
rightNarrow :: Double
-> Double
-> ( 𝕀 -> ( 𝕀, 𝕀 ) )
-> 𝕀
-> [ 𝕀 ]
rightNarrow ε_eq ε_bis ff' = right_narrow
where
right_narrow ( 𝕀 x_inf x_sup ) =
go x_inf ( 𝕀 ( if x_inf == x_sup then x_sup else prevFP x_sup ) x_sup )
go x_inf x_right =
let ( f_x_right, _f'_x_right ) = ff' x_right
in
if inf f_x_right <= 0 && sup f_x_right >= 0
then [ 𝕀 x_inf ( sup x_right ) ]
else
let x = 𝕀 x_inf ( inf x_right )
( _f_x, f'_x ) = ff' x
x's = do δ <- f_x_right ⊘ f'_x
let x_new = x_right - δ
map fst $ x_new ∩ x
in
if | null x's
-> []
| ( width x - sum ( map width x's ) ) < ε_eq
-> x's
| otherwise
-> do
x' <- x's
if sup x' - inf x' < ε_bis
then return x'
else right_narrow =<< NE.toList ( bisect x' )
--------------------------------------------------------------------------------
-- Adaptive shaving.
data AdaptiveShavingOptions
= AdaptiveShavingOptions
{ γ_init, σ_good, σ_bad, β_good, β_bad :: !Double
}
deriving stock Show
defaultAdaptiveShavingOptions :: AdaptiveShavingOptions
defaultAdaptiveShavingOptions =
AdaptiveShavingOptions
{ γ_init = 0.25
, σ_good = 0.25
, σ_bad = 0.75
, β_good = 1.5
, β_bad = 0.7
}
-- | Algorithm @lnar_sbc3ag@ (adaptive guessing) from the paper
-- "Box Consistency through Adaptive Shaving" (Goldsztejn & Goualard, 2010).
leftShave :: Double
-> Double
-> AdaptiveShavingOptions
-> ( 𝕀 -> ( 𝕀, 𝕀 ) )
-> 𝕀
-> [ 𝕀 ]
leftShave ε_eq ε_bis
( AdaptiveShavingOptions { γ_init, σ_good, σ_bad, β_good, β_bad } )
ff' i0 =
left_narrow γ_init i0
where
w0 = width i0
left_narrow :: Double -> 𝕀 -> [ 𝕀 ]
left_narrow γ i@( 𝕀 x_inf x_sup )
-- Stop if the box is too small.
| width i < ε_bis
= [ i ]
| otherwise
= go γ x_sup ( 𝕀 x_inf ( if x_inf == x_sup then x_inf else succFP x_inf ) )
go :: Double -> Double -> 𝕀 -> [ 𝕀 ]
go γ x_sup x_left =
let ( f_x_left, _f'_x_left ) = ff' x_left
in
if 0 ∈ f_x_left
-- Box-consistency achieved; finish.
then [ 𝕀 ( inf x_left ) x_sup ]
else
-- Otherwise, try to shave off a chunk on the left of the interval.
let x' = 𝕀 ( sup x_left ) x_sup
inf_guess = sup x_left
sup_guess = min x_sup ( inf_guess + γ * w0 ) -- * width x' )
guess = 𝕀 inf_guess sup_guess
-- NB: this always uses the initial width (to "avoid asymptotic behaviour" according to the paper)
( f_guess, f'_guess ) = ff' guess
x_minus_guess = 𝕀 ( min x_sup ( succFP $ sup guess ) ) x_sup
in if not ( 0 ∈ f_guess )
then
-- We successfully shaved "guess" off; go round again after removing it.
-- TODO: here we could go back to the top with a new "w" maybe?
left_narrow ( β_good * γ ) x_minus_guess
else
-- Do a Newton step to try to reduce the guess interval.
-- Starting from "guess", we get a collection of sub-intervals
-- "guesses'" refining where the function can be zero.
let guesses' :: [ 𝕀 ]
guesses' = do
δ <- f_x_left ⊘ f'_guess
let guess' = singleton ( inf guess ) - δ
map fst $ guess' ∩ guess
w_guess = width guess
w_guesses'
| null guesses'
= 0
| otherwise
= sup_guess - minimum ( map inf guesses' )
γ' | w_guesses' < σ_good * w_guess
-- Good improvement, try larger guesses in the future.
= β_good * γ
| w_guesses' > σ_bad * w_guess
-- Poor improvement, try smaller guesses in the future.
= β_bad * γ
| otherwise
-- Otherwise, keep the γ factor the same.
= γ
xs' = x_minus_guess : guesses'
in if ( width x' - sum ( map width xs' ) ) < ε_eq
then xs'
else left_narrow γ' =<< xs'
{-# INLINEABLE leftShave #-}
--------------------------------------------------------------------------------
-- Two-sided shaving.
-- | @sbc@ algorithm from the paper
--
-- "A Data-Parallel Algorithm to Reliably Solve Systems of Nonlinear Equations",
-- (Frédéric Goualard, Alexandre Goldsztejn, 2008).
sbc :: Double
-> ( 𝕀 -> ( 𝕀, 𝕀 ) )
-> 𝕀 -> [ 𝕀 ]
sbc ε ff' = go
where
go :: 𝕀 -> [ 𝕀 ]
go x@( 𝕀 x_l x_r )
| width x <= ε
= [ x ]
| otherwise
= let x_mid = 0.5 * ( x_l + x_r )
( left_done, left_todo )
| 0 ∈ ( fst $ ff' ( 𝕀 x_l x_lp ) )
= ( True, [ ] )
| not $ 0 ∈ ( fst $ ff' i_l )
= ( False, [ ] )
| otherwise
= let
xls = do
let l = 𝕀 x_lp x_lp
δ <- fst ( ff' l ) ⊘ snd ( ff' i_l )
map fst $ ( l - δ ) ∩ i_l
in ( False, xls )
where x_lp = min ( succFP x_l ) x_mid
i_l = 𝕀 x_lp x_mid
( right_done, right_todo )
| 0 ∈ ( fst $ ff' ( 𝕀 x_rm x_r ) )
= ( True, [ ] )
| not $ 0 ∈ ( fst $ ff' i_r )
= ( False, [ ] )
| otherwise
= let
xrs = do
let r = 𝕀 x_rm x_rm
δ <- fst ( ff' r ) ⊘ snd ( ff' i_r )
map fst $ ( r - δ ) ∩ i_r
in ( False, xrs )
where x_rm = max ( prevFP x_r ) x_mid
i_r = 𝕀 x_mid x_rm
in do let lefts' = if left_done
then [ 𝕀 x_l x_mid ]
else go =<< left_todo
rights' = if right_done
then [ 𝕀 x_mid x_r ]
else go =<< right_todo
lefts' ++ rights'
{-# INLINEABLE sbc #-}