-- SPDX-License-Identifier: MPL-2.0

{- | Checked, finite-difference differentiation and numerical quadrature for
Double callbacks. User callbacks must be total and side-effect free. Budgets
bound refinement; error estimates are numerical estimates, not rigorous bounds.
-}
module HMath.Calculus (
    Estimate (..),
    Tolerance (..),
    defaultTolerance,
    Difference (..),
    derivative,
    derivativeN,
    complexDerivative,
    richardson,
    ridders,
    limit,
    partialDerivative,
    gradient,
    jacobian,
    hessian,
    directionalDerivative,
    jvp,
    vjp,
    hvp,
    divergence,
    curl,
    laplacian,
    Quadrature (..),
    integrateComposite,
    romberg,
    gaussLegendre,
    gaussHermite,
    gaussLaguerre,
    gaussChebyshev,
    gaussKronrod,
    adaptiveIntegrate,
    clenshawCurtis,
    tanhSinh,
    integrateInfinite,
    integrateOscillatory,
    integrateND,
    monteCarlo,
    quasiMonteCarlo,
    cumulativeIntegral,
    pathIntegral,
    lineIntegral,
    surfaceIntegral,
    volumeIntegral,
) where

import Data.Complex (Complex ((:+)), imagPart, realPart)
import Data.List (transpose)
import HMath.Error (MathError (..))
import HMath.Matrix (Matrix, fromListChecked)

-- | Estimated value, absolute error estimate and callback evaluation count.
data Estimate a = Estimate
    { forall a. Estimate a -> a
estimateValue :: a
    -- ^ Returned numerical derivative/integral estimate.
    , forall a. Estimate a -> Double
estimateError :: Double
    -- ^ Absolute error estimate from the selected approximation/refinement rule; not a universal rigorous bound.
    , forall a. Estimate a -> Int
evaluations :: Int
    -- ^ Number of function evaluations used by the calculation.
    }
    deriving (Estimate a -> Estimate a -> Bool
(Estimate a -> Estimate a -> Bool)
-> (Estimate a -> Estimate a -> Bool) -> Eq (Estimate a)
forall a. Eq a => Estimate a -> Estimate a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => Estimate a -> Estimate a -> Bool
== :: Estimate a -> Estimate a -> Bool
$c/= :: forall a. Eq a => Estimate a -> Estimate a -> Bool
/= :: Estimate a -> Estimate a -> Bool
Eq, Int -> Estimate a -> ShowS
[Estimate a] -> ShowS
Estimate a -> String
(Int -> Estimate a -> ShowS)
-> (Estimate a -> String)
-> ([Estimate a] -> ShowS)
-> Show (Estimate a)
forall a. Show a => Int -> Estimate a -> ShowS
forall a. Show a => [Estimate a] -> ShowS
forall a. Show a => Estimate a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> Estimate a -> ShowS
showsPrec :: Int -> Estimate a -> ShowS
$cshow :: forall a. Show a => Estimate a -> String
show :: Estimate a -> String
$cshowList :: forall a. Show a => [Estimate a] -> ShowS
showList :: [Estimate a] -> ShowS
Show)

-- | Positive absolute tolerance, nonnegative relative tolerance and positive budget.
data Tolerance = Tolerance
    { Tolerance -> Double
absoluteTolerance :: Double
    -- ^ Requested finite nonnegative absolute error threshold.
    , Tolerance -> Double
relativeTolerance :: Double
    -- ^ Requested finite nonnegative error threshold relative to estimate magnitude.
    , Tolerance -> Int
maximumEvaluations :: Int
    -- ^ Positive caller-supplied function-evaluation budget.
    }
    deriving (Tolerance -> Tolerance -> Bool
(Tolerance -> Tolerance -> Bool)
-> (Tolerance -> Tolerance -> Bool) -> Eq Tolerance
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Tolerance -> Tolerance -> Bool
== :: Tolerance -> Tolerance -> Bool
$c/= :: Tolerance -> Tolerance -> Bool
/= :: Tolerance -> Tolerance -> Bool
Eq, Int -> Tolerance -> ShowS
[Tolerance] -> ShowS
Tolerance -> String
(Int -> Tolerance -> ShowS)
-> (Tolerance -> String)
-> ([Tolerance] -> ShowS)
-> Show Tolerance
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Tolerance -> ShowS
showsPrec :: Int -> Tolerance -> ShowS
$cshow :: Tolerance -> String
show :: Tolerance -> String
$cshowList :: [Tolerance] -> ShowS
showList :: [Tolerance] -> ShowS
Show)

-- | Default adaptive tolerance: 1e-10 absolute, 1e-8 relative, 100000 evaluations.
defaultTolerance :: Tolerance
defaultTolerance :: Tolerance
defaultTolerance = Double -> Double -> Int -> Tolerance
Tolerance Double
1e-10 Double
1e-8 Int
100000

-- | One-sided or central first differences.
data Difference = Forward | Backward | Central deriving (Difference -> Difference -> Bool
(Difference -> Difference -> Bool)
-> (Difference -> Difference -> Bool) -> Eq Difference
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Difference -> Difference -> Bool
== :: Difference -> Difference -> Bool
$c/= :: Difference -> Difference -> Bool
/= :: Difference -> Difference -> Bool
Eq, Int -> Difference -> ShowS
[Difference] -> ShowS
Difference -> String
(Int -> Difference -> ShowS)
-> (Difference -> String)
-> ([Difference] -> ShowS)
-> Show Difference
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Difference -> ShowS
showsPrec :: Int -> Difference -> ShowS
$cshow :: Difference -> String
show :: Difference -> String
$cshowList :: [Difference] -> ShowS
showList :: [Difference] -> ShowS
Show)

{- | Composite Newton-Cotes formulas. Simpson panels must be a multiple of 2
and ThreeEighths panels a multiple of 3.
-}
data Quadrature = Trapezoidal | Simpson | ThreeEighths deriving (Quadrature -> Quadrature -> Bool
(Quadrature -> Quadrature -> Bool)
-> (Quadrature -> Quadrature -> Bool) -> Eq Quadrature
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Quadrature -> Quadrature -> Bool
== :: Quadrature -> Quadrature -> Bool
$c/= :: Quadrature -> Quadrature -> Bool
/= :: Quadrature -> Quadrature -> Bool
Eq, Int -> Quadrature -> ShowS
[Quadrature] -> ShowS
Quadrature -> String
(Int -> Quadrature -> ShowS)
-> (Quadrature -> String)
-> ([Quadrature] -> ShowS)
-> Show Quadrature
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Quadrature -> ShowS
showsPrec :: Int -> Quadrature -> ShowS
$cshow :: Quadrature -> String
show :: Quadrature -> String
$cshowList :: [Quadrature] -> ShowS
showList :: [Quadrature] -> ShowS
Show)

finite :: Double -> Bool
finite :: Double -> Bool
finite Double
x = Bool -> Bool
not (Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x)
check :: Double -> Either MathError Double
check :: Double -> Either MathError Double
check Double
x | Double -> Bool
finite Double
x = Double -> Either MathError Double
forall a b. b -> Either a b
Right Double
x | Bool
otherwise = MathError -> Either MathError Double
forall a b. a -> Either a b
Left MathError
NonFiniteResult
require :: Bool -> String -> Either MathError ()
require :: Bool -> String -> Either MathError ()
require Bool
p String
message = if Bool
p then () -> Either MathError ()
forall a b. b -> Either a b
Right () else MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
message)
validateTolerance :: Tolerance -> Either MathError ()
validateTolerance :: Tolerance -> Either MathError ()
validateTolerance (Tolerance Double
a Double
r Int
n) =
    Bool -> String -> Either MathError ()
require
        (Double -> Bool
finite Double
a Bool -> Bool -> Bool
&& Double -> Bool
finite Double
r Bool -> Bool -> Bool
&& Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
r Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0)
        String
"invalid numerical tolerance or evaluation budget"
target :: Tolerance -> Double -> Double
target :: Tolerance -> Double -> Double
target Tolerance
t Double
x = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max (Tolerance -> Double
absoluteTolerance Tolerance
t) (Tolerance -> Double
relativeTolerance Tolerance
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Num a => a -> a
abs Double
x)
sample :: (Double -> Double) -> Double -> Either MathError Double
sample :: (Double -> Double) -> Double -> Either MathError Double
sample Double -> Double
f Double
x = Double -> Either MathError Double
check Double
x Either MathError Double
-> Either MathError Double -> Either MathError Double
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Double -> Either MathError Double
check (Double -> Double
f Double
x)

-- | First derivative with explicit positive finite step; rounded-away steps fail.
derivative :: Difference -> Double -> (Double -> Double) -> Double -> Either MathError Double
derivative :: Difference
-> Double
-> (Double -> Double)
-> Double
-> Either MathError Double
derivative Difference
method Double
h Double -> Double
f Double
x = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
h Bool -> Bool -> Bool
&& Double
h Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
x) String
"invalid differentiation step or point"
    let (Double
a, Double
b) = case Difference
method of Difference
Forward -> (Double
x, Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h); Difference
Backward -> (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
h, Double
x); Difference
Central -> (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
h, Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h)
    Bool -> String -> Either MathError ()
require (Double
a Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
b) String
"unrepresentable differentiation step"
    fa <- (Double -> Double) -> Double -> Either MathError Double
sample Double -> Double
f Double
a
    fb <- sample f b
    check ((fb - fa) / (b - a))

{- | Arbitrary nonnegative derivative order by recursively centered differences.
Work grows exponentially with the requested order; orders are limited to 12.
-}
derivativeN :: Int -> Double -> (Double -> Double) -> Double -> Either MathError Double
derivativeN :: Int
-> Double
-> (Double -> Double)
-> Double
-> Either MathError Double
derivativeN Int
order Double
h Double -> Double
f Double
x = do
    Bool -> String -> Either MathError ()
require
        (Int
order Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
order Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
12 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
h Bool -> Bool -> Bool
&& Double
h Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
x)
        String
