{-# LANGUAGE DeriveFoldable #-}
{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE DeriveTraversable #-}
{-# LANGUAGE PatternSynonyms #-}
{-# LANGUAGE ViewPatterns #-}

-- SPDX-License-Identifier: MPL-2.0

{- | Uniform iteration reports retain the best current estimate on an exhausted
budget. Input/domain failures use Either; successful execution need not mean
numerical convergence. Residuals are finite nonnegative absolute measures.
-}
module HMath.API.Convergence (
    StopReason (..),
    Convergence (..),
    IterationOptions (..),
    defaultIterationOptions,
    iterateUntil,
    requireConvergence,
    resultMaybe,
    fromMaybeError,
    pattern Converged,
) where

import HMath.Error

{- | Stagnation is opt-in and uses an exact-equality caller predicate; it does
not reinterpret a small update as proof of a small residual.
-}
data StopReason = ResidualSatisfied | IterationLimit | Stagnation deriving (StopReason -> StopReason -> Bool
(StopReason -> StopReason -> Bool)
-> (StopReason -> StopReason -> Bool) -> Eq StopReason
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: StopReason -> StopReason -> Bool
== :: StopReason -> StopReason -> Bool
$c/= :: StopReason -> StopReason -> Bool
/= :: StopReason -> StopReason -> Bool
Eq, Int -> StopReason -> ShowS
[StopReason] -> ShowS
StopReason -> String
(Int -> StopReason -> ShowS)
-> (StopReason -> String)
-> ([StopReason] -> ShowS)
-> Show StopReason
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> StopReason -> ShowS
showsPrec :: Int -> StopReason -> ShowS
$cshow :: StopReason -> String
show :: StopReason -> String
$cshowList :: [StopReason] -> ShowS
showList :: [StopReason] -> ShowS
Show)

{- | Iteration report retaining the estimate and residual even on a budget/stagnation stop.
Inspect stopReason before assuming convergence.
-}
data Convergence a = Convergence
    { forall a. Convergence a -> a
estimate :: a
    -- ^ Best current estimate, retained even when the iteration budget is exhausted.
    , forall a. Convergence a -> Int
iterationCount :: !Int
    -- ^ Number of completed updates; the initial estimate is iteration zero.
    , forall a. Convergence a -> Double
finalResidual :: !Double
    -- ^ Finite nonnegative absolute residual at the returned estimate.
    , forall a. Convergence a -> StopReason
stopReason :: !StopReason
    -- ^ ResidualSatisfied,IterationLimit,or Stagnation; successful execution does not itself imply convergence.
    }
    deriving (Convergence a -> Convergence a -> Bool
(Convergence a -> Convergence a -> Bool)
-> (Convergence a -> Convergence a -> Bool) -> Eq (Convergence a)
forall a. Eq a => Convergence a -> Convergence a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => Convergence a -> Convergence a -> Bool
== :: Convergence a -> Convergence a -> Bool
$c/= :: forall a. Eq a => Convergence a -> Convergence a -> Bool
/= :: Convergence a -> Convergence a -> Bool
Eq, Int -> Convergence a -> ShowS
[Convergence a] -> ShowS
Convergence a -> String
(Int -> Convergence a -> ShowS)
-> (Convergence a -> String)
-> ([Convergence a] -> ShowS)
-> Show (Convergence a)
forall a. Show a => Int -> Convergence a -> ShowS
forall a. Show a => [Convergence a] -> ShowS
forall a. Show a => Convergence a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> Convergence a -> ShowS
showsPrec :: Int -> Convergence a -> ShowS
$cshow :: forall a. Show a => Convergence a -> String
show :: Convergence a -> String
$cshowList :: forall a. Show a => [Convergence a] -> ShowS
showList :: [Convergence a] -> ShowS
Show, (forall a b. (a -> b) -> Convergence a -> Convergence b)
-> (forall a b. a -> Convergence b -> Convergence a)
-> Functor Convergence
forall a b. a -> Convergence b -> Convergence a
forall a b. (a -> b) -> Convergence a -> Convergence b
forall (f :: * -> *).
(forall a b. (a -> b) -> f a -> f b)
-> (forall a b. a -> f b -> f a) -> Functor f
$cfmap :: forall a b. (a -> b) -> Convergence a -> Convergence b
fmap :: forall a b. (a -> b) -> Convergence a -> Convergence b
$c<$ :: forall a b. a -> Convergence b -> Convergence a
<$ :: forall a b. a -> Convergence b -> Convergence a
Functor, (forall m. Monoid m => Convergence m -> m)
-> (forall m a. Monoid m => (a -> m) -> Convergence a -> m)
-> (forall m a. Monoid m => (a -> m) -> Convergence a -> m)
-> (forall a b. (a -> b -> b) -> b -> Convergence a -> b)
-> (forall a b. (a -> b -> b) -> b -> Convergence a -> b)
-> (forall b a. (b -> a -> b) -> b -> Convergence a -> b)
-> (forall b a. (b -> a -> b) -> b -> Convergence a -> b)
-> (forall a. (a -> a -> a) -> Convergence a -> a)
-> (forall a. (a -> a -> a) -> Convergence a -> a)
-> (forall a. Convergence a -> [a])
-> (forall a. Convergence a -> Bool)
-> (forall a. Convergence a -> Int)
-> (forall a. Eq a => a -> Convergence a -> Bool)
-> (forall a. Ord a => Convergence a -> a)
-> (forall a. Ord a => Convergence a -> a)
-> (forall a. Num a => Convergence a -> a)
-> (forall a. Num a => Convergence a -> a)
-> Foldable Convergence
forall a. Eq a => a -> Convergence a -> Bool
forall a. Num a => Convergence a -> a
forall a. Ord a => Convergence a -> a
forall m. Monoid m => Convergence m -> m
forall a. Convergence a -> Bool
forall a. Convergence a -> Int
forall a. Convergence a -> [a]
forall a. (a -> a -> a) -> Convergence a -> a
forall m a. Monoid m => (a -> m) -> Convergence a -> m
forall b a. (b -> a -> b) -> b -> Convergence a -> b
forall a b. (a -> b -> b) -> b -> Convergence a -> b
forall (t :: * -> *).
(forall m. Monoid m => t m -> m)
-> (forall m a. Monoid m => (a -> m) -> t a -> m)
-> (forall m a. Monoid m => (a -> m) -> t a -> m)
-> (forall a b. (a -> b -> b) -> b -> t a -> b)
-> (forall a b. (a -> b -> b) -> b -> t a -> b)
-> (forall b a. (b -> a -> b) -> b -> t a -> b)
-> (forall b a. (b -> a -> b) -> b -> t a -> b)
-> (forall a. (a -> a -> a) -> t a -> a)
-> (forall a. (a -> a -> a) -> t a -> a)
-> (forall a. t a -> [a])
-> (forall a. t a -> Bool)
-> (forall a. t a -> Int)
-> (forall a. Eq a => a -> t a -> Bool)
-> (forall a. Ord a => t a -> a)
-> (forall a. Ord a => t a -> a)
-> (forall a. Num a => t a -> a)
-> (forall a. Num a => t a -> a)
-> Foldable t
$cfold :: forall m. Monoid m => Convergence m -> m
fold :: forall m. Monoid m => Convergence m -> m
$cfoldMap :: forall m a. Monoid m => (a -> m) -> Convergence a -> m
foldMap :: forall m a. Monoid m => (a -> m) -> Convergence a -> m
$cfoldMap' :: forall m a. Monoid m => (a -> m) -> Convergence a -> m
foldMap' :: forall m a. Monoid m => (a -> m) -> Convergence a -> m
$cfoldr :: forall a b. (a -> b -> b) -> b -> Convergence a -> b
foldr :: forall a b. (a -> b -> b) -> b -> Convergence a -> b
$cfoldr' :: forall a b. (a -> b -> b) -> b -> Convergence a -> b
foldr' :: forall a b. (a -> b -> b) -> b -> Convergence a -> b
$cfoldl :: forall b a. (b -> a -> b) -> b -> Convergence a -> b
foldl :: forall b a. (b -> a -> b) -> b -> Convergence a -> b
$cfoldl' :: forall b a. (b -> a -> b) -> b -> Convergence a -> b
foldl' :: forall b a. (b -> a -> b) -> b -> Convergence a -> b
$cfoldr1 :: forall a. (a -> a -> a) -> Convergence a -> a
foldr1 :: forall a. (a -> a -> a) -> Convergence a -> a
$cfoldl1 :: forall a. (a -> a -> a) -> Convergence a -> a
foldl1 :: forall a. (a -> a -> a) -> Convergence a -> a
$ctoList :: forall a. Convergence a -> [a]
toList :: forall a. Convergence a -> [a]
$cnull :: forall a. Convergence a -> Bool
null :: forall a. Convergence a -> Bool
$clength :: forall a. Convergence a -> Int
length :: forall a. Convergence a -> Int
$celem :: forall a. Eq a => a -> Convergence a -> Bool
elem :: forall a. Eq a => a -> Convergence a -> Bool
$cmaximum :: forall a. Ord a => Convergence a -> a
maximum :: forall a. Ord a => Convergence a -> a
$cminimum :: forall a. Ord a => Convergence a -> a
minimum :: forall a. Ord a => Convergence a -> a
$csum :: forall a. Num a => Convergence a -> a
sum :: forall a. Num a => Convergence a -> a
$cproduct :: forall a. Num a => Convergence a -> a
product :: forall a. Num a => Convergence a -> a
Foldable, Functor Convergence
Foldable Convergence
(Functor Convergence, Foldable Convergence) =>
(forall (f :: * -> *) a b.
 Applicative f =>
 (a -> f b) -> Convergence a -> f (Convergence b))
-> (forall (f :: * -> *) a.
    Applicative f =>
    Convergence (f a) -> f (Convergence a))
-> (forall (m :: * -> *) a b.
    Monad m =>
    (a -> m b) -> Convergence a -> m (Convergence b))
-> (forall (m :: * -> *) a.
    Monad m =>
    Convergence (m a) -> m (Convergence a))
-> Traversable Convergence
forall (t :: * -> *).
(Functor t, Foldable t) =>
(forall (f :: * -> *) a b.
 Applicative f =>
 (a -> f b) -> t a -> f (t b))
-> (forall (f :: * -> *) a. Applicative f => t (f a) -> f (t a))
-> (forall (m :: * -> *) a b.
    Monad m =>
    (a -> m b) -> t a -> m (t b))
-> (forall (m :: * -> *) a. Monad m => t (m a) -> m (t a))
-> Traversable t
forall (m :: * -> *) a.
Monad m =>
Convergence (m a) -> m (Convergence a)
forall (f :: * -> *) a.
Applicative f =>
Convergence (f a) -> f (Convergence a)
forall (m :: * -> *) a b.
Monad m =>
(a -> m b) -> Convergence a -> m (Convergence b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> Convergence a -> f (Convergence b)
$ctraverse :: forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> Convergence a -> f (Convergence b)
traverse :: forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> Convergence a -> f (Convergence b)
$csequenceA :: forall (f :: * -> *) a.
Applicative f =>
Convergence (f a) -> f (Convergence a)
sequenceA :: forall (f :: * -> *) a.
Applicative f =>
Convergence (f a) -> f (Convergence a)
$cmapM :: forall (m :: * -> *) a b.
Monad m =>
(a -> m b) -> Convergence a -> m (Convergence b)
mapM :: forall (m :: * -> *) a b.
Monad m =>
(a -> m b) -> Convergence a -> m (Convergence b)
$csequence :: forall (m :: * -> *) a.
Monad m =>
Convergence (m a) -> m (Convergence a)
sequence :: forall (m :: * -> *) a.
Monad m =>
Convergence (m a) -> m (Convergence a)
Traversable)

{- | Caller-controlled iteration count and absolute/relative stopping tolerances, validated
by iterateUntil.
-}
data IterationOptions = IterationOptions
    { IterationOptions -> Int
maximumIterations :: !Int
    -- ^ Nonnegative maximum number of updates permitted by iterateUntil.
    , IterationOptions -> Double
absoluteTolerance :: !Double
    -- ^ Finite nonnegative absolute stopping threshold.
    , IterationOptions -> Double
relativeTolerance :: !Double
    -- ^ Finite nonnegative tolerance scaled by the initial residual.
    }
    deriving (IterationOptions -> IterationOptions -> Bool
(IterationOptions -> IterationOptions -> Bool)
-> (IterationOptions -> IterationOptions -> Bool)
-> Eq IterationOptions
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: IterationOptions -> IterationOptions -> Bool
== :: IterationOptions -> IterationOptions -> Bool
$c/= :: IterationOptions -> IterationOptions -> Bool
/= :: IterationOptions -> IterationOptions -> Bool
Eq, Int -> IterationOptions -> ShowS
[IterationOptions] -> ShowS
IterationOptions -> String
(Int -> IterationOptions -> ShowS)
-> (IterationOptions -> String)
-> ([IterationOptions] -> ShowS)
-> Show IterationOptions
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> IterationOptions -> ShowS
showsPrec :: Int -> IterationOptions -> ShowS
$cshow :: IterationOptions -> String
show :: IterationOptions -> String
$cshowList :: [IterationOptions] -> ShowS
showList :: [IterationOptions] -> ShowS
Show)

{- | Default iteration budget 1000, absolute tolerance 1e-12, relative tolerance 1e-10. Use
record updates to configure a solver; iterateUntil validates nonnegative finite
tolerances.
-}
defaultIterationOptions :: IterationOptions
defaultIterationOptions :: IterationOptions
defaultIterationOptions = Int -> Double -> Double -> IterationOptions
IterationOptions Int
1000 Double
1e-12 Double
1e-10

-- | Read-only successful view. The report still retains its residual and count.
pattern Converged :: a -> Convergence a
pattern $mConverged :: forall {r} {a}. Convergence a -> (a -> r) -> ((# #) -> r) -> r
Converged x <- Convergence x _ _ ResidualSatisfied

{- | Stop at residual <= atol + rtol * initialResidual. Comparisons are scaled
to avoid overflow in that sum. The initial value is tested at iteration zero.
-}
iterateUntil ::
    IterationOptions ->
    (a -> Either MathError a) ->
    (a -> Either MathError Double) ->
    (a -> a -> Bool) ->
    a ->
    Either MathError (Convergence a)
iterateUntil :: forall a.
IterationOptions
-> (a -> Either MathError a)
-> (a -> Either MathError Double)
-> (a -> a -> Bool)
-> a
-> Either MathError (Convergence a)
iterateUntil IterationOptions
options a -> Either MathError a
step a -> Either MathError Double
residual a -> a -> Bool
stagnant a
initial = do
    let limit :: Int
limit = IterationOptions -> Int
maximumIterations IterationOptions
options
        atol :: Double
atol = IterationOptions -> Double
absoluteTolerance IterationOptions
options
        rtol :: Double
rtol = IterationOptions -> Double
relativeTolerance IterationOptions
options
    if Int
limit Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 Bool -> Bool -> Bool
|| Bool -> Bool
not (Double -> Bool
forall {a}. RealFloat a => a -> Bool
finite Double
atol Bool -> Bool -> Bool
&& Double -> Bool
forall {a}. RealFloat a => a -> Bool
finite Double
rtol) Bool -> Bool -> Bool
|| Double
atol Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 Bool -> Bool -> Bool
|| Double
rtol Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0
        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"invalid convergence budget/tolerances")
        else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    initialResidual <- a -> Either MathError Double
residual a
initial Either MathError Double
-> (Double -> Either MathError Double) -> Either MathError Double
forall a b.
Either MathError a
-> (a -> Either MathError b) -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= Double -> Either MathError Double
forall {b}. RealFloat b => b -> Either MathError b
checkResidual
    loop initialResidual 0 initial initialResidual
  where
    finite :: a -> Bool
finite a
x = Bool -> Bool
not (a -> Bool
forall {a}. RealFloat a => a -> Bool
isNaN a
x Bool -> Bool -> Bool
|| a -> Bool
forall {a}. RealFloat a => a -> Bool
isInfinite a
x)
    checkResidual :: b -> Either MathError b
checkResidual b
x
        | b -> Bool
forall {a}. RealFloat a => a -> Bool
finite b
x Bool -> Bool -> Bool
&& b
x b -> b -> Bool
forall a. Ord a => a -> a -> Bool
>= b
0 = b -> Either MathError b
forall a b. b -> Either a b
Right b
x
        | Bool
otherwise = MathError -> Either MathError b
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"residual must be finite and nonnegative")
    satisfied :: Double -> Double -> Bool
satisfied Double
initialResidual Double
current
        | Double
current Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= IterationOptions -> Double
absoluteTolerance IterationOptions
options = Bool
True
        | Double
initialResidual Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 = Bool
False
        | Bool
otherwise = (Double
current Double -> Double -> Double
forall a. Num a => a -> a -> a
- IterationOptions -> Double
absoluteTolerance IterationOptions
options) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
initialResidual Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= IterationOptions -> Double
relativeTolerance IterationOptions
options
    loop :: Double -> Int -> a -> Double -> Either MathError (Convergence a)
loop Double
initialResidual Int
n a
x Double
r
        | Double -> Double -> Bool
satisfied Double
initialResidual Double
r = Convergence a -> Either MathError (Convergence a)
forall a b. b -> Either a b
Right (a -> Int -> Double -> StopReason -> Convergence a
forall a. a -> Int -> Double -> StopReason -> Convergence a
Convergence a
x Int
n Double
r StopReason
ResidualSatisfied)
        | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= IterationOptions -> Int
maximumIterations IterationOptions
options = Convergence a -> Either MathError (Convergence a)
forall a b. b -> Either a b
Right (a -> Int -> Double -> StopReason -> Convergence a
forall a. a -> Int -> Double -> StopReason -> Convergence a
Convergence a
x Int
n Double
r StopReason
IterationLimit)
        | Bool
otherwise = do
            y <- a -> Either MathError a
step a
x
            ry <- residual y >>= checkResidual
            if not (satisfied initialResidual ry) && stagnant x y
                then Right (Convergence y (n + 1) ry Stagnation)
                else loop initialResidual (n + 1) y ry

{- | Given operation name then convergence report, return its estimate only for
ResidualSatisfied. IterationLimit or Stagnation becomes Left (NonConvergence name
iterationCount), preserving the distinction between execution and convergence.
-}
requireConvergence :: String -> Convergence a -> Either MathError a
requireConvergence :: forall a. String -> Convergence a -> Either MathError a
requireConvergence String
name Convergence a
r
    | Convergence a -> StopReason
forall a. Convergence a -> StopReason
stopReason Convergence a
r StopReason -> StopReason -> Bool
forall a. Eq a => a -> a -> Bool
== StopReason
ResidualSatisfied = a -> Either MathError a
forall a b. b -> Either a b
Right (Convergence a -> a
forall a. Convergence a -> a
estimate Convergence a
r)
    | Bool
otherwise = MathError -> Either MathError a
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
name (Convergence a -> Int
forall a. Convergence a -> Int
iterationCount Convergence a
r))

-- | Explicitly discard error details for callers retaining a Maybe interface.
resultMaybe :: Either e a -> Maybe a
resultMaybe :: forall e a. Either e a -> Maybe a
resultMaybe = (e -> Maybe a) -> (a -> Maybe a) -> Either e a -> Maybe a
forall a c b. (a -> c) -> (b -> c) -> Either a b -> c
either (Maybe a -> e -> Maybe a
forall a b. a -> b -> a
const Maybe a
forall a. Maybe a
Nothing) a -> Maybe a
forall a. a -> Maybe a
Just

{- | Given error value then Maybe result, convert Nothing to Left error and Just x to Right
x. Example: fromMaybeError "missing" (Just 3) == Right 3.
-}
fromMaybeError :: e -> Maybe a -> Either e a
fromMaybeError :: forall e a. e -> Maybe a -> Either e a
fromMaybeError e
e = Either e a -> (a -> Either e a) -> Maybe a -> Either e a
forall b a. b -> (a -> b) -> Maybe a -> b
maybe (e -> Either e a
forall a b. a -> Either a b
Left e
e) a -> Either e a
forall a b. b -> Either a b
Right