-- SPDX-License-Identifier: MPL-2.0

{- | Checked reference interpolation and approximation. Samples are strictly
increasing, finite, and unique; no silent sorting or duplicate averaging.
Piecewise coefficients are ascending powers of (x-left knot), intervals are
left closed/right open except the final endpoint. Error results are sampled
estimates, not mathematical enclosures. Dense fitting is bounded to 1024
samples, degree 64, and 16 scattered dimensions.
-}
module HMath.Interpolation (
    Extrapolation (..),
    InterpolationMethod (..),
    interpolate,
    hermite,
    cubicHermite,
    SplineBoundary (..),
    PiecewisePolynomial,
    piecewisePolynomial,
    piecewiseKnots,
    evaluatePiecewise,
    differentiatePiecewise,
    integratePiecewise,
    cubicSpline,
    pchip,
    akima,
    smoothingSpline,
    BSpline,
    bSpline,
    evaluateBSpline,
    derivativeBSpline,
    insertKnot,
    bezier,
    rationalBezier,
    cardinalSpline,
    catmullRom,
    kochanekBartels,
    curveIntegral,
    arcLength,
    gridInterpolate,
    RadialKernel (..),
    RBF,
    fitRBF,
    evaluateRBF,
    inverseDistance,
    shepard,
    kriging,
    Approximation (..),
    polynomialFit,
    chebyshevApproximate,
    evaluateChebyshev,
    pade,
    rationalFit,
    remez,
    nonlinearFit,
    interpolationError,
) where

import qualified HMath.Calculus as C
import HMath.Error (MathError (..))
import qualified HMath.Matrix as M
import qualified HMath.Matrix.QR as QR
import HMath.Signal.Internal

-- | Reject out-of-domain queries, clamp to an endpoint, or evaluate the endpoint polynomial beyond its interval.
data Extrapolation = Reject | Clamp | Extrapolate deriving (Extrapolation -> Extrapolation -> Bool
(Extrapolation -> Extrapolation -> Bool)
-> (Extrapolation -> Extrapolation -> Bool) -> Eq Extrapolation
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Extrapolation -> Extrapolation -> Bool
== :: Extrapolation -> Extrapolation -> Bool
$c/= :: Extrapolation -> Extrapolation -> Bool
/= :: Extrapolation -> Extrapolation -> Bool
Eq, Int -> Extrapolation -> ShowS
[Extrapolation] -> ShowS
Extrapolation -> String
(Int -> Extrapolation -> ShowS)
-> (Extrapolation -> String)
-> ([Extrapolation] -> ShowS)
-> Show Extrapolation
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Extrapolation -> ShowS
showsPrec :: Int -> Extrapolation -> ShowS
$cshow :: Extrapolation -> String
show :: Extrapolation -> String
$cshowList :: [Extrapolation] -> ShowS
showList :: [Extrapolation] -> ShowS
Show)

-- | Independent scalar nearest, linear, Lagrange, Newton divided-difference or barycentric formula.
data InterpolationMethod = Nearest | Linear | Lagrange | Newton | Barycentric deriving (InterpolationMethod -> InterpolationMethod -> Bool
(InterpolationMethod -> InterpolationMethod -> Bool)
-> (InterpolationMethod -> InterpolationMethod -> Bool)
-> Eq InterpolationMethod
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: InterpolationMethod -> InterpolationMethod -> Bool
== :: InterpolationMethod -> InterpolationMethod -> Bool
$c/= :: InterpolationMethod -> InterpolationMethod -> Bool
/= :: InterpolationMethod -> InterpolationMethod -> Bool
Eq, Int -> InterpolationMethod -> ShowS
[InterpolationMethod] -> ShowS
InterpolationMethod -> String
(Int -> InterpolationMethod -> ShowS)
-> (InterpolationMethod -> String)
-> ([InterpolationMethod] -> ShowS)
-> Show InterpolationMethod
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> InterpolationMethod -> ShowS
showsPrec :: Int -> InterpolationMethod -> ShowS
$cshow :: InterpolationMethod -> String
show :: InterpolationMethod -> String
$cshowList :: [InterpolationMethod] -> ShowS
showList :: [InterpolationMethod] -> ShowS
Show)

-- | Natural zero curvature, supplied clamped endpoint slopes, or periodic value/derivative matching.
data SplineBoundary = Natural | Clamped Double Double | Periodic deriving (SplineBoundary -> SplineBoundary -> Bool
(SplineBoundary -> SplineBoundary -> Bool)
-> (SplineBoundary -> SplineBoundary -> Bool) -> Eq SplineBoundary
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: SplineBoundary -> SplineBoundary -> Bool
== :: SplineBoundary -> SplineBoundary -> Bool
$c/= :: SplineBoundary -> SplineBoundary -> Bool
/= :: SplineBoundary -> SplineBoundary -> Bool
Eq, Int -> SplineBoundary -> ShowS
[SplineBoundary] -> ShowS
SplineBoundary -> String
(Int -> SplineBoundary -> ShowS)
-> (SplineBoundary -> String)
-> ([SplineBoundary] -> ShowS)
-> Show SplineBoundary
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> SplineBoundary -> ShowS
showsPrec :: Int -> SplineBoundary -> ShowS
$cshow :: SplineBoundary -> String
show :: SplineBoundary -> String
$cshowList :: [SplineBoundary] -> ShowS
showList :: [SplineBoundary] -> ShowS
Show)

-- | Opaque ascending local polynomial coefficients and strictly increasing breakpoints.
data PiecewisePolynomial = PiecewisePolynomial [Double] [[Double]] deriving (PiecewisePolynomial -> PiecewisePolynomial -> Bool
(PiecewisePolynomial -> PiecewisePolynomial -> Bool)
-> (PiecewisePolynomial -> PiecewisePolynomial -> Bool)
-> Eq PiecewisePolynomial
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: PiecewisePolynomial -> PiecewisePolynomial -> Bool
== :: PiecewisePolynomial -> PiecewisePolynomial -> Bool
$c/= :: PiecewisePolynomial -> PiecewisePolynomial -> Bool
/= :: PiecewisePolynomial -> PiecewisePolynomial -> Bool
Eq, Int -> PiecewisePolynomial -> ShowS
[PiecewisePolynomial] -> ShowS
PiecewisePolynomial -> String
(Int -> PiecewisePolynomial -> ShowS)
-> (PiecewisePolynomial -> String)
-> ([PiecewisePolynomial] -> ShowS)
-> Show PiecewisePolynomial
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> PiecewisePolynomial -> ShowS
showsPrec :: Int -> PiecewisePolynomial -> ShowS
$cshow :: PiecewisePolynomial -> String
show :: PiecewisePolynomial -> String
$cshowList :: [PiecewisePolynomial] -> ShowS
showList :: [PiecewisePolynomial] -> ShowS
Show)

-- | Opaque rational B-spline with positive weights and a validated full knot vector.
data BSpline = BSpline Int [Double] [[Double]] [Double] deriving (BSpline -> BSpline -> Bool
(BSpline -> BSpline -> Bool)
-> (BSpline -> BSpline -> Bool) -> Eq BSpline
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: BSpline -> BSpline -> Bool
== :: BSpline -> BSpline -> Bool
$c/= :: BSpline -> BSpline -> Bool
/= :: BSpline -> BSpline -> Bool
Eq, Int -> BSpline -> ShowS
[BSpline] -> ShowS
BSpline -> String
(Int -> BSpline -> ShowS)
-> (BSpline -> String) -> ([BSpline] -> ShowS) -> Show BSpline
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> BSpline -> ShowS
showsPrec :: Int -> BSpline -> ShowS
$cshow :: BSpline -> String
show :: BSpline -> String
$cshowList :: [BSpline] -> ShowS
showList :: [BSpline] -> ShowS
Show)

-- | Radial interpolation kernels; the first three accept a positive inverse-length shape parameter.
data RadialKernel
    = Gaussian Double
    | Multiquadric Double
    | InverseMultiquadric Double
    | ThinPlate
    | Cubic
    deriving (RadialKernel -> RadialKernel -> Bool
(RadialKernel -> RadialKernel -> Bool)
-> (RadialKernel -> RadialKernel -> Bool) -> Eq RadialKernel
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: RadialKernel -> RadialKernel -> Bool
== :: RadialKernel -> RadialKernel -> Bool
$c/= :: RadialKernel -> RadialKernel -> Bool
/= :: RadialKernel -> RadialKernel -> Bool
Eq, Int -> RadialKernel -> ShowS
[RadialKernel] -> ShowS
RadialKernel -> String
(Int -> RadialKernel -> ShowS)
-> (RadialKernel -> String)
-> ([RadialKernel] -> ShowS)
-> Show RadialKernel
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> RadialKernel -> ShowS
showsPrec :: Int -> RadialKernel -> ShowS
$cshow :: RadialKernel -> String
show :: RadialKernel -> String
$cshowList :: [RadialKernel] -> ShowS
showList :: [RadialKernel] -> ShowS
Show)