"invalid derivative order or step"
    let go :: t -> Double -> Either MathError Double
go t
0 Double
p = (Double -> Double) -> Double -> Either MathError Double
sample Double -> Double
f Double
p
        go t
n Double
p = do
            Bool -> String -> Either MathError ()
require (Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
p Bool -> Bool -> Bool
&& Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
h Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
p) String
"rounded-away derivative step"
            a <- t -> Double -> Either MathError Double
go (t
n t -> t -> t
forall a. Num a => a -> a -> a
- t
1) (Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
h)
            b <- go (n - 1) (p + h)
            check ((b - a) / (2 * h))
    Int -> Double -> Either MathError Double
forall {t}. (Eq t, Num t) => t -> Double -> Either MathError Double
go Int
order Double
x

{- | Complex-step first derivative. The callback must have a holomorphic
extension and be real on the real axis.
-}
complexDerivative ::
    Double -> (Complex Double -> Complex Double) -> Double -> Either MathError Double
complexDerivative :: Double
-> (Complex Double -> Complex Double)
-> Double
-> Either MathError Double
complexDerivative Double
h Complex Double -> Complex Double
f Double
x = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
h Bool -> Bool -> Bool
&& Double
h Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
x) String
"invalid complex-step arguments"
    let value :: Complex Double
value = Complex Double -> Complex Double
f (Double
x Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
h)
    _ <- Double -> Either MathError Double
check (Complex Double -> Double
forall a. Complex a -> a
realPart Complex Double
value)
    check (imagPart value / h)

{- | Richardson extrapolation of a central first difference. Each row halves
the step. Evaluation budget is 2 per row; convergence uses consecutive diagonals.
-}
richardson ::
    Tolerance -> Double -> (Double -> Double) -> Double -> Either MathError (Estimate Double)
richardson :: Tolerance
-> Double
-> (Double -> Double)
-> Double
-> Either MathError (Estimate Double)
richardson Tolerance
tol Double
initial Double -> Double
f Double
x = do
    Tolerance -> Either MathError ()
validateTolerance Tolerance
tol
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
initial Bool -> Bool -> Bool
&& Double
initial Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"invalid Richardson step"
    let go :: Int
-> [Double]
-> Double
-> Maybe Double
-> Either MathError (Estimate Double)
go Int
k [Double]
previous Double
h Maybe Double
old = do
            Bool -> String -> Either MathError ()
require (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Tolerance -> Int
maximumEvaluations Tolerance
tol) String
"Richardson budget must permit another row"
            base <- Difference
-> Double
-> (Double -> Double)
-> Double
-> Either MathError Double
derivative Difference
Central Double
h Double -> Double
f Double
x
            row <- build 1 [base] previous
            let value = [Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
row
                err = Double -> (Double -> Double) -> Maybe Double -> Double
forall b a. b -> (a -> b) -> Maybe a -> b
maybe (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0) (Double -> Double
forall a. Num a => a -> a
abs (Double -> Double) -> (Double -> Double) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double
value Double -> Double -> Double
forall a. Num a => a -> a -> a
-)) Maybe Double
old
            if err <= target tol value
                then Right (Estimate value err (2 * (k + 1)))
                else
                    if k >= 20 || 2 * (k + 2) > maximumEvaluations tol
                        then Left (NonConvergence "Richardson" (2 * (k + 1)))
                        else go (k + 1) row (h / 2) (Just value)
        build :: Int -> [Double] -> [Double] -> Either MathError [Double]
build Int
_ [Double]
acc [] = [Double] -> Either MathError [Double]
forall a b. b -> Either a b
Right [Double]
acc
        build Int
j [Double]
acc (Double
p : [Double]
ps) = do
            let current :: Double
current = [Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
acc
            value <- Double -> Either MathError Double
check (Double
current Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
current Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
p) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
4 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1))
            build (j + 1 :: Int) (acc ++ [value]) ps
    Int
-> [Double]
-> Double
-> Maybe Double
-> Either MathError (Estimate Double)
go (Int
0 :: Int) [] Double
initial Maybe Double
forall a. Maybe a
Nothing

-- | Ridders extrapolation with factor 1.4 and a best-error stopping criterion.
ridders :: Tolerance -> Double -> (Double -> Double) -> Double -> Either MathError (Estimate Double)
ridders :: Tolerance
-> Double
-> (Double -> Double)
-> Double
-> Either MathError (Estimate Double)
ridders Tolerance
tol Double
initial Double -> Double
f Double
x = do
    Tolerance -> Either MathError ()
validateTolerance Tolerance
tol
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
initial Bool -> Bool -> Bool
&& Double
initial Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"invalid Ridders step"
    let go :: Int
-> [Double]
-> Double
-> (Double, Double)
-> Either MathError (Estimate Double)
go Int
k [Double]
previous Double
h (Double, Double)
best = do
            d <- Difference
-> Double
-> (Double -> Double)
-> Double
-> Either MathError Double
derivative Difference
Central Double
h Double -> Double
f Double
x
            row <- build 1 [d] previous
            let candidates = (Double -> Double -> Double -> (Double, Double))
-> [Double] -> [Double] -> [Double] -> [(Double, Double)]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 (\Double
v Double
a Double
b -> (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max (Double -> Double
forall a. Num a => a -> a
abs (Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a)) (Double -> Double
forall a. Num a => a -> a
abs (Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b)), Double
v)) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
row) [Double]
row [Double]
previous
                best' = ((Double, Double) -> (Double, Double) -> (Double, Double))
-> (Double, Double) -> [(Double, Double)] -> (Double, Double)
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (\(Double, Double)
p (Double, Double)
q -> if (Double, Double) -> Double
forall a b. (a, b) -> a
fst (Double, Double)
q Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< (Double, Double) -> Double
forall a b. (a, b) -> a
fst (Double, Double)
p then (Double, Double)
q else (Double, Double)
p) (Double, Double)
best [(Double, Double)]
candidates
                (err, value) = best'
            if err <= target tol value
                then Right (Estimate value err (2 * (k + 1)))
                else
                    if k >= 20 || 2 * (k + 2) > maximumEvaluations tol
                        then Left (NonConvergence "Ridders" (2 * (k + 1)))
                        else go (k + 1) row (h / 1.4) best'
        build :: Int -> [Double] -> [Double] -> Either MathError [Double]
build Int
_ [Double]
acc [] = [Double] -> Either MathError [Double]
forall a b. b -> Either a b
Right [Double]
acc
        build Int
j [Double]
acc (Double
p : [Double]
ps) = do
            let current :: Double
current = [Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
acc; factor :: Double
factor = Double
1.4 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
j)
            value <- Double -> Either MathError Double
check ((Double
factor Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
current Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
p) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
factor Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1))
            build (j + 1 :: Int) (acc ++ [value]) ps
    Bool -> String -> Either MathError ()
require (Tolerance -> Int
maximumEvaluations Tolerance
tol Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
4) String
"Ridders needs at least four evaluations"
    Int
-> [Double]
-> Double
-> (Double, Double)
-> Either MathError (Estimate Double)
go (Int
0 :: Int) [] Double
initial (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0, Double
0)

{- | One-sided numerical limit, approaching geometrically from the requested
side. Central compares the average of both sides; it cannot prove continuity.
-}
limit ::
    Difference ->
    Tolerance ->
    Double ->
    (Double -> Double) ->
    Double ->
    Either MathError (Estimate Double)
limit :: Difference
-> Tolerance
-> Double
-> (Double -> Double)
-> Double
-> Either MathError (Estimate Double)
limit Difference
side Tolerance
tol Double
initial Double -> Double
f Double
x = do
    Tolerance -> Either MathError ()
validateTolerance Tolerance
tol
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
initial Bool -> Bool -> Bool
&& Double
initial Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
x) String
"invalid limit arguments"
    let cost :: Int
cost = if Difference
side Difference -> Difference -> Bool
forall a. Eq a => a -> a -> Bool
== Difference
Central then Int
2 else Int
1
        go :: Int -> Double -> Maybe Double -> Either MathError (Estimate Double)
go Int
n Double
h Maybe Double
old = do
            Bool -> String -> Either MathError ()
require (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
x Bool -> Bool -> Bool
&& Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
h Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
x) String
"limit reached floating-point resolution"
            value <- case Difference
side of
                Difference
Forward -> (Double -> Double) -> Double -> Either MathError Double
sample Double -> Double
f (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h)
                Difference
Backward -> (Double -> Double) -> Double -> Either MathError Double
sample Double -> Double
f (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
h)
                Difference
Central -> do a <- (Double -> Double) -> Double -> Either MathError Double
sample Double -> Double
f (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
h); b <- sample f (x + h); check ((a + b) / 2)
            let err = Double -> (Double -> Double) -> Maybe Double -> Double
forall b a. b -> (a -> b) -> Maybe a -> b
maybe (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0) (Double -> Double
forall a. Num a => a -> a
abs (Double -> Double) -> (Double -> Double) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double
value Double -> Double -> Double
forall a. Num a => a -> a -> a
-)) Maybe Double
old
            if err <= target tol value
                then Right (Estimate value err (n + cost))
                else
                    if n + 2 * cost > maximumEvaluations tol
                        then Left (NonConvergence "limit" (n + cost))
                        else go (n + cost) (h / 2) (Just value)
    Bool -> String -> Either MathError ()
require (Tolerance -> Int
maximumEvaluations Tolerance
tol Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
cost) String
"limit budget is too small for one sample"
    Int -> Double -> Maybe Double -> Either MathError (Estimate Double)
go Int
0 Double
initial Maybe Double
forall a. Maybe a
Nothing

replace :: Int -> a -> [a] -> [a]
replace :: forall a. Int -> a -> [a] -> [a]
replace Int
i a
x = (Int -> a -> a) -> [Int] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Int
j a
v -> if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i then a
x else a
v) [Int
0 ..]

-- | Partial derivative at a zero-based coordinate.
partialDerivative :: Int -> Double -> ([Double] -> Double) -> [Double] -> Either MathError Double
partialDerivative :: Int
-> Double
-> ([Double] -> Double)
-> [Double]
-> Either MathError Double
partialDerivative Int
i Double
h [Double] -> Double
f [Double]
xs = do
    Bool -> String -> Either MathError ()
require (Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
xs) String
"invalid derivative coordinate or point"
    Difference
-> Double
-> (Double -> Double)
-> Double
-> Either MathError Double
derivative Difference
Central Double
h (\Double
x -> [Double] -> Double
f (Int -> Double -> [Double] -> [Double]
forall a. Int -> a -> [a] -> [a]
replace Int
i Double
x [Double]
xs)) ([Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)

-- | Gradient using one central difference per coordinate.
gradient :: Double -> ([Double] -> Double) -> [Double] -> Either MathError [Double]
gradient :: Double
-> ([Double] -> Double) -> [Double] -> Either MathError [Double]
gradient Double
h [Double] -> Double
f [Double]
xs = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
h Bool -> Bool -> Bool
&& Double
h Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
xs) String
"invalid gradient step or point"
    (Int -> Either MathError Double)
-> [Int] -> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse (\Int
i -> Int
-> Double
-> ([Double] -> Double)
-> [Double]
-> Either MathError Double
partialDerivative Int
i Double
h [Double] -> Double
f [Double]
xs) [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]

-- | Row-major Jacobian; callback output dimensions must remain constant.
jacobian :: Double -> ([Double] -> [Double]) -> [Double] -> Either MathError (Matrix Double)
jacobian :: Double
-> ([Double] -> [Double])
-> [Double]
-> Either MathError (Matrix Double)
jacobian Double
h [Double] -> [Double]
f [Double]
xs = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
h Bool -> Bool -> Bool
&& Double
h Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
xs) String
"invalid Jacobian step or point"
    let ys :: [Double]
ys = [Double] -> [Double]
f [Double]
xs; m :: Int
m = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys; n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs
    Bool -> String -> Either MathError ()
require ((Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
xs Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
ys) String
"nonfinite Jacobian point or result"
    columns <- (Int -> Either MathError [Double])
-> [Int] -> Either MathError [[Double]]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse (Int -> Int -> Either MathError [Double]
column Int
m) [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
    fromListChecked m n (concat (transpose columns))
  where
    column :: Int -> Int -> Either MathError [Double]
column Int
m Int
i = do
        Bool -> String -> Either MathError ()
require
            (Double -> Bool
finite Double
h Bool -> Bool -> Bool
&& Double
h Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Bool -> Bool -> Bool
&& [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
h Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)
            String
"invalid Jacobian step"
        left <- (Double -> Either MathError Double)
-> [Double] -> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse Double -> Either MathError Double
check (Int -> Double -> [Double] -> [Double]
forall a. Int -> a -> [a] -> [a]
replace Int
i ([Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
h) [Double]
xs)
        right <- traverse check (replace i (xs !! i + h) xs)
        let a = [Double] -> [Double]
f [Double]
left; b = [Double] -> [Double]
f [Double]
right
        require (length a == m && length b == m) "callback changed Jacobian output dimension"
        _ <- traverse check (a ++ b)
        traverse check (zipWith (\Double
u Double
v -> (Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
u) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
h)) a b)

-- | Hessian with centered diagonal and mixed differences.
hessian :: Double -> ([Double] -> Double) -> [Double] -> Either MathError (Matrix Double)
hessian :: Double
-> ([Double] -> Double)
-> [Double]
-> Either MathError (Matrix Double)
hessian Double
h [Double] -> Double
f [Double]
xs = do
    Bool -> String -> Either MathError ()
require
        (Double -> Bool
finite Double
h Bool -> Bool -> Bool
&& Double
h Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
xs Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Double
x -> Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
x Bool -> Bool -> Bool
&& Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
h Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
x) [Double]
xs)
        String
"invalid Hessian step or point"
    f0 <- Double -> Either MathError Double
check ([Double] -> Double
f [Double]
xs)
    values <- traverse (entry f0) [(i, j) | i <- [0 .. n - 1], j <- [0 .. n - 1]]
    fromListChecked n n values
  where
    n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs
    point :: Int -> Double -> Int -> Double -> [Double]
point Int
i Double
d Int
j Double
e = Int -> Double -> [Double] -> [Double]
forall a. Int -> a -> [a] -> [a]
replace Int
i ([Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d) (Int -> Double -> [Double] -> [Double]
forall a. Int -> a -> [a] -> [a]
replace Int
j ([Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
e) [Double]
xs)
    evaluatePoint :: [Double] -> Either MathError Double
evaluatePoint [Double]
p = (Double -> Either MathError Double)
-> [Double] -> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse Double -> Either MathError Double
check [Double]
p Either MathError [Double]
-> Either MathError Double -> Either MathError Double
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Double -> Either MathError Double
check ([Double] -> Double
f [Double]
p)
    entry :: Double -> (Int, Int) -> Either MathError Double
entry Double
f0 (Int
i, Int
j)
        | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j = do
            a <- [Double] -> Either MathError Double
evaluatePoint (Int -> Double -> [Double] -> [Double]
forall a. Int -> a -> [a] -> [a]
replace Int
i ([Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h) [Double]
xs)
            b <- evaluatePoint (replace i (xs !! i - h) xs)
            check ((a - 2 * f0 + b) / (h * h))
        | Bool
otherwise = do
            a <- [Double] -> Either MathError Double
evaluatePoint (Int -> Double -> Int -> Double -> [Double]
point Int
i Double
h Int
j Double
h)
            b <- evaluatePoint (point i h j (-h))
            c <- evaluatePoint (point i (-h) j h)
            d <- evaluatePoint (point i (-h) j (-h))
            check ((a - b - c + d) / (4 * h * h))

-- | Directional derivative; the direction is used as supplied, without normalization.
directionalDerivative ::
    Double -> ([Double] -> Double) -> [Double] -> [Double] -> Either MathError Double
directionalDerivative :: Double
-> ([Double] -> Double)
-> [Double]
-> [Double]
-> Either MathError Double
directionalDerivative Double
h [Double] -> Double
f [Double]
x [Double]
v = do
    Bool -> String -> Either MathError ()
require
        ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
v Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite ([Double]
x [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
v))
        String
"direction dimension mismatch or nonfinite direction"
    Difference
-> Double
-> (Double -> Double)
-> Double
-> Either MathError Double
derivative
        Difference
Central
        Double
h
        (\Double
t -> let p :: [Double]
p = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
a Double
b -> Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b) [Double]
x [Double]
v in if (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
p then [Double] -> Double
f [Double]
p else Double
0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0)
        Double
0

-- | Jacobian-vector product by directional central differences.
jvp :: Double -> ([Double] -> [Double]) -> [Double] -> [Double] -> Either MathError [Double]
jvp :: Double
-> ([Double] -> [Double])
-> [Double]
-> [Double]
-> Either MathError [Double]
jvp Double
h [Double] -> [Double]
f [Double]
x [Double]
v = do
    Bool -> String -> Either MathError ()
require ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
v Bool -> Bool -> Bool
&& Double -> Bool
finite Double
h Bool -> Bool -> Bool
&& Double
h Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite ([Double]
x [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
v)) String
"invalid JVP arguments"
    left <- (Double -> Either MathError Double)
-> [Double] -> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse Double -> Either MathError Double
check ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
u Double
w -> Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w) [Double]
x [Double]
v)
    right <- traverse check (zipWith (\Double
u Double
w -> Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w) x v)
    let a = [Double] -> [Double]
f [Double]
left; b = [Double] -> [Double]
f [Double]
right
    require (length a == length b && length a == length (f x)) "callback changed JVP output dimension"
    _ <- traverse check (a ++ b)
    traverse check (zipWith (\Double
u Double
w -> (Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
u) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
h)) a b)

-- | Vector-Jacobian product as the gradient of a scalar contraction.
vjp :: Double -> ([Double] -> [Double]) -> [Double] -> [Double] -> Either MathError [Double]
vjp :: Double
-> ([Double] -> [Double])
-> [Double]
-> [Double]
-> Either MathError [Double]
vjp Double
h [Double] -> [Double]
f [Double]
x [Double]
v = do
    Bool -> String -> Either MathError ()
require ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Double] -> [Double]
f [Double]
x) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
v Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
v) String
"VJP output dimension mismatch"
    let contracted :: [Double] -> Double
contracted [Double]
p = let y :: [Double]
y = [Double] -> [Double]
f [Double]
p in if [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
y Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
v then Double
0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0 else [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) [Double]
v [Double]
y)
    Double
-> ([Double] -> Double) -> [Double] -> Either MathError [Double]
gradient Double
h [Double] -> Double
contracted [Double]
x

-- | Hessian-vector product with four directional/coordinate samples per component.
hvp :: Double -> ([Double] -> Double) -> [Double] -> [Double] -> Either MathError [Double]
hvp :: Double
-> ([Double] -> Double)
-> [Double]
-> [Double]
-> Either MathError [Double]
hvp Double
h [Double] -> Double
f [Double]
x [Double]
v = do
    Bool -> String -> Either MathError ()
require ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
v Bool -> Bool -> Bool
&& Double -> Bool
finite Double
h Bool -> Bool -> Bool
&& Double
h Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite ([Double]
x [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
v)) String
"invalid HVP arguments"
    let at :: Double -> [Double]
at Double
sign = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
a Double
b -> Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
sign Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b) [Double]
x [Double]
v
    a <- Double
-> ([Double] -> Double) -> [Double] -> Either MathError [Double]
gradient Double
h [Double] -> Double
f (Double -> [Double]
at Double
1)
    b <- gradient h f (at (-1))
    traverse check (zipWith (\Double
u Double
w -> (Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
w) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
h)) a b)

-- | Divergence of a vector field with equal input and output dimensions.
divergence :: Double -> ([Double] -> [Double]) -> [Double] -> Either MathError Double
divergence :: Double
-> ([Double] -> [Double]) -> [Double] -> Either MathError Double
divergence Double
h [Double] -> [Double]
f [Double]
x = do
    Bool -> String -> Either MathError ()