-- | Opaque fitted radial model with an affine polynomial tail.
data RBF = RBF RadialKernel [[Double]] [Double] [Double] deriving (RBF -> RBF -> Bool
(RBF -> RBF -> Bool) -> (RBF -> RBF -> Bool) -> Eq RBF
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: RBF -> RBF -> Bool
== :: RBF -> RBF -> Bool
$c/= :: RBF -> RBF -> Bool
/= :: RBF -> RBF -> Bool
Eq, Int -> RBF -> ShowS
[RBF] -> ShowS
RBF -> String
(Int -> RBF -> ShowS)
-> (RBF -> String) -> ([RBF] -> ShowS) -> Show RBF
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> RBF -> ShowS
showsPrec :: Int -> RBF -> ShowS
$cshow :: RBF -> String
show :: RBF -> String
$cshowList :: [RBF] -> ShowS
showList :: [RBF] -> ShowS
Show)

-- | Coefficients or parameters, observed absolute residual, and iteration/evaluation count; error is an estimate.
data Approximation = Approximation
    { Approximation -> [Double]
approximationCoefficients :: [Double]
    -- ^ Fitted coefficients in the approximation routine's documented basis/order.
    , Approximation -> Double
approximationError :: Double
    -- ^ Reported approximation residual/error measure, not a general analytic bound.
    , Approximation -> Int
approximationIterations :: Int
    -- ^ Number of fitting/refinement iterations performed.
    }
    deriving (Approximation -> Approximation -> Bool
(Approximation -> Approximation -> Bool)
-> (Approximation -> Approximation -> Bool) -> Eq Approximation
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Approximation -> Approximation -> Bool
== :: Approximation -> Approximation -> Bool
$c/= :: Approximation -> Approximation -> Bool
/= :: Approximation -> Approximation -> Bool
Eq, Int -> Approximation -> ShowS
[Approximation] -> ShowS
Approximation -> String
(Int -> Approximation -> ShowS)
-> (Approximation -> String)
-> ([Approximation] -> ShowS)
-> Show Approximation
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Approximation -> ShowS
showsPrec :: Int -> Approximation -> ShowS
$cshow :: Approximation -> String
show :: Approximation -> String
$cshowList :: [Approximation] -> ShowS
showList :: [Approximation] -> ShowS
Show)

samples :: [(Double, Double)] -> Either MathError ([Double], [Double])
samples :: [(Double, Double)] -> Either MathError ([Double], [Double])
samples [(Double, Double)]
points = do
    let ([Double]
xs, [Double]
ys) = [(Double, Double)] -> ([Double], [Double])
forall a b. [(a, b)] -> ([a], [b])
unzip [(Double, Double)]
points
    Bool -> String -> Either MathError ()
require
        (Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
xs) Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
1024 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]
xs (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
xs)))
        String
"samples must be strictly increasing, unique and nonempty (maximum 1024)"
    _ <- [Double] -> Either MathError [Double]
vector ([Double]
xs [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
ys)
    Right (xs, ys)
pointPolicy :: Extrapolation -> [Double] -> Double -> Either MathError Double
pointPolicy :: Extrapolation -> [Double] -> Double -> Either MathError Double
pointPolicy Extrapolation
policy [Double]
xs Double
x = do
    _ <- Double -> Either MathError Double
checked Double
x
    let lo = [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0; hi = [Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
xs
    require (policy /= Reject || x >= lo && x <= hi) "point outside interpolation domain"
    Right (if policy == Clamp then max lo (min hi x) else x)

-- | Interpolate ordered (x,y) samples using the selected formula and extrapolation policy; nearest ties choose the earlier node.
interpolate ::
    InterpolationMethod -> Extrapolation -> [(Double, Double)] -> Double -> Either MathError Double
interpolate :: InterpolationMethod
-> Extrapolation
-> [(Double, Double)]
-> Double
-> Either MathError Double
interpolate InterpolationMethod
method Extrapolation
policy [(Double, Double)]
points Double
x0 = do
    (xs, ys) <- [(Double, Double)] -> Either MathError ([Double], [Double])
samples [(Double, Double)]
points
    x <- pointPolicy policy xs x0
    let n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs
        interval = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2) (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
0 ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((Double -> Bool) -> [Double] -> [Double]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
x) [Double]
xs) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1))
        polynomialValue = case InterpolationMethod
method of
            InterpolationMethod
Nearest -> (Double, Double) -> Double
forall a b. (a, b) -> b
snd (Double -> [(Double, Double)] -> (Double, Double)
forall {t :: * -> *} {a} {b}.
(Foldable t, Ord a, Num a) =>
a -> t (a, b) -> (a, b)
minimumByDistance Double
x [(Double, Double)]
points)
            InterpolationMethod
Linear ->
                if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1
                    then [Double]