require
        (Double -> Bool
finite Double
h Bool -> Bool -> Bool
&& Double
h Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
x Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Double] -> [Double]
f [Double]
x) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite ([Double] -> [Double]
f [Double]
x))
        String
"divergence requires finite matching dimensions and a positive step"
    values <-
        (Int -> Either MathError Double)
-> [Int] -> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse
            ( \Int
i ->
                Int
-> Double
-> ([Double] -> Double)
-> [Double]
-> Either MathError Double
partialDerivative
                    Int
i
                    Double
h
                    (\[Double]
p -> let y :: [Double]
y = [Double] -> [Double]
f [Double]
p in if [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
y Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
y then [Double]
y [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i else Double
0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0)
                    [Double]
x
            )
            [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
    check (sum values)

-- | Curl of a three-dimensional vector field.
curl :: Double -> ([Double] -> [Double]) -> [Double] -> Either MathError [Double]
curl :: Double
-> ([Double] -> [Double]) -> [Double] -> Either MathError [Double]
curl Double
h [Double] -> [Double]
f [Double]
x = do
    Bool -> String -> Either MathError ()
require ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
3 Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Double] -> [Double]
f [Double]
x) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
3) String
"curl requires three dimensions"
    let part :: Int -> Int -> Either MathError Double
part Int
i Int
j =
            Int
-> Double
-> ([Double] -> Double)
-> [Double]
-> Either MathError Double
partialDerivative
                Int
j
                Double
h
                (\[Double]
p -> let y :: [Double]
y = [Double] -> [Double]
f [Double]
p in if [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
y Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
3 Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
y then [Double]
y [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i else Double
0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0)
                [Double]
x
    a <- Int -> Int -> Either MathError Double
part Int
2 Int
1
    b <- part 1 2
    c <- part 0 2
    d <- part 2 0
    e <- part 1 0
    g <- part 0 1
    traverse check [a - b, c - d, e - g]

-- | Laplacian as the sum of centered second coordinate differences.
laplacian :: Double -> ([Double] -> Double) -> [Double] -> Either MathError Double
laplacian :: Double
-> ([Double] -> Double) -> [Double] -> Either MathError Double
laplacian Double
h [Double] -> Double
f [Double]
x = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
h Bool -> Bool -> Bool
&& Double
h Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
x) String
"invalid Laplacian arguments"
    values <- (Int -> Either MathError Double)
-> [Int] -> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse (\Int
i -> Int
-> Double
-> (Double -> Double)
-> Double
-> Either MathError Double
derivativeN Int
2 Double
h (\Double
p -> [Double] -> Double
f (Int -> Double -> [Double] -> [Double]
forall a. Int -> a -> [a] -> [a]
replace Int
i Double
p [Double]
x)) ([Double]
x [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)) [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
    check (sum values)

-- | Composite finite-interval integration with exact panel constraints.
integrateComposite ::
    Quadrature -> Int -> (Double -> Double) -> Double -> Double -> Either MathError Double
integrateComposite :: Quadrature
-> Int
-> (Double -> Double)
-> Double
-> Double
-> Either MathError Double
integrateComposite Quadrature
method Int
n Double -> Double
f Double
a Double
b = do
    let divisor :: Int
divisor = case Quadrature
method of Quadrature
Trapezoidal -> Int
1; Quadrature
Simpson -> Int
2; Quadrature
ThreeEighths -> Int
3
    Bool -> String -> Either MathError ()
require
        (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
divisor Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
a Bool -> Bool -> Bool
&& Double -> Bool
finite Double
b)
        String
"invalid integration panels or endpoints"
    if Double
a Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
b
        then Double -> Either MathError Double
forall a b. b -> Either a b
Right Double
0
        else do
            let h :: Double
h = (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
                weight :: Int -> a
weight Int
i = case Quadrature
method of
                    Quadrature
Trapezoidal -> if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 Bool -> Bool -> Bool
|| Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n then a
0.5 else a
1
                    Quadrature
Simpson -> if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 Bool -> Bool -> Bool
|| Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n then a
1 else if Int -> Bool
forall a. Integral a => a -> Bool
odd Int
i then a
4 else a
2
                    Quadrature
ThreeEighths -> if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 Bool -> Bool -> Bool
|| Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n then a
1 else if Int
i Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
3 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then a
2 else a
3
                factor :: Double
factor = case Quadrature
method of Quadrature
Trapezoidal -> Double
1; Quadrature
Simpson -> Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3; Quadrature
ThreeEighths -> Double
3 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
8
            values <- (Int -> Either MathError Double)
-> [Int] -> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse (\Int
i -> (Int -> Double
forall {a}. Fractional a => Int -> a
weight Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
*) (Double -> Double)
-> Either MathError Double -> Either MathError Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> (Double -> Double) -> Double -> Either MathError Double
sample Double -> Double
f (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
h)) [Int
0 .. Int
n]
            check (h * factor * sum values)

-- | Romberg integration using trapezoidal refinement and Richardson diagonals.
romberg :: Tolerance -> (Double -> Double) -> Double -> Double -> Either MathError (Estimate Double)
romberg :: Tolerance
-> (Double -> Double)
-> Double
-> Double
-> Either MathError (Estimate Double)
romberg Tolerance
tol Double -> Double
f Double
a Double
b = do
    Tolerance -> Either MathError ()
validateTolerance Tolerance
tol
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
a Bool -> Bool -> Bool
&& Double -> Bool
finite Double
b) String
"nonfinite integration endpoint"
    if Double
a Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
b
        then Estimate Double -> Either MathError (Estimate Double)
forall a b. b -> Either a b
Right (Double -> Double -> Int -> Estimate Double
forall a. a -> Double -> Int -> Estimate a
Estimate Double
0 Double
0 Int
0)
        else do
            Bool -> String -> Either MathError ()
require (Tolerance -> Int
maximumEvaluations Tolerance
tol Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
3) String
"Romberg needs at least three evaluations"
            fa <- (Double -> Double) -> Double -> Either MathError Double
sample Double -> Double
f Double
a
            fb <- sample f b
            let go t
level [Double]
previous Int
count = do
                    let n :: Int
n = Int
2 Int -> t -> Int
forall a b. (Num a, Integral b) => a -> b -> a
^ t
level; newCount :: Int
newCount = Int
n Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
                    if Int
count Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
newCount Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Tolerance -> Int
maximumEvaluations Tolerance
tol
                        then MathError -> Either MathError (Estimate Double)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"Romberg" Int
count)
                        else do
                            values <-
                                (Int -> Either MathError Double)
-> [Int] -> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse (\Int
i -> (Double -> Double) -> Double -> Either MathError Double
sample Double -> Double
f (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)) [Int
1 .. Int
newCount]
                            let first = Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
previous [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
values
                                row =
                                    (Double -> (Int, Double) -> Double)
-> Double -> [(Int, Double)] -> [Double]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl
                                        (\Double
current (Int
j, Double
p) -> Double
current Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
current Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
p) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
4 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1))
                                        Double
first
                                        ([Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
1 :: Int ..] [Double]
previous)
                                value = [Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
row
                                err = Double -> Double
forall a. Num a => a -> a
abs (Double
value Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
previous)
                            _ <- check value
                            if err <= target tol value
                                then Right (Estimate value err (count + newCount))
                                else
                                    if level >= 20
                                        then Left (NonConvergence "Romberg" (count + newCount))
                                        else go (level + 1) row (count + newCount)
            go (1 :: Int) [(b - a) * (fa + fb) / 2] 2

-- | n-point Gauss-Legendre quadrature; nodes are refined by Newton iteration.
gaussLegendre :: Int -> (Double -> Double) -> Double -> Double -> Either MathError Double
gaussLegendre :: Int
-> (Double -> Double)
-> Double
-> Double
-> Either MathError Double
gaussLegendre Int
n Double -> Double
f Double
a Double
b = do
    Bool -> String -> Either MathError ()
require (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
512 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
a Bool -> Bool -> Bool
&& Double -> Bool
finite Double
b) String
"invalid Gauss-Legendre order or bounds"
    nodes <- (Int -> Either MathError (Double, Double))
-> [Int] -> Either MathError [(Double, Double)]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse Int -> Either MathError (Double, Double)
forall {t} {a}.
(Floating t, Integral a, Ord t) =>
a -> Either MathError (t, t)
root [Int
1 .. Int
n]
    values <- traverse (\(Double
x, Double
w) -> (Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
*) (Double -> Double)
-> Either MathError Double -> Either MathError Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> (Double -> Double) -> Double -> Either MathError Double
sample Double -> Double
f ((Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x)) nodes
    check ((b - a) / 2 * sum values)
  where
    poly :: a -> (a, a)
poly a
x =
        ((a, a) -> Int -> (a, a)) -> (a, a) -> [Int] -> (a, a)
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl
            (\(a
p, a
q) Int
k -> (((a
2 a -> a -> a
forall a. Num a => a -> a -> a
* Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k a -> a -> a
forall a. Num a => a -> a -> a
- a
1) a -> a -> a
forall a. Num a => a -> a -> a
* a
x a -> a -> a
forall a. Num a => a -> a -> a
* a
p a -> a -> a
forall a. Num a => a -> a -> a
- (Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k a -> a -> a
forall a. Num a => a -> a -> a
- a
1) a -> a -> a
forall a. Num a => a -> a -> a
* a
q) a -> a -> a
forall a. Fractional a => a -> a -> a
/ Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k, a
p))
            (a
1, a
0)
            [Int
1 .. Int
n]
    root :: a -> Either MathError (t, t)
root a
i = Int -> t -> Either MathError (t, t)
forall {t}.
(Ord t, Fractional t) =>
Int -> t -> Either MathError (t, t)
go (Int
0 :: Int) (t -> t
forall a. Floating a => a -> a
cos (t
forall a. Floating a => a
pi t -> t -> t
forall a. Num a => a -> a -> a
* (a -> t
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
i t -> t -> t
forall a. Num a => a -> a -> a
- t
0.25) t -> t -> t
forall a. Fractional a => a -> a -> a
/ (Int -> t
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n t -> t -> t
forall a. Num a => a -> a -> a
+ t
0.5)))
    go :: Int -> t -> Either MathError (t, t)
go Int
k t
x =
        let (t
p, t
q) = t -> (t, t)
forall {a}. Fractional a => a -> (a, a)
poly t
x; d :: t
d = Int -> t
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n t -> t -> t
forall a. Num a => a -> a -> a
* (t
x t -> t -> t
forall a. Num a => a -> a -> a
* t
p t -> t -> t
forall a. Num a => a -> a -> a
- t
q) t -> t -> t
forall a. Fractional a => a -> a -> a
/ (t
x t -> t -> t
forall a. Num a => a -> a -> a
* t
x t -> t -> t
forall a. Num a => a -> a -> a
- t
1); next :: t
next = t
x t -> t -> t
forall a. Num a => a -> a -> a
- t
p t -> t -> t
forall a. Fractional a => a -> a -> a
/ t
d
         in if Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
100
                then MathError -> Either MathError (t, t)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"Gauss-Legendre node" Int
k)
                else
                    if t -> t
forall a. Num a => a -> a
abs (t
next t -> t -> t
forall a. Num a => a -> a -> a
- t
x) t -> t -> Bool
forall a. Ord a => a -> a -> Bool
<= t
4e-16 then (t, t) -> Either MathError (t, t)
forall a b. b -> Either a b
Right (t
next, t
2 t -> t -> t
forall a. Fractional a => a -> a -> a
/ ((t
1 t -> t -> t
forall a. Num a => a -> a -> a
- t
next t -> t -> t
forall a. Num a => a -> a -> a
* t
next) t -> t -> t
forall a. Num a => a -> a -> a
* t
d t -> t -> t
forall a. Num a => a -> a -> a
* t
d)) else Int -> t -> Either MathError (t, t)
go (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) t
next

{- | Weighted Gaussian quadrature on the real line for weight exp(-x^2).
Uses the symmetric Jacobi recurrence and a checked implicit QL eigensystem.
-}
gaussHermite :: Int -> (Double -> Double) -> Either MathError Double
gaussHermite :: Int -> (Double -> Double) -> Either MathError Double
gaussHermite Int
n = Int
-> [Double]
-> [Double]
-> Double
-> (Double -> Double)
-> Either MathError Double
gaussianRule Int
n (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
n Double
0) [Double -> Double
forall a. Floating a => a -> a
sqrt (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) | Int
k <- [Int
1 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] (Double -> Double
forall a. Floating a => a -> a
sqrt Double
forall a. Floating a => a
pi)

-- | Weighted Gaussian quadrature on [0,infinity) for weight exp(-x).
gaussLaguerre :: Int -> (Double -> Double) -> Either MathError Double
gaussLaguerre :: Int -> (Double -> Double) -> Either MathError Double
gaussLaguerre Int
n = Int
-> [Double]
-> [Double]
-> Double
-> (Double -> Double)
-> Either MathError Double
gaussianRule Int
n [Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1 | Int
k <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] ((Int -> Double) -> [Int] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral [Int
1 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]) Double
1

-- | Weighted first-kind Chebyshev rule on [-1,1], weight 1/sqrt(1-x^2).
gaussChebyshev :: Int -> (Double -> Double) -> Either MathError Double
gaussChebyshev :: Int -> (Double -> Double) -> Either MathError Double
gaussChebyshev Int
n Double -> Double
f = do
    Bool -> String -> Either MathError ()
require (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
100000) String
"invalid Gauss-Chebyshev order"
    values <- (Int -> Either MathError Double)
-> [Int] -> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse ((Double -> Double) -> Double -> Either MathError Double
sample Double -> Double
f (Double -> Either MathError Double)
-> (Int -> Double) -> Int -> Either MathError Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
forall a. Floating a => a -> a
cos (Double -> Double) -> (Int -> Double) -> Int -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (\Int
i -> Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.5) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)) [Int
1 .. Int
n]
    check (pi / fromIntegral n * sum values)

-- Small symmetric Jacobi eigensolver for quadrature construction; largest
-- off-diagonal rotations preserve eigenvectors and terminate by a scale check.
gaussianRule ::
    Int -> [Double] -> [Double] -> Double -> (Double -> Double) -> Either MathError Double
gaussianRule :: Int
-> [Double]
-> [Double]
-> Double
-> (Double -> Double)
-> Either MathError Double
gaussianRule Int
n [Double]
diagonal [Double]
off Double
mass Double -> Double
f = do
    Bool -> String -> Either MathError ()
require (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
64) String
"weighted Gaussian order must lie in [1,64]"
    let a :: [[Double]]
a =
            [ [ if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then [Double]
diagonal [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i else if Int -> Int
forall a. Num a => a -> a
abs (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
j) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 then [Double]
off [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
i Int
j else Double
0 | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
              ]
            | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
            ]
        unit :: [[Double]]
unit = [[if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then Double
1 else Double
0 | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        go :: Int
-> [[Double]]
-> [[Double]]
-> Either MathError ([[Double]], [[Double]])
go Int
k [[Double]]
matrix [[Double]]
vectors =
            let candidates :: [(Double, Int, Int)]
candidates = [(Double -> Double
forall a. Num a => a -> a
abs ([[Double]]
matrix [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j), Int
i, Int
j) | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
j <- [Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
                (Double
largest, Int
p, Int
q) = if [(Double, Int, Int)] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(Double, Int, Int)]
candidates then (Double
0, Int
0, Int
0) else [(Double, Int, Int)] -> (Double, Int, Int)
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [(Double, Int, Int)]
candidates
                scale :: Double
scale = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Double
1 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
abs [Double]
diagonal)
             in if Double
largest Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1e-14 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
scale
                    then ([[Double]], [[Double]])
-> Either MathError ([[Double]], [[Double]])
forall a b. b -> Either a b
Right ([[Double]]
matrix, [[Double]]
vectors)
                    else
                        if Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
100 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n
                            then MathError -> Either MathError ([[Double]], [[Double]])
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"Gaussian quadrature eigensystem" Int
k)
                            else
                                let app :: Double
app = [[Double]]
matrix [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
p [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
p
                                    aqq :: Double
aqq = [[Double]]
matrix [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
q [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
q
                                    apq :: Double
apq = [[Double]]
matrix [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
p [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
q
                                    theta :: Double
theta = Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
apq) (Double
aqq Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
app)
                                    c :: Double
c = Double -> Double
forall a. Floating a => a -> a
cos Double
theta
                                    s :: Double
s = Double -> Double
forall a. Floating a => a -> a
sin Double
theta
                                    rotate :: [Double] -> [Double]
rotate [Double]
row = Int -> Double -> [Double] -> [Double]
forall a. Int -> a -> [a] -> [a]
replace Int
q (Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
p) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
q)) (Int -> Double -> [Double] -> [Double]
forall a. Int -> a -> [a] -> [a]
replace Int
p (Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
p) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
q)) [Double]
row)
                                    columns :: [[Double]]
columns = ([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map [Double] -> [Double]
rotate [[Double]]
matrix
                                    matrix' :: [[Double]]
matrix' = [[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
transpose (([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map [Double] -> [Double]
rotate ([[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
transpose [[Double]]
columns))
                                    vectors' :: [[Double]]
vectors' = ([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map [Double] -> [Double]
rotate [[Double]]
vectors
                                 in Int
-> [[Double]]
-> [[Double]]
-> Either MathError ([[Double]], [[Double]])
go (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [[Double]]
matrix' [[Double]]
vectors'
    (eigen, vectors) <- Int
-> [[Double]]
-> [[Double]]
-> Either MathError ([[Double]], [[Double]])
go (Int
0 :: Int) [[Double]]
a [[Double]]
unit
    values <-
        traverse
            (\Int
i -> do y <- (Double -> Double) -> Double -> Either MathError Double
sample Double -> Double
f ([[Double]]
eigen [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i); check (mass * (vectors !! 0 !! i) ^ (2 :: Int) * y))
            [0 .. n - 1]
    check (sum values)

{- | Embedded 7-point Gauss / 15-point Kronrod estimate. Nodes and weights are
the standard Gauss-Kronrod pair used by QUADPACK dqk15 (Netlib).
-}
gaussKronrod :: (Double -> Double) -> Double -> Double -> Either MathError (Estimate Double)
gaussKronrod :: (Double -> Double)
-> Double -> Double -> Either MathError (Estimate Double)
gaussKronrod Double -> Double
f Double
a Double
b = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
a Bool -> Bool -> Bool
&& Double -> Bool
finite Double
b) String
"nonfinite quadrature bounds"
    if Double
a Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
b
        then Estimate Double -> Either MathError (Estimate Double)
forall a b. b -> Either a b
Right (Double -> Double -> Int -> Estimate Double
forall a. a -> Double -> Int -> Estimate a
Estimate Double
0 Double
0 Int
0)
        else do
            center <- (Double -> Double) -> Double -> Either MathError Double
sample Double -> Double
f ((Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)
            pairs <-
                traverse
                    ( \Double
x -> do
                        u <- (Double -> Double) -> Double -> Either MathError Double
sample Double -> Double
f ((Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x)
                        v <- sample f ((a + b) / 2 + (b - a) / 2 * x)
                        check (u + v)
                    )
                    nodes
            let kronrod = (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
wk0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
center Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) [Double]
weights [Double]
pairs))
                gauss = (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
wg0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
center Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) [Double]
gaussWeights [Double]
pairs))
            value <- check kronrod
            err <- check (abs (kronrod - gauss))
            Right (Estimate value err 15)
  where
    nodes :: [Double]
nodes =
        [ Double
0.9914553711208126
        , Double
0.9491079123427585
        , Double
0.8648644233597691
        , Double
0.7415311855993944
        , Double
0.5860872354676911
        , Double
0.4058451513773972
        , Double
0.2077849550078985
        ]
    weights :: [Double]
weights =
        [ Double
0.02293532201052922
        , Double
0.06309209262997855
        , Double
0.1047900103222502
        , Double
0.1406532597155259
        , Double
0.1690047266392679
        , Double
0.1903505780647854
        , Double
0.2044329400752989
        ]
    gaussWeights :: [Double]
gaussWeights = [Double
0, Double
0.1294849661688697, Double
0, Double
0.2797053914892767, Double
0, Double
0.3818300505051189, Double
0]
    wk0 :: Double
wk0 = Double
0.2094821410847278
    wg0 :: Double
wg0 = Double
0.4179591836734694

{- | Adaptive Gauss-Kronrod with global summed error and a strict evaluation
budget. The interval with the largest estimated error is bisected first.
-}
adaptiveIntegrate ::
    Tolerance -> (Double -> Double) -> Double -> Double -> Either MathError (Estimate Double)
adaptiveIntegrate :: Tolerance
-> (Double -> Double)
-> Double
-> Double
-> Either MathError (Estimate Double)
adaptiveIntegrate Tolerance
tol Double -> Double
f Double
a Double
b = do
    Tolerance -> Either MathError ()
validateTolerance Tolerance
tol
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
a Bool -> Bool -> Bool
&& Double -> Bool
finite Double
b) String
"adaptive bounds must be finite"
    if Double
a Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
b
        then Estimate Double -> Either MathError (Estimate Double)
forall a b. b -> Either a b
Right (Double -> Double -> Int -> Estimate Double
forall a. a -> Double -> Int -> Estimate a
Estimate Double
0 Double
0 Int
0)
        else do
            Bool -> String -> Either MathError ()
require (Tolerance -> Int
maximumEvaluations Tolerance
tol Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
15) String
"Gauss-Kronrod needs fifteen evaluations"
            first <- (Double -> Double)
-> Double -> Double -> Either MathError (Estimate Double)
gaussKronrod Double -> Double
f Double
a Double
b
            let go [(Double, Double, Estimate Double)]
panels Int
count = do
                    let value :: Double
value = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Estimate Double -> Double
forall a. Estimate a -> a
estimateValue Estimate Double
e | (Double
_, Double
_, Estimate Double
e) <- [(Double, Double, Estimate Double)]
panels]
                        err :: Double
err = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Estimate Double -> Double
forall a. Estimate a -> Double
estimateError Estimate Double
e | (Double
_, Double
_, Estimate Double
e) <- [(Double, Double, Estimate Double)]
panels]
                    _ <- Double -> Either MathError Double
check Double
value
                    if err <= target tol value
                        then Right (Estimate value err count)
                        else
                            if count + 30 > maximumEvaluations tol
                                then Left (NonConvergence "adaptive Gauss-Kronrod" count)
                                else do
                                    let index = (Double, Int) -> Int
forall a b. (a, b) -> b
snd ([(Double, Int)] -> (Double, Int)
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [(Estimate Double -> Double
forall a. Estimate a -> Double
estimateError Estimate Double
e, Int
i) | (Int
i, (Double
_, Double
_, Estimate Double
e)) <- [Int]
-> [(Double, Double, Estimate Double)]
-> [(Int, (Double, Double, Estimate Double))]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [(Double, Double, Estimate Double)]
panels])
                                        (left, right, _) = panels !! index
                                        middle = Double
left Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
right Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
left) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
                                    if middle == left || middle == right
                                        then Left (NonConvergence "adaptive interval resolution" count)
                                        else do
                                            l <- gaussKronrod f left middle
                                            r <- gaussKronrod f middle right
                                            go
                                                (take index panels ++ [(left, middle, l), (middle, right, r)] ++ drop (index + 1) panels)
                                                (count + 30)
            go [(a, b, first)] 15

{- | Clenshaw-Curtis quadrature on n cosine-spaced panels using exact moments
of the cosine interpolant. n must be positive and at most 4096.
-}
clenshawCurtis :: Int -> (Double -> Double) -> Double -> Double -> Either MathError Double
clenshawCurtis :: Int
-> (Double -> Double)
-> Double
-> Double
-> Either MathError Double
clenshawCurtis Int
n Double -> Double
f Double
a Double
b = do
    Bool -> String -> Either MathError ()
require (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
4096 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
a Bool -> Bool -> Bool
&& Double -> Bool
finite Double
b) String
"invalid Clenshaw-Curtis arguments"
    values <-
        (Int -> Either MathError Double)
-> [Int] -> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse
            (\Int
j -> (Double -> Double) -> Double -> Either MathError Double
sample Double -> Double
f ((Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
j Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)))
            [Int
0 .. Int
n]
    let coefficients =
            [ Double
2
                Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
                Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
                    [ (if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 Bool -> Bool -> Bool
|| Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n then Double
0.5 else Double
1)
                        Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
values [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j
                        Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
j) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)
                    | Int
j <- [Int
0 .. Int
n]
                    ]
            | Int
k <- [Int
0 .. Int
n]
            ]
        moment a
k = if a -> Bool
forall a. Integral a => a -> Bool
odd a
k then a
0 else a
2 a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
1 a -> a -> a
forall a. Num a => a -> a -> a
- a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
k a -> a -> a
forall a. Num a => a -> a -> a
* a
k))
        integral = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [(if Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 Bool -> Bool -> Bool
|| Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n then Double
0.5 else Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
coefficients [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall {a} {a}. (Integral a, Fractional a) => a -> a
moment Int
k | Int
k <- [Int
0 .. Int
n]]
    check ((b - a) / 2 * integral)

{- | Tanh-sinh quadrature with geometric step refinement on a finite interval.
Rounded-to-endpoint nodes are omitted, supporting integrable endpoint poles.
-}
tanhSinh ::
    Tolerance -> (Double -> Double) -> Double -> Double -> Either MathError (Estimate Double)
tanhSinh :: Tolerance
-> (Double -> Double)
-> Double
-> Double
-> Either MathError (Estimate Double)
tanhSinh Tolerance
tol Double -> Double
f Double
a Double
b = do
    Tolerance -> Either MathError ()
validateTolerance Tolerance
tol
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
a Bool -> Bool -> Bool
&& Double -> Bool
finite Double
b) String
"tanh-sinh bounds must be finite"
    if Double
a Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
b
        then Estimate Double -> Either MathError (Estimate Double)
forall a b. b -> Either a b
Right (Double -> Double -> Int -> Estimate Double
forall a. a -> Double -> Int -> Estimate a
Estimate Double
0 Double
0 Int
0)
        else do
            let rule :: Double -> Either MathError (Double, Int)
rule Double
h = do
                    let n :: Int
n = Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
ceiling (Double
3.5 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
h) :: Int
                    Bool -> String -> Either MathError ()
require (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Tolerance -> Int
maximumEvaluations Tolerance
tol) String
"tanh-sinh initial budget too small"
                    terms <- (Int -> Either MathError Double)
-> [Int] -> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse (Double -> Int -> Either MathError Double
forall {p}. Integral p => Double -> p -> Either MathError Double
node Double
h) [-Int
n .. Int
n]
                    value <- check (h * sum terms)
                    Right (value, 2 * n + 1)
                node :: Double -> p -> Either MathError Double
node Double
h p
k =
                    let t :: Double
t = p -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral p
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
h
                        u :: Double
u = Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sinh Double
t
                        z :: Double
z = Double -> Double
forall a. Floating a => a -> a
tanh Double
u
                        x :: Double
x = (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z
                        weight :: Double
weight = (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cosh Double
t Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double -> Double
forall a. Floating a => a -> a
cosh Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cosh Double
u)
                     in if Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
a Bool -> Bool -> Bool
|| Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
b then Double -> Either MathError Double
forall a b. b -> Either a b
Right Double
0 else do y <- (Double -> Double) -> Double -> Either MathError Double
sample Double -> Double
f Double
x; check (weight * y)
                go :: Double -> Maybe Double -> Int -> Either MathError (Estimate Double)
go Double
h Maybe Double
old Int
count = do
                    let n :: Int
n = Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
ceiling (Double
3.5 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
h) :: Int) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
                    if Int
count Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Tolerance -> Int
maximumEvaluations Tolerance
tol
                        then MathError -> Either MathError (Estimate Double)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"tanh-sinh" Int
count)
                        else do
                            (value, cost) <- Double -> Either MathError (Double, Int)
rule Double
h
                            let err = Double -> (Double -> Double) -> Maybe Double -> Double
forall b a. b -> (a -> b) -> Maybe a -> b
maybe (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0) (Double -> Double
forall a. Num a => a -> a
abs (Double -> Double) -> (Double -> Double) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double
value Double -> Double -> Double
forall a. Num a => a -> a -> a
-)) Maybe Double
old
                            if err <= target tol value
                                then Right (Estimate value err (count + cost))
                                else go (h / 2) (Just value) (count + cost)
            Double -> Maybe Double -> Int -> Either MathError (Estimate Double)
go Double
1 Maybe Double
forall a. Maybe a
Nothing Int
0

{- | Integrate on finite, half-infinite or doubly-infinite intervals by an open
coordinate transform followed by adaptive Gauss-Kronrod. Reversed bounds
change the sign. Divergent integrals have no guaranteed detection criterion.
-}
integrateInfinite ::
    Tolerance -> (Double -> Double) -> Double -> Double -> Either MathError (Estimate Double)
integrateInfinite :: Tolerance
-> (Double -> Double)
-> Double
-> Double
-> Either MathError (Estimate Double)
integrateInfinite Tolerance
tol Double -> Double
f Double
a Double
b
    | Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
a Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
b = MathError -> Either MathError (Estimate Double)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"NaN integration bound")
    | Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
b = (Estimate Double -> Estimate Double)
-> Either MathError (Estimate Double)
-> Either MathError (Estimate Double)
forall a b. (a -> b) -> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\Estimate Double
e -> Estimate Double
e{estimateValue = negate (estimateValue e)}) (Tolerance
-> (Double -> Double)
-> Double
-> Double
-> Either MathError (Estimate Double)
integrateInfinite Tolerance
tol Double -> Double
f Double
b Double
a)
    | Double -> Bool
finite Double
a Bool -> Bool -> Bool
&& Double -> Bool
finite Double
b = Tolerance
-> (Double -> Double)
-> Double
-> Double
-> Either MathError (Estimate Double)
adaptiveIntegrate Tolerance
tol Double -> Double
f Double
a Double
b
    | Double
a Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== (-Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0) Bool -> Bool -> Bool
&& Double
b Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0 =
        Tolerance
-> (Double -> Double)
-> Double
-> Double
-> Either MathError (Estimate Double)
adaptiveIntegrate Tolerance
tol (\Double
t -> let p :: Double
p = Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.5); c :: Double
c = Double -> Double
forall a. Floating a => a -> a
cos Double
p in Double -> Double
f (Double -> Double
forall a. Floating a => a -> a
tan Double
p) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c)) Double
0 Double
1
    | Double -> Bool