ys [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0
                    else
                        let i :: Int
i = Int
interval; t :: Double
t = (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ([Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)
                         in (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t) Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
ys [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
ys [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
            InterpolationMethod
Lagrange ->
                [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
                    [ [Double]
ys [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i 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
product [(Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ([Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
i]
                    | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                    ]
            InterpolationMethod
Newton ->
                let columns :: [[Double]]
columns =
                        ([Double] -> Int -> [Double]) -> [Double] -> [Int] -> [[Double]]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl
                            (\[Double]
v Int
k -> [([Double]
v [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
v [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ([Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
k) Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) | Int
i <- [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
v Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2]])
                            [Double]
ys
                            [Int
1 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                    -- Explicit divided-difference table, separate from Lagrange.
                    coeffs :: [Double]
coeffs = ([Double] -> Double) -> [[Double]] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) [[Double]]
columns
                 in [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [[Double]
coeffs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
k 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
product (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take Int
k ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
-) [Double]
xs)) | Int
k <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
            InterpolationMethod
Barycentric -> case Double -> [(Double, Double)] -> Maybe Double
forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup Double
x [(Double, Double)]
points of
                Just Double
y -> Double
y
                Maybe Double
Nothing ->
                    let weights :: [Double]
weights =
                            [Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Double
xi Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
xj | (Int
j, Double
xj) <- [Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [Double]
xs, Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
i] | (Int
i, Double
xi) <- [Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [Double]
xs]
                        terms :: [Double]
terms = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
w Double
xi -> Double
w Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
xi)) [Double]
weights [Double]
xs
                     in [Double] -> [Double] -> Double
dot [Double]
terms [Double]
ys Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
terms
    checked polynomialValue
  where
    minimumByDistance :: a -> t (a, b) -> (a, b)
minimumByDistance a
x = ((a, b) -> (a, b) -> (a, b)) -> t (a, b) -> (a, b)
forall a. (a -> a -> a) -> t a -> a
forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 (\(a, b)
a (a, b)
b -> if a -> a
forall a. Num a => a -> a
abs ((a, b) -> a
forall a b. (a, b) -> a
fst (a, b)
a a -> a -> a
forall a. Num a => a -> a -> a
- a
x) a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a -> a
forall a. Num a => a -> a
abs ((a, b) -> a
forall a b. (a, b) -> a
fst (a, b)
b a -> a -> a
forall a. Num a => a -> a -> a
- a
x) then (a, b)
a else (a, b)
b)

{- | Global confluent Hermite polynomial, distinct nodes and supplied first
derivatives. Uses the squared Lagrange cardinal functions.
-}
hermite :: [(Double, Double, Double)] -> Double -> Either MathError Double
hermite :: [(Double, Double, Double)] -> Double -> Either MathError Double
hermite [(Double, Double, Double)]
triples Double
x = do
    (xs, _) <- [(Double, Double)] -> Either MathError ([Double], [Double])
samples [(Double
a, Double
b) | (Double
a, Double
b, Double
_) <- [(Double, Double, Double)]
triples]
    _ <- vector (x : [d | (_, _, d) <- triples])
    let n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs
        term Int
i (a
_, Double
y, Double
d) =
            let xi :: Double
xi = [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i
                l :: Double
l = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [(Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
xi Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
i]
                slope :: Double
slope = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
xi Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
i]
             in ((Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
xi) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
slope) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
xi) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
l Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
l
    checked (sum (zipWith term [0 ..] triples))

-- | Validate breakpoints and nonempty coefficient rows. Each row uses ascending powers of displacement from its left knot.
piecewisePolynomial :: [Double] -> [[Double]] -> Either MathError PiecewisePolynomial
piecewisePolynomial :: [Double] -> [[Double]] -> Either MathError PiecewisePolynomial
piecewisePolynomial [Double]
knots [[Double]]
coefficients = do
    Bool -> String -> Either MathError ()
require
        ( [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
knots Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
2
            Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
knots Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
1024
            Bool -> Bool -> Bool
&& [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
coefficients 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]
knots Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
            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]
knots (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
knots))
            Bool -> Bool -> Bool
&& ([Double] -> Bool) -> [[Double]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\[Double]
r -> Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
r) Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
r Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
65) [[Double]]
coefficients
        )
        String
"invalid piecewise polynomial knots/coefficients"
    _ <- [Double] -> Either MathError [Double]
vector ([Double]
knots [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Double]]
coefficients)
    Right (PiecewisePolynomial knots coefficients)

-- | Return the immutable piecewise breakpoints, including both outer endpoints.
piecewiseKnots :: PiecewisePolynomial -> [Double]
piecewiseKnots :: PiecewisePolynomial -> [Double]
piecewiseKnots (PiecewisePolynomial [Double]
xs [[Double]]
_) = [Double]
xs

-- | Evaluate by Horner rule in the containing interval; use the final interval at the right endpoint.
evaluatePiecewise :: Extrapolation -> PiecewisePolynomial -> Double -> Either MathError Double
evaluatePiecewise :: Extrapolation
-> PiecewisePolynomial -> Double -> Either MathError Double
evaluatePiecewise Extrapolation
policy (PiecewisePolynomial [Double]
xs [[Double]]
coeffs) Double
x0 = do
    x <- Extrapolation -> [Double] -> Double -> Either MathError Double
pointPolicy Extrapolation
policy [Double]
xs Double
x0
    let i = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min ([[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
coeffs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
0 ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((Double -> Bool) -> [Double] -> [Double]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
x) [Double]
xs) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1))
    checked (polyEval (coeffs !! i) (x - xs !! i))

polyEval :: [Double] -> Double -> Double
polyEval :: [Double] -> Double -> Double
polyEval [Double]
cs Double
x = (Double -> Double -> Double) -> Double -> [Double] -> Double
forall a b. (a -> b -> b) -> b -> [a] -> b
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\Double
a Double
b -> Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b) Double
0 [Double]
cs

-- | Differentiate each local coefficient row without changing the breakpoints.
differentiatePiecewise :: PiecewisePolynomial -> Either MathError PiecewisePolynomial
differentiatePiecewise :: PiecewisePolynomial -> Either MathError PiecewisePolynomial
differentiatePiecewise (PiecewisePolynomial [Double]
xs [[Double]]
cs) = [Double] -> [[Double]] -> Either MathError PiecewisePolynomial
piecewisePolynomial [Double]
xs [if [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
r Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 then [Double
0] else (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
1 ..] (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
r) | [Double]
r <- [[Double]]
cs]

-- | Exact polynomial antiderivative across all crossed intervals; both limits must lie inside the domain.
integratePiecewise :: PiecewisePolynomial -> Double -> Double -> Either MathError Double
integratePiecewise :: PiecewisePolynomial -> Double -> Double -> Either MathError Double
integratePiecewise pp :: PiecewisePolynomial
pp@(PiecewisePolynomial [Double]
xs [[Double]]
cs) Double
a Double
b = do
    _ <- Extrapolation -> [Double] -> Double -> Either MathError Double
pointPolicy Extrapolation
Reject [Double]
xs Double
a
    _ <- pointPolicy Reject xs b
    if a > b
        then negate <$> integratePiecewise pp b a
        else
            checked
                ( sum
                    [ primitive i (min b (xs !! (i + 1))) - primitive i (max a (xs !! i))
                    | i <- [0 .. length cs - 1]
                    , max a (xs !! i) < min b (xs !! (i + 1))
                    ]
                )
  where
    primitive :: Int -> Double -> Double
primitive Int
i Double
x = [Double] -> Double -> Double
polyEval (Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Fractional a => a -> a -> a
(/) ([[Double]]
cs [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) [Double
1 ..]) (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)

{- | Piecewise cubic Hermite interpolation with supplied endpoint derivatives.
Triples are (coordinate,value,slope); first derivatives are continuous at
shared knots, while second derivatives need not agree.
-}
cubicHermite :: [(Double, Double, Double)] -> Either MathError PiecewisePolynomial
cubicHermite :: [(Double, Double, Double)] -> Either MathError PiecewisePolynomial
cubicHermite [(Double, Double, Double)]
triples = do
    (xs, ys) <- [(Double, Double)] -> Either MathError ([Double], [Double])
samples [(Double
x, Double
y) | (Double
x, Double
y, Double
_) <- [(Double, Double, Double)]
triples]
    let derivatives = [Double
d | (Double
_, Double
_, Double
d) <- [(Double, Double, Double)]
triples]
    require (length xs >= 2) "cubic Hermite needs two samples"
    _ <- vector derivatives
    hermitePieces xs ys derivatives

hermitePieces :: [Double] -> [Double] -> [Double] -> Either MathError PiecewisePolynomial
hermitePieces :: [Double]
-> [Double] -> [Double] -> Either MathError PiecewisePolynomial
hermitePieces [Double]
xs [Double]
ys [Double]
ds =
    [Double] -> [[Double]] -> Either MathError PiecewisePolynomial
piecewisePolynomial
        [Double]
xs
        [ let h :: Double
h = [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i; s :: Double
s = ([Double]
ys [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
ys [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
h
           in [ [Double]
ys [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i
              , [Double]
ds [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i
              , (Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
ds [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
ds [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
h
              , ([Double]
ds [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double]
ds [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
h)
              ]
        | Int
i <- [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
2]
        ]

{- | Cubic C2 interpolation. Periodic requires equal endpoint values and at least
three intervals; endpoint first and second derivatives agree.
-}
cubicSpline :: SplineBoundary -> [(Double, Double)] -> Either MathError PiecewisePolynomial
cubicSpline :: SplineBoundary
-> [(Double, Double)] -> Either MathError PiecewisePolynomial
cubicSpline SplineBoundary
boundary [(Double, Double)]
points = do
    (xs, ys) <- [(Double, Double)] -> Either MathError ([Double], [Double])
samples [(Double, Double)]
points
    let n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs; h = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
xs) [Double]
xs; slopes = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Fractional a => a -> a -> a
(/) ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
ys) [Double]
ys) [Double]
h
    require (n >= 2) "cubic spline needs two samples"
    second <- case boundary of
        SplineBoundary
Periodic -> do
            Bool -> String -> Either MathError ()
require (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
4 Bool -> Bool -> Bool
&& [Double]
ys [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0 Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
ys) String
"periodic cubic needs >=4 samples and equal endpoint values"
            let m :: Int
m = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
                row :: Int -> [Double]
row Int
i =
                    [ (if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i then Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! ((Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
m) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) else Double
0)
                        Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
m then [Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! ((Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
m) else Double
0)
                        Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
m then [Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i else Double
0)
                    | Int
j <- [Int
0 .. Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                    ]
            v <-
                [[Double]] -> [Double] -> Either MathError [Double]
solve
                    [Int -> [Double]
row Int
i | Int
i <- [Int
0 .. Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
                    [Double
6 Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
slopes [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
slopes [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! ((Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
m)) | Int
i <- [Int
0 .. Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
            Right (v ++ [v !! 0])
        SplineBoundary
_ -> do
            let first :: [Double]
first = case SplineBoundary
boundary of
                    Clamped Double
_ Double
_ -> [Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0, [Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0] [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2) Double
0
                    SplineBoundary
_ -> Double
1 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Double
0
                end :: [Double]
end = case SplineBoundary
boundary of
                    Clamped Double
_ Double
_ -> Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2) Double
0 [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [[Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
h, Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
h]
                    SplineBoundary
_ -> Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Double
0 [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double
1]
                middle :: Int -> [Double]
middle Int
i =
                    [ if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
                        then [Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
                        else if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i then Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) else if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 then [Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i else Double
0
                    | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                    ]
                left :: Double
left = case SplineBoundary
boundary of Clamped Double
d Double
_ -> Double
6 Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
slopes [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
d); SplineBoundary
_ -> Double
0
                right :: Double
right = case SplineBoundary
boundary of Clamped Double
_ Double
d -> Double
6 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
slopes); SplineBoundary
_ -> Double
0
            [[Double]] -> [Double] -> Either MathError [Double]
solve
                ([Double]
first [Double] -> [[Double]] -> [[Double]]
forall a. a -> [a] -> [a]
: [Int -> [Double]
middle Int
i | Int
i <- [Int
1 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2]] [[Double]] -> [[Double]] -> [[Double]]
forall a. [a] -> [a] -> [a]
++ [[Double]
end])
                (Double
left Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [Double
6 Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
slopes [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
slopes [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)) | Int
i <- [Int
1 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2]] [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double
right])
    piecewisePolynomial
        xs
        [ [ ys !! i
          , slopes !! i - h !! i * (2 * second !! i + second !! (i + 1)) / 6
          , second !! i / 2
          , (second !! (i + 1) - second !! i) / (6 * h !! i)
          ]
        | i <- [0 .. n - 2]
        ]

{- | Fritsch--Carlson/PCHIP weighted harmonic interior slopes and limited one-
sided endpoints preserve monotonicity on monotone data.
-}
pchip :: [(Double, Double)] -> Either MathError PiecewisePolynomial
pchip :: [(Double, Double)] -> Either MathError PiecewisePolynomial
pchip [(Double, Double)]
points = do
    (xs, ys) <- [(Double, Double)] -> Either MathError ([Double], [Double])
samples [(Double, Double)]
points
    require (length xs >= 2) "PCHIP needs two samples"
    let n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs
        h = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
xs) [Double]
xs
        s = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Fractional a => a -> a -> a
(/) ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
ys) [Double]
ys) [Double]
h
        endpoint a
h0 a
h1 a
s0 a
s1 =
            let d :: a
d = ((a
2 a -> a -> a
forall a. Num a => a -> a -> a
* a
h0 a -> a -> a
forall a. Num a => a -> a -> a
+ a
h1) a -> a -> a
forall a. Num a => a -> a -> a
* a
s0 a -> a -> a
forall a. Num a => a -> a -> a
- a
h0 a -> a -> a
forall a. Num a => a -> a -> a
* a
s1) a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
h0 a -> a -> a
forall a. Num a => a -> a -> a
+ a
h1)
             in if a
d a -> a -> a
forall a. Num a => a -> a -> a
* a
s0 a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0 then a
0 else if a
s0 a -> a -> a
forall a. Num a => a -> a -> a
* a
s1 a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 Bool -> Bool -> Bool
&& a -> a
forall a. Num a => a -> a
abs a
d a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
3 a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Num a => a -> a
abs a
s0 then a
3 a -> a -> a
forall a. Num a => a -> a -> a
* a
s0 else a
d
        interior Int
i =
            if [Double]
s [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
s [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0
                then Double
0
                else
                    let w1 :: Double
w1 = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1); w2 :: Double
w2 = [Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
                     in (Double
w1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
w2) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
w1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double]
s [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
w2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double]
s [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)
        ds =
            if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
2
                then [[Double]
s [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0, [Double]
s [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0]
                else
                    Double -> Double -> Double -> Double -> Double
forall {a}. (Ord a, Fractional a) => a -> a -> a -> a -> a
endpoint ([Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) ([Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
1) ([Double]
s [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) ([Double]
s [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
1)
                        Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [Int -> Double
interior Int
i | Int
i <- [Int
1 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2]]
                        [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double -> Double -> Double -> Double -> Double
forall {a}. (Ord a, Fractional a) => a -> a -> a -> a -> a
endpoint ([Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
h) ([Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
3)) ([Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
s) ([Double]
s [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
3))]
    hermitePieces xs ys ds

-- | Akima cubic interpolation using slope-difference weights and linearly extended endpoint secants; at least three samples.
akima :: [(Double, Double)] -> Either MathError PiecewisePolynomial
akima :: [(Double, Double)] -> Either MathError PiecewisePolynomial
akima [(Double, Double)]
points = do
    (xs, ys) <- [(Double, Double)] -> Either MathError ([Double], [Double])
samples [(Double, Double)]
points
    require (length xs >= 3) "Akima needs three samples"
    let s = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Fractional a => a -> a -> a
(/) ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
ys) [Double]
ys) ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
xs) [Double]
xs)
        ext =
            [Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
s [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
s [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
1, Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
s [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
s [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
1]
                [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
s
                [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
s Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
s [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
s Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2), Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
s Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
s [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
s Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2)]
        ds =
            [ let a :: Double
a = Double -> Double
forall a. Num a => a -> a
abs ([Double]
ext [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
3) Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
ext [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2)); b :: Double
b = Double -> Double
forall a. Num a => a -> a
abs ([Double]
ext [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
ext [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)
               in if Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0
                    then ([Double]
ext [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double]
ext [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
                    else (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
ext [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
ext [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b)
            | Int
i <- [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]
            ]
    hermitePieces xs ys ds

{- | Natural cubic smoothing spline minimizing sum (y_i-s(x_i))² plus lambda
times integral s''². Reinsch Q/R formulation, unit observation weights.
-}
smoothingSpline :: Double -> [(Double, Double)] -> Either MathError PiecewisePolynomial
smoothingSpline :: Double
-> [(Double, Double)] -> Either MathError PiecewisePolynomial
smoothingSpline Double
lambda [(Double, Double)]
points = do
    (xs, ys) <- [(Double, Double)] -> Either MathError ([Double], [Double])
samples [(Double, Double)]
points
    require (finite lambda && lambda >= 0 && length xs >= 3) "invalid smoothing spline"
    let n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs
        m = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2
        h = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
xs) [Double]
xs
        q =
            [ [ if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j
                    then Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j
                    else
                        if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 then -Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) else if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2 then Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) else Double
0
              | Int
j <- [Int
0 .. Int
m 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]
            ]
        r =
            [ [ if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then ([Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3 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]
h [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
i Int
j Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
6 else Double
0
              | Int
j <- [Int
0 .. Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
              ]
            | Int
i <- [Int
0 .. Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
            ]
        system = ([Double] -> [Double] -> [Double])
-> [[Double]] -> [[Double]] -> [[Double]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith ((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]]
r (([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double
lambda Double -> Double -> Double
forall a. Num a => a -> a -> a
*)) ([[Double]] -> [[Double]] -> [[Double]]
mm ([[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
tr [[Double]]
q) [[Double]]
q))
    gamma <- solve system (mv (tr q) ys)
    fitted <- vector (zipWith (-) ys (map (lambda *) (mv q gamma)))
    cubicSpline Natural (zip xs fitted)

{- | Nondecreasing full knot vector, degree 0..16, positive rational weights,
control points of equal dimension. Evaluation domain [k_p,k_n].
-}
bSpline :: Int -> [Double] -> [[Double]] -> [Double] -> Either MathError BSpline
bSpline :: Int
-> [Double] -> [[Double]] -> [Double] -> Either MathError BSpline
bSpline Int
degree [Double]
knots [[Double]]
controls [Double]
weights = do
    let n :: Int
n = [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
controls; d :: Int
d = if [[Double]] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [[Double]]
controls then Int
0 else [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([[Double]]
controls [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
0)
    Bool -> String -> Either MathError ()
require
        ( Int
degree Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0
            Bool -> Bool -> Bool
&& Int
degree Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
16
            Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
degree
            Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
1024
            Bool -> Bool -> Bool
&& Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0
            Bool -> Bool -> Bool
&& Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
16
            Bool -> Bool -> Bool
&& ([Double] -> Bool) -> [[Double]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
d) (Int -> Bool) -> ([Double] -> Int) -> [Double] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) [[Double]]
controls
            Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
weights Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n
            Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) [Double]
weights
            Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
knots Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
degree Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
            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]
knots (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
knots))
            Bool -> Bool -> Bool
&& [Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
degree Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< [Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
n
        )
        String
"invalid B-spline degree/knots/control points/positive weights"
    _ <- [Double] -> Either MathError [Double]
vector ([Double]
knots [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Double]]
controls [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
weights)
    Right (BSpline degree knots controls weights)

basis :: Int -> [Double] -> Int -> Double -> Double
basis :: Int -> [Double] -> Int -> Double -> Double
basis Int
p [Double]
knots Int
i Double
x
    | Int
p Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = if [Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
x Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< [Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) then Double
1 else Double
0
    | Bool
otherwise =
        let a :: Double
a = [Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
p) Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i; b :: Double
b = [Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
         in (if Double
a Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
0 else (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> [Double] -> Int -> Double -> Double
basis (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) [Double]
knots Int
i Double
x)
                Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (if Double
b Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
0 else ([Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> [Double] -> Int -> Double -> Double
basis (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) [Double]
knots (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double
x)

splineHomogeneous :: BSpline -> Double -> Either MathError ([Double], Double, [Double], Double)
splineHomogeneous :: BSpline
-> Double -> Either MathError ([Double], Double, [Double], Double)
splineHomogeneous (BSpline Int
p [Double]
knots [[Double]]
controls [Double]
ws) Double
x0 = do
    let n :: Int
n = [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
controls; hi :: Double
hi = [Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
n
    _ <- Extrapolation -> [Double] -> Double -> Either MathError Double
pointPolicy Extrapolation
Reject [[Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
p, Double
hi] Double
x0
    -- Evaluate the right endpoint using the exact limiting basis.
    let evalBasis Int
degree Int
i =
            if Double
x0 Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
hi
                then Int -> [Double] -> Int -> Double -> Double
basis Int
degree ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
negate ([Double] -> [Double]
forall a. [a] -> [a]
reverse [Double]
knots)) ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
knots Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
degree Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
i) (-Double
x0)
                else Int -> [Double] -> Int -> Double -> Double
basis Int
degree [Double]
knots Int
i Double
x0
        bs = [Int -> Int -> Double
evalBasis Int
p Int
i | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        db Int
i =
            if Int
p Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
                then Double
0
                else
                    let a :: Double
a = [Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
p) Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i; b :: Double
b = [Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
                     in (if Double
a Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
0 else Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
p Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Int -> Double
evalBasis (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
i)
                            Double -> Double -> Double
forall a. Num a => a -> a -> a
- (if Double
b Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
0 else Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
p Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Int -> Double
evalBasis (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1))
        numerator = [[Double]] -> [Double] -> [Double]
mv ([[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
tr [[Double]]
controls) ((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]
bs [Double]
ws)
        denominator = [Double] -> [Double] -> Double
dot [Double]
bs [Double]
ws
        dnumerator = [[Double]] -> [Double] -> [Double]
mv ([[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
tr [[Double]]
controls) ((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
(*) [Int -> Double
db Int
i | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] [Double]
ws)
        ddenominator = [Double] -> [Double] -> Double
dot [Int -> Double
db Int
i | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] [Double]
ws
    require (denominator > 0 && finite denominator) "zero/overflow rational spline denominator"
    _ <- vector (numerator ++ dnumerator ++ [ddenominator])
    Right (numerator, denominator, dnumerator, ddenominator)

-- | Evaluate rational Cox--de Boor basis functions and divide homogeneous coordinates by their positive total weight.
evaluateBSpline :: BSpline -> Double -> Either MathError [Double]
evaluateBSpline :: BSpline -> Double -> Either MathError [Double]
evaluateBSpline BSpline
s Double
x = do (v, w, _, _) <- BSpline
-> Double -> Either MathError ([Double], Double, [Double], Double)
splineHomogeneous BSpline
s Double
x; vector (map (/ w) v)

-- | Analytic first derivative of a rational B-spline using basis derivatives and the homogeneous quotient rule.
derivativeBSpline :: BSpline -> Double -> Either MathError [Double]
derivativeBSpline :: BSpline -> Double -> Either MathError [Double]
derivativeBSpline BSpline
s Double
x = do (v, w, dv, dw) <- BSpline
-> Double -> Either MathError ([Double], Double, [Double], Double)
splineHomogeneous BSpline
s Double
x; vector (zipWith (\Double
a Double
b -> (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
dw Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
w)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
w) dv v)

-- | One Boehm knot insertion in homogeneous coordinates; shape is preserved.
insertKnot :: BSpline -> Double -> Either MathError BSpline
insertKnot :: BSpline -> Double -> Either MathError BSpline
insertKnot (BSpline Int
p [Double]
knots [[Double]]
controls [Double]
ws) Double
u = do
    let n :: Int
n = [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
controls
        k :: Int
k = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((Double -> Bool) -> [Double] -> [Double]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
u) [Double]
knots) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
        s :: Int
s = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((Double -> Bool) -> [Double] -> [Double]
forall a. (a -> Bool) -> [a] -> [a]
filter (Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
u) [Double]
knots)
    Bool -> String -> Either MathError ()
require
        (Double -> Bool
finite Double
u Bool -> Bool -> Bool
&& Double
u Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> [Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
p Bool -> Bool -> Bool
&& Double
u Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< [Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
n Bool -> Bool -> Bool
&& Int
s Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
p)
        String
"knot insertion needs interior knot with multiplicity less than degree"
    let hc :: [[Double]]
hc = ([Double] -> Double -> [Double])
-> [[Double]] -> [Double] -> [[Double]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\[Double]
v Double
w -> (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
*) [Double]
v [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double
w]) [[Double]]
controls [Double]
ws
        at :: Int -> [Double]
at Int
i
            | Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
p = [[Double]]
hc [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i
            | Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
s Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 = [[Double]]
hc [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
            | Bool
otherwise =
                let a :: Double
a = (Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ([Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
p) Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
knots [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)
                 in (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
x Double
y -> Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
1 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]]
hc [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) ([[Double]]
hc [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1))
        next :: [[Double]]
next = (Int -> [Double]) -> [Int] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map Int -> [Double]
at [Int
0 .. Int
n]
        ww :: [Double]
ww = ([Double] -> Double) -> [[Double]] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map [Double] -> Double
forall a. HasCallStack => [a] -> a
last [[Double]]
next
        cc :: [[Double]]
cc = ([Double] -> Double -> [Double])
-> [[Double]] -> [Double] -> [[Double]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\[Double]
v Double
w -> (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
w) ([Double] -> [Double]
forall a. HasCallStack => [a] -> [a]
init [Double]
v)) [[Double]]
next [Double]
ww
    Int
-> [Double] -> [[Double]] -> [Double] -> Either MathError BSpline
bSpline Int
p (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [Double]
knots [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double
u] [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [Double]
knots) [[Double]]
cc [Double]
ww

-- | Positive-weight rational Bezier curve on [0,1], represented exactly as an open clamped B-spline (degree <=16).
rationalBezier :: [[Double]] -> [Double] -> Double -> Either MathError [Double]
rationalBezier :: [[Double]] -> [Double] -> Double -> Either MathError [Double]
rationalBezier [[Double]]
controls [Double]
ws Double
t = do
    Bool -> String -> Either MathError ()
require
        ([[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
controls Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
1 Bool -> Bool -> Bool
&& [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
controls Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
17 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
t Bool -> Bool -> Bool
&& Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 Bool -> Bool -> Bool
&& Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1)
        String
"invalid Bezier degree/parameter"
    spline <-
        Int
-> [Double] -> [[Double]] -> [Double] -> Either MathError BSpline
bSpline
            ([[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
controls Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
            (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate ([[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
controls) Double
0 [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate ([[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
controls) Double
1)
            [[Double]]
controls
            [Double]
ws
    evaluateBSpline spline t

-- | Polynomial Bezier curve with unit weights; controls are vectors with 1..16 coordinates.
bezier :: [[Double]] -> Double -> Either MathError [Double]
bezier :: [[Double]] -> Double -> Either MathError [Double]
bezier [[Double]]
controls = [[Double]] -> [Double] -> Double -> Either MathError [Double]
rationalBezier [[Double]]
controls (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate ([[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
controls) Double
1)

{- | Interpolating cubic Hermite segment sequence with Kochanek--Bartels
tension, continuity and bias in [-1,1]; endpoint secants are repeated.
-}
kochanekBartels ::
    Double -> Double -> Double -> [(Double, Double)] -> Either MathError PiecewisePolynomial
kochanekBartels :: Double
-> Double
-> Double
-> [(Double, Double)]
-> Either MathError PiecewisePolynomial
kochanekBartels Double
tension Double
continuity Double
bias [(Double, Double)]
points = do
    (xs, ys) <- [(Double, Double)] -> Either MathError ([Double], [Double])
samples [(Double, Double)]
points
    require
        (length xs >= 2 && all (\Double
x -> Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& Double -> Double
forall a. Num a => a -> a
abs Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1) [tension, continuity, bias])
        "invalid Kochanek-Bartels parameters"
    let n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs
        s = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Fractional a => a -> a -> a
(/) ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
ys) [Double]
ys) ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
xs) [Double]
xs)
        before Int
i = [Double]
s [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
0 (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
        after Int
i = [Double]
s [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int -> Int -> Int
forall a. Ord a => a -> a -> a
min (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2) Int
i
        out Int
i =
            (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
tension)
                Double -> Double -> Double
forall a. Num a => a -> a -> a
* ((Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
continuity) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
bias) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
before 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
continuity) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
bias) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
after Int
i)
                Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
        incoming Int
i =
            (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
tension)
                Double -> Double -> Double
forall a. Num a => a -> a -> a
* ((Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
continuity) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
bias) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
before 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
continuity) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
bias) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
after Int
i)
                Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
    piecewisePolynomial
        xs
        [ let h = [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i; sec = [Double]
s [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i; a = Int -> Double
out Int
i; b = Int -> Double
incoming (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
           in [ys !! i, a, (3 * sec - 2 * a - b) / h, (a + b - 2 * sec) / (h * h)]
        | i <- [0 .. n - 2]
        ]

-- | Kochanek--Bartels interpolation with zero continuity/bias and supplied tension in [-1,1].
cardinalSpline :: Double -> [(Double, Double)] -> Either MathError PiecewisePolynomial
cardinalSpline :: Double
-> [(Double, Double)] -> Either MathError PiecewisePolynomial
cardinalSpline Double
t = Double
-> Double
-> Double
-> [(Double, Double)]
-> Either MathError PiecewisePolynomial
kochanekBartels Double
t Double
0 Double
0

-- | Cardinal cubic interpolation with zero tension and repeated endpoint secants.
catmullRom :: [(Double, Double)] -> Either MathError PiecewisePolynomial
catmullRom :: [(Double, Double)] -> Either MathError PiecewisePolynomial
catmullRom = Double
-> [(Double, Double)] -> Either MathError PiecewisePolynomial
cardinalSpline Double
0

-- | Componentwise adaptive integral of a B-spline between in-domain limits, with estimated quadrature errors.
curveIntegral :: C.Tolerance -> BSpline -> Double -> Double -> Either MathError [C.Estimate Double]
curveIntegral :: Tolerance
-> BSpline
-> Double
-> Double
-> Either MathError [Estimate Double]
curveIntegral Tolerance
tolerance s :: BSpline
s@(BSpline Int
_ [Double]
_ [[Double]]
controls [Double]
_) Double
a Double
b = do
    _ <- BSpline -> Double -> Either MathError [Double]
evaluateBSpline BSpline
s Double
a
    _ <- evaluateBSpline s b
    mapM
        (\Int
i -> Tolerance
-> (Double -> Double)
-> Double
-> Double
-> Either MathError (Estimate Double)
C.adaptiveIntegrate Tolerance
tolerance (\Double
t -> (MathError -> Double)
-> ([Double] -> Double) -> Either MathError [Double] -> Double
forall a c b. (a -> c) -> (b -> c) -> Either a b -> c
either (Double -> MathError -> Double
forall a b. a -> b -> a
const (Double
0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0)) ([Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) (BSpline -> Double -> Either MathError [Double]
evaluateBSpline BSpline
s Double
t)) Double
a Double
b)
        [0 .. length (controls !! 0) - 1]

-- | Adaptive integral of the Euclidean analytic derivative norm; limits must increase and remain in the spline domain.
arcLength :: C.Tolerance -> BSpline -> Double -> Double -> Either MathError (C.Estimate Double)
arcLength :: Tolerance
-> BSpline
-> Double
-> Double
-> Either MathError (Estimate Double)
arcLength Tolerance
tolerance BSpline
spline Double
a Double
b = do
    Bool -> String -> Either MathError ()
require (Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
b) String
"arc length limits must be ordered"
    _ <- BSpline -> Double -> Either MathError [Double]
evaluateBSpline BSpline
spline Double
a
    _ <- evaluateBSpline spline b
    C.adaptiveIntegrate tolerance (\Double
t -> (MathError -> Double)
-> ([Double] -> Double) -> Either MathError [Double] -> Double
forall a c b. (a -> c) -> (b -> c) -> Either a b -> c
either (Double -> MathError -> Double
forall a b. a -> b -> a
const (Double
0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0)) [Double] -> Double
norm (BSpline -> Double -> Either MathError [Double]
derivativeBSpline BSpline
spline Double
t)) a b

{- | Rectilinear multilinear interpolation; values are row-major/last axis fast.
Nearest uses Euclidean distance to grid nodes; extrapolation policy per axis.
-}
gridInterpolate ::
    InterpolationMethod ->
    Extrapolation ->
    [[Double]] ->
    [Double] ->
    [Double] ->
    Either MathError Double
gridInterpolate :: InterpolationMethod
-> Extrapolation
-> [[Double]]
-> [Double]
-> [Double]
-> Either MathError Double
gridInterpolate InterpolationMethod
method Extrapolation
policy [[Double]]
axes [Double]
values [Double]
point = do
    let dims :: [Int]
dims = ([Double] -> Int) -> [[Double]] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
axes; d :: Int
d = [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
axes
    Bool -> String -> Either MathError ()
require
        ( Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0
            Bool -> Bool -> Bool
&& Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
16
            Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
point Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
d
            Bool -> Bool -> Bool
&& (Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
2) [Int]
dims
            Bool -> Bool -> Bool
&& [Integer] -> Integer
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ((Int -> Integer) -> [Int] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Integer
forall a. Integral a => a -> Integer
toInteger [Int]
dims) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
1048576
            Bool -> Bool -> Bool
&& [Int] -> Int
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Int]
dims 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]
values
            Bool -> Bool -> Bool
&& (InterpolationMethod
method InterpolationMethod -> InterpolationMethod -> Bool
forall a. Eq a => a -> a -> Bool
== InterpolationMethod
Linear Bool -> Bool -> Bool
|| InterpolationMethod
method InterpolationMethod -> InterpolationMethod -> Bool
forall a. Eq a => a -> a -> Bool
== InterpolationMethod
Nearest)
        )
        String
"invalid grid interpolation shape/method"
    ([Double] -> Either MathError ([Double], [Double]))
-> [[Double]] -> Either MathError ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\[Double]
a -> [(Double, Double)] -> Either MathError ([Double], [Double])
samples ([Double] -> [Double] -> [(Double, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
a (Double -> [Double]
forall a. a -> [a]
repeat Double
0))) [[Double]]
axes
    _ <- [Double] -> Either MathError [Double]
vector [Double]
values
    p <- sequence (zipWith (pointPolicy policy) axes point)
    let intervals = [Int -> Int -> Int
forall a. Ord a => a -> a -> a
min ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
a Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2) (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
0 ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((Double -> Bool) -> [Double] -> [Double]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
x) [Double]
a) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)) | ([Double]
a, Double
x) <- [[Double]] -> [Double] -> [([Double], Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [[Double]]
axes [Double]
p]
        ts = [(Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
a [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ([Double]
a [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
a [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) | ([Double]
a, Double
x, Int
i) <- [[Double]] -> [Double] -> [Int] -> [([Double], Double, Int)]
forall a b c. [a] -> [b] -> [c] -> [(a, b, c)]
zip3 [[Double]]
axes [Double]
p [Int]
intervals]
        strides = Int -> [Int] -> [Int]
forall a. Int -> [a] -> [a]
drop Int
1 ((Int -> Int -> Int) -> Int -> [Int] -> [Int]
forall a b. (a -> b -> b) -> b -> [a] -> [b]
scanr Int -> Int -> Int
forall a. Num a => a -> a -> a
(*) Int
1 [Int]
dims)
        corner [Int]
bits =
            let indices :: [Int]
indices = (Int -> Int -> Int) -> [Int] -> [Int] -> [Int]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Int -> Int -> Int
forall a. Num a => a -> a -> a
(+) [Int]
intervals [Int]
bits
                weight :: Double
weight = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ((Double -> Int -> Double) -> [Double] -> [Int] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
t Int
bit -> if Int
bit Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t else Double
t) [Double]
ts [Int]
bits)
             in Double
weight Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
values [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! [Int] -> Int
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Int -> Int -> Int) -> [Int] -> [Int] -> [Int]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Int -> Int -> Int
forall a. Num a => a -> a -> a
(*) [Int]
indices [Int]
strides))
    if method == Nearest
        then
            checked
                (values !! sum (zipWith (*) (zipWith (\Int
i Double
t -> Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ if Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0.5 then Int
0 else Int
1) intervals ts) strides))
        else checked (sum (map corner (sequence (replicate d [0, 1]))))

scattered :: [[Double]] -> [Double] -> Either MathError Int
scattered :: [[Double]] -> [Double] -> Either MathError Int
scattered [[Double]]
nodes [Double]
ys = do
    let n :: Int
n = [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
nodes; d :: Int
d = if [[Double]] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [[Double]]
nodes then Int
0 else [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([[Double]]
nodes [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
0)
    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
1024 Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n Bool -> Bool -> Bool
&& Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
d Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
16 Bool -> Bool -> Bool
&& ([Double] -> Bool) -> [[Double]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
d) (Int -> Bool) -> ([Double] -> Int) -> [Double] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) [[Double]]
nodes)
        String
"invalid scattered-data shape"
    _ <- [Double] -> Either MathError [Double]
vector ([[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Double]]
nodes [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
ys)
    require
        (and [nodes !! i /= nodes !! j | i <- [0 .. n - 1], j <- [0 .. i - 1]])
        "duplicate scattered nodes"
    Right d
distance :: [Double] -> [Double] -> Double
distance :: [Double] -> [Double] -> Double
distance [Double]
x [Double]
y = [Double] -> Double
norm ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Double]
x [Double]
y)
kernel :: RadialKernel -> Double -> Double
kernel :: RadialKernel -> Double -> Double
kernel RadialKernel
method Double
r = case RadialKernel
method of
    Gaussian Double
e -> Double -> Double
forall a. Floating a => a -> a
exp (-(Double
e Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
r) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int))
    Multiquadric Double
e -> Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
e Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
r) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int))
    InverseMultiquadric Double
e -> Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
e Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
r) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int))
    RadialKernel
ThinPlate -> if Double
r Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
0 else Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
r
    RadialKernel
Cubic -> Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
r

-- | Solve the augmented radial/affine system with nonnegative diagonal smoothing; singular or affinely degenerate systems fail.
fitRBF :: RadialKernel -> Double -> [[Double]] -> [Double] -> Either MathError RBF
fitRBF :: RadialKernel
-> Double -> [[Double]] -> [Double] -> Either MathError RBF
fitRBF RadialKernel
method Double
smoothing [[Double]]
nodes [Double]
ys = do
    d <- [[Double]] -> [Double] -> Either MathError Int
scattered [[Double]]
nodes [Double]
ys
    require (finite smoothing && smoothing >= 0) "invalid RBF smoothing"
    case method of
        Gaussian Double
e -> Double -> Either MathError ()
validScale Double
e
        Multiquadric Double
e -> Double -> Either MathError ()
validScale Double
e
        InverseMultiquadric Double
e -> Double -> Either MathError ()
validScale Double
e
        RadialKernel
_ -> () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    let n = [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
nodes
        p = ([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map (Double
1 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
:) [[Double]]
nodes
        a =
            [ [ RadialKernel -> Double -> Double
kernel RadialKernel
method ([Double] -> [Double] -> Double
distance [Double]
x [Double]
y) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then Double
smoothing else Double
0) | (Int
j, [Double]
y) <- [Int] -> [[Double]] -> [(Int, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [[Double]]
nodes
              ]
                [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [[Double]]
p [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i
            | (Int
i, [Double]
x) <- [Int] -> [[Double]] -> [(Int, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [[Double]]
nodes
            ]
                [[Double]] -> [[Double]] -> [[Double]]
forall a. [a] -> [a] -> [a]
++ [[Double]
row [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate (Int
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double
0 | [Double]
row <- [[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
tr [[Double]]
p]
    coeffs <- solve a (ys ++ replicate (d + 1) 0)
    Right (RBF method nodes (take n coeffs) (drop n coeffs))
  where
    validScale :: Double -> Either MathError ()
validScale Double
e = Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
e Bool -> Bool -> Bool
&& Double
e Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"RBF shape must be positive"

-- | Evaluate a fitted scattered model at one finite point with the original coordinate dimension.
evaluateRBF :: RBF -> [Double] -> Either MathError Double
evaluateRBF :: RBF -> [Double] -> Either MathError Double
evaluateRBF (RBF RadialKernel
method [[Double]]
nodes [Double]
weights [Double]
affine) [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
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([[Double]]
nodes [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
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
"RBF query shape/nonfinite"
    Double -> Either MathError Double
checked
        ([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
w [Double]
c -> Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
* RadialKernel -> Double -> Double
kernel RadialKernel
method ([Double] -> [Double] -> Double
distance [Double]
x [Double]
c)) [Double]
weights [[Double]]
nodes) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double] -> [Double] -> Double
dot [Double]
affine (Double
1 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [Double]
x))

-- | Normalized inverse-distance interpolation with positive power; coincident queries return the exact observation.
inverseDistance :: Double -> [[Double]] -> [Double] -> [Double] -> Either MathError Double
inverseDistance :: Double
-> [[Double]] -> [Double] -> [Double] -> Either MathError Double
inverseDistance Double
power [[Double]]
nodes [Double]
ys [Double]
point = do
    d <- [[Double]] -> [Double] -> Either MathError Int
scattered [[Double]]
nodes [Double]
ys
    require
        (finite power && power > 0 && length point == d && all finite point)
        "invalid inverse distance query/power"
    case lookup point (zip nodes ys) of
        Just Double
y -> Double -> Either MathError Double
forall a b. b -> Either a b
Right Double
y
        Maybe Double
Nothing ->
            let distances :: [Double]
distances = ([Double] -> Double) -> [[Double]] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ([Double] -> [Double] -> Double
distance [Double]
point) [[Double]]
nodes
                nearest :: Double
nearest = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Double]
distances
                weights :: [Double]
weights = (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\Double
r -> (Double
nearest Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
r) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
power) [Double]
distances
             in Double -> Either MathError Double
checked ([Double] -> [Double] -> Double
dot [Double]
weights [Double]
ys Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
weights)

-- | Classical inverse-square Shepard interpolation in the supplied Euclidean coordinates.
shepard :: [[Double]] -> [Double] -> [Double] -> Either MathError Double
shepard :: [[Double]] -> [Double] -> [Double] -> Either MathError Double
shepard = Double
-> [[Double]] -> [Double] -> [Double] -> Either MathError Double
inverseDistance Double
2

{- | Ordinary kriging with stationary exponential covariance exp(-distance/range),
optional independent nugget, unbiased constant mean. Returns mean and variance.
-}
kriging ::
    Double -> Double -> [[Double]] -> [Double] -> [Double] -> Either MathError (Double, Double)
kriging :: Double
-> Double
-> [[Double]]
-> [Double]
-> [Double]
-> Either MathError (Double, Double)
kriging Double
range Double
nugget [[Double]]
nodes [Double]
ys [Double]
point = do
    d <- [[Double]] -> [Double] -> Either MathError Int
scattered [[Double]]
nodes [Double]
ys
    require
        (finite range && range > 0 && finite nugget && nugget >= 0 && length point == d && all finite point)
        "invalid kriging parameters/query"
    let n = [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
nodes
        c = ([Double] -> Double) -> [[Double]] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (\[Double]
x -> Double -> Double
forall a. Floating a => a -> a
exp (-[Double] -> [Double] -> Double
distance [Double]
x [Double]
point Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
range)) [[Double]]
nodes
        a =
            [ [Double -> Double
forall a. Floating a => a -> a
exp (-[Double] -> [Double] -> Double
distance [Double]
x [Double]
y Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
range) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then Double
nugget else Double
0) | (Int
j, [Double]
y) <- [Int] -> [[Double]] -> [(Int, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [[Double]]
nodes]
                [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double
1]
            | (Int
i, [Double]
x) <- [Int] -> [[Double]] -> [(Int, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [[Double]]
nodes
            ]
                [[Double]] -> [[Double]] -> [[Double]]
forall a. [a] -> [a] -> [a]
++ [Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
n Double
1 [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double
0]]
    w <- solve a (c ++ [1])
    mean <- checked (dot (take n w) ys)
    variance <- checked (1 - dot (take n w) c - last w)
    require (variance >= -1e-10) "kriging covariance produced negative variance"
    Right (mean, max 0 variance)

-- | Full-column-rank polynomial least squares via pivoted QR, returning ascending coefficients and maximum observation error.
polynomialFit :: Int -> [(Double, Double)] -> Either MathError Approximation
polynomialFit :: Int -> [(Double, Double)] -> Either MathError Approximation
polynomialFit Int
degree [(Double, Double)]
points = do
    (xs, ys) <- [(Double, Double)] -> Either MathError ([Double], [Double])
samples [(Double, Double)]
points
    require (degree >= 0 && degree <= 64 && degree < length xs) "invalid polynomial fit degree"
    design <-
        M.fromListChecked (length xs) (degree + 1) (concat [[x ^ j | j <- [0 .. degree]] | x <- xs])
    cs <- QR.leastSquares design ys
    err <- interpolationError (polyEval cs) points
    Right (Approximation cs err 1)

-- | Chebyshev-root interpolation on a finite interval; report maximum error on a separate 513-point uniform validation grid.
chebyshevApproximate ::
    Int -> (Double -> Double) -> Double -> Double -> Either MathError Approximation
chebyshevApproximate :: Int
-> (Double -> Double)
-> Double
-> Double
-> Either MathError Approximation
chebyshevApproximate Int
degree Double -> Double
f Double
a Double
b = do
    Bool -> String -> Either MathError ()
require
        (Int
degree Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
degree Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
64 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
a Bool -> Bool -> Bool
&& Double -> Bool
finite Double
b Bool -> Bool -> Bool
&& Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
b Bool -> Bool -> Bool
&& Double -> Bool
finite (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a))
        String
"invalid Chebyshev interval/degree"
    let n :: Int
n = Int
degree Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
        theta :: a -> a
theta a
j = a
forall a. Floating a => a
pi a -> a -> a
forall a. Num a => a -> a -> a
* (a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
j a -> a -> a
forall a. Num a => a -> a -> a
+ a
0.5) a -> a -> a
forall a. Fractional a => a -> a -> a
/ Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
        xs :: [Double]
xs = [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
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
cos (Int -> Double
forall {a} {a}. (Floating a, Integral a) => a -> a
theta Int
j)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
    ys <- [Double] -> Either MathError [Double]
vector ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
f [Double]
xs)
    let cs =
            [ (if Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then Double
1 else 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 [[Double]
ys [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 (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall {a} {a}. (Floating a, Integral a) => a -> a
theta Int
j) | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
            | Int
k <- [Int
0 .. Int
degree]
            ]
        grid = [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
j Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
512 | Int
j <- [Int
0 .. Int
512 :: Int]]
    fitted <- mapM (evaluateChebyshev cs a b) grid
    _ <- vector grid
    truth <- vector (map f grid)
    Right (Approximation cs (maximum (map abs (zipWith (-) truth fitted))) n)

-- | Evaluate ascending Chebyshev coefficients after mapping the original interval to [-1,1].
evaluateChebyshev :: [Double] -> Double -> Double -> Double -> Either MathError Double
evaluateChebyshev :: [Double] -> Double -> Double -> Double -> Either MathError Double
evaluateChebyshev [Double]
cs Double
a Double
b Double
x = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
a Bool -> Bool -> Bool
&& Double -> Bool
finite Double
b Bool -> Bool -> Bool
&& Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
b Bool -> Bool -> Bool
&& Double -> Bool
finite Double
x Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
cs) String
"invalid Chebyshev evaluation"
    let z :: Double
z = Double
2 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. Fractional 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
1; ts :: [Double]
ts = Double
1 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: Double
z Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
u Double
v -> Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
u) [Double]
ts (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
ts)
    Double -> Either MathError Double
checked ([Double] -> [Double] -> Double
dot [Double]
cs (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
cs) [Double]
ts))

-- | Padé [m/n] from ascending Taylor coefficients, q_0=1.
pade :: Int -> Int -> [Double] -> Either MathError ([Double], [Double])
pade :: Int -> Int -> [Double] -> Either MathError ([Double], [Double])
pade Int
m Int
n [Double]
series = do
    Bool -> String -> Either MathError ()
require
        (Int
m 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
0 Bool -> Bool -> Bool
&& Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
64 Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
series Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
series)
        String
"invalid Pade orders/series"
    let a :: Int -> Double
a Int
k = if Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 then Double
0 else [Double]
series [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
k
    qs <-
        if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
            then [Double] -> Either MathError [Double]
forall a b. b -> Either a b
Right []
            else [[Double]] -> [Double] -> Either MathError [Double]
solve [[Int -> Double
a (Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
j) | Int
j <- [Int
1 .. Int
n]] | Int
i <- [Int
1 .. Int
n]] [-Int -> Double
a (Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
i) | Int
i <- [Int
1 .. Int
n]]
    let q = Double
1 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [Double]
qs; p = [[Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [[Double]
q [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
a (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
j) | Int
j <- [Int
0 .. Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
k Int
n]] | Int
k <- [Int
0 .. Int
m]]
    _ <- vector (p ++ q)
    Right (p, q)

{- | Linearized rational least squares y*q(x)-p(x)=0 with q_0=1;
poles at observations are rejected. Error is the actual rational residual.
-}
rationalFit :: Int -> Int -> [(Double, Double)] -> Either MathError ([Double], [Double], Double)
rationalFit :: Int
-> Int
-> [(Double, Double)]
-> Either MathError ([Double], [Double], Double)
rationalFit Int
m Int
n [(Double, Double)]
points = do
    (xs, ys) <- [(Double, Double)] -> Either MathError ([Double], [Double])
samples [(Double, Double)]
points
    require (m >= 0 && n >= 0 && m + n <= 64 && length xs >= m + n + 1) "invalid rational fit degrees"
    a <-
        M.fromListChecked
            (length xs)
            (m + n + 1)
            (concat [[x ^ j | j <- [0 .. m]] ++ [-y * x ^ j | j <- [1 .. n]] | (x, y) <- points])
    cs <- QR.leastSquares a ys
    let p = Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take (Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [Double]
cs; q = Double
1 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop (Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [Double]
cs
    require (all ((/= 0) . polyEval q) xs) "rational fit has an observation pole"
    err <- interpolationError (\Double
x -> [Double] -> Double -> Double
polyEval [Double]
p Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> Double -> Double
polyEval [Double]
q Double
x) points
    Right (p, q, err)

{- | Discrete Remez exchange on an explicit dense grid, 1..64 degree. Success
requires alternating extrema and stable equioscillation; returned error is
the maximum on this grid only, not a certified continuous minimax bound.
-}
remez ::
    Int ->
    Int ->
    Double ->
    Int ->
    (Double -> Double) ->
    Double ->
    Double ->
    Either MathError Approximation
remez :: Int
-> Int
-> Double
-> Int
-> (Double -> Double)
-> Double
-> Double
-> Either MathError Approximation
remez Int
degree Int
gridCount Double
tolerance Int
budget Double -> Double
f Double
a Double
b = do
    Bool -> String -> Either MathError ()
require
        ( Int
degree Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0
            Bool -> Bool -> Bool
&& Int
degree Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
32
            Bool -> Bool -> Bool
&& Int
gridCount Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
8 Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
degree Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2)
            Bool -> Bool -> Bool
&& Int
gridCount Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
65536
            Bool -> Bool -> Bool
&& Double -> Bool
finite Double
tolerance
            Bool -> Bool -> Bool
&& Double
tolerance Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0
            Bool -> Bool -> Bool
&& Int
budget Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0
            Bool -> Bool -> Bool
&& Double -> Bool
finite Double
a
            Bool -> Bool -> Bool
&& Double -> Bool
finite Double
b
            Bool -> Bool -> Bool
&& Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
b
            Bool -> Bool -> Bool
&& Double -> Bool
finite (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a)
        )
        String
"invalid Remez configuration"
    let grid :: [Double]
grid = [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
j Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
gridCount Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) | Int
j <- [Int
0 .. Int
gridCount Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
    _ <- [Double] -> Either MathError [Double]
vector [Double]
grid
    truth <- vector (map f grid)
    let initial =
            [ Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
round
                ( (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
gridCount Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) :: Double)
                    Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 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
degree Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)))
                    Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
                )
            | Int
j <- [Int
0 .. Int
degree Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1]
            ]
        go Int
iteration [Int]
indices = do
            let rows :: [[Double]]
rows =
                    [ [([Double]
grid [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
k | Int
k <- [Int
0 .. Int
degree]] [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [if Int -> Bool
forall a. Integral a => a -> Bool
even Int
j then Double
1 else -Double
1]
                    | (Int
j, Int
i) <- [Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [Int]
indices
                    ]
            cs <- [[Double]] -> [Double] -> Either MathError [Double]
solve [[Double]]
rows [[Double]
truth [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i | Int
i <- [Int]
indices]
            let p = Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take (Int
degree Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [Double]
cs
                err = [[Double]
truth [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double] -> Double -> Double
polyEval [Double]
p ([Double]
grid [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) | Int
i <- [Int
0 .. Int
gridCount Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
                local =
                    Int
0
                        Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: [ Int
i
                          | Int
i <- [Int
1 .. Int
gridCount Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2]
                          , Double -> Double
forall a. Num a => a -> a
abs ([Double]
err [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double -> Double
forall a. Num a => a -> a
abs ([Double]
err [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1))
                          , Double -> Double
forall a. Num a => a -> a
abs ([Double]
err [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double -> Double
forall a. Num a => a -> a
abs ([Double]
err [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1))
                          ]
                        [Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++ [Int
gridCount Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                alternate =
                    ([Int] -> Int -> [Int]) -> [Int] -> [Int] -> [Int]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl
                        ( \[Int]
acc Int
i -> case [Int] -> [Int]
forall a. [a] -> [a]
reverse [Int]
acc of
                            [] -> [Int
i]
                            Int
j : [Int]
_ ->
                                if [Double]
err [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
err [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0
                                    then if Double -> Double
forall a. Num a => a -> a
abs ([Double]
err [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double -> Double
forall a. Num a => a -> a
abs ([Double]
err [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) then [Int] -> [Int]
forall a. HasCallStack => [a] -> [a]
init [Int]
acc [Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++ [Int
i] else [Int]
acc
                                    else [Int]
acc [Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++ [Int
i]
                        )
                        []
                        [Int]
local
                trim [Int]
is
                    | [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
is Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
degree Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2 = [Int]
is
                    | Double -> Double
forall a. Num a => a -> a
abs ([Double]
err [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! ([Int]
is [Int] -> Int -> Int
forall a. HasCallStack => [a] -> Int -> a
!! Int
0)) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double -> Double
forall a. Num a => a -> a
abs ([Double]
err [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! [Int] -> Int
forall a. HasCallStack => [a] -> a
last [Int]
is) = [Int] -> [Int]
trim (Int -> [Int] -> [Int]
forall a. Int -> [a] -> [a]
drop Int
1 [Int]
is)
                    | Bool
otherwise = [Int] -> [Int]
trim ([Int] -> [Int]
forall a. HasCallStack => [a] -> [a]
init [Int]
is)
                next = [Int] -> [Int]
trim [Int]
alternate
                maximumError = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
abs [Double]
err)
                ripple = Double -> Double
forall a. Num a => a -> a
abs ([Double] -> Double
forall a. HasCallStack => [a] -> a
last [Double]
cs)
            _ <- vector (p ++ err)
            if maximumError <= ripple + tolerance || maximumError <= tolerance
                then Right (Approximation p maximumError iteration)
                else
                    if iteration >= budget || length next < degree + 2
                        then Left (NonConvergence "Remez exchange" iteration)
                        else go (iteration + 1) next
    go 1 initial

{- | Damped Gauss--Newton least squares with central-difference Jacobian, strict
residual/step convergence and finite total callbacks; returns fitted parameters.
-}
nonlinearFit ::
    Int ->
    Double ->
    (Double -> [Double] -> Double) ->
    [(Double, Double)] ->
    [Double] ->
    Either MathError Approximation
nonlinearFit :: Int
-> Double
-> (Double -> [Double] -> Double)
-> [(Double, Double)]
-> [Double]
-> Either MathError Approximation
nonlinearFit Int
budget Double
tolerance Double -> [Double] -> Double
model [(Double, Double)]
points [Double]
initial = do
    _ <- [(Double, Double)] -> Either MathError ([Double], [Double])
samples [(Double, Double)]
points
    require
        ( budget > 0
            && finite tolerance
            && tolerance > 0
            && not (null initial)
            && length initial <= 32
            && length points >= length initial
        )
        "invalid nonlinear fit"
    _ <- vector initial
    let residual [Double]
ps = [Double] -> Either MathError [Double]
vector [Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> [Double] -> Double
model Double
x [Double]
ps | (Double
x, Double
y) <- [(Double, Double)]
points]
        go Int
iteration Double
damping [Double]
ps = do
            r <- [Double] -> Either MathError [Double]
residual [Double]
ps
            columns <-
                mapM
                    ( \Int
j -> do
                        let h :: Double
h = Double
1e-5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1 (Double -> Double
forall a. Num a => a -> a
abs ([Double]
ps [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j))
                            plus :: [Double]
plus = Int -> Double -> [Double] -> [Double]
forall a. Int -> a -> [a] -> [a]
replace Int
j ([Double]
ps [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h) [Double]
ps
                            minus :: [Double]
minus = Int -> Double -> [Double] -> [Double]
forall a. Int -> a -> [a] -> [a]
replace Int
j ([Double]
ps [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
h) [Double]
ps
                        _ <- [Double] -> Either MathError [Double]
vector ([Double]
plus [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
minus)
                        vector [(model x plus - model x minus) / (2 * h) | (x, _) <- points]
                    )
                    [0 .. length ps - 1]
            let jac = [[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
tr [[Double]]
columns
                jtj = [[Double]] -> [[Double]] -> [[Double]]
mm [[Double]]
columns [[Double]]
jac
                a = ([Double] -> [Double] -> [Double])
-> [[Double]] -> [[Double]] -> [[Double]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith ((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]]
jtj (([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double
damping Double -> Double -> Double
forall a. Num a => a -> a -> a
*)) (Int -> [[Double]]
eye ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ps)))
            step <- solve a (mv columns r)
            next <- vector (zipWith (+) ps step)
            rr <- residual next
            let old = [Double] -> Double
norm [Double]
r; new = [Double] -> Double
norm [Double]
rr
            if old <= tolerance || new <= old && norm step <= tolerance * (1 + norm ps)
                then Right (Approximation (if old <= tolerance then ps else next) (min old new) iteration)
                else
                    if iteration >= budget
                        then Left (NonConvergence "nonlinear fit" iteration)
                        else
                            if new < old
                                then go (iteration + 1) (max 1e-15 (damping / 3)) next
                                else go (iteration + 1) (damping * 10) ps
    go 1 1e-3 initial

-- | Maximum absolute discrepancy at the supplied nonempty observations; nonfinite callback results fail.
interpolationError :: (Double -> Double) -> [(Double, Double)] -> Either MathError Double
interpolationError :: (Double -> Double) -> [(Double, Double)] -> Either MathError Double
interpolationError Double -> Double
f [(Double, Double)]
points = do
    Bool -> String -> Either MathError ()
require (Bool -> Bool
not ([(Double, Double)] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(Double, Double)]
points)) String
"error estimate needs observations"
    errors <- [Double] -> Either MathError [Double]
vector [Double -> Double
forall a. Num a => a -> a
abs (Double -> Double
f Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
y) | (Double
x, Double
y) <- [(Double, Double)]
points]
    Right (maximum errors)