finite Double
a Bool -> Bool -> Bool
&& Double
b Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0 =
        Tolerance
-> (Double -> Double)
-> Double
-> Double
-> Either MathError (Estimate Double)
adaptiveIntegrate Tolerance
tol (\Double
t -> Double -> Double
f (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
t Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) Double
0 Double
1
    | Double
a Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== (-Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0) Bool -> Bool -> Bool
&& Double -> Bool
finite Double
b =
        Tolerance
-> (Double -> Double)
-> Double
-> Double
-> Either MathError (Estimate Double)
adaptiveIntegrate Tolerance
tol (\Double
t -> Double -> Double
f (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) Double
0 Double
1
    | Bool
otherwise = MathError -> Either MathError (Estimate Double)
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"invalid infinite integration bounds")

{- | Integrate f(x)*cos(omega*x) or f(x)*sin(omega*x), splitting at least once
per half-period to avoid accepting an aliased embedded quadrature estimate.
-}
integrateOscillatory ::
    Bool ->
    Double ->
    Tolerance ->
    (Double -> Double) ->
    Double ->
    Double ->
    Either MathError (Estimate Double)
integrateOscillatory :: Bool
-> Double
-> Tolerance
-> (Double -> Double)
-> Double
-> Double
-> Either MathError (Estimate Double)
integrateOscillatory Bool
sine Double
omega Tolerance
tol Double -> Double
f Double
a Double
b = do
    Tolerance -> Either MathError ()
validateTolerance Tolerance
tol
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
omega Bool -> Bool -> Bool
&& Double -> Bool
finite Double
a Bool -> Bool -> Bool
&& Double -> Bool
finite Double
b) String
"invalid oscillatory integration arguments"
    let panelsDouble :: Integer
panelsDouble = Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
max Integer
1 (Double -> Integer
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
ceiling (Double -> Double
forall a. Num a => a -> a
abs (Double
omega Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
forall a. Floating a => a
pi) :: Integer)
    Bool -> String -> Either MathError ()
require
        (Integer
panelsDouble Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
15 Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Tolerance -> Int
maximumEvaluations Tolerance
tol))
        String
"oscillatory budget is too small for the frequency"
    let n :: Int
n = Integer -> Int
forall a. Num a => Integer -> a
fromInteger Integer
panelsDouble :: Int
        trig :: Double -> Double
trig = if Bool
sine then Double -> Double
forall a. Floating a => a -> a
sin else Double -> Double
forall a. Floating a => a -> a
cos
        local :: Tolerance
local =
            Tolerance
tol
                { absoluteTolerance = absoluteTolerance tol / fromIntegral n
                , maximumEvaluations = maximumEvaluations tol `div` n
                }
    estimates <-
        (Int -> Either MathError (Estimate Double))
-> [Int] -> Either MathError [Estimate Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse
            ( \Int
i ->
                Tolerance
-> (Double -> Double)
-> Double
-> Double
-> Either MathError (Estimate Double)
adaptiveIntegrate
                    Tolerance
local
                    (\Double
x -> Double -> Double
f Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
trig (Double
omega Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x))
                    (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)
                    (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)
            )
            [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
    value <- check (sum (map estimateValue estimates))
    Right (Estimate value (sum (map estimateError estimates)) (sum (map evaluations estimates)))

{- | Tensor-product midpoint integration over a box. The product of panel
counts is checked against the explicit budget before callback evaluation.
-}
integrateND :: Int -> Int -> ([Double] -> Double) -> [(Double, Double)] -> Either MathError Double
integrateND :: Int
-> Int
-> ([Double] -> Double)
-> [(Double, Double)]
-> Either MathError Double
integrateND Int
panels Int
budget [Double] -> Double
f [(Double, Double)]
bounds = do
    Bool -> String -> Either MathError ()
require
        (Int
panels Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
budget Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& ((Double, Double) -> Bool) -> [(Double, Double)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\(Double
a, Double
b) -> Double -> Bool
finite Double
a Bool -> Bool -> Bool
&& Double -> Bool
finite Double
b) [(Double, Double)]
bounds)
        String
"invalid multidimensional integration arguments"
    let count :: Integer
count = Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
panels Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ [(Double, Double)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(Double, Double)]
bounds
    Bool -> String -> Either MathError ()
require (Integer
count Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
budget) String
"multidimensional quadrature exceeds evaluation budget"
    let points :: [[Double]]
points =
            [[Double]] -> [[Double]]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
forall (m :: * -> *) a. Monad m => [m a] -> m [a]
sequence
                [ [Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.5) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
panels | Int
i <- [Int
0 .. Int
panels Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
                | (Double
a, Double
b) <- [(Double, Double)]
bounds
                ]
        volume :: Double
volume = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [(Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
panels | (Double
a, Double
b) <- [(Double, Double)]
bounds]
    values <- ([Double] -> Either MathError Double)
-> [[Double]] -> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse (\[Double]
p -> (Double -> Either MathError Double)
-> [Double] -> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse Double -> Either MathError Double
check [Double]
p Either MathError [Double]
-> Either MathError Double -> Either MathError Double
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Double -> Either MathError Double
check ([Double] -> Double
f [Double]
p)) [[Double]]
points
    check (volume * sum values)

{- | Reproducible Monte Carlo box integral with an explicit 64-bit seed.
Error is one estimated standard error from Welford's sample variance.
-}
monteCarlo ::
    Integer -> Int -> ([Double] -> Double) -> [(Double, Double)] -> Either MathError (Estimate Double)
monteCarlo :: Integer
-> Int
-> ([Double] -> Double)
-> [(Double, Double)]
-> Either MathError (Estimate Double)
monteCarlo Integer
seed Int
n [Double] -> Double
f [(Double, Double)]
bounds = do
    Bool -> String -> Either MathError ()
require
        (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
2 Bool -> Bool -> Bool
&& ((Double, Double) -> Bool) -> [(Double, Double)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\(Double
a, Double
b) -> Double -> Bool
finite Double
a Bool -> Bool -> Bool
&& Double -> Bool
finite Double
b) [(Double, Double)]
bounds)
        String
"Monte Carlo needs two samples and finite bounds"
    let next :: a -> a
next a
s = (a
6364136223846793005 a -> a -> a
forall a. Num a => a -> a -> a
* a
s a -> a -> a
forall a. Num a => a -> a -> a
+ a
1442695040888963407) a -> a -> a
forall a. Integral a => a -> a -> a
`mod` (a
2 a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
64 :: Int))
        point :: Integer -> (Integer, [Double])
point Integer
s =
            ((Integer, [Double]) -> (Double, Double) -> (Integer, [Double]))
-> (Integer, [Double]) -> [(Double, Double)] -> (Integer, [Double])
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl
                ( \(Integer
state, [Double]
acc) (Double
a, Double
b) ->
                    let t :: Integer
t = Integer -> Integer
forall {a}. Integral a => a -> a
next Integer
state; u :: Double
u = Integer -> Double
forall a. Num a => Integer -> a
fromInteger (Integer
t Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` Integer
2048) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
9007199254740992
                     in (Integer
t, [Double]
acc [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u])
                )
                (Integer
s, [])
                [(Double, Double)]
bounds
        go :: t
-> Integer
-> t
-> Double
-> Double
-> Either MathError (Double, Double)
go t
0 Integer
_ t
_ Double
mean Double
m2 = (Double, Double) -> Either MathError (Double, Double)
forall a b. b -> Either a b
Right (Double
mean, Double
m2)
        go t
remaining Integer
state t
count Double
mean Double
m2 = do
            let (Integer
state', [Double]
xs) = Integer -> (Integer, [Double])
point Integer
state
            _ <- (Double -> Either MathError Double)
-> [Double] -> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse Double -> Either MathError Double
check [Double]
xs
            y <- check (f xs)
            let count' = t
count t -> t -> t
forall a. Num a => a -> a -> a
+ t
1; delta = Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
mean; mean' = Double
mean Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
delta Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ t -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral t
count'
            go (remaining - 1) state' count' mean' (m2 + delta * (y - mean'))
    (mean, m2) <- Int
-> Integer
-> Int
-> Double
-> Double
-> Either MathError (Double, Double)
forall {t} {t}.
(Integral t, Num t, Eq t) =>
t
-> Integer
-> t
-> Double
-> Double
-> Either MathError (Double, Double)
go Int
n Integer
seed (Int
0 :: Int) Double
0 Double
0
    let volume = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a | (Double
a, Double
b) <- [(Double, Double)]
bounds]
    value <- check (volume * mean)
    err <- check (abs volume * sqrt (max 0 m2 / fromIntegral (n - 1) / fromIntegral n))
    Right (Estimate value err n)

{- | Deterministic Halton-box quadrature, up to 16 dimensions. No probabilistic
error bound is claimed; estimateError is the difference of the two halves.
-}
quasiMonteCarlo ::
    Int -> ([Double] -> Double) -> [(Double, Double)] -> Either MathError (Estimate Double)
quasiMonteCarlo :: Int
-> ([Double] -> Double)
-> [(Double, Double)]
-> Either MathError (Estimate Double)
quasiMonteCarlo Int
n [Double] -> Double
f [(Double, Double)]
bounds = do
    Bool -> String -> Either MathError ()
require
        (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
2 Bool -> Bool -> Bool
&& [(Double, Double)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(Double, Double)]
bounds Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
16 Bool -> Bool -> Bool
&& ((Double, Double) -> Bool) -> [(Double, Double)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\(Double
a, Double
b) -> Double -> Bool
finite Double
a Bool -> Bool -> Bool
&& Double -> Bool
finite Double
b) [(Double, Double)]
bounds)
        String
"invalid Halton quadrature arguments"
    let bases :: [Int]
bases = [Int
2, Int
3, Int
5, Int
7, Int
11, Int
13, Int
17, Int
19, Int
23, Int
29, Int
31, Int
37, Int
41, Int
43, Int
47, Int
53]
        radical :: a -> a -> t
radical a
base a
index = a -> t -> t -> t
forall {t}. Fractional t => a -> t -> t -> t
go a
index (t
1 t -> t -> t
forall a. Fractional a => a -> a -> a
/ a -> t
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
base) t
0
          where
            go :: a -> t -> t -> t
go a
0 t
_ t
acc = t
acc
            go a
k t
factor t
acc = a -> t -> t -> t
go (a
k a -> a -> a
forall a. Integral a => a -> a -> a
`div` a
base) (t
factor t -> t -> t
forall a. Fractional a => a -> a -> a
/ a -> t
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
base) (t
acc t -> t -> t
forall a. Num a => a -> a -> a
+ a -> t
forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
k a -> a -> a
forall a. Integral a => a -> a -> a
`mod` a
base) t -> t -> t
forall a. Num a => a -> a -> a
* t
factor)
        point :: Int -> [Double]
point Int
i = ((Double, Double) -> Int -> Double)
-> [(Double, Double)] -> [Int] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\(Double
a, Double
b) Int
base -> Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Int -> Double
forall {t} {a}. (Fractional t, Integral a) => a -> a -> t
radical Int
base Int
i) [(Double, Double)]
bounds [Int]
bases
        volume :: Double
volume = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a | (Double
a, Double
b) <- [(Double, Double)]
bounds]
    values <- (Int -> Either MathError Double)
-> [Int] -> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse (\Int
i -> let p :: [Double]
p = Int -> [Double]
point Int
i in (Double -> Either MathError Double)
-> [Double] -> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse Double -> Either MathError Double
check [Double]
p Either MathError [Double]
-> Either MathError Double -> Either MathError Double
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Double -> Either MathError Double
check ([Double] -> Double
f [Double]
p)) [Int
1 .. Int
n]
    let half = Int
n Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
        left = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take Int
half [Double]
values) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
half
        right = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
half [Double]
values) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
half)
    value <- check (volume * sum values / fromIntegral n)
    err <- check (abs (volume * (left - right)))
    Right (Estimate value err n)

-- | Cumulative trapezoidal integrals at strictly increasing sample coordinates.
cumulativeIntegral :: [(Double, Double)] -> Either MathError [Double]
cumulativeIntegral :: [(Double, Double)] -> Either MathError [Double]
cumulativeIntegral [] = [Double] -> Either MathError [Double]
forall a b. b -> Either a b
Right []
cumulativeIntegral [(Double, Double)]
values = do
    Bool -> String -> Either MathError ()
require
        ( ((Double, Double) -> Bool) -> [(Double, Double)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\(Double
x, Double
y) -> Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double -> Bool
finite Double
y) [(Double, Double)]
values
            Bool -> Bool -> Bool
&& [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and ((Double -> Double -> Bool) -> [Double] -> [Double] -> [Bool]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
(<) (((Double, Double) -> Double) -> [(Double, Double)] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double, Double) -> Double
forall a b. (a, b) -> a
fst [(Double, Double)]
values) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 (((Double, Double) -> Double) -> [(Double, Double)] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double, Double) -> Double
forall a b. (a, b) -> a
fst [(Double, Double)]
values)))
        )
        String
"samples must be finite with strictly increasing coordinates"
    increments <-
        (((Double, Double), (Double, Double)) -> Either MathError Double)
-> [((Double, Double), (Double, Double))]
-> Either MathError [Double]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse (\((Double
x, Double
y), (Double
u, Double
v)) -> Double -> Either MathError Double
check ((Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
v) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)) ([(Double, Double)]
-> [(Double, Double)] -> [((Double, Double), (Double, Double))]
forall a b. [a] -> [b] -> [(a, b)]
zip [(Double, Double)]
values (Int -> [(Double, Double)] -> [(Double, Double)]
forall a. Int -> [a] -> [a]
drop Int
1 [(Double, Double)]
values))
    traverse check (scanl (+) 0 increments)

-- | Scalar line integral f(r(t))*|r'(t)| along an explicitly parameterized path.
pathIntegral ::
    Tolerance ->
    (Double -> [Double]) ->
    (Double -> [Double]) ->
    ([Double] -> Double) ->
    Double ->
    Double ->
    Either MathError (Estimate Double)
pathIntegral :: Tolerance
-> (Double -> [Double])
-> (Double -> [Double])
-> ([Double] -> Double)
-> Double
-> Double
-> Either MathError (Estimate Double)
pathIntegral Tolerance
tol Double -> [Double]
curve Double -> [Double]
tangent [Double] -> Double
f Double
a Double
b = Tolerance
-> (Double -> Double)
-> Double
-> Double
-> Either MathError (Estimate Double)
adaptiveIntegrate Tolerance
tol Double -> Double
integrand Double
a Double
b
  where
    integrand :: Double -> Double
integrand Double
t =
        let p :: [Double]
p = Double -> [Double]
curve Double
t; v :: [Double]
v = Double -> [Double]
tangent Double
t
         in if [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
p Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
v Bool -> Bool -> Bool
|| Bool -> Bool
not ((Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite ([Double]
p [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
v))
                then Double
0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0
                else [Double] -> Double
f [Double]
p Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) [Double]
v))

-- | Oriented vector-field line integral along a parameterized path.
lineIntegral ::
    Tolerance ->
    (Double -> [Double]) ->
    (Double -> [Double]) ->
    ([Double] -> [Double]) ->
    Double ->
    Double ->
    Either MathError (Estimate Double)
lineIntegral :: Tolerance
-> (Double -> [Double])
-> (Double -> [Double])
-> ([Double] -> [Double])
-> Double
-> Double
-> Either MathError (Estimate Double)
lineIntegral Tolerance
tol Double -> [Double]
curve Double -> [Double]
tangent [Double] -> [Double]
field Double
a Double
b = Tolerance
-> (Double -> Double)
-> Double
-> Double
-> Either MathError (Estimate Double)
adaptiveIntegrate Tolerance
tol Double -> Double
integrand Double
a Double
b
  where
    integrand :: Double -> Double
integrand Double
t =
        let p :: [Double]
p = Double -> [Double]
curve Double
t; v :: [Double]
v = Double -> [Double]
tangent Double
t
         in if [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
p Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
v Bool -> Bool -> Bool
|| Bool -> Bool
not ((Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite ([Double]
p [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
v))
                then Double
0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0
                else
                    let w :: [Double]
w = [Double] -> [Double]
field [Double]
p
                     in if [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
v Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
w Bool -> Bool -> Bool
|| Bool -> Bool
not ((Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
w) then Double
0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0 else [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) [Double]
v [Double]
w)

-- | Scalar surface integral over a rectangular parameter domain in R3.
surfaceIntegral ::
    Int ->
    Int ->
    ([Double] -> [Double]) ->
    ([Double] -> [Double]) ->
    ([Double] -> [Double]) ->
    ([Double] -> Double) ->
    [(Double, Double)] ->
    Either MathError Double
surfaceIntegral :: Int
-> Int
-> ([Double] -> [Double])
-> ([Double] -> [Double])
-> ([Double] -> [Double])
-> ([Double] -> Double)
-> [(Double, Double)]
-> Either MathError Double
surfaceIntegral Int
panels Int
budget [Double] -> [Double]
surface [Double] -> [Double]
du [Double] -> [Double]
dv [Double] -> Double
f [(Double, Double)]
bounds = do
    Bool -> String -> Either MathError ()
require ([(Double, Double)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(Double, Double)]
bounds Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
2) String
"surface parameter domain must have two dimensions"
    let integrand :: [Double] -> Double
integrand [Double]
uv = case ([Double] -> [Double]
surface [Double]
uv, [Double] -> [Double]
du [Double]
uv, [Double] -> [Double]
dv [Double]
uv) of
            ([Double]
p, [Double
a, Double
b, Double
c], [Double
x, Double
y, Double
z])
                | [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
p Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
3 Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite ([Double]
p [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double
a, Double
b, Double
c, Double
x, Double
y, Double
z]) ->
                    [Double] -> Double
f [Double]
p
                        Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt ((Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int))
            ([Double], [Double], [Double])
_ -> Double
0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0
    Int
-> Int
-> ([Double] -> Double)
-> [(Double, Double)]
-> Either MathError Double
integrateND Int
panels Int
budget [Double] -> Double
integrand [(Double, Double)]
bounds

-- | Volume integral over a three-dimensional rectangular box.
volumeIntegral ::
    Int -> Int -> ([Double] -> Double) -> [(Double, Double)] -> Either MathError Double
volumeIntegral :: Int
-> Int
-> ([Double] -> Double)
-> [(Double, Double)]
-> Either MathError Double
volumeIntegral Int
panels Int
budget [Double] -> Double
f [(Double, Double)]
bounds =
    Bool -> String -> Either MathError ()
require ([(Double, Double)] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(Double, Double)]
bounds Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
3) String
"volume domain must have three dimensions"
        Either MathError ()
-> Either MathError Double -> Either MathError Double
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Int
-> Int
-> ([Double] -> Double)
-> [(Double, Double)]
-> Either MathError Double
integrateND Int
panels Int
budget [Double] -> Double
f [(Double, Double)]
bounds